diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 2256abfe2..71441de2b 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -231,11 +231,11 @@ subroutine shortwave_ccsm3 (aicen, vicen, & alvdrn, alidrn, & alvdfn, alidfn, & fswsfc, fswint, & - fswthrun, & - fswthrun_vdr, & - fswthrun_vdf, & - fswthrun_idr, & - fswthrun_idf, & + fswthrun, & + fswthrun_vdr, & + fswthrun_vdf, & + fswthrun_idr, & + fswthrun_idf, & fswpenl, & Iswabs, SSwabs, & albin, albsn, & @@ -271,15 +271,15 @@ subroutine shortwave_ccsm3 (aicen, vicen, & alidfn , & ! near-ir, diffuse, avg (fraction) fswsfc , & ! SW absorbed at ice/snow surface (W m-2) fswint , & ! SW absorbed in ice interior, below surface (W m-2) - fswthrun , & ! SW through ice to ocean (W m-2) + fswthrun , & ! SW through ice to ocean (W m-2) albin , & ! bare ice albedo albsn ! snow albedo 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), intent(inout) :: & coszen ! cosine(zenith angle) @@ -308,10 +308,10 @@ subroutine shortwave_ccsm3 (aicen, vicen, & ! needed for optional fswthrun arrays when passed as scalars real (kind=dbl_kind) :: & - l_fswthru_vdr , & ! vis dir SW through ice to ocean (W m-2) - l_fswthru_vdf , & ! vis dif SW through ice to ocean (W m-2) - l_fswthru_idr , & ! nir dir SW through ice to ocean (W m-2) - l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) + l_fswthru_vdr, & ! vis dir SW through ice to ocean (W m-2) + l_fswthru_vdf, & ! vis dif SW through ice to ocean (W m-2) + l_fswthru_idr, & ! nir dir SW through ice to ocean (W m-2) + l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) character(len=*),parameter :: subname='(shortwave_ccsm3)' @@ -324,7 +324,7 @@ subroutine shortwave_ccsm3 (aicen, vicen, & do n = 1, ncat - Sswabs(:,n) = c0 + Sswabs(:,n) = c0 alvdrni = albocn alidrni = albocn @@ -344,10 +344,10 @@ subroutine shortwave_ccsm3 (aicen, vicen, & albin(n) = c0 albsn(n) = c0 - fswsfc(n) = c0 - fswint(n) = c0 - fswthrun(n) = c0 - fswpenl(:,n) = c0 + fswsfc(n) = c0 + fswint(n) = c0 + fswthrun(n) = c0 + fswpenl(:,n) = c0 Iswabs (:,n) = c0 if (aicen(n) > puny) then @@ -417,20 +417,20 @@ subroutine shortwave_ccsm3 (aicen, vicen, & alidrns, alidfns, & fswsfc=fswsfc(n), & fswint=fswint(n), & - fswthru=fswthrun(n), & - fswthru_vdr=l_fswthru_vdr,& - fswthru_vdf=l_fswthru_vdf,& - fswthru_idr=l_fswthru_idr,& - fswthru_idf=l_fswthru_idf,& + fswthru=fswthrun(n), & + fswthru_vdr=l_fswthru_vdr, & + fswthru_vdf=l_fswthru_vdf, & + fswthru_idr=l_fswthru_idr, & + fswthru_idf=l_fswthru_idf, & fswpenl=fswpenl(:,n), & Iswabs=Iswabs(:,n)) if (icepack_warnings_aborted(subname)) return - if(present(fswthrun_vdr)) fswthrun_vdr(n) = l_fswthru_vdr - if(present(fswthrun_vdf)) fswthrun_vdf(n) = l_fswthru_vdf - if(present(fswthrun_idr)) fswthrun_idr(n) = l_fswthru_idr - if(present(fswthrun_idf)) fswthrun_idf(n) = l_fswthru_idf + if (present(fswthrun_vdr)) fswthrun_vdr(n) = l_fswthru_vdr + if (present(fswthrun_vdf)) fswthrun_vdf(n) = l_fswthru_vdf + if (present(fswthrun_idr)) fswthrun_idr(n) = l_fswthru_idr + if (present(fswthrun_idf)) fswthrun_idf(n) = l_fswthru_idf endif ! aicen > puny @@ -488,7 +488,7 @@ subroutine compute_albedos (aicen, vicen, & ! 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 @@ -509,8 +509,6 @@ subroutine compute_albedos (aicen, vicen, & character(len=*),parameter :: subname='(compute_albedos)' - fhtan = atan(ahmax*c4) - !----------------------------------------------------------------- ! Compute albedo for each thickness category. !----------------------------------------------------------------- @@ -519,6 +517,7 @@ subroutine compute_albedos (aicen, vicen, & hs = vsnon / aicen ! bare ice, thickness dependence + fhtan = atan(ahmax*c4) fh = min(atan(hi*c4)/fhtan,c1) albo = albocn*(c1-fh) alvdfni = albicev*fh + albo @@ -757,8 +756,8 @@ subroutine absorbed_solar (aicen, & ! Initialize !----------------------------------------------------------------- - trantop = c0 - tranbot = c0 + trantop = c0 + tranbot = c0 hs = vsnon / aicen @@ -873,13 +872,13 @@ subroutine run_dEdd(dt, & vsnon, Tsfcn, & alvln, apndn, & hpndn, ipndn, & - aeron, & + aeron, & trcrn_bgcsw, & TLAT, TLON, & calendar_type, & days_per_year, & nextsw_cday, yday, & - sec, & + sec, & swvdr, swvdf, & swidr, swidf, & coszen, fsnow, & @@ -937,15 +936,15 @@ subroutine run_dEdd(dt, & ipndn ! pond refrozen lid thickness (m) real(kind=dbl_kind), dimension(:,:), intent(in) :: & - aeron, & ! aerosols (kg/m^3) + aeron, & ! aerosols (kg/m^3) trcrn_bgcsw ! zaerosols (kg/m^3) + chlorophyll on shorthwave grid real(kind=dbl_kind), dimension(:), intent(inout) :: & - ffracn,& ! fraction of fsurfn used to melt ipond - dhsn ! depth difference for snow on sea ice and pond ice + ffracn, & ! fraction of fsurfn used to melt ipond + 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) @@ -976,14 +975,12 @@ subroutine run_dEdd(dt, & rsnow ! snow grain radius tracer (10^-6 m) logical (kind=log_kind), intent(in) :: & - l_print_point + l_print_point ! print diagnostic information logical (kind=log_kind), optional :: & initonly ! flag to indicate init only, default is false - ! local temporary variables - - ! other local variables + ! local variables ! snow variables for Delta-Eddington shortwave real (kind=dbl_kind) :: & fsn , & ! snow horizontal fraction @@ -1014,7 +1011,6 @@ subroutine run_dEdd(dt, & hmx , & ! maximum available snow infiltration equivalent depth dhs , & ! local difference in snow depth on sea ice and pond ice spn , & ! snow depth on refrozen pond (m) - rnslyr , & ! 1/nslyr tmp ! 0 or 1 ! needed for optional fswthrun arrays when passed as scalars @@ -1025,7 +1021,7 @@ subroutine run_dEdd(dt, & l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) logical (kind=log_kind) :: & - l_initonly ! local initonly value + l_initonly ! local initonly value real(kind=dbl_kind), dimension(nslyr) :: & l_rsnows ! snow grain radius tracer (10^-6 m) @@ -1058,15 +1054,15 @@ subroutine run_dEdd(dt, & hsn = c0 rhosnwn(:) = c0 rsnwn(:) = c0 - apeffn(n) = c0 ! for history - snowfracn(n) = c0 ! for history + apeffn(n) = c0 ! for history + snowfracn(n) = c0 ! for history if (aicen(n) > puny) then if (snwgrain) then l_rsnows(:) = rsnow(:,n) endif - call shortwave_dEdd_set_snow(R_snw, & + call shortwave_dEdd_set_snow(R_snw, & dT_mlt, rsnw_mlt, & aicen(n), vsnon(n), & Tsfcn(n), fsn, & @@ -1382,26 +1378,26 @@ subroutine shortwave_dEdd (coszen, & vsno ! volume of snow real (kind=dbl_kind) :: & - swdn , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j) - swab , & ! fswsfc(i,j)+fswint(i,j)+fswthru(i,j) - swalb ! (1.-swab/(swdn+.0001)) + swdn , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j) + swab , & ! fswsfc(i,j)+fswint(i,j)+fswthru(i,j) + swalb ! (1.-swab/(swdn+.0001)) ! 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)' !----------------------------------------------------------------------- - klev = nslyr + nilyr + 1 ! number of radiation layers - 1 - klevp = klev + 1 ! number of radiation interfaces - 1 - ! (0 layer is included also) + klev = nslyr + nilyr + 1 ! number of radiation layers - 1 + klevp = klev + 1 ! number of radiation interfaces - 1 + ! (0 layer is included also) - ! zero storage albedos and fluxes for accumulation over surface types: + ! set storage albedos and fluxes to zero for accumulation over surface types hstmp = c0 hi = c0 fi = c0 @@ -1425,12 +1421,12 @@ subroutine shortwave_dEdd (coszen, & if( swidr + swidf > puny ) then fnidr = swidr/(swidr+swidf) endif - albice = c0 - albsno = c0 - albpnd = c0 + albice = c0 + albsno = c0 + albpnd = c0 fswpenl(:) = c0 - Sswabs(:) = c0 - Iswabs(:) = c0 + Sswabs (:) = c0 + Iswabs (:) = c0 ! compute aerosol mass path @@ -1471,26 +1467,20 @@ subroutine shortwave_dEdd (coszen, & srftyp = 0 call compute_dEdd_3bd( & - klev, klevp, zbio, & - fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hstmp, rhosnw, rsnw, hi, hp, & - fi, aero_mp, avdrl, avdfl, & - aidrl, aidfl, & - fswsfc, fswint, & - fswthru, & - fswthru_vdr, & - fswthru_vdf, & - fswthru_idr, & - fswthru_idf, & - Sswabs, & - Iswabs, fswpenl ) + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hstmp, rhosnw, rsnw, hi, hp, & + fi, aero_mp, avdrl, avdfl, & + aidrl, aidfl, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) if (icepack_warnings_aborted(subname)) return - alvdr = alvdr + avdrl *fi - alvdf = alvdf + avdfl *fi - alidr = alidr + aidrl *fi - alidf = alidf + aidfl *fi + alvdr = alvdr + avdrl*fi + alvdf = alvdf + avdfl*fi + alidr = alidr + aidrl*fi + alidf = alidf + aidfl*fi ! for history albice = albice & + awtvdr*avdrl + awtidr*aidrl & @@ -1508,41 +1498,34 @@ subroutine shortwave_dEdd (coszen, & srftyp = 1 if (trim(shortwave) == 'dEdd_snicar_ad') then - call compute_dEdd_5bd(klev, klevp, & - zbio, & - fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hs, rhosnw, rsnw, hi, hp, & - fs, aero_mp, avdrl, avdfl, & - aidrl, aidfl, & - fswsfc, fswint, & - fswthru, Sswabs, & - Iswabs, fswpenl ) + call compute_dEdd_5bd( & + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hs, rhosnw, rsnw, hi, hp, & + fs, aero_mp, avdrl, avdfl, & + aidrl, aidfl, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) else !echmod - this can be combined with the 5bd call above, if we use module data - call compute_dEdd_3bd( & - klev, klevp, zbio, & - fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hs, rhosnw, rsnw, hi, hp, & - fs, aero_mp, avdrl, avdfl, & - aidrl, aidfl, & - fswsfc, fswint, & - fswthru, & - fswthru_vdr, & - fswthru_vdf, & - fswthru_idr, & - fswthru_idf, & - Sswabs, & - Iswabs, fswpenl ) + call compute_dEdd_3bd( & + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hs, rhosnw, rsnw, hi, hp, & + fs, aero_mp, avdrl, avdfl, & + aidrl, aidfl, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) endif if (icepack_warnings_aborted(subname)) return - alvdr = alvdr + avdrl *fs - alvdf = alvdf + avdfl *fs - alidr = alidr + aidrl *fs - alidf = alidf + aidfl *fs + alvdr = alvdr + avdrl*fs + alvdf = alvdf + avdfl*fs + alidr = alidr + aidrl*fs + alidf = alidf + aidfl*fs ! for history albsno = albsno & + awtvdr*avdrl + awtidr*aidrl & @@ -1564,27 +1547,21 @@ subroutine shortwave_dEdd (coszen, & ! calculate ponded ice srftyp = 2 - call compute_dEdd_3bd( & - klev, klevp, zbio, & - fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hs, rhosnw, rsnw, hi, hp, & - fp, aero_mp, avdrl, avdfl, & - aidrl, aidfl, & - fswsfc, fswint, & - fswthru, & - fswthru_vdr, & - fswthru_vdf, & - fswthru_idr, & - fswthru_idf, & - Sswabs, & - Iswabs, fswpenl ) + call compute_dEdd_3bd( & + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hs, rhosnw, rsnw, hi, hp, & + fp, aero_mp, avdrl, avdfl, & + aidrl, aidfl, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) if (icepack_warnings_aborted(subname)) return - alvdr = alvdr + avdrl *fp - alvdf = alvdf + avdfl *fp - alidr = alidr + aidrl *fp - alidf = alidf + aidfl *fp + alvdr = alvdr + avdrl*fp + alvdf = alvdf + avdfl*fp + alidr = alidr + aidrl*fp + alidf = alidf + aidfl*fp ! for history albpnd = albpnd & + awtvdr*avdrl + awtidr*aidrl & @@ -1680,21 +1657,15 @@ end subroutine shortwave_dEdd ! 2013: E Hunke merged with NCAR version ! 2022: E Hunke, T Craig moved data (now module data) - subroutine compute_dEdd_3bd( & - klev, klevp, zbio, & - fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hs, rhosnw, rsnw, hi, hp, & - fi, aero_mp, alvdr, alvdf, & - alidr, alidf, & - fswsfc, fswint, & - fswthru, & - fswthru_vdr, & - fswthru_vdf, & - fswthru_idr, & - fswthru_idf, & - Sswabs, & - Iswabs, fswpenl ) + subroutine compute_dEdd_3bd( & + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hs, rhosnw, rsnw, hi, hp, & + fi, aero_mp, alvdr, alvdf, & + alidr, alidf, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) integer (kind=int_kind), intent(in) :: & klev , & ! number of radiation layers - 1 @@ -1702,49 +1673,49 @@ subroutine compute_dEdd_3bd( & ! (0 layer is included also) real (kind=dbl_kind), intent(in) :: & - fnidr , & ! fraction of direct to total down flux in nir - coszen , & ! cosine solar zenith angle - swvdr , & ! shortwave down at surface, visible, direct (W/m^2) - 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) + fnidr , & ! fraction of direct to total down flux in nir + coszen, & ! cosine solar zenith angle + swvdr , & ! shortwave down at surface, visible, direct (W/m^2) + 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) :: & - srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) + srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) real (kind=dbl_kind), intent(in) :: & - hs ! snow thickness (m) + hs ! snow thickness (m) real (kind=dbl_kind), dimension (:), intent(in) :: & - rhosnw , & ! snow density in snow layer (kg/m3) - rsnw , & ! snow grain radius in snow layer (m) - zbio , & ! zaerosol + chla shortwave tracers kg/m^3 - aero_mp ! aerosol mass path in kg/m2 + rhosnw, & ! snow density in snow layer (kg/m3) + rsnw , & ! snow grain radius in snow layer (m) + zbio , & ! zaerosol + chla shortwave tracers kg/m^3 + aero_mp ! aerosol mass path in kg/m2 real (kind=dbl_kind), intent(in) :: & - hi , & ! ice thickness (m) - hp , & ! pond depth (m) - fi ! snow/bare ice fractional coverage (0 to 1) + 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) - 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) + 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) real (kind=dbl_kind), intent(inout) :: & - fswthru_vdr , & ! vis dir SW through snow/bare ice/ponded ice into ocean (W m-2) - 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) + fswthru_vdr, & ! vis dir SW through snow/bare ice/ponded ice into ocean (W m-2) + 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) + 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) !----------------------------------------------------------------------- ! @@ -1886,12 +1857,12 @@ subroutine compute_dEdd_3bd( & fthrul ! shortwave through to ice layers (W m-2) real (kind=dbl_kind), dimension (nspint_3bd) :: & - wghtns ! spectral weights + wghtns ! spectral weights 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 + 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 @@ -1927,53 +1898,54 @@ subroutine compute_dEdd_3bd( & ! ice and ponded ice IOPs, allowing for tuning ! modifications of the above "_mn" value real (kind=dbl_kind), dimension (nspint_3bd) :: & - ki_ssl , & ! Surface-scattering-layer ice extinction coefficient (/m) - wi_ssl , & ! Surface-scattering-layer ice single scattering albedo - gi_ssl , & ! Surface-scattering-layer ice asymmetry parameter - ki_dl , & ! Drained-layer ice extinction coefficient (/m) - wi_dl , & ! Drained-layer ice single scattering albedo - gi_dl , & ! Drained-layer ice asymmetry parameter - ki_int , & ! Interior-layer ice extinction coefficient (/m) - wi_int , & ! Interior-layer ice single scattering albedo - gi_int , & ! Interior-layer ice asymmetry parameter - ki_p_ssl , & ! Ice under pond srf scat layer extinction coefficient (/m) - wi_p_ssl , & ! Ice under pond srf scat layer single scattering albedo - gi_p_ssl , & ! Ice under pond srf scat layer asymmetry parameter - ki_p_int , & ! Ice under pond extinction coefficient (/m) - wi_p_int , & ! Ice under pond single scattering albedo - gi_p_int ! Ice under pond asymmetry parameter + ki_ssl , & ! Surface-scattering-layer ice extinction coefficient (/m) + wi_ssl , & ! Surface-scattering-layer ice single scattering albedo + gi_ssl , & ! Surface-scattering-layer ice asymmetry parameter + ki_dl , & ! Drained-layer ice extinction coefficient (/m) + wi_dl , & ! Drained-layer ice single scattering albedo + gi_dl , & ! Drained-layer ice asymmetry parameter + ki_int , & ! Interior-layer ice extinction coefficient (/m) + wi_int , & ! Interior-layer ice single scattering albedo + gi_int , & ! Interior-layer ice asymmetry parameter + ki_p_ssl, & ! Ice under pond srf scat layer extinction coefficient (/m) + wi_p_ssl, & ! Ice under pond srf scat layer single scattering albedo + gi_p_ssl, & ! Ice under pond srf scat layer asymmetry parameter + ki_p_int, & ! Ice under pond extinction coefficient (/m) + wi_p_int, & ! Ice under pond single scattering albedo + gi_p_int ! Ice under pond asymmetry parameter real (kind=dbl_kind), dimension(0:klev) :: & - dzk ! layer thickness + dzk ! layer thickness real (kind=dbl_kind) :: & - dz , & ! snow, sea ice or pond water layer thickness - 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 - ! conserve computed albedo for constant physical depth of sea - ! ice when the number of sea ice layers vary + dz , & ! snow, sea ice or pond water layer thickness + 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 + ! conserve computed albedo for constant physical depth of sea + ! ice when the number of sea ice layers vary + real (kind=dbl_kind) :: & - sig , & ! scattering coefficient for tuning - kabs , & ! absorption coefficient for tuning - sigp ! modified scattering coefficient for tuning + sig , & ! scattering coefficient for tuning + kabs , & ! absorption coefficient for tuning + sigp ! modified scattering coefficient for tuning real (kind=dbl_kind), dimension(nspint_3bd, 0:klev) :: & - kabs_chl , & ! absorption coefficient for chlorophyll (/m) - tzaer , & ! total aerosol extinction optical depth - wzaer , & ! total aerosol single scatter albedo - gzaer ! total aerosol asymmetry parameter + kabs_chl, & ! absorption coefficient for chlorophyll (/m) + tzaer , & ! total aerosol extinction optical depth + wzaer , & ! total aerosol single scatter albedo + gzaer ! total aerosol asymmetry parameter real (kind=dbl_kind) :: & - albodr , & ! spectral ocean albedo to direct rad - albodf ! spectral ocean albedo to diffuse rad + 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 real (kind=dbl_kind) :: & - sig_i , & ! ice scattering coefficient (/m) - sig_p , & ! pond scattering coefficient (/m) - kext ! weighted extinction coefficient (/m) + sig_i , & ! ice scattering coefficient (/m) + sig_p , & ! pond scattering coefficient (/m) + kext ! weighted extinction coefficient (/m) ! aerosol optical properties from Mark Flanner, 26 June 2008 ! order assumed: hydrophobic black carbon, hydrophilic black carbon, @@ -1984,26 +1956,26 @@ subroutine compute_dEdd_3bd( & ! and 1.19-5.0 micron in wavelength) integer (kind=int_kind) :: & - na , n ! aerosol index + na , n ! aerosol index real (kind=dbl_kind) :: & - taer , & ! total aerosol extinction optical depth - waer , & ! total aerosol single scatter albedo - gaer , & ! total aerosol asymmetry parameter - swdr , & ! shortwave down at surface, direct (W/m^2) - swdf , & ! shortwave down at surface, diffuse (W/m^2) - rnilyr , & ! real(nilyr) - rnslyr , & ! real(nslyr) - rns , & ! real(ns) - tmp_0, tmp_ks, tmp_kl ! temp variables + taer , & ! total aerosol extinction optical depth + waer , & ! total aerosol single scatter albedo + gaer , & ! total aerosol asymmetry parameter + swdr , & ! shortwave down at surface, direct (W/m^2) + swdf , & ! shortwave down at surface, diffuse (W/m^2) + rnilyr , & ! 1/real(nilyr) + rnslyr , & ! 1/real(nslyr) + rns , & ! real(ns) + tmp_0, tmp_ks, tmp_kl ! temporary variables integer(kind=int_kind), dimension(0:klev) :: & - k_bcini , & ! index - k_bcins , & - k_bcexs + k_bcini , & ! index + k_bcins , & ! = 2 hardwired + k_bcexs ! = 2 hardwired real(kind=dbl_kind):: & - tmp_gs, tmp1 ! temp variables + tmp_gs, tmp1 ! temporary variables real (kind=dbl_kind), parameter :: & fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max @@ -2015,13 +1987,13 @@ subroutine compute_dEdd_3bd( & fp_pnd = 2.00_dbl_kind, & ! ponded ice fraction of scat coeff for + stn dev in alb fm_pnd = 0.50_dbl_kind ! ponded ice fraction of scat coeff for - stn dev in alb - real (kind=dbl_kind), parameter :: & !chla-specific absorption coefficient - kchl_tab = 0.01 !0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1 - ! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014 - ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011) - ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014) - ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting) - !chlorophyll mass extinction cross section (m^2/mg chla) + real (kind=dbl_kind), parameter :: & ! chla-specific absorption coefficient + kchl_tab = p01 ! 0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1 + ! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014 + ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011) + ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014) + ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting) + ! chlorophyll mass extinction cross section (m^2/mg chla) character(len=*),parameter :: subname='(compute_dEdd_3bd)' @@ -2037,24 +2009,24 @@ subroutine compute_dEdd_3bd( & kii = nslyr + 1 ! initialize albedos and fluxes to 0 - fthrul = c0 - Iabs = c0 + fthrul = c0 + Iabs = c0 kabs_chl(:,:) = c0 - tzaer(:,:) = c0 - wzaer(:,:) = c0 - gzaer(:,:) = c0 - - avdr = c0 - avdf = c0 - aidr = c0 - aidf = c0 - fsfc = c0 - fint = c0 - fthru = c0 - fthruvdr = c0 - fthruvdf = c0 - fthruidr = c0 - fthruidf = c0 + tzaer (:,:) = c0 + wzaer (:,:) = c0 + gzaer (:,:) = c0 + + avdr = c0 + avdf = c0 + aidr = c0 + aidf = c0 + fsfc = c0 + fint = c0 + fthru = c0 + fthruvdr = c0 + fthruvdf = c0 + fthruidr = c0 + fthruidf = c0 ! spectral 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 @@ -2200,108 +2172,123 @@ subroutine compute_dEdd_3bd( & ! aerosols if (modal_aero) then - 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. - ! CICE adjusted snow grain radius rsnw to frsnw in the original 3-band - ! scheme, while for SNICAR the snow grain radius is used directly. - ksnow = k - min(k-1,0) - tmp_gs = frsnw(ksnow) - - ! grain size index - if (tmp_gs < 125._dbl_kind) then - tmp1 = tmp_gs/50._dbl_kind - k_bcini(k) = nint(tmp1) - elseif (tmp_gs < 175._dbl_kind) then - k_bcini(k) = 2 - else - tmp1 = (tmp_gs/250._dbl_kind) + c2 - k_bcini(k) = nint(tmp1) - endif - else ! use the largest snow grain size for ice - k_bcini(k) = 8 - endif - ! Set index corresponding to BC effective radius. Here, - ! asssume constant BC effective radius of 100nm - ! (corresponding to index 2) - k_bcins(k) = 2 - k_bcexs(k) = 2 - - ! check bounds: - if (k_bcini(k) < 1) k_bcini(k) = 1 - if (k_bcini(k) > 8) k_bcini(k) = 8 - if (k_bcins(k) < 1) k_bcins(k) = 1 - if (k_bcins(k) > 10) k_bcins(k) = 10 - if (k_bcexs(k) < 1) k_bcexs(k) = 1 - if (k_bcexs(k) > 10) k_bcexs(k) = 10 - - ! print ice radius index: - ! write(warnstr,*) subname, "MGFICE2:k, ice index= ",k, k_bcini(k) - ! call icepack_warnings_add(warnstr) - enddo ! k - - 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_3bd ! not weighted by aice - tzaer(ns,k) = tzaer(ns,k)+kaer_bc_3bd(ns,k_bcexs(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer(ns,k) = wzaer(ns,k)+kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer(ns,k) = gzaer(ns,k)+kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k))* & - gaer_bc_3bd(ns,k_bcexs(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k) - enddo ! nspint - enddo - elseif (n==2) then ! within-ice BC - do k = 0, klev - do ns = 1,nspint_3bd - tzaer(ns,k) = tzaer(ns,k)+kaer_bc_3bd(ns,k_bcins(k)) * & - bcenh_3bd(ns,k_bcins(k),k_bcini(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer(ns,k) = wzaer(ns,k)+kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer(ns,k) = gzaer(ns,k)+kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k))* & - gaer_bc_3bd(ns,k_bcins(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k) - enddo ! nspint - enddo - else ! dust - 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. + ! CICE adjusted snow grain radius rsnw to frsnw in the original 3-band + ! scheme, while for SNICAR the snow grain radius is used directly. + ksnow = k - min(k-1,0) + tmp_gs = frsnw(ksnow) + + ! grain size index + if (tmp_gs < 125._dbl_kind) then + tmp1 = tmp_gs/50._dbl_kind + k_bcini(k) = nint(tmp1) + elseif (tmp_gs < 175._dbl_kind) then + k_bcini(k) = 2 + else + tmp1 = (tmp_gs/250._dbl_kind) + c2 + k_bcini(k) = nint(tmp1) + endif + else ! use the largest snow grain size for ice + k_bcini(k) = 8 + endif + ! Set index corresponding to BC effective radius. Here, + ! asssume constant BC effective radius of 100nm + ! (corresponding to index 2) + k_bcins(k) = 2 ! hardwired + k_bcexs(k) = 2 + + ! check bounds + if (k_bcini(k) < 1) k_bcini(k) = 1 + if (k_bcini(k) > 8) k_bcini(k) = 8 +! if (k_bcins(k) < 1) k_bcins(k) = 1 ! hardwired +! if (k_bcins(k) > 10) k_bcins(k) = 10 +! if (k_bcexs(k) < 1) k_bcexs(k) = 1 +! if (k_bcexs(k) > 10) k_bcexs(k) = 10 + enddo ! k + + 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_3bd ! not weighted by aice + tzaer(ns,k) = tzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer(ns,k) = wzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer(ns,k) = gzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) & + * gaer_bc_3bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + enddo + enddo + elseif (n==2) then ! within-ice BC + do k = 0, klev + do ns = 1, nspint_3bd + tzaer(ns,k) = tzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcins(k)) & + * bcenh_3bd(ns,k_bcins(k),k_bcini(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer(ns,k) = wzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer(ns,k) = gzaer (ns,k) & + + kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) & + * gaer_bc_3bd(ns,k_bcins(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + enddo + enddo + else ! dust + do k = 0, klev do ns = 1,nspint_3bd ! not weighted by aice - tzaer(ns,k) = tzaer(ns,k)+kaer_3bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer(ns,k) = wzaer(ns,k)+kaer_3bd(ns,n)*waer_3bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer(ns,k) = gzaer(ns,k)+kaer_3bd(ns,n)*waer_3bd(ns,n)* & - gaer_3bd(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k) + tzaer(ns,k) = tzaer (ns,k) & + + kaer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer(ns,k) = wzaer (ns,k) & + + kaer_3bd(ns,n) & + * waer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer(ns,k) = gzaer (ns,k) & + + kaer_3bd(ns,n) & + * waer_3bd(ns,n) & + * gaer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) enddo ! nspint - enddo - endif !(n=1) + enddo ! k + endif ! n enddo ! n_zaero endif ! tr_zaero and dEdd_algae - else ! Bulk aerosol treatment - if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll - do n = 1,n_zaero ! multiply by aice? + else ! Bulk aerosol treatment + if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll + do n = 1, n_zaero ! multiply by aice? do k = 0, klev - do ns = 1,nspint_3bd ! not weighted by aice - tzaer(ns,k) = tzaer(ns,k)+kaer_3bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer(ns,k) = wzaer(ns,k)+kaer_3bd(ns,n)*waer_3bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer(ns,k) = gzaer(ns,k)+kaer_3bd(ns,n)*waer_3bd(ns,n)* & - gaer_3bd(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k) - enddo ! nspint - enddo - enddo - endif !tr_zaero - - endif ! modal_aero + do ns = 1, nspint_3bd ! not weighted by aice + tzaer(ns,k) = tzaer (ns,k) & + + kaer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer(ns,k) = wzaer (ns,k) & + + kaer_3bd(ns,n) & + * waer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer(ns,k) = gzaer (ns,k) & + + kaer_3bd(ns,n) & + * waer_3bd(ns,n) & + * gaer_3bd(ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + enddo ! nspint + enddo ! k + enddo ! n + endif ! tr_zaero + endif ! modal_aero !----------------------------------------------------------------------- @@ -2311,42 +2298,42 @@ subroutine compute_dEdd_3bd( & ! set optical properties of air/snow/pond overlying sea ice ! air - if( srftyp == 0 ) then + if (srftyp == 0 ) then do k=0,nslyr tau(k) = c0 w0(k) = c0 g(k) = c0 enddo ! snow - else if( srftyp == 1 ) then + elseif (srftyp == 1 ) then ! interpolate snow iops using input snow grain radius, ! snow density and tabular data - do k=0,nslyr + do k = 0, nslyr ! use top rsnw, rhosnw for snow ssl and rest of top layer ksnow = k - min(k-1,0) ! find snow iops using input snow density and snow grain radius: - if( frsnw(ksnow) < rsnw_tab(1) ) then - Qs = Qs_tab(ns,1) - ws = ws_tab(ns,1) - gs = gs_tab(ns,1) - else if( frsnw(ksnow) >= rsnw_tab(nmbrad_snw) ) then - Qs = Qs_tab(ns,nmbrad_snw) - ws = ws_tab(ns,nmbrad_snw) - gs = gs_tab(ns,nmbrad_snw) + if (frsnw(ksnow) < rsnw_tab(1)) then + Qs = Qs_tab(ns,1) + ws = ws_tab(ns,1) + gs = gs_tab(ns,1) + elseif (frsnw(ksnow) >= rsnw_tab(nmbrad_snw)) then + Qs = Qs_tab(ns,nmbrad_snw) + ws = ws_tab(ns,nmbrad_snw) + gs = gs_tab(ns,nmbrad_snw) else ! linear interpolation in rsnw - do nr=2,nmbrad_snw - if( rsnw_tab(nr-1) <= frsnw(ksnow) .and. & + do nr = 2, nmbrad_snw + if (rsnw_tab(nr-1) <= frsnw(ksnow) .and. & frsnw(ksnow) < rsnw_tab(nr)) then delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / & (rsnw_tab(nr) - rsnw_tab(nr-1)) Qs = Qs_tab(ns,nr-1)*(c1-delr) + & - Qs_tab(ns,nr)*delr + Qs_tab(ns,nr )* delr ws = ws_tab(ns,nr-1)*(c1-delr) + & - ws_tab(ns,nr)*delr + ws_tab(ns,nr )* delr gs = gs_tab(ns,nr-1)*(c1-delr) + & - gs_tab(ns,nr)*delr + gs_tab(ns,nr )* delr endif enddo ! nr endif @@ -2354,8 +2341,8 @@ subroutine compute_dEdd_3bd( & (4._dbl_kind*frsnw(ksnow)*1.0e-6_dbl_kind)) tau(k) = (ks + kabs_chl(ns,k))*dzk(k) - w0(k) = ks/(ks + kabs_chl(ns,k)) *ws - g(k) = gs + w0 (k) = ks/(ks + kabs_chl(ns,k)) * ws + g (k) = gs enddo ! k ! aerosol in snow @@ -2373,7 +2360,7 @@ subroutine compute_dEdd_3bd( & waer = c0 gaer = c0 - do na=1,4*n_aero,4 + do na = 1, 4*n_aero, 4 if (modal_aero) then if (na == 1) then ! interstitial BC taer = taer + aero_mp(na)*kaer_bc_3bd(ns,k_bcexs(k)) @@ -2383,13 +2370,13 @@ subroutine compute_dEdd_3bd( & *waer_bc_3bd(ns,k_bcexs(k)) & *gaer_bc_3bd(ns,k_bcexs(k)) elseif (na == 5) then ! within-ice BC - taer = taer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & - *bcenh_3bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & - *waer_bc_3bd(ns,k_bcins(k)) - gaer = gaer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & - *waer_bc_3bd(ns,k_bcins(k)) & - *gaer_bc_3bd(ns,k_bcins(k)) + taer = taer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & + * bcenh_3bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & + *waer_bc_3bd(ns,k_bcins(k)) + gaer = gaer + aero_mp(na)*kaer_bc_3bd(ns,k_bcins(k)) & + *waer_bc_3bd(ns,k_bcins(k)) & + *gaer_bc_3bd(ns,k_bcins(k)) else ! other species (dust) taer = taer + aero_mp(na)*kaer_3bd(ns,(1+(na-1)/4)) waer = waer + aero_mp(na)*kaer_3bd(ns,(1+(na-1)/4)) & @@ -2405,65 +2392,66 @@ subroutine compute_dEdd_3bd( & gaer = gaer + aero_mp(na)*kaer_3bd(ns,(1+(na-1)/4)) & *waer_3bd(ns,(1+(na-1)/4)) & *gaer_3bd(ns,(1+(na-1)/4)) - endif !modal_aero - enddo ! na + endif ! modal_aero + enddo ! na gaer = gaer/(waer+puny) waer = waer/(taer+puny) ! tcraig, why does the above section exist if taer=waer=gaer=0 below - do k=1,nslyr + do k = 1, nslyr taer = c0 waer = c0 gaer = c0 - do na=1,4*n_aero,4 + do na = 1, 4*n_aero, 4 if (modal_aero) then - if (na==1) then ! interstitial BC - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcexs(k)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k))*gaer_bc_3bd(ns,k_bcexs(k)) - elseif (na==5) then ! within-ice BC - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcins(k))*& - bcenh_3bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k))*gaer_bc_3bd(ns,k_bcins(k)) - - else ! other species (dust) - - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) - endif ! na==1 - + if (na==1) then ! interstitial BC + taer = taer + (aero_mp(na+1)/rnslyr) & !echmod: BUG should be *rnslyr + * kaer_bc_3bd(ns,k_bcexs(k)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) & + * gaer_bc_3bd(ns,k_bcexs(k)) + elseif (na==5) then ! within-ice BC + taer = taer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * bcenh_3bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) & + * gaer_bc_3bd(ns,k_bcins(k)) + else ! other species (dust) + taer = taer + (aero_mp(na+1)/rnslyr) & + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) + endif ! na else - taer = taer + & - (aero_mp(na+1)*rnslyr)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+1)*rnslyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+1)*rnslyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) + taer = taer + (aero_mp(na+1)*rnslyr) & !echmod NOTE * rnslyr not / + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+1)*rnslyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+1)*rnslyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) endif ! modal_aero enddo ! na - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & + w0 (k) = (w0(k)*tau(k) + waer) / & (tau(k) + taer) tau(k) = tau(k) + taer enddo ! k @@ -2474,8 +2462,8 @@ subroutine compute_dEdd_3bd( & dz = hp/(c1/rnslyr+c1) do k=0,nslyr tau(k) = kw(ns)*dz - w0(k) = ww(ns) - g(k) = gw(ns) + w0 (k) = ww(ns) + g (k) = gw(ns) ! no aerosol in pond enddo ! k endif ! srftyp @@ -2483,40 +2471,40 @@ subroutine compute_dEdd_3bd( & ! set optical properties of sea ice ! bare or snow-covered sea ice layers - if( srftyp <= 1 ) then + if (srftyp <= 1) then ! ssl k = kii - tau(k) = (ki_ssl(ns)+kabs_chl(ns,k))*dzk(k) - w0(k) = ki_ssl(ns)/(ki_ssl(ns) + kabs_chl(ns,k))*wi_ssl(ns) - g(k) = gi_ssl(ns) + tau(k) = (ki_ssl(ns) + kabs_chl(ns,k)) * dzk(k) + w0 (k) = ki_ssl(ns)/(ki_ssl(ns) + kabs_chl(ns,k)) * wi_ssl(ns) + g (k) = gi_ssl(ns) ! dl k = kii + 1 ! scale dz for dl relative to 4 even-layer-thickness 1.5m case fs = p25/rnilyr - tau(k) = (ki_dl(ns) + kabs_chl(ns,k)) *dzk(k)*fs - w0(k) = ki_dl(ns)/(ki_dl(ns) + kabs_chl(ns,k)) *wi_dl(ns) - g(k) = gi_dl(ns) + tau(k) = (ki_dl(ns) + kabs_chl(ns,k)) * dzk(k) * fs + w0 (k) = ki_dl(ns)/(ki_dl(ns) + kabs_chl(ns,k)) * wi_dl(ns) + g (k) = gi_dl(ns) ! int above lowest layer if (kii+2 <= klev-1) then do k = kii+2, klev-1 - tau(k) = (ki_int(ns) + kabs_chl(ns,k))*dzk(k) - w0(k) = ki_int(ns)/(ki_int(ns) + kabs_chl(ns,k)) *wi_int(ns) - g(k) = gi_int(ns) + tau(k) = (ki_int(ns) + kabs_chl(ns,k)) * dzk(k) + w0 (k) = ki_int(ns)/(ki_int(ns) + kabs_chl(ns,k)) * wi_int(ns) + g (k) = gi_int(ns) enddo endif ! lowest layer k = klev ! add algae to lowest sea ice layer, visible only: kabs = ki_int(ns)*(c1-wi_int(ns)) - if( ns == 1 ) then + 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) endif - sig = ki_int(ns)*wi_int(ns) - tau(k) = (kabs+sig)*dzk(k) - w0(k) = (sig/(sig+kabs)) - g(k) = gi_int(ns) + sig = ki_int(ns) * wi_int(ns) + tau(k) = (kabs+sig) * dzk(k) + w0 (k) = sig/(sig+kabs) + g (k) = gi_int(ns) ! aerosol in sea ice if (tr_zaero .and. dEdd_algae) then do k = kii, klev @@ -2534,151 +2522,158 @@ subroutine compute_dEdd_3bd( & do na=1,4*n_aero,4 if (modal_aero) then if (na==1) then ! interstitial BC - taer = taer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcexs(k)) - waer = waer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k)) - gaer = gaer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k))*gaer_bc_3bd(ns,k_bcexs(k)) + taer = taer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcexs(k)) + waer = waer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) + gaer = gaer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) & + * gaer_bc_3bd(ns,k_bcexs(k)) elseif (na==5) then ! within-ice BC - taer = taer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcins(k))* & - bcenh_3bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k)) - gaer = gaer + & - aero_mp(na+2)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k))*gaer_bc_3bd(ns,k_bcins(k)) + taer = taer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * bcenh_3bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) + gaer = gaer + aero_mp(na+2) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) & + * gaer_bc_3bd(ns,k_bcins(k)) else ! other species (dust) - taer = taer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) + taer = taer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) endif else ! bulk - taer = taer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na+2)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) + taer = taer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na+2) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) endif ! modal_aero enddo ! na - - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & + w0 (k) = (w0(k)*tau(k) + waer) / & (tau(k) + taer) tau(k) = tau(k) + taer do k = kii+1, klev taer = c0 waer = c0 gaer = c0 - do na=1,4*n_aero,4 + do na = 1, 4*n_aero, 4 if (modal_aero) then if (na==1) then ! interstitial BC - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcexs(k)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcexs(k))* & - waer_bc_3bd(ns,k_bcexs(k))*gaer_bc_3bd(ns,k_bcexs(k)) + taer = taer + (aero_mp(na+3)/rnilyr) & !echmod: BUG should be *rnilyr + * kaer_bc_3bd(ns,k_bcexs(k)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_3bd(ns,k_bcexs(k)) & + * waer_bc_3bd(ns,k_bcexs(k)) & + * gaer_bc_3bd(ns,k_bcexs(k)) elseif (na==5) then ! within-ice BC - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcins(k))* & - bcenh_3bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_3bd(ns,k_bcins(k))* & - waer_bc_3bd(ns,k_bcins(k))*gaer_bc_3bd(ns,k_bcins(k)) + taer = taer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * bcenh_3bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_3bd(ns,k_bcins(k)) & + * waer_bc_3bd(ns,k_bcins(k)) & + * gaer_bc_3bd(ns,k_bcins(k)) else ! other species (dust) - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) + taer = taer + (aero_mp(na+3)/rnilyr) & + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) endif else ! bulk - taer = taer + & - (aero_mp(na+3)*rnilyr)*kaer_3bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+3)*rnilyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+3)*rnilyr)*kaer_3bd(ns,(1+(na-1)/4))* & - waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4)) + taer = taer + (aero_mp(na+3)*rnilyr) & !echmod NOTE * rnilyr, not / + * kaer_3bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+3)*rnilyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+3)*rnilyr) & + * kaer_3bd(ns,(1+(na-1)/4)) & + * waer_3bd(ns,(1+(na-1)/4)) & + * gaer_3bd(ns,(1+(na-1)/4)) endif ! modal_aero enddo ! na - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & + w0 (k) = (w0(k)*tau(k) + waer) / & (tau(k) + taer) tau(k) = tau(k) + taer enddo ! k - endif ! tr_aero + endif ! tr_aero else ! srftyp == 2 ! sea ice layers under ponds k = kii tau(k) = ki_p_ssl(ns)*dzk(k) - w0(k) = wi_p_ssl(ns) - g(k) = gi_p_ssl(ns) + w0 (k) = wi_p_ssl(ns) + g (k) = gi_p_ssl(ns) k = kii + 1 tau(k) = ki_p_int(ns)*dzk(k) - w0(k) = wi_p_int(ns) - g(k) = gi_p_int(ns) + w0 (k) = wi_p_int(ns) + g (k) = gi_p_int(ns) if (kii+2 <= klev) then do k = kii+2, klev tau(k) = ki_p_int(ns)*dzk(k) - w0(k) = wi_p_int(ns) - g(k) = gi_p_int(ns) + w0 (k) = wi_p_int(ns) + g (k) = gi_p_int(ns) enddo ! k endif ! adjust pond iops if pond depth within specified range if( hpmin <= hp .and. hp <= hp0 ) then k = kii - sig_i = ki_ssl(ns)*wi_ssl(ns) - sig_p = ki_p_ssl(ns)*wi_p_ssl(ns) - sig = sig_i + (sig_p-sig_i)*(hp/hp0) - kext = sig + ki_p_ssl(ns)*(c1-wi_p_ssl(ns)) + sig_i = ki_ssl (ns) * wi_ssl (ns) + sig_p = ki_p_ssl(ns) * wi_p_ssl(ns) + sig = sig_i + (sig_p-sig_i) * (hp/hp0) + kext = sig + ki_p_ssl(ns) * (c1-wi_p_ssl(ns)) tau(k) = kext*dzk(k) - w0(k) = sig/kext - g(k) = gi_p_int(ns) + w0 (k) = sig/kext + g (k) = gi_p_int(ns) k = kii + 1 ! scale dz for dl relative to 4 even-layer-thickness 1.5m case fs = p25/rnilyr - sig_i = ki_dl(ns)*wi_dl(ns)*fs - sig_p = ki_p_int(ns)*wi_p_int(ns) - sig = sig_i + (sig_p-sig_i)*(hp/hp0) - kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns)) + sig_i = ki_dl (ns) * wi_dl (ns) * fs + sig_p = ki_p_int(ns) * wi_p_int(ns) + sig = sig_i + (sig_p-sig_i) * (hp/hp0) + kext = sig + ki_p_int(ns) * (c1-wi_p_int(ns)) tau(k) = kext*dzk(k) - w0(k) = sig/kext - g(k) = gi_p_int(ns) + w0 (k) = sig/kext + g (k) = gi_p_int(ns) if (kii+2 <= klev) then do k = kii+2, klev - sig_i = ki_int(ns)*wi_int(ns) - sig_p = ki_p_int(ns)*wi_p_int(ns) - sig = sig_i + (sig_p-sig_i)*(hp/hp0) - kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns)) + sig_i = ki_int (ns) * wi_int (ns) + sig_p = ki_p_int(ns) * wi_p_int(ns) + sig = sig_i + (sig_p-sig_i) * (hp/hp0) + kext = sig + ki_p_int(ns) * (c1-wi_p_int(ns)) tau(k) = kext*dzk(k) - w0(k) = sig/kext - g(k) = gi_p_int(ns) + w0 (k) = sig/kext + g (k) = gi_p_int(ns) enddo ! k endif endif ! small pond depth transition to bare sea ice @@ -2697,8 +2692,8 @@ subroutine compute_dEdd_3bd( & ! 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, & + call solution_dEdd ( & + coszen, srftyp, klev, klevp, & tau, w0, g, albodr, albodf, & trndir, trntdr, trndif, rupdir, rupdif, & rdndif) @@ -2717,7 +2712,7 @@ subroutine compute_dEdd_3bd( & do k = 0, klevp ! interface scattering - refk = c1/(c1 - rdndif(k)*rupdif(k)) + refk = c1/(c1 - rdndif(k)*rupdif(k)) ! dir tran ref from below times interface scattering, plus diff ! tran and ref from below times interface scattering ! fdirup(k) = (trndir(k)*rupdir(k) + & @@ -2745,14 +2740,11 @@ subroutine compute_dEdd_3bd( & ! 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) - + 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 @@ -2769,9 +2761,9 @@ subroutine compute_dEdd_3bd( & fthruvdf = fthruvdf + dfdif(klevp)*swdf ! if snow covered ice, set snow internal absorption; else, Sabs=0 - if( srftyp == 1 ) then + if (srftyp == 1) then ki = 0 - do k=1,nslyr + do k = 1, nslyr ! skip snow SSL, since SSL absorption included in the surface ! absorption fsfc above km = k @@ -2785,14 +2777,14 @@ subroutine compute_dEdd_3bd( & ! complex indexing to insure proper absorptions for sea ice ki = 0 - do k=nslyr+2,nslyr+1+nilyr + do k = nslyr+2, nslyr+1+nilyr ! for bare ice, DL absorption for sea ice layer 1 km = k kp = km + 1 ! modify for top sea ice layer for snow over sea ice - if( srftyp == 1 ) then + if (srftyp == 1) then ! must add SSL and DL absorption for sea ice layer 1 - if( k == nslyr+2 ) then + if (k == nslyr+2) then km = k - 1 kp = km + 2 endif @@ -2833,14 +2825,14 @@ subroutine compute_dEdd_3bd( & fthruidf = fthruidf + dfdif(klevp)*swdf*wghtns(ns) ! if snow covered ice, set snow internal absorption; else, Sabs=0 - if( srftyp == 1 ) then + if (srftyp == 1) then ki = 0 - do k=1,nslyr + do k = 1, nslyr ! skip snow SSL, since SSL absorption included in the surface ! absorption fsfc above - km = k - kp = km + 1 - ki = ki + 1 + km = k + kp = km + 1 + ki = ki + 1 Sabs(ki) = Sabs(ki) & + (dfdir(km)*swdr + dfdif(km)*swdf & - (dfdir(kp)*swdr + dfdif(kp)*swdf)) & @@ -2850,14 +2842,14 @@ subroutine compute_dEdd_3bd( & ! complex indexing to insure proper absorptions for sea ice ki = 0 - do k=nslyr+2,nslyr+1+nilyr + do k = nslyr+2, nslyr+1+nilyr ! for bare ice, DL absorption for sea ice layer 1 km = k kp = km + 1 ! modify for top sea ice layer for snow over sea ice - if( srftyp == 1 ) then + if (srftyp == 1) then ! must add SSL and DL absorption for sea ice layer 1 - if( k == nslyr+2 ) then + if (k == nslyr+2) then km = k - 1 kp = km + 2 endif @@ -2868,15 +2860,13 @@ subroutine compute_dEdd_3bd( & - (dfdir(kp)*swdr + dfdif(kp)*swdf)) & * wghtns(ns) enddo ! k + endif ! ns + enddo ! ns: end spectral loop - endif ! ns = 1, ns > 1 - - enddo ! ns: end spectral loop - - alvdr = avdr - alvdf = avdf - alidr = aidr - alidf = aidf + alvdr = avdr + alvdf = avdf + alidr = aidr + alidf = aidf ! accumulate fluxes over bare sea ice fswsfc = fswsfc + fsfc *fi @@ -2889,17 +2879,14 @@ subroutine compute_dEdd_3bd( & do k = 1, nslyr Sswabs(k) = Sswabs(k) + Sabs(k)*fi - enddo ! k + enddo do k = 1, nilyr Iswabs(k) = Iswabs(k) + Iabs(k)*fi - ! 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 + enddo + fswpenl(nilyr+1) = fswpenl(nilyr+1) + fthrul(nilyr+1)*fi #ifdef UNDEPRECATE_0LAYER !---------------------------------------------------------------- @@ -2929,8 +2916,8 @@ end subroutine compute_dEdd_3bd ! ! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version - subroutine solution_dEdd & - (coszen, srftyp, klev, klevp, & + subroutine solution_dEdd ( & + coszen, srftyp, klev, klevp, & tau, w0, g, albodr, albodf, & trndir, trntdr, trndif, rupdir, rupdif, & rdndif) @@ -3289,19 +3276,15 @@ subroutine solution_dEdd & ! the fresnel (refractive) layer, always taken to be above ! the present layer k (i.e. be the top interface): - rintfc = c1 / (c1-Rf_dif_b*rdif_a(k)) - tdir(k) = Tf_dir_a*tdir(k) + & - Tf_dir_a*rdir(k) * & - Rf_dif_b*rintfc*tdif_a(k) - rdir(k) = Rf_dir_a + & - Tf_dir_a*rdir(k) * & - rintfc*Tf_dif_b - rdif_a(k) = Rf_dif_a + & - Tf_dif_a*rdif_a(k) * & - rintfc*Tf_dif_b - rdif_b(k) = rdif_b(k) + & - tdif_b(k)*Rf_dif_b * & - rintfc*tdif_a(k) + rintfc = c1 / (c1-Rf_dif_b*rdif_a(k)) + tdir (k) = Tf_dir_a*tdir(k) & + + Tf_dir_a*rdir(k) * Rf_dif_b*rintfc*tdif_a(k) + rdir (k) = Rf_dir_a & + + Tf_dir_a*rdir (k) * rintfc*Tf_dif_b + rdif_a(k) = Rf_dif_a & + + Tf_dif_a*rdif_a(k) * rintfc*Tf_dif_b + rdif_b(k) = rdif_b(k) & + + tdif_b(k)*Rf_dif_b * rintfc*tdif_a(k) tdif_a(k) = tdif_a(k)*rintfc*Tf_dif_a tdif_b(k) = tdif_b(k)*rintfc*Tf_dif_b @@ -3337,13 +3320,13 @@ subroutine solution_dEdd & ! \\\\\\\ ocean \\\\\\\ trndir(k+1) = trndir(k)*trnlay(k) - refkm1 = c1/(c1 - rdndif(k)*rdif_a(k)) - tdrrdir = trndir(k)*rdir(k) - tdndif = trntdr(k) - trndir(k) - trntdr(k+1) = trndir(k)*tdir(k) + & - (tdndif + tdrrdir*rdndif(k))*refkm1*tdif_a(k) - rdndif(k+1) = rdif_b(k) + & - (tdif_b(k)*rdndif(k)*refkm1*tdif_a(k)) + refkm1 = c1/(c1 - rdndif(k)*rdif_a(k)) + tdrrdir = trndir(k)*rdir(k) + tdndif = trntdr(k) - trndir(k) + trntdr(k+1) = trndir(k)*tdir(k) & + + (tdndif + tdrrdir*rdndif(k))*refkm1*tdif_a(k) + rdndif(k+1) = rdif_b(k) & + + (tdif_b(k)*rdndif(k)*refkm1*tdif_a(k)) trndif(k+1) = trndif(k)*refkm1*tdif_a(k) enddo ! k end main level loop @@ -3365,7 +3348,7 @@ subroutine solution_dEdd & do k=klev,0,-1 ! interface scattering - refkp1 = c1/( c1 - rdif_b(k)*rupdif(k+1)) + refkp1 = c1/( c1 - rdif_b(k)*rupdif(k+1)) ! dir from top layer plus exp tran ref from lower layer, interface ! scattered and tran thru top layer from below, plus diff tran ref ! from lower layer with interface scattering tran thru top from below @@ -3387,7 +3370,7 @@ end subroutine solution_dEdd ! author: Bruce P. Briegleb, NCAR ! 2013: E Hunke merged with NCAR version - subroutine shortwave_dEdd_set_snow(R_snw, & + subroutine shortwave_dEdd_set_snow(R_snw, & dT_mlt, rsnw_mlt, & aice, vsno, & Tsfc, fs, & @@ -3429,8 +3412,8 @@ subroutine shortwave_dEdd_set_snow(R_snw, & real (kind=dbl_kind), parameter :: & ! units for the following are 1.e-6 m (micro-meters) - rsnw_nonmelt = 500._dbl_kind, & ! nonmelt snow grain radius - rsnw_sig = 250._dbl_kind ! assumed sigma for snow grain radius + rsnw_nonmelt = 500._dbl_kind, & ! nonmelt snow grain radius + rsnw_sig = 250._dbl_kind ! assumed sigma for snow grain radius character(len=*),parameter :: subname='(shortwave_dEdd_set_snow)' @@ -3518,19 +3501,15 @@ subroutine shortwave_dEdd_set_pond(Tsfc, & end subroutine shortwave_dEdd_set_pond -! End Delta-Eddington shortwave method - !======================================================================= ! ! authors Nicole Jeffery, LANL - subroutine compute_shortwave_trcr( & - bgcN, zaero, & - trcrn_bgcsw, & - sw_grid, hin, & - hbri, & - i_grid, & - skl_bgc, z_tracers ) + subroutine compute_shortwave_trcr(bgcN, zaero, & + trcrn_bgcsw, sw_grid, & + hin, hbri, & + i_grid, skl_bgc, & + z_tracers) real (kind=dbl_kind), dimension (:), intent(in) :: & bgcN , & ! Nit tracer @@ -3540,16 +3519,16 @@ subroutine compute_shortwave_trcr( & trcrn_bgcsw ! ice on shortwave grid tracers real (kind=dbl_kind), dimension (:), intent(in) :: & - sw_grid , & ! - i_grid ! CICE bio grid + sw_grid , & ! + i_grid ! CICE bio grid real(kind=dbl_kind), intent(in) :: & - hin , & ! CICE ice thickness - hbri ! brine height + hin , & ! CICE ice thickness + hbri ! brine height logical (kind=log_kind), intent(in) :: & - skl_bgc , & ! skeletal layer bgc - z_tracers ! zbgc + skl_bgc , & ! skeletal layer bgc + z_tracers ! zbgc ! local variables @@ -3593,8 +3572,8 @@ subroutine compute_shortwave_trcr( & do k = 1, nblyr+1 do n = 1, n_algae - trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) + & - R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1 + k-1) + trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) & + + R_chl2N(n) * F_abs_chl(n) * bgcN(nt_bgc_N(n)-nt_bgc_N(1) + k) enddo ! n enddo ! k @@ -3668,6 +3647,7 @@ subroutine compute_shortwave_trcr( & enddo endif + end subroutine compute_shortwave_trcr !======================================================================= @@ -3759,21 +3739,21 @@ subroutine icepack_prep_radiation(aice, aicen, & ! Scale absorbed solar radiation for change in net shortwave !----------------------------------------------------------------- - fswsfcn(n) = scale_factor*fswsfcn (n) - fswintn(n) = scale_factor*fswintn (n) - fswthrun(n) = scale_factor*fswthrun(n) - if (present(fswthrun_vdr)) fswthrun_vdr(n) = scale_factor*fswthrun_vdr(n) - if (present(fswthrun_vdf)) fswthrun_vdf(n) = scale_factor*fswthrun_vdf(n) - if (present(fswthrun_idr)) fswthrun_idr(n) = scale_factor*fswthrun_idr(n) - if (present(fswthrun_idf)) fswthrun_idf(n) = scale_factor*fswthrun_idf(n) + fswsfcn(n) = scale_factor * fswsfcn (n) + fswintn(n) = scale_factor * fswintn (n) + fswthrun(n) = scale_factor * fswthrun(n) + if (present(fswthrun_vdr)) fswthrun_vdr(n) = scale_factor * fswthrun_vdr(n) + if (present(fswthrun_vdf)) fswthrun_vdf(n) = scale_factor * fswthrun_vdf(n) + if (present(fswthrun_idr)) fswthrun_idr(n) = scale_factor * fswthrun_idr(n) + if (present(fswthrun_idf)) fswthrun_idf(n) = scale_factor * fswthrun_idf(n) do k = 1,nilyr+1 - fswpenln(k,n) = scale_factor*fswpenln(k,n) + fswpenln(k,n) = scale_factor * fswpenln(k,n) enddo !k do k=1,nslyr - Sswabsn(k,n) = scale_factor*Sswabsn(k,n) + Sswabsn (k,n) = scale_factor * Sswabsn (k,n) enddo do k=1,nilyr - Iswabsn(k,n) = scale_factor*Iswabsn(k,n) + Iswabsn (k,n) = scale_factor * Iswabsn (k,n) enddo endif @@ -4081,16 +4061,16 @@ subroutine icepack_step_radiation (dt, & ! Initialize for safety do n = 1, ncat - alvdrn(n) = c0 - alidrn(n) = c0 - alvdfn(n) = c0 - alidfn(n) = c0 - fswsfcn(n) = c0 - fswintn(n) = c0 + alvdrn (n) = c0 + alidrn (n) = c0 + alvdfn (n) = c0 + alidfn (n) = c0 + fswsfcn (n) = c0 + fswintn (n) = c0 fswthrun(n) = c0 enddo ! ncat - Iswabsn(:,:) = c0 - Sswabsn(:,:) = c0 + Iswabsn (:,:) = c0 + Sswabsn (:,:) = c0 endif ! calc_Tsfc @@ -4128,7 +4108,7 @@ real(kind=dbl_kind) function n(uu,et) real(kind=dbl_kind), intent(in) :: uu, et - n = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et) + n = ((uu+c1)*(uu+c1)/et) - ((uu-c1)*(uu-c1)*et) end function n @@ -4212,61 +4192,65 @@ end function asys ! for a unified treatment of cryospheric surfaces, The Cryosphere, 13, ! 2325-2343, https://doi.org/10.5194/tc-13-2325-2019, 2019. - subroutine compute_dEdd_5bd (klev, klevp, & - zbio, fnidr, coszen, & - swvdr, swvdf, swidr, swidf, srftyp, & - hs, rhosnw, rsnw, hi, hp, & - fi, aero_mp, alvdr, alvdf, & - alidr, alidf, & - fswsfc, fswint, & - fswthru, Sswabs, & - Iswabs, fswpenl ) + subroutine compute_dEdd_5bd( & + klev, klevp, zbio, fnidr, coszen, & + swvdr, swvdf, swidr, swidf, srftyp, & + hs, rhosnw, rsnw, hi, hp, & + fi, aero_mp, alvdr, alvdf, & + alidr, alidf, fswsfc, fswint, fswthru, & + fswthru_vdr, fswthru_vdf, & + fswthru_idr, fswthru_idf, & + Sswabs, Iswabs, fswpenl ) integer (kind=int_kind), intent(in) :: & - klev , & ! number of radiation layers - 1 - klevp ! number of radiation interfaces - 1 - ! (0 layer is included also) - - ! dEdd tuning parameters, set in namelist + klev , & ! number of radiation layers - 1 + klevp ! number of radiation interfaces - 1 + ! (0 layer is included also) real (kind=dbl_kind), intent(in) :: & - fnidr , & ! fraction of direct to total down flux in nir - coszen , & ! cosine solar zenith angle - swvdr , & ! shortwave down at surface, visible, direct (W/m^2) - 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) + fnidr , & ! fraction of direct to total down flux in nir + coszen, & ! cosine solar zenith angle + swvdr , & ! shortwave down at surface, visible, direct (W/m^2) + 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) :: & - srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) + srftyp ! surface type over ice: (0=air, 1=snow, 2=pond) real (kind=dbl_kind), intent(in) :: & - hs ! snow thickness (m) + hs ! snow thickness (m) real (kind=dbl_kind), dimension (:), intent(in) :: & - rhosnw , & ! snow density in snow layer (kg/m3) - rsnw , & ! snow grain radius in snow layer (m) - zbio , & ! zaerosol + chla shortwave tracers kg/m^3 - aero_mp ! aerosol mass path in kg/m2 + rhosnw, & ! snow density in snow layer (kg/m3) + rsnw , & ! snow grain radius in snow layer (m) + zbio , & ! zaerosol + chla shortwave tracers kg/m^3 + aero_mp ! aerosol mass path in kg/m2 real (kind=dbl_kind), intent(in) :: & - hi , & ! ice thickness (m) - hp , & ! pond depth (m) - fi ! snow/bare ice fractional coverage (0 to 1) + 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) - 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) + 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) + + real (kind=dbl_kind), intent(inout) :: & + fswthru_vdr, & ! vis dir SW through snow/bare ice/ponded ice into ocean (W m-2) + 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) + 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) !----------------------------------------------------------------------- ! Set up optical property profiles, based on snow, sea ice and ponded @@ -4391,7 +4375,11 @@ subroutine compute_dEdd_5bd (klev, klevp, & 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) + fthru , & ! shortwave through snow/bare ice/ponded ice to ocean (W/m^2) + fthruvdr, & ! vis dir shortwave through snow/bare ice/ponded ice to ocean (W/m^2) + fthruvdf, & ! vis dif shortwave through snow/bare ice/ponded ice to ocean (W/m^2) + 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) :: & Sabs ! shortwave absorbed in snow layer (W m-2) @@ -4403,11 +4391,11 @@ subroutine compute_dEdd_5bd (klev, klevp, & fthrul ! shortwave through to ice layers (W m-2) real (kind=dbl_kind), parameter :: & - cp67 = 0.67_dbl_kind , & ! nir band weight parameter - cp33 = 0.33_dbl_kind , & ! nir band weight parameter - cp78 = 0.78_dbl_kind , & ! nir band weight parameter - cp22 = 0.22_dbl_kind , & ! nir band weight parameter - cp01 = 0.01_dbl_kind ! for ocean visible albedo + cp67 = 0.67_dbl_kind, & ! nir band weight parameter + cp33 = 0.33_dbl_kind, & ! nir band weight parameter + cp78 = 0.78_dbl_kind, & ! nir band weight parameter + cp22 = 0.22_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 @@ -4437,35 +4425,35 @@ subroutine compute_dEdd_5bd (klev, klevp, & ws , & ! Snow single scattering albedo gs ! Snow asymmetry parameter - real (kind=dbl_kind), dimension(nslyr) :: & - frsnw ! snow grain radius in snow layer * adjustment factor (m) +! real (kind=dbl_kind), dimension(nslyr) :: & +! frsnw ! snow grain radius in snow layer * adjustment factor (m) real (kind=dbl_kind), dimension(0:klev) :: & - dzk ! layer thickness + dzk ! layer thickness real (kind=dbl_kind) :: & - dz , & ! snow, sea ice or pond water layer thickness - 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 - ! conserve computed albedo for constant physical depth of sea - ! ice when the number of sea ice layers vary - real (kind=dbl_kind) :: & - sig , & ! scattering coefficient for tuning - kabs , & ! absorption coefficient for tuning - sigp ! modified scattering coefficient for tuning + dz , & ! snow, sea ice or pond water layer thickness + 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 + ! conserve computed albedo for constant physical depth of sea + ! ice when the number of sea ice layers vary + real (kind=dbl_kind) :: & + sig , & ! scattering coefficient for tuning + kabs , & ! absorption coefficient for tuning + sigp ! modified scattering coefficient for tuning real (kind=dbl_kind) :: & - albodr , & ! spectral ocean albedo to direct rad - albodf ! spectral ocean albedo to diffuse rad + 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 real (kind=dbl_kind) :: & - sig_i , & ! ice scattering coefficient (/m) - sig_p , & ! pond scattering coefficient (/m) - kext ! weighted extinction coefficient (/m) + sig_i , & ! ice scattering coefficient (/m) + sig_p , & ! pond scattering coefficient (/m) + kext ! weighted extinction coefficient (/m) ! aerosol optical properties from Mark Flanner, 26 June 2008 ! order assumed: hydrophobic black carbon, hydrophilic black carbon, @@ -4476,30 +4464,30 @@ subroutine compute_dEdd_5bd (klev, klevp, & ! and 1.19-5.0 micron in wavelength) integer (kind=int_kind) :: & - na , n ! aerosol index + na , n ! aerosol index real (kind=dbl_kind) :: & - taer , & ! total aerosol extinction optical depth - waer , & ! total aerosol single scatter albedo - gaer , & ! total aerosol asymmetry parameter - swdr , & ! shortwave down at surface, direct (W/m^2) - swdf , & ! shortwave down at surface, diffuse (W/m^2) - rnilyr , & ! real(nilyr) - rnslyr , & ! real(nslyr) - rns , & ! real(ns) - tmp_0, tmp_ks, tmp_kl ! temp variables + taer , & ! total aerosol extinction optical depth + waer , & ! total aerosol single scatter albedo + gaer , & ! total aerosol asymmetry parameter + swdr , & ! shortwave down at surface, direct (W/m^2) + swdf , & ! shortwave down at surface, diffuse (W/m^2) + rnilyr , & ! 1/real(nilyr) + rnslyr , & ! 1/real(nslyr) + rns , & ! real(ns) + tmp_0, tmp_ks, tmp_kl ! temporary variables integer(kind=int_kind), dimension(0:klev) :: & - k_bcini , & - k_bcins , & - k_bcexs + k_bcini , & ! index + k_bcins , & ! = 2 hardwired + k_bcexs ! = 2 hardwired real(kind=dbl_kind):: & - tmp_gs, tmp1 ! temp variables + tmp_gs, tmp1 ! temporary variables real (kind=dbl_kind), parameter :: & - fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max - fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min + fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max + fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min ! tuning parameters ! ice and pond scat coeff fractional change for +- one-sigma in albedo fp_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for + stn dev in alb @@ -4507,32 +4495,29 @@ subroutine compute_dEdd_5bd (klev, klevp, & fp_pnd = 2.00_dbl_kind, & ! ponded ice fraction of scat coeff for + stn dev in alb fm_pnd = 0.50_dbl_kind ! ponded ice fraction of scat coeff for - stn dev in alb - real (kind=dbl_kind), parameter :: & !chla-specific absorption coefficient - kchl_tab = 0.01 !0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1 - ! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014 - ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011) - ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014) - ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting) - !chlorophyll mass extinction cross section (m^2/mg chla) - - character(len=char_len_long) :: & - warning ! warning message + real (kind=dbl_kind), parameter :: & ! chla-specific absorption coefficient + kchl_tab = p01 ! 0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1 + ! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014 + ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011) + ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014) + ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting) + ! chlorophyll mass extinction cross section (m^2/mg chla) real (kind=dbl_kind), dimension (nspint_5bd) :: & - wghtns_5bd_dfs, & ! spectral weights for diffuse incident - wghtns_5bd_drc ! spectral weights for direct incident + wghtns_5bd_dfs , & ! spectral weights for diffuse incident + wghtns_5bd_drc ! spectral weights for direct incident ! snow grain single-scattering properties for ! direct (drc) and diffuse (dfs) shortwave incidents ! local variable names, point to table data ! TODO use variable names in ice_shortwave_data directly real (kind=dbl_kind), pointer, dimension(:,:) :: & ! Model SNICAR snow SSP - asm_prm_ice_drc , & ! snow asymmetry factor (cos(theta)) - asm_prm_ice_dfs , & ! snow asymmetry factor (cos(theta)) - ss_alb_ice_drc , & ! snow single scatter albedo (fraction) - ss_alb_ice_dfs , & ! snow single scatter albedo (fraction) - ext_cff_mss_ice_drc , & ! snow mass extinction cross section (m2/kg) - ext_cff_mss_ice_dfs ! snow mass extinction cross section (m2/kg) + asm_prm_ice_drc , & ! snow asymmetry factor (cos(theta)) + asm_prm_ice_dfs , & ! snow asymmetry factor (cos(theta)) + ss_alb_ice_drc , & ! snow single scatter albedo (fraction) + ss_alb_ice_dfs , & ! snow single scatter albedo (fraction) + ext_cff_mss_ice_drc , & ! snow mass extinction cross section (m2/kg) + ext_cff_mss_ice_dfs ! snow mass extinction cross section (m2/kg) ! FUTURE-WORK: update 5-band sea ice iops when avalible real (kind=dbl_kind), dimension (nspint_5bd) :: & ! for ice only @@ -4555,7 +4540,7 @@ subroutine compute_dEdd_5bd (klev, klevp, & ! index integer (kind=int_kind) :: & - nsky !sky = 1 (2) for direct (diffuse) downward SW incident + nsky ! sky = 1 (2) for direct (diffuse) downward SW incident ! temporary variables used to assign variables for direct/diffuse incident ! based on snicar 5 band IOPs @@ -4565,7 +4550,7 @@ subroutine compute_dEdd_5bd (klev, klevp, & rupdir_snicar , & ! reflectivity to direct radiation for layers below rupdif_snicar ! reflectivity to diffuse radiation for layers above - ! solar zenith angle parameterizations + ! solar zenith angle parameters real (kind=dbl_kind), parameter :: & sza_a0 = 0.085730_dbl_kind , & sza_a1 = -0.630883_dbl_kind , & @@ -4573,7 +4558,7 @@ subroutine compute_dEdd_5bd (klev, klevp, & sza_b0 = 1.467291_dbl_kind , & sza_b1 = -3.338043_dbl_kind , & sza_b2 = 6.807489_dbl_kind , & - mu_75 = 0.2588_dbl_kind ! cosine of 75 degree + mu_75 = 0.2588_dbl_kind ! cos(75 degrees) real (kind=dbl_kind) :: & sza_c1 , & ! parameter for high sza adjustment @@ -4581,7 +4566,7 @@ subroutine compute_dEdd_5bd (klev, klevp, & sza_factor , & ! parameter for high sza adjustment mu0 - character(len=*),parameter :: subname='(compute_dEdd_3bd)' + character(len=*),parameter :: subname='(compute_dEdd_5bd)' !----------------------------------------------------------------------- ! Initialize and tune bare ice/ponded ice iops @@ -4603,20 +4588,24 @@ subroutine compute_dEdd_5bd (klev, klevp, & kii = nslyr + 1 ! initialize albedos and fluxes to 0 - fthrul = c0 - Iabs = c0 + fthrul = c0 + Iabs = c0 kabs_chl_5bd(:,:) = c0 - tzaer_5bd(:,:) = c0 - wzaer_5bd(:,:) = c0 - gzaer_5bd(:,:) = c0 - - avdr = c0 - avdf = c0 - aidr = c0 - aidf = c0 - fsfc = c0 - fint = c0 - fthru = c0 + tzaer_5bd (:,:) = c0 + wzaer_5bd (:,:) = c0 + gzaer_5bd (:,:) = c0 + + avdr = c0 + avdf = c0 + aidr = c0 + aidf = c0 + fsfc = c0 + fint = c0 + fthru = c0 + fthruvdr = c0 + fthruvdf = c0 + fthruidr = c0 + fthruidf = c0 ! spectral weights - 3 bands ! this section of code is kept for future mearge between 5band and 3 band @@ -4626,29 +4615,26 @@ subroutine compute_dEdd_5bd (klev, klevp, & ! 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) +! wghtns(1) = c1 +! wghtns(2) = cp67 + (cp78-cp67)*(c1-fnidr) ! wghtns(3) = cp33 + (cp22-cp33)*(c1-fnidr) - !wghtns(3) = c1 - wghtns(2) +! wghtns(3) = c1 - wghtns(2) ! spectral weights - 5 bands ! direct beam incident ! add-local-variable - wghtns_5bd_drc(1) = 1._dbl_kind - wghtns_5bd_drc(2) = 0.49352158521175_dbl_kind!0.49352_dbl_kind!0.50_dbl_kind - wghtns_5bd_drc(3) = 0.18099494230665_dbl_kind!0.18100_dbl_kind!0.18_dbl_kind - wghtns_5bd_drc(4) = 0.12094898498813_dbl_kind!0.12095_dbl_kind!0.12_dbl_kind ! + wghtns_5bd_drc(1) = c1 + wghtns_5bd_drc(2) = 0.49352158521175_dbl_kind + wghtns_5bd_drc(3) = 0.18099494230665_dbl_kind + wghtns_5bd_drc(4) = 0.12094898498813_dbl_kind wghtns_5bd_drc(5) = c1-(wghtns_5bd_drc(2)+wghtns_5bd_drc(3)+wghtns_5bd_drc(4)) - !wghtns_5bd_drc(5) = 0.20453448749347_dbl_kind!0.20453_dbl_kind!0.20_dbl_kind ! ! diffuse incident - wghtns_5bd_dfs(1) = 1._dbl_kind - wghtns_5bd_dfs(2) = 0.58581507618433_dbl_kind!0.58582_dbl_kind!0.59_dbl_kind ! - wghtns_5bd_dfs(3) = 0.20156903770812_dbl_kind!0.20157_dbl_kind!0.20_dbl_kind ! - wghtns_5bd_dfs(4) = 0.10917889346386_dbl_kind!0.10918_dbl_kind!0.11_dbl_kind ! + wghtns_5bd_dfs(1) = c1 + wghtns_5bd_dfs(2) = 0.58581507618433_dbl_kind + wghtns_5bd_dfs(3) = 0.20156903770812_dbl_kind + wghtns_5bd_dfs(4) = 0.10917889346386_dbl_kind wghtns_5bd_dfs(5) = c1-(wghtns_5bd_dfs(2)+wghtns_5bd_dfs(3)+wghtns_5bd_dfs(4)) - !wghtns_5bd_dfs(5) = 0.10343699264369_dbl_kind!0.10343_dbl_kind!0.10_dbl_kind ! - do k = 1, nslyr !frsnw(k) = (fr_max*fnidr + fr_min*(c1-fnidr))*rsnw(k) @@ -4695,37 +4681,37 @@ subroutine compute_dEdd_5bd (klev, klevp, & ! Note: the albedo change becomes non-linear for R values > +1 or < -1 if( R_ice >= c0 ) then do ns = 1, nspint_5bd - sigp = ki_ssl_mn_5bd(ns)*wi_ssl_mn_5bd(ns)*(c1+fp_ice*R_ice) + sigp = ki_ssl_mn_5bd(ns)*wi_ssl_mn_5bd(ns)*(c1+fp_ice*R_ice) ki_ssl_5bd(ns) = sigp+ki_ssl_mn_5bd(ns)*(c1-wi_ssl_mn_5bd(ns)) wi_ssl_5bd(ns) = sigp/ki_ssl_5bd(ns) gi_ssl_5bd(ns) = gi_ssl_mn_5bd(ns) - sigp = ki_dl_mn_5bd(ns)*wi_dl_mn_5bd(ns)*(c1+fp_ice*R_ice) + sigp = ki_dl_mn_5bd(ns)*wi_dl_mn_5bd(ns)*(c1+fp_ice*R_ice) ki_dl_5bd(ns) = sigp+ki_dl_mn_5bd(ns)*(c1-wi_dl_mn_5bd(ns)) wi_dl_5bd(ns) = sigp/ki_dl_5bd(ns) gi_dl_5bd(ns) = gi_dl_mn_5bd(ns) - sigp = ki_int_mn_5bd(ns)*wi_int_mn_5bd(ns)*(c1+fp_ice*R_ice) + sigp = ki_int_mn_5bd(ns)*wi_int_mn_5bd(ns)*(c1+fp_ice*R_ice) ki_int_5bd(ns) = sigp+ki_int_mn_5bd(ns)*(c1-wi_int_mn_5bd(ns)) wi_int_5bd(ns) = sigp/ki_int_5bd(ns) gi_int_5bd(ns) = gi_int_mn_5bd(ns) enddo else !if( R_ice < c0 ) then do ns = 1, nspint_5bd - sigp = ki_ssl_mn_5bd(ns)*wi_ssl_mn_5bd(ns)*(c1+fm_ice*R_ice) - sigp = max(sigp, c0) + sigp = ki_ssl_mn_5bd(ns)*wi_ssl_mn_5bd(ns)*(c1+fm_ice*R_ice) + sigp = max(sigp, c0) ki_ssl_5bd(ns) = sigp+ki_ssl_mn_5bd(ns)*(c1-wi_ssl_mn_5bd(ns)) wi_ssl_5bd(ns) = sigp/ki_ssl_5bd(ns) gi_ssl_5bd(ns) = gi_ssl_mn_5bd(ns) - sigp = ki_dl_mn_5bd(ns)*wi_dl_mn_5bd(ns)*(c1+fm_ice*R_ice) - sigp = max(sigp, c0) + sigp = ki_dl_mn_5bd(ns)*wi_dl_mn_5bd(ns)*(c1+fm_ice*R_ice) + sigp = max(sigp, c0) ki_dl_5bd(ns) = sigp+ki_dl_mn_5bd(ns)*(c1-wi_dl_mn_5bd(ns)) wi_dl_5bd(ns) = sigp/ki_dl_5bd(ns) gi_dl_5bd(ns) = gi_dl_mn_5bd(ns) - sigp = ki_int_mn_5bd(ns)*wi_int_mn_5bd(ns)*(c1+fm_ice*R_ice) - sigp = max(sigp, c0) + sigp = ki_int_mn_5bd(ns)*wi_int_mn_5bd(ns)*(c1+fm_ice*R_ice) + sigp = max(sigp, c0) ki_int_5bd(ns) = sigp+ki_int_mn_5bd(ns)*(c1-wi_int_mn_5bd(ns)) wi_int_5bd(ns) = sigp/ki_int_5bd(ns) gi_int_5bd(ns) = gi_int_mn_5bd(ns) @@ -4742,750 +4728,737 @@ subroutine compute_dEdd_5bd (klev, klevp, & else k = klev kabs_chl_5bd(1,k) = kalg*(0.50_dbl_kind/dzk(k)) - !print *, 'aerosol, k, kabs_chl_5bd(1,k)', k, kabs_chl_5bd(1,k) endif -!mgf++ if (modal_aero) then - 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. - ! CICE adjusted snow grain radius rsnw to frsnw, while for - ! SNICAR there is no need, the tmp_gs is therefore calculated - ! differently from code in subroutine compute_dEdd - ksnow = k - min(k-1,0) - tmp_gs = rsnw(ksnow) ! use rsnw not frsnw - - ! get grain size index: - ! works for 25 < snw_rds < 1625 um: - if (tmp_gs < 125) then - tmp1 = tmp_gs/50 - k_bcini(k) = nint(tmp1) - elseif (tmp_gs < 175) then - k_bcini(k) = 2 - else - tmp1 = (tmp_gs/250)+2 - k_bcini(k) = nint(tmp1) - endif - else ! use the largest snow grain size for ice - k_bcini(k) = 8 - endif - ! Set index corresponding to BC effective radius. Here, - ! asssume constant BC effective radius of 100nm - ! (corresponding to index 2) - k_bcins(k) = 2 - k_bcexs(k) = 2 - - ! check bounds: - if (k_bcini(k) < 1) k_bcini(k) = 1 - if (k_bcini(k) > 8) k_bcini(k) = 8 - if (k_bcins(k) < 1) k_bcins(k) = 1 - if (k_bcins(k) > 10) k_bcins(k) = 10 - if (k_bcexs(k) < 1) k_bcexs(k) = 1 - if (k_bcexs(k) > 10) k_bcexs(k) = 10 - - ! print ice radius index: - ! write(warning,*) "MGFICE2:k, ice index= ",k, k_bcini(k) - ! call add_warning(warning) - enddo ! k - ! assign the aerosol index - - 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_5bd ! not weighted by aice - tzaer_5bd(ns,k) = tzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcexs(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer_5bd(ns,k) = wzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer_5bd(ns,k) = gzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))* & - gaer_bc_5bd(ns,k_bcexs(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k) - enddo ! nspint - enddo - elseif (n==2) then ! within-ice BC - do k = 0, klev - do ns = 1,nspint_5bd - tzaer_5bd(ns,k) = tzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcins(k)) * & - bcenh_5bd(ns,k_bcins(k),k_bcini(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer_5bd(ns,k) = wzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer_5bd(ns,k) = gzaer_5bd(ns,k)+kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))* & - gaer_bc_5bd(ns,k_bcins(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k) - enddo ! nspint - enddo - else ! dust - do k = 0, klev - do ns = 1,nspint_5bd ! not weighted by aice - tzaer_5bd(ns,k) = tzaer_5bd(ns,k)+kaer_5bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer_5bd(ns,k) = wzaer_5bd(ns,k)+kaer_5bd(ns,n)*waer_5bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer_5bd(ns,k) = gzaer_5bd(ns,k)+kaer_5bd(ns,n)*waer_5bd(ns,n)* & - gaer_5bd(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k) + 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. + ! CICE adjusted snow grain radius rsnw to frsnw, while for + ! SNICAR there is no need, the tmp_gs is therefore calculated + ! differently from code in subroutine compute_dEdd + ksnow = k - min(k-1,0) + tmp_gs = rsnw(ksnow) ! use rsnw not frsnw + + ! grain size index + ! works for 25 < snw_rds < 1625 um: + if (tmp_gs < 125._dbl_kind) then + tmp1 = tmp_gs/50._dbl_kind + k_bcini(k) = nint(tmp1) + elseif (tmp_gs < 175._dbl_kind) then + k_bcini(k) = 2 + else + tmp1 = (tmp_gs/250._dbl_kind) + c2 + k_bcini(k) = nint(tmp1) + endif + else ! use the largest snow grain size for ice + k_bcini(k) = 8 + endif + ! Set index corresponding to BC effective radius. Here, + ! asssume constant BC effective radius of 100nm + ! (corresponding to index 2) + k_bcins(k) = 2 ! hardwired + k_bcexs(k) = 2 ! hardwired + + ! check bounds + if (k_bcini(k) < 1) k_bcini(k) = 1 + if (k_bcini(k) > 8) k_bcini(k) = 8 +! if (k_bcins(k) < 1) k_bcins(k) = 1 ! hardwired +! if (k_bcins(k) > 10) k_bcins(k) = 10 +! if (k_bcexs(k) < 1) k_bcexs(k) = 1 +! if (k_bcexs(k) > 10) k_bcexs(k) = 10 + enddo ! k + + 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_5bd ! not weighted by aice + tzaer_5bd(ns,k) = tzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer_5bd(ns,k) = wzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer_5bd(ns,k) = gzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) & + * gaer_bc_5bd(ns,k_bcexs(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + enddo + enddo + elseif (n==2) then ! within-ice BC + do k = 0, klev + do ns = 1, nspint_5bd + tzaer_5bd(ns,k) = tzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcins(k)) & + * bcenh_5bd(ns,k_bcins(k),k_bcini(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer_5bd(ns,k) = wzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer_5bd(ns,k) = gzaer_5bd (ns,k) & + + kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) & + * gaer_bc_5bd(ns,k_bcins(k)) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + enddo + enddo + else ! dust + do k = 0, klev + do ns = 1, nspint_5bd ! not weighted by aice + tzaer_5bd(ns,k) = tzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer_5bd(ns,k) = wzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * waer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer_5bd(ns,k) = gzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * waer_5bd (ns,n) & + * gaer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) enddo ! nspint - enddo - endif !(n=1) - enddo ! n_zaero - endif ! tr_zaero and dEdd_algae - - else ! Bulk aerosol treatment - if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll - do n = 1,n_zaero ! multiply by aice? + enddo ! k + endif ! n + enddo ! n_zaero + endif ! tr_zaero and dEdd_algae + + else ! Bulk aerosol treatment + if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll + do n = 1, n_zaero ! multiply by aice? do k = 0, klev - do ns = 1,nspint_5bd ! not weighted by aice - tzaer_5bd(ns,k) = tzaer_5bd(ns,k)+kaer_5bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - wzaer_5bd(ns,k) = wzaer_5bd(ns,k)+kaer_5bd(ns,n)*waer_5bd(ns,n)* & - zbio(nlt_zaero_sw(n)+k)*dzk(k) - gzaer_5bd(ns,k) = gzaer_5bd(ns,k)+kaer_5bd(ns,n)*waer_5bd(ns,n)* & - gaer_5bd(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k) + do ns = 1, nspint_5bd ! not weighted by aice + tzaer_5bd(ns,k) = tzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + wzaer_5bd(ns,k) = wzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * waer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) + gzaer_5bd(ns,k) = gzaer_5bd(ns,k) & + + kaer_5bd (ns,n) & + * waer_5bd (ns,n) & + * gaer_5bd (ns,n) & + * zbio(nlt_zaero_sw(n)+k) * dzk(k) enddo ! nspint - enddo - enddo - endif !tr_zaero - - endif ! modal_aero - + enddo ! k + enddo ! n + endif ! tr_zaero + endif ! modal_aero !----------------------------------------------------------------------- + ! begin spectral loop do ns = 1, nspint_5bd - ! for snow-covered sea ice, comput 5 bands - !if( srftyp == 1 ) then + + ! for snow-covered sea ice, compute 5 bands + !if( srftyp == 1 ) then ! SNICAR-AD major changes ! 1. loop through 5bands: do ns = 1, nspint_5bd based on nsky ! 2. use snow grain size rsnow, not scaled frsnw ! 3. replace $IOPs_tab with $IOPs_snicar ! 4. replace wghtns with wghtns_5bd - do nsky = 1,2 ! loop for both direct beam and diffuse beam - if (nsky == 1) then ! direc incident - do k=0,nslyr - ! use top rsnw, rhosnw for snow ssl and rest of top layer - ksnow = k - min(k-1,0) - if (rsnw(ksnow) <= rsnw_snicar_min) then - ks = ext_cff_mss_ice_drc(ns,1) - ws = ss_alb_ice_drc(ns,1) - gs = asm_prm_ice_drc(ns,1) - elseif (rsnw(ksnow) >= rsnw_snicar_max) then - ks = ext_cff_mss_ice_drc(ns,nmbrad_snicar) - ws = ss_alb_ice_drc(ns,nmbrad_snicar) - gs = asm_prm_ice_drc(ns,nmbrad_snicar) - elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then - nr = ceiling(rsnw(ksnow)) - 30 + 1 + do nsky = 1, 2 ! loop for both direct beam and diffuse beam + if (nsky == 1) then ! direct incident + do k = 0, nslyr + ! use top rsnw, rhosnw for snow ssl and rest of top layer + ksnow = k - min(k-1,0) + if (rsnw(ksnow) <= rsnw_snicar_min) then + ks = ext_cff_mss_ice_drc(ns,1) + ws = ss_alb_ice_drc (ns,1) + gs = asm_prm_ice_drc (ns,1) + elseif (rsnw(ksnow) >= rsnw_snicar_max) then + ks = ext_cff_mss_ice_drc(ns,nmbrad_snicar) + ws = ss_alb_ice_drc (ns,nmbrad_snicar) + gs = asm_prm_ice_drc (ns,nmbrad_snicar) + elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then + ! radius = 30 --> nr = 1 in SNICAR table ! NOTE + nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30 ks = ext_cff_mss_ice_drc(ns,nr) - ws = ss_alb_ice_drc(ns,nr) - gs = asm_prm_ice_drc(ns,nr) - else ! linear interpolation in rsnw - ! radius = 30 --> nr = 1 in SNICAR table - nr = ceiling(rsnw(ksnow)) - 30 + 1 - delr = (rsnw(ksnow) - floor(rsnw(ksnow))) / & - (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow))) - ks = ext_cff_mss_ice_drc(ns,nr-1)*(delr) + & - ext_cff_mss_ice_drc(ns,nr)*(c1-delr) - ws = ss_alb_ice_drc(ns,nr-1)*(delr) + & - ss_alb_ice_drc(ns,nr)*(c1-delr) - gs = asm_prm_ice_drc(ns,nr-1)*(delr) + & - asm_prm_ice_drc(ns,nr)*(c1-delr) - endif - tau(k) = (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k))*dzk(k) - w0(k) = (ks*rhosnw(ksnow))/(ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k)) * ws - g(k) = gs - - !write(warning, *) "sky, k, tau, w0, g =", nsky, k, tau(k), w0(k), g(k) - !write(warning, *) "ns, ks, kabs_chl_5bd(ns,k), ", ns, ks, kabs_chl_5bd(ns,k) - !call add_warning(warning) - !print *, "rsnw(ksnow)", rsnw(ksnow) - !print *, "sky, k, tau, w0, g =", nsky, k, tau(k), w0(k), g(k) - !print *, "ns, ks, kabs_chl_5bd(ns,k), ",ns, ks, kabs_chl_5bd(ns,k) - + ws = ss_alb_ice_drc (ns,nr) + gs = asm_prm_ice_drc (ns,nr) + else ! linear interpolation in rsnw + nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30 + delr = (rsnw(ksnow) - floor(rsnw(ksnow))) & ! hardwired delta radius = 1 in table + / (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow))) ! denom always = 1? + ks = ext_cff_mss_ice_drc(ns,nr-1)* delr & ! echmod: BUG delr factors opposite 3bd case + + ext_cff_mss_ice_drc(ns,nr )*(c1-delr) ! echmod: and nsky=2 case below + ws = ss_alb_ice_drc (ns,nr-1)* delr & + + ss_alb_ice_drc (ns,nr )*(c1-delr) + gs = asm_prm_ice_drc (ns,nr-1)* delr & + + asm_prm_ice_drc (ns,nr )*(c1-delr) + endif + tau(k) = (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k))*dzk(k) + w0 (k) = ks*rhosnw(ksnow) / (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k)) * ws + g (k) = gs enddo ! k - elseif (nsky == 2) then ! diffuse incident - do k=0,nslyr - ! use top rsnw, rhosnw for snow ssl and rest of top layer - ksnow = k - min(k-1,0) - if (rsnw(ksnow) < rsnw_snicar_min) then - ks = ext_cff_mss_ice_dfs(ns,1) - ws = ss_alb_ice_dfs(ns,1) - gs = asm_prm_ice_dfs(ns,1) - elseif (rsnw(ksnow) > rsnw_snicar_max) then - ks = ext_cff_mss_ice_dfs(ns,nmbrad_snicar) - ws = ss_alb_ice_dfs(ns,nmbrad_snicar) - gs = asm_prm_ice_dfs(ns,nmbrad_snicar) - elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then - nr = ceiling(rsnw(ksnow)) - 30 + 1 - ks = ext_cff_mss_ice_dfs(ns,nr) - ws = ss_alb_ice_dfs(ns,nr) - gs = asm_prm_ice_dfs(ns,nr) + elseif (nsky == 2) then ! diffuse incident + do k = 0, nslyr + ! use top rsnw, rhosnw for snow ssl and rest of top layer + ksnow = k - min(k-1,0) + if (rsnw(ksnow) < rsnw_snicar_min) then + ks = ext_cff_mss_ice_dfs(ns,1) + ws = ss_alb_ice_dfs (ns,1) + gs = asm_prm_ice_dfs (ns,1) + elseif (rsnw(ksnow) > rsnw_snicar_max) then + ks = ext_cff_mss_ice_dfs(ns,nmbrad_snicar) + ws = ss_alb_ice_dfs (ns,nmbrad_snicar) + gs = asm_prm_ice_dfs (ns,nmbrad_snicar) + elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then + ! radius = 30 --> nr = 1 in SNICAR table ! NOTE + nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30 + ks = ext_cff_mss_ice_dfs(ns,nr) + ws = ss_alb_ice_dfs (ns,nr) + gs = asm_prm_ice_dfs (ns,nr) else ! linear interpolation in rsnw - ! radius = 30 --> nr = 1 in SNICAR table - nr = ceiling(rsnw(ksnow)) - 30 + 1 - delr = (rsnw(ksnow) - floor(rsnw(ksnow))) / & - (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow))) - ks = ext_cff_mss_ice_dfs(ns,nr-1)*(c1-delr) + & - ext_cff_mss_ice_dfs(ns,nr)*delr - ws = ss_alb_ice_dfs(ns,nr-1)*(c1-delr) + & - ss_alb_ice_dfs(ns,nr)*delr - gs = asm_prm_ice_dfs(ns,nr-1)*(c1-delr) + & - asm_prm_ice_dfs(ns,nr)*delr - endif - tau(k) = (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k))*dzk(k) - w0(k) = (ks*rhosnw(ksnow))/(ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k)) * ws - g(k) = gs - - !write(warning, *) "sky, k, tau, w0, g =", nsky, k, tau(k), w0(k), g(k) - !write(warning, *) "ns, ks, kabs_chl_5bd(ns,k), ", ns, ks, kabs_chl_5bd(ns,k) - !call add_warning(warning) + nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30 + delr = (rsnw(ksnow) - floor(rsnw(ksnow))) & ! hardwired delta radius = 1 in table + / (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow))) + ks = ext_cff_mss_ice_dfs(ns,nr-1)*(c1-delr) & + + ext_cff_mss_ice_dfs(ns,nr )* delr + ws = ss_alb_ice_dfs (ns,nr-1)*(c1-delr) & + + ss_alb_ice_dfs (ns,nr )* delr + gs = asm_prm_ice_dfs (ns,nr-1)*(c1-delr) & + + asm_prm_ice_dfs (ns,nr )* delr + endif + tau(k) = (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k))*dzk(k) + w0 (k) = ks*rhosnw(ksnow) / (ks*rhosnw(ksnow) + kabs_chl_5bd(ns,k)) * ws + g (k) = gs enddo ! k - endif ! end if nsky for snow IOPs assignment - !------------------------------------------------------------------------------ - - !aerosol in snow - if (tr_zaero .and. dEdd_algae) then - do k = 0,nslyr - g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / & - (w0(k)*tau(k) + wzaer_5bd(ns,k)) - w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / & - (tau(k) + tzaer_5bd(ns,k)) - tau(k) = tau(k) + tzaer_5bd(ns,k) - enddo - elseif (tr_aero) then + endif ! nsky for snow IOPs + + !------------------------------------------------------------------------------ + + ! aerosol in snow + if (tr_zaero .and. dEdd_algae) then + do k = 0,nslyr + g (k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / & + (w0(k)*tau(k) + wzaer_5bd(ns,k)) + w0 (k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / & + (tau(k) + tzaer_5bd(ns,k)) + tau(k) = tau(k) + tzaer_5bd(ns,k) + enddo + elseif (tr_aero) then k = 0 ! snow SSL taer = c0 waer = c0 gaer = c0 - do na=1,4*n_aero,4 -! mgf++ + do na = 1, 4*n_aero, 4 if (modal_aero) then - if (na == 1) then - !interstitial BC - taer = taer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k)) - waer = waer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k)) - gaer = gaer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))*gaer_bc_5bd(ns,k_bcexs(k)) - elseif (na == 5)then - !within-ice BC - taer = taer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k))* & - bcenh_5bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k)) - gaer = gaer + & - aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))*gaer_bc_5bd(ns,k_bcins(k)) - else - ! other species (dust) - taer = taer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif + if (na == 1) then ! interstitial BC + taer = taer + aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k)) + waer = waer + aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k)) & + *waer_bc_5bd(ns,k_bcexs(k)) + gaer = gaer + aero_mp(na)*kaer_bc_5bd(ns,k_bcexs(k)) & + *waer_bc_5bd(ns,k_bcexs(k)) & + *gaer_bc_5bd(ns,k_bcexs(k)) + elseif (na == 5) then ! within-ice BC + taer = taer + aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k)) & + * bcenh_5bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k)) & + *waer_bc_5bd(ns,k_bcins(k)) + gaer = gaer + aero_mp(na)*kaer_bc_5bd(ns,k_bcins(k)) & + *waer_bc_5bd(ns,k_bcins(k)) & + *gaer_bc_5bd(ns,k_bcins(k)) + else ! other species (dust) + taer = taer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) & + *waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) & + *waer_5bd(ns,(1+(na-1)/4)) & + *gaer_5bd(ns,(1+(na-1)/4)) + endif else - taer = taer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif !modal_aero -!mgf-- - enddo ! na + taer = taer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) & + *waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na)*kaer_5bd(ns,(1+(na-1)/4)) & + *waer_5bd(ns,(1+(na-1)/4)) & + *gaer_5bd(ns,(1+(na-1)/4)) + endif ! modal_aero + enddo ! na gaer = gaer/(waer+puny) waer = waer/(taer+puny) ! tcraig, again why does the above exist if taer=waer=gaer=0 below - do k=1,nslyr + do k = 1, nslyr taer = c0 waer = c0 gaer = c0 - do na=1,4*n_aero,4 + do na = 1, 4*n_aero, 4 if (modal_aero) then -!mgf++ - if (na==1) then - ! interstitial BC - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcexs(k)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))*gaer_bc_5bd(ns,k_bcexs(k)) - elseif (na==5) then - ! within-ice BC - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcins(k))*& - bcenh_5bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))*gaer_bc_5bd(ns,k_bcins(k)) - - else - ! other species (dust) - taer = taer + & - (aero_mp(na+1)/rnslyr)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+1)/rnslyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+1)/rnslyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif !(na==1) - + if (na==1) then ! interstitial BC + taer = taer + (aero_mp(na+1)/rnslyr) & !echmod: BUG should be *rnslyr + * kaer_bc_5bd(ns,k_bcexs(k)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) & + * gaer_bc_5bd(ns,k_bcexs(k)) + elseif (na==5) then ! within-ice BC + taer = taer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * bcenh_5bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) & + * gaer_bc_5bd(ns,k_bcins(k)) + else ! other species (dust) + taer = taer + (aero_mp(na+1)/rnslyr) & + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+1)/rnslyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+1)/rnslyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) + endif ! na else - taer = taer + & - (aero_mp(na+1)*rnslyr)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+1)*rnslyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+1)*rnslyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) + taer = taer + (aero_mp(na+1)*rnslyr) & !echmod NOTE * rnslyr not / + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+1)*rnslyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+1)*rnslyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) endif ! modal_aero -!mgf-- enddo ! na - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & + w0 (k) = (w0(k)*tau(k) + waer) / & (tau(k) + taer) tau(k) = tau(k) + taer enddo ! k endif ! tr_aero - ! set optical properties of sea ice + ! set optical properties of sea ice - ! bare or snow-covered sea ice layers - !if( srftyp <= 1 ) then - ! ssl - k = kii - tau(k) = (ki_ssl_5bd(ns)+kabs_chl_5bd(ns,k))*dzk(k) - w0(k) = ki_ssl_5bd(ns)/(ki_ssl_5bd(ns) + kabs_chl_5bd(ns,k))*wi_ssl_5bd(ns) - g(k) = gi_ssl_5bd(ns) - ! dl - k = kii + 1 - ! scale dz for dl relative to 4 even-layer-thickness 1.5m case - fs = p25/rnilyr - tau(k) = (ki_dl_5bd(ns) + kabs_chl_5bd(ns,k)) *dzk(k)*fs - w0(k) = ki_dl_5bd(ns)/(ki_dl_5bd(ns) + kabs_chl_5bd(ns,k)) *wi_dl_5bd(ns) - g(k) = gi_dl_5bd(ns) - ! int above lowest layer - if (kii+2 <= klev-1) then - do k = kii+2, klev-1 - tau(k) = (ki_int_5bd(ns) + kabs_chl_5bd(ns,k))*dzk(k) - w0(k) = ki_int_5bd(ns)/(ki_int_5bd(ns) + kabs_chl_5bd(ns,k)) *wi_int_5bd(ns) - g(k) = gi_int_5bd(ns) - enddo - endif - ! lowest layer - k = klev - ! add algae to lowest sea ice layer, visible only: - kabs = ki_int_5bd(ns)*(c1-wi_int_5bd(ns)) - 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_5bd(ns,k) - endif - sig = ki_int_5bd(ns)*wi_int_5bd(ns) - tau(k) = (kabs+sig)*dzk(k) - w0(k) = (sig/(sig+kabs)) - g(k) = gi_int_5bd(ns) - ! aerosol in sea ice - if (tr_zaero .and. dEdd_algae) then - do k = kii, klev - g(k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / & - (w0(k)*tau(k) + wzaer_5bd(ns,k)) - w0(k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / & - (tau(k) + tzaer_5bd(ns,k)) - tau(k) = tau(k) + tzaer_5bd(ns,k) - enddo - elseif (tr_aero) then - k = kii ! sea ice SSL + ! bare or snow-covered sea ice layers + !if (srftyp <= 1) then + ! ssl + k = kii + tau(k) = (ki_ssl_5bd(ns) + kabs_chl_5bd(ns,k)) * dzk(k) + w0 (k) = ki_ssl_5bd(ns)/(ki_ssl_5bd(ns) + kabs_chl_5bd(ns,k)) * wi_ssl_5bd(ns) + g (k) = gi_ssl_5bd(ns) + ! dl + k = kii + 1 + ! scale dz for dl relative to 4 even-layer-thickness 1.5m case + fs = p25/rnilyr + tau(k) = (ki_dl_5bd(ns) + kabs_chl_5bd(ns,k)) * dzk(k) * fs + w0 (k) = ki_dl_5bd(ns)/(ki_dl_5bd(ns) + kabs_chl_5bd(ns,k)) * wi_dl_5bd(ns) + g (k) = gi_dl_5bd(ns) + ! int above lowest layer + if (kii+2 <= klev-1) then + do k = kii+2, klev-1 + tau(k) = (ki_int_5bd(ns) + kabs_chl_5bd(ns,k)) * dzk(k) + w0 (k) = ki_int_5bd(ns)/(ki_int_5bd(ns) + kabs_chl_5bd(ns,k)) * wi_int_5bd(ns) + g (k) = gi_int_5bd(ns) + enddo + endif + ! lowest layer + k = klev + ! add algae to lowest sea ice layer, visible only: + kabs = ki_int_5bd(ns)*(c1-wi_int_5bd(ns)) + 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_5bd(ns,k) + endif + sig = ki_int_5bd(ns)*wi_int_5bd(ns) + tau(k) = (kabs+sig) * dzk(k) + w0 (k) = sig/(sig+kabs) + g (k) = gi_int_5bd(ns) + ! aerosol in sea ice + if (tr_zaero .and. dEdd_algae) then + do k = kii, klev + g (k) = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / & + (w0(k)*tau(k) + wzaer_5bd(ns,k)) + w0 (k) = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / & + (tau(k) + tzaer_5bd(ns,k)) + tau(k) = tau(k) + tzaer_5bd(ns,k) + enddo + elseif (tr_aero) then + k = kii ! sea ice SSL + taer = c0 + waer = c0 + gaer = c0 + do na = 1, 4*n_aero, 4 + if (modal_aero) then + if (na==1) then ! interstitial BC + taer = taer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcexs(k)) + waer = waer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) + gaer = gaer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) & + * gaer_bc_5bd(ns,k_bcexs(k)) + elseif (na==5) then ! within-ice BC + taer = taer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * bcenh_5bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) + gaer = gaer + aero_mp(na+2) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) & + * gaer_bc_5bd(ns,k_bcins(k)) + else ! other species (dust) + taer = taer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) + endif + else ! bulk + taer = taer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + aero_mp(na+2) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) + endif ! modal_aero + enddo ! na + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & + (w0(k)*tau(k) + waer) + w0 (k) = (w0(k)*tau(k) + waer) / & + (tau(k) + taer) + tau(k) = tau(k) + taer + do k = kii+1, klev taer = c0 waer = c0 gaer = c0 - do na=1,4*n_aero,4 - !mgf++ + do na = 1, 4*n_aero, 4 if (modal_aero) then - if (na==1) then - ! interstitial BC - taer = taer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcexs(k)) - waer = waer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k)) - gaer = gaer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))*gaer_bc_5bd(ns,k_bcexs(k)) - elseif (na==5) then - ! within-ice BC - taer = taer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcins(k))* & - bcenh_5bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k)) - gaer = gaer + & - aero_mp(na+2)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))*gaer_bc_5bd(ns,k_bcins(k)) - else - ! other species (dust) - taer = taer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) + if (na==1) then ! interstitial BC + taer = taer + (aero_mp(na+3)/rnilyr) & !echmod: BUG should be *rnilyr + * kaer_bc_5bd(ns,k_bcexs(k)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_5bd(ns,k_bcexs(k)) & + * waer_bc_5bd(ns,k_bcexs(k)) & + * gaer_bc_5bd(ns,k_bcexs(k)) + elseif (na==5) then ! within-ice BC + taer = taer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * bcenh_5bd(ns,k_bcins(k),k_bcini(k)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_bc_5bd(ns,k_bcins(k)) & + * waer_bc_5bd(ns,k_bcins(k)) & + * gaer_bc_5bd(ns,k_bcins(k)) + else ! other species (dust) + taer = taer + (aero_mp(na+3)/rnilyr) & + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+3)/rnilyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+3)/rnilyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) endif - else !bulk - taer = taer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - aero_mp(na+2)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif ! modal_aero - !mgf-- - enddo ! na - - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & + else !bulk + taer = taer + (aero_mp(na+3)*rnilyr) & !echmod NOTE * rnilyr, not / + * kaer_5bd(ns,(1+(na-1)/4)) + waer = waer + (aero_mp(na+3)*rnilyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) + gaer = gaer + (aero_mp(na+3)*rnilyr) & + * kaer_5bd(ns,(1+(na-1)/4)) & + * waer_5bd(ns,(1+(na-1)/4)) & + * gaer_5bd(ns,(1+(na-1)/4)) + endif ! modal_aero + enddo ! na + g (k) = (g(k)*w0(k)*tau(k) + gaer) / & (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & + w0 (k) = (w0(k)*tau(k) + waer) / & (tau(k) + taer) tau(k) = tau(k) + taer - do k = kii+1, klev - taer = c0 - waer = c0 - gaer = c0 - do na=1,4*n_aero,4 - !mgf++ - if (modal_aero) then - if (na==1) then - ! interstitial BC - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcexs(k)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcexs(k))* & - waer_bc_5bd(ns,k_bcexs(k))*gaer_bc_5bd(ns,k_bcexs(k)) - elseif (na==5) then - ! within-ice BC - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcins(k))* & - bcenh_5bd(ns,k_bcins(k),k_bcini(k)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_bc_5bd(ns,k_bcins(k))* & - waer_bc_5bd(ns,k_bcins(k))*gaer_bc_5bd(ns,k_bcins(k)) - - else - ! other species (dust) - taer = taer + & - (aero_mp(na+3)/rnilyr)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+3)/rnilyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+3)/rnilyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif - else !bulk - - taer = taer + & - (aero_mp(na+3)*rnilyr)*kaer_5bd(ns,(1+(na-1)/4)) - waer = waer + & - (aero_mp(na+3)*rnilyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4)) - gaer = gaer + & - (aero_mp(na+3)*rnilyr)*kaer_5bd(ns,(1+(na-1)/4))* & - waer_5bd(ns,(1+(na-1)/4))*gaer_5bd(ns,(1+(na-1)/4)) - endif ! modal_aero - !mgf-- - enddo ! na - g(k) = (g(k)*w0(k)*tau(k) + gaer) / & - (w0(k)*tau(k) + waer) - w0(k) = (w0(k)*tau(k) + waer) / & - (tau(k) + taer) - tau(k) = tau(k) + taer - enddo ! k - endif ! tr_aero + enddo ! k + endif ! tr_aero + ! --------------------------------------------------------------------------- - ! set reflectivities for ocean underlying sea ice - ! if ns == 1 (visible), albedo is 0.1, else, albedo is zero - 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 - ! reflectivities and transmissivities for each layer; then, - ! combine the layers going downwards accounting for multiple - ! 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, & - tau, w0, g, albodr, albodf, & - trndir, trntdr, trndif, rupdir, rupdif, & - rdndif) - ! 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 - ! combined layer properties at each interface: - ! - ! layers interface - ! - ! --------------------- k - ! k - ! --------------------- - - do k = 0, klevp - ! interface scattering - refk = c1/(c1 - rdndif(k)*rupdif(k)) - ! dir tran ref from below times interface scattering, plus diff - ! tran and ref from below times interface scattering - ! fdirup(k) = (trndir(k)*rupdir(k) + & - ! (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 - ! fdirdn(k) = trndir(k) + (trntdr(k) & - ! - trndir(k) + trndir(k) & - ! *rupdir(k)*rdndif(k))*refk - ! diffuse tran ref from below times interface scattering - ! fdifup(k) = trndif(k)*rupdif(k)*refk - ! diffuse tran times interface scattering - ! fdifdn(k) = trndif(k)*refk - - ! dfdir = fdirdn - fdirup - dfdir(k) = trndir(k) & - + (trntdr(k)-trndir(k)) * (c1 - rupdif(k)) * refk & - - trndir(k)*rupdir(k) * (c1 - rdndif(k)) * refk - if (dfdir(k) < puny) dfdir(k) = c0 !echmod necessary? - ! dfdif = fdifdn - fdifup - dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk - if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary? - enddo ! k - - ! note that because the snow IOPs for diffuse and direct incidents - ! are different, the snow albedo needs to be calculated twice for - ! direct incident and diffuse incident respectively - if (nsky == 1) then ! direc beam (keep the direct beam results) - do k = 0, klevp - dfdir_snicar(k) = dfdir(k) - rupdir_snicar(k) = rupdir(k) - enddo - elseif (nsky == 2) then ! diffuse (keep the diffuse incident results) - do k = 0, klevp + ! set reflectivities for ocean underlying sea ice + ! if ns == 1 (visible), albedo is 0.1, else, albedo is zero + 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 + ! reflectivities and transmissivities for each layer; then, + ! combine the layers going downwards accounting for multiple + ! 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, & + tau, w0, g, albodr, albodf, & + trndir, trntdr, trndif, rupdir, rupdif, & + 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 + ! combined layer properties at each interface: + ! + ! layers interface + ! + ! --------------------- k + ! k + ! --------------------- + + do k = 0, klevp + ! interface scattering + refk = c1/(c1 - rdndif(k)*rupdif(k)) + ! dir tran ref from below times interface scattering, plus diff + ! tran and ref from below times interface scattering + ! fdirup(k) = (trndir(k)*rupdir(k) + & + ! (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 + ! fdirdn(k) = trndir(k) + (trntdr(k) & + ! - trndir(k) + trndir(k) & + ! *rupdir(k)*rdndif(k))*refk + ! diffuse tran ref from below times interface scattering + ! fdifup(k) = trndif(k)*rupdif(k)*refk + ! diffuse tran times interface scattering + ! fdifdn(k) = trndif(k)*refk + + ! dfdir = fdirdn - fdirup + dfdir(k) = trndir(k) & + + (trntdr(k)-trndir(k)) * (c1 - rupdif(k)) * refk & + - trndir(k)*rupdir(k) * (c1 - rdndif(k)) * refk + if (dfdir(k) < puny) dfdir(k) = c0 !echmod necessary? + ! dfdif = fdifdn - fdifup + dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk + if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary? + enddo ! k + + ! note that because the snow IOPs for diffuse and direct incidents + ! are different, the snow albedo needs to be calculated twice for + ! direct incident and diffuse incident respectively + if (nsky == 1) then ! direct beam (keep the direct beam results) + do k = 0, klevp + dfdir_snicar(k) = dfdir(k) + rupdir_snicar(k) = rupdir(k) + enddo + elseif (nsky == 2) then ! diffuse (keep the diffuse incident results) + do k = 0, klevp dfdif_snicar(k) = dfdif(k) rupdif_snicar(k) = rupdif(k) - enddo - endif - enddo ! end direct/diffuse incident nsky - - ! 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_snicar(0) - avdf = rupdif_snicar(0) - tmp_0 = dfdir_snicar(0 )*swdr + dfdif_snicar(0 )*swdf - tmp_ks = dfdir_snicar(ksrf )*swdr + dfdif_snicar(ksrf )*swdf - tmp_kl = dfdir_snicar(klevp)*swdr + dfdif_snicar(klevp)*swdf - - ! for layer biology: save visible only - do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering - fthrul(k-nslyr-1) = dfdir_snicar(k)*swdr + dfdif_snicar(k)*swdf - enddo + enddo + endif + enddo ! end direct/diffuse nsky loop ------------------------------------ + + ! 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_snicar(0) + avdf = rupdif_snicar(0) + tmp_0 = dfdir_snicar(0 )*swdr + dfdif_snicar(0 )*swdf + tmp_ks = dfdir_snicar(ksrf )*swdr + dfdif_snicar(ksrf )*swdf + tmp_kl = dfdir_snicar(klevp)*swdr + dfdif_snicar(klevp)*swdf + + ! for layer biology: save visible only + do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering + fthrul(k-nslyr-1) = dfdir_snicar(k)*swdr + dfdif_snicar(k)*swdf + enddo - fsfc = fsfc + tmp_0 - tmp_ks - fint = fint + tmp_ks - tmp_kl - fthru = fthru + tmp_kl + fsfc = fsfc + tmp_0 - tmp_ks + fint = fint + tmp_ks - tmp_kl + fthru = fthru + tmp_kl + fthruvdr = fthruvdr + dfdir_snicar(klevp)*swdr + fthruvdf = fthruvdf + dfdif_snicar(klevp)*swdf - ! if snow covered ice, set snow internal absorption; else, Sabs=0 - if( srftyp == 1 ) then + ! if snow covered ice, set snow internal absorption; else, Sabs=0 + if (srftyp == 1) then + ki = 0 + do k = 1, nslyr + ! skip snow SSL, since SSL absorption included in the surface + ! absorption fsfc above + km = k + kp = km + 1 + ki = ki + 1 + Sabs(ki) = Sabs(ki) & + + dfdir_snicar(km)*swdr + dfdif_snicar(km)*swdf & + - (dfdir_snicar(kp)*swdr + dfdif_snicar(kp)*swdf) + enddo ! k + endif + + ! complex indexing to insure proper absorptions for sea ice ki = 0 - do k=1,nslyr - ! skip snow SSL, since SSL absorption included in the surface - ! absorption fsfc above - km = k - kp = km + 1 - ki = ki + 1 - Sabs(ki) = Sabs(ki) & + do k = nslyr+2, nslyr+1+nilyr + ! for bare ice, DL absorption for sea ice layer 1 + km = k + kp = km + 1 + ! modify for top sea ice layer for snow over sea ice + if (srftyp == 1) then + ! must add SSL and DL absorption for sea ice layer 1 + if (k == nslyr+2) then + km = k - 1 + kp = km + 2 + endif + endif + ki = ki + 1 + Iabs(ki) = Iabs(ki) & + dfdir_snicar(km)*swdr + dfdif_snicar(km)*swdf & - (dfdir_snicar(kp)*swdr + dfdif_snicar(kp)*swdf) enddo ! k - endif - ! complex indexing to insure proper absorptions for sea ice - ki = 0 - do k=nslyr+2,nslyr+1+nilyr - ! for bare ice, DL absorption for sea ice layer 1 - km = k - kp = km + 1 - ! modify for top sea ice layer for snow over sea ice - if( srftyp == 1 ) then - ! must add SSL and DL absorption for sea ice layer 1 - if( k == nslyr+2 ) then - km = k - 1 - kp = km + 2 - endif - endif - ki = ki + 1 - Iabs(ki) = Iabs(ki) & - + dfdir_snicar(km)*swdr + dfdif_snicar(km)*swdf & - - (dfdir_snicar(kp)*swdr + dfdif_snicar(kp)*swdf) - enddo ! k + else ! ns > 1, near IR - else !if(ns > 1) then ! near IR - - swdr = swidr - swdf = swidf - - ! let fr2(3,4,5) = alb_2(3,4,5)*swd*wght2(3,4,5) - ! the ns=2(3,4,5) reflected fluxes respectively, - ! where alb_2(3,4,5) are the band - ! albedos, swd = nir incident shortwave flux, and wght2(3,4,5) are - ! the 2(3,4,5) band weights. thus, the total reflected flux is: - ! fr = fr2 + fr3 + fr4 + fr5 - ! = alb_2*swd*wght2 + alb_3*swd*wght3 + alb_4*swd*wght4 + alb_5*swd*wght5 - ! hence, the 2,3,4,5 nir band albedo is - ! alb = fr/swd = alb_2*wght2 + alb_3*wght3 + alb_4*wght4 + alb_5*wght5 - - aidr = aidr + rupdir_snicar(0)*wghtns_5bd_drc(ns) - aidf = aidf + rupdif_snicar(0)*wghtns_5bd_dfs(ns) - - tmp_0 = dfdir_snicar(0 )*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(0 )*swdf*wghtns_5bd_dfs(ns) - tmp_ks = dfdir_snicar(ksrf )*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(ksrf )*swdf*wghtns_5bd_dfs(ns) - tmp_kl = dfdir_snicar(klevp)*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(klevp)*swdf*wghtns_5bd_dfs(ns) - - fsfc = fsfc + tmp_0 - tmp_ks - fint = fint + tmp_ks - tmp_kl - fthru = fthru + tmp_kl - - ! if snow covered ice, set snow internal absorption; else, Sabs=0 - if( srftyp == 1 ) then - ki = 0 - do k=1,nslyr - ! skip snow SSL, since SSL absorption included in the surface - ! absorption fsfc above - km = k - kp = km + 1 - ki = ki + 1 - Sabs(ki) = Sabs(ki) & - + dfdir_snicar(km)*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(km)*swdf*wghtns_5bd_dfs(ns) & - -(dfdir_snicar(kp)*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(kp)*swdf*wghtns_5bd_dfs(ns)) + swdr = swidr + swdf = swidf + ! let fr2(3,4,5) = alb_2(3,4,5)*swd*wght2(3,4,5) + ! the ns=2(3,4,5) reflected fluxes respectively, + ! where alb_2(3,4,5) are the band + ! albedos, swd = nir incident shortwave flux, and wght2(3,4,5) are + ! the 2(3,4,5) band weights. thus, the total reflected flux is: + ! fr = fr2 + fr3 + fr4 + fr5 + ! = alb_2*swd*wght2 + alb_3*swd*wght3 + alb_4*swd*wght4 + alb_5*swd*wght5 + ! hence, the 2,3,4,5 nir band albedo is + ! alb = fr/swd = alb_2*wght2 + alb_3*wght3 + alb_4*wght4 + alb_5*wght5 + + aidr = aidr + rupdir_snicar(0)*wghtns_5bd_drc(ns) + aidf = aidf + rupdif_snicar(0)*wghtns_5bd_dfs(ns) + + tmp_0 = dfdir_snicar(0 )*swdr*wghtns_5bd_drc(ns) & + + dfdif_snicar(0 )*swdf*wghtns_5bd_dfs(ns) + tmp_ks = dfdir_snicar(ksrf )*swdr*wghtns_5bd_drc(ns) & + + dfdif_snicar(ksrf )*swdf*wghtns_5bd_dfs(ns) + tmp_kl = dfdir_snicar(klevp)*swdr*wghtns_5bd_drc(ns) & + + dfdif_snicar(klevp)*swdf*wghtns_5bd_dfs(ns) + + fsfc = fsfc + tmp_0 - tmp_ks + fint = fint + tmp_ks - tmp_kl + fthru = fthru + tmp_kl + fthruidr = fthruidr + dfdir_snicar(klevp)*swdr*wghtns_5bd_drc(ns) + fthruidf = fthruidf + dfdif_snicar(klevp)*swdf*wghtns_5bd_dfs(ns) + + ! if snow covered ice, set snow internal absorption; else, Sabs=0 + if (srftyp == 1) then + ki = 0 + do k = 1, nslyr + ! skip snow SSL, since SSL absorption included in the surface + ! absorption fsfc above + km = k + kp = km + 1 + ki = ki + 1 + Sabs(ki) = Sabs(ki) & + + dfdir_snicar(km)*swdr*wghtns_5bd_drc(ns) & + + dfdif_snicar(km)*swdf*wghtns_5bd_dfs(ns) & + - dfdir_snicar(kp)*swdr*wghtns_5bd_drc(ns) & + - dfdif_snicar(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 - ! for bare ice, DL absorption for sea ice layer 1 - km = k - kp = km + 1 - ! modify for top sea ice layer for snow over sea ice - if( srftyp == 1 ) then - ! must add SSL and DL absorption for sea ice layer 1 - if( k == nslyr+2 ) then - km = k - 1 - kp = km + 2 + ! complex indexing to insure proper absorptions for sea ice + ki = 0 + do k = nslyr+2, nslyr+1+nilyr + ! for bare ice, DL absorption for sea ice layer 1 + km = k + kp = km + 1 + ! modify for top sea ice layer for snow over sea ice + if (srftyp == 1) then + ! must add SSL and DL absorption for sea ice layer 1 + if (k == nslyr+2) then + km = k - 1 + kp = km + 2 + endif endif - endif - ki = ki + 1 - Iabs(ki) = Iabs(ki) & - + dfdir_snicar(km)*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(km)*swdf*wghtns_5bd_dfs(ns) & - -(dfdir_snicar(kp)*swdr*wghtns_5bd_drc(ns) & - + dfdif_snicar(kp)*swdf*wghtns_5bd_dfs(ns)) - enddo ! k - endif ! ns = 1, ns > 1 - enddo ! end spectral loop ns - - - ! accumulate fluxes over bare sea ice + ki = ki + 1 + Iabs(ki) = Iabs(ki) & + + dfdir_snicar(km)*swdr*wghtns_5bd_drc(ns) & + + dfdif_snicar(km)*swdf*wghtns_5bd_dfs(ns) & + - dfdir_snicar(kp)*swdr*wghtns_5bd_drc(ns) & + - dfdif_snicar(kp)*swdf*wghtns_5bd_dfs(ns) + enddo ! k + endif ! ns + enddo ! ns: end spectral loop ! solar zenith angle parameterization - ! calculate the scaling factor for NIR direct albedo if SZA>75 degree + ! calculate the scaling factor for NIR direct albedo if SZA>75 degrees sza_factor = c1 - if( srftyp == 1 ) then - mu0 = max(coszen,p01) + if (srftyp == 1) then + mu0 = max(coszen, p01) if (mu0 < mu_75) then sza_c1 = sza_a0 + sza_a1 * mu0 + sza_a2 * mu0**2 sza_c0 = sza_b0 + sza_b1 * mu0 + sza_b2 * mu0**2 - sza_factor = sza_c1 * (log10(rsnw(1)) - 6.0) + sza_c0 - endif + sza_factor = sza_c1 * (log10(rsnw(1)) - 6.0_dbl_kind) + sza_c0 + endif endif - alvdr = avdr - alvdf = avdf - alidr = aidr * sza_factor !sza factor is always larger than or equal to 1 - alidf = aidf + alvdr = avdr + alvdf = avdf + alidr = aidr * sza_factor !sza factor is always larger than or equal to 1 + alidf = aidf + + ! accumulate fluxes over bare sea ice ! note that we assume the reduced NIR energy absorption by snow ! due to corrected snow albedo is absorbed by the snow single ! scattering layer only - this is generally true if snow SSL >= 2 cm ! by the default model set up: ! if snow_depth >= 8 cm, SSL = 4 cm, satisfy - ! esle if snow_depth >= 4 cm, SSL = snow_depth/2 >= 2 cm, satisfy - ! esle snow_depth < 4 cm, SSL = snow_depth/2, may overcool SSL layer + ! else if snow_depth >= 4 cm, SSL = snow_depth/2 >= 2 cm, satisfy + ! else snow_depth < 4 cm, SSL = snow_depth/2, may overcool SSL layer fswsfc = fswsfc + (fsfc- (sza_factor-c1)*aidr*swidr)*fi fswint = fswint + fint *fi fswthru = fswthru + fthru*fi - + fswthru_vdr = fswthru_vdr + fthruvdr*fi + fswthru_vdf = fswthru_vdf + fthruvdf*fi + fswthru_idr = fswthru_idr + fthruidr*fi + fswthru_idf = fswthru_idf + fthruidf*fi do k = 1, nslyr Sswabs(k) = Sswabs(k) + Sabs(k)*fi - enddo ! k + enddo do k = 1, nilyr Iswabs(k) = Iswabs(k) + Iabs(k)*fi - ! 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 + fswpenl(k) = fswpenl(k) + fthrul(k)*fi + enddo + fswpenl(nilyr+1) = fswpenl(nilyr+1) + fthrul(nilyr+1)*fi +#ifdef UNDEPRECATE_0LAYER !---------------------------------------------------------------- ! if ice has zero heat capacity, no SW can be absorbed ! in the ice/snow interior, so add to surface absorption. @@ -5503,7 +5476,7 @@ subroutine compute_dEdd_5bd (klev, klevp, & ! Sswabs(1) = c0 ! endif ! heat_capacity - +#endif end subroutine compute_dEdd_5bd !=======================================================================