From d2f38709bdd62660ab47f13b2b6b25cfb083e649 Mon Sep 17 00:00:00 2001 From: apcraig Date: Wed, 24 Aug 2022 14:08:24 -0600 Subject: [PATCH] clean up trailing blanks --- columnphysics/icepack_shortwave.F90 | 550 ++++++++++++++-------------- 1 file changed, 275 insertions(+), 275 deletions(-) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 0dd424304..737e5b01e 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -4,26 +4,26 @@ ! snow over ice, bare ice and ponded ice. ! ! Presently, two methods are included: -! (1) CCSM3 -! (2) Delta-Eddington +! (1) CCSM3 +! (2) Delta-Eddington ! as two distinct routines. ! Either can be called from the ice driver. ! ! The Delta-Eddington method is described here: ! -! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple -! Scattering Parameterization for Solar Radiation in the Sea Ice -! Component of the Community Climate System Model, NCAR Technical +! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple +! Scattering Parameterization for Solar Radiation in the Sea Ice +! Component of the Community Climate System Model, NCAR Technical ! Note NCAR/TN-472+STR February 2007 ! ! name: originally ice_albedo ! ! authors: Bruce P. Briegleb, NCAR ! Elizabeth C. Hunke and William H. Lipscomb, LANL -! 2005, WHL: Moved absorbed_solar from icepack_therm_vertical to this +! 2005, WHL: Moved absorbed_solar from icepack_therm_vertical to this ! module and changed name from ice_albedo ! 2006, WHL: Added Delta Eddington routines from Bruce Briegleb -! 2006, ECH: Changed data statements in Delta Eddington routines (no +! 2006, ECH: Changed data statements in Delta Eddington routines (no ! longer hardwired) ! Converted to free source form (F90) ! 2007, BPB: Completely updated Delta-Eddington code, so that: @@ -89,7 +89,7 @@ module icepack_shortwave hpmin = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m) hp0 = 0.200_dbl_kind ! pond depth below which transition to bare ice - real (kind=dbl_kind), parameter :: & + real (kind=dbl_kind), parameter :: & exp_argmax = c10 ! maximum argument of exponential !======================================================================= @@ -339,19 +339,19 @@ subroutine shortwave_ccsm3 (aicen, vicen, & alidrni = albocn alvdfni = albocn alidfni = albocn - + alvdrns = albocn alidrns = albocn alvdfns = albocn alidfns = albocn - + alvdrn(n) = albocn alidrn(n) = albocn alvdfn(n) = albocn alidfn(n) = albocn - + albin(n) = c0 - albsn(n) = c0 + albsn(n) = c0 fswsfc(n) = c0 fswint(n) = c0 @@ -502,13 +502,13 @@ subroutine compute_albedos (aicen, vicen, & alidrn , & ! near-ir, direct, avg (fraction) alvdfn , & ! visible, diffuse, avg (fraction) alidfn , & ! near-ir, diffuse, avg (fraction) - albin , & ! bare ice - albsn ! snow + albin , & ! bare ice + albsn ! snow ! local variables real (kind=dbl_kind), parameter :: & - dT_melt = c1 , & ! change in temp to give dalb_mlt + dT_melt = c1 , & ! change in temp to give dalb_mlt ! albedo change dalb_mlt = -0.075_dbl_kind, & ! albedo change per dT_melt change ! in temp for ice @@ -536,7 +536,7 @@ subroutine compute_albedos (aicen, vicen, & !----------------------------------------------------------------- hi = vicen / aicen - hs = vsnon / aicen + hs = vsnon / aicen ! bare ice, thickness dependence fh = min(atan(hi*c4)/fhtan,c1) @@ -590,9 +590,9 @@ subroutine compute_albedos (aicen, vicen, & ! save ice and snow albedos (for history) albin = awtvdr*alvdrni + awtidr*alidrni & - + awtvdf*alvdfni + awtidf*alidfni + + awtvdf*alvdfni + awtidf*alidfni albsn = awtvdr*alvdrns + awtidr*alidrns & - + awtvdf*alvdfns + awtidf*alidfns + + awtvdf*alvdfns + awtidf*alidfns end subroutine compute_albedos @@ -628,8 +628,8 @@ subroutine constant_albedos (aicen, & alidrn , & ! near-ir, direct, avg (fraction) alvdfn , & ! visible, diffuse, avg (fraction) alidfn , & ! near-ir, diffuse, avg (fraction) - albin , & ! bare ice - albsn ! snow + albin , & ! bare ice + albsn ! snow ! local variables @@ -685,9 +685,9 @@ subroutine constant_albedos (aicen, & ! save ice and snow albedos (for history) albin = awtvdr*alvdrni + awtidr*alidrni & - + awtvdf*alvdfni + awtidf*alidfni + + awtvdf*alvdfni + awtidf*alidfni albsn = awtvdr*alvdrns + awtidr*alidrns & - + awtvdf*alvdfns + awtidf*alidfns + + awtvdf*alvdfns + awtidf*alidfns end subroutine constant_albedos @@ -724,7 +724,7 @@ subroutine absorbed_solar (nilyr, aicen, & logical(kind=log_kind), intent(in) :: & heat_capacity ! if true, ice has nonzero heat capacity #endif - integer (kind=int_kind), intent(in) :: & + integer (kind=int_kind), intent(in) :: & nilyr ! number of ice layers real (kind=dbl_kind), intent(in) :: & @@ -828,10 +828,10 @@ subroutine absorbed_solar (nilyr, aicen, & ! no penetrating radiation in near IR ! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir -! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir +! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir fswpen = fswpenvdr + fswpenvdf - + fswsfc = swabs - fswpen trantop = c1 ! transmittance at top of ice @@ -875,7 +875,7 @@ subroutine absorbed_solar (nilyr, aicen, & ! if zero-layer model (no heat capacity), no SW is absorbed in ice ! interior, so add to surface absorption !---------------------------------------------------------------- - + if (.not. heat_capacity) then ! SW absorbed at snow/ice surface @@ -893,7 +893,7 @@ end subroutine absorbed_solar !======================================================================= ! Begin Delta-Eddington shortwave method -! Compute initial data for Delta-Eddington method, specifically, +! Compute initial data for Delta-Eddington method, specifically, ! the approximate exponential look-up table. ! ! author: Bruce P. Briegleb, NCAR @@ -912,7 +912,7 @@ subroutine run_dEdd(dt, ncat, & #ifdef UNDEPRECATE_0LAYER heat_capacity, & #endif - tlat, tlon, & + tlat, tlon, & calendar_type, & days_per_year, & nextsw_cday, yday, & @@ -979,12 +979,12 @@ subroutine run_dEdd(dt, ncat, & kaer_3bd, & ! aerosol mass extinction cross section (m2/kg) waer_3bd, & ! aerosol single scatter albedo (fraction) gaer_3bd ! aerosol asymmetry parameter (cos(theta)) - + real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment kaer_bc_3bd, & ! aerosol mass extinction cross section (m2/kg) waer_bc_3bd, & ! aerosol single scatter albedo (fraction) gaer_bc_3bd ! aerosol asymmetry parameter (cos(theta)) - + real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment bcenh_3bd ! BC absorption enhancement factor @@ -1058,7 +1058,7 @@ subroutine run_dEdd(dt, ncat, & dhsn ! depth difference for snow on sea ice and pond ice real(kind=dbl_kind), intent(inout) :: & - coszen ! cosine solar zenith angle, < 0 for sun below horizon + coszen ! cosine solar zenith angle, < 0 for sun below horizon real(kind=dbl_kind), dimension(:), intent(inout) :: & alvdrn, & ! visible direct albedo (fraction) @@ -1067,23 +1067,23 @@ subroutine run_dEdd(dt, ncat, & alidfn, & ! near-ir diffuse albedo (fraction) fswsfcn, & ! SW absorbed at ice/snow surface (W m-2) fswintn, & ! SW absorbed in ice interior, below surface (W m-2) - fswthrun, & ! SW through ice to ocean (W/m^2) - albicen, & ! albedo bare ice - albsnon, & ! albedo snow - albpndn, & ! albedo pond + fswthrun, & ! SW through ice to ocean (W/m^2) + albicen, & ! albedo bare ice + albsnon, & ! albedo snow + albpndn, & ! albedo pond apeffn, & ! effective pond area used for radiation calculation snowfracn ! snow fraction on each category used for radiation real(kind=dbl_kind), dimension(:), intent(out), optional :: & - fswthrun_vdr, & ! vis dir SW through ice to ocean (W/m^2) - fswthrun_vdf, & ! vis dif SW through ice to ocean (W/m^2) - fswthrun_idr, & ! nir dir SW through ice to ocean (W/m^2) - fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + fswthrun_vdr, & ! vis dir SW through ice to ocean (W/m^2) + fswthrun_vdf, & ! vis dif SW through ice to ocean (W/m^2) + fswthrun_idr, & ! nir dir SW through ice to ocean (W/m^2) + fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real(kind=dbl_kind), dimension(:,:), intent(inout) :: & rsnow , & ! snow grain radius tracer (10^-6 m) Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) - Iswabsn , & ! SW radiation absorbed in ice layers (W m-2) + Iswabsn , & ! SW radiation absorbed in ice layers (W m-2) fswpenln ! visible SW entering ice layers (W m-2) logical (kind=log_kind), intent(in) :: & @@ -1200,7 +1200,7 @@ subroutine run_dEdd(dt, ncat, & if (tr_pond_cesm) then ! fraction of ice area fpn = apndn(n) - ! pond depth over fraction fpn + ! pond depth over fraction fpn hpn = hpndn(n) ! snow infiltration if (hsn >= hs_min .and. hs0 > puny) then @@ -1250,25 +1250,25 @@ subroutine run_dEdd(dt, ncat, & spn = hsnlvl - dhs ! snow depth on pond ice if (.not. linitonly .and. ipn*spn < puny) dhs = c0 dhsn(n) = dhs ! save: constant until reset to 0 - + ! not using ipn assumes that lid ice is perfectly clear ! if (ipn <= 0.3_dbl_kind) then - + ! fraction of ice area - fpn = apndn(n) * alvln(n) + fpn = apndn(n) * alvln(n) ! pond depth over fraction fpn hpn = hpndn(n) - + ! reduce effective pond area absorbing surface heat flux ! due to flux already having been used to melt pond ice fpn = (c1 - ffracn(n)) * fpn - + ! taper pond area with snow on pond ice if (dhs > puny .and. spn >= puny .and. hs1 > puny) then asnow = min(spn/hs1, c1) fpn = (c1 - asnow) * fpn endif - + ! infiltrate snow hp = hpn if (hp > puny) then @@ -1303,7 +1303,7 @@ subroutine run_dEdd(dt, ncat, & fpn = c0 endif if (apndn(n) > puny) then - hpn = hpndn(n) + hpn = hpndn(n) else fpn = c0 hpn = c0 @@ -1313,7 +1313,7 @@ subroutine run_dEdd(dt, ncat, & if (hpn < hpmin) fpn = c0 ! If ponds are present snow fraction reduced to - ! non-ponded part dEdd scheme + ! non-ponded part dEdd scheme fsn = min(fsn, c1-fpn) apeffn(n) = fpn @@ -1324,12 +1324,12 @@ subroutine run_dEdd(dt, ncat, & fsn, fpn, & hpn) if (icepack_warnings_aborted(subname)) return - + apeffn(n) = fpn ! for history fpn = c0 hpn = c0 endif ! pond type - + snowfracn(n) = fsn ! for history call shortwave_dEdd(dEdd_algae, & @@ -1408,33 +1408,33 @@ subroutine run_dEdd(dt, ncat, & deallocate(l_fswthrun_idf) end subroutine run_dEdd - + !======================================================================= ! -! Compute snow/bare ice/ponded ice shortwave albedos, absorbed and transmitted +! Compute snow/bare ice/ponded ice shortwave albedos, absorbed and transmitted ! flux using the Delta-Eddington solar radiation method as described in: ! ! A Delta-Eddington Multiple Scattering Parameterization for Solar Radiation ! in the Sea Ice Component of the Community Climate System Model ! B.P.Briegleb and B.Light NCAR/TN-472+STR February 2007 ! -! Compute shortwave albedos and fluxes for three surface types: -! snow over ice, bare ice and ponded ice. -! -! Albedos and fluxes are output for later use by thermodynamic routines. -! Invokes three calls to compute_dEdd, which sets inherent optical properties -! appropriate for the surface type. Within compute_dEdd, a call to solution_dEdd +! Compute shortwave albedos and fluxes for three surface types: +! snow over ice, bare ice and ponded ice. +! +! Albedos and fluxes are output for later use by thermodynamic routines. +! Invokes three calls to compute_dEdd, which sets inherent optical properties +! appropriate for the surface type. Within compute_dEdd, a call to solution_dEdd ! evaluates the Delta-Eddington solution. The final albedos and fluxes are then -! evaluated in compute_dEdd. Albedos and fluxes are transferred to output in +! evaluated in compute_dEdd. Albedos and fluxes are transferred to output in ! this routine. ! ! NOTE regarding albedo diagnostics: This method yields zero albedo values ! if there is no incoming solar and thus the albedo diagnostics are masked -! out when the sun is below the horizon. To estimate albedo from the history +! out when the sun is below the horizon. To estimate albedo from the history ! output (post-processing), compute ice albedo using ! (1 - albedo)*swdn = swabs. -ECH ! -! author: Bruce P. Briegleb, NCAR +! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version ! subroutine shortwave_dEdd (dEdd_algae, & @@ -1445,7 +1445,7 @@ subroutine shortwave_dEdd (dEdd_algae, & coszen, & #endif aice, vice, & - hs, fs, & + hs, fs, & rhosnw, rsnw, & fp, hp, & aero, & @@ -1497,7 +1497,7 @@ subroutine shortwave_dEdd (dEdd_algae, & #endif dEdd_algae, & ! .true. use prognostic chla in dEdd modal_aero ! .true. use modal aerosol treatment - + real (kind=dbl_kind), dimension(:,:), intent(in) :: & kaer_3bd, & ! aerosol mass extinction cross section (m2/kg) waer_3bd, & ! aerosol single scatter albedo (fraction) @@ -1507,7 +1507,7 @@ subroutine shortwave_dEdd (dEdd_algae, & kaer_bc_3bd, & ! aerosol mass extinction cross section (m2/kg) waer_bc_3bd, & ! aerosol single scatter albedo (fraction) gaer_bc_3bd ! aerosol asymmetry parameter (cos(theta)) - + real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment bcenh_3bd ! BC absorption enhancement factor @@ -1541,8 +1541,8 @@ subroutine shortwave_dEdd (dEdd_algae, & kalg , & ! algae absorption coefficient R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo - aice , & ! concentration of ice - vice , & ! volume of ice + aice , & ! concentration of ice + vice , & ! volume of ice hs , & ! snow depth fs ! horizontal coverage of snow @@ -1553,19 +1553,19 @@ subroutine shortwave_dEdd (dEdd_algae, & zbio ! shortwave tracers (zaero+chla) real (kind=dbl_kind), intent(in) :: & - fp , & ! pond fractional coverage (0 to 1) - hp , & ! pond depth (m) + fp , & ! pond fractional coverage (0 to 1) + hp , & ! pond depth (m) swvdr , & ! sw down, visible, direct (W/m^2) swvdf , & ! sw down, visible, diffuse (W/m^2) swidr , & ! sw down, near IR, direct (W/m^2) swidf ! sw down, near IR, diffuse (W/m^2) real (kind=dbl_kind), intent(inout) :: & - coszen , & ! cosine of solar zenith angle - alvdr , & ! visible, direct, albedo (fraction) - alvdf , & ! visible, diffuse, albedo (fraction) - alidr , & ! near-ir, direct, albedo (fraction) - alidf , & ! near-ir, diffuse, albedo (fraction) + coszen , & ! cosine of solar zenith angle + alvdr , & ! visible, direct, albedo (fraction) + alvdf , & ! visible, diffuse, albedo (fraction) + alidr , & ! near-ir, direct, albedo (fraction) + alidf , & ! near-ir, diffuse, albedo (fraction) fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2) fswint , & ! SW interior absorption (below surface, above ocean,W m-2) fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2) @@ -1575,16 +1575,16 @@ subroutine shortwave_dEdd (dEdd_algae, & fswthru_vdf , & ! vis dif SW through snow/bare ice/ponded ice into ocean (W m-2) fswthru_idr , & ! nir dir SW through snow/bare ice/ponded ice into ocean (W m-2) fswthru_idf ! nir dif SW through snow/bare ice/ponded ice into ocean (W m-2) - + real (kind=dbl_kind), dimension (:), intent(inout) :: & fswpenl , & ! visible SW entering ice layers (W m-2) Sswabs , & ! SW absorbed in snow layer (W m-2) Iswabs ! SW absorbed in ice layer (W m-2) real (kind=dbl_kind), intent(out) :: & - albice , & ! bare ice albedo, for history - albsno , & ! snow albedo, for history - albpnd ! pond albedo, for history + albice , & ! bare ice albedo, for history + albsno , & ! snow albedo, for history + albpnd ! pond albedo, for history logical (kind=log_kind) , intent(in) :: & l_print_point @@ -1604,7 +1604,7 @@ subroutine shortwave_dEdd (dEdd_algae, & integer (kind=int_kind) :: & nspint , & ! number of solar spectral bands srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) - + integer (kind=int_kind) :: & k , & ! level index na , & ! aerosol index @@ -1613,7 +1613,7 @@ subroutine shortwave_dEdd (dEdd_algae, & ! (0 layer is included also) real (kind=dbl_kind) :: & - vsno ! volume of snow + vsno ! volume of snow real (kind=dbl_kind) :: & swdn , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j) @@ -1622,10 +1622,10 @@ subroutine shortwave_dEdd (dEdd_algae, & ! for history real (kind=dbl_kind) :: & - avdrl , & ! visible, direct, albedo (fraction) - avdfl , & ! visible, diffuse, albedo (fraction) - aidrl , & ! near-ir, direct, albedo (fraction) - aidfl ! near-ir, diffuse, albedo (fraction) + avdrl , & ! visible, direct, albedo (fraction) + avdfl , & ! visible, diffuse, albedo (fraction) + aidrl , & ! near-ir, direct, albedo (fraction) + aidfl ! near-ir, diffuse, albedo (fraction) character(len=*),parameter :: subname='(shortwave_dEdd)' @@ -1667,7 +1667,7 @@ subroutine shortwave_dEdd (dEdd_algae, & Iswabs(:) = c0 ! compute aerosol mass path - + aero_mp(:) = c0 if( tr_aero ) then ! check 4 layers for each aerosol, a snow SSL, snow below SSL, @@ -1689,9 +1689,9 @@ subroutine shortwave_dEdd (dEdd_algae, & enddo ! na endif ! if aerosols - ! compute shortwave radiation accounting for snow/ice (both snow over + ! compute shortwave radiation accounting for snow/ice (both snow over ! ice and bare ice) and ponded ice (if any): - + ! sea ice points with sun above horizon netsw = swvdr + swidr + swvdf + swidf if (netsw > puny) then ! sun above horizon @@ -1730,7 +1730,7 @@ subroutine shortwave_dEdd (dEdd_algae, & Iswabs, fswpenl, & l_use_snicar) if (icepack_warnings_aborted(subname)) return - + alvdr = alvdr + avdrl *fi alvdf = alvdf + avdfl *fi alidr = alidr + aidrl *fi @@ -1738,10 +1738,10 @@ subroutine shortwave_dEdd (dEdd_algae, & ! for history albice = albice & + awtvdr*avdrl + awtidr*aidrl & - + awtvdf*avdfl + awtidf*aidfl + + awtvdf*avdfl + awtidf*aidfl endif endif - + ! sea ice points with sun above horizon netsw = swvdr + swidr + swvdf + swidf if (netsw > puny) then ! sun above horizon @@ -1808,7 +1808,7 @@ subroutine shortwave_dEdd (dEdd_algae, & l_use_snicar) endif if (icepack_warnings_aborted(subname)) return - + alvdr = alvdr + avdrl *fs alvdf = alvdf + avdfl *fs alidr = alidr + aidrl *fs @@ -1816,7 +1816,7 @@ subroutine shortwave_dEdd (dEdd_algae, & ! for history albsno = albsno & + awtvdr*avdrl + awtidr*aidrl & - + awtvdf*avdfl + awtidf*aidfl + + awtvdf*avdfl + awtidf*aidfl endif endif @@ -1830,7 +1830,7 @@ subroutine shortwave_dEdd (dEdd_algae, & ! if nonzero pond fraction and sufficient pond depth ! if( fp > puny .and. hp > hpmin ) then if (fp > puny) then - + ! calculate ponded ice srftyp = 2 @@ -1846,7 +1846,7 @@ subroutine shortwave_dEdd (dEdd_algae, & kaer_3bd, waer_3bd, gaer_3bd, & kaer_bc_3bd, waer_bc_3bd, gaer_bc_3bd, & bcenh_3bd, modal_aero, kalg, & - swvdr, swvdf, swidr, swidf, srftyp, & + swvdr, swvdf, swidr, swidf, srftyp, & hs, rhosnw, rsnw, hi, hp, & fp, aero_mp, avdrl, avdfl, & aidrl, aidfl, & @@ -1860,7 +1860,7 @@ subroutine shortwave_dEdd (dEdd_algae, & Iswabs, fswpenl, & l_use_snicar) if (icepack_warnings_aborted(subname)) return - + alvdr = alvdr + avdrl *fp alvdf = alvdf + avdfl *fp alidr = alidr + aidrl *fp @@ -1868,7 +1868,7 @@ subroutine shortwave_dEdd (dEdd_algae, & ! for history albpnd = albpnd & + awtvdr*avdrl + awtidr*aidrl & - + awtvdf*avdfl + awtidf*aidfl + + awtvdf*avdfl + awtidf*aidfl endif endif @@ -1928,7 +1928,7 @@ subroutine shortwave_dEdd (dEdd_algae, & swab = fswsfc+fswint+fswthru swalb = (1.-swab/(swdn+.0001)) write(warnstr,*) subname, ' swdn swab swalb = ',swdn,swab,swalb - do k = 1, nslyr + do k = 1, nslyr write(warnstr,*) subname, ' snow layer k = ', k, & ' rhosnw = ', & rhosnw(k), & @@ -1936,12 +1936,12 @@ subroutine shortwave_dEdd (dEdd_algae, & rsnw(k) call icepack_warnings_add(warnstr) enddo - do k = 1, nslyr + do k = 1, nslyr write(warnstr,*) subname, ' snow layer k = ', k, & ' Sswabs(k) = ', Sswabs(k) call icepack_warnings_add(warnstr) enddo - do k = 1, nilyr + do k = 1, nilyr write(warnstr,*) subname, ' sea ice layer k = ', k, & ' Iswabs(k) = ', Iswabs(k) call icepack_warnings_add(warnstr) @@ -1953,10 +1953,10 @@ end subroutine shortwave_dEdd !======================================================================= ! -! Evaluate snow/ice/ponded ice inherent optical properties (IOPs), and +! Evaluate snow/ice/ponded ice inherent optical properties (IOPs), and ! then calculate the multiple scattering solution by calling solution_dEdd. ! -! author: Bruce P. Briegleb, NCAR +! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version ! 2018: Cheng Dang merged with SNICAR 5-band snow and aersols IOPs, UC Irvine @@ -2009,7 +2009,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & klev , & ! number of radiation layers - 1 klevp ! number of radiation interfaces - 1 ! (0 layer is included also) - + logical (kind=log_kind), intent(in) :: & #ifdef UNDEPRECATE_0LAYER heat_capacity,& ! if true, ice has nonzero heat capacity @@ -2021,7 +2021,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg) waer_bc_tab, & ! aerosol single scatter albedo (fraction) gaer_bc_tab ! aerosol asymmetry parameter (cos(theta)) - + real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment bcenh ! BC absorption enhancement factor @@ -2030,11 +2030,11 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo R_pnd ! ponded ice tuning parameter; +1 > 1sig increase in albedo - real (kind=dbl_kind), dimension(:,:), intent(in) :: & + real (kind=dbl_kind), dimension(:,:), intent(in) :: & kaer_tab, & ! aerosol mass extinction cross section (m2/kg) waer_tab, & ! aerosol single scatter albedo (fraction) gaer_tab ! aerosol asymmetry parameter (cos(theta)) - + real (kind=dbl_kind), intent(in) :: & kalg , & ! algae absorption coefficient fnidr , & ! fraction of direct to total down flux in nir @@ -2043,11 +2043,11 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & swvdf , & ! shortwave down at surface, visible, diffuse (W/m^2) swidr , & ! shortwave down at surface, near IR, direct (W/m^2) swidf ! shortwave down at surface, near IR, diffuse (W/m^2) - + integer (kind=int_kind), intent(in) :: & nspint , & ! number of solar spectral bands srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) - + real (kind=dbl_kind), intent(in) :: & hs ! snow thickness (m) @@ -2061,12 +2061,12 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & hi , & ! ice thickness (m) hp , & ! pond depth (m) fi ! snow/bare ice fractional coverage (0 to 1) - + real (kind=dbl_kind), intent(inout) :: & - alvdr , & ! visible, direct, albedo (fraction) - alvdf , & ! visible, diffuse, albedo (fraction) - alidr , & ! near-ir, direct, albedo (fraction) - alidf , & ! near-ir, diffuse, albedo (fraction) + alvdr , & ! visible, direct, albedo (fraction) + alvdf , & ! visible, diffuse, albedo (fraction) + alidr , & ! near-ir, direct, albedo (fraction) + alidf , & ! near-ir, diffuse, albedo (fraction) fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2) fswint , & ! SW interior absorption (below surface, above ocean,W m-2) fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2) @@ -2076,7 +2076,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & fswthru_vdf , & ! vis dif SW through snow/bare ice/ponded ice into ocean (W m-2) fswthru_idr , & ! nir dir SW through snow/bare ice/ponded ice into ocean (W m-2) fswthru_idf ! nir dif SW through snow/bare ice/ponded ice into ocean (W m-2) - + real (kind=dbl_kind), dimension (:), intent(inout) :: & fswpenl , & ! visible SW entering ice layers (W m-2) Sswabs , & ! SW absorbed in snow layer (W m-2) @@ -2097,12 +2097,12 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & !----------------------------------------------------------------------- ! -! Set up optical property profiles, based on snow, sea ice and ponded +! Set up optical property profiles, based on snow, sea ice and ponded ! ice IOPs from: ! -! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple -! Scattering Parameterization for Solar Radiation in the Sea Ice -! Component of the Community Climate System Model, NCAR Technical +! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple +! Scattering Parameterization for Solar Radiation in the Sea Ice +! Component of the Community Climate System Model, NCAR Technical ! Note NCAR/TN-472+STR February 2007 ! ! Computes column Delta-Eddington radiation solution for specific @@ -2111,7 +2111,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! Divides solar spectrum into 3 intervals: 0.2-0.7, 0.7-1.19, and ! 1.19-5.0 micro-meters. The latter two are added (using an assumed ! partition of incident shortwave in the 0.7-5.0 micro-meter band between -! the 0.7-1.19 and 1.19-5.0 micro-meter band) to give the final output +! the 0.7-1.19 and 1.19-5.0 micro-meter band) to give the final output ! of 0.2-0.7 visible and 0.7-5.0 near-infrared albedos and fluxes. ! ! Specifies vertical layer optical properties based on input snow depth, @@ -2126,14 +2126,14 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! Please read the following; otherwise, there is 99.9% chance you ! will be confused about indices at some point in time........ :) ! -! CICE4.0 snow treatment has one snow layer above the sea ice. This +! CICE4.0 snow treatment has one snow layer above the sea ice. This ! snow layer has finite heat capacity, so that surface absorption must ! be distinguished from internal. The Delta-Eddington solar radiation ! thus adds extra surface scattering layers to both snow and sea ice. ! Note that in the following, we assume a fixed vertical layer structure -! for the radiation calculation. In other words, we always have the -! structure shown below for one snow and four sea ice layers, but for -! ponded ice the pond fills "snow" layer 1 over the sea ice, and for +! for the radiation calculation. In other words, we always have the +! structure shown below for one snow and four sea ice layers, but for +! ponded ice the pond fills "snow" layer 1 over the sea ice, and for ! bare sea ice the top layers over sea ice are treated as transparent air. ! ! SSL = surface scattering layer for either snow or sea ice @@ -2187,7 +2187,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! for surface heating, and that absorbed in the sea ice DL is ! used for sea ice layer 1 heating. ! -! Basically, vertical profiles of the layer extinction optical depth (tau), +! Basically, vertical profiles of the layer extinction optical depth (tau), ! single scattering albedo (w0) and asymmetry parameter (g) are required over ! the klev+1 layers, where klev+1 = 2 + nslyr + nilyr. All of the surface type ! information and snow/ice iop properties are evaulated in this routine, so @@ -2210,16 +2210,16 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ksnow , & ! level index for snow density and grain size kii ! level starting index for sea ice (nslyr+1) - integer (kind=int_kind), parameter :: & + integer (kind=int_kind), parameter :: & nmbrad = 32 ! number of snow grain radii in tables - - real (kind=dbl_kind) :: & + + real (kind=dbl_kind) :: & avdr , & ! visible albedo, direct (fraction) avdf , & ! visible albedo, diffuse (fraction) aidr , & ! near-ir albedo, direct (fraction) aidf ! near-ir albedo, diffuse (fraction) - - real (kind=dbl_kind) :: & + + real (kind=dbl_kind) :: & fsfc , & ! shortwave absorbed at snow/bare ice/ponded ice surface (W m-2) fint , & ! shortwave absorbed in interior (W m-2) fthru , & ! shortwave through snow/bare ice/ponded ice to ocean (W/m^2) @@ -2228,28 +2228,28 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & fthruidr, & ! nir dir shortwave through snow/bare ice/ponded ice to ocean (W/m^2) fthruidf ! nir dif shortwave through snow/bare ice/ponded ice to ocean (W/m^2) - real (kind=dbl_kind), dimension(nslyr) :: & + real (kind=dbl_kind), dimension(nslyr) :: & Sabs ! shortwave absorbed in snow layer (W m-2) - real (kind=dbl_kind), dimension(nilyr) :: & + real (kind=dbl_kind), dimension(nilyr) :: & Iabs ! shortwave absorbed in ice layer (W m-2) - - real (kind=dbl_kind), dimension(nilyr+1) :: & + + real (kind=dbl_kind), dimension(nilyr+1) :: & fthrul ! shortwave through to ice layers (W m-2) real (kind=dbl_kind), dimension (nspint_3bd) :: & wghtns ! spectral weights - - real (kind=dbl_kind), parameter :: & + + real (kind=dbl_kind), parameter :: & cp67 = 0.67_dbl_kind , & ! nir band weight parameter cp78 = 0.78_dbl_kind , & ! nir band weight parameter cp01 = 0.01_dbl_kind ! for ocean visible albedo - + real (kind=dbl_kind), dimension (0:klev) :: & tau , & ! layer extinction optical depth w0 , & ! layer single scattering albedo g ! layer asymmetry parameter - + ! following arrays are defined at model interfaces; 0 is the top of the ! layer above the sea ice; klevp is the sea ice/ocean interface. real (kind=dbl_kind), dimension (0:klevp) :: & @@ -2259,7 +2259,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & rupdir , & ! reflectivity to direct radiation for layers below rupdif , & ! reflectivity to diffuse radiation for layers below rdndif ! reflectivity to diffuse radiation for layers above - + real (kind=dbl_kind), dimension (0:klevp) :: & dfdir , & ! down-up flux at interface due to direct beam at top surface dfdif ! down-up flux at interface due to diffuse beam at top surface @@ -2273,7 +2273,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ws , & ! Snow single scattering albedo gs ! Snow asymmetry parameter - real (kind=dbl_kind), dimension(nslyr) :: & + real (kind=dbl_kind), dimension(nslyr) :: & frsnw ! snow grain radius in snow layer * adjustment factor (m) ! ice and ponded ice IOPs, allowing for tuning @@ -2303,7 +2303,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & dz_ssl , & ! snow or sea ice surface scattering layer thickness fs ! scaling factor to reduce (nilyr<4) or increase (nilyr>4) DL ! extinction coefficient to maintain DL optical depth constant - ! with changing number of sea ice layers, to approximately + ! with changing number of sea ice layers, to approximately ! conserve computed albedo for constant physical depth of sea ! ice when the number of sea ice layers vary real (kind=dbl_kind) :: & @@ -2320,8 +2320,8 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & real (kind=dbl_kind) :: & albodr , & ! spectral ocean albedo to direct rad albodf ! spectral ocean albedo to diffuse rad - - ! for melt pond transition to bare sea ice for small pond depths + + ! for melt pond transition to bare sea ice for small pond depths real (kind=dbl_kind) :: & sig_i , & ! ice scattering coefficient (/m) sig_p , & ! pond scattering coefficient (/m) @@ -2631,7 +2631,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & kii = nslyr + 1 ! initialize albedos and fluxes to 0 - fthrul = c0 + fthrul = c0 Iabs = c0 kabs_chl(:,:) = c0 tzaer(:,:) = c0 @@ -2649,7 +2649,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & fthruvdf = c0 fthruidr = c0 fthruidf = c0 - + if (nspint == nspint_3bd) then ki_ssl_mn = ki_ssl_mn_3bd wi_ssl_mn = wi_ssl_mn_3bd @@ -2676,10 +2676,10 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! spectral weights if (nspint == nspint_3bd) then - ! weights 2 (0.7-1.19 micro-meters) and 3 (1.19-5.0 micro-meters) - ! are chosen based on 1D calculations using ratio of direct to total - ! near-infrared solar (0.7-5.0 micro-meter) which indicates clear/cloudy - ! conditions: more cloud, the less 1.19-5.0 relative to the + ! weights 2 (0.7-1.19 micro-meters) and 3 (1.19-5.0 micro-meters) + ! are chosen based on 1D calculations using ratio of direct to total + ! near-infrared solar (0.7-5.0 micro-meter) which indicates clear/cloudy + ! conditions: more cloud, the less 1.19-5.0 relative to the ! 0.7-1.19 micro-meter due to cloud absorption. wghtns(1) = c1 wghtns(2) = cp67 + (cp78-cp67)*(c1-fnidr) @@ -2723,7 +2723,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & dzk(k) = dz enddo endif - + ! ice dz = hi*rnilyr ! empirical reduction in sea ice ssl thickness for ice thinner than 1.5m; @@ -2735,7 +2735,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! set sea ice ssl thickness to half top layer if sea ice thin enough !ech: note this is highly resolution dependent! dz_ssl = min(dz_ssl, dz/c2) - + dzk(kii) = dz_ssl dzk(kii+1) = dz - dz_ssl if (kii+2 <= klev) then @@ -2840,12 +2840,12 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & enddo else k = klev - kabs_chl(1,k) = kalg*(0.50_dbl_kind/dzk(k)) + kabs_chl(1,k) = kalg*(0.50_dbl_kind/dzk(k)) endif ! kabs_chl ! aerosols if (modal_aero) then - do k=0,klev + do k=0,klev if (k < nslyr+1) then ! define indices for snow layer ! use top rsnw, rhosnw for snow ssl and rest of top layer ! Cheng: note that aerosol IOPs are related to snow grain radius. @@ -2857,7 +2857,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & elseif (l_use_snicar .and. nspint == nspint_5bd) then tmp_gs = rsnw(ksnow) endif - + ! grain size index: works for 25 < snw_rds < 1625 um: if (tmp_gs < 125._dbl_kind) then tmp1 = tmp_gs/50._dbl_kind @@ -2890,8 +2890,8 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! call icepack_warnings_add(warnstr) enddo ! k - if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll - do n = 1,n_zaero + if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll + do n = 1,n_zaero if (n == 1) then ! interstitial BC do k = 0, klev do ns = 1,nspint ! not weighted by aice @@ -2907,7 +2907,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & enddo elseif (n==2) then ! within-ice BC do k = 0, klev - do ns = 1,nspint + do ns = 1,nspint tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k)) * & bcenh(ns,k_bcins(k),k_bcini(k))* & zbio(nlt_zaero_sw(n)+k)*dzk(k) @@ -2948,20 +2948,20 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & enddo ! nspint enddo enddo - endif !tr_zaero + endif !tr_zaero endif ! modal_aero !----------------------------------------------------------------------- - + ! begin spectral loop !echmod - split this loop for efficiency, if possible (move conditionals outside of the loop) do ns = 1, nspint - + ! set optical properties of air/snow/pond overlying sea ice ! air if( srftyp == 0 ) then - do k=0,nslyr + do k=0,nslyr tau(k) = c0 w0(k) = c0 g(k) = c0 @@ -3178,9 +3178,9 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & gaer = gaer + & (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* & waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k)) - + else ! other species (dust) - + taer = taer + & (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4)) waer = waer + & @@ -3222,9 +3222,9 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! no aerosol in pond enddo ! k endif ! srftyp - + ! set optical properties of sea ice - + ! bare or snow-covered sea ice layers if( srftyp <= 1 ) then ! ssl @@ -3254,7 +3254,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & if( ns == 1 ) then ! total layer absorption optical depth fixed at value ! of kalg*0.50m, independent of actual layer thickness - kabs = kabs + kabs_chl(ns,k) + kabs = kabs + kabs_chl(ns,k) endif sig = ki_int(ns)*wi_int(ns) tau(k) = (kabs+sig)*dzk(k) @@ -3262,7 +3262,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & g(k) = gi_int(ns) ! aerosol in sea ice if (tr_zaero .and. dEdd_algae) then - do k = kii, klev + do k = kii, klev gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny) wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny) g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) & @@ -3271,7 +3271,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & / (tau(k) + tzaer(ns,k)) tau(k) = tau(k) + tzaer(ns,k) enddo - elseif (tr_aero) then + elseif (tr_aero) then k = kii ! sea ice SSL taer = c0 waer = c0 @@ -3296,7 +3296,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & waer_bc_tab(ns,k_bcins(k)) gaer = gaer + & aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* & - waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k)) + waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k)) else ! other species (dust) taer = taer + & aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4)) @@ -3431,40 +3431,40 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & enddo ! k endif endif ! small pond depth transition to bare sea ice - endif ! srftyp - + endif ! srftyp + ! set reflectivities for ocean underlying sea ice rns = real(ns-1, kind=dbl_kind) albodr = cp01 * (c1 - min(rns, c1)) albodf = cp01 * (c1 - min(rns, c1)) - + ! layer input properties now completely specified: tau, w0, g, - ! albodr, albodf; now compute the Delta-Eddington solution + ! albodr, albodf; now compute the Delta-Eddington solution ! reflectivities and transmissivities for each layer; then, ! combine the layers going downwards accounting for multiple - ! scattering between layers, and finally start from the + ! scattering between layers, and finally start from the ! underlying ocean and combine successive layers upwards to ! the surface; see comments in solution_dEdd for more details. - + call solution_dEdd & (coszen, srftyp, klev, klevp, nslyr, & tau, w0, g, albodr, albodf, & trndir, trntdr, trndif, rupdir, rupdif, & - rdndif) + rdndif) if (icepack_warnings_aborted(subname)) return ! the interface reflectivities and transmissivities required ! to evaluate interface fluxes are returned from solution_dEdd; - ! now compute up and down fluxes for each interface, using the + ! now compute up and down fluxes for each interface, using the ! combined layer properties at each interface: ! ! layers interface ! ! --------------------- k ! k - ! --------------------- - - do k = 0, klevp + ! --------------------- + + do k = 0, klevp ! interface scattering refk = c1/(c1 - rdndif(k)*rupdif(k)) ! dir tran ref from below times interface scattering, plus diff @@ -3473,7 +3473,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! (trntdr(k)-trndir(k)) & ! *rupdif(k))*refk ! dir tran plus total diff trans times interface scattering plus - ! dir tran with up dir ref and down dif ref times interface scattering + ! dir tran with up dir ref and down dif ref times interface scattering ! fdirdn(k) = trndir(k) + (trntdr(k) & ! - trndir(k) + trndir(k) & ! *rupdir(k)*rdndif(k))*refk @@ -3490,8 +3490,8 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! dfdif = fdifdn - fdifup dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary? - enddo ! k - + enddo ! k + !echmod - is this necessary? ! ! note that because the snow IOPs for diffuse and direct incidents ! ! are different, the snow albedo needs to be calculated twice for @@ -3511,14 +3511,14 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ! calculate final surface albedos and fluxes- ! all absorbed flux above ksrf is included in surface absorption - + if( ns == 1) then ! visible - + swdr = swvdr swdf = swvdf avdr = rupdir(0) avdf = rupdif(0) - + tmp_0 = dfdir(0 )*swdr + dfdif(0 )*swdf tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf @@ -3527,13 +3527,13 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering fthrul(k-nslyr-1) = dfdir(k)*swdr + dfdif(k)*swdf enddo - + fsfc = fsfc + tmp_0 - tmp_ks fint = fint + tmp_ks - tmp_kl fthru = fthru + tmp_kl fthruvdr = fthruvdr + dfdir(klevp)*swdr fthruvdf = fthruvdf + dfdif(klevp)*swdf - + ! if snow covered ice, set snow internal absorption; else, Sabs=0 if( srftyp == 1 ) then ki = 0 @@ -3553,7 +3553,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & ki = 0 do k=nslyr+2,nslyr+1+nilyr ! for bare ice, DL absorption for sea ice layer 1 - km = k + km = k kp = km + 1 ! modify for top sea ice layer for snow over sea ice if( srftyp == 1 ) then @@ -3568,9 +3568,9 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & + dfdir(km)*swdr + dfdif(km)*swdf & - (dfdir(kp)*swdr + dfdif(kp)*swdf) enddo ! k - + else ! ns > 1, near IR - + swdr = swidr swdf = swidf @@ -3681,7 +3681,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & + (dfdif(km)-dfdir(kp))*swdf*wghtns_5bd_dfs(ns) enddo ! k endif - + ! complex indexing to insure proper absorptions for sea ice ki = 0 do k=nslyr+2,nslyr+1+nilyr @@ -3705,7 +3705,7 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & endif ! nspint endif ! ns = 1, ns > 1 - + enddo ! end spectral loop ns ! solar zenith angle parameterization @@ -3746,29 +3746,29 @@ subroutine compute_dEdd(nilyr, nslyr, nspint, & do k = 1, nslyr Sswabs(k) = Sswabs(k) + Sabs(k)*fi enddo ! k - + do k = 1, nilyr Iswabs(k) = Iswabs(k) + Iabs(k)*fi - - ! bgc layer + + ! bgc layer fswpenl(k) = fswpenl(k) + fthrul(k)* fi if (k == nilyr) then fswpenl(k+1) = fswpenl(k+1) + fthrul(k+1)*fi endif enddo ! k - + #ifdef UNDEPRECATE_0LAYER !---------------------------------------------------------------- - ! if ice has zero heat capacity, no SW can be absorbed + ! if ice has zero heat capacity, no SW can be absorbed ! in the ice/snow interior, so add to surface absorption. ! Note: nilyr = nslyr = 1 for this case !---------------------------------------------------------------- if (.not. heat_capacity) then - + ! SW absorbed at snow/ice surface fswsfc = fswsfc + Iswabs(1) + Sswabs(1) - + ! SW absorbed in ice interior fswint = c0 Iswabs(1) = c0 @@ -3780,7 +3780,7 @@ end subroutine compute_dEdd !======================================================================= ! -! Given input vertical profiles of optical properties, evaluate the +! Given input vertical profiles of optical properties, evaluate the ! monochromatic Delta-Eddington solution. ! ! author: Bruce P. Briegleb, NCAR @@ -3800,16 +3800,16 @@ subroutine solution_dEdd & klevp , & ! number of radiation interfaces - 1 ! (0 layer is included also) nslyr ! number of snow layers - + real (kind=dbl_kind), dimension(0:klev), intent(in) :: & tau , & ! layer extinction optical depth w0 , & ! layer single scattering albedo g ! layer asymmetry parameter - + real (kind=dbl_kind), intent(in) :: & albodr , & ! ocean albedo to direct rad albodf ! ocean albedo to diffuse rad - + ! following arrays are defined at model interfaces; 0 is the top of the ! layer above the sea ice; klevp is the sea ice/ocean interface. real (kind=dbl_kind), dimension (0:klevp), intent(out) :: & @@ -3841,7 +3841,7 @@ subroutine solution_dEdd & ! Assumes monochromatic (spectrally uniform) properties across a band ! for the input optical parameters. ! -! If total transmission of the direct beam to the interface above a particular +! If total transmission of the direct beam to the interface above a particular ! layer is less than trmin, then no further Delta-Eddington solutions are ! evaluated for layers below. ! @@ -3849,16 +3849,16 @@ subroutine solution_dEdd & ! ! First, we assume that radiation is refracted when entering either ! sea ice at the base of the surface scattering layer, or water (i.e. melt -! pond); we assume that radiation does not refract when entering snow, nor -! upon entering sea ice from a melt pond, nor upon entering the underlying +! pond); we assume that radiation does not refract when entering snow, nor +! upon entering sea ice from a melt pond, nor upon entering the underlying ! ocean from sea ice. ! ! To handle refraction, we define a "fresnel" layer, which physically -! is of neglible thickness and is non-absorbing, which can be combined to -! any sea ice layer or top of melt pond. The fresnel layer accounts for +! is of neglible thickness and is non-absorbing, which can be combined to +! any sea ice layer or top of melt pond. The fresnel layer accounts for ! refraction of direct beam and associated reflection and transmission for -! solar radiation. A fresnel layer is combined with the top of a melt pond -! or to the surface scattering layer of sea ice if no melt pond lies over it. +! solar radiation. A fresnel layer is combined with the top of a melt pond +! or to the surface scattering layer of sea ice if no melt pond lies over it. ! ! Some caution must be exercised for the fresnel layer, because any layer ! to which it is combined is no longer a homogeneous layer, as are all other @@ -3875,11 +3875,11 @@ subroutine solution_dEdd & integer (kind=int_kind) :: & kfrsnl ! radiation interface index for fresnel layer - + ! following variables are defined for each layer; 0 refers to the top - ! layer. In general we must distinguish directions above and below in + ! layer. In general we must distinguish directions above and below in ! the diffuse reflectivity and transmissivity, as layers are not assumed - ! to be homogeneous (apart from the single layer Delta-Edd solutions); + ! to be homogeneous (apart from the single layer Delta-Edd solutions); ! the direct is always from above. real (kind=dbl_kind), dimension (0:klev) :: & rdir , & ! layer reflectivity to direct radiation @@ -3890,14 +3890,14 @@ subroutine solution_dEdd & tdif_b , & ! layer transmission to diffuse radiation from below trnlay ! solar beam transm for layer (direct beam only) - integer (kind=int_kind) :: & + integer (kind=int_kind) :: & k ! level index - + real (kind=dbl_kind), parameter :: & trmin = 0.001_dbl_kind ! minimum total transmission allowed ! total transmission is that due to the direct beam; i.e. it includes ! both the directly transmitted solar beam and the diffuse downwards - ! transmitted radiation resulting from scattering out of the direct beam + ! transmitted radiation resulting from scattering out of the direct beam real (kind=dbl_kind) :: & tautot , & ! layer optical depth wtot , & ! layer single scattering albedo @@ -3909,9 +3909,9 @@ subroutine solution_dEdd & rintfc , & ! reflection (multiple) at an interface refkp1 , & ! interface multiple scattering for k+1 refkm1 , & ! interface multiple scattering for k-1 - tdrrdir , & ! direct tran times layer direct ref + tdrrdir , & ! direct tran times layer direct ref tdndif ! total down diffuse = tot tran - direct tran - + ! perpendicular and parallel relative to plane of incidence and scattering real (kind=dbl_kind) :: & R1 , & ! perpendicular polarization reflection amplitude @@ -3924,21 +3924,21 @@ subroutine solution_dEdd & Rf_dif_b , & ! fresnel reflection to diff radiation from below Tf_dif_a , & ! fresnel transmission to diff radiation from above Tf_dif_b ! fresnel transmission to diff radiation from below - + ! refractive index for sea ice, water; pre-computed, band-independent, ! diffuse fresnel reflectivities - real (kind=dbl_kind), parameter :: & + real (kind=dbl_kind), parameter :: & refindx = 1.310_dbl_kind , & ! refractive index of sea ice (water also) cp063 = 0.063_dbl_kind , & ! diffuse fresnel reflectivity from above cp455 = 0.455_dbl_kind ! diffuse fresnel reflectivity from below - + real (kind=dbl_kind) :: & mu0 , & ! cosine solar zenith angle incident mu0nij ! cosine solar zenith angle in medium below fresnel level - + real (kind=dbl_kind) :: & mu0n ! cosine solar zenith angle in medium - + real (kind=dbl_kind) :: & alp , & ! temporary for alpha gam , & ! temporary for agamm @@ -3984,7 +3984,7 @@ subroutine solution_dEdd & !----------------------------------------------------------------------- - do k = 0, klevp + do k = 0, klevp trndir(k) = c0 trntdr(k) = c0 trndif(k) = c0 @@ -3993,20 +3993,20 @@ subroutine solution_dEdd & rdndif(k) = c0 enddo - ! initialize top interface of top layer + ! initialize top interface of top layer trndir(0) = c1 trntdr(0) = c1 trndif(0) = c1 rdndif(0) = c0 - ! mu0 is cosine solar zenith angle above the fresnel level; make + ! mu0 is cosine solar zenith angle above the fresnel level; make ! sure mu0 is large enough for stable and meaningful radiation ! solution: .01 is like sun just touching horizon with its lower edge mu0 = max(coszen,p01) ! mu0n is cosine solar zenith angle used to compute the layer ! Delta-Eddington solution; it is initially computed to be the - ! value below the fresnel level, i.e. the cosine solar zenith + ! value below the fresnel level, i.e. the cosine solar zenith ! angle below the fresnel level for the refracted solar beam: mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx))) @@ -4017,7 +4017,7 @@ subroutine solution_dEdd & ! at base of sea ice SSL (and top of the sea ice DL); the ! snow SSL counts for one, then the number of snow layers, ! then the sea ice SSL which also counts for one: - if( srftyp < 2 ) kfrsnl = nslyr + 2 + if( srftyp < 2 ) kfrsnl = nslyr + 2 ! proceed down one layer at a time; if the total transmission to ! the interface just above a given layer is less than trmin, then no @@ -4037,7 +4037,7 @@ subroutine solution_dEdd & ! compute next layer Delta-eddington solution only if total transmission ! of radiation to the interface just above the layer exceeds trmin. - + if (trntdr(k) > trmin ) then ! calculation over layers with penetrating radiation @@ -4076,7 +4076,7 @@ subroutine solution_dEdd & amg = alp - gam rdir(k) = apg*rdif_a(k) + amg*(tdif_a(k)*trnlay(k) - c1) tdir(k) = apg*tdif_a(k) + (amg* rdif_a(k)-apg+c1)*trnlay(k) - + ! recalculate rdif,tdif using direct angular integration over rdir,tdir, ! since Delta-Eddington rdif formula is not well-behaved (it is usually ! biased low and can even be negative); use ngmax angles and gaussian @@ -4103,20 +4103,20 @@ subroutine solution_dEdd & enddo ! ng rdif_a(k) = smr/swt tdif_a(k) = smt/swt - + ! homogeneous layer rdif_b(k) = rdif_a(k) tdif_b(k) = tdif_a(k) - ! add fresnel layer to top of desired layer if either - ! air or snow overlies ice; we ignore refraction in ice + ! add fresnel layer to top of desired layer if either + ! air or snow overlies ice; we ignore refraction in ice ! if a melt pond overlies it: if( k == kfrsnl ) then ! compute fresnel reflection and transmission amplitudes ! for two polarizations: 1=perpendicular and 2=parallel to ! the plane containing incident, reflected and refracted rays. - R1 = (mu0 - refindx*mu0n) / & + R1 = (mu0 - refindx*mu0n) / & (mu0 + refindx*mu0n) R2 = (refindx*mu0 - mu0n) / & (refindx*mu0 + mu0n) @@ -4142,7 +4142,7 @@ subroutine solution_dEdd & Rf_dif_b = cp455 Tf_dif_b = c1 - Rf_dif_b - ! the k = kfrsnl layer properties are updated to combined + ! the k = kfrsnl layer properties are updated to combined ! the fresnel (refractive) layer, always taken to be above ! the present layer k (i.e. be the top interface): @@ -4168,31 +4168,31 @@ subroutine solution_dEdd & endif ! k = kfrsnl endif ! trntdr(k) > trmin - + ! initialize current layer properties to zero; only if total ! transmission to the top interface of the current layer exceeds the ! minimum, will these values be computed below: ! Calculate the solar beam transmission, total transmission, and - ! reflectivity for diffuse radiation from below at interface k, + ! reflectivity for diffuse radiation from below at interface k, ! the top of the current layer k: ! ! layers interface - ! - ! --------------------- k-1 + ! + ! --------------------- k-1 ! k-1 ! --------------------- k ! k - ! --------------------- + ! --------------------- ! For k = klevp ! note that we ignore refraction between sea ice and underlying ocean: ! ! layers interface ! - ! --------------------- k-1 + ! --------------------- k-1 ! k-1 ! --------------------- k ! \\\\\\\ ocean \\\\\\\ - + trndir(k+1) = trndir(k)*trnlay(k) refkm1 = c1/(c1 - rdndif(k)*rdif_a(k)) tdrrdir = trndir(k)*rdir(k) @@ -4205,8 +4205,8 @@ subroutine solution_dEdd & enddo ! k end main level loop - ! compute reflectivity to direct and diffuse radiation for layers - ! below by adding succesive layers starting from the underlying + ! compute reflectivity to direct and diffuse radiation for layers + ! below by adding succesive layers starting from the underlying ! ocean and working upwards: ! ! layers interface @@ -4218,7 +4218,7 @@ subroutine solution_dEdd & ! --------------------- rupdir(klevp) = albodr - rupdif(klevp) = albodf + rupdif(klevp) = albodf do k=klev,0,-1 ! interface scattering @@ -4238,10 +4238,10 @@ end subroutine solution_dEdd !======================================================================= ! -! Set snow horizontal coverage, density and grain radius diagnostically +! Set snow horizontal coverage, density and grain radius diagnostically ! for the Delta-Eddington solar radiation method. ! -! author: Bruce P. Briegleb, NCAR +! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version subroutine shortwave_dEdd_set_snow(nslyr, R_snw, & @@ -4252,7 +4252,7 @@ subroutine shortwave_dEdd_set_snow(nslyr, R_snw, & rhosnw, rsnw, & rsnow) - integer (kind=int_kind), intent(in) :: & + integer (kind=int_kind), intent(in) :: & nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & @@ -4263,7 +4263,7 @@ subroutine shortwave_dEdd_set_snow(nslyr, R_snw, & real (kind=dbl_kind), intent(in) :: & aice , & ! concentration of ice vsno , & ! volume of snow - Tsfc , & ! surface temperature + Tsfc , & ! surface temperature hs0 ! snow depth for transition to bare sea ice (m) real (kind=dbl_kind), intent(inout) :: & @@ -4298,12 +4298,12 @@ subroutine shortwave_dEdd_set_snow(nslyr, R_snw, & ! set snow horizontal fraction hs = vsno / aice - + if (hs >= hs_min) then fs = c1 if (hs0 > puny) fs = min(hs/hs0, c1) endif - + if (snwgrain) then ! use snow grain tracer do ks = 1, nslyr @@ -4341,7 +4341,7 @@ end subroutine shortwave_dEdd_set_snow ! Set pond fraction and depth diagnostically for ! the Delta-Eddington solar radiation method. ! -! author: Bruce P. Briegleb, NCAR +! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version subroutine shortwave_dEdd_set_pond(Tsfc, & @@ -4375,7 +4375,7 @@ subroutine shortwave_dEdd_set_pond(Tsfc, & ! pond fp = 0.3_dbl_kind*fT*(c1-fs) hp = 0.3_dbl_kind*fT*(c1-fs) - + end subroutine shortwave_dEdd_set_pond ! End Delta-Eddington shortwave method @@ -4392,13 +4392,13 @@ subroutine compute_shortwave_trcr(nslyr, & nilyr, nblyr, & i_grid, & skl_bgc, z_tracers ) - + integer (kind=int_kind), intent(in) :: & nslyr ! number of snow layers integer (kind=int_kind), intent(in) :: & nblyr , & ! number of bio layers - nilyr ! number of ice layers + nilyr ! number of ice layers real (kind=dbl_kind), dimension (:), intent(in) :: & bgcN , & ! Nit tracer @@ -4408,9 +4408,9 @@ subroutine compute_shortwave_trcr(nslyr, & trcrn_bgcsw ! ice on shortwave grid tracers real (kind=dbl_kind), dimension (:), intent(in) :: & - sw_grid , & ! + sw_grid , & ! i_grid ! CICE bio grid - + real(kind=dbl_kind), intent(in) :: & hin , & ! CICE ice thickness hbri ! brine height @@ -4431,7 +4431,7 @@ subroutine compute_shortwave_trcr(nslyr, & icegrid ! correct for large ice surface layers real (kind=dbl_kind):: & - top_conc ! 1% (min_bgc) of surface concentration + top_conc ! 1% (min_bgc) of surface concentration ! when hin > hbri: just used in sw calculation character(len=*),parameter :: subname='(compute_shortwave_trcr)' @@ -4446,7 +4446,7 @@ subroutine compute_shortwave_trcr(nslyr, & do k = 1,nilyr+1 icegrid(k) = sw_grid(k) - enddo + enddo if (sw_grid(1)*hin*c2 > hi_ssl) then icegrid(1) = hi_ssl/c2/hin endif @@ -4465,7 +4465,7 @@ subroutine compute_shortwave_trcr(nslyr, & R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1 + k-1) enddo ! n enddo ! k - + top_conc = trtmp0(nt_bgc_N(1))*min_bgc call remap_zbgc (nilyr+1, & nt_bgc_N(1), & @@ -4474,7 +4474,7 @@ subroutine compute_shortwave_trcr(nslyr, & 1, nblyr+1, & hin, hbri, & icegrid(1:nilyr+1), & - i_grid(1:nblyr+1), top_conc ) + i_grid(1:nblyr+1), top_conc ) if (icepack_warnings_aborted(subname)) return do k = 1, nilyr+1 @@ -4533,7 +4533,7 @@ subroutine compute_shortwave_trcr(nslyr, & + F_abs_chl(nn)*R_chl2N(nn) & * bgcN(nt_bgc_N(nn)-nt_bgc_N(1)+1)*sk_l/hin & * real(nilyr,kind=dbl_kind) - enddo + enddo endif end subroutine compute_shortwave_trcr @@ -4604,10 +4604,10 @@ subroutine icepack_prep_radiation (ncat, nilyr, nslyr, & ! local variables integer (kind=int_kind) :: & - k , & ! vertical index + k , & ! vertical index n ! thickness category index - real (kind=dbl_kind) :: netsw + real (kind=dbl_kind) :: netsw character(len=*),parameter :: subname='(icepack_prep_radiation)' @@ -4736,14 +4736,14 @@ subroutine icepack_step_radiation (dt, ncat, & yday ! day of the year real (kind=dbl_kind), intent(inout) :: & - coszen ! cosine solar zenith angle, < 0 for sun below horizon + coszen ! cosine solar zenith angle, < 0 for sun below horizon real (kind=dbl_kind), dimension (:), intent(in) :: & igrid ! biology vertical interface points - + real (kind=dbl_kind), dimension (:), intent(in) :: & swgrid ! grid for ice tracers used in dEdd scheme - + real (kind=dbl_kind), dimension(:,:), intent(in) :: & kaer_3bd, & ! aerosol mass extinction cross section (m2/kg) waer_3bd, & ! aerosol single scatter albedo (fraction) @@ -4766,7 +4766,7 @@ subroutine icepack_step_radiation (dt, ncat, & apndn , & ! pond area fraction hpndn , & ! pond depth (m) ipndn , & ! pond refrozen lid thickness (m) - fbri ! brine fraction + fbri ! brine fraction real(kind=dbl_kind), dimension(:,:), intent(in) :: & aeron , & ! aerosols (kg/m^3) @@ -4788,9 +4788,9 @@ subroutine icepack_step_radiation (dt, ncat, & dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond ! albedo components for history - albicen , & ! bare ice - albsnon , & ! snow - albpndn , & ! pond + albicen , & ! bare ice + albsnon , & ! snow + albpndn , & ! pond apeffn ! effective pond area used for radiation calculation real (kind=dbl_kind), dimension(:), intent(inout), optional :: & @@ -4893,7 +4893,7 @@ subroutine icepack_step_radiation (dt, ncat, & if (calc_Tsfc) then if (trim(shortwave(1:4)) == 'dEdd') then ! delta Eddington - + call run_dEdd(dt, ncat, & dEdd_algae, & nilyr, nslyr, & @@ -4948,7 +4948,7 @@ subroutine icepack_step_radiation (dt, ncat, & l_use_snicar=use_snicar, & initonly=linitonly) if (icepack_warnings_aborted(subname)) return - + elseif (trim(shortwave(1:4)) == 'ccsm') then call shortwave_ccsm3(aicen, vicen, & @@ -4995,7 +4995,7 @@ subroutine icepack_step_radiation (dt, ncat, & if (tr_pond_topo) then do n = 1, ncat - apeffn(n) = c0 + apeffn(n) = c0 if (aicen(n) > puny) then ! Lid effective if thicker than hp1 if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then @@ -5006,7 +5006,7 @@ subroutine icepack_step_radiation (dt, ncat, & if (apndn(n) < puny) apeffn(n) = c0 endif enddo ! ncat - + endif ! tr_pond_topo ! Initialize for safety @@ -5120,7 +5120,7 @@ real(kind=dbl_kind) function asys(gg,f) asys = (gg - f)/(c1 - f) end function asys - + !======================================================================= end module icepack_shortwave