LCOV - code coverage report
Current view: top level - icepack/columnphysics - icepack_shortwave.F90 (source / functions) Hit Total Coverage
Test: 200617-180449:aec9683041:7:first,base,travis,decomp,reprosum,io,quick Lines: 1131 1554 72.78 %
Date: 2020-06-17 18:05:09 Functions: 20 21 95.24 %

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! The albedo and absorbed/transmitted flux parameterizations for
       4             : ! snow over ice, bare ice and ponded ice.
       5             : !
       6             : ! Presently, two methods are included:
       7             : !   (1) CCSM3 
       8             : !   (2) Delta-Eddington 
       9             : ! as two distinct routines.
      10             : ! Either can be called from the ice driver.
      11             : !
      12             : ! The Delta-Eddington method is described here:
      13             : !
      14             : ! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple 
      15             : !    Scattering Parameterization for Solar Radiation in the Sea Ice 
      16             : !    Component of the Community Climate System Model, NCAR Technical 
      17             : !    Note  NCAR/TN-472+STR  February 2007
      18             : !
      19             : ! name: originally ice_albedo
      20             : !
      21             : ! authors:  Bruce P. Briegleb, NCAR
      22             : !           Elizabeth C. Hunke and William H. Lipscomb, LANL
      23             : ! 2005, WHL: Moved absorbed_solar from icepack_therm_vertical to this 
      24             : !            module and changed name from ice_albedo
      25             : ! 2006, WHL: Added Delta Eddington routines from Bruce Briegleb
      26             : ! 2006, ECH: Changed data statements in Delta Eddington routines (no 
      27             : !            longer hardwired)
      28             : !            Converted to free source form (F90)
      29             : ! 2007, BPB: Completely updated Delta-Eddington code, so that:
      30             : !            (1) multiple snow layers enabled (i.e. nslyr > 1)
      31             : !            (2) included SSL for snow surface absorption
      32             : !            (3) added Sswabs for internal snow layer absorption
      33             : !            (4) variable sea ice layers allowed (i.e. not hardwired)
      34             : !            (5) updated all inherent optical properties
      35             : !            (6) included algae absorption for sea ice lowest layer
      36             : !            (7) very complete internal documentation included
      37             : ! 2007, ECH: Improved efficiency
      38             : ! 2008, BPB: Added aerosols to Delta Eddington code
      39             : ! 2013, ECH: merged with NCAR version, cleaned up
      40             : 
      41             :       module icepack_shortwave
      42             : 
      43             :       use icepack_kinds
      44             :       use icepack_parameters, only: c0, c1, c1p5, c2, c3, c4, c10
      45             :       use icepack_parameters, only: p01, p1, p15, p25, p5, p75, puny
      46             :       use icepack_parameters, only: albocn, Timelt, snowpatch, awtvdr, awtidr, awtvdf, awtidf
      47             :       use icepack_parameters, only: kappav, hs_min, rhofresh, rhos, nspint
      48             :       use icepack_parameters, only: hi_ssl, hs_ssl, min_bgc, sk_l
      49             :       use icepack_parameters, only: z_tracers, skl_bgc, calc_tsfc, shortwave, kalg, heat_capacity
      50             :       use icepack_parameters, only: r_ice, r_pnd, r_snw, dt_mlt, rsnw_mlt, hs0, hs1, hp1
      51             :       use icepack_parameters, only: pndaspect, albedo_type, albicev, albicei, albsnowv, albsnowi, ahmax
      52             :       use icepack_tracers,    only: ntrcr, nbtrcr_sw
      53             :       use icepack_tracers,    only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
      54             :       use icepack_tracers,    only: tr_bgc_N, tr_aero
      55             :       use icepack_tracers,    only: nt_bgc_N, nt_zaero, tr_bgc_N
      56             :       use icepack_tracers,    only: tr_zaero, nlt_chl_sw, nlt_zaero_sw
      57             :       use icepack_tracers,    only: n_algae, n_aero, n_zaero
      58             :       use icepack_warnings,   only: warnstr, icepack_warnings_add
      59             :       use icepack_warnings,   only: icepack_warnings_setabort, icepack_warnings_aborted
      60             : 
      61             :       use icepack_zbgc_shared,only: R_chl2N, F_abs_chl
      62             :       use icepack_zbgc_shared,only: remap_zbgc
      63             :       use icepack_orbital, only: compute_coszen
      64             : 
      65             : 
      66             :       implicit none
      67             : 
      68             :       private
      69             :       public :: run_dEdd, &
      70             :                 shortwave_ccsm3, &
      71             :                 compute_shortwave_trcr, &
      72             :                 icepack_prep_radiation, &
      73             :                 icepack_step_radiation
      74             : 
      75             :       real (kind=dbl_kind), parameter :: &
      76             :          hpmin  = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
      77             :          hp0    = 0.200_dbl_kind    ! pond depth below which transition to bare ice
      78             : 
      79             :       real (kind=dbl_kind), parameter :: & 
      80             :          exp_argmax = c10    ! maximum argument of exponential
      81             : 
      82             : !=======================================================================
      83             : 
      84             :       contains
      85             : 
      86             : !=======================================================================
      87             : !
      88             : ! Driver for basic solar radiation from CCSM3.  Albedos and absorbed solar.
      89             : 
      90    92741504 :       subroutine shortwave_ccsm3 (aicen,    vicen,    &
      91    46370752 :                                   vsnon,    Tsfcn,    &
      92             :                                   swvdr,    swvdf,    &
      93             :                                   swidr,    swidf,    &
      94             :                                   heat_capacity,      &
      95     5796344 :                                   albedo_type,        &
      96             :                                   albicev,  albicei,  &
      97             :                                   albsnowv, albsnowi, &
      98             :                                   ahmax,              &
      99    46370752 :                                   alvdrn,   alidrn,   &
     100    46370752 :                                   alvdfn,   alidfn,   &
     101    46370752 :                                   fswsfc,   fswint,   &
     102    46370752 :                                   fswthru,  fswpenl,  &
     103    92741504 :                                   Iswabs,   SSwabs,   &
     104    46370752 :                                   albin,    albsn,    &
     105             :                                   coszen,   ncat,     &
     106             :                                   nilyr)
     107             : 
     108             :       integer (kind=int_kind), intent(in) :: &
     109             :          nilyr    , & ! number of ice layers
     110             :          ncat         ! number of ice thickness categories
     111             : 
     112             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     113             :          aicen    , & ! concentration of ice per category
     114             :          vicen    , & ! volume of ice per category
     115             :          vsnon    , & ! volume of ice per category
     116             :          Tsfcn        ! surface temperature
     117             : 
     118             :       real (kind=dbl_kind), intent(in) :: &
     119             :          swvdr    , & ! sw down, visible, direct  (W/m^2)
     120             :          swvdf    , & ! sw down, visible, diffuse (W/m^2)
     121             :          swidr    , & ! sw down, near IR, direct  (W/m^2)
     122             :          swidf        ! sw down, near IR, diffuse (W/m^2)
     123             : 
     124             :       ! baseline albedos for ccsm3 shortwave, set in namelist
     125             :       real (kind=dbl_kind), intent(in) :: &
     126             :          albicev , & ! visible ice albedo for h > ahmax
     127             :          albicei , & ! near-ir ice albedo for h > ahmax
     128             :          albsnowv, & ! cold snow albedo, visible
     129             :          albsnowi, & ! cold snow albedo, near IR
     130             :          ahmax       ! thickness above which ice albedo is constant (m)
     131             : 
     132             :       logical(kind=log_kind), intent(in) :: &
     133             :          heat_capacity! if true, ice has nonzero heat capacity
     134             : 
     135             :       character (len=char_len), intent(in) :: &
     136             :          albedo_type  ! albedo parameterization, 'ccsm3' or 'constant'
     137             : 
     138             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     139             :          alvdrn   , & ! visible, direct, avg   (fraction)
     140             :          alidrn   , & ! near-ir, direct, avg   (fraction)
     141             :          alvdfn   , & ! visible, diffuse, avg  (fraction)
     142             :          alidfn   , & ! near-ir, diffuse, avg  (fraction)
     143             :          fswsfc   , & ! SW absorbed at ice/snow surface (W m-2)
     144             :          fswint   , & ! SW absorbed in ice interior, below surface (W m-2)
     145             :          fswthru  , & ! SW through ice to ocean (W m-2)
     146             :          albin    , & ! bare ice albedo
     147             :          albsn        ! snow albedo
     148             : 
     149             :       real (kind=dbl_kind), intent(inout) :: &
     150             :          coszen       ! cosine(zenith angle)
     151             : 
     152             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
     153             :          fswpenl  , & ! SW entering ice layers (W m-2)
     154             :          Iswabs   , & ! SW absorbed in particular layer (W m-2)
     155             :          Sswabs       ! SW absorbed in particular layer (W m-2)
     156             : 
     157             :       ! local variables
     158             : 
     159             :       integer (kind=int_kind) :: &
     160             :          n                  ! thickness category index
     161             : 
     162             :       ! ice and snow albedo for each category
     163             : 
     164             :       real (kind=dbl_kind) :: &
     165     5796344 :          alvdrni, & ! visible, direct, ice    (fraction)
     166     5796344 :          alidrni, & ! near-ir, direct, ice    (fraction)
     167     5796344 :          alvdfni, & ! visible, diffuse, ice   (fraction)
     168     5796344 :          alidfni, & ! near-ir, diffuse, ice   (fraction)
     169     5796344 :          alvdrns, & ! visible, direct, snow   (fraction)
     170     5796344 :          alidrns, & ! near-ir, direct, snow   (fraction)
     171     5796344 :          alvdfns, & ! visible, diffuse, snow  (fraction)
     172     5796344 :          alidfns    ! near-ir, diffuse, snow  (fraction)
     173             : 
     174             :       character(len=*),parameter :: subname='(shortwave_ccsm3)'
     175             : 
     176             :       !-----------------------------------------------------------------
     177             :       ! Solar radiation: albedo and absorbed shortwave
     178             :       !-----------------------------------------------------------------
     179             : 
     180             :       ! For basic shortwave, set coszen to a constant between 0 and 1.
     181    46370752 :       coszen = p5 ! sun above the horizon
     182             : 
     183   185483008 :       do n = 1, ncat
     184             : 
     185   278224512 :          Sswabs(:,n) = c0
     186             : 
     187   139112256 :       alvdrni = albocn
     188   139112256 :       alidrni = albocn
     189   139112256 :       alvdfni = albocn
     190   139112256 :       alidfni = albocn
     191             :       
     192   139112256 :       alvdrns = albocn
     193   139112256 :       alidrns = albocn
     194   139112256 :       alvdfns = albocn
     195   139112256 :       alidfns = albocn
     196             :       
     197   139112256 :       alvdrn(n) = albocn
     198   139112256 :       alidrn(n) = albocn
     199   139112256 :       alvdfn(n) = albocn
     200   139112256 :       alidfn(n) = albocn
     201             :       
     202   139112256 :       albin(n) = c0
     203   139112256 :       albsn(n) = c0    
     204             : 
     205   139112256 :       fswsfc(n)  = c0
     206   139112256 :       fswint(n)  = c0
     207   139112256 :       fswthru(n) = c0
     208   556449024 :       fswpenl(:,n)  = c0
     209   417336768 :       Iswabs (:,n) = c0
     210             : 
     211   185483008 :       if (aicen(n) > puny) then
     212             : 
     213             :       !-----------------------------------------------------------------
     214             :       ! Compute albedos for ice and snow.
     215             :       !-----------------------------------------------------------------
     216             : 
     217    42678992 :          if (trim(albedo_type) == 'constant') then
     218             : 
     219     4680160 :             call constant_albedos (aicen(n),             &
     220     4680160 :                                    vsnon(n),             &
     221     4680160 :                                    Tsfcn(n),             &
     222             :                                    alvdrni,    alidrni,  &
     223             :                                    alvdfni,    alidfni,  &
     224             :                                    alvdrns,    alidrns,  &
     225             :                                    alvdfns,    alidfns,  &
     226     4680160 :                                    alvdrn(n),            &
     227     4680160 :                                    alidrn(n),            &
     228     4680160 :                                    alvdfn(n),            &
     229     4680160 :                                    alidfn(n),            &
     230     4680160 :                                    albin(n),             &
     231    37441280 :                                    albsn(n))
     232    37441280 :             if (icepack_warnings_aborted(subname)) return
     233             : 
     234     5237712 :          elseif (trim(albedo_type) == 'ccsm3') then
     235             : 
     236      654714 :             call compute_albedos (aicen(n),             &
     237      654714 :                                   vicen(n),             &
     238      654714 :                                   vsnon(n),             &
     239      654714 :                                   Tsfcn(n),             &
     240             :                                   albicev,    albicei,  &
     241             :                                   albsnowv,   albsnowi, &
     242             :                                   ahmax,                &
     243             :                                   alvdrni,    alidrni,  &
     244             :                                   alvdfni,    alidfni,  &
     245             :                                   alvdrns,    alidrns,  &
     246             :                                   alvdfns,    alidfns,  &
     247      654714 :                                   alvdrn(n),            &
     248      654714 :                                   alidrn(n),            &
     249      654714 :                                   alvdfn(n),            &
     250      654714 :                                   alidfn(n),            &
     251      654714 :                                   albin(n),             &
     252     5237712 :                                   albsn(n))
     253     5237712 :             if (icepack_warnings_aborted(subname)) return
     254             : 
     255             :          else
     256             : 
     257           0 :             call icepack_warnings_add(subname//' ERROR: albedo_type '//trim(albedo_type)//' unknown')
     258           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     259           0 :             return
     260             : 
     261             :          endif
     262             : 
     263             :       !-----------------------------------------------------------------
     264             :       ! Compute solar radiation absorbed in ice and penetrating to ocean.
     265             :       !-----------------------------------------------------------------
     266             : 
     267             :          call absorbed_solar  (heat_capacity,        &
     268             :                                nilyr,                &
     269     5334874 :                                aicen(n),             &
     270     5334874 :                                vicen(n),             &
     271     5334874 :                                vsnon(n),             &
     272             :                                swvdr,      swvdf,    &
     273             :                                swidr,      swidf,    &
     274             :                                alvdrni,    alvdfni,  &
     275             :                                alidrni,    alidfni,  &
     276             :                                alvdrns,    alvdfns,  &
     277             :                                alidrns,    alidfns,  &
     278     5334874 :                                fswsfc(n),            &
     279     5334874 :                                fswint(n),            &
     280     5334874 :                                fswthru(n),           &
     281           0 :                                fswpenl(:,n),         &
     282    42678992 :                                Iswabs(:,n))
     283    42678992 :          if (icepack_warnings_aborted(subname)) return
     284             : 
     285             :       endif ! aicen > puny
     286             : 
     287             :       enddo                  ! ncat
     288             : 
     289             :       end subroutine shortwave_ccsm3
     290             : 
     291             : !=======================================================================
     292             : !
     293             : ! Compute albedos for each thickness category
     294             : 
     295     5237712 :       subroutine compute_albedos (aicen,    vicen,    &
     296             :                                   vsnon,    Tsfcn,    &
     297             :                                   albicev,  albicei,  &
     298             :                                   albsnowv, albsnowi, &
     299             :                                   ahmax,              &
     300             :                                   alvdrni,  alidrni,  &
     301             :                                   alvdfni,  alidfni,  &
     302             :                                   alvdrns,  alidrns,  &
     303             :                                   alvdfns,  alidfns,  &
     304             :                                   alvdrn,   alidrn,   &
     305             :                                   alvdfn,   alidfn,   &
     306             :                                   albin,    albsn)
     307             : 
     308             :       real (kind=dbl_kind), intent(in) :: &
     309             :          aicen   , & ! concentration of ice per category
     310             :          vicen   , & ! volume of ice per category
     311             :          vsnon   , & ! volume of ice per category
     312             :          Tsfcn       ! surface temperature
     313             : 
     314             :       ! baseline albedos for ccsm3 shortwave, set in namelist
     315             :       real (kind=dbl_kind), intent(in) :: &
     316             :          albicev , & ! visible ice albedo for h > ahmax
     317             :          albicei , & ! near-ir ice albedo for h > ahmax
     318             :          albsnowv, & ! cold snow albedo, visible
     319             :          albsnowi, & ! cold snow albedo, near IR
     320             :          ahmax       ! thickness above which ice albedo is constant (m)
     321             : 
     322             :       real (kind=dbl_kind), intent(out) :: &
     323             :          alvdrni  , & ! visible, direct, ice   (fraction)
     324             :          alidrni  , & ! near-ir, direct, ice   (fraction)
     325             :          alvdfni  , & ! visible, diffuse, ice  (fraction)
     326             :          alidfni  , & ! near-ir, diffuse, ice  (fraction)
     327             :          alvdrns  , & ! visible, direct, snow  (fraction)
     328             :          alidrns  , & ! near-ir, direct, snow  (fraction)
     329             :          alvdfns  , & ! visible, diffuse, snow (fraction)
     330             :          alidfns  , & ! near-ir, diffuse, snow (fraction)
     331             :          alvdrn   , & ! visible, direct, avg   (fraction)
     332             :          alidrn   , & ! near-ir, direct, avg   (fraction)
     333             :          alvdfn   , & ! visible, diffuse, avg  (fraction)
     334             :          alidfn   , & ! near-ir, diffuse, avg  (fraction)
     335             :          albin    , & ! bare ice 
     336             :          albsn        ! snow 
     337             : 
     338             :       ! local variables
     339             : 
     340             :       real (kind=dbl_kind), parameter :: &
     341             :          dT_melt    = c1          , & ! change in temp to give dalb_mlt 
     342             :                                      ! albedo change
     343             :          dalb_mlt  = -0.075_dbl_kind, & ! albedo change per dT_melt change
     344             :                                      ! in temp for ice
     345             :          dalb_mltv = -p1         , & ! albedo vis change per dT_melt change
     346             :                                      ! in temp for snow
     347             :          dalb_mlti = -p15            ! albedo nir change per dT_melt change
     348             :                                      ! in temp for snow
     349             : 
     350             :       real (kind=dbl_kind) :: &
     351      654714 :          hi  , & ! ice thickness  (m)
     352      654714 :          hs  , & ! snow thickness  (m)
     353      654714 :          albo, & ! effective ocean albedo, function of ice thickness
     354      654714 :          fh  , & ! piecewise linear function of thickness
     355      654714 :          fT  , & ! piecewise linear function of surface temperature
     356      654714 :          dTs , & ! difference of Tsfc and Timelt
     357      654714 :          fhtan,& ! factor used in albedo dependence on ice thickness
     358      654714 :          asnow   ! fractional area of snow cover
     359             : 
     360             :       character(len=*),parameter :: subname='(compute_albedos)'
     361             : 
     362     5237712 :       fhtan = atan(ahmax*c4)
     363             : 
     364             :       !-----------------------------------------------------------------
     365             :       ! Compute albedo for each thickness category.
     366             :       !-----------------------------------------------------------------
     367             : 
     368     5237712 :          hi = vicen / aicen
     369     5237712 :          hs = vsnon / aicen            
     370             : 
     371             :          ! bare ice, thickness dependence
     372     5237712 :          fh = min(atan(hi*c4)/fhtan,c1)
     373     5237712 :          albo = albocn*(c1-fh)
     374     5237712 :          alvdfni = albicev*fh + albo
     375     5237712 :          alidfni = albicei*fh + albo
     376             : 
     377             :          ! bare ice, temperature dependence
     378     5237712 :          dTs = Timelt - Tsfcn
     379     5237712 :          fT = min(dTs/dT_melt-c1,c0)
     380     5237712 :          alvdfni = alvdfni - dalb_mlt*fT
     381     5237712 :          alidfni = alidfni - dalb_mlt*fT
     382             : 
     383             :          ! avoid negative albedos for thin, bare, melting ice
     384     5237712 :          alvdfni = max (alvdfni, albocn)
     385     5237712 :          alidfni = max (alidfni, albocn)
     386             : 
     387     5237712 :          if (hs > puny) then
     388             : 
     389     4870288 :             alvdfns = albsnowv
     390     4870288 :             alidfns = albsnowi
     391             : 
     392             :             ! snow on ice, temperature dependence
     393     4870288 :             alvdfns = alvdfns - dalb_mltv*fT
     394     4870288 :             alidfns = alidfns - dalb_mlti*fT
     395             : 
     396             :          endif                  ! hs > puny
     397             : 
     398             :          ! direct albedos (same as diffuse for now)
     399     5237712 :          alvdrni = alvdfni
     400     5237712 :          alidrni = alidfni
     401     5237712 :          alvdrns = alvdfns
     402     5237712 :          alidrns = alidfns
     403             : 
     404             :          ! fractional area of snow cover
     405     5237712 :          if (hs > puny) then
     406     4870288 :             asnow = hs / (hs + snowpatch)
     407             :          else
     408      367424 :             asnow = c0
     409             :          endif
     410             : 
     411             :          ! combine ice and snow albedos (for coupler)
     412             :          alvdfn = alvdfni*(c1-asnow) + &
     413     5237712 :                   alvdfns*asnow
     414             :          alidfn = alidfni*(c1-asnow) + &
     415     5237712 :                   alidfns*asnow
     416             :          alvdrn = alvdrni*(c1-asnow) + &
     417     5237712 :                   alvdrns*asnow
     418             :          alidrn = alidrni*(c1-asnow) + &
     419     5237712 :                   alidrns*asnow
     420             : 
     421             :          ! save ice and snow albedos (for history)
     422             :          albin = awtvdr*alvdrni + awtidr*alidrni &
     423     5237712 :                + awtvdf*alvdfni + awtidf*alidfni 
     424             :          albsn = awtvdr*alvdrns + awtidr*alidrns &
     425     5237712 :                + awtvdf*alvdfns + awtidf*alidfns 
     426             : 
     427     5237712 :       end subroutine compute_albedos
     428             : 
     429             : !=======================================================================
     430             : !
     431             : ! Compute albedos for each thickness category
     432             : 
     433    37441280 :       subroutine constant_albedos (aicen,              &
     434             :                                    vsnon,    Tsfcn,    &
     435             :                                    alvdrni,  alidrni,  &
     436             :                                    alvdfni,  alidfni,  &
     437             :                                    alvdrns,  alidrns,  &
     438             :                                    alvdfns,  alidfns,  &
     439             :                                    alvdrn,   alidrn,   &
     440             :                                    alvdfn,   alidfn,   &
     441             :                                    albin,    albsn)
     442             : 
     443             :       real (kind=dbl_kind), intent(in) :: &
     444             :          aicen   , & ! concentration of ice per category
     445             :          vsnon   , & ! volume of ice per category
     446             :          Tsfcn       ! surface temperature
     447             : 
     448             :       real (kind=dbl_kind), intent(out) :: &
     449             :          alvdrni  , & ! visible, direct, ice   (fraction)
     450             :          alidrni  , & ! near-ir, direct, ice   (fraction)
     451             :          alvdfni  , & ! visible, diffuse, ice  (fraction)
     452             :          alidfni  , & ! near-ir, diffuse, ice  (fraction)
     453             :          alvdrns  , & ! visible, direct, snow  (fraction)
     454             :          alidrns  , & ! near-ir, direct, snow  (fraction)
     455             :          alvdfns  , & ! visible, diffuse, snow (fraction)
     456             :          alidfns  , & ! near-ir, diffuse, snow (fraction)
     457             :          alvdrn   , & ! visible, direct, avg   (fraction)
     458             :          alidrn   , & ! near-ir, direct, avg   (fraction)
     459             :          alvdfn   , & ! visible, diffuse, avg  (fraction)
     460             :          alidfn   , & ! near-ir, diffuse, avg  (fraction)
     461             :          albin    , & ! bare ice 
     462             :          albsn        ! snow 
     463             : 
     464             :       ! local variables
     465             : 
     466             :       real (kind=dbl_kind), parameter :: &
     467             :          warmice  = 0.68_dbl_kind, &
     468             :          coldice  = 0.70_dbl_kind, &
     469             :          warmsnow = 0.77_dbl_kind, &
     470             :          coldsnow = 0.81_dbl_kind
     471             : 
     472             :       real (kind=dbl_kind) :: &
     473     4680160 :          hs      ! snow thickness  (m)
     474             : 
     475             :       character(len=*),parameter :: subname='(constant_albedos)'
     476             : 
     477             :       !-----------------------------------------------------------------
     478             :       ! Compute albedo for each thickness category.
     479             :       !-----------------------------------------------------------------
     480             : 
     481    37441280 :          hs = vsnon / aicen
     482             : 
     483    37441280 :          if (hs > puny) then
     484             :             ! snow, temperature dependence
     485    17955200 :             if (Tsfcn >= -c2*puny) then
     486           0 :                alvdfn = warmsnow
     487           0 :                alidfn = warmsnow
     488             :             else
     489    17955200 :                alvdfn = coldsnow
     490    17955200 :                alidfn = coldsnow
     491             :             endif
     492             :          else      ! hs < puny
     493             :             ! bare ice, temperature dependence
     494    19486080 :             if (Tsfcn >= -c2*puny) then
     495           0 :                alvdfn = warmice
     496           0 :                alidfn = warmice
     497             :             else
     498    19486080 :                alvdfn = coldice
     499    19486080 :                alidfn = coldice
     500             :             endif
     501             :          endif                  ! hs > puny
     502             : 
     503             :          ! direct albedos (same as diffuse for now)
     504    37441280 :          alvdrn  = alvdfn
     505    37441280 :          alidrn  = alidfn
     506             : 
     507    37441280 :          alvdrni = alvdrn
     508    37441280 :          alidrni = alidrn
     509    37441280 :          alvdrns = alvdrn
     510    37441280 :          alidrns = alidrn
     511    37441280 :          alvdfni = alvdfn
     512    37441280 :          alidfni = alidfn
     513    37441280 :          alvdfns = alvdfn
     514    37441280 :          alidfns = alidfn
     515             : 
     516             :          ! save ice and snow albedos (for history)
     517             :          albin = awtvdr*alvdrni + awtidr*alidrni &
     518    37441280 :                + awtvdf*alvdfni + awtidf*alidfni 
     519             :          albsn = awtvdr*alvdrns + awtidr*alidrns &
     520    37441280 :                + awtvdf*alvdfns + awtidf*alidfns 
     521             : 
     522    37441280 :       end subroutine constant_albedos
     523             : 
     524             : !=======================================================================
     525             : !
     526             : ! Compute solar radiation absorbed in ice and penetrating to ocean
     527             : !
     528             : ! authors William H. Lipscomb, LANL
     529             : !         C. M. Bitz, UW
     530             : 
     531    42678992 :       subroutine absorbed_solar (heat_capacity,      &
     532             :                                  nilyr,    aicen,    &
     533             :                                  vicen,    vsnon,    &
     534             :                                  swvdr,    swvdf,    &
     535             :                                  swidr,    swidf,    &
     536             :                                  alvdrni,  alvdfni,  &
     537             :                                  alidrni,  alidfni,  &
     538             :                                  alvdrns,  alvdfns,  &
     539             :                                  alidrns,  alidfns,  &
     540             :                                  fswsfc,   fswint,   &
     541    42678992 :                                  fswthru,  fswpenl,  &
     542    42678992 :                                  Iswabs)
     543             : 
     544             :       logical(kind=log_kind), intent(in) :: &
     545             :          heat_capacity   ! if true, ice has nonzero heat capacity
     546             : 
     547             :       integer (kind=int_kind), intent(in) :: & 
     548             :          nilyr           ! number of ice layers
     549             : 
     550             :       real (kind=dbl_kind), intent(in) :: &
     551             :          aicen       , & ! fractional ice area
     552             :          vicen       , & ! ice volume
     553             :          vsnon       , & ! snow volume
     554             :          swvdr       , & ! sw down, visible, direct  (W/m^2)
     555             :          swvdf       , & ! sw down, visible, diffuse (W/m^2)
     556             :          swidr       , & ! sw down, near IR, direct  (W/m^2)
     557             :          swidf       , & ! sw down, near IR, diffuse (W/m^2)
     558             :          alvdrni     , & ! visible, direct albedo,ice
     559             :          alidrni     , & ! near-ir, direct albedo,ice
     560             :          alvdfni     , & ! visible, diffuse albedo,ice
     561             :          alidfni     , & ! near-ir, diffuse albedo,ice
     562             :          alvdrns     , & ! visible, direct albedo, snow
     563             :          alidrns     , & ! near-ir, direct albedo, snow
     564             :          alvdfns     , & ! visible, diffuse albedo, snow
     565             :          alidfns         ! near-ir, diffuse albedo, snow
     566             : 
     567             :       real (kind=dbl_kind), intent(out):: &
     568             :          fswsfc      , & ! SW absorbed at ice/snow surface (W m-2)
     569             :          fswint      , & ! SW absorbed in ice interior, below surface (W m-2)
     570             :          fswthru         ! SW through ice to ocean (W m-2)
     571             : 
     572             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
     573             :          Iswabs      , & ! SW absorbed in particular layer (W m-2)
     574             :          fswpenl         ! visible SW entering ice layers (W m-2)
     575             : 
     576             :       ! local variables
     577             : 
     578             :       real (kind=dbl_kind), parameter :: &
     579             :          i0vis = 0.70_dbl_kind  ! fraction of penetrating solar rad (visible)
     580             : 
     581             :       integer (kind=int_kind) :: &
     582             :          k               ! ice layer index
     583             : 
     584             :       real (kind=dbl_kind) :: &
     585     5334874 :          fswpen      , & ! SW penetrating beneath surface (W m-2)
     586     5334874 :          trantop     , & ! transmitted frac of penetrating SW at layer top
     587     5334874 :          tranbot         ! transmitted frac of penetrating SW at layer bot
     588             : 
     589             :       real (kind=dbl_kind) :: &
     590     5334874 :          swabs       , & ! net SW down at surface (W m-2)
     591     5334874 :          swabsv      , & ! swabs in vis (wvlngth < 700nm)  (W/m^2)
     592     5334874 :          swabsi      , & ! swabs in nir (wvlngth > 700nm)  (W/m^2)
     593     5334874 :          fswpenvdr   , & ! penetrating SW, vis direct
     594     5334874 :          fswpenvdf   , & ! penetrating SW, vis diffuse
     595     5334874 :          hi          , & ! ice thickness (m)
     596     5334874 :          hs          , & ! snow thickness (m)
     597     5334874 :          hilyr       , & ! ice layer thickness
     598     5334874 :          asnow           ! fractional area of snow cover
     599             : 
     600             :       character(len=*),parameter :: subname='(absorbed_solar)'
     601             : 
     602             :       !-----------------------------------------------------------------
     603             :       ! Initialize
     604             :       !-----------------------------------------------------------------
     605             : 
     606    42678992 :       trantop = c0
     607    42678992 :       tranbot = c0
     608             : 
     609    42678992 :          hs  = vsnon / aicen
     610             : 
     611             :       !-----------------------------------------------------------------
     612             :       ! Fractional snow cover
     613             :       !-----------------------------------------------------------------
     614    42678992 :          if (hs > puny) then
     615    22825488 :             asnow = hs / (hs + snowpatch)
     616             :          else
     617    19853504 :             asnow = c0
     618             :          endif
     619             : 
     620             :       !-----------------------------------------------------------------
     621             :       ! Shortwave flux absorbed at surface, absorbed internally,
     622             :       !  and penetrating to mixed layer.
     623             :       ! This parameterization assumes that all IR is absorbed at the
     624             :       !  surface; only visible is absorbed in the ice interior or
     625             :       !  transmitted to the ocean.
     626             :       !-----------------------------------------------------------------
     627             : 
     628             :          swabsv  = swvdr * ( (c1-alvdrni)*(c1-asnow) &
     629             :                            + (c1-alvdrns)*asnow ) &
     630             :                  + swvdf * ( (c1-alvdfni)*(c1-asnow) &
     631    42678992 :                            + (c1-alvdfns)*asnow )
     632             : 
     633             :          swabsi  = swidr * ( (c1-alidrni)*(c1-asnow) &
     634             :                            + (c1-alidrns)*asnow ) &
     635             :                  + swidf * ( (c1-alidfni)*(c1-asnow) &
     636    42678992 :                            + (c1-alidfns)*asnow )
     637             : 
     638    42678992 :          swabs   = swabsv + swabsi
     639             : 
     640    42678992 :          fswpenvdr = swvdr * (c1-alvdrni) * (c1-asnow) * i0vis
     641    42678992 :          fswpenvdf = swvdf * (c1-alvdfni) * (c1-asnow) * i0vis
     642             : 
     643             :           ! no penetrating radiation in near IR
     644             : !         fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir
     645             : !         fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir  
     646             : 
     647    42678992 :          fswpen = fswpenvdr + fswpenvdf
     648             :                       
     649    42678992 :          fswsfc = swabs - fswpen
     650             : 
     651    42678992 :          trantop = c1  ! transmittance at top of ice
     652             : 
     653             :       !-----------------------------------------------------------------
     654             :       ! penetrating SW absorbed in each ice layer
     655             :       !-----------------------------------------------------------------
     656             : 
     657   116784256 :          do k = 1, nilyr
     658             : 
     659    74105264 :             hi  = vicen / aicen
     660    74105264 :             hilyr = hi / real(nilyr,kind=dbl_kind)
     661             : 
     662    74105264 :             tranbot = exp (-kappav * hilyr * real(k,kind=dbl_kind))
     663    74105264 :             Iswabs(k) = fswpen * (trantop-tranbot)
     664             : 
     665             :             ! bottom of layer k = top of layer k+1
     666    74105264 :             trantop = tranbot
     667             : 
     668             :             ! bgc layer model
     669   116784256 :             if (k == 1) then   ! surface flux
     670    42678992 :                fswpenl(k)   = fswpen
     671    42678992 :                fswpenl(k+1) = fswpen * tranbot
     672             :             else
     673    31426272 :                fswpenl(k+1) = fswpen * tranbot
     674             :             endif
     675             :          enddo                     ! nilyr
     676             : 
     677             :          ! SW penetrating thru ice into ocean
     678    42678992 :          fswthru = fswpen * tranbot
     679             : 
     680             :          ! SW absorbed in ice interior
     681    42678992 :          fswint  = fswpen - fswthru
     682             : 
     683             :       !----------------------------------------------------------------
     684             :       ! if zero-layer model (no heat capacity), no SW is absorbed in ice
     685             :       ! interior, so add to surface absorption
     686             :       !----------------------------------------------------------------
     687             :          
     688    42678992 :          if (.not. heat_capacity) then
     689             : 
     690             :             ! SW absorbed at snow/ice surface
     691    37441280 :             fswsfc = fswsfc + fswint
     692             : 
     693             :             ! SW absorbed in ice interior (nilyr = 1)
     694    37441280 :             fswint    = c0
     695    37441280 :             Iswabs(1) = c0
     696             : 
     697             :          endif                       ! heat_capacity
     698             : 
     699    42678992 :       end subroutine absorbed_solar
     700             : 
     701             : ! End ccsm3 shortwave method
     702             : !=======================================================================
     703             : ! Begin Delta-Eddington shortwave method
     704             : 
     705             : ! Compute initial data for Delta-Eddington method, specifically, 
     706             : ! the approximate exponential look-up table.
     707             : !
     708             : ! author:  Bruce P. Briegleb, NCAR
     709             : ! 2011 ECH modified for melt pond tracers
     710             : ! 2013 ECH merged with NCAR version
     711             : 
     712   655811942 :       subroutine run_dEdd(dt,       ncat,      &
     713             :                           dEdd_algae,          &
     714             :                           nilyr,    nslyr,     &
     715  1311623884 :                           aicen,    vicen,     &
     716   655811942 :                           vsnon,    Tsfcn,     &
     717  1311623884 :                           alvln,    apndn,     &
     718  1311623884 :                           hpndn,    ipndn,     &
     719   655811942 :                           aeron,    kalg,      &
     720   655811942 :                           trcrn_bgcsw,         &
     721             :                           heat_capacity,       &
     722             :                           tlat,     tlon,      & 
     723           0 :                           calendar_type,       &
     724             :                           days_per_year,       &
     725             :                           nextsw_cday,   yday, &
     726             :                           sec,      R_ice,     &
     727             :                           R_pnd,    R_snw,     &
     728             :                           dT_mlt,   rsnw_mlt,  &
     729             :                           hs0,      hs1,  hp1, &
     730             :                           pndaspect,           &
     731   655811942 :                           kaer_tab, waer_tab,  &
     732   655811942 :                           gaer_tab,            &
     733   655811942 :                           kaer_bc_tab,         &
     734   655811942 :                           waer_bc_tab,         &
     735   655811942 :                           gaer_bc_tab,         &
     736   655811942 :                           bcenh,               &
     737             :                           modal_aero,          &
     738             :                           swvdr,    swvdf,     &
     739             :                           swidr,    swidf,     &
     740             :                           coszen,   fsnow,     &
     741   655811942 :                           alvdrn,   alvdfn,    &
     742   655811942 :                           alidrn,   alidfn,    &
     743   655811942 :                           fswsfcn,  fswintn,   &
     744  1311623884 :                           fswthrun, fswpenln,  &
     745   655811942 :                           Sswabsn,  Iswabsn,   &
     746   655811942 :                           albicen,  albsnon,   &
     747  1311623884 :                           albpndn,  apeffn,    &
     748   655811942 :                           snowfracn,           &
     749   655811942 :                           dhsn,     ffracn,    &
     750             :                           l_print_point,       &
     751             :                           initonly)
     752             : 
     753             :       integer (kind=int_kind), intent(in) :: &
     754             :          ncat   , & ! number of ice thickness categories
     755             :          nilyr  , & ! number of ice layers
     756             :          nslyr      ! number of snow layers
     757             : 
     758             :       logical(kind=log_kind), intent(in) :: &
     759             :          heat_capacity,& ! if true, ice has nonzero heat capacity
     760             :          dEdd_algae,   & ! .true. use prognostic chla in dEdd
     761             :          modal_aero      ! .true. use modal aerosol treatment
     762             : 
     763             :       ! dEdd tuning parameters, set in namelist
     764             :       real (kind=dbl_kind), intent(in) :: &
     765             :          R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
     766             :          R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
     767             :          R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
     768             :          dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
     769             :          rsnw_mlt, & ! maximum melting snow grain radius (10^-6 m)
     770             :          hs0      , & ! snow depth for transition to bare sea ice (m)
     771             :          pndaspect, & ! ratio of pond depth to pond fraction
     772             :          hs1      , & ! tapering parameter for snow on pond ice
     773             :          hp1      , & ! critical parameter for pond ice thickness
     774             :          kalg         ! algae absorption coefficient
     775             : 
     776             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & 
     777             :          kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
     778             :          waer_tab, & ! aerosol single scatter albedo (fraction)
     779             :          gaer_tab    ! aerosol asymmetry parameter (cos(theta))
     780             :    
     781             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
     782             :          kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
     783             :          waer_bc_tab, & ! aerosol single scatter albedo (fraction)
     784             :          gaer_bc_tab    ! aerosol asymmetry parameter (cos(theta))
     785             :    
     786             :       real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
     787             :          bcenh          ! BC absorption enhancement factor
     788             : 
     789             :       character (len=char_len), intent(in) :: &
     790             :          calendar_type       ! differentiates Gregorian from other calendars
     791             : 
     792             :       integer (kind=int_kind), intent(in) :: &
     793             :          days_per_year, &    ! number of days in one year
     794             :          sec                 ! elapsed seconds into date
     795             : 
     796             :       real (kind=dbl_kind), intent(in) :: &
     797             :          nextsw_cday     , & ! julian day of next shortwave calculation
     798             :          yday                ! day of the year
     799             : 
     800             :       real(kind=dbl_kind), intent(in) :: &
     801             :            dt,    & ! time step (s)
     802             :            tlat,  & ! latitude of temp pts (radians)
     803             :            tlon,  & ! longitude of temp pts (radians)
     804             :            swvdr, & ! sw down, visible, direct  (W/m^2)
     805             :            swvdf, & ! sw down, visible, diffuse (W/m^2)
     806             :            swidr, & ! sw down, near IR, direct  (W/m^2)
     807             :            swidf, & ! sw down, near IR, diffuse (W/m^2)
     808             :            fsnow    ! snowfall rate (kg/m^2 s)
     809             : 
     810             :       real(kind=dbl_kind), dimension(:), intent(in) :: &
     811             :            aicen, & ! concentration of ice
     812             :            vicen, & ! volume per unit area of ice (m)
     813             :            vsnon, & ! volume per unit area of snow (m)
     814             :            Tsfcn, & ! surface temperature (deg C)
     815             :            alvln, & ! level-ice area fraction
     816             :            apndn, & ! pond area fraction
     817             :            hpndn, & ! pond depth (m)
     818             :            ipndn    ! pond refrozen lid thickness (m)
     819             : 
     820             :       real(kind=dbl_kind), dimension(:,:), intent(in) :: &
     821             :            aeron, & ! aerosols (kg/m^3)
     822             :            trcrn_bgcsw ! zaerosols (kg/m^3) + chlorophyll on shorthwave grid
     823             : 
     824             :       real(kind=dbl_kind), dimension(:), intent(inout) :: &
     825             :            ffracn,& ! fraction of fsurfn used to melt ipond
     826             :            dhsn     ! depth difference for snow on sea ice and pond ice
     827             : 
     828             :       real(kind=dbl_kind), intent(inout) :: &
     829             :            coszen   ! cosine solar zenith angle, < 0 for sun below horizon 
     830             : 
     831             :       real(kind=dbl_kind), dimension(:), intent(inout) :: &
     832             :            alvdrn,   & ! visible direct albedo (fraction)
     833             :            alvdfn,   & ! near-ir direct albedo (fraction)
     834             :            alidrn,   & ! visible diffuse albedo (fraction)
     835             :            alidfn,   & ! near-ir diffuse albedo (fraction)
     836             :            fswsfcn,  & ! SW absorbed at ice/snow surface (W m-2)
     837             :            fswintn,  & ! SW absorbed in ice interior, below surface (W m-2)
     838             :            fswthrun, & ! SW through ice to ocean (W/m^2) 
     839             :            albicen,  & ! albedo bare ice 
     840             :            albsnon,  & ! albedo snow 
     841             :            albpndn,  & ! albedo pond 
     842             :            apeffn,   & ! effective pond area used for radiation calculation
     843             :            snowfracn   ! snow fraction on each category used for radiation
     844             : 
     845             :       real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
     846             :            Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
     847             :            Iswabsn , & ! SW radiation absorbed in ice layers (W m-2) 
     848             :            fswpenln    ! visible SW entering ice layers (W m-2)
     849             : 
     850             :       logical (kind=log_kind), intent(in) :: &
     851             :            l_print_point
     852             : 
     853             :       logical (kind=log_kind), optional :: &
     854             :            initonly    ! flag to indicate init only, default is false
     855             : 
     856             :       ! local temporary variables
     857             : 
     858             :       ! other local variables
     859             :       ! snow variables for Delta-Eddington shortwave
     860             :       real (kind=dbl_kind) :: &
     861    36756038 :          fsn         , & ! snow horizontal fraction
     862    36756038 :          hsn             ! snow depth (m)
     863             : 
     864             :       real (kind=dbl_kind), dimension (nslyr) :: &
     865  1385135960 :          rhosnwn     , & ! snow density (kg/m3)
     866   766080056 :          rsnwn           ! snow grain radius (micrometers)
     867             : 
     868             :       ! pond variables for Delta-Eddington shortwave
     869             :       real (kind=dbl_kind) :: &
     870    36756038 :          fpn         , & ! pond fraction of ice cover
     871    36756038 :          hpn             ! actual pond depth (m)
     872             : 
     873             :       integer (kind=int_kind) :: &
     874             :          n               ! thickness category index
     875             : 
     876             :       real (kind=dbl_kind) :: &
     877    36756038 :          ipn         , & ! refrozen pond ice thickness (m), mean over ice fraction
     878    36756038 :          hp          , & ! pond depth
     879    36756038 :          hs          , & ! snow depth
     880    36756038 :          asnow       , & ! fractional area of snow cover
     881    36756038 :          rp          , & ! volume fraction of retained melt water to total liquid content
     882    36756038 :          hmx         , & ! maximum available snow infiltration equivalent depth
     883    36756038 :          dhs         , & ! local difference in snow depth on sea ice and pond ice
     884    36756038 :          spn         , & ! snow depth on refrozen pond (m)
     885    36756038 :          tmp             ! 0 or 1
     886             : 
     887             :       logical (kind=log_kind) :: &
     888             :          linitonly       ! local initonly value
     889             : 
     890             :       character(len=*),parameter :: subname='(run_dEdd)'
     891             : 
     892   655811942 :       linitonly = .false.
     893   655811942 :       if (present(initonly)) then
     894   655811942 :          linitonly = initonly
     895             :       endif
     896             : 
     897             :       ! cosine of the zenith angle
     898             : #ifdef CESMCOUPLED
     899             :       call compute_coszen (tlat,          tlon, &
     900             :                            yday,  sec, coszen,  &
     901             :                            days_per_year, nextsw_cday, calendar_type)
     902             : #else
     903             :       call compute_coszen (tlat,          tlon, &
     904   655811942 :                            yday,  sec, coszen)
     905             : #endif
     906   655811942 :       if (icepack_warnings_aborted(subname)) return
     907             : 
     908  3934871652 :       do n = 1, ncat
     909             : 
     910             :       ! note that rhoswn, rsnw, fp, hp and Sswabs ARE NOT dimensioned with ncat
     911             :       ! BPB 19 Dec 2006
     912             : 
     913             :          ! set snow properties
     914  3279059710 :          fsn        = c0
     915  3279059710 :          hsn        = c0
     916  6558119420 :          rhosnwn(:) = c0
     917  6558119420 :          rsnwn(:)   = c0
     918  3279059710 :          apeffn(n)    = c0 ! for history
     919  3279059710 :          snowfracn(n) = c0 ! for history
     920             : 
     921  3934871652 :          if (aicen(n) > puny) then
     922             : 
     923             :             call shortwave_dEdd_set_snow(nslyr,      R_snw,    &
     924             :                                          dT_mlt,     rsnw_mlt, &
     925    45664357 :                                          aicen(n),   vsnon(n), &
     926    45664357 :                                          Tsfcn(n),   fsn,      &
     927             :                                          hs0,        hsn,      &
     928   705443757 :                                          rhosnwn,    rsnwn)    
     929   705443757 :             if (icepack_warnings_aborted(subname)) return
     930             : 
     931             :             ! set pond properties
     932   705443757 :             if (tr_pond_cesm) then
     933             :                ! fraction of ice area
     934    22906571 :                fpn = apndn(n)
     935             :                ! pond depth over fraction fpn 
     936    22906571 :                hpn = hpndn(n)
     937             :                ! snow infiltration
     938    22906571 :                if (hsn >= hs_min .and. hs0 > puny) then
     939           0 :                   asnow = min(hsn/hs0, c1) ! delta-Eddington formulation
     940           0 :                   fpn = (c1 - asnow) * fpn
     941           0 :                   hpn = pndaspect * fpn
     942             :                endif
     943             :                ! Zero out fraction of thin ponds for radiation only
     944    22906571 :                if (hpn < hpmin) fpn = c0
     945    22906571 :                fsn = min(fsn, c1-fpn)
     946    22906571 :                apeffn(n) = fpn ! for history
     947   682537186 :             elseif (tr_pond_lvl) then
     948   605913432 :                fpn = c0  ! fraction of ice covered in pond
     949   605913432 :                hpn = c0  ! pond depth over fpn
     950             :                ! refrozen pond lid thickness avg over ice
     951             :                ! allow snow to cover pond ice
     952   605913432 :                ipn = alvln(n) * apndn(n) * ipndn(n)
     953   605913432 :                dhs = dhsn(n) ! snow depth difference, sea ice - pond
     954             :                if (.not. linitonly .and. ipn > puny .and. &
     955   605913432 :                     dhs < puny .and. fsnow*dt > hs_min) &
     956    10014488 :                     dhs = hsn - fsnow*dt ! initialize dhs>0
     957   605913432 :                spn = hsn - dhs   ! snow depth on pond ice
     958   605913432 :                if (.not. linitonly .and. ipn*spn < puny) dhs = c0
     959   605913432 :                dhsn(n) = dhs ! save: constant until reset to 0
     960             :                   
     961             :                ! not using ipn assumes that lid ice is perfectly clear
     962             :                ! if (ipn <= 0.3_dbl_kind) then
     963             :                
     964             :                ! fraction of ice area
     965   605913432 :                fpn = apndn(n) * alvln(n) 
     966             :                ! pond depth over fraction fpn
     967   605913432 :                hpn = hpndn(n)
     968             :                
     969             :                ! reduce effective pond area absorbing surface heat flux
     970             :                ! due to flux already having been used to melt pond ice
     971   605913432 :                fpn = (c1 - ffracn(n)) * fpn
     972             :                
     973             :                ! taper pond area with snow on pond ice
     974   605913432 :                if (dhs > puny .and. spn >= puny .and. hs1 > puny) then
     975    19477996 :                   asnow = min(spn/hs1, c1)
     976    19477996 :                   fpn = (c1 - asnow) * fpn
     977             :                endif
     978             :                
     979             :                ! infiltrate snow
     980   605913432 :                hp = hpn
     981   605913432 :                if (hp > puny) then
     982   382596642 :                   hs = hsn
     983   382596642 :                   rp = rhofresh*hp/(rhofresh*hp + rhos*hs)
     984   382596642 :                   if (rp < p15) then
     985   286207314 :                      fpn = c0
     986   286207314 :                      hpn = c0
     987             :                   else
     988    96389328 :                      hmx = hs*(rhofresh - rhos)/rhofresh
     989    96389328 :                      tmp = max(c0, sign(c1, hp-hmx)) ! 1 if hp>=hmx, else 0
     990             :                      hp = (rhofresh*hp + rhos*hs*tmp) &
     991    96389328 :                           / (rhofresh    - rhos*(c1-tmp))
     992    96389328 :                      hsn = hs - hp*fpn*(c1-tmp)
     993    96389328 :                      hpn = hp * tmp
     994    96389328 :                      fpn = fpn * tmp
     995             :                   endif
     996             :                endif ! hp > puny
     997             : 
     998             :                ! Zero out fraction of thin ponds for radiation only
     999   605913432 :                if (hpn < hpmin) fpn = c0
    1000   605913432 :                fsn = min(fsn, c1-fpn)
    1001             : 
    1002             :                ! endif    ! masking by lid ice
    1003   605913432 :                apeffn(n) = fpn ! for history
    1004             : 
    1005    76623754 :             elseif (tr_pond_topo) then
    1006             :                ! Lid effective if thicker than hp1
    1007           0 :                if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
    1008           0 :                   fpn = apndn(n)
    1009             :                else
    1010           0 :                   fpn = c0
    1011             :                endif
    1012           0 :                if (apndn(n) > puny) then
    1013           0 :                   hpn = hpndn(n) 
    1014             :                else
    1015           0 :                   fpn = c0
    1016           0 :                   hpn = c0
    1017             :                endif
    1018             : 
    1019             :                ! Zero out fraction of thin ponds for radiation only
    1020           0 :                if (hpn < hpmin) fpn = c0
    1021             : 
    1022             :                ! If ponds are present snow fraction reduced to
    1023             :                ! non-ponded part dEdd scheme 
    1024           0 :                fsn = min(fsn, c1-fpn)
    1025             : 
    1026           0 :                apeffn(n) = fpn
    1027             :             else
    1028    76623754 :                fpn = c0
    1029    76623754 :                hpn = c0
    1030    13444925 :                call shortwave_dEdd_set_pond(Tsfcn(n),   &
    1031             :                                             fsn, fpn,   &
    1032    76623754 :                                             hpn)
    1033    76623754 :                if (icepack_warnings_aborted(subname)) return
    1034             :                
    1035    76623754 :                apeffn(n) = fpn ! for history
    1036    76623754 :                fpn = c0
    1037    76623754 :                hpn = c0
    1038             :             endif ! pond type
    1039             :             
    1040   705443757 :          snowfracn(n) = fsn ! for history
    1041             : 
    1042             :          call shortwave_dEdd(dEdd_algae,                    &
    1043             :                              nslyr,         nilyr,          &
    1044             :                              coszen,        heat_capacity,  &
    1045    45664357 :                              aicen(n),      vicen(n),       &
    1046             :                              hsn,           fsn,            &
    1047           0 :                              rhosnwn,       rsnwn,          &
    1048             :                              fpn,           hpn,            &
    1049           0 :                              aeron(:,n),                    &
    1050             :                              R_ice,         R_pnd,          &
    1051           0 :                              kaer_tab,      waer_tab,       &
    1052           0 :                              gaer_tab,                      &
    1053           0 :                              kaer_bc_tab,                   &
    1054           0 :                              waer_bc_tab,                   &
    1055           0 :                              gaer_bc_tab,                   &
    1056           0 :                              bcenh,         modal_aero,     &   
    1057             :                              kalg,                          &
    1058             :                              swvdr,         swvdf,          &
    1059             :                              swidr,         swidf,          &
    1060    45664357 :                              alvdrn(n),     alvdfn(n),      &
    1061    45664357 :                              alidrn(n),     alidfn(n),      &
    1062    45664357 :                              fswsfcn(n),    fswintn(n),     &
    1063    45664357 :                              fswthrun(n),                   &
    1064           0 :                              Sswabsn(:,n),                  &
    1065           0 :                              Iswabsn(:,n),                  &
    1066    45664357 :                              albicen(n),                    &
    1067    45664357 :                              albsnon(n),    albpndn(n),     &
    1068           0 :                              fswpenln(:,n),                 &
    1069           0 :                              trcrn_bgcsw(:,n),              &
    1070   796772471 :                              l_print_point)
    1071   705443757 :          if (icepack_warnings_aborted(subname)) return
    1072             : 
    1073             :          endif ! aicen > puny
    1074             : 
    1075             :       enddo  ! ncat
    1076             : 
    1077             :       end subroutine run_dEdd
    1078             :  
    1079             : !=======================================================================
    1080             : !
    1081             : !   Compute snow/bare ice/ponded ice shortwave albedos, absorbed and transmitted 
    1082             : !   flux using the Delta-Eddington solar radiation method as described in:
    1083             : !
    1084             : !   A Delta-Eddington Multiple Scattering Parameterization for Solar Radiation
    1085             : !        in the Sea Ice Component of the Community Climate System Model
    1086             : !            B.P.Briegleb and B.Light   NCAR/TN-472+STR  February 2007
    1087             : !
    1088             : !   Compute shortwave albedos and fluxes for three surface types: 
    1089             : !   snow over ice, bare ice and ponded ice. 
    1090             : !   
    1091             : !   Albedos and fluxes are output for later use by thermodynamic routines. 
    1092             : !   Invokes three calls to compute_dEdd, which sets inherent optical properties 
    1093             : !   appropriate for the surface type. Within compute_dEdd, a call to solution_dEdd 
    1094             : !   evaluates the Delta-Eddington solution. The final albedos and fluxes are then
    1095             : !   evaluated in compute_dEdd. Albedos and fluxes are transferred to output in 
    1096             : !   this routine.
    1097             : !
    1098             : !   NOTE regarding albedo diagnostics:  This method yields zero albedo values
    1099             : !   if there is no incoming solar and thus the albedo diagnostics are masked
    1100             : !   out when the sun is below the horizon.  To estimate albedo from the history 
    1101             : !   output (post-processing), compute ice albedo using
    1102             : !   (1 - albedo)*swdn = swabs. -ECH
    1103             : !
    1104             : ! author:  Bruce P. Briegleb, NCAR 
    1105             : !   2013:  E Hunke merged with NCAR version
    1106             : !
    1107   705443757 :       subroutine shortwave_dEdd  (dEdd_algae,            &
    1108             :                                   nslyr,    nilyr,       &
    1109             :                                   coszen,   heat_capacity,&
    1110             :                                   aice,     vice,        &
    1111             :                                   hs,       fs,          & 
    1112   705443757 :                                   rhosnw,   rsnw,        &
    1113             :                                   fp,       hp,          &
    1114   705443757 :                                   aero,                  &
    1115             :                                   R_ice,    R_pnd,       &
    1116   705443757 :                                   kaer_tab, waer_tab,    &
    1117   705443757 :                                   gaer_tab,              &
    1118   705443757 :                                   kaer_bc_tab,           &
    1119   705443757 :                                   waer_bc_tab,           &
    1120   705443757 :                                   gaer_bc_tab,           &
    1121   705443757 :                                   bcenh,  modal_aero,    &   
    1122             :                                   kalg,                  &
    1123             :                                   swvdr,    swvdf,       &
    1124             :                                   swidr,    swidf,       &
    1125             :                                   alvdr,    alvdf,       &
    1126             :                                   alidr,    alidf,       &
    1127             :                                   fswsfc,   fswint,      &
    1128   705443757 :                                   fswthru,  Sswabs,      &
    1129   705443757 :                                   Iswabs,   albice,      &
    1130             :                                   albsno,   albpnd,      &
    1131  1410887514 :                                   fswpenl,  zbio,        &
    1132             :                                   l_print_point)
    1133             : 
    1134             :       integer (kind=int_kind), intent(in) :: &
    1135             :          nilyr   , & ! number of ice layers
    1136             :          nslyr       ! number of snow layers
    1137             : 
    1138             :       logical (kind=log_kind), intent(in) :: &
    1139             :          heat_capacity, & ! if true, ice has nonzero heat capacity
    1140             :          dEdd_algae,    & ! .true. use prognostic chla in dEdd
    1141             :          modal_aero       ! .true. use modal aerosol treatment
    1142             :  
    1143             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
    1144             :          kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
    1145             :          waer_bc_tab, & ! aerosol single scatter albedo (fraction)
    1146             :          gaer_bc_tab    ! aerosol asymmetry parameter (cos(theta))
    1147             :    
    1148             :       real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
    1149             :          bcenh          ! BC absorption enhancement factor
    1150             : 
    1151             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & 
    1152             :          kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
    1153             :          waer_tab, & ! aerosol single scatter albedo (fraction)
    1154             :          gaer_tab    ! aerosol asymmetry parameter (cos(theta))
    1155             : 
    1156             :       real (kind=dbl_kind), intent(in) :: &
    1157             :          kalg    , & ! algae absorption coefficient
    1158             :          R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
    1159             :          R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
    1160             :          aice    , & ! concentration of ice 
    1161             :          vice    , & ! volume of ice 
    1162             :          hs      , & ! snow depth
    1163             :          fs          ! horizontal coverage of snow
    1164             : 
    1165             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1166             :          rhosnw  , & ! density in snow layer (kg/m3)
    1167             :          rsnw    , & ! grain radius in snow layer (m)
    1168             :          aero    , & ! aerosol tracers
    1169             :          zbio        ! shortwave tracers (zaero+chla)
    1170             : 
    1171             :       real (kind=dbl_kind), intent(in) :: &
    1172             :          fp      , & ! pond fractional coverage (0 to 1) 
    1173             :          hp      , & ! pond depth (m) 
    1174             :          swvdr   , & ! sw down, visible, direct  (W/m^2)
    1175             :          swvdf   , & ! sw down, visible, diffuse (W/m^2)
    1176             :          swidr   , & ! sw down, near IR, direct  (W/m^2)
    1177             :          swidf       ! sw down, near IR, diffuse (W/m^2)
    1178             : 
    1179             :       real (kind=dbl_kind), intent(inout) :: &
    1180             :          coszen  , & ! cosine of solar zenith angle 
    1181             :          alvdr   , & ! visible, direct, albedo (fraction) 
    1182             :          alvdf   , & ! visible, diffuse, albedo (fraction) 
    1183             :          alidr   , & ! near-ir, direct, albedo (fraction) 
    1184             :          alidf   , & ! near-ir, diffuse, albedo (fraction) 
    1185             :          fswsfc  , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
    1186             :          fswint  , & ! SW interior absorption (below surface, above ocean,W m-2)
    1187             :          fswthru     ! SW through snow/bare ice/ponded ice into ocean (W m-2)
    1188             :  
    1189             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1190             :          fswpenl , & ! visible SW entering ice layers (W m-2)
    1191             :          Sswabs  , & ! SW absorbed in snow layer (W m-2)
    1192             :          Iswabs      ! SW absorbed in ice layer (W m-2)
    1193             : 
    1194             :       real (kind=dbl_kind), intent(out) :: &
    1195             :          albice  , & ! bare ice albedo, for history  
    1196             :          albsno  , & ! snow albedo, for history  
    1197             :          albpnd      ! pond albedo, for history  
    1198             : 
    1199             :       logical (kind=log_kind) , intent(in) :: &
    1200             :          l_print_point
    1201             : 
    1202             :       ! local variables
    1203             : 
    1204             :       real (kind=dbl_kind) :: &
    1205    45664357 :          netsw    , & ! net shortwave
    1206    45664357 :          fnidr    , & ! fraction of direct to total down surface flux in nir
    1207    45664357 :          hstmp    , & ! snow thickness (set to 0 for bare ice case)
    1208    45664357 :          hi       , & ! ice thickness (all sea ice layers, m)
    1209    45664357 :          fi           ! snow/bare ice fractional coverage (0 to 1)
    1210             : 
    1211             :       real (kind=dbl_kind), dimension (4*n_aero) :: &
    1212   933089790 :          aero_mp      ! aerosol mass path in kg/m2
    1213             : 
    1214             :       integer (kind=int_kind) :: &
    1215             :          srftyp       ! surface type over ice: (0=air, 1=snow, 2=pond)
    1216             :  
    1217             :       integer (kind=int_kind) :: &
    1218             :          k        , & ! level index
    1219             :          na       , & ! aerosol index
    1220             :          klev     , & ! number of radiation layers - 1
    1221             :          klevp        ! number of radiation interfaces - 1
    1222             :                       ! (0 layer is included also)
    1223             : 
    1224             :       real (kind=dbl_kind) :: &
    1225    45664357 :          vsno         ! volume of snow 
    1226             : 
    1227             :       real (kind=dbl_kind) :: &
    1228    45664357 :          swdn  , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j)
    1229    45664357 :          swab  , & ! fswsfc(i,j)+fswint(i,j)+fswthru(i,j)
    1230    45664357 :          swalb     ! (1.-swab/(swdn+.0001))
    1231             : 
    1232             :       ! for history
    1233             :       real (kind=dbl_kind) :: &
    1234    45664357 :          avdrl   , & ! visible, direct, albedo (fraction) 
    1235    45664357 :          avdfl   , & ! visible, diffuse, albedo (fraction) 
    1236    45664357 :          aidrl   , & ! near-ir, direct, albedo (fraction) 
    1237    45664357 :          aidfl       ! near-ir, diffuse, albedo (fraction) 
    1238             : 
    1239             :       character(len=*),parameter :: subname='(shortwave_dEdd)'
    1240             : 
    1241             : !-----------------------------------------------------------------------
    1242             : 
    1243   705443757 :          klev    = nslyr + nilyr + 1   ! number of radiation layers - 1
    1244   705443757 :          klevp   = klev  + 1           ! number of radiation interfaces - 1
    1245             :                                        ! (0 layer is included also)
    1246             : 
    1247             :       ! zero storage albedos and fluxes for accumulation over surface types:
    1248   705443757 :       hstmp    = c0
    1249   705443757 :       hi       = c0
    1250   705443757 :       fi       = c0
    1251   705443757 :       alvdr    = c0
    1252   705443757 :       alvdf    = c0
    1253   705443757 :       alidr    = c0
    1254   705443757 :       alidf    = c0
    1255   705443757 :       avdrl    = c0
    1256   705443757 :       avdfl    = c0
    1257   705443757 :       aidrl    = c0
    1258   705443757 :       aidfl    = c0
    1259   705443757 :       fswsfc   = c0
    1260   705443757 :       fswint   = c0
    1261   705443757 :       fswthru  = c0
    1262             :       ! compute fraction of nir down direct to total over all points:
    1263   705443757 :       fnidr = c0
    1264   705443757 :       if( swidr + swidf > puny ) then
    1265   335099540 :          fnidr = swidr/(swidr+swidf)
    1266             :       endif
    1267   705443757 :       albice    = c0
    1268   705443757 :       albsno    = c0
    1269   705443757 :       albpnd    = c0
    1270  5889318405 :       fswpenl(:) = c0
    1271  1410887514 :       Sswabs(:)  = c0
    1272  5183874648 :       Iswabs(:)  = c0
    1273             : 
    1274             :       ! compute aerosol mass path
    1275             :       
    1276  3291009201 :          aero_mp(:) = c0
    1277   705443757 :          if( tr_aero ) then
    1278             :             ! check 4 layers for each aerosol, a snow SSL, snow below SSL,
    1279             :             ! sea ice SSL, and sea ice below SSL, in that order.
    1280    22372076 :             if (size(aero) < 4*n_aero) then
    1281           0 :                call icepack_warnings_add(subname//' ERROR: size(aero) too small')
    1282           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1283           0 :                return
    1284             :             endif
    1285    44744152 :             do na = 1, 4*n_aero, 4
    1286    22372076 :                vsno = hs * aice
    1287    22372076 :                netsw = swvdr + swidr + swvdf + swidf
    1288    44744152 :                if (netsw > puny) then ! sun above horizon
    1289    15245921 :                   aero_mp(na  ) = aero(na  )*vsno
    1290    15245921 :                   aero_mp(na+1) = aero(na+1)*vsno
    1291    15245921 :                   aero_mp(na+2) = aero(na+2)*vice
    1292    15245921 :                   aero_mp(na+3) = aero(na+3)*vice
    1293             :                endif                  ! aice > 0 and netsw > 0
    1294             :             enddo      ! na
    1295             :          endif      ! if aerosols
    1296             : 
    1297             :          ! compute shortwave radiation accounting for snow/ice (both snow over 
    1298             :          ! ice and bare ice) and ponded ice (if any):
    1299             :          
    1300             :          ! sea ice points with sun above horizon
    1301   705443757 :          netsw = swvdr + swidr + swvdf + swidf
    1302   705443757 :          if (netsw > puny) then ! sun above horizon
    1303   335099540 :             coszen = max(puny,coszen)
    1304             :             ! evaluate sea ice thickness and fraction
    1305   335099540 :             hi  = vice / aice
    1306   335099540 :             fi  = c1 - fs - fp
    1307             :             ! bare sea ice points
    1308   335099540 :             if(fi > c0) then
    1309             :                ! calculate bare sea ice
    1310             : 
    1311    49079912 :                srftyp = 0
    1312             :                call compute_dEdd(nilyr,       nslyr,   klev,   klevp,   & 
    1313           0 :                       zbio,      dEdd_algae,                            &
    1314             :                       heat_capacity,          fnidr,   coszen,          &
    1315             :                       R_ice,     R_pnd,                                 &
    1316           0 :                       kaer_tab,    waer_tab,    gaer_tab,               &
    1317           0 :                       kaer_bc_tab, waer_bc_tab, gaer_bc_tab,            &
    1318           0 :                       bcenh,     modal_aero,  kalg,                     &
    1319             :                       swvdr,     swvdf,       swidr,   swidf,  srftyp,  &
    1320           0 :                       hstmp,     rhosnw,      rsnw,    hi,     hp,      &
    1321           0 :                       fi,        aero_mp,     avdrl,   avdfl,           &
    1322             :                       aidrl,     aidfl,                                 &
    1323             :                       fswsfc,    fswint,                                &
    1324           0 :                       fswthru,   Sswabs,                                &
    1325    49079912 :                       Iswabs,    fswpenl)
    1326    49079912 :                if (icepack_warnings_aborted(subname)) return
    1327             :                
    1328    49079912 :                alvdr   = alvdr   + avdrl *fi
    1329    49079912 :                alvdf   = alvdf   + avdfl *fi
    1330    49079912 :                alidr   = alidr   + aidrl *fi
    1331    49079912 :                alidf   = alidf   + aidfl *fi
    1332             :                ! for history
    1333             :                albice = albice &
    1334             :                       + awtvdr*avdrl + awtidr*aidrl &
    1335    49079912 :                       + awtvdf*avdfl + awtidf*aidfl 
    1336             :             endif
    1337             :          endif
    1338             :          
    1339             :          ! sea ice points with sun above horizon
    1340   705443757 :          netsw = swvdr + swidr + swvdf + swidf
    1341   705443757 :          if (netsw > puny) then ! sun above horizon
    1342   335099540 :             coszen = max(puny,coszen)
    1343             :             ! snow-covered sea ice points
    1344   335099540 :             if(fs > c0) then
    1345             :                ! calculate snow covered sea ice
    1346             : 
    1347   292471682 :                srftyp = 1
    1348             :                call compute_dEdd(nilyr,       nslyr,   klev,   klevp,   & 
    1349           0 :                       zbio,      dEdd_algae,                            &
    1350             :                       heat_capacity,          fnidr,   coszen,          &
    1351             :                       R_ice,     R_pnd,                                 &
    1352           0 :                       kaer_tab,    waer_tab,    gaer_tab,               &
    1353           0 :                       kaer_bc_tab, waer_bc_tab, gaer_bc_tab,            &
    1354           0 :                       bcenh,     modal_aero,  kalg,                     &
    1355             :                       swvdr,     swvdf,       swidr,   swidf,  srftyp,  &
    1356           0 :                       hs,        rhosnw,      rsnw,    hi,     hp,      &
    1357           0 :                       fs,        aero_mp,     avdrl,   avdfl,           &
    1358             :                       aidrl,     aidfl,                                 &
    1359             :                       fswsfc,    fswint,                                &
    1360           0 :                       fswthru,   Sswabs,                                &
    1361   292471682 :                       Iswabs,    fswpenl)
    1362   292471682 :                if (icepack_warnings_aborted(subname)) return
    1363             :                
    1364   292471682 :                alvdr   = alvdr   + avdrl *fs
    1365   292471682 :                alvdf   = alvdf   + avdfl *fs
    1366   292471682 :                alidr   = alidr   + aidrl *fs
    1367   292471682 :                alidf   = alidf   + aidfl *fs
    1368             :                ! for history
    1369             :                albsno = albsno &
    1370             :                       + awtvdr*avdrl + awtidr*aidrl &
    1371   292471682 :                       + awtvdf*avdfl + awtidf*aidfl 
    1372             :             endif
    1373             :          endif
    1374             :          
    1375   705443757 :          hi = c0
    1376             : 
    1377             :          ! sea ice points with sun above horizon
    1378   705443757 :          netsw = swvdr + swidr + swvdf + swidf
    1379   705443757 :          if (netsw > puny) then ! sun above horizon
    1380   335099540 :             coszen = max(puny,coszen)
    1381   335099540 :             hi  = vice / aice
    1382             :             ! if nonzero pond fraction and sufficient pond depth
    1383             :             ! if( fp > puny .and. hp > hpmin ) then
    1384   335099540 :             if (fp > puny) then
    1385             :                
    1386             :                ! calculate ponded ice
    1387             : 
    1388    44295577 :                srftyp = 2
    1389             :                call compute_dEdd(nilyr,       nslyr,   klev,   klevp,   & 
    1390           0 :                       zbio,      dEdd_algae,                            &
    1391             :                       heat_capacity,          fnidr,   coszen,          &
    1392             :                       R_ice,     R_pnd,                                 &
    1393           0 :                       kaer_tab,    waer_tab,    gaer_tab,               &
    1394           0 :                       kaer_bc_tab, waer_bc_tab, gaer_bc_tab,            &
    1395           0 :                       bcenh,     modal_aero,  kalg,                     &
    1396             :                       swvdr,     swvdf,       swidr,   swidf,  srftyp,  & 
    1397           0 :                       hs,        rhosnw,      rsnw,    hi,     hp,      &
    1398           0 :                       fp,        aero_mp,     avdrl,   avdfl,           &
    1399             :                       aidrl,     aidfl,                                 &
    1400             :                       fswsfc,    fswint,                                &
    1401           0 :                       fswthru,   Sswabs,                                &
    1402    44295577 :                       Iswabs,    fswpenl)
    1403    44295577 :                if (icepack_warnings_aborted(subname)) return
    1404             :                
    1405    44295577 :                alvdr   = alvdr   + avdrl *fp
    1406    44295577 :                alvdf   = alvdf   + avdfl *fp
    1407    44295577 :                alidr   = alidr   + aidrl *fp
    1408    44295577 :                alidf   = alidf   + aidfl *fp
    1409             :                ! for history
    1410             :                albpnd = albpnd &
    1411             :                       + awtvdr*avdrl + awtidr*aidrl &
    1412    44295577 :                       + awtvdf*avdfl + awtidf*aidfl 
    1413             :             endif
    1414             :          endif
    1415             : 
    1416             :          ! if no incoming shortwave, set albedos to 1
    1417   705443757 :          netsw = swvdr + swidr + swvdf + swidf
    1418   705443757 :          if (netsw <= puny) then ! sun above horizon
    1419   370344217 :             alvdr = c1
    1420   370344217 :             alvdf = c1
    1421   370344217 :             alidr = c1
    1422   370344217 :             alidf = c1
    1423             :          endif
    1424             : 
    1425   705443757 :       if (l_print_point .and. netsw > puny) then
    1426             : 
    1427           0 :          write(warnstr,*) subname, ' printing point'
    1428           0 :          call icepack_warnings_add(warnstr)
    1429           0 :          write(warnstr,*) subname, ' coszen = ', &
    1430           0 :                             coszen
    1431           0 :          call icepack_warnings_add(warnstr)
    1432           0 :          write(warnstr,*) subname, ' swvdr  swvdf = ', &
    1433           0 :                             swvdr,swvdf
    1434           0 :          call icepack_warnings_add(warnstr)
    1435           0 :          write(warnstr,*) subname, ' swidr  swidf = ', &
    1436           0 :                             swidr,swidf
    1437           0 :          call icepack_warnings_add(warnstr)
    1438           0 :          write(warnstr,*) subname, ' aice = ', &
    1439           0 :                             aice
    1440           0 :          call icepack_warnings_add(warnstr)
    1441           0 :          write(warnstr,*) subname, ' hs = ', &
    1442           0 :                             hs
    1443           0 :          call icepack_warnings_add(warnstr)
    1444           0 :          write(warnstr,*) subname, ' hp = ', &
    1445           0 :                             hp
    1446           0 :          call icepack_warnings_add(warnstr)
    1447           0 :          write(warnstr,*) subname, ' fs = ', &
    1448           0 :                             fs
    1449           0 :          call icepack_warnings_add(warnstr)
    1450           0 :          write(warnstr,*) subname, ' fi = ', &
    1451           0 :                             fi
    1452           0 :          call icepack_warnings_add(warnstr)
    1453           0 :          write(warnstr,*) subname, ' fp = ', &
    1454           0 :                             fp
    1455           0 :          call icepack_warnings_add(warnstr)
    1456           0 :          write(warnstr,*) subname, ' hi = ', &
    1457           0 :                             hi
    1458           0 :          call icepack_warnings_add(warnstr)
    1459           0 :          write(warnstr,*) subname, ' alvdr  alvdf = ', &
    1460           0 :                             alvdr,alvdf
    1461           0 :          call icepack_warnings_add(warnstr)
    1462           0 :          write(warnstr,*) subname, ' alidr  alidf = ', &
    1463           0 :                             alidr,alidf
    1464           0 :          call icepack_warnings_add(warnstr)
    1465           0 :          write(warnstr,*) subname, ' fswsfc fswint fswthru = ', &
    1466           0 :                             fswsfc,fswint,fswthru
    1467           0 :          call icepack_warnings_add(warnstr)
    1468           0 :          swdn  = swvdr+swvdf+swidr+swidf
    1469           0 :          swab  = fswsfc+fswint+fswthru
    1470           0 :          swalb = (1.-swab/(swdn+.0001))
    1471           0 :          write(warnstr,*) subname, ' swdn swab swalb = ',swdn,swab,swalb
    1472           0 :          do k = 1, nslyr               
    1473           0 :             write(warnstr,*) subname, ' snow layer k    = ', k, &
    1474           0 :                              ' rhosnw = ', &
    1475           0 :                                rhosnw(k), &
    1476           0 :                              ' rsnw = ', &
    1477           0 :                                rsnw(k)
    1478           0 :             call icepack_warnings_add(warnstr)
    1479             :          enddo
    1480           0 :          do k = 1, nslyr               
    1481           0 :             write(warnstr,*) subname, ' snow layer k    = ', k, &
    1482           0 :                              ' Sswabs(k)       = ', Sswabs(k)
    1483           0 :             call icepack_warnings_add(warnstr)
    1484             :          enddo
    1485           0 :          do k = 1, nilyr               
    1486           0 :             write(warnstr,*) subname, ' sea ice layer k = ', k, &
    1487           0 :                              ' Iswabs(k)       = ', Iswabs(k)
    1488           0 :             call icepack_warnings_add(warnstr)
    1489             :          enddo
    1490             : 
    1491             :       endif  ! l_print_point .and. coszen > .01
    1492             : 
    1493             :       end subroutine shortwave_dEdd
    1494             : 
    1495             : !=======================================================================
    1496             : !
    1497             : ! Evaluate snow/ice/ponded ice inherent optical properties (IOPs), and 
    1498             : ! then calculate the multiple scattering solution by calling solution_dEdd.
    1499             : !
    1500             : ! author:  Bruce P. Briegleb, NCAR 
    1501             : !   2013:  E Hunke merged with NCAR version
    1502             : 
    1503   385847171 :       subroutine compute_dEdd (nilyr,    nslyr,    klev,  klevp,  &
    1504   385847171 :                     zbio,     dEdd_algae,                         &
    1505             :                     heat_capacity,       fnidr,    coszen,        &
    1506             :                     R_ice,     R_pnd,                             &
    1507   385847171 :                     kaer_tab,    waer_tab,         gaer_tab,      &
    1508   385847171 :                     kaer_bc_tab, waer_bc_tab,      gaer_bc_tab,   &
    1509   385847171 :                     bcenh,     modal_aero,         kalg,          &
    1510             :                     swvdr,     swvdf,    swidr,    swidf, srftyp, &
    1511   771694342 :                     hs,        rhosnw,   rsnw,     hi,    hp,     &
    1512   385847171 :                     fi,        aero_mp,  alvdr,    alvdf,         &
    1513             :                                alidr,    alidf,         &
    1514             :                                fswsfc,   fswint,        &
    1515   385847171 :                                fswthru,  Sswabs,        &
    1516   385847171 :                                Iswabs,   fswpenl)
    1517             : 
    1518             :       integer (kind=int_kind), intent(in) :: &
    1519             :          nilyr , & ! number of ice layers
    1520             :          nslyr , & ! number of snow layers
    1521             :          klev  , & ! number of radiation layers - 1
    1522             :          klevp     ! number of radiation interfaces - 1
    1523             :                    ! (0 layer is included also)
    1524             :  
    1525             :       logical (kind=log_kind), intent(in) :: &
    1526             :          heat_capacity,& ! if true, ice has nonzero heat capacity
    1527             :          dEdd_algae,   & ! .true. use prognostic chla in dEdd
    1528             :          modal_aero      ! .true. use modal aerosol treatment
    1529             :  
    1530             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
    1531             :          kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
    1532             :          waer_bc_tab, & ! aerosol single scatter albedo (fraction)
    1533             :          gaer_bc_tab    ! aerosol asymmetry parameter (cos(theta))
    1534             :    
    1535             :       real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
    1536             :          bcenh          ! BC absorption enhancement factor
    1537             : 
    1538             :       ! dEdd tuning parameters, set in namelist
    1539             :       real (kind=dbl_kind), intent(in) :: &
    1540             :          R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
    1541             :          R_pnd     ! ponded ice tuning parameter; +1 > 1sig increase in albedo
    1542             : 
    1543             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & 
    1544             :          kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
    1545             :          waer_tab, & ! aerosol single scatter albedo (fraction)
    1546             :          gaer_tab    ! aerosol asymmetry parameter (cos(theta))
    1547             :    
    1548             :       real (kind=dbl_kind), intent(in) :: &
    1549             :          kalg    , & ! algae absorption coefficient
    1550             :          fnidr   , & ! fraction of direct to total down flux in nir
    1551             :          coszen  , & ! cosine solar zenith angle
    1552             :          swvdr   , & ! shortwave down at surface, visible, direct  (W/m^2)
    1553             :          swvdf   , & ! shortwave down at surface, visible, diffuse (W/m^2)
    1554             :          swidr   , & ! shortwave down at surface, near IR, direct  (W/m^2)
    1555             :          swidf       ! shortwave down at surface, near IR, diffuse (W/m^2)
    1556             :  
    1557             :       integer (kind=int_kind), intent(in) :: &
    1558             :          srftyp      ! surface type over ice: (0=air, 1=snow, 2=pond)
    1559             :  
    1560             :       real (kind=dbl_kind), intent(in) :: &
    1561             :          hs          ! snow thickness (m)
    1562             : 
    1563             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1564             :          rhosnw  , & ! snow density in snow layer (kg/m3)
    1565             :          rsnw    , & ! snow grain radius in snow layer (m)
    1566             :          zbio    , & ! zaerosol + chla shortwave tracers kg/m^3
    1567             :          aero_mp     ! aerosol mass path in kg/m2
    1568             : 
    1569             :       real (kind=dbl_kind), intent(in) :: &
    1570             :          hi      , & ! ice thickness (m)
    1571             :          hp      , & ! pond depth (m)
    1572             :          fi          ! snow/bare ice fractional coverage (0 to 1)
    1573             :  
    1574             :       real (kind=dbl_kind), intent(inout) :: &
    1575             :          alvdr   , & ! visible, direct, albedo (fraction) 
    1576             :          alvdf   , & ! visible, diffuse, albedo (fraction) 
    1577             :          alidr   , & ! near-ir, direct, albedo (fraction) 
    1578             :          alidf   , & ! near-ir, diffuse, albedo (fraction) 
    1579             :          fswsfc  , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
    1580             :          fswint  , & ! SW interior absorption (below surface, above ocean,W m-2)
    1581             :          fswthru     ! SW through snow/bare ice/ponded ice into ocean (W m-2)
    1582             :  
    1583             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1584             :          fswpenl , & ! visible SW entering ice layers (W m-2)
    1585             :          Sswabs  , & ! SW absorbed in snow layer (W m-2)
    1586             :          Iswabs      ! SW absorbed in ice layer (W m-2)
    1587             : 
    1588             : !-----------------------------------------------------------------------
    1589             : !
    1590             : ! Set up optical property profiles, based on snow, sea ice and ponded 
    1591             : ! ice IOPs from:
    1592             : !
    1593             : ! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple 
    1594             : !    Scattering Parameterization for Solar Radiation in the Sea Ice 
    1595             : !    Component of the Community Climate System Model, NCAR Technical 
    1596             : !    Note  NCAR/TN-472+STR  February 2007
    1597             : !
    1598             : ! Computes column Delta-Eddington radiation solution for specific
    1599             : ! surface type: either snow over sea ice, bare sea ice, or ponded sea ice.
    1600             : !
    1601             : ! Divides solar spectrum into 3 intervals: 0.2-0.7, 0.7-1.19, and
    1602             : ! 1.19-5.0 micro-meters. The latter two are added (using an assumed
    1603             : ! partition of incident shortwave in the 0.7-5.0 micro-meter band between
    1604             : ! the 0.7-1.19 and 1.19-5.0 micro-meter band) to give the final output 
    1605             : ! of 0.2-0.7 visible and 0.7-5.0 near-infrared albedos and fluxes.
    1606             : !
    1607             : ! Specifies vertical layer optical properties based on input snow depth,
    1608             : ! density and grain radius, along with ice and pond depths, then computes
    1609             : ! layer by layer Delta-Eddington reflectivity, transmissivity and combines
    1610             : ! layers (done by calling routine solution_dEdd). Finally, surface albedos
    1611             : ! and internal fluxes/flux divergences are evaluated.
    1612             : !
    1613             : !  Description of the level and layer index conventions. This is
    1614             : !  for the standard case of one snow layer and four sea ice layers.
    1615             : !
    1616             : !  Please read the following; otherwise, there is 99.9% chance you
    1617             : !  will be confused about indices at some point in time........ :)
    1618             : !
    1619             : !  CICE4.0 snow treatment has one snow layer above the sea ice. This 
    1620             : !  snow layer has finite heat capacity, so that surface absorption must
    1621             : !  be distinguished from internal. The Delta-Eddington solar radiation
    1622             : !  thus adds extra surface scattering layers to both snow and sea ice.
    1623             : !  Note that in the following, we assume a fixed vertical layer structure
    1624             : !  for the radiation calculation. In other words, we always have the 
    1625             : !  structure shown below for one snow and four sea ice layers, but for 
    1626             : !  ponded ice the pond fills "snow" layer 1 over the sea ice, and for 
    1627             : !  bare sea ice the top layers over sea ice are treated as transparent air.
    1628             : !
    1629             : !  SSL = surface scattering layer for either snow or sea ice
    1630             : !  DL  = drained layer for sea ice immediately under sea ice SSL
    1631             : !  INT = interior layers for sea ice below the drained layer.
    1632             : !
    1633             : !  Notice that the radiation level starts with 0 at the top. Thus,
    1634             : !  the total number radiation layers is klev+1, where klev is the
    1635             : !  sum of nslyr, the number of CCSM snow layers, and nilyr, the
    1636             : !  number of CCSM sea ice layers, plus the sea ice SSL:
    1637             : !  klev = 1 + nslyr + nilyr
    1638             : !
    1639             : !  For the standard case illustrated below, nslyr=1, nilyr=4,
    1640             : !  and klev=6, with the number of layer interfaces klevp=klev+1.
    1641             : !  Layer interfaces are the surfaces on which reflectivities,
    1642             : !  transmissivities and fluxes are evaluated.
    1643             : !
    1644             : !  CCSM3 Sea Ice Model            Delta-Eddington Solar Radiation
    1645             : !                                     Layers and Interfaces
    1646             : !                             Layer Index             Interface Index
    1647             : !    ---------------------            ---------------------  0
    1648             : !                                  0  \\\   snow SSL    \\\
    1649             : !       snow layer 1                  ---------------------  1
    1650             : !                                  1    rest of snow layer
    1651             : !    +++++++++++++++++++++            +++++++++++++++++++++  2
    1652             : !                                  2  \\\ sea ice SSL   \\\
    1653             : !      sea ice layer 1                ---------------------  3
    1654             : !                                  3      sea ice  DL
    1655             : !    ---------------------            ---------------------  4
    1656             : !
    1657             : !      sea ice layer 2             4      sea ice INT
    1658             : !
    1659             : !    ---------------------            ---------------------  5
    1660             : !
    1661             : !      sea ice layer 3             5      sea ice INT
    1662             : !
    1663             : !    ---------------------            ---------------------  6
    1664             : !
    1665             : !      sea ice layer 4             6      sea ice INT
    1666             : !
    1667             : !    ---------------------            ---------------------  7
    1668             : !
    1669             : ! When snow lies over sea ice, the radiation absorbed in the
    1670             : ! snow SSL is used for surface heating, and that in the rest
    1671             : ! of the snow layer for its internal heating. For sea ice in
    1672             : ! this case, all of the radiant heat absorbed in both the
    1673             : ! sea ice SSL and the DL are used for sea ice layer 1 heating.
    1674             : !
    1675             : ! When pond lies over sea ice, and for bare sea ice, all of the
    1676             : ! radiant heat absorbed within and above the sea ice SSL is used
    1677             : ! for surface heating, and that absorbed in the sea ice DL is
    1678             : ! used for sea ice layer 1 heating.
    1679             : !
    1680             : ! Basically, vertical profiles of the layer extinction optical depth (tau), 
    1681             : ! single scattering albedo (w0) and asymmetry parameter (g) are required over
    1682             : ! the klev+1 layers, where klev+1 = 2 + nslyr + nilyr. All of the surface type
    1683             : ! information and snow/ice iop properties are evaulated in this routine, so
    1684             : ! the tau,w0,g profiles can be passed to solution_dEdd for multiple scattering
    1685             : ! evaluation. Snow, bare ice and ponded ice iops are contained in data arrays
    1686             : ! in this routine.
    1687             : !
    1688             : !-----------------------------------------------------------------------
    1689             : 
    1690             :       ! local variables
    1691             : 
    1692             :       integer (kind=int_kind) :: &
    1693             :          k       , & ! level index
    1694             :          ns      , & ! spectral index
    1695             :          nr      , & ! index for grain radius tables
    1696             :          ki      , & ! index for internal absorption
    1697             :          km      , & ! k starting index for snow, sea ice internal absorption
    1698             :          kp      , & ! k+1 or k+2 index for snow, sea ice internal absorption
    1699             :          ksrf    , & ! level index for surface absorption
    1700             :          ksnow   , & ! level index for snow density and grain size
    1701             :          kii         ! level starting index for sea ice (nslyr+1)
    1702             : 
    1703             :       integer (kind=int_kind), parameter :: & 
    1704             :          nmbrad  = 32        ! number of snow grain radii in tables
    1705             :  
    1706             :       real (kind=dbl_kind) :: & 
    1707    24277435 :          avdr    , & ! visible albedo, direct   (fraction)
    1708    24277435 :          avdf    , & ! visible albedo, diffuse  (fraction)
    1709    24277435 :          aidr    , & ! near-ir albedo, direct   (fraction)
    1710    24277435 :          aidf        ! near-ir albedo, diffuse  (fraction)
    1711             :  
    1712             :       real (kind=dbl_kind) :: & 
    1713    24277435 :          fsfc    , & ! shortwave absorbed at snow/bare ice/ponded ice surface (W m-2)
    1714    24277435 :          fint    , & ! shortwave absorbed in interior (W m-2)
    1715    24277435 :          fthru       ! shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
    1716             : 
    1717             :       real (kind=dbl_kind), dimension(nslyr) :: & 
    1718   771694342 :          Sabs        ! shortwave absorbed in snow layer (W m-2)
    1719             : 
    1720             :       real (kind=dbl_kind), dimension(nilyr) :: & 
    1721   917358952 :          Iabs        ! shortwave absorbed in ice layer (W m-2)
    1722             :  
    1723             :       real (kind=dbl_kind), dimension(nilyr+1) :: & 
    1724   941636387 :          fthrul      ! shortwave through to ice layers (W m-2)
    1725             : 
    1726             :       real (kind=dbl_kind), dimension (nspint) :: &
    1727    97109740 :          wghtns              ! spectral weights
    1728             :  
    1729             :       real (kind=dbl_kind), parameter :: & 
    1730             :          cp67    = 0.67_dbl_kind   , & ! nir band weight parameter
    1731             :          cp78    = 0.78_dbl_kind   , & ! nir band weight parameter
    1732             :          cp01    = 0.01_dbl_kind       ! for ocean visible albedo
    1733             :  
    1734             :       real (kind=dbl_kind), dimension (0:klev) :: &
    1735   990191257 :          tau     , & ! layer extinction optical depth
    1736   990191257 :          w0      , & ! layer single scattering albedo
    1737   990191257 :          g           ! layer asymmetry parameter
    1738             :  
    1739             :       ! following arrays are defined at model interfaces; 0 is the top of the
    1740             :       ! layer above the sea ice; klevp is the sea ice/ocean interface.
    1741             :       real (kind=dbl_kind), dimension (0:klevp) :: &
    1742  1014468692 :          trndir  , & ! solar beam down transmission from top
    1743  1014468692 :          trntdr  , & ! total transmission to direct beam for layers above
    1744  1014468692 :          trndif  , & ! diffuse transmission to diffuse beam for layers above
    1745  1014468692 :          rupdir  , & ! reflectivity to direct radiation for layers below
    1746  1014468692 :          rupdif  , & ! reflectivity to diffuse radiation for layers below
    1747  1014468692 :          rdndif      ! reflectivity to diffuse radiation for layers above
    1748             :  
    1749             :       real (kind=dbl_kind), dimension (0:klevp) :: &
    1750  1014468692 :          dfdir   , & ! down-up flux at interface due to direct beam at top surface
    1751  1063023562 :          dfdif       ! down-up flux at interface due to diffuse beam at top surface
    1752             : 
    1753             :       real (kind=dbl_kind) :: &
    1754    24277435 :          refk    , & ! interface k multiple scattering term
    1755    24277435 :          delr    , & ! snow grain radius interpolation parameter
    1756             :       ! inherent optical properties (iop) for snow
    1757    24277435 :          Qs      , & ! Snow extinction efficiency
    1758    24277435 :          ks      , & ! Snow extinction coefficient (/m)
    1759    24277435 :          ws      , & ! Snow single scattering albedo
    1760    24277435 :          gs          ! Snow asymmetry parameter
    1761             : 
    1762             :       real (kind=dbl_kind), dimension(nslyr) :: & 
    1763   771694342 :          frsnw       ! snow grain radius in snow layer * adjustment factor (m)
    1764             : 
    1765             :       ! actual used ice and ponded ice IOPs, allowing for tuning 
    1766             :       ! modifications of the above "_mn" value
    1767             :       real (kind=dbl_kind), dimension (nspint) :: &
    1768    97109740 :          ki_ssl       , & ! Surface-scattering-layer ice extinction coefficient (/m)
    1769    97109740 :          wi_ssl       , & ! Surface-scattering-layer ice single scattering albedo
    1770    97109740 :          gi_ssl       , & ! Surface-scattering-layer ice asymmetry parameter
    1771    97109740 :          ki_dl        , & ! Drained-layer ice extinction coefficient (/m)
    1772    97109740 :          wi_dl        , & ! Drained-layer ice single scattering albedo
    1773    97109740 :          gi_dl        , & ! Drained-layer ice asymmetry parameter
    1774    97109740 :          ki_int       , & ! Interior-layer ice extinction coefficient (/m)
    1775    97109740 :          wi_int       , & ! Interior-layer ice single scattering albedo
    1776    97109740 :          gi_int       , & ! Interior-layer ice asymmetry parameter
    1777    97109740 :          ki_p_ssl     , & ! Ice under pond srf scat layer extinction coefficient (/m)
    1778    97109740 :          wi_p_ssl     , & ! Ice under pond srf scat layer single scattering albedo
    1779    97109740 :          gi_p_ssl     , & ! Ice under pond srf scat layer asymmetry parameter
    1780    97109740 :          ki_p_int     , & ! Ice under pond extinction coefficient (/m)
    1781    97109740 :          wi_p_int     , & ! Ice under pond single scattering albedo
    1782    97109740 :          gi_p_int         ! Ice under pond asymmetry parameter
    1783             : 
    1784             :       real (kind=dbl_kind), dimension(0:klev) :: &
    1785   990191257 :          dzk              ! layer thickness
    1786             : 
    1787             :       real (kind=dbl_kind) :: &
    1788    24277435 :          dz           , & ! snow, sea ice or pond water layer thickness
    1789    24277435 :          dz_ssl       , & ! snow or sea ice surface scattering layer thickness
    1790    24277435 :          fs               ! scaling factor to reduce (nilyr<4) or increase (nilyr>4) DL
    1791             :                           ! extinction coefficient to maintain DL optical depth constant
    1792             :                           ! with changing number of sea ice layers, to approximately 
    1793             :                           ! conserve computed albedo for constant physical depth of sea
    1794             :                           ! ice when the number of sea ice layers vary
    1795             :       real (kind=dbl_kind) :: &
    1796    24277435 :          sig          , & ! scattering coefficient for tuning
    1797    24277435 :          kabs         , & ! absorption coefficient for tuning
    1798    24277435 :          sigp             ! modified scattering coefficient for tuning
    1799             : 
    1800             :       real (kind=dbl_kind), dimension(nspint, 0:klev) :: &
    1801  1718514307 :          kabs_chl    , & ! absorption coefficient for chlorophyll (/m)
    1802  1718514307 :          tzaer       , & ! total aerosol extinction optical depth
    1803  1718514307 :          wzaer       , & ! total aerosol single scatter albedo
    1804  1718514307 :          gzaer           ! total aerosol asymmetry parameter
    1805             : 
    1806             :       real (kind=dbl_kind) :: &
    1807    24277435 :          albodr       , & ! spectral ocean albedo to direct rad
    1808    24277435 :          albodf           ! spectral ocean albedo to diffuse rad
    1809             :       
    1810             :       ! for melt pond transition to bare sea ice for small pond depths 
    1811             :       real (kind=dbl_kind) :: &
    1812    24277435 :          sig_i        , & ! ice scattering coefficient (/m)
    1813    24277435 :          sig_p        , & ! pond scattering coefficient (/m)
    1814    24277435 :          kext             ! weighted extinction coefficient (/m)
    1815             : 
    1816             :       ! aerosol optical properties from Mark Flanner, 26 June 2008
    1817             :       ! order assumed: hydrophobic black carbon, hydrophilic black carbon,
    1818             :       ! four dust aerosols by particle size range:
    1819             :       ! dust1(.05-0.5 micron), dust2(0.5-1.25 micron),
    1820             :       ! dust3(1.25-2.5 micron), dust4(2.5-5.0 micron)
    1821             :       ! spectral bands same as snow/sea ice: (0.3-0.7 micron, 0.7-1.19 micron
    1822             :       ! and 1.19-5.0 micron in wavelength)
    1823             : 
    1824             :       integer (kind=int_kind) :: &
    1825             :          na , n                        ! aerosol index
    1826             : 
    1827             :       real (kind=dbl_kind) :: &
    1828    24277435 :          taer                     , & ! total aerosol extinction optical depth
    1829    24277435 :          waer                     , & ! total aerosol single scatter albedo
    1830    24277435 :          gaer                     , & ! total aerosol asymmetry parameter
    1831    24277435 :          swdr                     , & ! shortwave down at surface, direct  (W/m^2)
    1832    24277435 :          swdf                     , & ! shortwave down at surface, diffuse (W/m^2)
    1833    24277435 :          rnilyr                   , & ! real(nilyr)
    1834    24277435 :          rnslyr                   , & ! real(nslyr)
    1835    24277435 :          rns                      , & ! real(ns)
    1836    48554870 :          tmp_0, tmp_ks, tmp_kl        ! temp variables
    1837             : 
    1838             :       integer(kind=int_kind), dimension(0:klev) :: &
    1839   771694342 :          k_bcini        , & !
    1840   771694342 :          k_bcins        , &
    1841   410124606 :          k_bcexs
    1842             :       real(kind=dbl_kind)::  &
    1843    24277435 :           tmp_gs, tmp1               ! temp variables
    1844             : 
    1845             :       ! snow grain radii (micro-meters) for table
    1846             :       real (kind=dbl_kind), dimension(nmbrad), parameter :: &
    1847             :          rsnw_tab = (/ &   ! snow grain radius for each table entry (micro-meters)
    1848             :           5._dbl_kind,    7._dbl_kind,   10._dbl_kind,   15._dbl_kind, &
    1849             :          20._dbl_kind,   30._dbl_kind,   40._dbl_kind,   50._dbl_kind, &
    1850             :          65._dbl_kind,   80._dbl_kind,  100._dbl_kind,  120._dbl_kind, &
    1851             :         140._dbl_kind,  170._dbl_kind,  200._dbl_kind,  240._dbl_kind, &
    1852             :         290._dbl_kind,  350._dbl_kind,  420._dbl_kind,  500._dbl_kind, &
    1853             :         570._dbl_kind,  660._dbl_kind,  760._dbl_kind,  870._dbl_kind, &
    1854             :        1000._dbl_kind, 1100._dbl_kind, 1250._dbl_kind, 1400._dbl_kind, &
    1855             :        1600._dbl_kind, 1800._dbl_kind, 2000._dbl_kind, 2500._dbl_kind/)
    1856             : 
    1857             :       ! snow extinction efficiency (unitless)
    1858             :       real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
    1859             :          Qs_tab = reshape((/ &
    1860             :           2.131798_dbl_kind,  2.187756_dbl_kind,  2.267358_dbl_kind, &
    1861             :           2.104499_dbl_kind,  2.148345_dbl_kind,  2.236078_dbl_kind, &
    1862             :           2.081580_dbl_kind,  2.116885_dbl_kind,  2.175067_dbl_kind, &
    1863             :           2.062595_dbl_kind,  2.088937_dbl_kind,  2.130242_dbl_kind, &
    1864             :           2.051403_dbl_kind,  2.072422_dbl_kind,  2.106610_dbl_kind, &
    1865             :           2.039223_dbl_kind,  2.055389_dbl_kind,  2.080586_dbl_kind, &
    1866             :           2.032383_dbl_kind,  2.045751_dbl_kind,  2.066394_dbl_kind, &
    1867             :           2.027920_dbl_kind,  2.039388_dbl_kind,  2.057224_dbl_kind, &
    1868             :           2.023444_dbl_kind,  2.033137_dbl_kind,  2.048055_dbl_kind, &
    1869             :           2.020412_dbl_kind,  2.028840_dbl_kind,  2.041874_dbl_kind, &
    1870             :           2.017608_dbl_kind,  2.024863_dbl_kind,  2.036046_dbl_kind, &
    1871             :           2.015592_dbl_kind,  2.022021_dbl_kind,  2.031954_dbl_kind, &
    1872             :           2.014083_dbl_kind,  2.019887_dbl_kind,  2.028853_dbl_kind, &
    1873             :           2.012368_dbl_kind,  2.017471_dbl_kind,  2.025353_dbl_kind, &
    1874             :           2.011092_dbl_kind,  2.015675_dbl_kind,  2.022759_dbl_kind, &
    1875             :           2.009837_dbl_kind,  2.013897_dbl_kind,  2.020168_dbl_kind, &
    1876             :           2.008668_dbl_kind,  2.012252_dbl_kind,  2.017781_dbl_kind, &
    1877             :           2.007627_dbl_kind,  2.010813_dbl_kind,  2.015678_dbl_kind, &
    1878             :           2.006764_dbl_kind,  2.009577_dbl_kind,  2.013880_dbl_kind, &
    1879             :           2.006037_dbl_kind,  2.008520_dbl_kind,  2.012382_dbl_kind, &
    1880             :           2.005528_dbl_kind,  2.007807_dbl_kind,  2.011307_dbl_kind, &
    1881             :           2.005025_dbl_kind,  2.007079_dbl_kind,  2.010280_dbl_kind, &
    1882             :           2.004562_dbl_kind,  2.006440_dbl_kind,  2.009333_dbl_kind, &
    1883             :           2.004155_dbl_kind,  2.005898_dbl_kind,  2.008523_dbl_kind, &
    1884             :           2.003794_dbl_kind,  2.005379_dbl_kind,  2.007795_dbl_kind, &
    1885             :           2.003555_dbl_kind,  2.005041_dbl_kind,  2.007329_dbl_kind, &
    1886             :           2.003264_dbl_kind,  2.004624_dbl_kind,  2.006729_dbl_kind, &
    1887             :           2.003037_dbl_kind,  2.004291_dbl_kind,  2.006230_dbl_kind, &
    1888             :           2.002776_dbl_kind,  2.003929_dbl_kind,  2.005700_dbl_kind, &
    1889             :           2.002590_dbl_kind,  2.003627_dbl_kind,  2.005276_dbl_kind, &
    1890             :           2.002395_dbl_kind,  2.003391_dbl_kind,  2.004904_dbl_kind, &
    1891             :           2.002071_dbl_kind,  2.002922_dbl_kind,  2.004241_dbl_kind/), &
    1892             :           (/nspint,nmbrad/))
    1893             : 
    1894             :       ! snow single scattering albedo (unitless)
    1895             :       real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
    1896             :         ws_tab = reshape((/ &
    1897             :          0.9999994_dbl_kind,  0.9999673_dbl_kind,  0.9954589_dbl_kind, &
    1898             :          0.9999992_dbl_kind,  0.9999547_dbl_kind,  0.9938576_dbl_kind, &
    1899             :          0.9999990_dbl_kind,  0.9999382_dbl_kind,  0.9917989_dbl_kind, &
    1900             :          0.9999985_dbl_kind,  0.9999123_dbl_kind,  0.9889724_dbl_kind, &
    1901             :          0.9999979_dbl_kind,  0.9998844_dbl_kind,  0.9866190_dbl_kind, &
    1902             :          0.9999970_dbl_kind,  0.9998317_dbl_kind,  0.9823021_dbl_kind, &
    1903             :          0.9999960_dbl_kind,  0.9997800_dbl_kind,  0.9785269_dbl_kind, &
    1904             :          0.9999951_dbl_kind,  0.9997288_dbl_kind,  0.9751601_dbl_kind, &
    1905             :          0.9999936_dbl_kind,  0.9996531_dbl_kind,  0.9706974_dbl_kind, &
    1906             :          0.9999922_dbl_kind,  0.9995783_dbl_kind,  0.9667577_dbl_kind, &
    1907             :          0.9999903_dbl_kind,  0.9994798_dbl_kind,  0.9621007_dbl_kind, &
    1908             :          0.9999885_dbl_kind,  0.9993825_dbl_kind,  0.9579541_dbl_kind, &
    1909             :          0.9999866_dbl_kind,  0.9992862_dbl_kind,  0.9541924_dbl_kind, &
    1910             :          0.9999838_dbl_kind,  0.9991434_dbl_kind,  0.9490959_dbl_kind, &
    1911             :          0.9999810_dbl_kind,  0.9990025_dbl_kind,  0.9444940_dbl_kind, &
    1912             :          0.9999772_dbl_kind,  0.9988171_dbl_kind,  0.9389141_dbl_kind, &
    1913             :          0.9999726_dbl_kind,  0.9985890_dbl_kind,  0.9325819_dbl_kind, &
    1914             :          0.9999670_dbl_kind,  0.9983199_dbl_kind,  0.9256405_dbl_kind, &
    1915             :          0.9999605_dbl_kind,  0.9980117_dbl_kind,  0.9181533_dbl_kind, &
    1916             :          0.9999530_dbl_kind,  0.9976663_dbl_kind,  0.9101540_dbl_kind, &
    1917             :          0.9999465_dbl_kind,  0.9973693_dbl_kind,  0.9035031_dbl_kind, &
    1918             :          0.9999382_dbl_kind,  0.9969939_dbl_kind,  0.8953134_dbl_kind, &
    1919             :          0.9999289_dbl_kind,  0.9965848_dbl_kind,  0.8865789_dbl_kind, &
    1920             :          0.9999188_dbl_kind,  0.9961434_dbl_kind,  0.8773350_dbl_kind, &
    1921             :          0.9999068_dbl_kind,  0.9956323_dbl_kind,  0.8668233_dbl_kind, &
    1922             :          0.9998975_dbl_kind,  0.9952464_dbl_kind,  0.8589990_dbl_kind, &
    1923             :          0.9998837_dbl_kind,  0.9946782_dbl_kind,  0.8476493_dbl_kind, &
    1924             :          0.9998699_dbl_kind,  0.9941218_dbl_kind,  0.8367318_dbl_kind, &
    1925             :          0.9998515_dbl_kind,  0.9933966_dbl_kind,  0.8227881_dbl_kind, &
    1926             :          0.9998332_dbl_kind,  0.9926888_dbl_kind,  0.8095131_dbl_kind, &
    1927             :          0.9998148_dbl_kind,  0.9919968_dbl_kind,  0.7968620_dbl_kind, &
    1928             :          0.9997691_dbl_kind,  0.9903277_dbl_kind,  0.7677887_dbl_kind/), &
    1929             :          (/nspint,nmbrad/))
    1930             : 
    1931             :       ! snow asymmetry parameter (unitless)
    1932             :       real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
    1933             :          gs_tab = reshape((/ &
    1934             :           0.859913_dbl_kind,  0.848003_dbl_kind,  0.824415_dbl_kind, &
    1935             :           0.867130_dbl_kind,  0.858150_dbl_kind,  0.848445_dbl_kind, &
    1936             :           0.873381_dbl_kind,  0.867221_dbl_kind,  0.861714_dbl_kind, &
    1937             :           0.878368_dbl_kind,  0.874879_dbl_kind,  0.874036_dbl_kind, &
    1938             :           0.881462_dbl_kind,  0.879661_dbl_kind,  0.881299_dbl_kind, &
    1939             :           0.884361_dbl_kind,  0.883903_dbl_kind,  0.890184_dbl_kind, &
    1940             :           0.885937_dbl_kind,  0.886256_dbl_kind,  0.895393_dbl_kind, &
    1941             :           0.886931_dbl_kind,  0.887769_dbl_kind,  0.899072_dbl_kind, &
    1942             :           0.887894_dbl_kind,  0.889255_dbl_kind,  0.903285_dbl_kind, &
    1943             :           0.888515_dbl_kind,  0.890236_dbl_kind,  0.906588_dbl_kind, &
    1944             :           0.889073_dbl_kind,  0.891127_dbl_kind,  0.910152_dbl_kind, &
    1945             :           0.889452_dbl_kind,  0.891750_dbl_kind,  0.913100_dbl_kind, &
    1946             :           0.889730_dbl_kind,  0.892213_dbl_kind,  0.915621_dbl_kind, &
    1947             :           0.890026_dbl_kind,  0.892723_dbl_kind,  0.918831_dbl_kind, &
    1948             :           0.890238_dbl_kind,  0.893099_dbl_kind,  0.921540_dbl_kind, &
    1949             :           0.890441_dbl_kind,  0.893474_dbl_kind,  0.924581_dbl_kind, &
    1950             :           0.890618_dbl_kind,  0.893816_dbl_kind,  0.927701_dbl_kind, &
    1951             :           0.890762_dbl_kind,  0.894123_dbl_kind,  0.930737_dbl_kind, &
    1952             :           0.890881_dbl_kind,  0.894397_dbl_kind,  0.933568_dbl_kind, &
    1953             :           0.890975_dbl_kind,  0.894645_dbl_kind,  0.936148_dbl_kind, &
    1954             :           0.891035_dbl_kind,  0.894822_dbl_kind,  0.937989_dbl_kind, &
    1955             :           0.891097_dbl_kind,  0.895020_dbl_kind,  0.939949_dbl_kind, &
    1956             :           0.891147_dbl_kind,  0.895212_dbl_kind,  0.941727_dbl_kind, &
    1957             :           0.891189_dbl_kind,  0.895399_dbl_kind,  0.943339_dbl_kind, &
    1958             :           0.891225_dbl_kind,  0.895601_dbl_kind,  0.944915_dbl_kind, &
    1959             :           0.891248_dbl_kind,  0.895745_dbl_kind,  0.945950_dbl_kind, &
    1960             :           0.891277_dbl_kind,  0.895951_dbl_kind,  0.947288_dbl_kind, &
    1961             :           0.891299_dbl_kind,  0.896142_dbl_kind,  0.948438_dbl_kind, &
    1962             :           0.891323_dbl_kind,  0.896388_dbl_kind,  0.949762_dbl_kind, &
    1963             :           0.891340_dbl_kind,  0.896623_dbl_kind,  0.950916_dbl_kind, &
    1964             :           0.891356_dbl_kind,  0.896851_dbl_kind,  0.951945_dbl_kind, &
    1965             :           0.891386_dbl_kind,  0.897399_dbl_kind,  0.954156_dbl_kind/), &
    1966             :           (/nspint,nmbrad/))
    1967             : 
    1968             :       ! inherent optical property (iop) arrays for ice and ponded ice
    1969             :       ! mn = specified mean (or base) value
    1970             :       ! ki = extinction coefficient (/m)
    1971             :       ! wi = single scattering albedo
    1972             :       ! gi = asymmetry parameter
    1973             : 
    1974             :       ! ice surface scattering layer (ssl) iops
    1975             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    1976             :          ki_ssl_mn = (/ 1000.1_dbl_kind, 1003.7_dbl_kind, 7042._dbl_kind/), &
    1977             :          wi_ssl_mn = (/ .9999_dbl_kind,  .9963_dbl_kind,  .9088_dbl_kind/), &
    1978             :          gi_ssl_mn = (/  .94_dbl_kind,     .94_dbl_kind,    .94_dbl_kind/)
    1979             : 
    1980             :       ! ice drained layer (dl) iops
    1981             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    1982             :          ki_dl_mn = (/ 100.2_dbl_kind, 107.7_dbl_kind,  1309._dbl_kind /), &
    1983             :          wi_dl_mn = (/ .9980_dbl_kind,  .9287_dbl_kind, .0305_dbl_kind /), &
    1984             :          gi_dl_mn = (/ .94_dbl_kind,     .94_dbl_kind,    .94_dbl_kind /)
    1985             : 
    1986             :       ! ice interior layer (int) iops
    1987             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    1988             :          ki_int_mn = (/  20.2_dbl_kind,  27.7_dbl_kind,  1445._dbl_kind /), &
    1989             :          wi_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind,  .0277_dbl_kind /), &
    1990             :          gi_int_mn = (/ .94_dbl_kind,    .94_dbl_kind,     .94_dbl_kind /)
    1991             : 
    1992             :       ! ponded ice surface scattering layer (ssl) iops
    1993             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    1994             :          ki_p_ssl_mn = (/ 70.2_dbl_kind,  77.7_dbl_kind,  1309._dbl_kind/), &
    1995             :          wi_p_ssl_mn = (/ .9972_dbl_kind, .9009_dbl_kind, .0305_dbl_kind/), &
    1996             :          gi_p_ssl_mn = (/ .94_dbl_kind,   .94_dbl_kind,   .94_dbl_kind  /)
    1997             : 
    1998             :       ! ponded ice interior layer (int) iops
    1999             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    2000             :          ki_p_int_mn = (/  20.2_dbl_kind,  27.7_dbl_kind, 1445._dbl_kind/), &
    2001             :          wi_p_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind/), &
    2002             :          gi_p_int_mn = (/ .94_dbl_kind,   .94_dbl_kind,   .94_dbl_kind  /)
    2003             : 
    2004             :       ! inherent optical property (iop) arrays for pond water and underlying ocean
    2005             :       ! kw = Pond water extinction coefficient (/m)
    2006             :       ! ww = Pond water single scattering albedo
    2007             :       ! gw = Pond water asymmetry parameter
    2008             :       real (kind=dbl_kind), dimension (nspint), parameter :: &
    2009             :          kw = (/ 0.20_dbl_kind,   12.0_dbl_kind,   729._dbl_kind /), &
    2010             :          ww = (/ 0.00_dbl_kind,   0.00_dbl_kind,   0.00_dbl_kind /), &
    2011             :          gw = (/ 0.00_dbl_kind,   0.00_dbl_kind,   0.00_dbl_kind /)
    2012             : 
    2013             :       real (kind=dbl_kind), parameter :: &
    2014             :          rhoi   = 917.0_dbl_kind,& ! pure ice mass density (kg/m3)
    2015             :          fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max
    2016             :          fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min
    2017             :       ! tuning parameters
    2018             :       ! ice and pond scat coeff fractional change for +- one-sigma in albedo
    2019             :          fp_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for + stn dev in alb
    2020             :          fm_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for - stn dev in alb
    2021             :          fp_pnd = 2.00_dbl_kind, & ! ponded ice fraction of scat coeff for + stn dev in alb
    2022             :          fm_pnd = 0.50_dbl_kind    ! ponded ice fraction of scat coeff for - stn dev in alb
    2023             : 
    2024             :       real (kind=dbl_kind),  parameter :: &   !chla-specific absorption coefficient
    2025             :          kchl_tab = 0.01 !0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1
    2026             :                          ! found values of 0.006 to 0.023 m^2/ mg  (676 nm)  Neukermans 2014
    2027             :                          ! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011)
    2028             :                          ! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014)
    2029             :                          ! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting)
    2030             :                          !chlorophyll mass extinction cross section (m^2/mg chla)
    2031             : 
    2032             :       character(len=*),parameter :: subname='(compute_dEdd)'
    2033             : 
    2034             : !-----------------------------------------------------------------------
    2035             : ! Initialize and tune bare ice/ponded ice iops
    2036             : 
    2037  4244318881 :       k_bcini(:) = c0
    2038  4244318881 :       k_bcins(:) = c0
    2039  4244318881 :       k_bcexs(:) = c0
    2040             : 
    2041   385847171 :       rnilyr = c1/real(nilyr,kind=dbl_kind)
    2042   385847171 :       rnslyr = c1/real(nslyr,kind=dbl_kind)
    2043   385847171 :       kii = nslyr + 1
    2044             : 
    2045             :       ! initialize albedos and fluxes to 0
    2046  3472624539 :       fthrul = c0                
    2047  3086777368 :       Iabs = c0
    2048 15819734011 :       kabs_chl(:,:) = c0
    2049 15819734011 :       tzaer(:,:) = c0
    2050 15819734011 :       wzaer(:,:) = c0
    2051 15819734011 :       gzaer(:,:) = c0
    2052             : 
    2053   385847171 :       avdr   = c0
    2054   385847171 :       avdf   = c0
    2055   385847171 :       aidr   = c0
    2056   385847171 :       aidf   = c0
    2057   385847171 :       fsfc   = c0
    2058   385847171 :       fint   = c0
    2059   385847171 :       fthru  = c0
    2060             :  
    2061             :       ! spectral weights
    2062             :       ! weights 2 (0.7-1.19 micro-meters) and 3 (1.19-5.0 micro-meters) 
    2063             :       ! are chosen based on 1D calculations using ratio of direct to total 
    2064             :       ! near-infrared solar (0.7-5.0 micro-meter) which indicates clear/cloudy 
    2065             :       ! conditions: more cloud, the less 1.19-5.0 relative to the 
    2066             :       ! 0.7-1.19 micro-meter due to cloud absorption.
    2067   385847171 :       wghtns(1) = c1
    2068   385847171 :       wghtns(2) = cp67 + (cp78-cp67)*(c1-fnidr)
    2069             : !      wghtns(3) = cp33 + (cp22-cp33)*(c1-fnidr)
    2070   385847171 :       wghtns(3) = c1 - wghtns(2)
    2071             : 
    2072             :       ! find snow grain adjustment factor, dependent upon clear/overcast sky
    2073             :       ! estimate. comparisons with SNICAR show better agreement with DE when
    2074             :       ! this factor is included (clear sky near 1 and overcast near 0.8 give
    2075             :       ! best agreement).  Multiply by rnsw here for efficiency.
    2076   771694342 :       do k = 1, nslyr
    2077   385847171 :          frsnw(k) = (fr_max*fnidr + fr_min*(c1-fnidr))*rsnw(k)
    2078   771694342 :          Sabs(k) = c0
    2079             :       enddo
    2080             : 
    2081             :       ! layer thicknesses
    2082             :       ! snow
    2083   385847171 :       dz = hs*rnslyr
    2084             :       ! for small enough snow thickness, ssl thickness half of top snow layer
    2085             : !ech: note this is highly resolution dependent!
    2086   385847171 :       dzk(0) = min(hs_ssl, dz/c2)
    2087   385847171 :       dzk(1) = dz - dzk(0)
    2088   385847171 :       if (nslyr > 1) then
    2089           0 :          do k = 2, nslyr
    2090           0 :             dzk(k) = dz
    2091             :          enddo
    2092             :       endif
    2093             :       
    2094             :       ! ice
    2095   385847171 :       dz = hi*rnilyr
    2096             :       ! empirical reduction in sea ice ssl thickness for ice thinner than 1.5m;
    2097             :       ! factor of 30 gives best albedo comparison with limited observations
    2098   385847171 :       dz_ssl = hi_ssl
    2099             : !ech: note hardwired parameters
    2100             : !         if( hi < 1.5_dbl_kind ) dz_ssl = hi/30._dbl_kind
    2101   385847171 :       dz_ssl = min(hi_ssl, hi/30._dbl_kind)
    2102             :       ! set sea ice ssl thickness to half top layer if sea ice thin enough
    2103             : !ech: note this is highly resolution dependent!
    2104   385847171 :       dz_ssl = min(dz_ssl, dz/c2)
    2105             :       
    2106   385847171 :       dzk(kii)   = dz_ssl
    2107   385847171 :       dzk(kii+1) = dz - dz_ssl
    2108   385847171 :       if (kii+2 <= klev) then
    2109  2700930197 :          do k = kii+2, klev
    2110  2700930197 :             dzk(k) = dz
    2111             :          enddo
    2112             :       endif
    2113             : 
    2114             :       ! adjust sea ice iops with tuning parameters; tune only the
    2115             :       ! scattering coefficient by factors of R_ice, R_pnd, where
    2116             :       ! R values of +1 correspond approximately to +1 sigma changes in albedo, and
    2117             :       ! R values of -1 correspond approximately to -1 sigma changes in albedo
    2118             :       ! Note: the albedo change becomes non-linear for R values > +1 or < -1
    2119   385847171 :       if( R_ice >= c0 ) then
    2120  1543388684 :         do ns = 1, nspint
    2121  1157541513 :           sigp       = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fp_ice*R_ice)
    2122  1157541513 :           ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
    2123  1157541513 :           wi_ssl(ns) = sigp/ki_ssl(ns)
    2124  1157541513 :           gi_ssl(ns) = gi_ssl_mn(ns)
    2125             : 
    2126  1157541513 :           sigp       = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fp_ice*R_ice)
    2127  1157541513 :           ki_dl(ns)  = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
    2128  1157541513 :           wi_dl(ns)  = sigp/ki_dl(ns)
    2129  1157541513 :           gi_dl(ns)  = gi_dl_mn(ns)
    2130             : 
    2131  1157541513 :           sigp       = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fp_ice*R_ice)
    2132  1157541513 :           ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
    2133  1157541513 :           wi_int(ns) = sigp/ki_int(ns)
    2134  1543388684 :           gi_int(ns) = gi_int_mn(ns)
    2135             :         enddo
    2136             :       else !if( R_ice < c0 ) then
    2137           0 :         do ns = 1, nspint
    2138           0 :           sigp       = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fm_ice*R_ice)
    2139           0 :           sigp       = max(sigp, c0)
    2140           0 :           ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
    2141           0 :           wi_ssl(ns) = sigp/ki_ssl(ns)
    2142           0 :           gi_ssl(ns) = gi_ssl_mn(ns)
    2143             : 
    2144           0 :           sigp       = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fm_ice*R_ice)
    2145           0 :           sigp       = max(sigp, c0)
    2146           0 :           ki_dl(ns)  = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
    2147           0 :           wi_dl(ns)  = sigp/ki_dl(ns)
    2148           0 :           gi_dl(ns)  = gi_dl_mn(ns)
    2149             : 
    2150           0 :           sigp       = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fm_ice*R_ice)
    2151           0 :           sigp       = max(sigp, c0)
    2152           0 :           ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
    2153           0 :           wi_int(ns) = sigp/ki_int(ns)
    2154           0 :           gi_int(ns) = gi_int_mn(ns)
    2155             :         enddo
    2156             :       endif          ! adjust ice iops
    2157             : 
    2158             :       ! adjust ponded ice iops with tuning parameters
    2159   385847171 :       if( R_pnd >= c0 ) then
    2160  1543388684 :         do ns = 1, nspint
    2161  1157541513 :           sigp         = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fp_pnd*R_pnd)
    2162  1157541513 :           ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
    2163  1157541513 :           wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
    2164  1157541513 :           gi_p_ssl(ns) = gi_p_ssl_mn(ns)
    2165             : 
    2166  1157541513 :           sigp         = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fp_pnd*R_pnd)
    2167  1157541513 :           ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
    2168  1157541513 :           wi_p_int(ns) = sigp/ki_p_int(ns)
    2169  1543388684 :           gi_p_int(ns) = gi_p_int_mn(ns)
    2170             :         enddo
    2171             :       else !if( R_pnd < c0 ) then
    2172           0 :         do ns = 1, nspint
    2173           0 :           sigp         = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fm_pnd*R_pnd)
    2174           0 :           sigp         = max(sigp, c0)
    2175           0 :           ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
    2176           0 :           wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
    2177           0 :           gi_p_ssl(ns) = gi_p_ssl_mn(ns)
    2178             : 
    2179           0 :           sigp         = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fm_pnd*R_pnd)
    2180           0 :           sigp         = max(sigp, c0)
    2181           0 :           ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
    2182           0 :           wi_p_int(ns) = sigp/ki_p_int(ns)
    2183           0 :           gi_p_int(ns) = gi_p_int_mn(ns)
    2184             :         enddo
    2185             :       endif            ! adjust ponded ice iops
    2186             : 
    2187             :       ! use srftyp to determine interface index of surface absorption
    2188   385847171 :       if (srftyp == 1) then
    2189             :          ! snow covered sea ice
    2190   292471682 :          ksrf = 1
    2191             :       else
    2192             :          ! bare sea ice or ponded ice
    2193    93375489 :          ksrf = nslyr + 2 
    2194             :       endif
    2195             : 
    2196   385847171 :       if (tr_bgc_N .and. dEdd_algae) then ! compute kabs_chl for chlorophyll
    2197           0 :           do k = 0, klev
    2198           0 :              kabs_chl(1,k) = kchl_tab*zbio(nlt_chl_sw+k)
    2199             :           enddo
    2200             :       else
    2201   385847171 :             k = klev
    2202   385847171 :             kabs_chl(1,k) = kalg*(0.50_dbl_kind/dzk(k)) 
    2203             :       endif        ! kabs_chl
    2204             : 
    2205             : !mgf++
    2206   385847171 :       if (modal_aero) then
    2207           0 :            do k=0,klev   
    2208           0 :               if (k < nslyr+1) then ! define indices for snow layer
    2209             :                  ! use top rsnw, rhosnw for snow ssl and rest of top layer
    2210           0 :                  ksnow = k - min(k-1,0)
    2211           0 :                  tmp_gs = frsnw(ksnow)
    2212             :                 
    2213             :                  ! get grain size index:
    2214             :                  ! works for 25 < snw_rds < 1625 um:
    2215           0 :                  if (tmp_gs < 125) then
    2216           0 :                    tmp1 = tmp_gs/50
    2217           0 :                    k_bcini(k) = nint(tmp1)
    2218           0 :                  elseif (tmp_gs < 175) then
    2219           0 :                    k_bcini(k) = 2
    2220             :                  else
    2221           0 :                    tmp1 = (tmp_gs/250)+2
    2222           0 :                    k_bcini(k) = nint(tmp1)
    2223             :                  endif
    2224             :               else                  ! use the largest snow grain size for ice
    2225           0 :                  k_bcini(k) = 8
    2226             :               endif
    2227             :               ! Set index corresponding to BC effective radius.  Here,
    2228             :               ! asssume constant BC effective radius of 100nm
    2229             :               ! (corresponding to index 2)
    2230           0 :               k_bcins(k) = 2
    2231           0 :               k_bcexs(k) = 2
    2232             : 
    2233             :               ! check bounds:
    2234           0 :               if (k_bcini(k) < 1)  k_bcini(k) = 1
    2235           0 :               if (k_bcini(k) > 8)  k_bcini(k) = 8
    2236           0 :               if (k_bcins(k) < 1)  k_bcins(k) = 1
    2237           0 :               if (k_bcins(k) > 10) k_bcins(k) = 10
    2238           0 :               if (k_bcexs(k) < 1)  k_bcexs(k) = 1
    2239           0 :               if (k_bcexs(k) > 10) k_bcexs(k) = 10
    2240             : 
    2241             :               ! print ice radius index:
    2242             :               ! write(warnstr,*) subname, "MGFICE2:k, ice index= ",k,  k_bcini(k)
    2243             :               ! call icepack_warnings_add(warnstr)
    2244             :             enddo   ! k
    2245             : 
    2246           0 :         if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll  
    2247           0 :         do n = 1,n_zaero        
    2248           0 :            if (n == 1) then ! interstitial BC
    2249           0 :             do k = 0, klev
    2250           0 :                do ns = 1,nspint   ! not weighted by aice
    2251           0 :                   tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
    2252           0 :                                 zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2253           0 :                   wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
    2254           0 :                                 waer_bc_tab(ns,k_bcexs(k))* &
    2255           0 :                                 zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2256           0 :                   gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
    2257           0 :                                 waer_bc_tab(ns,k_bcexs(k))* &
    2258           0 :                                 gaer_bc_tab(ns,k_bcexs(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2259             :                enddo  ! nspint
    2260             :             enddo
    2261           0 :            elseif (n==2) then ! within-ice BC
    2262           0 :             do k = 0, klev
    2263           0 :                do ns = 1,nspint  
    2264           0 :                   tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))  * &
    2265           0 :                                 bcenh(ns,k_bcins(k),k_bcini(k))* &
    2266           0 :                                 zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2267           0 :                   wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
    2268           0 :                                 waer_bc_tab(ns,k_bcins(k))* &
    2269           0 :                                 zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2270           0 :                   gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
    2271           0 :                                 waer_bc_tab(ns,k_bcins(k))* &
    2272           0 :                                 gaer_bc_tab(ns,k_bcins(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2273             :                enddo  ! nspint
    2274             :             enddo
    2275             :            else                ! dust
    2276           0 :             do k = 0, klev
    2277           0 :                do ns = 1,nspint   ! not weighted by aice
    2278           0 :                   tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
    2279           0 :                                    zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2280           0 :                   wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
    2281           0 :                                    zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2282           0 :                   gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
    2283           0 :                                    gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2284             :                enddo  ! nspint
    2285             :             enddo
    2286             :            endif      !(n=1)
    2287             :         enddo         ! n_zaero
    2288             :         endif         ! tr_zaero and dEdd_algae
    2289             : 
    2290             :      else  ! Bulk aerosol treatment
    2291   385847171 :         if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll
    2292           0 :         do n = 1,n_zaero          ! multiply by aice?
    2293           0 :             do k = 0, klev
    2294           0 :                do ns = 1,nspint   ! not weighted by aice
    2295           0 :                   tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
    2296           0 :                                    zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2297           0 :                   wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
    2298           0 :                                    zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2299           0 :                   gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
    2300           0 :                                    gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
    2301             :                enddo  ! nspint
    2302             :             enddo
    2303             :         enddo
    2304             :         endif  !tr_zaero  
    2305             : 
    2306             :      endif  ! modal_aero
    2307             : 
    2308             : !-----------------------------------------------------------------------
    2309             :  
    2310             :       ! begin spectral loop
    2311  1543388684 :       do ns = 1, nspint
    2312             :          
    2313             :          ! set optical properties of air/snow/pond overlying sea ice
    2314             :          ! air
    2315  1157541513 :          if( srftyp == 0 ) then
    2316   441719208 :             do k=0,nslyr 
    2317   294479472 :                tau(k) = c0
    2318   294479472 :                w0(k)  = c0
    2319   441719208 :                g(k)   = c0
    2320             :             enddo
    2321             :             ! snow
    2322  1010301777 :          else if( srftyp == 1 ) then
    2323             :             ! interpolate snow iops using input snow grain radius,
    2324             :             ! snow density and tabular data
    2325  2632245138 :             do k=0,nslyr
    2326             :                ! use top rsnw, rhosnw for snow ssl and rest of top layer
    2327  1754830092 :                ksnow = k - min(k-1,0)
    2328             :                ! find snow iops using input snow density and snow grain radius:
    2329  1754830092 :                if( frsnw(ksnow) < rsnw_tab(1) ) then
    2330           0 :                   Qs     = Qs_tab(ns,1)
    2331           0 :                   ws     = ws_tab(ns,1)
    2332           0 :                   gs     = gs_tab(ns,1)
    2333  1754830092 :                else if( frsnw(ksnow) >= rsnw_tab(nmbrad) ) then
    2334           0 :                   Qs     = Qs_tab(ns,nmbrad)
    2335           0 :                   ws     = ws_tab(ns,nmbrad)
    2336           0 :                   gs     = gs_tab(ns,nmbrad)
    2337             :                else
    2338             :                   ! linear interpolation in rsnw
    2339 56154562944 :                   do nr=2,nmbrad
    2340 58034769540 :                      if( rsnw_tab(nr-1) <= frsnw(ksnow) .and. &
    2341  5389866780 :                          frsnw(ksnow) < rsnw_tab(nr)) then
    2342   234518496 :                         delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / &
    2343  1872089340 :                                (rsnw_tab(nr) - rsnw_tab(nr-1))
    2344   117259248 :                         Qs   = Qs_tab(ns,nr-1)*(c1-delr) + &
    2345  1872089340 :                                Qs_tab(ns,nr)*delr
    2346   117259248 :                         ws   = ws_tab(ns,nr-1)*(c1-delr) + &
    2347  1872089340 :                                ws_tab(ns,nr)*delr
    2348   117259248 :                         gs   = gs_tab(ns,nr-1)*(c1-delr) + &
    2349  1872089340 :                                gs_tab(ns,nr)*delr
    2350             :                      endif
    2351             :                   enddo       ! nr
    2352             :                endif
    2353   117259248 :                ks = Qs*((rhosnw(ksnow)/rhoi)*3._dbl_kind / &
    2354  1754830092 :                        (4._dbl_kind*frsnw(ksnow)*1.0e-6_dbl_kind))
    2355             : 
    2356  1754830092 :                tau(k) = (ks + kabs_chl(ns,k))*dzk(k)
    2357  1754830092 :                w0(k)  = ks/(ks + kabs_chl(ns,k)) *ws
    2358  2632245138 :                g(k)   = gs
    2359             :             enddo       ! k
    2360             : 
    2361             : 
    2362             :             ! aerosol in snow
    2363   877415046 :              if (tr_zaero .and. dEdd_algae) then 
    2364           0 :                do k = 0,nslyr
    2365           0 :                   gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
    2366           0 :                   wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
    2367           0 :                   g(k)   = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
    2368           0 :                                   (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
    2369           0 :                   w0(k)  = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
    2370           0 :                                    (tau(k) + tzaer(ns,k))
    2371           0 :                   tau(k) = tau(k) + tzaer(ns,k)
    2372             :                enddo
    2373   877415046 :              elseif (tr_aero) then
    2374    41202300 :                   k = 0  ! snow SSL
    2375    41202300 :                taer = c0
    2376    41202300 :                waer = c0
    2377    41202300 :                gaer = c0
    2378             : 
    2379    82404600 :                do na=1,4*n_aero,4
    2380             : ! mgf++
    2381    82404600 :                if (modal_aero) then
    2382           0 :                   if (na == 1) then  
    2383             :                   !interstitial BC
    2384             :                      taer = taer + &
    2385           0 :                           aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))
    2386             :                      waer = waer + &
    2387           0 :                           aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
    2388           0 :                           waer_bc_tab(ns,k_bcexs(k))
    2389             :                      gaer = gaer + &
    2390           0 :                           aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
    2391           0 :                            waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
    2392           0 :                   elseif (na == 5)then  
    2393             :                   !within-ice BC
    2394             :                       taer = taer + &
    2395           0 :                            aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
    2396           0 :                            bcenh(ns,k_bcins(k),k_bcini(k))
    2397             :                       waer = waer + &
    2398           0 :                            aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
    2399           0 :                            waer_bc_tab(ns,k_bcins(k))
    2400             :                       gaer = gaer + &
    2401           0 :                            aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
    2402           0 :                            waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
    2403             :                    else
    2404             :                       ! other species (dust)
    2405             :                       taer = taer + &
    2406           0 :                            aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
    2407             :                       waer = waer + &
    2408           0 :                            aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
    2409           0 :                            waer_tab(ns,(1+(na-1)/4))
    2410             :                       gaer = gaer + &
    2411           0 :                            aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
    2412           0 :                            waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2413             :                    endif
    2414             :                else
    2415             :                   taer = taer + &
    2416    41202300 :                        aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
    2417             :                   waer = waer + &
    2418      690870 :                        aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
    2419    41547735 :                        waer_tab(ns,(1+(na-1)/4))
    2420             :                   gaer = gaer + &
    2421      690870 :                        aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
    2422    41547735 :                        waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2423             :                endif  !modal_aero
    2424             : !mgf--
    2425             :                enddo       ! na
    2426    41202300 :                gaer = gaer/(waer+puny)
    2427    41202300 :                waer = waer/(taer+puny)
    2428             : 
    2429    82404600 :                do k=1,nslyr
    2430    41202300 :                   taer = c0
    2431    41202300 :                   waer = c0
    2432    41202300 :                   gaer = c0
    2433    82404600 :                   do na=1,4*n_aero,4
    2434    82404600 :                   if (modal_aero) then
    2435             : !mgf++
    2436           0 :                     if (na==1) then
    2437             :                       ! interstitial BC
    2438             :                       taer = taer + &
    2439           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))
    2440             :                       waer = waer + &
    2441           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
    2442           0 :                            waer_bc_tab(ns,k_bcexs(k))
    2443             :                       gaer = gaer + &
    2444           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
    2445           0 :                            waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
    2446           0 :                     elseif (na==5) then
    2447             :                       ! within-ice BC
    2448             :                       taer = taer + &
    2449           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))*&
    2450           0 :                            bcenh(ns,k_bcins(k),k_bcini(k))
    2451             :                       waer = waer + &
    2452           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
    2453           0 :                            waer_bc_tab(ns,k_bcins(k))
    2454             :                       gaer = gaer + &
    2455           0 :                            (aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
    2456           0 :                            waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
    2457             :                       
    2458             :                     else
    2459             :                       ! other species (dust)
    2460             :                       taer = taer + &
    2461           0 :                            (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))
    2462             :                       waer = waer + &
    2463           0 :                            (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2464           0 :                            waer_tab(ns,(1+(na-1)/4))
    2465             :                       gaer = gaer + &
    2466           0 :                            (aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2467           0 :                            waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2468             :                     endif   !(na==1)
    2469             : 
    2470             :                   else
    2471             :                      taer = taer + &
    2472    41202300 :                             (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))
    2473             :                      waer = waer + &
    2474      690870 :                             (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2475    41893170 :                             waer_tab(ns,(1+(na-1)/4))
    2476             :                      gaer = gaer + &
    2477      690870 :                             (aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2478    41893170 :                             waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2479             :                   endif       ! modal_aero
    2480             : !mgf--
    2481             :                   enddo       ! na
    2482    41202300 :                   gaer = gaer/(waer+puny)
    2483    41202300 :                   waer = waer/(taer+puny)
    2484      345435 :                   g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
    2485    41202300 :                            (w0(k)*tau(k) + waer*taer)
    2486      345435 :                   w0(k)  = (w0(k)*tau(k) + waer*taer) / &
    2487    41202300 :                            (tau(k) + taer)
    2488    82404600 :                   tau(k) = tau(k) + taer
    2489             :                enddo       ! k
    2490             :             endif     ! tr_aero
    2491             : 
    2492             :             ! pond
    2493             :          else !if( srftyp == 2 ) then
    2494             :             ! pond water layers evenly spaced
    2495   132886731 :             dz = hp/(c1/rnslyr+c1)
    2496   398660193 :             do k=0,nslyr
    2497   265773462 :                tau(k) = kw(ns)*dz
    2498   265773462 :                w0(k)  = ww(ns)
    2499   398660193 :                g(k)   = gw(ns)
    2500             :                ! no aerosol in pond
    2501             :             enddo       ! k
    2502             :          endif        ! srftyp
    2503             :          
    2504             :          ! set optical properties of sea ice
    2505             :          
    2506             :          ! bare or snow-covered sea ice layers
    2507  1157541513 :          if( srftyp <= 1 ) then
    2508             :             ! ssl
    2509  1024654782 :             k = kii
    2510  1024654782 :             tau(k) = (ki_ssl(ns)+kabs_chl(ns,k))*dzk(k)
    2511  1024654782 :             w0(k)  = ki_ssl(ns)/(ki_ssl(ns) + kabs_chl(ns,k))*wi_ssl(ns)
    2512  1024654782 :             g(k)   = gi_ssl(ns)
    2513             :             ! dl
    2514  1024654782 :             k = kii + 1
    2515             :             ! scale dz for dl relative to 4 even-layer-thickness 1.5m case
    2516  1024654782 :             fs = p25/rnilyr
    2517  1024654782 :             tau(k) = (ki_dl(ns) + kabs_chl(ns,k)) *dzk(k)*fs
    2518  1024654782 :             w0(k)  = ki_dl(ns)/(ki_dl(ns) + kabs_chl(ns,k)) *wi_dl(ns)
    2519  1024654782 :             g(k)   = gi_dl(ns)
    2520             :             ! int above lowest layer
    2521  1024654782 :             if (kii+2 <= klev-1) then
    2522  6147928692 :                do k = kii+2, klev-1
    2523  5123273910 :                   tau(k) = (ki_int(ns) + kabs_chl(ns,k))*dzk(k)
    2524  5123273910 :                   w0(k)  = ki_int(ns)/(ki_int(ns) + kabs_chl(ns,k)) *wi_int(ns)
    2525  6147928692 :                   g(k)   = gi_int(ns)
    2526             :                enddo
    2527             :             endif
    2528             :             ! lowest layer
    2529  1024654782 :             k = klev
    2530             :             ! add algae to lowest sea ice layer, visible only:
    2531  1024654782 :             kabs = ki_int(ns)*(c1-wi_int(ns))
    2532  1024654782 :             if( ns == 1 ) then
    2533             :                ! total layer absorption optical depth fixed at value
    2534             :                ! of kalg*0.50m, independent of actual layer thickness
    2535   341551594 :                kabs = kabs + kabs_chl(ns,k) 
    2536             :             endif
    2537  1024654782 :             sig        = ki_int(ns)*wi_int(ns)
    2538  1024654782 :             tau(k) = (kabs+sig)*dzk(k)
    2539  1024654782 :             w0(k)  = (sig/(sig+kabs))
    2540  1024654782 :             g(k)   = gi_int(ns)
    2541             :             ! aerosol in sea ice
    2542  1024654782 :             if (tr_zaero .and. dEdd_algae) then
    2543           0 :                do k = kii, klev                  
    2544           0 :                   gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
    2545           0 :                   wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
    2546           0 :                   g(k)   = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
    2547           0 :                                   (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
    2548           0 :                   w0(k)  = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
    2549           0 :                                    (tau(k) + tzaer(ns,k))
    2550           0 :                   tau(k) = tau(k) + tzaer(ns,k)
    2551             :                enddo
    2552  1024654782 :             elseif (tr_aero) then 
    2553    47385600 :                k = kii   ! sea ice SSL
    2554    47385600 :                taer = c0
    2555    47385600 :                waer = c0
    2556    47385600 :                gaer = c0
    2557    94771200 :                do na=1,4*n_aero,4
    2558             : !mgf++
    2559    94771200 :                if (modal_aero) then
    2560           0 :                   if (na==1) then
    2561             :                   ! interstitial BC
    2562             :                      taer = taer + &
    2563           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))
    2564             :                      waer = waer + &
    2565           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
    2566           0 :                           waer_bc_tab(ns,k_bcexs(k))
    2567             :                      gaer = gaer + &
    2568           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
    2569           0 :                           waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
    2570           0 :                   elseif (na==5) then
    2571             :                   ! within-ice BC
    2572             :                      taer = taer + &
    2573           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
    2574           0 :                           bcenh(ns,k_bcins(k),k_bcini(k))
    2575             :                      waer = waer + &
    2576           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
    2577           0 :                           waer_bc_tab(ns,k_bcins(k))
    2578             :                      gaer = gaer + &
    2579           0 :                           aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
    2580           0 :                           waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))    
    2581             :                   else
    2582             :                   ! other species (dust)
    2583             :                      taer = taer + &
    2584           0 :                           aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
    2585             :                      waer = waer + &
    2586           0 :                           aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
    2587           0 :                           waer_tab(ns,(1+(na-1)/4))
    2588             :                      gaer = gaer + &
    2589           0 :                           aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
    2590           0 :                           waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2591             :                   endif
    2592             :                else      !bulk
    2593             :                   taer = taer + &
    2594    47385600 :                          aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
    2595             :                   waer = waer + &
    2596      708966 :                          aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
    2597    48094566 :                          waer_tab(ns,(1+(na-1)/4))
    2598             :                   gaer = gaer + &
    2599      708966 :                          aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
    2600    48094566 :                          waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2601             :                 endif     ! modal_aero
    2602             : !mgf--
    2603             :                enddo      ! na
    2604             : 
    2605    47385600 :                gaer = gaer/(waer+puny)
    2606    47385600 :                waer = waer/(taer+puny)
    2607      354483 :                g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
    2608    47385600 :                     (w0(k)*tau(k) + waer*taer)
    2609      354483 :                w0(k)  = (w0(k)*tau(k) + waer*taer) / &
    2610    47385600 :                     (tau(k) + taer)
    2611    47385600 :                tau(k) = tau(k) + taer
    2612   379084800 :                do k = kii+1, klev
    2613   331699200 :                   taer = c0
    2614   331699200 :                   waer = c0
    2615   331699200 :                   gaer = c0
    2616   663398400 :                   do na=1,4*n_aero,4
    2617             : !mgf++
    2618   663398400 :                   if (modal_aero) then
    2619           0 :                      if (na==1) then
    2620             :                         ! interstitial BC
    2621             :                         taer = taer + &
    2622           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))
    2623             :                         waer = waer + &
    2624           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
    2625           0 :                              waer_bc_tab(ns,k_bcexs(k))
    2626             :                         gaer = gaer + &
    2627           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
    2628           0 :                              waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
    2629           0 :                      elseif (na==5) then
    2630             :                         ! within-ice BC
    2631             :                         taer = taer + &
    2632           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
    2633           0 :                              bcenh(ns,k_bcins(k),k_bcini(k))
    2634             :                         waer = waer + &
    2635           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
    2636           0 :                              waer_bc_tab(ns,k_bcins(k))
    2637             :                         gaer = gaer + &
    2638           0 :                              (aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
    2639           0 :                              waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
    2640             :                         
    2641             :                      else
    2642             :                         ! other species (dust)
    2643             :                         taer = taer + &
    2644           0 :                              (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))
    2645             :                         waer = waer + &
    2646           0 :                              (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2647           0 :                              waer_tab(ns,(1+(na-1)/4))
    2648             :                         gaer = gaer + &
    2649           0 :                              (aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2650           0 :                              waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2651             :                      endif
    2652             :                   else       !bulk
    2653             : 
    2654             :                      taer = taer + &
    2655   331699200 :                           (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))
    2656             :                      waer = waer + &
    2657     4962762 :                           (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2658   336661962 :                           waer_tab(ns,(1+(na-1)/4))
    2659             :                      gaer = gaer + &
    2660     4962762 :                           (aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
    2661   336661962 :                           waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
    2662             :                   endif       ! modal_aero
    2663             : !mgf--
    2664             :                   enddo       ! na
    2665   331699200 :                   gaer = gaer/(waer+puny)
    2666   331699200 :                   waer = waer/(taer+puny)
    2667     2481381 :                   g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
    2668   331699200 :                            (w0(k)*tau(k) + waer*taer)
    2669     2481381 :                   w0(k)  = (w0(k)*tau(k) + waer*taer) / &
    2670   331699200 :                            (tau(k) + taer)
    2671   379084800 :                   tau(k) = tau(k) + taer
    2672             :                enddo ! k
    2673             :             endif ! tr_aero
    2674             : 
    2675             :             ! sea ice layers under ponds
    2676             :          else !if( srftyp == 2 ) then
    2677   132886731 :             k = kii
    2678   132886731 :             tau(k) = ki_p_ssl(ns)*dzk(k)
    2679   132886731 :             w0(k)  = wi_p_ssl(ns)
    2680   132886731 :             g(k)   = gi_p_ssl(ns)
    2681   132886731 :             k = kii + 1
    2682   132886731 :             tau(k) = ki_p_int(ns)*dzk(k)
    2683   132886731 :             w0(k)  = wi_p_int(ns)
    2684   132886731 :             g(k)   = gi_p_int(ns)
    2685   132886731 :             if (kii+2 <= klev) then
    2686   930207117 :                do k = kii+2, klev
    2687   797320386 :                   tau(k) = ki_p_int(ns)*dzk(k)
    2688   797320386 :                   w0(k)  = wi_p_int(ns)
    2689   930207117 :                   g(k)   = gi_p_int(ns)
    2690             :                enddo       ! k
    2691             :             endif
    2692             :             ! adjust pond iops if pond depth within specified range
    2693   132886731 :             if( hpmin <= hp .and. hp <= hp0 ) then
    2694   128569785 :                k = kii
    2695   128569785 :                sig_i      = ki_ssl(ns)*wi_ssl(ns)
    2696   128569785 :                sig_p      = ki_p_ssl(ns)*wi_p_ssl(ns)
    2697   128569785 :                sig        = sig_i + (sig_p-sig_i)*(hp/hp0)
    2698   128569785 :                kext       = sig + ki_p_ssl(ns)*(c1-wi_p_ssl(ns))
    2699   128569785 :                tau(k) = kext*dzk(k)
    2700   128569785 :                w0(k) = sig/kext
    2701   128569785 :                g(k)  = gi_p_int(ns)
    2702   128569785 :                k = kii + 1
    2703             :                ! scale dz for dl relative to 4 even-layer-thickness 1.5m case
    2704   128569785 :                fs = p25/rnilyr
    2705   128569785 :                sig_i      = ki_dl(ns)*wi_dl(ns)*fs
    2706   128569785 :                sig_p      = ki_p_int(ns)*wi_p_int(ns)
    2707   128569785 :                sig        = sig_i + (sig_p-sig_i)*(hp/hp0)
    2708   128569785 :                kext       = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
    2709   128569785 :                tau(k) = kext*dzk(k)
    2710   128569785 :                w0(k) = sig/kext
    2711   128569785 :                g(k)  = gi_p_int(ns)
    2712   128569785 :                if (kii+2 <= klev) then
    2713   899988495 :                   do k = kii+2, klev
    2714   771418710 :                      sig_i      = ki_int(ns)*wi_int(ns)
    2715   771418710 :                      sig_p      = ki_p_int(ns)*wi_p_int(ns)
    2716   771418710 :                      sig        = sig_i + (sig_p-sig_i)*(hp/hp0)
    2717   771418710 :                      kext       = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
    2718   771418710 :                      tau(k) = kext*dzk(k)
    2719   771418710 :                      w0(k) = sig/kext
    2720   899988495 :                      g(k)  = gi_p_int(ns)
    2721             :                   enddo       ! k
    2722             :                endif
    2723             :             endif        ! small pond depth transition to bare sea ice
    2724             :          endif         ! srftyp  
    2725             :          
    2726             :          ! set reflectivities for ocean underlying sea ice
    2727  1157541513 :          rns = real(ns-1, kind=dbl_kind)
    2728  1157541513 :          albodr = cp01 * (c1 - min(rns, c1))
    2729  1157541513 :          albodf = cp01 * (c1 - min(rns, c1))
    2730             :          
    2731             :          ! layer input properties now completely specified: tau, w0, g,
    2732             :          ! albodr, albodf; now compute the Delta-Eddington solution 
    2733             :          ! reflectivities and transmissivities for each layer; then,
    2734             :          ! combine the layers going downwards accounting for multiple
    2735             :          ! scattering between layers, and finally start from the 
    2736             :          ! underlying ocean and combine successive layers upwards to
    2737             :          ! the surface; see comments in solution_dEdd for more details.
    2738             :          
    2739             :          call solution_dEdd &
    2740             :                (coszen,     srftyp,     klev,       klevp,      nslyr,     &
    2741             :                 tau,        w0,         g,          albodr,     albodf,    &
    2742             :                 trndir,     trntdr,     trndif,     rupdir,     rupdif,    &
    2743  1157541513 :                 rdndif)   
    2744  1157541513 :          if (icepack_warnings_aborted(subname)) return
    2745             : 
    2746             :          ! the interface reflectivities and transmissivities required
    2747             :          ! to evaluate interface fluxes are returned from solution_dEdd;
    2748             :          ! now compute up and down fluxes for each interface, using the 
    2749             :          ! combined layer properties at each interface:
    2750             :          !
    2751             :          !              layers       interface
    2752             :          !
    2753             :          !       ---------------------  k
    2754             :          !                 k
    2755             :          !       --------------------- 
    2756             :          
    2757 13890498156 :          do k = 0, klevp 
    2758             :             ! interface scattering
    2759 12732956643 :             refk          = c1/(c1 - rdndif(k)*rupdif(k))
    2760             :             ! dir tran ref from below times interface scattering, plus diff
    2761             :             ! tran and ref from below times interface scattering
    2762             :             ! fdirup(k) = (trndir(k)*rupdir(k) + &
    2763             :             !                 (trntdr(k)-trndir(k))  &
    2764             :             !                 *rupdif(k))*refk
    2765             :             ! dir tran plus total diff trans times interface scattering plus
    2766             :             ! dir tran with up dir ref and down dif ref times interface scattering 
    2767             :             ! fdirdn(k) = trndir(k) + (trntdr(k) &
    2768             :             !               - trndir(k) + trndir(k)  &
    2769             :             !               *rupdir(k)*rdndif(k))*refk
    2770             :             ! diffuse tran ref from below times interface scattering
    2771             :             ! fdifup(k) = trndif(k)*rupdif(k)*refk
    2772             :             ! diffuse tran times interface scattering
    2773             :             ! fdifdn(k) = trndif(k)*refk
    2774             : 
    2775             :             ! dfdir = fdirdn - fdirup
    2776   801155355 :             dfdir(k) = trndir(k) &
    2777   801155355 :                         + (trntdr(k)-trndir(k)) * (c1 - rupdif(k)) * refk &
    2778 12732956643 :                         -  trndir(k)*rupdir(k)  * (c1 - rdndif(k)) * refk
    2779 12732956643 :             if (dfdir(k) < puny) dfdir(k) = c0 !echmod necessary?
    2780             :             ! dfdif = fdifdn - fdifup
    2781 12732956643 :             dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk
    2782 13890498156 :             if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary?
    2783             :          enddo       ! k 
    2784             :          
    2785             :          ! calculate final surface albedos and fluxes-
    2786             :          ! all absorbed flux above ksrf is included in surface absorption
    2787             :          
    2788  1543388684 :          if( ns == 1) then      ! visible
    2789             :             
    2790   385847171 :             swdr = swvdr
    2791   385847171 :             swdf = swvdf
    2792   385847171 :             avdr  = rupdir(0)
    2793   385847171 :             avdf  = rupdif(0)
    2794             :             
    2795   385847171 :             tmp_0  = dfdir(0    )*swdr + dfdif(0    )*swdf
    2796   385847171 :             tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
    2797   385847171 :             tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
    2798             : 
    2799             :             ! for layer biology: save visible only
    2800  3472624539 :             do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering
    2801  3472624539 :                fthrul(k-nslyr-1) = dfdir(k)*swdr + dfdif(k)*swdf
    2802             :             enddo
    2803             :             
    2804   385847171 :             fsfc  = fsfc  + tmp_0  - tmp_ks
    2805   385847171 :             fint  = fint  + tmp_ks - tmp_kl
    2806   385847171 :             fthru = fthru + tmp_kl
    2807             :             
    2808             :             ! if snow covered ice, set snow internal absorption; else, Sabs=0
    2809   385847171 :             if( srftyp == 1 ) then
    2810   292471682 :                ki = 0
    2811   584943364 :                do k=1,nslyr
    2812             :                   ! skip snow SSL, since SSL absorption included in the surface
    2813             :                   ! absorption fsfc above
    2814   292471682 :                   km  = k
    2815   292471682 :                   kp  = km + 1
    2816   292471682 :                   ki  = ki + 1
    2817    19543208 :                   Sabs(ki) = Sabs(ki) &
    2818    19543208 :                            +  dfdir(km)*swdr + dfdif(km)*swdf &
    2819   584943364 :                            - (dfdir(kp)*swdr + dfdif(kp)*swdf)
    2820             :                enddo       ! k
    2821             :             endif
    2822             :             
    2823             :             ! complex indexing to insure proper absorptions for sea ice
    2824   385847171 :             ki = 0
    2825  3086777368 :             do k=nslyr+2,nslyr+1+nilyr
    2826             :                ! for bare ice, DL absorption for sea ice layer 1
    2827  2700930197 :                km = k  
    2828  2700930197 :                kp = km + 1
    2829             :                ! modify for top sea ice layer for snow over sea ice
    2830  2700930197 :                if( srftyp == 1 ) then
    2831             :                   ! must add SSL and DL absorption for sea ice layer 1
    2832  2047301774 :                   if( k == nslyr+2 ) then
    2833   292471682 :                      km = k  - 1
    2834   292471682 :                      kp = km + 2
    2835             :                   endif
    2836             :                endif
    2837  2700930197 :                ki = ki + 1
    2838   169942045 :                Iabs(ki) = Iabs(ki) &
    2839   169942045 :                         +  dfdir(km)*swdr + dfdif(km)*swdf &
    2840  3086777368 :                         - (dfdir(kp)*swdr + dfdif(kp)*swdf)
    2841             :             enddo       ! k
    2842             :             
    2843             :          else !if(ns > 1) then  ! near IR
    2844             :             
    2845   771694342 :             swdr = swidr
    2846   771694342 :             swdf = swidf
    2847             : 
    2848             :             ! let fr1 = alb_1*swd*wght1 and fr2 = alb_2*swd*wght2 be the ns=2,3
    2849             :             ! reflected fluxes respectively, where alb_1, alb_2 are the band
    2850             :             ! albedos, swd = nir incident shortwave flux, and wght1, wght2 are
    2851             :             ! the 2,3 band weights. thus, the total reflected flux is:
    2852             :             ! fr = fr1 + fr2 = alb_1*swd*wght1 + alb_2*swd*wght2  hence, the
    2853             :             ! 2,3 nir band albedo is alb = fr/swd = alb_1*wght1 + alb_2*wght2
    2854             : 
    2855   771694342 :             aidr   = aidr + rupdir(0)*wghtns(ns)
    2856   771694342 :             aidf   = aidf + rupdif(0)*wghtns(ns)
    2857             : 
    2858   771694342 :             tmp_0  = dfdir(0    )*swdr + dfdif(0    )*swdf
    2859   771694342 :             tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
    2860   771694342 :             tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
    2861             : 
    2862   771694342 :             tmp_0  = tmp_0  * wghtns(ns)
    2863   771694342 :             tmp_ks = tmp_ks * wghtns(ns)
    2864   771694342 :             tmp_kl = tmp_kl * wghtns(ns)
    2865             : 
    2866   771694342 :             fsfc  = fsfc  + tmp_0  - tmp_ks
    2867   771694342 :             fint  = fint  + tmp_ks - tmp_kl
    2868   771694342 :             fthru = fthru + tmp_kl
    2869             : 
    2870             :             ! if snow covered ice, set snow internal absorption; else, Sabs=0
    2871   771694342 :             if( srftyp == 1 ) then
    2872   584943364 :                ki = 0
    2873  1169886728 :                do k=1,nslyr
    2874             :                   ! skip snow SSL, since SSL absorption included in the surface
    2875             :                   ! absorption fsfc above
    2876   584943364 :                   km  = k
    2877   584943364 :                   kp  = km + 1
    2878   584943364 :                   ki  = ki + 1
    2879    39086416 :                   Sabs(ki) = Sabs(ki) &
    2880    39086416 :                            + (dfdir(km)*swdr + dfdif(km)*swdf   &
    2881    39086416 :                            - (dfdir(kp)*swdr + dfdif(kp)*swdf)) & 
    2882  1169886728 :                            * wghtns(ns)
    2883             :                enddo       ! k
    2884             :             endif
    2885             :             
    2886             :             ! complex indexing to insure proper absorptions for sea ice
    2887   771694342 :             ki = 0
    2888  6173554736 :             do k=nslyr+2,nslyr+1+nilyr
    2889             :                ! for bare ice, DL absorption for sea ice layer 1
    2890  5401860394 :                km = k  
    2891  5401860394 :                kp = km + 1
    2892             :                ! modify for top sea ice layer for snow over sea ice
    2893  5401860394 :                if( srftyp == 1 ) then
    2894             :                   ! must add SSL and DL absorption for sea ice layer 1
    2895  4094603548 :                   if( k == nslyr+2 ) then
    2896   584943364 :                      km = k  - 1
    2897   584943364 :                      kp = km + 2
    2898             :                   endif
    2899             :                endif
    2900  5401860394 :                ki = ki + 1
    2901   339884090 :                Iabs(ki) = Iabs(ki) &
    2902   339884090 :                         + (dfdir(km)*swdr + dfdif(km)*swdf &
    2903   339884090 :                         - (dfdir(kp)*swdr + dfdif(kp)*swdf)) &
    2904  6173554736 :                         * wghtns(ns)
    2905             :             enddo       ! k
    2906             :             
    2907             :          endif        ! ns = 1, ns > 1
    2908             :          
    2909             :       enddo         ! end spectral loop  ns
    2910             : 
    2911             :       ! accumulate fluxes over bare sea ice
    2912   385847171 :       alvdr   = avdr
    2913   385847171 :       alvdf   = avdf
    2914   385847171 :       alidr   = aidr
    2915   385847171 :       alidf   = aidf
    2916   385847171 :       fswsfc  = fswsfc  + fsfc *fi
    2917   385847171 :       fswint  = fswint  + fint *fi
    2918   385847171 :       fswthru = fswthru + fthru*fi
    2919             : 
    2920   771694342 :       do k = 1, nslyr
    2921   771694342 :          Sswabs(k) = Sswabs(k) + Sabs(k)*fi
    2922             :       enddo                     ! k
    2923             :       
    2924  3086777368 :       do k = 1, nilyr
    2925  2700930197 :          Iswabs(k) = Iswabs(k) + Iabs(k)*fi
    2926             :          
    2927             :          ! bgc layer 
    2928  2700930197 :          fswpenl(k) = fswpenl(k) + fthrul(k)* fi
    2929  3086777368 :          if (k == nilyr) then
    2930   385847171 :             fswpenl(k+1) = fswpenl(k+1) + fthrul(k+1)*fi
    2931             :          endif
    2932             :       enddo                     ! k
    2933             :       
    2934             :       !----------------------------------------------------------------
    2935             :       ! if ice has zero heat capacity, no SW can be absorbed 
    2936             :       ! in the ice/snow interior, so add to surface absorption.
    2937             :       ! Note: nilyr = nslyr = 1 for this case
    2938             :       !----------------------------------------------------------------
    2939             : 
    2940   385847171 :       if (.not. heat_capacity) then
    2941             :          
    2942             :          ! SW absorbed at snow/ice surface
    2943           0 :          fswsfc = fswsfc + Iswabs(1) + Sswabs(1)
    2944             :          
    2945             :          ! SW absorbed in ice interior
    2946           0 :          fswint   = c0
    2947           0 :          Iswabs(1) = c0
    2948           0 :          Sswabs(1) = c0
    2949             : 
    2950             :       endif                       ! heat_capacity
    2951             : 
    2952             :       end subroutine compute_dEdd
    2953             : 
    2954             : !=======================================================================
    2955             : !
    2956             : ! Given input vertical profiles of optical properties, evaluate the 
    2957             : ! monochromatic Delta-Eddington solution.
    2958             : !
    2959             : ! author:  Bruce P. Briegleb, NCAR
    2960             : !   2013:  E Hunke merged with NCAR version
    2961  1157541513 :       subroutine solution_dEdd                                 &
    2962             :             (coszen,     srftyp,    klev,      klevp,  nslyr,  &
    2963  1157541513 :              tau,        w0,        g,         albodr, albodf, &
    2964  1157541513 :              trndir,     trntdr,    trndif,    rupdir, rupdif, &
    2965  1157541513 :              rdndif)
    2966             : 
    2967             :       real (kind=dbl_kind), intent(in) :: &
    2968             :          coszen      ! cosine solar zenith angle
    2969             : 
    2970             :       integer (kind=int_kind), intent(in) :: &
    2971             :          srftyp   , & ! surface type over ice: (0=air, 1=snow, 2=pond)
    2972             :          klev     , & ! number of radiation layers - 1
    2973             :          klevp    , & ! number of radiation interfaces - 1
    2974             :                       ! (0 layer is included also)
    2975             :          nslyr        ! number of snow layers
    2976             :  
    2977             :       real (kind=dbl_kind), dimension(0:klev), intent(in) :: &
    2978             :          tau     , & ! layer extinction optical depth
    2979             :          w0      , & ! layer single scattering albedo
    2980             :          g           ! layer asymmetry parameter
    2981             :  
    2982             :       real (kind=dbl_kind), intent(in) :: &
    2983             :          albodr  , & ! ocean albedo to direct rad
    2984             :          albodf      ! ocean albedo to diffuse rad
    2985             :  
    2986             :       ! following arrays are defined at model interfaces; 0 is the top of the
    2987             :       ! layer above the sea ice; klevp is the sea ice/ocean interface.
    2988             :       real (kind=dbl_kind), dimension (0:klevp), intent(out) :: &
    2989             :          trndir  , & ! solar beam down transmission from top
    2990             :          trntdr  , & ! total transmission to direct beam for layers above
    2991             :          trndif  , & ! diffuse transmission to diffuse beam for layers above
    2992             :          rupdir  , & ! reflectivity to direct radiation for layers below
    2993             :          rupdif  , & ! reflectivity to diffuse radiation for layers below
    2994             :          rdndif      ! reflectivity to diffuse radiation for layers above
    2995             : 
    2996             : !-----------------------------------------------------------------------
    2997             : !
    2998             : ! Delta-Eddington solution for snow/air/pond over sea ice
    2999             : !
    3000             : ! Generic solution for a snow/air/pond input column of klev+1 layers,
    3001             : ! with srftyp determining at what interface fresnel refraction occurs.
    3002             : !
    3003             : ! Computes layer reflectivities and transmissivities, from the top down
    3004             : ! to the lowest interface using the Delta-Eddington solutions for each
    3005             : ! layer; combines layers from top down to lowest interface, and from the
    3006             : ! lowest interface (underlying ocean) up to the top of the column.
    3007             : !
    3008             : ! Note that layer diffuse reflectivity and transmissivity are computed
    3009             : ! by integrating the direct over several gaussian angles. This is
    3010             : ! because the diffuse reflectivity expression sometimes is negative,
    3011             : ! but the direct reflectivity is always well-behaved. We assume isotropic
    3012             : ! radiation in the upward and downward hemispheres for this integration.
    3013             : !
    3014             : ! Assumes monochromatic (spectrally uniform) properties across a band
    3015             : ! for the input optical parameters.
    3016             : !
    3017             : ! If total transmission of the direct beam to the interface above a particular 
    3018             : ! layer is less than trmin, then no further Delta-Eddington solutions are
    3019             : ! evaluated for layers below.
    3020             : !
    3021             : ! The following describes how refraction is handled in the calculation.
    3022             : !
    3023             : ! First, we assume that radiation is refracted when entering either
    3024             : ! sea ice at the base of the surface scattering layer, or water (i.e. melt
    3025             : ! pond); we assume that radiation does not refract when entering snow, nor 
    3026             : ! upon entering sea ice from a melt pond, nor upon entering the underlying 
    3027             : ! ocean from sea ice.
    3028             : !
    3029             : ! To handle refraction, we define a "fresnel" layer, which physically
    3030             : ! is of neglible thickness and is non-absorbing, which can be combined to 
    3031             : ! any sea ice layer or top of melt pond. The fresnel layer accounts for 
    3032             : ! refraction of direct beam and associated reflection and transmission for
    3033             : ! solar radiation. A fresnel layer is combined with the top of a melt pond 
    3034             : ! or to the surface scattering layer of sea ice if no melt pond lies over it. 
    3035             : !
    3036             : ! Some caution must be exercised for the fresnel layer, because any layer
    3037             : ! to which it is combined is no longer a homogeneous layer, as are all other
    3038             : ! individual layers. For all other layers for example, the direct and diffuse
    3039             : ! reflectivities/transmissivities (R/T) are the same for radiation above or
    3040             : ! below the layer. This is the meaning of homogeneous! But for the fresnel
    3041             : ! layer this is not so. Thus, the R/T for this layer must be distinguished
    3042             : ! for radiation above from that from radiation below. For generality, we
    3043             : ! treat all layers to be combined as inhomogeneous.
    3044             : !
    3045             : !-----------------------------------------------------------------------
    3046             : 
    3047             :       ! local variables
    3048             : 
    3049             :       integer (kind=int_kind) :: &
    3050             :          kfrsnl      ! radiation interface index for fresnel layer
    3051             :  
    3052             :       ! following variables are defined for each layer; 0 refers to the top
    3053             :       ! layer. In general we must distinguish directions above and below in 
    3054             :       ! the diffuse reflectivity and transmissivity, as layers are not assumed
    3055             :       ! to be homogeneous (apart from the single layer Delta-Edd solutions); 
    3056             :       ! the direct is always from above.
    3057             :       real (kind=dbl_kind), dimension (0:klev) :: &
    3058  2970573771 :          rdir    , & ! layer reflectivity to direct radiation
    3059  2970573771 :          rdif_a  , & ! layer reflectivity to diffuse radiation from above
    3060  2970573771 :          rdif_b  , & ! layer reflectivity to diffuse radiation from below
    3061  2970573771 :          tdir    , & ! layer transmission to direct radiation (solar beam + diffuse)
    3062  2970573771 :          tdif_a  , & ! layer transmission to diffuse radiation from above
    3063  2970573771 :          tdif_b  , & ! layer transmission to diffuse radiation from below
    3064  2970573771 :          trnlay      ! solar beam transm for layer (direct beam only)
    3065             : 
    3066             :       integer (kind=int_kind) :: & 
    3067             :          k           ! level index
    3068             :  
    3069             :       real (kind=dbl_kind), parameter :: &
    3070             :          trmin = 0.001_dbl_kind   ! minimum total transmission allowed
    3071             :       ! total transmission is that due to the direct beam; i.e. it includes
    3072             :       ! both the directly transmitted solar beam and the diffuse downwards
    3073             :       ! transmitted radiation resulting from scattering out of the direct beam 
    3074             :       real (kind=dbl_kind) :: &
    3075    72832305 :          tautot   , & ! layer optical depth
    3076    72832305 :          wtot     , & ! layer single scattering albedo
    3077    72832305 :          gtot     , & ! layer asymmetry parameter
    3078    72832305 :          ftot     , & ! layer forward scattering fraction
    3079    72832305 :          ts       , & ! layer scaled extinction optical depth
    3080    72832305 :          ws       , & ! layer scaled single scattering albedo
    3081    72832305 :          gs       , & ! layer scaled asymmetry parameter
    3082    72832305 :          rintfc   , & ! reflection (multiple) at an interface
    3083    72832305 :          refkp1   , & ! interface multiple scattering for k+1
    3084    72832305 :          refkm1   , & ! interface multiple scattering for k-1
    3085    72832305 :          tdrrdir  , & ! direct tran times layer direct ref 
    3086    72832305 :          tdndif       ! total down diffuse = tot tran - direct tran
    3087             :  
    3088             :       ! perpendicular and parallel relative to plane of incidence and scattering
    3089             :       real (kind=dbl_kind) :: &
    3090    72832305 :          R1       , & ! perpendicular polarization reflection amplitude
    3091    72832305 :          R2       , & ! parallel polarization reflection amplitude
    3092    72832305 :          T1       , & ! perpendicular polarization transmission amplitude
    3093    72832305 :          T2       , & ! parallel polarization transmission amplitude
    3094    72832305 :          Rf_dir_a , & ! fresnel reflection to direct radiation
    3095    72832305 :          Tf_dir_a , & ! fresnel transmission to direct radiation
    3096    72832305 :          Rf_dif_a , & ! fresnel reflection to diff radiation from above
    3097    72832305 :          Rf_dif_b , & ! fresnel reflection to diff radiation from below
    3098    72832305 :          Tf_dif_a , & ! fresnel transmission to diff radiation from above
    3099    72832305 :          Tf_dif_b     ! fresnel transmission to diff radiation from below
    3100             :  
    3101             :       ! refractive index for sea ice, water; pre-computed, band-independent,
    3102             :       ! diffuse fresnel reflectivities
    3103             :       real (kind=dbl_kind), parameter :: & 
    3104             :          refindx = 1.310_dbl_kind  , & ! refractive index of sea ice (water also)
    3105             :          cp063   = 0.063_dbl_kind  , & ! diffuse fresnel reflectivity from above
    3106             :          cp455   = 0.455_dbl_kind      ! diffuse fresnel reflectivity from below
    3107             :  
    3108             :       real (kind=dbl_kind) :: &
    3109    72832305 :          mu0      , & ! cosine solar zenith angle incident
    3110    72832305 :          mu0nij       ! cosine solar zenith angle in medium below fresnel level
    3111             :  
    3112             :       real (kind=dbl_kind) :: &
    3113    72832305 :          mu0n         ! cosine solar zenith angle in medium
    3114             :  
    3115             :       real (kind=dbl_kind) :: &
    3116    72832305 :          alp      , & ! temporary for alpha
    3117    72832305 :          gam      , & ! temporary for agamm
    3118    72832305 :          lm       , & ! temporary for el
    3119    72832305 :          mu       , & ! temporary for gauspt
    3120    72832305 :          ne       , & ! temporary for n
    3121    72832305 :          ue       , & ! temporary for u
    3122    72832305 :          extins   , & ! extinction
    3123    72832305 :          amg      , & ! alp - gam
    3124    72832305 :          apg          ! alp + gam
    3125             : 
    3126             :       integer (kind=int_kind), parameter :: &
    3127             :          ngmax = 8    ! number of gaussian angles in hemisphere
    3128             : 
    3129             :       real (kind=dbl_kind), dimension (ngmax), parameter :: &
    3130             :          gauspt     & ! gaussian angles (radians)
    3131             :             = (/ .9894009_dbl_kind,  .9445750_dbl_kind, &
    3132             :                  .8656312_dbl_kind,  .7554044_dbl_kind, &
    3133             :                  .6178762_dbl_kind,  .4580168_dbl_kind, &
    3134             :                  .2816036_dbl_kind,  .0950125_dbl_kind/), &
    3135             :          gauswt     & ! gaussian weights
    3136             :             = (/ .0271525_dbl_kind,  .0622535_dbl_kind, &
    3137             :                  .0951585_dbl_kind,  .1246290_dbl_kind, &
    3138             :                  .1495960_dbl_kind,  .1691565_dbl_kind, &
    3139             :                  .1826034_dbl_kind,  .1894506_dbl_kind/)
    3140             : 
    3141             :       integer (kind=int_kind) :: &
    3142             :          ng           ! gaussian integration index
    3143             : 
    3144             :       real (kind=dbl_kind) :: &
    3145    72832305 :          gwt      , & ! gaussian weight
    3146    72832305 :          swt      , & ! sum of weights
    3147    72832305 :          trn      , & ! layer transmission
    3148    72832305 :          rdr      , & ! rdir for gaussian integration
    3149    72832305 :          tdr      , & ! tdir for gaussian integration
    3150    72832305 :          smr      , & ! accumulator for rdif gaussian integration
    3151    72832305 :          smt          ! accumulator for tdif gaussian integration
    3152             : 
    3153             :       real (kind=dbl_kind) :: &
    3154    72832305 :          exp_min                    ! minimum exponential value
    3155             : 
    3156             :       character(len=*),parameter :: subname='(solution_dEdd)'
    3157             : 
    3158             : !-----------------------------------------------------------------------
    3159             : 
    3160 13890498156 :       do k = 0, klevp 
    3161 12732956643 :          trndir(k) = c0
    3162 12732956643 :          trntdr(k) = c0
    3163 12732956643 :          trndif(k) = c0
    3164 12732956643 :          rupdir(k) = c0
    3165 12732956643 :          rupdif(k) = c0
    3166 13890498156 :          rdndif(k) = c0
    3167             :       enddo
    3168             : 
    3169             :       ! initialize top interface of top layer 
    3170  1157541513 :       trndir(0) =   c1
    3171  1157541513 :       trntdr(0) =   c1
    3172  1157541513 :       trndif(0) =   c1
    3173  1157541513 :       rdndif(0) =   c0
    3174             : 
    3175             :       ! mu0 is cosine solar zenith angle above the fresnel level; make 
    3176             :       ! sure mu0 is large enough for stable and meaningful radiation
    3177             :       ! solution: .01 is like sun just touching horizon with its lower edge
    3178  1157541513 :       mu0  = max(coszen,p01)
    3179             : 
    3180             :       ! mu0n is cosine solar zenith angle used to compute the layer
    3181             :       ! Delta-Eddington solution; it is initially computed to be the
    3182             :       ! value below the fresnel level, i.e. the cosine solar zenith 
    3183             :       ! angle below the fresnel level for the refracted solar beam:
    3184  1157541513 :       mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
    3185             : 
    3186             :       ! compute level of fresnel refraction
    3187             :       ! if ponded sea ice, fresnel level is the top of the pond.
    3188  1157541513 :       kfrsnl = 0
    3189             :       ! if snow over sea ice or bare sea ice, fresnel level is
    3190             :       ! at base of sea ice SSL (and top of the sea ice DL); the
    3191             :       ! snow SSL counts for one, then the number of snow layers,
    3192             :       ! then the sea ice SSL which also counts for one:
    3193  1157541513 :       if( srftyp < 2 ) kfrsnl = nslyr + 2 
    3194             : 
    3195             :       ! proceed down one layer at a time; if the total transmission to
    3196             :       ! the interface just above a given layer is less than trmin, then no
    3197             :       ! Delta-Eddington computation for that layer is done.
    3198             : 
    3199             :       ! begin main level loop
    3200 12732956643 :       do k = 0, klev
    3201             : 
    3202             :          ! initialize all layer apparent optical properties to 0
    3203 11575415130 :          rdir  (k) = c0
    3204 11575415130 :          rdif_a(k) = c0
    3205 11575415130 :          rdif_b(k) = c0
    3206 11575415130 :          tdir  (k) = c0
    3207 11575415130 :          tdif_a(k) = c0
    3208 11575415130 :          tdif_b(k) = c0
    3209 11575415130 :          trnlay(k) = c0
    3210             : 
    3211             :          ! compute next layer Delta-eddington solution only if total transmission
    3212             :          ! of radiation to the interface just above the layer exceeds trmin.
    3213             :          
    3214 11575415130 :          if (trntdr(k) > trmin ) then
    3215             : 
    3216             :             ! calculation over layers with penetrating radiation
    3217             : 
    3218  5535215076 :             tautot  = tau(k)
    3219  5535215076 :             wtot    = w0(k)
    3220  5535215076 :             gtot    = g(k)
    3221  5535215076 :             ftot    = gtot*gtot
    3222             : 
    3223  5535215076 :             ts   = taus(wtot,ftot,tautot)
    3224  5535215076 :             ws   = omgs(wtot,ftot)
    3225  5535215076 :             gs   = asys(gtot,ftot)
    3226  5535215076 :             lm   = el(ws,gs)
    3227  5535215076 :             ue   = u(ws,gs,lm)
    3228             : 
    3229  5535215076 :             mu0n = mu0nij
    3230             :             ! if level k is above fresnel level and the cell is non-pond, use the
    3231             :             ! non-refracted beam instead
    3232  5535215076 :             if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
    3233             : 
    3234  5535215076 :             exp_min = min(exp_argmax,lm*ts)
    3235  5535215076 :             extins = exp(-exp_min)
    3236  5535215076 :             ne = n(ue,extins)
    3237             : 
    3238             :             ! first calculation of rdif, tdif using Delta-Eddington formulas
    3239             : !            rdif_a(k) = (ue+c1)*(ue-c1)*(c1/extins - extins)/ne
    3240  5535215076 :             rdif_a(k) = (ue**2-c1)*(c1/extins - extins)/ne
    3241  5535215076 :             tdif_a(k) = c4*ue/ne
    3242             : 
    3243             :             ! evaluate rdir,tdir for direct beam
    3244  5535215076 :             exp_min = min(exp_argmax,ts/mu0n)
    3245  5535215076 :             trnlay(k) = exp(-exp_min)
    3246  5535215076 :             alp = alpha(ws,mu0n,gs,lm)
    3247  5535215076 :             gam = agamm(ws,mu0n,gs,lm)
    3248  5535215076 :             apg = alp + gam
    3249  5535215076 :             amg = alp - gam
    3250  5535215076 :             rdir(k) = apg*rdif_a(k) +  amg*(tdif_a(k)*trnlay(k) - c1)
    3251  5535215076 :             tdir(k) = apg*tdif_a(k) + (amg* rdif_a(k)-apg+c1)*trnlay(k)
    3252             :             
    3253             :             ! recalculate rdif,tdif using direct angular integration over rdir,tdir,
    3254             :             ! since Delta-Eddington rdif formula is not well-behaved (it is usually
    3255             :             ! biased low and can even be negative); use ngmax angles and gaussian
    3256             :             ! integration for most accuracy:
    3257  5535215076 :             R1 = rdif_a(k) ! use R1 as temporary
    3258  5535215076 :             T1 = tdif_a(k) ! use T1 as temporary
    3259  5535215076 :             swt = c0
    3260  5535215076 :             smr = c0
    3261  5535215076 :             smt = c0
    3262 49816935684 :             do ng=1,ngmax
    3263 44281720608 :                mu  = gauspt(ng)
    3264 44281720608 :                gwt = gauswt(ng)
    3265 44281720608 :                swt = swt + mu*gwt
    3266 44281720608 :                exp_min = min(exp_argmax,ts/mu)
    3267 44281720608 :                trn = exp(-exp_min)
    3268 44281720608 :                alp = alpha(ws,mu,gs,lm)
    3269 44281720608 :                gam = agamm(ws,mu,gs,lm)
    3270 44281720608 :                apg = alp + gam
    3271 44281720608 :                amg = alp - gam
    3272 44281720608 :                rdr = apg*R1 + amg*T1*trn - amg
    3273 44281720608 :                tdr = apg*T1 + amg*R1*trn - apg*trn + trn
    3274 44281720608 :                smr = smr + mu*rdr*gwt
    3275 49816935684 :                smt = smt + mu*tdr*gwt
    3276             :             enddo      ! ng
    3277  5535215076 :             rdif_a(k) = smr/swt
    3278  5535215076 :             tdif_a(k) = smt/swt
    3279             :             
    3280             :             ! homogeneous layer
    3281  5535215076 :             rdif_b(k) = rdif_a(k)
    3282  5535215076 :             tdif_b(k) = tdif_a(k)
    3283             : 
    3284             :             ! add fresnel layer to top of desired layer if either 
    3285             :             ! air or snow overlies ice; we ignore refraction in ice 
    3286             :             ! if a melt pond overlies it:
    3287             : 
    3288  5535215076 :             if( k == kfrsnl ) then
    3289             :                ! compute fresnel reflection and transmission amplitudes
    3290             :                ! for two polarizations: 1=perpendicular and 2=parallel to
    3291             :                ! the plane containing incident, reflected and refracted rays.
    3292             :                R1 = (mu0 - refindx*mu0n) / & 
    3293   608086511 :                     (mu0 + refindx*mu0n)
    3294             :                R2 = (refindx*mu0 - mu0n) / &
    3295   608086511 :                     (refindx*mu0 + mu0n)
    3296             :                T1 = c2*mu0 / &
    3297   608086511 :                     (mu0 + refindx*mu0n)
    3298             :                T2 = c2*mu0 / &
    3299   608086511 :                     (refindx*mu0 + mu0n)
    3300             : 
    3301             :                ! unpolarized light for direct beam
    3302   608086511 :                Rf_dir_a = p5 * (R1*R1 + R2*R2)
    3303   608086511 :                Tf_dir_a = p5 * (T1*T1 + T2*T2)*refindx*mu0n/mu0
    3304             : 
    3305             :                ! precalculated diffuse reflectivities and transmissivities
    3306             :                ! for incident radiation above and below fresnel layer, using
    3307             :                ! the direct albedos and accounting for complete internal
    3308             :                ! reflection from below; precalculated because high order
    3309             :                ! number of gaussian points (~256) is required for convergence:
    3310             : 
    3311             :                ! above
    3312   608086511 :                Rf_dif_a = cp063
    3313   608086511 :                Tf_dif_a = c1 - Rf_dif_a
    3314             :                ! below
    3315   608086511 :                Rf_dif_b = cp455
    3316   608086511 :                Tf_dif_b = c1 - Rf_dif_b
    3317             : 
    3318             :                ! the k = kfrsnl layer properties are updated to combined 
    3319             :                ! the fresnel (refractive) layer, always taken to be above
    3320             :                ! the present layer k (i.e. be the top interface):
    3321             : 
    3322   608086511 :                rintfc   = c1 / (c1-Rf_dif_b*rdif_a(k))
    3323    34671652 :                tdir(k)   = Tf_dir_a*tdir(k) + &
    3324    34671652 :                     Tf_dir_a*rdir(k) * &
    3325   608086511 :                     Rf_dif_b*rintfc*tdif_a(k)
    3326    34671652 :                rdir(k)   = Rf_dir_a + &
    3327    34671652 :                     Tf_dir_a*rdir(k) * &
    3328   608086511 :                     rintfc*Tf_dif_b
    3329    34671652 :                rdif_a(k) = Rf_dif_a + &
    3330    34671652 :                     Tf_dif_a*rdif_a(k) * &
    3331   608086511 :                     rintfc*Tf_dif_b
    3332    34671652 :                rdif_b(k) = rdif_b(k) + &
    3333    34671652 :                     tdif_b(k)*Rf_dif_b * &
    3334   608086511 :                     rintfc*tdif_a(k)
    3335   608086511 :                tdif_a(k) = tdif_a(k)*rintfc*Tf_dif_a
    3336   608086511 :                tdif_b(k) = tdif_b(k)*rintfc*Tf_dif_b
    3337             : 
    3338             :                ! update trnlay to include fresnel transmission
    3339   608086511 :                trnlay(k) = Tf_dir_a*trnlay(k)
    3340             : 
    3341             :             endif      ! k = kfrsnl
    3342             : 
    3343             :          endif ! trntdr(k) > trmin
    3344             :          
    3345             :          ! initialize current layer properties to zero; only if total
    3346             :          ! transmission to the top interface of the current layer exceeds the
    3347             :          ! minimum, will these values be computed below:
    3348             :          ! Calculate the solar beam transmission, total transmission, and
    3349             :          ! reflectivity for diffuse radiation from below at interface k, 
    3350             :          ! the top of the current layer k:
    3351             :          !
    3352             :          !              layers       interface
    3353             :          !         
    3354             :          !       ---------------------  k-1 
    3355             :          !                k-1
    3356             :          !       ---------------------  k
    3357             :          !                 k
    3358             :          !       ---------------------  
    3359             :          !       For k = klevp
    3360             :          ! note that we ignore refraction between sea ice and underlying ocean:
    3361             :          !
    3362             :          !              layers       interface
    3363             :          !
    3364             :          !       ---------------------  k-1 
    3365             :          !                k-1
    3366             :          !       ---------------------  k
    3367             :          !       \\\\\\\ ocean \\\\\\\
    3368             :          
    3369 11575415130 :          trndir(k+1) = trndir(k)*trnlay(k)
    3370 11575415130 :          refkm1         = c1/(c1 - rdndif(k)*rdif_a(k))
    3371 11575415130 :          tdrrdir        = trndir(k)*rdir(k)
    3372 11575415130 :          tdndif         = trntdr(k) - trndir(k)
    3373   728323050 :          trntdr(k+1) = trndir(k)*tdir(k) + &
    3374 12303738180 :               (tdndif + tdrrdir*rdndif(k))*refkm1*tdif_a(k)
    3375   728323050 :          rdndif(k+1) = rdif_b(k) + &
    3376 12303738180 :               (tdif_b(k)*rdndif(k)*refkm1*tdif_a(k))
    3377 12732956643 :          trndif(k+1) = trndif(k)*refkm1*tdif_a(k)
    3378             : 
    3379             :       enddo       ! k   end main level loop
    3380             : 
    3381             :       ! compute reflectivity to direct and diffuse radiation for layers 
    3382             :       ! below by adding succesive layers starting from the underlying 
    3383             :       ! ocean and working upwards:
    3384             :       !
    3385             :       !              layers       interface
    3386             :       !
    3387             :       !       ---------------------  k
    3388             :       !                 k
    3389             :       !       ---------------------  k+1
    3390             :       !                k+1
    3391             :       !       ---------------------
    3392             : 
    3393  1157541513 :       rupdir(klevp) = albodr
    3394  1157541513 :       rupdif(klevp) = albodf 
    3395             : 
    3396 12732956643 :       do k=klev,0,-1
    3397             :          ! interface scattering
    3398 11575415130 :          refkp1        = c1/( c1 - rdif_b(k)*rupdif(k+1))
    3399             :          ! dir from top layer plus exp tran ref from lower layer, interface
    3400             :          ! scattered and tran thru top layer from below, plus diff tran ref
    3401             :          ! from lower layer with interface scattering tran thru top from below
    3402   728323050 :          rupdir(k) = rdir(k) &
    3403  1456646100 :               + (        trnlay(k)  *rupdir(k+1) &
    3404 12303738180 :               +  (tdir(k)-trnlay(k))*rupdif(k+1))*refkp1*tdif_b(k)
    3405             :          ! dif from top layer from above, plus dif tran upwards reflected and
    3406             :          ! interface scattered which tran top from below
    3407 12732956643 :          rupdif(k) = rdif_a(k) + tdif_a(k)*rupdif(k+1)*refkp1*tdif_b(k)
    3408             :       enddo       ! k
    3409             : 
    3410  1157541513 :       end subroutine solution_dEdd
    3411             : 
    3412             : !=======================================================================
    3413             : !
    3414             : !   Set snow horizontal coverage, density and grain radius diagnostically 
    3415             : !   for the Delta-Eddington solar radiation method.
    3416             : !
    3417             : ! author:  Bruce P. Briegleb, NCAR 
    3418             : !   2013:  E Hunke merged with NCAR version
    3419             : 
    3420   705443757 :       subroutine shortwave_dEdd_set_snow(nslyr,    R_snw,    &
    3421             :                                          dT_mlt,   rsnw_mlt, &
    3422             :                                          aice,     vsno,     &
    3423             :                                          Tsfc,     fs,       &
    3424             :                                          hs0,      hs,       &
    3425   705443757 :                                          rhosnw,   rsnw)
    3426             : 
    3427             :       integer (kind=int_kind), intent(in) :: & 
    3428             :          nslyr      ! number of snow layers
    3429             : 
    3430             :       real (kind=dbl_kind), intent(in) :: &
    3431             :          R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
    3432             :          dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
    3433             :          rsnw_mlt  ! maximum melting snow grain radius (10^-6 m)
    3434             : 
    3435             :       real (kind=dbl_kind), intent(in) :: &
    3436             :          aice   , & ! concentration of ice
    3437             :          vsno   , & ! volume of snow
    3438             :          Tsfc   , & ! surface temperature 
    3439             :          hs0        ! snow depth for transition to bare sea ice (m)
    3440             : 
    3441             :      real (kind=dbl_kind), intent(inout) :: &
    3442             :          fs     , & ! horizontal coverage of snow
    3443             :          hs         ! snow depth
    3444             : 
    3445             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
    3446             :          rhosnw , & ! density in snow layer (kg/m3)
    3447             :          rsnw       ! grain radius in snow layer (micro-meters)
    3448             : 
    3449             :       ! local variables
    3450             : 
    3451             :       integer (kind=int_kind) :: &
    3452             :          ks           ! snow vertical index
    3453             : 
    3454             :       real (kind=dbl_kind) :: &
    3455    45664357 :          fT  , & ! piecewise linear function of surface temperature
    3456    45664357 :          dTs , & ! difference of Tsfc and Timelt
    3457    45664357 :          rsnw_nm ! actual used nonmelt snow grain radius (micro-meters)
    3458             : 
    3459             :       real (kind=dbl_kind), parameter :: &
    3460             :          ! units for the following are 1.e-6 m (micro-meters)
    3461             :          rsnw_fresh    =  100._dbl_kind, & ! freshly-fallen snow grain radius 
    3462             :          rsnw_nonmelt  =  500._dbl_kind, & ! nonmelt snow grain radius
    3463             :          rsnw_sig      =  250._dbl_kind    ! assumed sigma for snow grain radius
    3464             : 
    3465             :       character(len=*),parameter :: subname='(shortwave_dEdd_set_snow)'
    3466             : 
    3467             : !-----------------------------------------------------------------------
    3468             : 
    3469             :       ! set snow horizontal fraction
    3470   705443757 :       hs = vsno / aice
    3471             :       
    3472   705443757 :       if (hs >= hs_min) then
    3473   574334283 :          fs = c1
    3474   574334283 :          if (hs0 > puny) fs = min(hs/hs0, c1)
    3475             :       endif
    3476             :       
    3477             :       ! bare ice, temperature dependence
    3478   705443757 :       dTs = Timelt - Tsfc
    3479   705443757 :       fT  = -min(dTs/dT_mlt-c1,c0)
    3480             :       ! tune nonmelt snow grain radius if desired: note that
    3481             :       ! the sign is negative so that if R_snw is 1, then the
    3482             :       ! snow grain radius is reduced and thus albedo increased.
    3483   705443757 :       rsnw_nm = rsnw_nonmelt - R_snw*rsnw_sig
    3484   705443757 :       rsnw_nm = max(rsnw_nm, rsnw_fresh)
    3485   705443757 :       rsnw_nm = min(rsnw_nm, rsnw_mlt) 
    3486  1410887514 :       do ks = 1, nslyr
    3487             :          ! snow density ccsm3 constant value
    3488   705443757 :          rhosnw(ks) = rhos
    3489             :          ! snow grain radius between rsnw_nonmelt and rsnw_mlt
    3490   705443757 :          rsnw(ks) = rsnw_nm + (rsnw_mlt-rsnw_nm)*fT
    3491   705443757 :          rsnw(ks) = max(rsnw(ks), rsnw_fresh)
    3492  1410887514 :          rsnw(ks) = min(rsnw(ks), rsnw_mlt)
    3493             :       enddo        ! ks
    3494             : 
    3495   705443757 :       end subroutine shortwave_dEdd_set_snow
    3496             : 
    3497             : !=======================================================================
    3498             : !
    3499             : !   Set pond fraction and depth diagnostically for
    3500             : !   the Delta-Eddington solar radiation method.
    3501             : !
    3502             : ! author:  Bruce P. Briegleb, NCAR 
    3503             : !   2013:  E Hunke merged with NCAR version
    3504             : 
    3505    76623754 :       subroutine shortwave_dEdd_set_pond(Tsfc,               &
    3506             :                                          fs,       fp,       &
    3507             :                                          hp)
    3508             : 
    3509             :       real (kind=dbl_kind), intent(in) :: &
    3510             :          Tsfc   , & ! surface temperature
    3511             :          fs         ! horizontal coverage of snow
    3512             : 
    3513             :       real (kind=dbl_kind), intent(out) :: &
    3514             :          fp     , & ! pond fractional coverage (0 to 1)
    3515             :          hp         ! pond depth (m)
    3516             : 
    3517             :       ! local variables
    3518             : 
    3519             :       real (kind=dbl_kind) :: &
    3520    13444925 :          fT  , & ! piecewise linear function of surface temperature
    3521    13444925 :          dTs     ! difference of Tsfc and Timelt
    3522             : 
    3523             :       real (kind=dbl_kind), parameter :: &
    3524             :          dT_pnd = c1   ! change in temp for pond fraction and depth
    3525             : 
    3526             :       character(len=*),parameter :: subname='(shortwave_dEdd_set_pond)'
    3527             : 
    3528             : !-----------------------------------------------------------------------
    3529             : 
    3530             :       ! bare ice, temperature dependence
    3531    76623754 :       dTs = Timelt - Tsfc
    3532    76623754 :       fT  = -min(dTs/dT_pnd-c1,c0)
    3533             :       ! pond
    3534    76623754 :       fp = 0.3_dbl_kind*fT*(c1-fs)
    3535    76623754 :       hp = 0.3_dbl_kind*fT*(c1-fs)
    3536             :       
    3537    76623754 :       end subroutine shortwave_dEdd_set_pond
    3538             : 
    3539             : ! End Delta-Eddington shortwave method
    3540             : 
    3541             : !=======================================================================
    3542             : !
    3543             : ! authors     Nicole Jeffery, LANL
    3544             : 
    3545           0 :       subroutine compute_shortwave_trcr(nslyr,               &
    3546           0 :                                     bgcN,         zaero,     &
    3547           0 :                                     trcrn_bgcsw,             &
    3548           0 :                                     sw_grid,      hin,       &
    3549             :                                     hbri,                    &
    3550             :                                     nilyr,        nblyr,     &
    3551           0 :                                     i_grid,                  &
    3552             :                                     skl_bgc,      z_tracers  )
    3553             :       
    3554             :       integer (kind=int_kind), intent(in) :: &
    3555             :          nslyr          ! number of snow layers
    3556             : 
    3557             :       integer (kind=int_kind), intent(in) :: &
    3558             :          nblyr      , & ! number of bio layers
    3559             :          nilyr          ! number of ice layers 
    3560             : 
    3561             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    3562             :          bgcN       , & ! Nit tracer
    3563             :          zaero          ! zaero tracer
    3564             : 
    3565             :       real (kind=dbl_kind), dimension (:), intent(out):: &
    3566             :          trcrn_bgcsw    ! ice on shortwave grid tracers
    3567             : 
    3568             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    3569             :          sw_grid     , & ! 
    3570             :          i_grid          ! CICE bio grid
    3571             :          
    3572             :       real(kind=dbl_kind), intent(in) :: &
    3573             :          hin         , & ! CICE ice thickness
    3574             :          hbri            ! brine height
    3575             : 
    3576             :       logical (kind=log_kind), intent(in) :: &
    3577             :          skl_bgc     , & ! skeletal layer bgc
    3578             :          z_tracers       ! zbgc
    3579             : 
    3580             :       !  local variables
    3581             : 
    3582             :       integer (kind=int_kind) :: k, n, nn
    3583             : 
    3584             :       real (kind=dbl_kind), dimension (ntrcr+2) :: &
    3585           0 :          trtmp0, &      ! temporary, remapped tracers
    3586           0 :          trtmp
    3587             : 
    3588             :       real (kind=dbl_kind), dimension (nilyr+1):: &
    3589           0 :          icegrid        ! correct for large ice surface layers
    3590             : 
    3591             :       real (kind=dbl_kind):: &
    3592           0 :          top_conc       ! 1% (min_bgc) of surface concentration 
    3593             :                         ! when hin > hbri:  just used in sw calculation
    3594             : 
    3595             :       character(len=*),parameter :: subname='(compute_shortwave_trcr)'
    3596             : 
    3597             :       !-----------------------------------------------------------------
    3598             :       ! Compute aerosols and algal chlorophyll on shortwave grid
    3599             :       !-----------------------------------------------------------------
    3600             : 
    3601           0 :       trtmp0(:) = c0
    3602           0 :       trtmp(:) = c0
    3603           0 :       trcrn_bgcsw(:) = c0
    3604             : 
    3605           0 :       do k = 1,nilyr+1
    3606           0 :          icegrid(k) = sw_grid(k)
    3607             :       enddo    
    3608           0 :       if (sw_grid(1)*hin*c2 > hi_ssl) then
    3609           0 :          icegrid(1) = hi_ssl/c2/hin
    3610             :       endif
    3611             : 
    3612           0 :       if (z_tracers) then
    3613           0 :       if (tr_bgc_N)  then
    3614           0 :          if (size(bgcN) < n_algae*(nblyr+3)) then
    3615           0 :             call icepack_warnings_add(subname//' ERROR: size(bgcN) too small')
    3616           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    3617           0 :             return
    3618             :          endif
    3619             : 
    3620           0 :          do k = 1, nblyr+1
    3621           0 :             do n = 1, n_algae
    3622           0 :                trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) + &
    3623           0 :                                 R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1 + k-1)
    3624             :             enddo ! n
    3625             :          enddo    ! k
    3626             :  
    3627           0 :          top_conc = trtmp0(nt_bgc_N(1))*min_bgc
    3628             :          call remap_zbgc (nilyr+1, &
    3629             :                           nt_bgc_N(1),                &
    3630           0 :                           trtmp0(1:ntrcr  ),          &
    3631           0 :                           trtmp (1:ntrcr+2),          &
    3632             :                           1,                 nblyr+1, &
    3633             :                           hin,               hbri,    &
    3634           0 :                           icegrid(1:nilyr+1),         &
    3635           0 :                           i_grid(1:nblyr+1), top_conc ) 
    3636           0 :          if (icepack_warnings_aborted(subname)) return
    3637             : 
    3638           0 :          do k = 1, nilyr+1
    3639           0 :             trcrn_bgcsw(nlt_chl_sw+nslyr+k) = trtmp(nt_bgc_N(1) + k-1)
    3640             :          enddo       ! k
    3641             : 
    3642           0 :          do n = 1, n_algae   ! snow contribution
    3643           0 :             trcrn_bgcsw(nlt_chl_sw)= trcrn_bgcsw(nlt_chl_sw) &
    3644           0 :                      + R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+1)
    3645             :                               ! snow surface layer
    3646           0 :             trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) = &
    3647           0 :                      trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) &
    3648           0 :                      + R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+2)
    3649             :                               ! only 1 snow layer in zaero
    3650             :          enddo ! n
    3651             :       endif    ! tr_bgc_N
    3652             : 
    3653           0 :       if (tr_zaero) then
    3654           0 :          if (size(zaero) < n_zaero*(nblyr+3)) then
    3655           0 :             call icepack_warnings_add(subname//' ERROR: size(zaero) too small')
    3656           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    3657           0 :             return
    3658             :          endif
    3659             : 
    3660           0 :          do n = 1, n_zaero
    3661             : 
    3662           0 :             trtmp0(:) = c0
    3663           0 :             trtmp(:) = c0
    3664             : 
    3665           0 :             do k = 1, nblyr+1
    3666           0 :                trtmp0(nt_zaero(n) + k-1) = zaero(nt_zaero(n)-nt_zaero(1)+1+k-1)
    3667             :             enddo
    3668             : 
    3669           0 :             top_conc = trtmp0(nt_zaero(n))*min_bgc
    3670             :             call remap_zbgc (nilyr+1, &
    3671           0 :                              nt_zaero(n),                &
    3672           0 :                              trtmp0(1:ntrcr  ),          &
    3673           0 :                              trtmp (1:ntrcr+2),          &
    3674             :                              1,                 nblyr+1, &
    3675             :                              hin,               hbri,    &
    3676           0 :                              icegrid(1:nilyr+1),         &
    3677           0 :                              i_grid(1:nblyr+1), top_conc )
    3678           0 :             if (icepack_warnings_aborted(subname)) return
    3679             : 
    3680           0 :             do k = 1,nilyr+1
    3681           0 :                trcrn_bgcsw(nlt_zaero_sw(n)+nslyr+k) = trtmp(nt_zaero(n) + k-1)
    3682             :             enddo
    3683           0 :             trcrn_bgcsw(nlt_zaero_sw(n))= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+1) !snow ssl
    3684           0 :             trcrn_bgcsw(nlt_zaero_sw(n)+1:nlt_zaero_sw(n)+nslyr)= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+2)
    3685             :          enddo ! n
    3686             :       endif    ! tr_zaero
    3687           0 :       elseif (skl_bgc) then
    3688             : 
    3689           0 :          do nn = 1,n_algae
    3690           0 :             trcrn_bgcsw(nbtrcr_sw) = trcrn_bgcsw(nbtrcr_sw) &
    3691           0 :                                 + F_abs_chl(nn)*R_chl2N(nn) &
    3692           0 :                                 * bgcN(nt_bgc_N(nn)-nt_bgc_N(1)+1)*sk_l/hin &
    3693           0 :                                 * real(nilyr,kind=dbl_kind)
    3694             :          enddo 
    3695             : 
    3696             :       endif
    3697             :       end subroutine compute_shortwave_trcr
    3698             : 
    3699             : !=======================================================================
    3700             : !autodocument_start icepack_prep_radiation
    3701             : ! Scales radiation fields computed on the previous time step.
    3702             : !
    3703             : ! authors: Elizabeth Hunke, LANL
    3704             : 
    3705   957276648 :       subroutine icepack_prep_radiation (ncat, nilyr, nslyr,   &
    3706   957276648 :                                         aice,        aicen,    &
    3707             :                                         swvdr,       swvdf,    &
    3708             :                                         swidr,       swidf,    &
    3709             :                                         alvdr_ai,    alvdf_ai, &
    3710             :                                         alidr_ai,    alidf_ai, &
    3711             :                                         scale_factor,          &
    3712   957276648 :                                         fswsfcn,     fswintn,  &
    3713   957276648 :                                         fswthrun,    fswpenln, &
    3714   957276648 :                                         Sswabsn,     Iswabsn)
    3715             : 
    3716             :       integer (kind=int_kind), intent(in) :: &
    3717             :          ncat    , & ! number of ice thickness categories
    3718             :          nilyr   , & ! number of ice layers
    3719             :          nslyr       ! number of snow layers
    3720             : 
    3721             :       real (kind=dbl_kind), intent(in) :: &
    3722             :          aice        , & ! ice area fraction
    3723             :          swvdr       , & ! sw down, visible, direct  (W/m^2)
    3724             :          swvdf       , & ! sw down, visible, diffuse (W/m^2)
    3725             :          swidr       , & ! sw down, near IR, direct  (W/m^2)
    3726             :          swidf       , & ! sw down, near IR, diffuse (W/m^2)
    3727             :          ! grid-box-mean albedos aggregated over categories (if calc_Tsfc)
    3728             :          alvdr_ai    , & ! visible, direct   (fraction)
    3729             :          alidr_ai    , & ! near-ir, direct   (fraction)
    3730             :          alvdf_ai    , & ! visible, diffuse  (fraction)
    3731             :          alidf_ai        ! near-ir, diffuse  (fraction)
    3732             : 
    3733             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    3734             :          aicen           ! ice area fraction in each category
    3735             : 
    3736             :       real (kind=dbl_kind), intent(inout) :: &
    3737             :          scale_factor    ! shortwave scaling factor, ratio new:old
    3738             : 
    3739             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    3740             :          fswsfcn     , & ! SW absorbed at ice/snow surface (W m-2)
    3741             :          fswintn     , & ! SW absorbed in ice interior, below surface (W m-2)
    3742             :          fswthrun        ! SW through ice to ocean (W/m^2)
    3743             : 
    3744             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    3745             :          fswpenln    , & ! visible SW entering ice layers (W m-2)
    3746             :          Iswabsn     , & ! SW radiation absorbed in ice layers (W m-2)
    3747             :          Sswabsn         ! SW radiation absorbed in snow layers (W m-2)
    3748             : 
    3749             : !autodocument_end
    3750             : 
    3751             :       ! local variables
    3752             : 
    3753             :       integer (kind=int_kind) :: &
    3754             :          k           , & ! vertical index       
    3755             :          n               ! thickness category index
    3756             : 
    3757    57986832 :       real (kind=dbl_kind) :: netsw 
    3758             : 
    3759             :       character(len=*),parameter :: subname='(icepack_prep_radiation)'
    3760             : 
    3761             :       !-----------------------------------------------------------------
    3762             :       ! Compute netsw scaling factor (new netsw / old netsw)
    3763             :       !-----------------------------------------------------------------
    3764             : 
    3765   957276648 :          if (aice > c0 .and. scale_factor > puny) then
    3766             :             netsw = swvdr*(c1 - alvdr_ai) &
    3767             :                   + swvdf*(c1 - alvdf_ai) &
    3768             :                   + swidr*(c1 - alidr_ai) &
    3769    75168313 :                   + swidf*(c1 - alidf_ai)
    3770    75168313 :             scale_factor = netsw / scale_factor
    3771             :          else
    3772   882108335 :             scale_factor = c1
    3773             :          endif
    3774             : 
    3775  5611698288 :          do n = 1, ncat
    3776             : 
    3777  5611698288 :             if (aicen(n) > puny) then
    3778             : 
    3779             :       !-----------------------------------------------------------------
    3780             :       ! Scale absorbed solar radiation for change in net shortwave
    3781             :       !-----------------------------------------------------------------
    3782             : 
    3783   744838517 :                fswsfcn(n)  = scale_factor*fswsfcn (n)
    3784   744838517 :                fswintn(n)  = scale_factor*fswintn (n)
    3785   744838517 :                fswthrun(n) = scale_factor*fswthrun(n)
    3786  6023820495 :                do k = 1,nilyr+1
    3787  6023820495 :                   fswpenln(k,n) = scale_factor*fswpenln(k,n)
    3788             :                enddo       !k
    3789  1489677034 :                do k=1,nslyr
    3790  1489677034 :                   Sswabsn(k,n) = scale_factor*Sswabsn(k,n)
    3791             :                enddo
    3792  5278981978 :                do k=1,nilyr
    3793  5278981978 :                   Iswabsn(k,n) = scale_factor*Iswabsn(k,n)
    3794             :                enddo
    3795             : 
    3796             :             endif
    3797             :          enddo                  ! ncat
    3798             : 
    3799   957276648 :       end subroutine icepack_prep_radiation
    3800             : 
    3801             : !=======================================================================
    3802             : !autodocument_start icepack_step_radiation
    3803             : ! Computes radiation fields
    3804             : !
    3805             : ! authors: William H. Lipscomb, LANL
    3806             : !          David Bailey, NCAR
    3807             : !          Elizabeth C. Hunke, LANL
    3808             : 
    3809   725368070 :       subroutine icepack_step_radiation (dt,       ncat,     &
    3810             :                                         nblyr,               &
    3811             :                                         nilyr,    nslyr,     &
    3812             :                                         dEdd_algae,          &
    3813   725368070 :                                         swgrid,   igrid,     &
    3814   725368070 :                                         fbri,                &
    3815   725368070 :                                         aicen,    vicen,     &
    3816   725368070 :                                         vsnon,    Tsfcn,     &
    3817   725368070 :                                         alvln,    apndn,     &
    3818   725368070 :                                         hpndn,    ipndn,     &
    3819   725368070 :                                         aeron,               &
    3820   725368070 :                                         bgcNn,    zaeron,    &
    3821   725368070 :                                         trcrn_bgcsw,         &
    3822             :                                         TLAT,     TLON,      &
    3823           0 :                                         calendar_type,       &
    3824             :                                         days_per_year,       &
    3825             :                                         nextsw_cday,         &
    3826             :                                         yday,     sec,       &
    3827   725368070 :                                         kaer_tab, waer_tab,  &
    3828   725368070 :                                         gaer_tab,            &
    3829   725368070 :                                         kaer_bc_tab,         &
    3830   725368070 :                                         waer_bc_tab,         &
    3831   725368070 :                                         gaer_bc_tab,         &
    3832   725368070 :                                         bcenh,               &
    3833             :                                         modal_aero,          &
    3834             :                                         swvdr,    swvdf,     &
    3835             :                                         swidr,    swidf,     &
    3836             :                                         coszen,   fsnow,     &
    3837  1450736140 :                                         alvdrn,   alvdfn,    &
    3838  1450736140 :                                         alidrn,   alidfn,    &
    3839   725368070 :                                         fswsfcn,  fswintn,   &
    3840   725368070 :                                         fswthrun, fswpenln,  &
    3841   725368070 :                                         Sswabsn,  Iswabsn,   &
    3842   725368070 :                                         albicen,  albsnon,   &
    3843   725368070 :                                         albpndn,  apeffn,    &
    3844   725368070 :                                         snowfracn,           &
    3845   725368070 :                                         dhsn,     ffracn,    &
    3846             :                                         l_print_point, &
    3847             :                                         initonly)
    3848             : 
    3849             :       integer (kind=int_kind), intent(in) :: &
    3850             :          ncat      , & ! number of ice thickness categories
    3851             :          nilyr     , & ! number of ice layers
    3852             :          nslyr     , & ! number of snow layers
    3853             :          nblyr         ! number of bgc layers
    3854             : 
    3855             :       real (kind=dbl_kind), intent(in) :: &
    3856             :          dt        , & ! time step (s)
    3857             :          swvdr     , & ! sw down, visible, direct  (W/m^2)
    3858             :          swvdf     , & ! sw down, visible, diffuse (W/m^2)
    3859             :          swidr     , & ! sw down, near IR, direct  (W/m^2)
    3860             :          swidf     , & ! sw down, near IR, diffuse (W/m^2)
    3861             :          fsnow     , & ! snowfall rate (kg/m^2 s)
    3862             :          TLAT, TLON    ! latitude and longitude (radian)
    3863             : 
    3864             :       character (len=char_len), intent(in) :: &
    3865             :          calendar_type       ! differentiates Gregorian from other calendars
    3866             : 
    3867             :       integer (kind=int_kind), intent(in) :: &
    3868             :          days_per_year, &    ! number of days in one year
    3869             :          sec                 ! elapsed seconds into date
    3870             : 
    3871             :       real (kind=dbl_kind), intent(in) :: &
    3872             :          nextsw_cday     , & ! julian day of next shortwave calculation
    3873             :          yday                ! day of the year
    3874             : 
    3875             :       real (kind=dbl_kind), intent(inout) :: &
    3876             :          coszen        ! cosine solar zenith angle, < 0 for sun below horizon 
    3877             : 
    3878             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    3879             :          igrid              ! biology vertical interface points
    3880             :  
    3881             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    3882             :          swgrid                ! grid for ice tracers used in dEdd scheme
    3883             :         
    3884             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & 
    3885             :          kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
    3886             :          waer_tab, & ! aerosol single scatter albedo (fraction)
    3887             :          gaer_tab    ! aerosol asymmetry parameter (cos(theta))
    3888             : 
    3889             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: & 
    3890             :          kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
    3891             :          waer_bc_tab, & ! aerosol single scatter albedo (fraction)
    3892             :          gaer_bc_tab    ! aerosol asymmetry parameter (cos(theta))
    3893             : 
    3894             :       real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & 
    3895             :          bcenh 
    3896             : 
    3897             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    3898             :          aicen     , & ! ice area fraction in each category
    3899             :          vicen     , & ! ice volume in each category (m)
    3900             :          vsnon     , & ! snow volume in each category (m)
    3901             :          Tsfcn     , & ! surface temperature (deg C)
    3902             :          alvln     , & ! level-ice area fraction
    3903             :          apndn     , & ! pond area fraction
    3904             :          hpndn     , & ! pond depth (m)
    3905             :          ipndn     , & ! pond refrozen lid thickness (m)
    3906             :          fbri           ! brine fraction 
    3907             : 
    3908             :       real(kind=dbl_kind), dimension(:,:), intent(in) :: &
    3909             :          aeron     , & ! aerosols (kg/m^3)
    3910             :          bgcNn     , & ! bgc Nit tracers
    3911             :          zaeron        ! bgcz aero tracers
    3912             : 
    3913             :       real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
    3914             :          trcrn_bgcsw   ! zaerosols (kg/m^3) and chla (mg/m^3)
    3915             : 
    3916             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    3917             :          alvdrn    , & ! visible, direct  albedo (fraction)
    3918             :          alidrn    , & ! near-ir, direct   (fraction)
    3919             :          alvdfn    , & ! visible, diffuse  (fraction)
    3920             :          alidfn    , & ! near-ir, diffuse  (fraction)
    3921             :          fswsfcn   , & ! SW absorbed at ice/snow surface (W m-2)
    3922             :          fswintn   , & ! SW absorbed in ice interior, below surface (W m-2)
    3923             :          fswthrun  , & ! SW through ice to ocean (W/m^2)
    3924             :          snowfracn , & ! snow fraction on each category
    3925             :          dhsn      , & ! depth difference for snow on sea ice and pond ice
    3926             :          ffracn    , & ! fraction of fsurfn used to melt ipond
    3927             :                        ! albedo components for history
    3928             :          albicen   , & ! bare ice 
    3929             :          albsnon   , & ! snow 
    3930             :          albpndn   , & ! pond 
    3931             :          apeffn        ! effective pond area used for radiation calculation
    3932             : 
    3933             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    3934             :          fswpenln  , & ! visible SW entering ice layers (W m-2)
    3935             :          Iswabsn   , & ! SW radiation absorbed in ice layers (W m-2)
    3936             :          Sswabsn       ! SW radiation absorbed in snow layers (W m-2)
    3937             : 
    3938             :       logical (kind=log_kind), intent(in) :: &
    3939             :          l_print_point, & ! flag for printing diagnostics
    3940             :          dEdd_algae   , & ! .true. use prognostic chla in dEdd
    3941             :          modal_aero       ! .true. use modal aerosol optical treatment
    3942             : 
    3943             :       logical (kind=log_kind), optional :: &
    3944             :          initonly         ! flag to indicate init only, default is false
    3945             : 
    3946             : !autodocument_end
    3947             : 
    3948             :       ! local variables
    3949             : 
    3950             :       integer (kind=int_kind) :: &
    3951             :          n                  ! thickness category index
    3952             : 
    3953             :       logical (kind=log_kind) :: &
    3954             :          linitonly       ! local flag for initonly
    3955             : 
    3956             :       real(kind=dbl_kind) :: &
    3957    45450554 :         hin,         & ! Ice thickness (m)
    3958    45450554 :         hbri           ! brine thickness (m)
    3959             : 
    3960             :       character(len=*),parameter :: subname='(icepack_step_radiation)'
    3961             : 
    3962   725368070 :         hin = c0
    3963   725368070 :         hbri = c0
    3964   725368070 :         linitonly = .false.
    3965   725368070 :         if (present(initonly)) then
    3966     2893598 :            linitonly = initonly
    3967             :         endif
    3968             : 
    3969             :          ! Initialize
    3970  4282652292 :          do n = 1, ncat
    3971  3557284222 :             alvdrn  (n) = c0
    3972  3557284222 :             alidrn  (n) = c0
    3973  3557284222 :             alvdfn  (n) = c0
    3974  3557284222 :             alidfn  (n) = c0
    3975  3557284222 :             fswsfcn (n) = c0
    3976  3557284222 :             fswintn (n) = c0
    3977  4282652292 :             fswthrun(n) = c0
    3978             :          enddo   ! ncat
    3979 31331295668 :          fswpenln (:,:) = c0
    3980 27774011446 :          Iswabsn  (:,:) = c0
    3981  7839936514 :          Sswabsn  (:,:) = c0
    3982  4282652292 :          trcrn_bgcsw(:,:) = c0
    3983             : 
    3984             :          ! Interpolate z-shortwave tracers to shortwave grid
    3985   725368070 :          if (dEdd_algae) then
    3986           0 :          do n = 1, ncat
    3987           0 :               if (aicen(n) .gt. puny) then
    3988           0 :                  hin = vicen(n)/aicen(n)
    3989           0 :                  hbri= fbri(n)*hin
    3990             :                  call compute_shortwave_trcr(nslyr,           &
    3991           0 :                                      bgcNn(:,n),              &
    3992           0 :                                      zaeron(:,n),             &
    3993           0 :                                      trcrn_bgcsw(:,n),        &
    3994           0 :                                      swgrid,       hin,       &
    3995             :                                      hbri,                    &
    3996             :                                      nilyr,        nblyr,     &
    3997           0 :                                      igrid,                   &
    3998           0 :                                      skl_bgc,      z_tracers  )
    3999           0 :                  if (icepack_warnings_aborted(subname)) return
    4000             :               endif
    4001             :          enddo
    4002             :          endif
    4003             : 
    4004   725368070 :          if (calc_Tsfc) then
    4005   702182694 :          if (trim(shortwave) == 'dEdd') then ! delta Eddington
    4006             :             
    4007             :             call run_dEdd(dt,           ncat,           &
    4008             :                           dEdd_algae,                   &
    4009             :                           nilyr,        nslyr,          &
    4010           0 :                           aicen,        vicen,          &
    4011           0 :                           vsnon,        Tsfcn,          &
    4012           0 :                           alvln,        apndn,          &
    4013           0 :                           hpndn,        ipndn,          &
    4014           0 :                           aeron,        kalg,           &
    4015           0 :                           trcrn_bgcsw,                  &
    4016             :                           heat_capacity,                &
    4017             :                           TLAT,         TLON,           &
    4018             :                           calendar_type,days_per_year,  &
    4019             :                           nextsw_cday,  yday,           &
    4020             :                           sec,          R_ice,          &
    4021             :                           R_pnd,        R_snw,          &
    4022             :                           dT_mlt,       rsnw_mlt,       &
    4023             :                           hs0,          hs1,            &
    4024             :                           hp1,          pndaspect,      &
    4025           0 :                           kaer_tab,     waer_tab,       &
    4026           0 :                           gaer_tab,                     &
    4027           0 :                           kaer_bc_tab,                  &
    4028           0 :                           waer_bc_tab,                  &
    4029           0 :                           gaer_bc_tab,                  &
    4030           0 :                           bcenh,                        &
    4031             :                           modal_aero,                   &
    4032             :                           swvdr,        swvdf,          &
    4033             :                           swidr,        swidf,          &
    4034             :                           coszen,       fsnow,          &
    4035           0 :                           alvdrn,       alvdfn,         &
    4036           0 :                           alidrn,       alidfn,         &
    4037           0 :                           fswsfcn,      fswintn,        &
    4038           0 :                           fswthrun,     fswpenln,       &
    4039           0 :                           Sswabsn,      Iswabsn,        &
    4040           0 :                           albicen,      albsnon,        &
    4041           0 :                           albpndn,      apeffn,         &
    4042           0 :                           snowfracn,                    &
    4043           0 :                           dhsn,         ffracn,         &
    4044             :                           l_print_point,                &
    4045   655811942 :                           linitonly)
    4046   655811942 :             if (icepack_warnings_aborted(subname)) return
    4047             :  
    4048    46370752 :          elseif (trim(shortwave) == 'ccsm3') then
    4049             : 
    4050           0 :             call shortwave_ccsm3(aicen,      vicen,      &
    4051           0 :                                  vsnon,                  &
    4052           0 :                                  Tsfcn,                  &
    4053             :                                  swvdr,      swvdf,      &
    4054             :                                  swidr,      swidf,      &
    4055             :                                  heat_capacity,          &
    4056             :                                  albedo_type,            &
    4057             :                                  albicev,    albicei,    &
    4058             :                                  albsnowv,   albsnowi,   &
    4059             :                                  ahmax,                  &
    4060           0 :                                  alvdrn,     alidrn,     &
    4061           0 :                                  alvdfn,     alidfn,     &
    4062           0 :                                  fswsfcn,    fswintn,    &
    4063           0 :                                  fswthrun,               &
    4064           0 :                                  fswpenln,               &
    4065           0 :                                  Iswabsn,                &
    4066           0 :                                  Sswabsn,                &
    4067           0 :                                  albicen,    albsnon,    &
    4068             :                                  coszen,     ncat,       &
    4069    46370752 :                                  nilyr)
    4070    46370752 :             if (icepack_warnings_aborted(subname)) return
    4071             : 
    4072             :          else
    4073             : 
    4074           0 :             call icepack_warnings_add(subname//' ERROR: shortwave '//trim(shortwave)//' unknown')
    4075           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    4076           0 :             return
    4077             : 
    4078             :          endif   ! shortwave
    4079             : 
    4080             :       else    ! .not. calc_Tsfc
    4081             : 
    4082             :       ! Calculate effective pond area for HadGEM
    4083             : 
    4084    23185376 :          if (tr_pond_topo) then
    4085   162297632 :             do n = 1, ncat
    4086   139112256 :                apeffn(n) = c0 
    4087   162297632 :                if (aicen(n) > puny) then
    4088             :                ! Lid effective if thicker than hp1
    4089    30050578 :                  if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
    4090       11140 :                     apeffn(n) = apndn(n)
    4091             :                  else
    4092    30039438 :                     apeffn(n) = c0
    4093             :                  endif
    4094    30050578 :                  if (apndn(n) < puny) apeffn(n) = c0
    4095             :                endif
    4096             :             enddo  ! ncat
    4097             :  
    4098             :          endif ! tr_pond_topo
    4099             : 
    4100             :          ! Initialize for safety
    4101   162297632 :          do n = 1, ncat
    4102   139112256 :             alvdrn(n) = c0
    4103   139112256 :             alidrn(n) = c0
    4104   139112256 :             alvdfn(n) = c0
    4105   139112256 :             alidfn(n) = c0
    4106   139112256 :             fswsfcn(n) = c0
    4107   139112256 :             fswintn(n) = c0
    4108   162297632 :             fswthrun(n) = c0
    4109             :          enddo   ! ncat
    4110  1136083424 :          Iswabsn(:,:) = c0
    4111   301409888 :          Sswabsn(:,:) = c0
    4112             : 
    4113             :       endif    ! calc_Tsfc
    4114             : 
    4115             :       end subroutine icepack_step_radiation
    4116             : 
    4117             :       ! Delta-Eddington solution expressions
    4118             : 
    4119             : !=======================================================================
    4120             : 
    4121 49816935684 :       real(kind=dbl_kind) function alpha(w,uu,gg,e)
    4122             : 
    4123             :       real(kind=dbl_kind), intent(in) :: w, uu, gg, e
    4124             : 
    4125 49816935684 :       alpha = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
    4126             : 
    4127 49816935684 :       end function alpha
    4128             : 
    4129             : !=======================================================================
    4130             : 
    4131 49816935684 :       real(kind=dbl_kind) function agamm(w,uu,gg,e)
    4132             : 
    4133             :       real(kind=dbl_kind), intent(in) :: w, uu, gg, e
    4134             : 
    4135 49816935684 :       agamm = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
    4136             : 
    4137 49816935684 :       end function agamm
    4138             : 
    4139             : !=======================================================================
    4140             : 
    4141  5535215076 :       real(kind=dbl_kind) function n(uu,et)
    4142             : 
    4143             :       real(kind=dbl_kind), intent(in) :: uu, et
    4144             : 
    4145  5535215076 :       n = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
    4146             : 
    4147  5535215076 :       end function n
    4148             : 
    4149             : !=======================================================================
    4150             : 
    4151  5535215076 :       real(kind=dbl_kind) function u(w,gg,e)
    4152             : 
    4153             :       real(kind=dbl_kind), intent(in) :: w, gg, e
    4154             : 
    4155  5535215076 :       u = c1p5*(c1 - w*gg)/e
    4156             : 
    4157  5535215076 :       end function u
    4158             : 
    4159             : !=======================================================================
    4160             : 
    4161  5535215076 :       real(kind=dbl_kind) function el(w,gg)
    4162             : 
    4163             :       real(kind=dbl_kind), intent(in) :: w, gg
    4164             : 
    4165  5535215076 :       el = sqrt(c3*(c1-w)*(c1 - w*gg))
    4166             : 
    4167  5535215076 :       end function el
    4168             : 
    4169             : !=======================================================================
    4170             : 
    4171  5535215076 :       real(kind=dbl_kind) function taus(w,f,t)
    4172             : 
    4173             :       real(kind=dbl_kind), intent(in) :: w, f, t
    4174             : 
    4175  5535215076 :       taus = (c1 - w*f)*t
    4176             : 
    4177  5535215076 :       end function taus
    4178             : 
    4179             : !=======================================================================
    4180             : 
    4181  5535215076 :       real(kind=dbl_kind) function omgs(w,f)
    4182             : 
    4183             :       real(kind=dbl_kind), intent(in) :: w, f
    4184             : 
    4185  5535215076 :       omgs = (c1 - f)*w/(c1 - w*f)
    4186             : 
    4187  5535215076 :       end function omgs
    4188             : 
    4189             : !=======================================================================
    4190             : 
    4191  5535215076 :       real(kind=dbl_kind) function asys(gg,f)
    4192             : 
    4193             :       real(kind=dbl_kind), intent(in) :: gg, f
    4194             : 
    4195  5535215076 :       asys = (gg - f)/(c1 - f)
    4196             : 
    4197  5535215076 :       end function asys
    4198             :  
    4199             : !=======================================================================
    4200             : 
    4201             :       end module icepack_shortwave
    4202             : 
    4203             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd