!=======================================================================
!
! The albedo and absorbed/transmitted flux parameterizations for
! snow over ice, bare ice and ponded ice.
!
! Presently, two methods are included:
! (1) CCSM3
! (2) Delta-Eddington
! as two distinct routines.
! Either can be called from the ice driver.
!
! The Delta-Eddington method is described here:
!
! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple
! Scattering Parameterization for Solar Radiation in the Sea Ice
! Component of the Community Climate System Model, NCAR Technical
! Note NCAR/TN-472+STR February 2007
!
! name: originally ice_albedo
!
! authors: Bruce P. Briegleb, NCAR
! Elizabeth C. Hunke and William H. Lipscomb, LANL
! 2005, WHL: Moved absorbed_solar from icepack_therm_vertical to this
! module and changed name from ice_albedo
! 2006, WHL: Added Delta Eddington routines from Bruce Briegleb
! 2006, ECH: Changed data statements in Delta Eddington routines (no
! longer hardwired)
! Converted to free source form (F90)
! 2007, BPB: Completely updated Delta-Eddington code, so that:
! (1) multiple snow layers enabled (i.e. nslyr > 1)
! (2) included SSL for snow surface absorption
! (3) added Sswabs for internal snow layer absorption
! (4) variable sea ice layers allowed (i.e. not hardwired)
! (5) updated all inherent optical properties
! (6) included algae absorption for sea ice lowest layer
! (7) very complete internal documentation included
! 2007, ECH: Improved efficiency
! 2008, BPB: Added aerosols to Delta Eddington code
! 2013, ECH: merged with NCAR version, cleaned up
module icepack_shortwave
use icepack_kinds
use icepack_parameters, only: c0, c1, c1p5, c2, c3, c4, c10
use icepack_parameters, only: p01, p1, p15, p25, p5, p75, puny
use icepack_parameters, only: albocn, Timelt, snowpatch, awtvdr, awtidr, awtvdf, awtidf
use icepack_parameters, only: kappav, hs_min, rhofresh, rhos, nspint
use icepack_parameters, only: hi_ssl, hs_ssl, min_bgc, sk_l
use icepack_parameters, only: z_tracers, skl_bgc, calc_tsfc, shortwave, kalg, heat_capacity
use icepack_parameters, only: r_ice, r_pnd, r_snw, dt_mlt, rsnw_mlt, hs0, hs1, hp1
use icepack_parameters, only: pndaspect, albedo_type, albicev, albicei, albsnowv, albsnowi, ahmax
use icepack_tracers, only: ntrcr, nbtrcr_sw
use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
use icepack_tracers, only: tr_bgc_N, tr_aero
use icepack_tracers, only: nt_bgc_N, nt_zaero, tr_bgc_N
use icepack_tracers, only: tr_zaero, nlt_chl_sw, nlt_zaero_sw
use icepack_tracers, only: n_algae, n_aero, n_zaero
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_zbgc_shared,only: R_chl2N, F_abs_chl
use icepack_zbgc_shared,only: remap_zbgc
use icepack_orbital, only: compute_coszen
implicit none
private
public :: run_dEdd, &
shortwave_ccsm3, &
compute_shortwave_trcr, &
icepack_prep_radiation, &
icepack_step_radiation
real (kind=dbl_kind), parameter :: &
hpmin = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
hp0 = 0.200_dbl_kind ! pond depth below which transition to bare ice
real (kind=dbl_kind), parameter :: &
exp_argmax = c10 ! maximum argument of exponential
!=======================================================================
contains
!=======================================================================
!
! Driver for basic solar radiation from CCSM3. Albedos and absorbed solar.
subroutine shortwave_ccsm3 (aicen, vicen, &
vsnon, Tsfcn, &
swvdr, swvdf, &
swidr, swidf, &
heat_capacity, &
albedo_type, &
albicev, albicei, &
albsnowv, albsnowi, &
ahmax, &
alvdrn, alidrn, &
alvdfn, alidfn, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
fswpenl, &
Iswabs, SSwabs, &
albin, albsn, &
coszen, ncat, &
nilyr)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
ncat ! number of ice thickness categories
real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen , & ! concentration of ice per category
vicen , & ! volume of ice per category
vsnon , & ! volume of ice per category
Tsfcn ! surface temperature
real (kind=dbl_kind), intent(in) :: &
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf ! sw down, near IR, diffuse (W/m^2)
! baseline albedos for ccsm3 shortwave, set in namelist
real (kind=dbl_kind), intent(in) :: &
albicev , & ! visible ice albedo for h > ahmax
albicei , & ! near-ir ice albedo for h > ahmax
albsnowv, & ! cold snow albedo, visible
albsnowi, & ! cold snow albedo, near IR
ahmax ! thickness above which ice albedo is constant (m)
logical(kind=log_kind), intent(in) :: &
heat_capacity! if true, ice has nonzero heat capacity
character (len=char_len), intent(in) :: &
albedo_type ! albedo parameterization, 'ccsm3' or 'constant'
real (kind=dbl_kind), dimension (:), intent(inout) :: &
alvdrn , & ! visible, direct, avg (fraction)
alidrn , & ! near-ir, direct, avg (fraction)
alvdfn , & ! visible, diffuse, avg (fraction)
alidfn , & ! near-ir, diffuse, avg (fraction)
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
fswint , & ! SW absorbed in ice interior, below surface (W m-2)
fswthru , & ! SW through ice to ocean (W m-2)
albin , & ! bare ice albedo
albsn ! snow albedo
real (kind=dbl_kind), dimension (:), intent(out), optional :: &
fswthru_vdr , & ! vis dir SW through ice to ocean (W m-2)
fswthru_vdf , & ! vis dif SW through ice to ocean (W m-2)
fswthru_idr , & ! nir dir SW through ice to ocean (W m-2)
fswthru_idf ! nir dif SW through ice to ocean (W m-2)
real (kind=dbl_kind), intent(inout) :: &
coszen ! cosine(zenith angle)
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
fswpenl , & ! SW entering ice layers (W m-2)
Iswabs , & ! SW absorbed in particular layer (W m-2)
Sswabs ! SW absorbed in particular layer (W m-2)
! local variables
integer (kind=int_kind) :: &
n ! thickness category index
! ice and snow albedo for each category
real (kind=dbl_kind) :: &
alvdrni, & ! visible, direct, ice (fraction)
alidrni, & ! near-ir, direct, ice (fraction)
alvdfni, & ! visible, diffuse, ice (fraction)
alidfni, & ! near-ir, diffuse, ice (fraction)
alvdrns, & ! visible, direct, snow (fraction)
alidrns, & ! near-ir, direct, snow (fraction)
alvdfns, & ! visible, diffuse, snow (fraction)
alidfns ! near-ir, diffuse, snow (fraction)
real (kind=dbl_kind), dimension(:), allocatable :: &
l_fswthru_vdr , & ! vis dir SW through ice to ocean (W m-2)
l_fswthru_vdf , & ! vis dif SW through ice to ocean (W m-2)
l_fswthru_idr , & ! nir dir SW through ice to ocean (W m-2)
l_fswthru_idf ! nir dif SW through ice to ocean (W m-2)
character(len=*),parameter :: subname='(shortwave_ccsm3)'
!-----------------------------------------------------------------
! Solar radiation: albedo and absorbed shortwave
!-----------------------------------------------------------------
allocate(l_fswthru_vdr(ncat))
allocate(l_fswthru_vdf(ncat))
allocate(l_fswthru_idr(ncat))
allocate(l_fswthru_idf(ncat))
! For basic shortwave, set coszen to a constant between 0 and 1.
coszen = p5 ! sun above the horizon
do n = 1, ncat
Sswabs(:,n) = c0
alvdrni = albocn
alidrni = albocn
alvdfni = albocn
alidfni = albocn
alvdrns = albocn
alidrns = albocn
alvdfns = albocn
alidfns = albocn
alvdrn(n) = albocn
alidrn(n) = albocn
alvdfn(n) = albocn
alidfn(n) = albocn
albin(n) = c0
albsn(n) = c0
fswsfc(n) = c0
fswint(n) = c0
fswthru(n) = c0
fswpenl(:,n) = c0
Iswabs (:,n) = c0
if (aicen(n) > puny) then
!-----------------------------------------------------------------
! Compute albedos for ice and snow.
!-----------------------------------------------------------------
if (trim(albedo_type) == 'constant') then
call constant_albedos (aicen(n), &
vsnon(n), &
Tsfcn(n), &
alvdrni, alidrni, &
alvdfni, alidfni, &
alvdrns, alidrns, &
alvdfns, alidfns, &
alvdrn(n), &
alidrn(n), &
alvdfn(n), &
alidfn(n), &
albin(n), &
albsn(n))
if (icepack_warnings_aborted(subname)) return
elseif (trim(albedo_type) == 'ccsm3') then
call compute_albedos (aicen(n), &
vicen(n), &
vsnon(n), &
Tsfcn(n), &
albicev, albicei, &
albsnowv, albsnowi, &
ahmax, &
alvdrni, alidrni, &
alvdfni, alidfni, &
alvdrns, alidrns, &
alvdfns, alidfns, &
alvdrn(n), &
alidrn(n), &
alvdfn(n), &
alidfn(n), &
albin(n), &
albsn(n))
if (icepack_warnings_aborted(subname)) return
else
call icepack_warnings_add(subname//' ERROR: albedo_type '//trim(albedo_type)//' unknown')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
!-----------------------------------------------------------------
! Compute solar radiation absorbed in ice and penetrating to ocean.
!-----------------------------------------------------------------
call absorbed_solar (heat_capacity, &
nilyr, &
aicen(n), &
vicen(n), &
vsnon(n), &
swvdr, swvdf, &
swidr, swidf, &
alvdrni, alvdfni, &
alidrni, alidfni, &
alvdrns, alvdfns, &
alidrns, alidfns, &
fswsfc=fswsfc(n), &
fswint=fswint(n), &
fswthru=fswthru(n), &
fswthru_vdr=l_fswthru_vdr(n),&
fswthru_vdf=l_fswthru_vdf(n),&
fswthru_idr=l_fswthru_idr(n),&
fswthru_idf=l_fswthru_idf(n),&
fswpenl=fswpenl(:,n), &
Iswabs=Iswabs(:,n))
if (icepack_warnings_aborted(subname)) return
endif ! aicen > puny
enddo ! ncat
if(present(fswthru_vdr)) fswthru_vdr = l_fswthru_vdr
if(present(fswthru_vdf)) fswthru_vdf = l_fswthru_vdf
if(present(fswthru_idr)) fswthru_idr = l_fswthru_idr
if(present(fswthru_idf)) fswthru_idf = l_fswthru_idf
deallocate(l_fswthru_vdr)
deallocate(l_fswthru_vdf)
deallocate(l_fswthru_idr)
deallocate(l_fswthru_idf)
end subroutine shortwave_ccsm3
!=======================================================================
!
! Compute albedos for each thickness category
subroutine compute_albedos (aicen, vicen, &
vsnon, Tsfcn, &
albicev, albicei, &
albsnowv, albsnowi, &
ahmax, &
alvdrni, alidrni, &
alvdfni, alidfni, &
alvdrns, alidrns, &
alvdfns, alidfns, &
alvdrn, alidrn, &
alvdfn, alidfn, &
albin, albsn)
real (kind=dbl_kind), intent(in) :: &
aicen , & ! concentration of ice per category
vicen , & ! volume of ice per category
vsnon , & ! volume of ice per category
Tsfcn ! surface temperature
! baseline albedos for ccsm3 shortwave, set in namelist
real (kind=dbl_kind), intent(in) :: &
albicev , & ! visible ice albedo for h > ahmax
albicei , & ! near-ir ice albedo for h > ahmax
albsnowv, & ! cold snow albedo, visible
albsnowi, & ! cold snow albedo, near IR
ahmax ! thickness above which ice albedo is constant (m)
real (kind=dbl_kind), intent(out) :: &
alvdrni , & ! visible, direct, ice (fraction)
alidrni , & ! near-ir, direct, ice (fraction)
alvdfni , & ! visible, diffuse, ice (fraction)
alidfni , & ! near-ir, diffuse, ice (fraction)
alvdrns , & ! visible, direct, snow (fraction)
alidrns , & ! near-ir, direct, snow (fraction)
alvdfns , & ! visible, diffuse, snow (fraction)
alidfns , & ! near-ir, diffuse, snow (fraction)
alvdrn , & ! visible, direct, avg (fraction)
alidrn , & ! near-ir, direct, avg (fraction)
alvdfn , & ! visible, diffuse, avg (fraction)
alidfn , & ! near-ir, diffuse, avg (fraction)
albin , & ! bare ice
albsn ! snow
! local variables
real (kind=dbl_kind), parameter :: &
dT_melt = c1 , & ! change in temp to give dalb_mlt
! albedo change
dalb_mlt = -0.075_dbl_kind, & ! albedo change per dT_melt change
! in temp for ice
dalb_mltv = -p1 , & ! albedo vis change per dT_melt change
! in temp for snow
dalb_mlti = -p15 ! albedo nir change per dT_melt change
! in temp for snow
real (kind=dbl_kind) :: &
hi , & ! ice thickness (m)
hs , & ! snow thickness (m)
albo, & ! effective ocean albedo, function of ice thickness
fh , & ! piecewise linear function of thickness
fT , & ! piecewise linear function of surface temperature
dTs , & ! difference of Tsfc and Timelt
fhtan,& ! factor used in albedo dependence on ice thickness
asnow ! fractional area of snow cover
character(len=*),parameter :: subname='(compute_albedos)'
fhtan = atan(ahmax*c4)
!-----------------------------------------------------------------
! Compute albedo for each thickness category.
!-----------------------------------------------------------------
hi = vicen / aicen
hs = vsnon / aicen
! bare ice, thickness dependence
fh = min(atan(hi*c4)/fhtan,c1)
albo = albocn*(c1-fh)
alvdfni = albicev*fh + albo
alidfni = albicei*fh + albo
! bare ice, temperature dependence
dTs = Timelt - Tsfcn
fT = min(dTs/dT_melt-c1,c0)
alvdfni = alvdfni - dalb_mlt*fT
alidfni = alidfni - dalb_mlt*fT
! avoid negative albedos for thin, bare, melting ice
alvdfni = max (alvdfni, albocn)
alidfni = max (alidfni, albocn)
if (hs > puny) then
alvdfns = albsnowv
alidfns = albsnowi
! snow on ice, temperature dependence
alvdfns = alvdfns - dalb_mltv*fT
alidfns = alidfns - dalb_mlti*fT
endif ! hs > puny
! direct albedos (same as diffuse for now)
alvdrni = alvdfni
alidrni = alidfni
alvdrns = alvdfns
alidrns = alidfns
! fractional area of snow cover
if (hs > puny) then
asnow = hs / (hs + snowpatch)
else
asnow = c0
endif
! combine ice and snow albedos (for coupler)
alvdfn = alvdfni*(c1-asnow) + &
alvdfns*asnow
alidfn = alidfni*(c1-asnow) + &
alidfns*asnow
alvdrn = alvdrni*(c1-asnow) + &
alvdrns*asnow
alidrn = alidrni*(c1-asnow) + &
alidrns*asnow
! save ice and snow albedos (for history)
albin = awtvdr*alvdrni + awtidr*alidrni &
+ awtvdf*alvdfni + awtidf*alidfni
albsn = awtvdr*alvdrns + awtidr*alidrns &
+ awtvdf*alvdfns + awtidf*alidfns
end subroutine compute_albedos
!=======================================================================
!
! Compute albedos for each thickness category
subroutine constant_albedos (aicen, &
vsnon, Tsfcn, &
alvdrni, alidrni, &
alvdfni, alidfni, &
alvdrns, alidrns, &
alvdfns, alidfns, &
alvdrn, alidrn, &
alvdfn, alidfn, &
albin, albsn)
real (kind=dbl_kind), intent(in) :: &
aicen , & ! concentration of ice per category
vsnon , & ! volume of ice per category
Tsfcn ! surface temperature
real (kind=dbl_kind), intent(out) :: &
alvdrni , & ! visible, direct, ice (fraction)
alidrni , & ! near-ir, direct, ice (fraction)
alvdfni , & ! visible, diffuse, ice (fraction)
alidfni , & ! near-ir, diffuse, ice (fraction)
alvdrns , & ! visible, direct, snow (fraction)
alidrns , & ! near-ir, direct, snow (fraction)
alvdfns , & ! visible, diffuse, snow (fraction)
alidfns , & ! near-ir, diffuse, snow (fraction)
alvdrn , & ! visible, direct, avg (fraction)
alidrn , & ! near-ir, direct, avg (fraction)
alvdfn , & ! visible, diffuse, avg (fraction)
alidfn , & ! near-ir, diffuse, avg (fraction)
albin , & ! bare ice
albsn ! snow
! local variables
real (kind=dbl_kind), parameter :: &
warmice = 0.68_dbl_kind, &
coldice = 0.70_dbl_kind, &
warmsnow = 0.77_dbl_kind, &
coldsnow = 0.81_dbl_kind
real (kind=dbl_kind) :: &
hs ! snow thickness (m)
character(len=*),parameter :: subname='(constant_albedos)'
!-----------------------------------------------------------------
! Compute albedo for each thickness category.
!-----------------------------------------------------------------
hs = vsnon / aicen
if (hs > puny) then
! snow, temperature dependence
if (Tsfcn >= -c2*puny) then
alvdfn = warmsnow
alidfn = warmsnow
else
alvdfn = coldsnow
alidfn = coldsnow
endif
else ! hs < puny
! bare ice, temperature dependence
if (Tsfcn >= -c2*puny) then
alvdfn = warmice
alidfn = warmice
else
alvdfn = coldice
alidfn = coldice
endif
endif ! hs > puny
! direct albedos (same as diffuse for now)
alvdrn = alvdfn
alidrn = alidfn
alvdrni = alvdrn
alidrni = alidrn
alvdrns = alvdrn
alidrns = alidrn
alvdfni = alvdfn
alidfni = alidfn
alvdfns = alvdfn
alidfns = alidfn
! save ice and snow albedos (for history)
albin = awtvdr*alvdrni + awtidr*alidrni &
+ awtvdf*alvdfni + awtidf*alidfni
albsn = awtvdr*alvdrns + awtidr*alidrns &
+ awtvdf*alvdfns + awtidf*alidfns
end subroutine constant_albedos
!=======================================================================
!
! Compute solar radiation absorbed in ice and penetrating to ocean
!
! authors William H. Lipscomb, LANL
! C. M. Bitz, UW
subroutine absorbed_solar (heat_capacity, &
nilyr, aicen, &
vicen, vsnon, &
swvdr, swvdf, &
swidr, swidf, &
alvdrni, alvdfni, &
alidrni, alidfni, &
alvdrns, alvdfns, &
alidrns, alidfns, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
fswpenl, &
Iswabs)
logical(kind=log_kind), intent(in) :: &
heat_capacity ! if true, ice has nonzero heat capacity
integer (kind=int_kind), intent(in) :: &
nilyr ! number of ice layers
real (kind=dbl_kind), intent(in) :: &
aicen , & ! fractional ice area
vicen , & ! ice volume
vsnon , & ! snow volume
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf , & ! sw down, near IR, diffuse (W/m^2)
alvdrni , & ! visible, direct albedo,ice
alidrni , & ! near-ir, direct albedo,ice
alvdfni , & ! visible, diffuse albedo,ice
alidfni , & ! near-ir, diffuse albedo,ice
alvdrns , & ! visible, direct albedo, snow
alidrns , & ! near-ir, direct albedo, snow
alvdfns , & ! visible, diffuse albedo, snow
alidfns ! near-ir, diffuse albedo, snow
real (kind=dbl_kind), intent(out):: &
fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
fswint , & ! SW absorbed in ice interior, below surface (W m-2)
fswthru ! SW through ice to ocean (W m-2)
real (kind=dbl_kind), intent(out) :: &
fswthru_vdr , & ! vis dir SW through ice to ocean (W m-2)
fswthru_vdf , & ! vis dif SW through ice to ocean (W m-2)
fswthru_idr , & ! nir dir SW through ice to ocean (W m-2)
fswthru_idf ! nir dif SW through ice to ocean (W m-2)
real (kind=dbl_kind), dimension (:), intent(out) :: &
Iswabs , & ! SW absorbed in particular layer (W m-2)
fswpenl ! visible SW entering ice layers (W m-2)
! local variables
real (kind=dbl_kind), parameter :: &
i0vis = 0.70_dbl_kind ! fraction of penetrating solar rad (visible)
integer (kind=int_kind) :: &
k ! ice layer index
real (kind=dbl_kind) :: &
fswpen , & ! SW penetrating beneath surface (W m-2)
trantop , & ! transmitted frac of penetrating SW at layer top
tranbot ! transmitted frac of penetrating SW at layer bot
real (kind=dbl_kind) :: &
swabs , & ! net SW down at surface (W m-2)
swabsv , & ! swabs in vis (wvlngth < 700nm) (W/m^2)
swabsi , & ! swabs in nir (wvlngth > 700nm) (W/m^2)
fswpenvdr , & ! penetrating SW, vis direct
fswpenvdf , & ! penetrating SW, vis diffuse
hi , & ! ice thickness (m)
hs , & ! snow thickness (m)
hilyr , & ! ice layer thickness
asnow ! fractional area of snow cover
character(len=*),parameter :: subname='(absorbed_solar)'
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
trantop = c0
tranbot = c0
hs = vsnon / aicen
!-----------------------------------------------------------------
! Fractional snow cover
!-----------------------------------------------------------------
if (hs > puny) then
asnow = hs / (hs + snowpatch)
else
asnow = c0
endif
!-----------------------------------------------------------------
! Shortwave flux absorbed at surface, absorbed internally,
! and penetrating to mixed layer.
! This parameterization assumes that all IR is absorbed at the
! surface; only visible is absorbed in the ice interior or
! transmitted to the ocean.
!-----------------------------------------------------------------
swabsv = swvdr * ( (c1-alvdrni)*(c1-asnow) &
+ (c1-alvdrns)*asnow ) &
+ swvdf * ( (c1-alvdfni)*(c1-asnow) &
+ (c1-alvdfns)*asnow )
swabsi = swidr * ( (c1-alidrni)*(c1-asnow) &
+ (c1-alidrns)*asnow ) &
+ swidf * ( (c1-alidfni)*(c1-asnow) &
+ (c1-alidfns)*asnow )
swabs = swabsv + swabsi
fswpenvdr = swvdr * (c1-alvdrni) * (c1-asnow) * i0vis
fswpenvdf = swvdf * (c1-alvdfni) * (c1-asnow) * i0vis
! no penetrating radiation in near IR
! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir
! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir
fswpen = fswpenvdr + fswpenvdf
fswsfc = swabs - fswpen
trantop = c1 ! transmittance at top of ice
!-----------------------------------------------------------------
! penetrating SW absorbed in each ice layer
!-----------------------------------------------------------------
do k = 1, nilyr
hi = vicen / aicen
hilyr = hi / real(nilyr,kind=dbl_kind)
tranbot = exp (-kappav * hilyr * real(k,kind=dbl_kind))
Iswabs(k) = fswpen * (trantop-tranbot)
! bottom of layer k = top of layer k+1
trantop = tranbot
! bgc layer model
if (k == 1) then ! surface flux
fswpenl(k) = fswpen
fswpenl(k+1) = fswpen * tranbot
else
fswpenl(k+1) = fswpen * tranbot
endif
enddo ! nilyr
! SW penetrating thru ice into ocean
fswthru = fswpen * tranbot
fswthru_vdr = fswpenvdr * tranbot
fswthru_vdf = fswpenvdf * tranbot
fswthru_idr = c0
fswthru_idf = c0
! SW absorbed in ice interior
fswint = fswpen - fswthru
!----------------------------------------------------------------
! if zero-layer model (no heat capacity), no SW is absorbed in ice
! interior, so add to surface absorption
!----------------------------------------------------------------
if (.not. heat_capacity) then
! SW absorbed at snow/ice surface
fswsfc = fswsfc + fswint
! SW absorbed in ice interior (nilyr = 1)
fswint = c0
Iswabs(1) = c0
endif ! heat_capacity
end subroutine absorbed_solar
! End ccsm3 shortwave method
!=======================================================================
! Begin Delta-Eddington shortwave method
! Compute initial data for Delta-Eddington method, specifically,
! the approximate exponential look-up table.
!
! author: Bruce P. Briegleb, NCAR
! 2011 ECH modified for melt pond tracers
! 2013 ECH merged with NCAR version
subroutine run_dEdd(dt, ncat, &
dEdd_algae, &
nilyr, nslyr, &
aicen, vicen, &
vsnon, Tsfcn, &
alvln, apndn, &
hpndn, ipndn, &
aeron, kalg, &
trcrn_bgcsw, &
heat_capacity, &
tlat, tlon, &
calendar_type, &
days_per_year, &
nextsw_cday, yday, &
sec, R_ice, &
R_pnd, R_snw, &
dT_mlt, rsnw_mlt, &
hs0, hs1, hp1, &
pndaspect, &
kaer_tab, waer_tab, &
gaer_tab, &
kaer_bc_tab, &
waer_bc_tab, &
gaer_bc_tab, &
bcenh, &
modal_aero, &
swvdr, swvdf, &
swidr, swidf, &
coszen, fsnow, &
alvdrn, alvdfn, &
alidrn, alidfn, &
fswsfcn, fswintn, &
fswthrun, &
fswthrun_vdr, &
fswthrun_vdf, &
fswthrun_idr, &
fswthrun_idf, &
fswpenln, &
Sswabsn, Iswabsn, &
albicen, albsnon, &
albpndn, apeffn, &
snowfracn, &
dhsn, ffracn, &
l_print_point, &
initonly)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of ice thickness categories
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical(kind=log_kind), intent(in) :: &
heat_capacity,& ! if true, ice has nonzero heat capacity
dEdd_algae, & ! .true. use prognostic chla in dEdd
modal_aero ! .true. use modal aerosol treatment
! dEdd tuning parameters, set in namelist
real (kind=dbl_kind), intent(in) :: &
R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
rsnw_mlt, & ! maximum melting snow grain radius (10^-6 m)
hs0 , & ! snow depth for transition to bare sea ice (m)
pndaspect, & ! ratio of pond depth to pond fraction
hs1 , & ! tapering parameter for snow on pond ice
hp1 , & ! critical parameter for pond ice thickness
kalg ! algae absorption coefficient
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_tab, & ! aerosol single scatter albedo (fraction)
gaer_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_bc_tab, & ! aerosol single scatter albedo (fraction)
gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
bcenh ! BC absorption enhancement factor
character (len=char_len), intent(in) :: &
calendar_type ! differentiates Gregorian from other calendars
integer (kind=int_kind), intent(in) :: &
days_per_year, & ! number of days in one year
sec ! elapsed seconds into date
real (kind=dbl_kind), intent(in) :: &
nextsw_cday , & ! julian day of next shortwave calculation
yday ! day of the year
real(kind=dbl_kind), intent(in) :: &
dt, & ! time step (s)
tlat, & ! latitude of temp pts (radians)
tlon, & ! longitude of temp pts (radians)
swvdr, & ! sw down, visible, direct (W/m^2)
swvdf, & ! sw down, visible, diffuse (W/m^2)
swidr, & ! sw down, near IR, direct (W/m^2)
swidf, & ! sw down, near IR, diffuse (W/m^2)
fsnow ! snowfall rate (kg/m^2 s)
real(kind=dbl_kind), dimension(:), intent(in) :: &
aicen, & ! concentration of ice
vicen, & ! volume per unit area of ice (m)
vsnon, & ! volume per unit area of snow (m)
Tsfcn, & ! surface temperature (deg C)
alvln, & ! level-ice area fraction
apndn, & ! pond area fraction
hpndn, & ! pond depth (m)
ipndn ! pond refrozen lid thickness (m)
real(kind=dbl_kind), dimension(:,:), intent(in) :: &
aeron, & ! aerosols (kg/m^3)
trcrn_bgcsw ! zaerosols (kg/m^3) + chlorophyll on shorthwave grid
real(kind=dbl_kind), dimension(:), intent(inout) :: &
ffracn,& ! fraction of fsurfn used to melt ipond
dhsn ! depth difference for snow on sea ice and pond ice
real(kind=dbl_kind), intent(inout) :: &
coszen ! cosine solar zenith angle, < 0 for sun below horizon
real(kind=dbl_kind), dimension(:), intent(inout) :: &
alvdrn, & ! visible direct albedo (fraction)
alvdfn, & ! near-ir direct albedo (fraction)
alidrn, & ! visible diffuse albedo (fraction)
alidfn, & ! near-ir diffuse albedo (fraction)
fswsfcn, & ! SW absorbed at ice/snow surface (W m-2)
fswintn, & ! SW absorbed in ice interior, below surface (W m-2)
fswthrun, & ! SW through ice to ocean (W/m^2)
albicen, & ! albedo bare ice
albsnon, & ! albedo snow
albpndn, & ! albedo pond
apeffn, & ! effective pond area used for radiation calculation
snowfracn ! snow fraction on each category used for radiation
real(kind=dbl_kind), dimension(:), intent(out), optional :: &
fswthrun_vdr, & ! vis dir SW through ice to ocean (W/m^2)
fswthrun_vdf, & ! vis dif SW through ice to ocean (W/m^2)
fswthrun_idr, & ! nir dir SW through ice to ocean (W/m^2)
fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
Sswabsn , & ! SW radiation absorbed in snow layers (W m-2)
Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
fswpenln ! visible SW entering ice layers (W m-2)
logical (kind=log_kind), intent(in) :: &
l_print_point
logical (kind=log_kind), optional :: &
initonly ! flag to indicate init only, default is false
! local temporary variables
! other local variables
! snow variables for Delta-Eddington shortwave
real (kind=dbl_kind) :: &
fsn , & ! snow horizontal fraction
hsn ! snow depth (m)
real (kind=dbl_kind), dimension (nslyr) :: &
rhosnwn , & ! snow density (kg/m3)
rsnwn ! snow grain radius (micrometers)
! pond variables for Delta-Eddington shortwave
real (kind=dbl_kind) :: &
fpn , & ! pond fraction of ice cover
hpn ! actual pond depth (m)
integer (kind=int_kind) :: &
n ! thickness category index
real (kind=dbl_kind) :: &
ipn , & ! refrozen pond ice thickness (m), mean over ice fraction
hp , & ! pond depth
hs , & ! snow depth
asnow , & ! fractional area of snow cover
rp , & ! volume fraction of retained melt water to total liquid content
hmx , & ! maximum available snow infiltration equivalent depth
dhs , & ! local difference in snow depth on sea ice and pond ice
spn , & ! snow depth on refrozen pond (m)
tmp ! 0 or 1
logical (kind=log_kind) :: &
linitonly ! local initonly value
real (kind=dbl_kind), dimension(:), allocatable :: &
l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W m-2)
l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W m-2)
l_fswthrun_idr , & ! nir dir SW through ice to ocean (W m-2)
l_fswthrun_idf ! nir dif SW through ice to ocean (W m-2)
character(len=*),parameter :: subname='(run_dEdd)'
allocate(l_fswthrun_vdr(ncat))
allocate(l_fswthrun_vdf(ncat))
allocate(l_fswthrun_idr(ncat))
allocate(l_fswthrun_idf(ncat))
linitonly = .false.
if (present(initonly)) then
linitonly = initonly
endif
! cosine of the zenith angle
#ifdef CESMCOUPLED
call compute_coszen (tlat, tlon, &
yday, sec, coszen, &
days_per_year, nextsw_cday, calendar_type)
#else
call compute_coszen (tlat, tlon, &
yday, sec, coszen)
#endif
if (icepack_warnings_aborted(subname)) return
do n = 1, ncat
! note that rhoswn, rsnw, fp, hp and Sswabs ARE NOT dimensioned with ncat
! BPB 19 Dec 2006
! set snow properties
fsn = c0
hsn = c0
rhosnwn(:) = c0
rsnwn(:) = c0
apeffn(n) = c0 ! for history
snowfracn(n) = c0 ! for history
if (aicen(n) > puny) then
call shortwave_dEdd_set_snow(nslyr, R_snw, &
dT_mlt, rsnw_mlt, &
aicen(n), vsnon(n), &
Tsfcn(n), fsn, &
hs0, hsn, &
rhosnwn, rsnwn)
if (icepack_warnings_aborted(subname)) return
! set pond properties
if (tr_pond_cesm) then
! fraction of ice area
fpn = apndn(n)
! pond depth over fraction fpn
hpn = hpndn(n)
! snow infiltration
if (hsn >= hs_min .and. hs0 > puny) then
asnow = min(hsn/hs0, c1) ! delta-Eddington formulation
fpn = (c1 - asnow) * fpn
hpn = pndaspect * fpn
endif
! Zero out fraction of thin ponds for radiation only
if (hpn < hpmin) fpn = c0
fsn = min(fsn, c1-fpn)
apeffn(n) = fpn ! for history
elseif (tr_pond_lvl) then
fpn = c0 ! fraction of ice covered in pond
hpn = c0 ! pond depth over fpn
! refrozen pond lid thickness avg over ice
! allow snow to cover pond ice
ipn = alvln(n) * apndn(n) * ipndn(n)
dhs = dhsn(n) ! snow depth difference, sea ice - pond
if (.not. linitonly .and. ipn > puny .and. &
dhs < puny .and. fsnow*dt > hs_min) &
dhs = hsn - fsnow*dt ! initialize dhs>0
spn = hsn - dhs ! snow depth on pond ice
if (.not. linitonly .and. ipn*spn < puny) dhs = c0
dhsn(n) = dhs ! save: constant until reset to 0
! not using ipn assumes that lid ice is perfectly clear
! if (ipn <= 0.3_dbl_kind) then
! fraction of ice area
fpn = apndn(n) * alvln(n)
! pond depth over fraction fpn
hpn = hpndn(n)
! reduce effective pond area absorbing surface heat flux
! due to flux already having been used to melt pond ice
fpn = (c1 - ffracn(n)) * fpn
! taper pond area with snow on pond ice
if (dhs > puny .and. spn >= puny .and. hs1 > puny) then
asnow = min(spn/hs1, c1)
fpn = (c1 - asnow) * fpn
endif
! infiltrate snow
hp = hpn
if (hp > puny) then
hs = hsn
rp = rhofresh*hp/(rhofresh*hp + rhos*hs)
if (rp < p15) then
fpn = c0
hpn = c0
else
hmx = hs*(rhofresh - rhos)/rhofresh
tmp = max(c0, sign(c1, hp-hmx)) ! 1 if hp>=hmx, else 0
hp = (rhofresh*hp + rhos*hs*tmp) &
/ (rhofresh - rhos*(c1-tmp))
hsn = hs - hp*fpn*(c1-tmp)
hpn = hp * tmp
fpn = fpn * tmp
endif
endif ! hp > puny
! Zero out fraction of thin ponds for radiation only
if (hpn < hpmin) fpn = c0
fsn = min(fsn, c1-fpn)
! endif ! masking by lid ice
apeffn(n) = fpn ! for history
elseif (tr_pond_topo) then
! Lid effective if thicker than hp1
if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
fpn = apndn(n)
else
fpn = c0
endif
if (apndn(n) > puny) then
hpn = hpndn(n)
else
fpn = c0
hpn = c0
endif
! Zero out fraction of thin ponds for radiation only
if (hpn < hpmin) fpn = c0
! If ponds are present snow fraction reduced to
! non-ponded part dEdd scheme
fsn = min(fsn, c1-fpn)
apeffn(n) = fpn
else
fpn = c0
hpn = c0
call shortwave_dEdd_set_pond(Tsfcn(n), &
fsn, fpn, &
hpn)
if (icepack_warnings_aborted(subname)) return
apeffn(n) = fpn ! for history
fpn = c0
hpn = c0
endif ! pond type
snowfracn(n) = fsn ! for history
call shortwave_dEdd(dEdd_algae, &
nslyr, nilyr, &
coszen, heat_capacity, &
aicen(n), vicen(n), &
hsn, fsn, &
rhosnwn, rsnwn, &
fpn, hpn, &
aeron(:,n), &
R_ice, R_pnd, &
kaer_tab, waer_tab, &
gaer_tab, &
kaer_bc_tab, &
waer_bc_tab, &
gaer_bc_tab, &
bcenh, modal_aero, &
kalg, &
swvdr, swvdf, &
swidr, swidf, &
alvdrn(n), alvdfn(n), &
alidrn(n), alidfn(n), &
fswsfcn(n), fswintn(n), &
fswthru=fswthrun(n), &
fswthru_vdr=l_fswthrun_vdr(n), &
fswthru_vdf=l_fswthrun_vdf(n), &
fswthru_idr=l_fswthrun_idr(n), &
fswthru_idf=l_fswthrun_idf(n), &
Sswabs=Sswabsn(:,n), &
Iswabs=Iswabsn(:,n), &
albice=albicen(n), &
albsno=albsnon(n), &
albpnd=albpndn(n), &
fswpenl=fswpenln(:,n), &
zbio=trcrn_bgcsw(:,n), &
l_print_point=l_print_point)
if (icepack_warnings_aborted(subname)) return
endif ! aicen > puny
enddo ! ncat
if(present(fswthrun_vdr)) fswthrun_vdr = l_fswthrun_vdr
if(present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf
if(present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr
if(present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf
deallocate(l_fswthrun_vdr)
deallocate(l_fswthrun_vdf)
deallocate(l_fswthrun_idr)
deallocate(l_fswthrun_idf)
end subroutine run_dEdd
!=======================================================================
!
! Compute snow/bare ice/ponded ice shortwave albedos, absorbed and transmitted
! flux using the Delta-Eddington solar radiation method as described in:
!
! A Delta-Eddington Multiple Scattering Parameterization for Solar Radiation
! in the Sea Ice Component of the Community Climate System Model
! B.P.Briegleb and B.Light NCAR/TN-472+STR February 2007
!
! Compute shortwave albedos and fluxes for three surface types:
! snow over ice, bare ice and ponded ice.
!
! Albedos and fluxes are output for later use by thermodynamic routines.
! Invokes three calls to compute_dEdd, which sets inherent optical properties
! appropriate for the surface type. Within compute_dEdd, a call to solution_dEdd
! evaluates the Delta-Eddington solution. The final albedos and fluxes are then
! evaluated in compute_dEdd. Albedos and fluxes are transferred to output in
! this routine.
!
! NOTE regarding albedo diagnostics: This method yields zero albedo values
! if there is no incoming solar and thus the albedo diagnostics are masked
! out when the sun is below the horizon. To estimate albedo from the history
! output (post-processing), compute ice albedo using
! (1 - albedo)*swdn = swabs. -ECH
!
! author: Bruce P. Briegleb, NCAR
! 2013: E Hunke merged with NCAR version
!
subroutine shortwave_dEdd (dEdd_algae, &
nslyr, nilyr, &
coszen, heat_capacity,&
aice, vice, &
hs, fs, &
rhosnw, rsnw, &
fp, hp, &
aero, &
R_ice, R_pnd, &
kaer_tab, waer_tab, &
gaer_tab, &
kaer_bc_tab, &
waer_bc_tab, &
gaer_bc_tab, &
bcenh, modal_aero, &
kalg, &
swvdr, swvdf, &
swidr, swidf, &
alvdr, alvdf, &
alidr, alidf, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
Sswabs, &
Iswabs, albice, &
albsno, albpnd, &
fswpenl, zbio, &
l_print_point)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr ! number of snow layers
logical (kind=log_kind), intent(in) :: &
heat_capacity, & ! if true, ice has nonzero heat capacity
dEdd_algae, & ! .true. use prognostic chla in dEdd
modal_aero ! .true. use modal aerosol treatment
real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_bc_tab, & ! aerosol single scatter albedo (fraction)
gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
bcenh ! BC absorption enhancement factor
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_tab, & ! aerosol single scatter albedo (fraction)
gaer_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), intent(in) :: &
kalg , & ! algae absorption coefficient
R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo
aice , & ! concentration of ice
vice , & ! volume of ice
hs , & ! snow depth
fs ! horizontal coverage of snow
real (kind=dbl_kind), dimension (:), intent(in) :: &
rhosnw , & ! density in snow layer (kg/m3)
rsnw , & ! grain radius in snow layer (m)
aero , & ! aerosol tracers
zbio ! shortwave tracers (zaero+chla)
real (kind=dbl_kind), intent(in) :: &
fp , & ! pond fractional coverage (0 to 1)
hp , & ! pond depth (m)
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf ! sw down, near IR, diffuse (W/m^2)
real (kind=dbl_kind), intent(inout) :: &
coszen , & ! cosine of solar zenith angle
alvdr , & ! visible, direct, albedo (fraction)
alvdf , & ! visible, diffuse, albedo (fraction)
alidr , & ! near-ir, direct, albedo (fraction)
alidf , & ! near-ir, diffuse, albedo (fraction)
fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
fswint , & ! SW interior absorption (below surface, above ocean,W m-2)
fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2)
real (kind=dbl_kind), intent(out) :: &
fswthru_vdr , & ! vis dir SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_vdf , & ! vis dif SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_idr , & ! nir dir SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_idf ! nir dif SW through snow/bare ice/ponded ice into ocean (W m-2)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
fswpenl , & ! visible SW entering ice layers (W m-2)
Sswabs , & ! SW absorbed in snow layer (W m-2)
Iswabs ! SW absorbed in ice layer (W m-2)
real (kind=dbl_kind), intent(out) :: &
albice , & ! bare ice albedo, for history
albsno , & ! snow albedo, for history
albpnd ! pond albedo, for history
logical (kind=log_kind) , intent(in) :: &
l_print_point
! local variables
real (kind=dbl_kind) :: &
netsw , & ! net shortwave
fnidr , & ! fraction of direct to total down surface flux in nir
hstmp , & ! snow thickness (set to 0 for bare ice case)
hi , & ! ice thickness (all sea ice layers, m)
fi ! snow/bare ice fractional coverage (0 to 1)
real (kind=dbl_kind), dimension (4*n_aero) :: &
aero_mp ! aerosol mass path in kg/m2
integer (kind=int_kind) :: &
srftyp ! surface type over ice: (0=air, 1=snow, 2=pond)
integer (kind=int_kind) :: &
k , & ! level index
na , & ! aerosol index
klev , & ! number of radiation layers - 1
klevp ! number of radiation interfaces - 1
! (0 layer is included also)
real (kind=dbl_kind) :: &
vsno ! volume of snow
real (kind=dbl_kind) :: &
swdn , & ! swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j)
swab , & ! fswsfc(i,j)+fswint(i,j)+fswthru(i,j)
swalb ! (1.-swab/(swdn+.0001))
! for history
real (kind=dbl_kind) :: &
avdrl , & ! visible, direct, albedo (fraction)
avdfl , & ! visible, diffuse, albedo (fraction)
aidrl , & ! near-ir, direct, albedo (fraction)
aidfl ! near-ir, diffuse, albedo (fraction)
character(len=*),parameter :: subname='(shortwave_dEdd)'
!-----------------------------------------------------------------------
klev = nslyr + nilyr + 1 ! number of radiation layers - 1
klevp = klev + 1 ! number of radiation interfaces - 1
! (0 layer is included also)
! zero storage albedos and fluxes for accumulation over surface types:
hstmp = c0
hi = c0
fi = c0
alvdr = c0
alvdf = c0
alidr = c0
alidf = c0
avdrl = c0
avdfl = c0
aidrl = c0
aidfl = c0
fswsfc = c0
fswint = c0
fswthru = c0
fswthru_vdr = c0
fswthru_vdf = c0
fswthru_idr = c0
fswthru_idf = c0
! compute fraction of nir down direct to total over all points:
fnidr = c0
if( swidr + swidf > puny ) then
fnidr = swidr/(swidr+swidf)
endif
albice = c0
albsno = c0
albpnd = c0
fswpenl(:) = c0
Sswabs(:) = c0
Iswabs(:) = c0
! compute aerosol mass path
aero_mp(:) = c0
if( tr_aero ) then
! check 4 layers for each aerosol, a snow SSL, snow below SSL,
! sea ice SSL, and sea ice below SSL, in that order.
if (size(aero) < 4*n_aero) then
call icepack_warnings_add(subname//' ERROR: size(aero) too small')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
do na = 1, 4*n_aero, 4
vsno = hs * aice
netsw = swvdr + swidr + swvdf + swidf
if (netsw > puny) then ! sun above horizon
aero_mp(na ) = aero(na )*vsno
aero_mp(na+1) = aero(na+1)*vsno
aero_mp(na+2) = aero(na+2)*vice
aero_mp(na+3) = aero(na+3)*vice
endif ! aice > 0 and netsw > 0
enddo ! na
endif ! if aerosols
! compute shortwave radiation accounting for snow/ice (both snow over
! ice and bare ice) and ponded ice (if any):
! sea ice points with sun above horizon
netsw = swvdr + swidr + swvdf + swidf
if (netsw > puny) then ! sun above horizon
coszen = max(puny,coszen)
! evaluate sea ice thickness and fraction
hi = vice / aice
fi = c1 - fs - fp
! bare sea ice points
if(fi > c0) then
! calculate bare sea ice
srftyp = 0
call compute_dEdd(nilyr, nslyr, klev, klevp, &
zbio, dEdd_algae, &
heat_capacity, fnidr, coszen, &
R_ice, R_pnd, &
kaer_tab, waer_tab, gaer_tab, &
kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
bcenh, modal_aero, kalg, &
swvdr, swvdf, swidr, swidf, srftyp, &
hstmp, rhosnw, rsnw, hi, hp, &
fi, aero_mp, avdrl, avdfl, &
aidrl, aidfl, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
Sswabs, &
Iswabs, fswpenl)
if (icepack_warnings_aborted(subname)) return
alvdr = alvdr + avdrl *fi
alvdf = alvdf + avdfl *fi
alidr = alidr + aidrl *fi
alidf = alidf + aidfl *fi
! for history
albice = albice &
+ awtvdr*avdrl + awtidr*aidrl &
+ awtvdf*avdfl + awtidf*aidfl
endif
endif
! sea ice points with sun above horizon
netsw = swvdr + swidr + swvdf + swidf
if (netsw > puny) then ! sun above horizon
coszen = max(puny,coszen)
! snow-covered sea ice points
if(fs > c0) then
! calculate snow covered sea ice
srftyp = 1
call compute_dEdd(nilyr, nslyr, klev, klevp, &
zbio, dEdd_algae, &
heat_capacity, fnidr, coszen, &
R_ice, R_pnd, &
kaer_tab, waer_tab, gaer_tab, &
kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
bcenh, modal_aero, kalg, &
swvdr, swvdf, swidr, swidf, srftyp, &
hs, rhosnw, rsnw, hi, hp, &
fs, aero_mp, avdrl, avdfl, &
aidrl, aidfl, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
Sswabs, &
Iswabs, fswpenl)
if (icepack_warnings_aborted(subname)) return
alvdr = alvdr + avdrl *fs
alvdf = alvdf + avdfl *fs
alidr = alidr + aidrl *fs
alidf = alidf + aidfl *fs
! for history
albsno = albsno &
+ awtvdr*avdrl + awtidr*aidrl &
+ awtvdf*avdfl + awtidf*aidfl
endif
endif
hi = c0
! sea ice points with sun above horizon
netsw = swvdr + swidr + swvdf + swidf
if (netsw > puny) then ! sun above horizon
coszen = max(puny,coszen)
hi = vice / aice
! if nonzero pond fraction and sufficient pond depth
! if( fp > puny .and. hp > hpmin ) then
if (fp > puny) then
! calculate ponded ice
srftyp = 2
call compute_dEdd(nilyr, nslyr, klev, klevp, &
zbio, dEdd_algae, &
heat_capacity, fnidr, coszen, &
R_ice, R_pnd, &
kaer_tab, waer_tab, gaer_tab, &
kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
bcenh, modal_aero, kalg, &
swvdr, swvdf, swidr, swidf, srftyp, &
hs, rhosnw, rsnw, hi, hp, &
fp, aero_mp, avdrl, avdfl, &
aidrl, aidfl, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
Sswabs, &
Iswabs, fswpenl)
if (icepack_warnings_aborted(subname)) return
alvdr = alvdr + avdrl *fp
alvdf = alvdf + avdfl *fp
alidr = alidr + aidrl *fp
alidf = alidf + aidfl *fp
! for history
albpnd = albpnd &
+ awtvdr*avdrl + awtidr*aidrl &
+ awtvdf*avdfl + awtidf*aidfl
endif
endif
! if no incoming shortwave, set albedos to 1
netsw = swvdr + swidr + swvdf + swidf
if (netsw <= puny) then ! sun above horizon
alvdr = c1
alvdf = c1
alidr = c1
alidf = c1
endif
if (l_print_point .and. netsw > puny) then
write(warnstr,*) subname, ' printing point'
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' coszen = ', &
coszen
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' swvdr swvdf = ', &
swvdr,swvdf
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' swidr swidf = ', &
swidr,swidf
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' aice = ', &
aice
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' hs = ', &
hs
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' hp = ', &
hp
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' fs = ', &
fs
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' fi = ', &
fi
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' fp = ', &
fp
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' hi = ', &
hi
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' alvdr alvdf = ', &
alvdr,alvdf
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' alidr alidf = ', &
alidr,alidf
call icepack_warnings_add(warnstr)
write(warnstr,*) subname, ' fswsfc fswint fswthru = ', &
fswsfc,fswint,fswthru
call icepack_warnings_add(warnstr)
swdn = swvdr+swvdf+swidr+swidf
swab = fswsfc+fswint+fswthru
swalb = (1.-swab/(swdn+.0001))
write(warnstr,*) subname, ' swdn swab swalb = ',swdn,swab,swalb
do k = 1, nslyr
write(warnstr,*) subname, ' snow layer k = ', k, &
' rhosnw = ', &
rhosnw(k), &
' rsnw = ', &
rsnw(k)
call icepack_warnings_add(warnstr)
enddo
do k = 1, nslyr
write(warnstr,*) subname, ' snow layer k = ', k, &
' Sswabs(k) = ', Sswabs(k)
call icepack_warnings_add(warnstr)
enddo
do k = 1, nilyr
write(warnstr,*) subname, ' sea ice layer k = ', k, &
' Iswabs(k) = ', Iswabs(k)
call icepack_warnings_add(warnstr)
enddo
endif ! l_print_point .and. coszen > .01
end subroutine shortwave_dEdd
!=======================================================================
!
! Evaluate snow/ice/ponded ice inherent optical properties (IOPs), and
! then calculate the multiple scattering solution by calling solution_dEdd.
!
! author: Bruce P. Briegleb, NCAR
! 2013: E Hunke merged with NCAR version
subroutine compute_dEdd (nilyr, nslyr, klev, klevp, &
zbio, dEdd_algae, &
heat_capacity, fnidr, coszen, &
R_ice, R_pnd, &
kaer_tab, waer_tab, gaer_tab, &
kaer_bc_tab, waer_bc_tab, gaer_bc_tab, &
bcenh, modal_aero, kalg, &
swvdr, swvdf, swidr, swidf, srftyp, &
hs, rhosnw, rsnw, hi, hp, &
fi, aero_mp, alvdr, alvdf, &
alidr, alidf, &
fswsfc, fswint, &
fswthru, &
fswthru_vdr, &
fswthru_vdf, &
fswthru_idr, &
fswthru_idf, &
Sswabs, &
Iswabs, fswpenl)
integer (kind=int_kind), intent(in) :: &
nilyr , & ! number of ice layers
nslyr , & ! number of snow layers
klev , & ! number of radiation layers - 1
klevp ! number of radiation interfaces - 1
! (0 layer is included also)
logical (kind=log_kind), intent(in) :: &
heat_capacity,& ! if true, ice has nonzero heat capacity
dEdd_algae, & ! .true. use prognostic chla in dEdd
modal_aero ! .true. use modal aerosol treatment
real (kind=dbl_kind), dimension(:,:), intent(in) :: & ! Modal aerosol treatment
kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_bc_tab, & ! aerosol single scatter albedo (fraction)
gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & ! Modal aerosol treatment
bcenh ! BC absorption enhancement factor
! dEdd tuning parameters, set in namelist
real (kind=dbl_kind), intent(in) :: &
R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo
R_pnd ! ponded ice tuning parameter; +1 > 1sig increase in albedo
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_tab, & ! aerosol single scatter albedo (fraction)
gaer_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), intent(in) :: &
kalg , & ! algae absorption coefficient
fnidr , & ! fraction of direct to total down flux in nir
coszen , & ! cosine solar zenith angle
swvdr , & ! shortwave down at surface, visible, direct (W/m^2)
swvdf , & ! shortwave down at surface, visible, diffuse (W/m^2)
swidr , & ! shortwave down at surface, near IR, direct (W/m^2)
swidf ! shortwave down at surface, near IR, diffuse (W/m^2)
integer (kind=int_kind), intent(in) :: &
srftyp ! surface type over ice: (0=air, 1=snow, 2=pond)
real (kind=dbl_kind), intent(in) :: &
hs ! snow thickness (m)
real (kind=dbl_kind), dimension (:), intent(in) :: &
rhosnw , & ! snow density in snow layer (kg/m3)
rsnw , & ! snow grain radius in snow layer (m)
zbio , & ! zaerosol + chla shortwave tracers kg/m^3
aero_mp ! aerosol mass path in kg/m2
real (kind=dbl_kind), intent(in) :: &
hi , & ! ice thickness (m)
hp , & ! pond depth (m)
fi ! snow/bare ice fractional coverage (0 to 1)
real (kind=dbl_kind), intent(inout) :: &
alvdr , & ! visible, direct, albedo (fraction)
alvdf , & ! visible, diffuse, albedo (fraction)
alidr , & ! near-ir, direct, albedo (fraction)
alidf , & ! near-ir, diffuse, albedo (fraction)
fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
fswint , & ! SW interior absorption (below surface, above ocean,W m-2)
fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2)
real (kind=dbl_kind), intent(inout) :: &
fswthru_vdr , & ! vis dir SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_vdf , & ! vis dif SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_idr , & ! nir dir SW through snow/bare ice/ponded ice into ocean (W m-2)
fswthru_idf ! nir dif SW through snow/bare ice/ponded ice into ocean (W m-2)
real (kind=dbl_kind), dimension (:), intent(inout) :: &
fswpenl , & ! visible SW entering ice layers (W m-2)
Sswabs , & ! SW absorbed in snow layer (W m-2)
Iswabs ! SW absorbed in ice layer (W m-2)
!-----------------------------------------------------------------------
!
! Set up optical property profiles, based on snow, sea ice and ponded
! ice IOPs from:
!
! Briegleb, B. P., and B. Light (2007): A Delta-Eddington Multiple
! Scattering Parameterization for Solar Radiation in the Sea Ice
! Component of the Community Climate System Model, NCAR Technical
! Note NCAR/TN-472+STR February 2007
!
! Computes column Delta-Eddington radiation solution for specific
! surface type: either snow over sea ice, bare sea ice, or ponded sea ice.
!
! Divides solar spectrum into 3 intervals: 0.2-0.7, 0.7-1.19, and
! 1.19-5.0 micro-meters. The latter two are added (using an assumed
! partition of incident shortwave in the 0.7-5.0 micro-meter band between
! the 0.7-1.19 and 1.19-5.0 micro-meter band) to give the final output
! of 0.2-0.7 visible and 0.7-5.0 near-infrared albedos and fluxes.
!
! Specifies vertical layer optical properties based on input snow depth,
! density and grain radius, along with ice and pond depths, then computes
! layer by layer Delta-Eddington reflectivity, transmissivity and combines
! layers (done by calling routine solution_dEdd). Finally, surface albedos
! and internal fluxes/flux divergences are evaluated.
!
! Description of the level and layer index conventions. This is
! for the standard case of one snow layer and four sea ice layers.
!
! Please read the following; otherwise, there is 99.9% chance you
! will be confused about indices at some point in time........ :)
!
! CICE4.0 snow treatment has one snow layer above the sea ice. This
! snow layer has finite heat capacity, so that surface absorption must
! be distinguished from internal. The Delta-Eddington solar radiation
! thus adds extra surface scattering layers to both snow and sea ice.
! Note that in the following, we assume a fixed vertical layer structure
! for the radiation calculation. In other words, we always have the
! structure shown below for one snow and four sea ice layers, but for
! ponded ice the pond fills "snow" layer 1 over the sea ice, and for
! bare sea ice the top layers over sea ice are treated as transparent air.
!
! SSL = surface scattering layer for either snow or sea ice
! DL = drained layer for sea ice immediately under sea ice SSL
! INT = interior layers for sea ice below the drained layer.
!
! Notice that the radiation level starts with 0 at the top. Thus,
! the total number radiation layers is klev+1, where klev is the
! sum of nslyr, the number of CCSM snow layers, and nilyr, the
! number of CCSM sea ice layers, plus the sea ice SSL:
! klev = 1 + nslyr + nilyr
!
! For the standard case illustrated below, nslyr=1, nilyr=4,
! and klev=6, with the number of layer interfaces klevp=klev+1.
! Layer interfaces are the surfaces on which reflectivities,
! transmissivities and fluxes are evaluated.
!
! CCSM3 Sea Ice Model Delta-Eddington Solar Radiation
! Layers and Interfaces
! Layer Index Interface Index
! --------------------- --------------------- 0
! 0 \\\ snow SSL \\\
! snow layer 1 --------------------- 1
! 1 rest of snow layer
! +++++++++++++++++++++ +++++++++++++++++++++ 2
! 2 \\\ sea ice SSL \\\
! sea ice layer 1 --------------------- 3
! 3 sea ice DL
! --------------------- --------------------- 4
!
! sea ice layer 2 4 sea ice INT
!
! --------------------- --------------------- 5
!
! sea ice layer 3 5 sea ice INT
!
! --------------------- --------------------- 6
!
! sea ice layer 4 6 sea ice INT
!
! --------------------- --------------------- 7
!
! When snow lies over sea ice, the radiation absorbed in the
! snow SSL is used for surface heating, and that in the rest
! of the snow layer for its internal heating. For sea ice in
! this case, all of the radiant heat absorbed in both the
! sea ice SSL and the DL are used for sea ice layer 1 heating.
!
! When pond lies over sea ice, and for bare sea ice, all of the
! radiant heat absorbed within and above the sea ice SSL is used
! for surface heating, and that absorbed in the sea ice DL is
! used for sea ice layer 1 heating.
!
! Basically, vertical profiles of the layer extinction optical depth (tau),
! single scattering albedo (w0) and asymmetry parameter (g) are required over
! the klev+1 layers, where klev+1 = 2 + nslyr + nilyr. All of the surface type
! information and snow/ice iop properties are evaulated in this routine, so
! the tau,w0,g profiles can be passed to solution_dEdd for multiple scattering
! evaluation. Snow, bare ice and ponded ice iops are contained in data arrays
! in this routine.
!
!-----------------------------------------------------------------------
! local variables
integer (kind=int_kind) :: &
k , & ! level index
ns , & ! spectral index
nr , & ! index for grain radius tables
ki , & ! index for internal absorption
km , & ! k starting index for snow, sea ice internal absorption
kp , & ! k+1 or k+2 index for snow, sea ice internal absorption
ksrf , & ! level index for surface absorption
ksnow , & ! level index for snow density and grain size
kii ! level starting index for sea ice (nslyr+1)
integer (kind=int_kind), parameter :: &
nmbrad = 32 ! number of snow grain radii in tables
real (kind=dbl_kind) :: &
avdr , & ! visible albedo, direct (fraction)
avdf , & ! visible albedo, diffuse (fraction)
aidr , & ! near-ir albedo, direct (fraction)
aidf ! near-ir albedo, diffuse (fraction)
real (kind=dbl_kind) :: &
fsfc , & ! shortwave absorbed at snow/bare ice/ponded ice surface (W m-2)
fint , & ! shortwave absorbed in interior (W m-2)
fthru , & ! shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
fthruvdr, & ! vis dir shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
fthruvdf, & ! vis dif shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
fthruidr, & ! nir dir shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
fthruidf ! nir dif shortwave through snow/bare ice/ponded ice to ocean (W/m^2)
real (kind=dbl_kind), dimension(nslyr) :: &
Sabs ! shortwave absorbed in snow layer (W m-2)
real (kind=dbl_kind), dimension(nilyr) :: &
Iabs ! shortwave absorbed in ice layer (W m-2)
real (kind=dbl_kind), dimension(nilyr+1) :: &
fthrul ! shortwave through to ice layers (W m-2)
real (kind=dbl_kind), dimension (nspint) :: &
wghtns ! spectral weights
real (kind=dbl_kind), parameter :: &
cp67 = 0.67_dbl_kind , & ! nir band weight parameter
cp78 = 0.78_dbl_kind , & ! nir band weight parameter
cp01 = 0.01_dbl_kind ! for ocean visible albedo
real (kind=dbl_kind), dimension (0:klev) :: &
tau , & ! layer extinction optical depth
w0 , & ! layer single scattering albedo
g ! layer asymmetry parameter
! following arrays are defined at model interfaces; 0 is the top of the
! layer above the sea ice; klevp is the sea ice/ocean interface.
real (kind=dbl_kind), dimension (0:klevp) :: &
trndir , & ! solar beam down transmission from top
trntdr , & ! total transmission to direct beam for layers above
trndif , & ! diffuse transmission to diffuse beam for layers above
rupdir , & ! reflectivity to direct radiation for layers below
rupdif , & ! reflectivity to diffuse radiation for layers below
rdndif ! reflectivity to diffuse radiation for layers above
real (kind=dbl_kind), dimension (0:klevp) :: &
dfdir , & ! down-up flux at interface due to direct beam at top surface
dfdif ! down-up flux at interface due to diffuse beam at top surface
real (kind=dbl_kind) :: &
refk , & ! interface k multiple scattering term
delr , & ! snow grain radius interpolation parameter
! inherent optical properties (iop) for snow
Qs , & ! Snow extinction efficiency
ks , & ! Snow extinction coefficient (/m)
ws , & ! Snow single scattering albedo
gs ! Snow asymmetry parameter
real (kind=dbl_kind), dimension(nslyr) :: &
frsnw ! snow grain radius in snow layer * adjustment factor (m)
! actual used ice and ponded ice IOPs, allowing for tuning
! modifications of the above "_mn" value
real (kind=dbl_kind), dimension (nspint) :: &
ki_ssl , & ! Surface-scattering-layer ice extinction coefficient (/m)
wi_ssl , & ! Surface-scattering-layer ice single scattering albedo
gi_ssl , & ! Surface-scattering-layer ice asymmetry parameter
ki_dl , & ! Drained-layer ice extinction coefficient (/m)
wi_dl , & ! Drained-layer ice single scattering albedo
gi_dl , & ! Drained-layer ice asymmetry parameter
ki_int , & ! Interior-layer ice extinction coefficient (/m)
wi_int , & ! Interior-layer ice single scattering albedo
gi_int , & ! Interior-layer ice asymmetry parameter
ki_p_ssl , & ! Ice under pond srf scat layer extinction coefficient (/m)
wi_p_ssl , & ! Ice under pond srf scat layer single scattering albedo
gi_p_ssl , & ! Ice under pond srf scat layer asymmetry parameter
ki_p_int , & ! Ice under pond extinction coefficient (/m)
wi_p_int , & ! Ice under pond single scattering albedo
gi_p_int ! Ice under pond asymmetry parameter
real (kind=dbl_kind), dimension(0:klev) :: &
dzk ! layer thickness
real (kind=dbl_kind) :: &
dz , & ! snow, sea ice or pond water layer thickness
dz_ssl , & ! snow or sea ice surface scattering layer thickness
fs ! scaling factor to reduce (nilyr<4) or increase (nilyr>4) DL
! extinction coefficient to maintain DL optical depth constant
! with changing number of sea ice layers, to approximately
! conserve computed albedo for constant physical depth of sea
! ice when the number of sea ice layers vary
real (kind=dbl_kind) :: &
sig , & ! scattering coefficient for tuning
kabs , & ! absorption coefficient for tuning
sigp ! modified scattering coefficient for tuning
real (kind=dbl_kind), dimension(nspint, 0:klev) :: &
kabs_chl , & ! absorption coefficient for chlorophyll (/m)
tzaer , & ! total aerosol extinction optical depth
wzaer , & ! total aerosol single scatter albedo
gzaer ! total aerosol asymmetry parameter
real (kind=dbl_kind) :: &
albodr , & ! spectral ocean albedo to direct rad
albodf ! spectral ocean albedo to diffuse rad
! for melt pond transition to bare sea ice for small pond depths
real (kind=dbl_kind) :: &
sig_i , & ! ice scattering coefficient (/m)
sig_p , & ! pond scattering coefficient (/m)
kext ! weighted extinction coefficient (/m)
! aerosol optical properties from Mark Flanner, 26 June 2008
! order assumed: hydrophobic black carbon, hydrophilic black carbon,
! four dust aerosols by particle size range:
! dust1(.05-0.5 micron), dust2(0.5-1.25 micron),
! dust3(1.25-2.5 micron), dust4(2.5-5.0 micron)
! spectral bands same as snow/sea ice: (0.3-0.7 micron, 0.7-1.19 micron
! and 1.19-5.0 micron in wavelength)
integer (kind=int_kind) :: &
na , n ! aerosol index
real (kind=dbl_kind) :: &
taer , & ! total aerosol extinction optical depth
waer , & ! total aerosol single scatter albedo
gaer , & ! total aerosol asymmetry parameter
swdr , & ! shortwave down at surface, direct (W/m^2)
swdf , & ! shortwave down at surface, diffuse (W/m^2)
rnilyr , & ! real(nilyr)
rnslyr , & ! real(nslyr)
rns , & ! real(ns)
tmp_0, tmp_ks, tmp_kl ! temp variables
integer(kind=int_kind), dimension(0:klev) :: &
k_bcini , & !
k_bcins , &
k_bcexs
real(kind=dbl_kind):: &
tmp_gs, tmp1 ! temp variables
! snow grain radii (micro-meters) for table
real (kind=dbl_kind), dimension(nmbrad), parameter :: &
rsnw_tab = (/ & ! snow grain radius for each table entry (micro-meters)
5._dbl_kind, 7._dbl_kind, 10._dbl_kind, 15._dbl_kind, &
20._dbl_kind, 30._dbl_kind, 40._dbl_kind, 50._dbl_kind, &
65._dbl_kind, 80._dbl_kind, 100._dbl_kind, 120._dbl_kind, &
140._dbl_kind, 170._dbl_kind, 200._dbl_kind, 240._dbl_kind, &
290._dbl_kind, 350._dbl_kind, 420._dbl_kind, 500._dbl_kind, &
570._dbl_kind, 660._dbl_kind, 760._dbl_kind, 870._dbl_kind, &
1000._dbl_kind, 1100._dbl_kind, 1250._dbl_kind, 1400._dbl_kind, &
1600._dbl_kind, 1800._dbl_kind, 2000._dbl_kind, 2500._dbl_kind/)
! snow extinction efficiency (unitless)
real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
Qs_tab = reshape((/ &
2.131798_dbl_kind, 2.187756_dbl_kind, 2.267358_dbl_kind, &
2.104499_dbl_kind, 2.148345_dbl_kind, 2.236078_dbl_kind, &
2.081580_dbl_kind, 2.116885_dbl_kind, 2.175067_dbl_kind, &
2.062595_dbl_kind, 2.088937_dbl_kind, 2.130242_dbl_kind, &
2.051403_dbl_kind, 2.072422_dbl_kind, 2.106610_dbl_kind, &
2.039223_dbl_kind, 2.055389_dbl_kind, 2.080586_dbl_kind, &
2.032383_dbl_kind, 2.045751_dbl_kind, 2.066394_dbl_kind, &
2.027920_dbl_kind, 2.039388_dbl_kind, 2.057224_dbl_kind, &
2.023444_dbl_kind, 2.033137_dbl_kind, 2.048055_dbl_kind, &
2.020412_dbl_kind, 2.028840_dbl_kind, 2.041874_dbl_kind, &
2.017608_dbl_kind, 2.024863_dbl_kind, 2.036046_dbl_kind, &
2.015592_dbl_kind, 2.022021_dbl_kind, 2.031954_dbl_kind, &
2.014083_dbl_kind, 2.019887_dbl_kind, 2.028853_dbl_kind, &
2.012368_dbl_kind, 2.017471_dbl_kind, 2.025353_dbl_kind, &
2.011092_dbl_kind, 2.015675_dbl_kind, 2.022759_dbl_kind, &
2.009837_dbl_kind, 2.013897_dbl_kind, 2.020168_dbl_kind, &
2.008668_dbl_kind, 2.012252_dbl_kind, 2.017781_dbl_kind, &
2.007627_dbl_kind, 2.010813_dbl_kind, 2.015678_dbl_kind, &
2.006764_dbl_kind, 2.009577_dbl_kind, 2.013880_dbl_kind, &
2.006037_dbl_kind, 2.008520_dbl_kind, 2.012382_dbl_kind, &
2.005528_dbl_kind, 2.007807_dbl_kind, 2.011307_dbl_kind, &
2.005025_dbl_kind, 2.007079_dbl_kind, 2.010280_dbl_kind, &
2.004562_dbl_kind, 2.006440_dbl_kind, 2.009333_dbl_kind, &
2.004155_dbl_kind, 2.005898_dbl_kind, 2.008523_dbl_kind, &
2.003794_dbl_kind, 2.005379_dbl_kind, 2.007795_dbl_kind, &
2.003555_dbl_kind, 2.005041_dbl_kind, 2.007329_dbl_kind, &
2.003264_dbl_kind, 2.004624_dbl_kind, 2.006729_dbl_kind, &
2.003037_dbl_kind, 2.004291_dbl_kind, 2.006230_dbl_kind, &
2.002776_dbl_kind, 2.003929_dbl_kind, 2.005700_dbl_kind, &
2.002590_dbl_kind, 2.003627_dbl_kind, 2.005276_dbl_kind, &
2.002395_dbl_kind, 2.003391_dbl_kind, 2.004904_dbl_kind, &
2.002071_dbl_kind, 2.002922_dbl_kind, 2.004241_dbl_kind/), &
(/nspint,nmbrad/))
! snow single scattering albedo (unitless)
real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
ws_tab = reshape((/ &
0.9999994_dbl_kind, 0.9999673_dbl_kind, 0.9954589_dbl_kind, &
0.9999992_dbl_kind, 0.9999547_dbl_kind, 0.9938576_dbl_kind, &
0.9999990_dbl_kind, 0.9999382_dbl_kind, 0.9917989_dbl_kind, &
0.9999985_dbl_kind, 0.9999123_dbl_kind, 0.9889724_dbl_kind, &
0.9999979_dbl_kind, 0.9998844_dbl_kind, 0.9866190_dbl_kind, &
0.9999970_dbl_kind, 0.9998317_dbl_kind, 0.9823021_dbl_kind, &
0.9999960_dbl_kind, 0.9997800_dbl_kind, 0.9785269_dbl_kind, &
0.9999951_dbl_kind, 0.9997288_dbl_kind, 0.9751601_dbl_kind, &
0.9999936_dbl_kind, 0.9996531_dbl_kind, 0.9706974_dbl_kind, &
0.9999922_dbl_kind, 0.9995783_dbl_kind, 0.9667577_dbl_kind, &
0.9999903_dbl_kind, 0.9994798_dbl_kind, 0.9621007_dbl_kind, &
0.9999885_dbl_kind, 0.9993825_dbl_kind, 0.9579541_dbl_kind, &
0.9999866_dbl_kind, 0.9992862_dbl_kind, 0.9541924_dbl_kind, &
0.9999838_dbl_kind, 0.9991434_dbl_kind, 0.9490959_dbl_kind, &
0.9999810_dbl_kind, 0.9990025_dbl_kind, 0.9444940_dbl_kind, &
0.9999772_dbl_kind, 0.9988171_dbl_kind, 0.9389141_dbl_kind, &
0.9999726_dbl_kind, 0.9985890_dbl_kind, 0.9325819_dbl_kind, &
0.9999670_dbl_kind, 0.9983199_dbl_kind, 0.9256405_dbl_kind, &
0.9999605_dbl_kind, 0.9980117_dbl_kind, 0.9181533_dbl_kind, &
0.9999530_dbl_kind, 0.9976663_dbl_kind, 0.9101540_dbl_kind, &
0.9999465_dbl_kind, 0.9973693_dbl_kind, 0.9035031_dbl_kind, &
0.9999382_dbl_kind, 0.9969939_dbl_kind, 0.8953134_dbl_kind, &
0.9999289_dbl_kind, 0.9965848_dbl_kind, 0.8865789_dbl_kind, &
0.9999188_dbl_kind, 0.9961434_dbl_kind, 0.8773350_dbl_kind, &
0.9999068_dbl_kind, 0.9956323_dbl_kind, 0.8668233_dbl_kind, &
0.9998975_dbl_kind, 0.9952464_dbl_kind, 0.8589990_dbl_kind, &
0.9998837_dbl_kind, 0.9946782_dbl_kind, 0.8476493_dbl_kind, &
0.9998699_dbl_kind, 0.9941218_dbl_kind, 0.8367318_dbl_kind, &
0.9998515_dbl_kind, 0.9933966_dbl_kind, 0.8227881_dbl_kind, &
0.9998332_dbl_kind, 0.9926888_dbl_kind, 0.8095131_dbl_kind, &
0.9998148_dbl_kind, 0.9919968_dbl_kind, 0.7968620_dbl_kind, &
0.9997691_dbl_kind, 0.9903277_dbl_kind, 0.7677887_dbl_kind/), &
(/nspint,nmbrad/))
! snow asymmetry parameter (unitless)
real (kind=dbl_kind), dimension (nspint,nmbrad), parameter :: &
gs_tab = reshape((/ &
0.859913_dbl_kind, 0.848003_dbl_kind, 0.824415_dbl_kind, &
0.867130_dbl_kind, 0.858150_dbl_kind, 0.848445_dbl_kind, &
0.873381_dbl_kind, 0.867221_dbl_kind, 0.861714_dbl_kind, &
0.878368_dbl_kind, 0.874879_dbl_kind, 0.874036_dbl_kind, &
0.881462_dbl_kind, 0.879661_dbl_kind, 0.881299_dbl_kind, &
0.884361_dbl_kind, 0.883903_dbl_kind, 0.890184_dbl_kind, &
0.885937_dbl_kind, 0.886256_dbl_kind, 0.895393_dbl_kind, &
0.886931_dbl_kind, 0.887769_dbl_kind, 0.899072_dbl_kind, &
0.887894_dbl_kind, 0.889255_dbl_kind, 0.903285_dbl_kind, &
0.888515_dbl_kind, 0.890236_dbl_kind, 0.906588_dbl_kind, &
0.889073_dbl_kind, 0.891127_dbl_kind, 0.910152_dbl_kind, &
0.889452_dbl_kind, 0.891750_dbl_kind, 0.913100_dbl_kind, &
0.889730_dbl_kind, 0.892213_dbl_kind, 0.915621_dbl_kind, &
0.890026_dbl_kind, 0.892723_dbl_kind, 0.918831_dbl_kind, &
0.890238_dbl_kind, 0.893099_dbl_kind, 0.921540_dbl_kind, &
0.890441_dbl_kind, 0.893474_dbl_kind, 0.924581_dbl_kind, &
0.890618_dbl_kind, 0.893816_dbl_kind, 0.927701_dbl_kind, &
0.890762_dbl_kind, 0.894123_dbl_kind, 0.930737_dbl_kind, &
0.890881_dbl_kind, 0.894397_dbl_kind, 0.933568_dbl_kind, &
0.890975_dbl_kind, 0.894645_dbl_kind, 0.936148_dbl_kind, &
0.891035_dbl_kind, 0.894822_dbl_kind, 0.937989_dbl_kind, &
0.891097_dbl_kind, 0.895020_dbl_kind, 0.939949_dbl_kind, &
0.891147_dbl_kind, 0.895212_dbl_kind, 0.941727_dbl_kind, &
0.891189_dbl_kind, 0.895399_dbl_kind, 0.943339_dbl_kind, &
0.891225_dbl_kind, 0.895601_dbl_kind, 0.944915_dbl_kind, &
0.891248_dbl_kind, 0.895745_dbl_kind, 0.945950_dbl_kind, &
0.891277_dbl_kind, 0.895951_dbl_kind, 0.947288_dbl_kind, &
0.891299_dbl_kind, 0.896142_dbl_kind, 0.948438_dbl_kind, &
0.891323_dbl_kind, 0.896388_dbl_kind, 0.949762_dbl_kind, &
0.891340_dbl_kind, 0.896623_dbl_kind, 0.950916_dbl_kind, &
0.891356_dbl_kind, 0.896851_dbl_kind, 0.951945_dbl_kind, &
0.891386_dbl_kind, 0.897399_dbl_kind, 0.954156_dbl_kind/), &
(/nspint,nmbrad/))
! inherent optical property (iop) arrays for ice and ponded ice
! mn = specified mean (or base) value
! ki = extinction coefficient (/m)
! wi = single scattering albedo
! gi = asymmetry parameter
! ice surface scattering layer (ssl) iops
real (kind=dbl_kind), dimension (nspint), parameter :: &
ki_ssl_mn = (/ 1000.1_dbl_kind, 1003.7_dbl_kind, 7042._dbl_kind/), &
wi_ssl_mn = (/ .9999_dbl_kind, .9963_dbl_kind, .9088_dbl_kind/), &
gi_ssl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind/)
! ice drained layer (dl) iops
real (kind=dbl_kind), dimension (nspint), parameter :: &
ki_dl_mn = (/ 100.2_dbl_kind, 107.7_dbl_kind, 1309._dbl_kind /), &
wi_dl_mn = (/ .9980_dbl_kind, .9287_dbl_kind, .0305_dbl_kind /), &
gi_dl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
! ice interior layer (int) iops
real (kind=dbl_kind), dimension (nspint), parameter :: &
ki_int_mn = (/ 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind /), &
wi_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind /), &
gi_int_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
! ponded ice surface scattering layer (ssl) iops
real (kind=dbl_kind), dimension (nspint), parameter :: &
ki_p_ssl_mn = (/ 70.2_dbl_kind, 77.7_dbl_kind, 1309._dbl_kind/), &
wi_p_ssl_mn = (/ .9972_dbl_kind, .9009_dbl_kind, .0305_dbl_kind/), &
gi_p_ssl_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
! ponded ice interior layer (int) iops
real (kind=dbl_kind), dimension (nspint), parameter :: &
ki_p_int_mn = (/ 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind/), &
wi_p_int_mn = (/ .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind/), &
gi_p_int_mn = (/ .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /)
! inherent optical property (iop) arrays for pond water and underlying ocean
! kw = Pond water extinction coefficient (/m)
! ww = Pond water single scattering albedo
! gw = Pond water asymmetry parameter
real (kind=dbl_kind), dimension (nspint), parameter :: &
kw = (/ 0.20_dbl_kind, 12.0_dbl_kind, 729._dbl_kind /), &
ww = (/ 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /), &
gw = (/ 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /)
real (kind=dbl_kind), parameter :: &
rhoi = 917.0_dbl_kind,& ! pure ice mass density (kg/m3)
fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max
fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min
! tuning parameters
! ice and pond scat coeff fractional change for +- one-sigma in albedo
fp_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for + stn dev in alb
fm_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for - stn dev in alb
fp_pnd = 2.00_dbl_kind, & ! ponded ice fraction of scat coeff for + stn dev in alb
fm_pnd = 0.50_dbl_kind ! ponded ice fraction of scat coeff for - stn dev in alb
real (kind=dbl_kind), parameter :: & !chla-specific absorption coefficient
kchl_tab = 0.01 !0.0023-0.0029 Perovich 1993, also 0.0067 m^2 (mg Chl)^-1
! found values of 0.006 to 0.023 m^2/ mg (676 nm) Neukermans 2014
! and averages over the 300-700nm of 0.0075 m^2/mg in ice Fritsen (2011)
! at 440nm values as high as 0.2 m^2/mg in under ice bloom (Balch 2014)
! Grenfell 1991 uses 0.004 (m^2/mg) which is (0.0078 * spectral weighting)
!chlorophyll mass extinction cross section (m^2/mg chla)
character(len=*),parameter :: subname='(compute_dEdd)'
!-----------------------------------------------------------------------
! Initialize and tune bare ice/ponded ice iops
k_bcini(:) = c0
k_bcins(:) = c0
k_bcexs(:) = c0
rnilyr = c1/real(nilyr,kind=dbl_kind)
rnslyr = c1/real(nslyr,kind=dbl_kind)
kii = nslyr + 1
! initialize albedos and fluxes to 0
fthrul = c0
Iabs = c0
kabs_chl(:,:) = c0
tzaer(:,:) = c0
wzaer(:,:) = c0
gzaer(:,:) = c0
avdr = c0
avdf = c0
aidr = c0
aidf = c0
fsfc = c0
fint = c0
fthru = c0
fthruvdr = c0
fthruvdf = c0
fthruidr = c0
fthruidf = c0
! spectral weights
! weights 2 (0.7-1.19 micro-meters) and 3 (1.19-5.0 micro-meters)
! are chosen based on 1D calculations using ratio of direct to total
! near-infrared solar (0.7-5.0 micro-meter) which indicates clear/cloudy
! conditions: more cloud, the less 1.19-5.0 relative to the
! 0.7-1.19 micro-meter due to cloud absorption.
wghtns(1) = c1
wghtns(2) = cp67 + (cp78-cp67)*(c1-fnidr)
! wghtns(3) = cp33 + (cp22-cp33)*(c1-fnidr)
wghtns(3) = c1 - wghtns(2)
! find snow grain adjustment factor, dependent upon clear/overcast sky
! estimate. comparisons with SNICAR show better agreement with DE when
! this factor is included (clear sky near 1 and overcast near 0.8 give
! best agreement). Multiply by rnsw here for efficiency.
do k = 1, nslyr
frsnw(k) = (fr_max*fnidr + fr_min*(c1-fnidr))*rsnw(k)
Sabs(k) = c0
enddo
! layer thicknesses
! snow
dz = hs*rnslyr
! for small enough snow thickness, ssl thickness half of top snow layer
!ech: note this is highly resolution dependent!
dzk(0) = min(hs_ssl, dz/c2)
dzk(1) = dz - dzk(0)
if (nslyr > 1) then
do k = 2, nslyr
dzk(k) = dz
enddo
endif
! ice
dz = hi*rnilyr
! empirical reduction in sea ice ssl thickness for ice thinner than 1.5m;
! factor of 30 gives best albedo comparison with limited observations
dz_ssl = hi_ssl
!ech: note hardwired parameters
! if( hi < 1.5_dbl_kind ) dz_ssl = hi/30._dbl_kind
dz_ssl = min(hi_ssl, hi/30._dbl_kind)
! set sea ice ssl thickness to half top layer if sea ice thin enough
!ech: note this is highly resolution dependent!
dz_ssl = min(dz_ssl, dz/c2)
dzk(kii) = dz_ssl
dzk(kii+1) = dz - dz_ssl
if (kii+2 <= klev) then
do k = kii+2, klev
dzk(k) = dz
enddo
endif
! adjust sea ice iops with tuning parameters; tune only the
! scattering coefficient by factors of R_ice, R_pnd, where
! R values of +1 correspond approximately to +1 sigma changes in albedo, and
! R values of -1 correspond approximately to -1 sigma changes in albedo
! Note: the albedo change becomes non-linear for R values > +1 or < -1
if( R_ice >= c0 ) then
do ns = 1, nspint
sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fp_ice*R_ice)
ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
wi_ssl(ns) = sigp/ki_ssl(ns)
gi_ssl(ns) = gi_ssl_mn(ns)
sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fp_ice*R_ice)
ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
wi_dl(ns) = sigp/ki_dl(ns)
gi_dl(ns) = gi_dl_mn(ns)
sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fp_ice*R_ice)
ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
wi_int(ns) = sigp/ki_int(ns)
gi_int(ns) = gi_int_mn(ns)
enddo
else !if( R_ice < c0 ) then
do ns = 1, nspint
sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fm_ice*R_ice)
sigp = max(sigp, c0)
ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
wi_ssl(ns) = sigp/ki_ssl(ns)
gi_ssl(ns) = gi_ssl_mn(ns)
sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fm_ice*R_ice)
sigp = max(sigp, c0)
ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
wi_dl(ns) = sigp/ki_dl(ns)
gi_dl(ns) = gi_dl_mn(ns)
sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fm_ice*R_ice)
sigp = max(sigp, c0)
ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
wi_int(ns) = sigp/ki_int(ns)
gi_int(ns) = gi_int_mn(ns)
enddo
endif ! adjust ice iops
! adjust ponded ice iops with tuning parameters
if( R_pnd >= c0 ) then
do ns = 1, nspint
sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fp_pnd*R_pnd)
ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
gi_p_ssl(ns) = gi_p_ssl_mn(ns)
sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fp_pnd*R_pnd)
ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
wi_p_int(ns) = sigp/ki_p_int(ns)
gi_p_int(ns) = gi_p_int_mn(ns)
enddo
else !if( R_pnd < c0 ) then
do ns = 1, nspint
sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fm_pnd*R_pnd)
sigp = max(sigp, c0)
ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
gi_p_ssl(ns) = gi_p_ssl_mn(ns)
sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fm_pnd*R_pnd)
sigp = max(sigp, c0)
ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
wi_p_int(ns) = sigp/ki_p_int(ns)
gi_p_int(ns) = gi_p_int_mn(ns)
enddo
endif ! adjust ponded ice iops
! use srftyp to determine interface index of surface absorption
if (srftyp == 1) then
! snow covered sea ice
ksrf = 1
else
! bare sea ice or ponded ice
ksrf = nslyr + 2
endif
if (tr_bgc_N .and. dEdd_algae) then ! compute kabs_chl for chlorophyll
do k = 0, klev
kabs_chl(1,k) = kchl_tab*zbio(nlt_chl_sw+k)
enddo
else
k = klev
kabs_chl(1,k) = kalg*(0.50_dbl_kind/dzk(k))
endif ! kabs_chl
!mgf++
if (modal_aero) then
do k=0,klev
if (k < nslyr+1) then ! define indices for snow layer
! use top rsnw, rhosnw for snow ssl and rest of top layer
ksnow = k - min(k-1,0)
tmp_gs = frsnw(ksnow)
! get grain size index:
! works for 25 < snw_rds < 1625 um:
if (tmp_gs < 125) then
tmp1 = tmp_gs/50
k_bcini(k) = nint(tmp1)
elseif (tmp_gs < 175) then
k_bcini(k) = 2
else
tmp1 = (tmp_gs/250)+2
k_bcini(k) = nint(tmp1)
endif
else ! use the largest snow grain size for ice
k_bcini(k) = 8
endif
! Set index corresponding to BC effective radius. Here,
! asssume constant BC effective radius of 100nm
! (corresponding to index 2)
k_bcins(k) = 2
k_bcexs(k) = 2
! check bounds:
if (k_bcini(k) < 1) k_bcini(k) = 1
if (k_bcini(k) > 8) k_bcini(k) = 8
if (k_bcins(k) < 1) k_bcins(k) = 1
if (k_bcins(k) > 10) k_bcins(k) = 10
if (k_bcexs(k) < 1) k_bcexs(k) = 1
if (k_bcexs(k) > 10) k_bcexs(k) = 10
! print ice radius index:
! write(warnstr,*) subname, "MGFICE2:k, ice index= ",k, k_bcini(k)
! call icepack_warnings_add(warnstr)
enddo ! k
if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll
do n = 1,n_zaero
if (n == 1) then ! interstitial BC
do k = 0, klev
do ns = 1,nspint ! not weighted by aice
tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))* &
gaer_bc_tab(ns,k_bcexs(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
enddo ! nspint
enddo
elseif (n==2) then ! within-ice BC
do k = 0, klev
do ns = 1,nspint
tzaer(ns,k) = tzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k)) * &
bcenh(ns,k_bcins(k),k_bcini(k))* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
wzaer(ns,k) = wzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
gzaer(ns,k) = gzaer(ns,k)+kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))* &
gaer_bc_tab(ns,k_bcins(k))*zbio(nlt_zaero_sw(n)+k)*dzk(k)
enddo ! nspint
enddo
else ! dust
do k = 0, klev
do ns = 1,nspint ! not weighted by aice
tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
enddo ! nspint
enddo
endif !(n=1)
enddo ! n_zaero
endif ! tr_zaero and dEdd_algae
else ! Bulk aerosol treatment
if (tr_zaero .and. dEdd_algae) then ! compute kzaero for chlorophyll
do n = 1,n_zaero ! multiply by aice?
do k = 0, klev
do ns = 1,nspint ! not weighted by aice
tzaer(ns,k) = tzaer(ns,k)+kaer_tab(ns,n)* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
wzaer(ns,k) = wzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
zbio(nlt_zaero_sw(n)+k)*dzk(k)
gzaer(ns,k) = gzaer(ns,k)+kaer_tab(ns,n)*waer_tab(ns,n)* &
gaer_tab(ns,n)*zbio(nlt_zaero_sw(n)+k)*dzk(k)
enddo ! nspint
enddo
enddo
endif !tr_zaero
endif ! modal_aero
!-----------------------------------------------------------------------
! begin spectral loop
do ns = 1, nspint
! set optical properties of air/snow/pond overlying sea ice
! air
if( srftyp == 0 ) then
do k=0,nslyr
tau(k) = c0
w0(k) = c0
g(k) = c0
enddo
! snow
else if( srftyp == 1 ) then
! interpolate snow iops using input snow grain radius,
! snow density and tabular data
do k=0,nslyr
! use top rsnw, rhosnw for snow ssl and rest of top layer
ksnow = k - min(k-1,0)
! find snow iops using input snow density and snow grain radius:
if( frsnw(ksnow) < rsnw_tab(1) ) then
Qs = Qs_tab(ns,1)
ws = ws_tab(ns,1)
gs = gs_tab(ns,1)
else if( frsnw(ksnow) >= rsnw_tab(nmbrad) ) then
Qs = Qs_tab(ns,nmbrad)
ws = ws_tab(ns,nmbrad)
gs = gs_tab(ns,nmbrad)
else
! linear interpolation in rsnw
do nr=2,nmbrad
if( rsnw_tab(nr-1) <= frsnw(ksnow) .and. &
frsnw(ksnow) < rsnw_tab(nr)) then
delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / &
(rsnw_tab(nr) - rsnw_tab(nr-1))
Qs = Qs_tab(ns,nr-1)*(c1-delr) + &
Qs_tab(ns,nr)*delr
ws = ws_tab(ns,nr-1)*(c1-delr) + &
ws_tab(ns,nr)*delr
gs = gs_tab(ns,nr-1)*(c1-delr) + &
gs_tab(ns,nr)*delr
endif
enddo ! nr
endif
ks = Qs*((rhosnw(ksnow)/rhoi)*3._dbl_kind / &
(4._dbl_kind*frsnw(ksnow)*1.0e-6_dbl_kind))
tau(k) = (ks + kabs_chl(ns,k))*dzk(k)
w0(k) = ks/(ks + kabs_chl(ns,k)) *ws
g(k) = gs
enddo ! k
! aerosol in snow
if (tr_zaero .and. dEdd_algae) then
do k = 0,nslyr
gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
(w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
(tau(k) + tzaer(ns,k))
tau(k) = tau(k) + tzaer(ns,k)
enddo
elseif (tr_aero) then
k = 0 ! snow SSL
taer = c0
waer = c0
gaer = c0
do na=1,4*n_aero,4
! mgf++
if (modal_aero) then
if (na == 1) then
!interstitial BC
taer = taer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))
waer = waer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))
gaer = gaer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
elseif (na == 5)then
!within-ice BC
taer = taer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
bcenh(ns,k_bcins(k),k_bcini(k))
waer = waer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))
gaer = gaer + &
aero_mp(na)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
else
! other species (dust)
taer = taer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif
else
taer = taer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
aero_mp(na)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif !modal_aero
!mgf--
enddo ! na
gaer = gaer/(waer+puny)
waer = waer/(taer+puny)
do k=1,nslyr
taer = c0
waer = c0
gaer = c0
do na=1,4*n_aero,4
if (modal_aero) then
!mgf++
if (na==1) then
! interstitial BC
taer = taer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))
waer = waer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))
gaer = gaer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
elseif (na==5) then
! within-ice BC
taer = taer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))*&
bcenh(ns,k_bcins(k),k_bcini(k))
waer = waer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))
gaer = gaer + &
(aero_mp(na+1)/rnslyr)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
else
! other species (dust)
taer = taer + &
(aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
(aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
(aero_mp(na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif !(na==1)
else
taer = taer + &
(aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
(aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
(aero_mp(na+1)*rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif ! modal_aero
!mgf--
enddo ! na
gaer = gaer/(waer+puny)
waer = waer/(taer+puny)
g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
(w0(k)*tau(k) + waer*taer)
w0(k) = (w0(k)*tau(k) + waer*taer) / &
(tau(k) + taer)
tau(k) = tau(k) + taer
enddo ! k
endif ! tr_aero
! pond
else !if( srftyp == 2 ) then
! pond water layers evenly spaced
dz = hp/(c1/rnslyr+c1)
do k=0,nslyr
tau(k) = kw(ns)*dz
w0(k) = ww(ns)
g(k) = gw(ns)
! no aerosol in pond
enddo ! k
endif ! srftyp
! set optical properties of sea ice
! bare or snow-covered sea ice layers
if( srftyp <= 1 ) then
! ssl
k = kii
tau(k) = (ki_ssl(ns)+kabs_chl(ns,k))*dzk(k)
w0(k) = ki_ssl(ns)/(ki_ssl(ns) + kabs_chl(ns,k))*wi_ssl(ns)
g(k) = gi_ssl(ns)
! dl
k = kii + 1
! scale dz for dl relative to 4 even-layer-thickness 1.5m case
fs = p25/rnilyr
tau(k) = (ki_dl(ns) + kabs_chl(ns,k)) *dzk(k)*fs
w0(k) = ki_dl(ns)/(ki_dl(ns) + kabs_chl(ns,k)) *wi_dl(ns)
g(k) = gi_dl(ns)
! int above lowest layer
if (kii+2 <= klev-1) then
do k = kii+2, klev-1
tau(k) = (ki_int(ns) + kabs_chl(ns,k))*dzk(k)
w0(k) = ki_int(ns)/(ki_int(ns) + kabs_chl(ns,k)) *wi_int(ns)
g(k) = gi_int(ns)
enddo
endif
! lowest layer
k = klev
! add algae to lowest sea ice layer, visible only:
kabs = ki_int(ns)*(c1-wi_int(ns))
if( ns == 1 ) then
! total layer absorption optical depth fixed at value
! of kalg*0.50m, independent of actual layer thickness
kabs = kabs + kabs_chl(ns,k)
endif
sig = ki_int(ns)*wi_int(ns)
tau(k) = (kabs+sig)*dzk(k)
w0(k) = (sig/(sig+kabs))
g(k) = gi_int(ns)
! aerosol in sea ice
if (tr_zaero .and. dEdd_algae) then
do k = kii, klev
gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
g(k) = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) / &
(w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k))
w0(k) = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) / &
(tau(k) + tzaer(ns,k))
tau(k) = tau(k) + tzaer(ns,k)
enddo
elseif (tr_aero) then
k = kii ! sea ice SSL
taer = c0
waer = c0
gaer = c0
do na=1,4*n_aero,4
!mgf++
if (modal_aero) then
if (na==1) then
! interstitial BC
taer = taer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))
waer = waer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))
gaer = gaer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
elseif (na==5) then
! within-ice BC
taer = taer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
bcenh(ns,k_bcins(k),k_bcini(k))
waer = waer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))
gaer = gaer + &
aero_mp(na+2)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
else
! other species (dust)
taer = taer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif
else !bulk
taer = taer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
aero_mp(na+2)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif ! modal_aero
!mgf--
enddo ! na
gaer = gaer/(waer+puny)
waer = waer/(taer+puny)
g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
(w0(k)*tau(k) + waer*taer)
w0(k) = (w0(k)*tau(k) + waer*taer) / &
(tau(k) + taer)
tau(k) = tau(k) + taer
do k = kii+1, klev
taer = c0
waer = c0
gaer = c0
do na=1,4*n_aero,4
!mgf++
if (modal_aero) then
if (na==1) then
! interstitial BC
taer = taer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))
waer = waer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))
gaer = gaer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcexs(k))* &
waer_bc_tab(ns,k_bcexs(k))*gaer_bc_tab(ns,k_bcexs(k))
elseif (na==5) then
! within-ice BC
taer = taer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
bcenh(ns,k_bcins(k),k_bcini(k))
waer = waer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))
gaer = gaer + &
(aero_mp(na+3)/rnilyr)*kaer_bc_tab(ns,k_bcins(k))* &
waer_bc_tab(ns,k_bcins(k))*gaer_bc_tab(ns,k_bcins(k))
else
! other species (dust)
taer = taer + &
(aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
(aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
(aero_mp(na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif
else !bulk
taer = taer + &
(aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))
waer = waer + &
(aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))
gaer = gaer + &
(aero_mp(na+3)*rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
endif ! modal_aero
!mgf--
enddo ! na
gaer = gaer/(waer+puny)
waer = waer/(taer+puny)
g(k) = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
(w0(k)*tau(k) + waer*taer)
w0(k) = (w0(k)*tau(k) + waer*taer) / &
(tau(k) + taer)
tau(k) = tau(k) + taer
enddo ! k
endif ! tr_aero
! sea ice layers under ponds
else !if( srftyp == 2 ) then
k = kii
tau(k) = ki_p_ssl(ns)*dzk(k)
w0(k) = wi_p_ssl(ns)
g(k) = gi_p_ssl(ns)
k = kii + 1
tau(k) = ki_p_int(ns)*dzk(k)
w0(k) = wi_p_int(ns)
g(k) = gi_p_int(ns)
if (kii+2 <= klev) then
do k = kii+2, klev
tau(k) = ki_p_int(ns)*dzk(k)
w0(k) = wi_p_int(ns)
g(k) = gi_p_int(ns)
enddo ! k
endif
! adjust pond iops if pond depth within specified range
if( hpmin <= hp .and. hp <= hp0 ) then
k = kii
sig_i = ki_ssl(ns)*wi_ssl(ns)
sig_p = ki_p_ssl(ns)*wi_p_ssl(ns)
sig = sig_i + (sig_p-sig_i)*(hp/hp0)
kext = sig + ki_p_ssl(ns)*(c1-wi_p_ssl(ns))
tau(k) = kext*dzk(k)
w0(k) = sig/kext
g(k) = gi_p_int(ns)
k = kii + 1
! scale dz for dl relative to 4 even-layer-thickness 1.5m case
fs = p25/rnilyr
sig_i = ki_dl(ns)*wi_dl(ns)*fs
sig_p = ki_p_int(ns)*wi_p_int(ns)
sig = sig_i + (sig_p-sig_i)*(hp/hp0)
kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
tau(k) = kext*dzk(k)
w0(k) = sig/kext
g(k) = gi_p_int(ns)
if (kii+2 <= klev) then
do k = kii+2, klev
sig_i = ki_int(ns)*wi_int(ns)
sig_p = ki_p_int(ns)*wi_p_int(ns)
sig = sig_i + (sig_p-sig_i)*(hp/hp0)
kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
tau(k) = kext*dzk(k)
w0(k) = sig/kext
g(k) = gi_p_int(ns)
enddo ! k
endif
endif ! small pond depth transition to bare sea ice
endif ! srftyp
! set reflectivities for ocean underlying sea ice
rns = real(ns-1, kind=dbl_kind)
albodr = cp01 * (c1 - min(rns, c1))
albodf = cp01 * (c1 - min(rns, c1))
! layer input properties now completely specified: tau, w0, g,
! albodr, albodf; now compute the Delta-Eddington solution
! reflectivities and transmissivities for each layer; then,
! combine the layers going downwards accounting for multiple
! scattering between layers, and finally start from the
! underlying ocean and combine successive layers upwards to
! the surface; see comments in solution_dEdd for more details.
call solution_dEdd &
(coszen, srftyp, klev, klevp, nslyr, &
tau, w0, g, albodr, albodf, &
trndir, trntdr, trndif, rupdir, rupdif, &
rdndif)
if (icepack_warnings_aborted(subname)) return
! the interface reflectivities and transmissivities required
! to evaluate interface fluxes are returned from solution_dEdd;
! now compute up and down fluxes for each interface, using the
! combined layer properties at each interface:
!
! layers interface
!
! --------------------- k
! k
! ---------------------
do k = 0, klevp
! interface scattering
refk = c1/(c1 - rdndif(k)*rupdif(k))
! dir tran ref from below times interface scattering, plus diff
! tran and ref from below times interface scattering
! fdirup(k) = (trndir(k)*rupdir(k) + &
! (trntdr(k)-trndir(k)) &
! *rupdif(k))*refk
! dir tran plus total diff trans times interface scattering plus
! dir tran with up dir ref and down dif ref times interface scattering
! fdirdn(k) = trndir(k) + (trntdr(k) &
! - trndir(k) + trndir(k) &
! *rupdir(k)*rdndif(k))*refk
! diffuse tran ref from below times interface scattering
! fdifup(k) = trndif(k)*rupdif(k)*refk
! diffuse tran times interface scattering
! fdifdn(k) = trndif(k)*refk
! dfdir = fdirdn - fdirup
dfdir(k) = trndir(k) &
+ (trntdr(k)-trndir(k)) * (c1 - rupdif(k)) * refk &
- trndir(k)*rupdir(k) * (c1 - rdndif(k)) * refk
if (dfdir(k) < puny) dfdir(k) = c0 !echmod necessary?
! dfdif = fdifdn - fdifup
dfdif(k) = trndif(k) * (c1 - rupdif(k)) * refk
if (dfdif(k) < puny) dfdif(k) = c0 !echmod necessary?
enddo ! k
! calculate final surface albedos and fluxes-
! all absorbed flux above ksrf is included in surface absorption
if( ns == 1) then ! visible
swdr = swvdr
swdf = swvdf
avdr = rupdir(0)
avdf = rupdif(0)
tmp_0 = dfdir(0 )*swdr + dfdif(0 )*swdf
tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
! for layer biology: save visible only
do k = nslyr+2, klevp ! Start at DL layer of ice after SSL scattering
fthrul(k-nslyr-1) = dfdir(k)*swdr + dfdif(k)*swdf
enddo
fsfc = fsfc + tmp_0 - tmp_ks
fint = fint + tmp_ks - tmp_kl
fthru = fthru + tmp_kl
fthruvdr = fthruvdr + dfdir(klevp)*swdr
fthruvdf = fthruvdf + dfdif(klevp)*swdf
! if snow covered ice, set snow internal absorption; else, Sabs=0
if( srftyp == 1 ) then
ki = 0
do k=1,nslyr
! skip snow SSL, since SSL absorption included in the surface
! absorption fsfc above
km = k
kp = km + 1
ki = ki + 1
Sabs(ki) = Sabs(ki) &
+ dfdir(km)*swdr + dfdif(km)*swdf &
- (dfdir(kp)*swdr + dfdif(kp)*swdf)
enddo ! k
endif
! complex indexing to insure proper absorptions for sea ice
ki = 0
do k=nslyr+2,nslyr+1+nilyr
! for bare ice, DL absorption for sea ice layer 1
km = k
kp = km + 1
! modify for top sea ice layer for snow over sea ice
if( srftyp == 1 ) then
! must add SSL and DL absorption for sea ice layer 1
if( k == nslyr+2 ) then
km = k - 1
kp = km + 2
endif
endif
ki = ki + 1
Iabs(ki) = Iabs(ki) &
+ dfdir(km)*swdr + dfdif(km)*swdf &
- (dfdir(kp)*swdr + dfdif(kp)*swdf)
enddo ! k
else !if(ns > 1) then ! near IR
swdr = swidr
swdf = swidf
! let fr1 = alb_1*swd*wght1 and fr2 = alb_2*swd*wght2 be the ns=2,3
! reflected fluxes respectively, where alb_1, alb_2 are the band
! albedos, swd = nir incident shortwave flux, and wght1, wght2 are
! the 2,3 band weights. thus, the total reflected flux is:
! fr = fr1 + fr2 = alb_1*swd*wght1 + alb_2*swd*wght2 hence, the
! 2,3 nir band albedo is alb = fr/swd = alb_1*wght1 + alb_2*wght2
aidr = aidr + rupdir(0)*wghtns(ns)
aidf = aidf + rupdif(0)*wghtns(ns)
tmp_0 = dfdir(0 )*swdr + dfdif(0 )*swdf
tmp_ks = dfdir(ksrf )*swdr + dfdif(ksrf )*swdf
tmp_kl = dfdir(klevp)*swdr + dfdif(klevp)*swdf
tmp_0 = tmp_0 * wghtns(ns)
tmp_ks = tmp_ks * wghtns(ns)
tmp_kl = tmp_kl * wghtns(ns)
fsfc = fsfc + tmp_0 - tmp_ks
fint = fint + tmp_ks - tmp_kl
fthru = fthru + tmp_kl
fthruidr = fthruidr + dfdir(klevp)*swdr*wghtns(ns)
fthruidf = fthruidf + dfdif(klevp)*swdf*wghtns(ns)
! if snow covered ice, set snow internal absorption; else, Sabs=0
if( srftyp == 1 ) then
ki = 0
do k=1,nslyr
! skip snow SSL, since SSL absorption included in the surface
! absorption fsfc above
km = k
kp = km + 1
ki = ki + 1
Sabs(ki) = Sabs(ki) &
+ (dfdir(km)*swdr + dfdif(km)*swdf &
- (dfdir(kp)*swdr + dfdif(kp)*swdf)) &
* wghtns(ns)
enddo ! k
endif
! complex indexing to insure proper absorptions for sea ice
ki = 0
do k=nslyr+2,nslyr+1+nilyr
! for bare ice, DL absorption for sea ice layer 1
km = k
kp = km + 1
! modify for top sea ice layer for snow over sea ice
if( srftyp == 1 ) then
! must add SSL and DL absorption for sea ice layer 1
if( k == nslyr+2 ) then
km = k - 1
kp = km + 2
endif
endif
ki = ki + 1
Iabs(ki) = Iabs(ki) &
+ (dfdir(km)*swdr + dfdif(km)*swdf &
- (dfdir(kp)*swdr + dfdif(kp)*swdf)) &
* wghtns(ns)
enddo ! k
endif ! ns = 1, ns > 1
enddo ! end spectral loop ns
! accumulate fluxes over bare sea ice
alvdr = avdr
alvdf = avdf
alidr = aidr
alidf = aidf
fswsfc = fswsfc + fsfc *fi
fswint = fswint + fint *fi
fswthru = fswthru + fthru*fi
fswthru_vdr = fswthru_vdr + fthruvdr*fi
fswthru_vdf = fswthru_vdf + fthruvdf*fi
fswthru_idr = fswthru_idr + fthruidr*fi
fswthru_idf = fswthru_idf + fthruidf*fi
do k = 1, nslyr
Sswabs(k) = Sswabs(k) + Sabs(k)*fi
enddo ! k
do k = 1, nilyr
Iswabs(k) = Iswabs(k) + Iabs(k)*fi
! bgc layer
fswpenl(k) = fswpenl(k) + fthrul(k)* fi
if (k == nilyr) then
fswpenl(k+1) = fswpenl(k+1) + fthrul(k+1)*fi
endif
enddo ! k
!----------------------------------------------------------------
! if ice has zero heat capacity, no SW can be absorbed
! in the ice/snow interior, so add to surface absorption.
! Note: nilyr = nslyr = 1 for this case
!----------------------------------------------------------------
if (.not. heat_capacity) then
! SW absorbed at snow/ice surface
fswsfc = fswsfc + Iswabs(1) + Sswabs(1)
! SW absorbed in ice interior
fswint = c0
Iswabs(1) = c0
Sswabs(1) = c0
endif ! heat_capacity
end subroutine compute_dEdd
!=======================================================================
!
! Given input vertical profiles of optical properties, evaluate the
! monochromatic Delta-Eddington solution.
!
! author: Bruce P. Briegleb, NCAR
! 2013: E Hunke merged with NCAR version
subroutine solution_dEdd &
(coszen, srftyp, klev, klevp, nslyr, &
tau, w0, g, albodr, albodf, &
trndir, trntdr, trndif, rupdir, rupdif, &
rdndif)
real (kind=dbl_kind), intent(in) :: &
coszen ! cosine solar zenith angle
integer (kind=int_kind), intent(in) :: &
srftyp , & ! surface type over ice: (0=air, 1=snow, 2=pond)
klev , & ! number of radiation layers - 1
klevp , & ! number of radiation interfaces - 1
! (0 layer is included also)
nslyr ! number of snow layers
real (kind=dbl_kind), dimension(0:klev), intent(in) :: &
tau , & ! layer extinction optical depth
w0 , & ! layer single scattering albedo
g ! layer asymmetry parameter
real (kind=dbl_kind), intent(in) :: &
albodr , & ! ocean albedo to direct rad
albodf ! ocean albedo to diffuse rad
! following arrays are defined at model interfaces; 0 is the top of the
! layer above the sea ice; klevp is the sea ice/ocean interface.
real (kind=dbl_kind), dimension (0:klevp), intent(out) :: &
trndir , & ! solar beam down transmission from top
trntdr , & ! total transmission to direct beam for layers above
trndif , & ! diffuse transmission to diffuse beam for layers above
rupdir , & ! reflectivity to direct radiation for layers below
rupdif , & ! reflectivity to diffuse radiation for layers below
rdndif ! reflectivity to diffuse radiation for layers above
!-----------------------------------------------------------------------
!
! Delta-Eddington solution for snow/air/pond over sea ice
!
! Generic solution for a snow/air/pond input column of klev+1 layers,
! with srftyp determining at what interface fresnel refraction occurs.
!
! Computes layer reflectivities and transmissivities, from the top down
! to the lowest interface using the Delta-Eddington solutions for each
! layer; combines layers from top down to lowest interface, and from the
! lowest interface (underlying ocean) up to the top of the column.
!
! Note that layer diffuse reflectivity and transmissivity are computed
! by integrating the direct over several gaussian angles. This is
! because the diffuse reflectivity expression sometimes is negative,
! but the direct reflectivity is always well-behaved. We assume isotropic
! radiation in the upward and downward hemispheres for this integration.
!
! Assumes monochromatic (spectrally uniform) properties across a band
! for the input optical parameters.
!
! If total transmission of the direct beam to the interface above a particular
! layer is less than trmin, then no further Delta-Eddington solutions are
! evaluated for layers below.
!
! The following describes how refraction is handled in the calculation.
!
! First, we assume that radiation is refracted when entering either
! sea ice at the base of the surface scattering layer, or water (i.e. melt
! pond); we assume that radiation does not refract when entering snow, nor
! upon entering sea ice from a melt pond, nor upon entering the underlying
! ocean from sea ice.
!
! To handle refraction, we define a "fresnel" layer, which physically
! is of neglible thickness and is non-absorbing, which can be combined to
! any sea ice layer or top of melt pond. The fresnel layer accounts for
! refraction of direct beam and associated reflection and transmission for
! solar radiation. A fresnel layer is combined with the top of a melt pond
! or to the surface scattering layer of sea ice if no melt pond lies over it.
!
! Some caution must be exercised for the fresnel layer, because any layer
! to which it is combined is no longer a homogeneous layer, as are all other
! individual layers. For all other layers for example, the direct and diffuse
! reflectivities/transmissivities (R/T) are the same for radiation above or
! below the layer. This is the meaning of homogeneous! But for the fresnel
! layer this is not so. Thus, the R/T for this layer must be distinguished
! for radiation above from that from radiation below. For generality, we
! treat all layers to be combined as inhomogeneous.
!
!-----------------------------------------------------------------------
! local variables
integer (kind=int_kind) :: &
kfrsnl ! radiation interface index for fresnel layer
! following variables are defined for each layer; 0 refers to the top
! layer. In general we must distinguish directions above and below in
! the diffuse reflectivity and transmissivity, as layers are not assumed
! to be homogeneous (apart from the single layer Delta-Edd solutions);
! the direct is always from above.
real (kind=dbl_kind), dimension (0:klev) :: &
rdir , & ! layer reflectivity to direct radiation
rdif_a , & ! layer reflectivity to diffuse radiation from above
rdif_b , & ! layer reflectivity to diffuse radiation from below
tdir , & ! layer transmission to direct radiation (solar beam + diffuse)
tdif_a , & ! layer transmission to diffuse radiation from above
tdif_b , & ! layer transmission to diffuse radiation from below
trnlay ! solar beam transm for layer (direct beam only)
integer (kind=int_kind) :: &
k ! level index
real (kind=dbl_kind), parameter :: &
trmin = 0.001_dbl_kind ! minimum total transmission allowed
! total transmission is that due to the direct beam; i.e. it includes
! both the directly transmitted solar beam and the diffuse downwards
! transmitted radiation resulting from scattering out of the direct beam
real (kind=dbl_kind) :: &
tautot , & ! layer optical depth
wtot , & ! layer single scattering albedo
gtot , & ! layer asymmetry parameter
ftot , & ! layer forward scattering fraction
ts , & ! layer scaled extinction optical depth
ws , & ! layer scaled single scattering albedo
gs , & ! layer scaled asymmetry parameter
rintfc , & ! reflection (multiple) at an interface
refkp1 , & ! interface multiple scattering for k+1
refkm1 , & ! interface multiple scattering for k-1
tdrrdir , & ! direct tran times layer direct ref
tdndif ! total down diffuse = tot tran - direct tran
! perpendicular and parallel relative to plane of incidence and scattering
real (kind=dbl_kind) :: &
R1 , & ! perpendicular polarization reflection amplitude
R2 , & ! parallel polarization reflection amplitude
T1 , & ! perpendicular polarization transmission amplitude
T2 , & ! parallel polarization transmission amplitude
Rf_dir_a , & ! fresnel reflection to direct radiation
Tf_dir_a , & ! fresnel transmission to direct radiation
Rf_dif_a , & ! fresnel reflection to diff radiation from above
Rf_dif_b , & ! fresnel reflection to diff radiation from below
Tf_dif_a , & ! fresnel transmission to diff radiation from above
Tf_dif_b ! fresnel transmission to diff radiation from below
! refractive index for sea ice, water; pre-computed, band-independent,
! diffuse fresnel reflectivities
real (kind=dbl_kind), parameter :: &
refindx = 1.310_dbl_kind , & ! refractive index of sea ice (water also)
cp063 = 0.063_dbl_kind , & ! diffuse fresnel reflectivity from above
cp455 = 0.455_dbl_kind ! diffuse fresnel reflectivity from below
real (kind=dbl_kind) :: &
mu0 , & ! cosine solar zenith angle incident
mu0nij ! cosine solar zenith angle in medium below fresnel level
real (kind=dbl_kind) :: &
mu0n ! cosine solar zenith angle in medium
real (kind=dbl_kind) :: &
alp , & ! temporary for alpha
gam , & ! temporary for agamm
lm , & ! temporary for el
mu , & ! temporary for gauspt
ne , & ! temporary for n
ue , & ! temporary for u
extins , & ! extinction
amg , & ! alp - gam
apg ! alp + gam
integer (kind=int_kind), parameter :: &
ngmax = 8 ! number of gaussian angles in hemisphere
real (kind=dbl_kind), dimension (ngmax), parameter :: &
gauspt & ! gaussian angles (radians)
= (/ .9894009_dbl_kind, .9445750_dbl_kind, &
.8656312_dbl_kind, .7554044_dbl_kind, &
.6178762_dbl_kind, .4580168_dbl_kind, &
.2816036_dbl_kind, .0950125_dbl_kind/), &
gauswt & ! gaussian weights
= (/ .0271525_dbl_kind, .0622535_dbl_kind, &
.0951585_dbl_kind, .1246290_dbl_kind, &
.1495960_dbl_kind, .1691565_dbl_kind, &
.1826034_dbl_kind, .1894506_dbl_kind/)
integer (kind=int_kind) :: &
ng ! gaussian integration index
real (kind=dbl_kind) :: &
gwt , & ! gaussian weight
swt , & ! sum of weights
trn , & ! layer transmission
rdr , & ! rdir for gaussian integration
tdr , & ! tdir for gaussian integration
smr , & ! accumulator for rdif gaussian integration
smt ! accumulator for tdif gaussian integration
real (kind=dbl_kind) :: &
exp_min ! minimum exponential value
character(len=*),parameter :: subname='(solution_dEdd)'
!-----------------------------------------------------------------------
do k = 0, klevp
trndir(k) = c0
trntdr(k) = c0
trndif(k) = c0
rupdir(k) = c0
rupdif(k) = c0
rdndif(k) = c0
enddo
! initialize top interface of top layer
trndir(0) = c1
trntdr(0) = c1
trndif(0) = c1
rdndif(0) = c0
! mu0 is cosine solar zenith angle above the fresnel level; make
! sure mu0 is large enough for stable and meaningful radiation
! solution: .01 is like sun just touching horizon with its lower edge
mu0 = max(coszen,p01)
! mu0n is cosine solar zenith angle used to compute the layer
! Delta-Eddington solution; it is initially computed to be the
! value below the fresnel level, i.e. the cosine solar zenith
! angle below the fresnel level for the refracted solar beam:
mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
! compute level of fresnel refraction
! if ponded sea ice, fresnel level is the top of the pond.
kfrsnl = 0
! if snow over sea ice or bare sea ice, fresnel level is
! at base of sea ice SSL (and top of the sea ice DL); the
! snow SSL counts for one, then the number of snow layers,
! then the sea ice SSL which also counts for one:
if( srftyp < 2 ) kfrsnl = nslyr + 2
! proceed down one layer at a time; if the total transmission to
! the interface just above a given layer is less than trmin, then no
! Delta-Eddington computation for that layer is done.
! begin main level loop
do k = 0, klev
! initialize all layer apparent optical properties to 0
rdir (k) = c0
rdif_a(k) = c0
rdif_b(k) = c0
tdir (k) = c0
tdif_a(k) = c0
tdif_b(k) = c0
trnlay(k) = c0
! compute next layer Delta-eddington solution only if total transmission
! of radiation to the interface just above the layer exceeds trmin.
if (trntdr(k) > trmin ) then
! calculation over layers with penetrating radiation
tautot = tau(k)
wtot = w0(k)
gtot = g(k)
ftot = gtot*gtot
ts = taus(wtot,ftot,tautot)
ws = omgs(wtot,ftot)
gs = asys(gtot,ftot)
lm = el(ws,gs)
ue = u(ws,gs,lm)
mu0n = mu0nij
! if level k is above fresnel level and the cell is non-pond, use the
! non-refracted beam instead
if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
exp_min = min(exp_argmax,lm*ts)
extins = exp(-exp_min)
ne = n(ue,extins)
! first calculation of rdif, tdif using Delta-Eddington formulas
! rdif_a(k) = (ue+c1)*(ue-c1)*(c1/extins - extins)/ne
rdif_a(k) = (ue**2-c1)*(c1/extins - extins)/ne
tdif_a(k) = c4*ue/ne
! evaluate rdir,tdir for direct beam
exp_min = min(exp_argmax,ts/mu0n)
trnlay(k) = exp(-exp_min)
alp = alpha(ws,mu0n,gs,lm)
gam = agamm(ws,mu0n,gs,lm)
apg = alp + gam
amg = alp - gam
rdir(k) = apg*rdif_a(k) + amg*(tdif_a(k)*trnlay(k) - c1)
tdir(k) = apg*tdif_a(k) + (amg* rdif_a(k)-apg+c1)*trnlay(k)
! recalculate rdif,tdif using direct angular integration over rdir,tdir,
! since Delta-Eddington rdif formula is not well-behaved (it is usually
! biased low and can even be negative); use ngmax angles and gaussian
! integration for most accuracy:
R1 = rdif_a(k) ! use R1 as temporary
T1 = tdif_a(k) ! use T1 as temporary
swt = c0
smr = c0
smt = c0
do ng=1,ngmax
mu = gauspt(ng)
gwt = gauswt(ng)
swt = swt + mu*gwt
exp_min = min(exp_argmax,ts/mu)
trn = exp(-exp_min)
alp = alpha(ws,mu,gs,lm)
gam = agamm(ws,mu,gs,lm)
apg = alp + gam
amg = alp - gam
rdr = apg*R1 + amg*T1*trn - amg
tdr = apg*T1 + amg*R1*trn - apg*trn + trn
smr = smr + mu*rdr*gwt
smt = smt + mu*tdr*gwt
enddo ! ng
rdif_a(k) = smr/swt
tdif_a(k) = smt/swt
! homogeneous layer
rdif_b(k) = rdif_a(k)
tdif_b(k) = tdif_a(k)
! add fresnel layer to top of desired layer if either
! air or snow overlies ice; we ignore refraction in ice
! if a melt pond overlies it:
if( k == kfrsnl ) then
! compute fresnel reflection and transmission amplitudes
! for two polarizations: 1=perpendicular and 2=parallel to
! the plane containing incident, reflected and refracted rays.
R1 = (mu0 - refindx*mu0n) / &
(mu0 + refindx*mu0n)
R2 = (refindx*mu0 - mu0n) / &
(refindx*mu0 + mu0n)
T1 = c2*mu0 / &
(mu0 + refindx*mu0n)
T2 = c2*mu0 / &
(refindx*mu0 + mu0n)
! unpolarized light for direct beam
Rf_dir_a = p5 * (R1*R1 + R2*R2)
Tf_dir_a = p5 * (T1*T1 + T2*T2)*refindx*mu0n/mu0
! precalculated diffuse reflectivities and transmissivities
! for incident radiation above and below fresnel layer, using
! the direct albedos and accounting for complete internal
! reflection from below; precalculated because high order
! number of gaussian points (~256) is required for convergence:
! above
Rf_dif_a = cp063
Tf_dif_a = c1 - Rf_dif_a
! below
Rf_dif_b = cp455
Tf_dif_b = c1 - Rf_dif_b
! the k = kfrsnl layer properties are updated to combined
! the fresnel (refractive) layer, always taken to be above
! the present layer k (i.e. be the top interface):
rintfc = c1 / (c1-Rf_dif_b*rdif_a(k))
tdir(k) = Tf_dir_a*tdir(k) + &
Tf_dir_a*rdir(k) * &
Rf_dif_b*rintfc*tdif_a(k)
rdir(k) = Rf_dir_a + &
Tf_dir_a*rdir(k) * &
rintfc*Tf_dif_b
rdif_a(k) = Rf_dif_a + &
Tf_dif_a*rdif_a(k) * &
rintfc*Tf_dif_b
rdif_b(k) = rdif_b(k) + &
tdif_b(k)*Rf_dif_b * &
rintfc*tdif_a(k)
tdif_a(k) = tdif_a(k)*rintfc*Tf_dif_a
tdif_b(k) = tdif_b(k)*rintfc*Tf_dif_b
! update trnlay to include fresnel transmission
trnlay(k) = Tf_dir_a*trnlay(k)
endif ! k = kfrsnl
endif ! trntdr(k) > trmin
! initialize current layer properties to zero; only if total
! transmission to the top interface of the current layer exceeds the
! minimum, will these values be computed below:
! Calculate the solar beam transmission, total transmission, and
! reflectivity for diffuse radiation from below at interface k,
! the top of the current layer k:
!
! layers interface
!
! --------------------- k-1
! k-1
! --------------------- k
! k
! ---------------------
! For k = klevp
! note that we ignore refraction between sea ice and underlying ocean:
!
! layers interface
!
! --------------------- k-1
! k-1
! --------------------- k
! \\\\\\\ ocean \\\\\\\
trndir(k+1) = trndir(k)*trnlay(k)
refkm1 = c1/(c1 - rdndif(k)*rdif_a(k))
tdrrdir = trndir(k)*rdir(k)
tdndif = trntdr(k) - trndir(k)
trntdr(k+1) = trndir(k)*tdir(k) + &
(tdndif + tdrrdir*rdndif(k))*refkm1*tdif_a(k)
rdndif(k+1) = rdif_b(k) + &
(tdif_b(k)*rdndif(k)*refkm1*tdif_a(k))
trndif(k+1) = trndif(k)*refkm1*tdif_a(k)
enddo ! k end main level loop
! compute reflectivity to direct and diffuse radiation for layers
! below by adding succesive layers starting from the underlying
! ocean and working upwards:
!
! layers interface
!
! --------------------- k
! k
! --------------------- k+1
! k+1
! ---------------------
rupdir(klevp) = albodr
rupdif(klevp) = albodf
do k=klev,0,-1
! interface scattering
refkp1 = c1/( c1 - rdif_b(k)*rupdif(k+1))
! dir from top layer plus exp tran ref from lower layer, interface
! scattered and tran thru top layer from below, plus diff tran ref
! from lower layer with interface scattering tran thru top from below
rupdir(k) = rdir(k) &
+ ( trnlay(k) *rupdir(k+1) &
+ (tdir(k)-trnlay(k))*rupdif(k+1))*refkp1*tdif_b(k)
! dif from top layer from above, plus dif tran upwards reflected and
! interface scattered which tran top from below
rupdif(k) = rdif_a(k) + tdif_a(k)*rupdif(k+1)*refkp1*tdif_b(k)
enddo ! k
end subroutine solution_dEdd
!=======================================================================
!
! Set snow horizontal coverage, density and grain radius diagnostically
! for the Delta-Eddington solar radiation method.
!
! author: Bruce P. Briegleb, NCAR
! 2013: E Hunke merged with NCAR version
subroutine shortwave_dEdd_set_snow(nslyr, R_snw, &
dT_mlt, rsnw_mlt, &
aice, vsno, &
Tsfc, fs, &
hs0, hs, &
rhosnw, rsnw)
integer (kind=int_kind), intent(in) :: &
nslyr ! number of snow layers
real (kind=dbl_kind), intent(in) :: &
R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo
dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C)
rsnw_mlt ! maximum melting snow grain radius (10^-6 m)
real (kind=dbl_kind), intent(in) :: &
aice , & ! concentration of ice
vsno , & ! volume of snow
Tsfc , & ! surface temperature
hs0 ! snow depth for transition to bare sea ice (m)
real (kind=dbl_kind), intent(inout) :: &
fs , & ! horizontal coverage of snow
hs ! snow depth
real (kind=dbl_kind), dimension (:), intent(out) :: &
rhosnw , & ! density in snow layer (kg/m3)
rsnw ! grain radius in snow layer (micro-meters)
! local variables
integer (kind=int_kind) :: &
ks ! snow vertical index
real (kind=dbl_kind) :: &
fT , & ! piecewise linear function of surface temperature
dTs , & ! difference of Tsfc and Timelt
rsnw_nm ! actual used nonmelt snow grain radius (micro-meters)
real (kind=dbl_kind), parameter :: &
! units for the following are 1.e-6 m (micro-meters)
rsnw_fresh = 100._dbl_kind, & ! freshly-fallen snow grain radius
rsnw_nonmelt = 500._dbl_kind, & ! nonmelt snow grain radius
rsnw_sig = 250._dbl_kind ! assumed sigma for snow grain radius
character(len=*),parameter :: subname='(shortwave_dEdd_set_snow)'
!-----------------------------------------------------------------------
! set snow horizontal fraction
hs = vsno / aice
if (hs >= hs_min) then
fs = c1
if (hs0 > puny) fs = min(hs/hs0, c1)
endif
! bare ice, temperature dependence
dTs = Timelt - Tsfc
fT = -min(dTs/dT_mlt-c1,c0)
! tune nonmelt snow grain radius if desired: note that
! the sign is negative so that if R_snw is 1, then the
! snow grain radius is reduced and thus albedo increased.
rsnw_nm = rsnw_nonmelt - R_snw*rsnw_sig
rsnw_nm = max(rsnw_nm, rsnw_fresh)
rsnw_nm = min(rsnw_nm, rsnw_mlt)
do ks = 1, nslyr
! snow density ccsm3 constant value
rhosnw(ks) = rhos
! snow grain radius between rsnw_nonmelt and rsnw_mlt
rsnw(ks) = rsnw_nm + (rsnw_mlt-rsnw_nm)*fT
rsnw(ks) = max(rsnw(ks), rsnw_fresh)
rsnw(ks) = min(rsnw(ks), rsnw_mlt)
enddo ! ks
end subroutine shortwave_dEdd_set_snow
!=======================================================================
!
! Set pond fraction and depth diagnostically for
! the Delta-Eddington solar radiation method.
!
! author: Bruce P. Briegleb, NCAR
! 2013: E Hunke merged with NCAR version
subroutine shortwave_dEdd_set_pond(Tsfc, &
fs, fp, &
hp)
real (kind=dbl_kind), intent(in) :: &
Tsfc , & ! surface temperature
fs ! horizontal coverage of snow
real (kind=dbl_kind), intent(out) :: &
fp , & ! pond fractional coverage (0 to 1)
hp ! pond depth (m)
! local variables
real (kind=dbl_kind) :: &
fT , & ! piecewise linear function of surface temperature
dTs ! difference of Tsfc and Timelt
real (kind=dbl_kind), parameter :: &
dT_pnd = c1 ! change in temp for pond fraction and depth
character(len=*),parameter :: subname='(shortwave_dEdd_set_pond)'
!-----------------------------------------------------------------------
! bare ice, temperature dependence
dTs = Timelt - Tsfc
fT = -min(dTs/dT_pnd-c1,c0)
! pond
fp = 0.3_dbl_kind*fT*(c1-fs)
hp = 0.3_dbl_kind*fT*(c1-fs)
end subroutine shortwave_dEdd_set_pond
! End Delta-Eddington shortwave method
!=======================================================================
!
! authors Nicole Jeffery, LANL
subroutine compute_shortwave_trcr(nslyr, &
bgcN, zaero, &
trcrn_bgcsw, &
sw_grid, hin, &
hbri, &
nilyr, nblyr, &
i_grid, &
skl_bgc, z_tracers )
integer (kind=int_kind), intent(in) :: &
nslyr ! number of snow layers
integer (kind=int_kind), intent(in) :: &
nblyr , & ! number of bio layers
nilyr ! number of ice layers
real (kind=dbl_kind), dimension (:), intent(in) :: &
bgcN , & ! Nit tracer
zaero ! zaero tracer
real (kind=dbl_kind), dimension (:), intent(out):: &
trcrn_bgcsw ! ice on shortwave grid tracers
real (kind=dbl_kind), dimension (:), intent(in) :: &
sw_grid , & !
i_grid ! CICE bio grid
real(kind=dbl_kind), intent(in) :: &
hin , & ! CICE ice thickness
hbri ! brine height
logical (kind=log_kind), intent(in) :: &
skl_bgc , & ! skeletal layer bgc
z_tracers ! zbgc
! local variables
integer (kind=int_kind) :: k, n, nn
real (kind=dbl_kind), dimension (ntrcr+2) :: &
trtmp0, & ! temporary, remapped tracers
trtmp
real (kind=dbl_kind), dimension (nilyr+1):: &
icegrid ! correct for large ice surface layers
real (kind=dbl_kind):: &
top_conc ! 1% (min_bgc) of surface concentration
! when hin > hbri: just used in sw calculation
character(len=*),parameter :: subname='(compute_shortwave_trcr)'
!-----------------------------------------------------------------
! Compute aerosols and algal chlorophyll on shortwave grid
!-----------------------------------------------------------------
trtmp0(:) = c0
trtmp(:) = c0
trcrn_bgcsw(:) = c0
do k = 1,nilyr+1
icegrid(k) = sw_grid(k)
enddo
if (sw_grid(1)*hin*c2 > hi_ssl) then
icegrid(1) = hi_ssl/c2/hin
endif
if (z_tracers) then
if (tr_bgc_N) then
if (size(bgcN) < n_algae*(nblyr+3)) then
call icepack_warnings_add(subname//' ERROR: size(bgcN) too small')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
do k = 1, nblyr+1
do n = 1, n_algae
trtmp0(nt_bgc_N(1) + k-1) = trtmp0(nt_bgc_N(1) + k-1) + &
R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1 + k-1)
enddo ! n
enddo ! k
top_conc = trtmp0(nt_bgc_N(1))*min_bgc
call remap_zbgc (nilyr+1, &
nt_bgc_N(1), &
trtmp0(1:ntrcr ), &
trtmp (1:ntrcr+2), &
1, nblyr+1, &
hin, hbri, &
icegrid(1:nilyr+1), &
i_grid(1:nblyr+1), top_conc )
if (icepack_warnings_aborted(subname)) return
do k = 1, nilyr+1
trcrn_bgcsw(nlt_chl_sw+nslyr+k) = trtmp(nt_bgc_N(1) + k-1)
enddo ! k
do n = 1, n_algae ! snow contribution
trcrn_bgcsw(nlt_chl_sw)= trcrn_bgcsw(nlt_chl_sw) &
+ R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+1)
! snow surface layer
trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) = &
trcrn_bgcsw(nlt_chl_sw+1:nlt_chl_sw+nslyr) &
+ R_chl2N(n)*F_abs_chl(n)*bgcN(nt_bgc_N(n)-nt_bgc_N(1)+1+nblyr+2)
! only 1 snow layer in zaero
enddo ! n
endif ! tr_bgc_N
if (tr_zaero) then
if (size(zaero) < n_zaero*(nblyr+3)) then
call icepack_warnings_add(subname//' ERROR: size(zaero) too small')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
do n = 1, n_zaero
trtmp0(:) = c0
trtmp(:) = c0
do k = 1, nblyr+1
trtmp0(nt_zaero(n) + k-1) = zaero(nt_zaero(n)-nt_zaero(1)+1+k-1)
enddo
top_conc = trtmp0(nt_zaero(n))*min_bgc
call remap_zbgc (nilyr+1, &
nt_zaero(n), &
trtmp0(1:ntrcr ), &
trtmp (1:ntrcr+2), &
1, nblyr+1, &
hin, hbri, &
icegrid(1:nilyr+1), &
i_grid(1:nblyr+1), top_conc )
if (icepack_warnings_aborted(subname)) return
do k = 1,nilyr+1
trcrn_bgcsw(nlt_zaero_sw(n)+nslyr+k) = trtmp(nt_zaero(n) + k-1)
enddo
trcrn_bgcsw(nlt_zaero_sw(n))= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+1) !snow ssl
trcrn_bgcsw(nlt_zaero_sw(n)+1:nlt_zaero_sw(n)+nslyr)= zaero(nt_zaero(n)-nt_zaero(1)+1+nblyr+2)
enddo ! n
endif ! tr_zaero
elseif (skl_bgc) then
do nn = 1,n_algae
trcrn_bgcsw(nbtrcr_sw) = trcrn_bgcsw(nbtrcr_sw) &
+ F_abs_chl(nn)*R_chl2N(nn) &
* bgcN(nt_bgc_N(nn)-nt_bgc_N(1)+1)*sk_l/hin &
* real(nilyr,kind=dbl_kind)
enddo
endif
end subroutine compute_shortwave_trcr
!=======================================================================
!autodocument_start icepack_prep_radiation
! Scales radiation fields computed on the previous time step.
!
! authors: Elizabeth Hunke, LANL
subroutine icepack_prep_radiation (ncat, nilyr, nslyr, &
aice, aicen, &
swvdr, swvdf, &
swidr, swidf, &
alvdr_ai, alvdf_ai, &
alidr_ai, alidf_ai, &
scale_factor, &
fswsfcn, fswintn, &
fswthrun, &
fswthrun_vdr, &
fswthrun_vdf, &
fswthrun_idr, &
fswthrun_idf, &
fswpenln, &
Sswabsn, Iswabsn)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of ice thickness categories
nilyr , & ! number of ice layers
nslyr ! number of snow layers
real (kind=dbl_kind), intent(in) :: &
aice , & ! ice area fraction
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf , & ! sw down, near IR, diffuse (W/m^2)
! grid-box-mean albedos aggregated over categories (if calc_Tsfc)
alvdr_ai , & ! visible, direct (fraction)
alidr_ai , & ! near-ir, direct (fraction)
alvdf_ai , & ! visible, diffuse (fraction)
alidf_ai ! near-ir, diffuse (fraction)
real (kind=dbl_kind), dimension(:), intent(in) :: &
aicen ! ice area fraction in each category
real (kind=dbl_kind), intent(inout) :: &
scale_factor ! shortwave scaling factor, ratio new:old
real (kind=dbl_kind), dimension(:), intent(inout) :: &
fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
fswthrun ! SW through ice to ocean (W/m^2)
real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
fswpenln , & ! visible SW entering ice layers (W m-2)
Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
Sswabsn ! SW radiation absorbed in snow layers (W m-2)
!autodocument_end
! local variables
integer (kind=int_kind) :: &
k , & ! vertical index
n ! thickness category index
real (kind=dbl_kind) :: netsw
character(len=*),parameter :: subname='(icepack_prep_radiation)'
!-----------------------------------------------------------------
! Compute netsw scaling factor (new netsw / old netsw)
!-----------------------------------------------------------------
if (aice > c0 .and. scale_factor > puny) then
netsw = swvdr*(c1 - alvdr_ai) &
+ swvdf*(c1 - alvdf_ai) &
+ swidr*(c1 - alidr_ai) &
+ swidf*(c1 - alidf_ai)
scale_factor = netsw / scale_factor
else
scale_factor = c1
endif
do n = 1, ncat
if (aicen(n) > puny) then
!-----------------------------------------------------------------
! Scale absorbed solar radiation for change in net shortwave
!-----------------------------------------------------------------
fswsfcn(n) = scale_factor*fswsfcn (n)
fswintn(n) = scale_factor*fswintn (n)
fswthrun(n) = scale_factor*fswthrun(n)
if (present(fswthrun_vdr)) fswthrun_vdr(n) = scale_factor*fswthrun_vdr(n)
if (present(fswthrun_vdf)) fswthrun_vdf(n) = scale_factor*fswthrun_vdf(n)
if (present(fswthrun_idr)) fswthrun_idr(n) = scale_factor*fswthrun_idr(n)
if (present(fswthrun_idf)) fswthrun_idf(n) = scale_factor*fswthrun_idf(n)
do k = 1,nilyr+1
fswpenln(k,n) = scale_factor*fswpenln(k,n)
enddo !k
do k=1,nslyr
Sswabsn(k,n) = scale_factor*Sswabsn(k,n)
enddo
do k=1,nilyr
Iswabsn(k,n) = scale_factor*Iswabsn(k,n)
enddo
endif
enddo ! ncat
end subroutine icepack_prep_radiation
!=======================================================================
!autodocument_start icepack_step_radiation
! Computes radiation fields
!
! authors: William H. Lipscomb, LANL
! David Bailey, NCAR
! Elizabeth C. Hunke, LANL
subroutine icepack_step_radiation (dt, ncat, &
nblyr, &
nilyr, nslyr, &
dEdd_algae, &
swgrid, igrid, &
fbri, &
aicen, vicen, &
vsnon, Tsfcn, &
alvln, apndn, &
hpndn, ipndn, &
aeron, &
bgcNn, zaeron, &
trcrn_bgcsw, &
TLAT, TLON, &
calendar_type, &
days_per_year, &
nextsw_cday, &
yday, sec, &
kaer_tab, waer_tab, &
gaer_tab, &
kaer_bc_tab, &
waer_bc_tab, &
gaer_bc_tab, &
bcenh, &
modal_aero, &
swvdr, swvdf, &
swidr, swidf, &
coszen, fsnow, &
alvdrn, alvdfn, &
alidrn, alidfn, &
fswsfcn, fswintn, &
fswthrun, &
fswthrun_vdr, &
fswthrun_vdf, &
fswthrun_idr, &
fswthrun_idf, &
fswpenln, &
Sswabsn, Iswabsn, &
albicen, albsnon, &
albpndn, apeffn, &
snowfracn, &
dhsn, ffracn, &
l_print_point, &
initonly)
integer (kind=int_kind), intent(in) :: &
ncat , & ! number of ice thickness categories
nilyr , & ! number of ice layers
nslyr , & ! number of snow layers
nblyr ! number of bgc layers
real (kind=dbl_kind), intent(in) :: &
dt , & ! time step (s)
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf , & ! sw down, near IR, diffuse (W/m^2)
fsnow , & ! snowfall rate (kg/m^2 s)
TLAT, TLON ! latitude and longitude (radian)
character (len=char_len), intent(in) :: &
calendar_type ! differentiates Gregorian from other calendars
integer (kind=int_kind), intent(in) :: &
days_per_year, & ! number of days in one year
sec ! elapsed seconds into date
real (kind=dbl_kind), intent(in) :: &
nextsw_cday , & ! julian day of next shortwave calculation
yday ! day of the year
real (kind=dbl_kind), intent(inout) :: &
coszen ! cosine solar zenith angle, < 0 for sun below horizon
real (kind=dbl_kind), dimension (:), intent(in) :: &
igrid ! biology vertical interface points
real (kind=dbl_kind), dimension (:), intent(in) :: &
swgrid ! grid for ice tracers used in dEdd scheme
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_tab, & ! aerosol single scatter albedo (fraction)
gaer_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:), intent(in) :: &
kaer_bc_tab, & ! aerosol mass extinction cross section (m2/kg)
waer_bc_tab, & ! aerosol single scatter albedo (fraction)
gaer_bc_tab ! aerosol asymmetry parameter (cos(theta))
real (kind=dbl_kind), dimension(:,:,:), intent(in) :: &
bcenh
real (kind=dbl_kind), dimension(:), intent(in) :: &
aicen , & ! ice area fraction in each category
vicen , & ! ice volume in each category (m)
vsnon , & ! snow volume in each category (m)
Tsfcn , & ! surface temperature (deg C)
alvln , & ! level-ice area fraction
apndn , & ! pond area fraction
hpndn , & ! pond depth (m)
ipndn , & ! pond refrozen lid thickness (m)
fbri ! brine fraction
real(kind=dbl_kind), dimension(:,:), intent(in) :: &
aeron , & ! aerosols (kg/m^3)
bgcNn , & ! bgc Nit tracers
zaeron ! bgcz aero tracers
real(kind=dbl_kind), dimension(:,:), intent(inout) :: &
trcrn_bgcsw ! zaerosols (kg/m^3) and chla (mg/m^3)
real (kind=dbl_kind), dimension(:), intent(inout) :: &
alvdrn , & ! visible, direct albedo (fraction)
alidrn , & ! near-ir, direct (fraction)
alvdfn , & ! visible, diffuse (fraction)
alidfn , & ! near-ir, diffuse (fraction)
fswsfcn , & ! SW absorbed at ice/snow surface (W m-2)
fswintn , & ! SW absorbed in ice interior, below surface (W m-2)
fswthrun , & ! SW through ice to ocean (W/m^2)
snowfracn , & ! snow fraction on each category
dhsn , & ! depth difference for snow on sea ice and pond ice
ffracn , & ! fraction of fsurfn used to melt ipond
! albedo components for history
albicen , & ! bare ice
albsnon , & ! snow
albpndn , & ! pond
apeffn ! effective pond area used for radiation calculation
real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
fswpenln , & ! visible SW entering ice layers (W m-2)
Iswabsn , & ! SW radiation absorbed in ice layers (W m-2)
Sswabsn ! SW radiation absorbed in snow layers (W m-2)
logical (kind=log_kind), intent(in) :: &
l_print_point, & ! flag for printing diagnostics
dEdd_algae , & ! .true. use prognostic chla in dEdd
modal_aero ! .true. use modal aerosol optical treatment
logical (kind=log_kind), optional :: &
initonly ! flag to indicate init only, default is false
!autodocument_end
! local variables
integer (kind=int_kind) :: &
n ! thickness category index
logical (kind=log_kind) :: &
linitonly ! local flag for initonly
real(kind=dbl_kind) :: &
hin, & ! Ice thickness (m)
hbri ! brine thickness (m)
real (kind=dbl_kind), dimension(:), allocatable :: &
l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2)
l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2)
l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2)
l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2)
character(len=*),parameter :: subname='(icepack_step_radiation)'
allocate(l_fswthrun_vdr(ncat))
allocate(l_fswthrun_vdf(ncat))
allocate(l_fswthrun_idr(ncat))
allocate(l_fswthrun_idf(ncat))
hin = c0
hbri = c0
linitonly = .false.
if (present(initonly)) then
linitonly = initonly
endif
! Initialize
do n = 1, ncat
alvdrn (n) = c0
alidrn (n) = c0
alvdfn (n) = c0
alidfn (n) = c0
fswsfcn (n) = c0
fswintn (n) = c0
fswthrun(n) = c0
enddo ! ncat
fswpenln (:,:) = c0
Iswabsn (:,:) = c0
Sswabsn (:,:) = c0
trcrn_bgcsw(:,:) = c0
! Interpolate z-shortwave tracers to shortwave grid
if (dEdd_algae) then
do n = 1, ncat
if (aicen(n) .gt. puny) then
hin = vicen(n)/aicen(n)
hbri= fbri(n)*hin
call compute_shortwave_trcr(nslyr, &
bgcNn(:,n), &
zaeron(:,n), &
trcrn_bgcsw(:,n), &
swgrid, hin, &
hbri, &
nilyr, nblyr, &
igrid, &
skl_bgc, z_tracers )
if (icepack_warnings_aborted(subname)) return
endif
enddo
endif
if (calc_Tsfc) then
if (trim(shortwave) == 'dEdd') then ! delta Eddington
call run_dEdd(dt, ncat, &
dEdd_algae, &
nilyr, nslyr, &
aicen, vicen, &
vsnon, Tsfcn, &
alvln, apndn, &
hpndn, ipndn, &
aeron, kalg, &
trcrn_bgcsw, &
heat_capacity, &
TLAT, TLON, &
calendar_type,days_per_year, &
nextsw_cday, yday, &
sec, R_ice, &
R_pnd, R_snw, &
dT_mlt, rsnw_mlt, &
hs0, hs1, &
hp1, pndaspect, &
kaer_tab, waer_tab, &
gaer_tab, &
kaer_bc_tab, &
waer_bc_tab, &
gaer_bc_tab, &
bcenh, &
modal_aero, &
swvdr, swvdf, &
swidr, swidf, &
coszen, fsnow, &
alvdrn, alvdfn, &
alidrn, alidfn, &
fswsfcn, fswintn, &
fswthrun=fswthrun, &
fswthrun_vdr=l_fswthrun_vdr, &
fswthrun_vdf=l_fswthrun_vdf, &
fswthrun_idr=l_fswthrun_idr, &
fswthrun_idf=l_fswthrun_idf, &
fswpenln=fswpenln, &
Sswabsn=Sswabsn, &
Iswabsn=Iswabsn, &
albicen=albicen, &
albsnon=albsnon, &
albpndn=albpndn, &
apeffn=apeffn, &
snowfracn=snowfracn, &
dhsn=dhsn, &
ffracn=ffracn, &
l_print_point=l_print_point, &
initonly=linitonly)
if (icepack_warnings_aborted(subname)) return
elseif (trim(shortwave) == 'ccsm3') then
call shortwave_ccsm3(aicen, vicen, &
vsnon, &
Tsfcn, &
swvdr, swvdf, &
swidr, swidf, &
heat_capacity, &
albedo_type, &
albicev, albicei, &
albsnowv, albsnowi, &
ahmax, &
alvdrn, alidrn, &
alvdfn, alidfn, &
fswsfcn, fswintn, &
fswthru=fswthrun, &
fswthru_vdr=l_fswthrun_vdr,&
fswthru_vdf=l_fswthrun_vdf,&
fswthru_idr=l_fswthrun_idr,&
fswthru_idf=l_fswthrun_idf,&
fswpenl=fswpenln, &
Iswabs=Iswabsn, &
Sswabs=Sswabsn, &
albin=albicen, &
albsn=albsnon, &
coszen=coszen, &
ncat=ncat, &
nilyr=nilyr)
if (icepack_warnings_aborted(subname)) return
else
call icepack_warnings_add(subname//' ERROR: shortwave '//trim(shortwave)//' unknown')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif ! shortwave
else ! .not. calc_Tsfc
! Calculate effective pond area for HadGEM
if (tr_pond_topo) then
do n = 1, ncat
apeffn(n) = c0
if (aicen(n) > puny) then
! Lid effective if thicker than hp1
if (apndn(n)*aicen(n) > puny .and. ipndn(n) < hp1) then
apeffn(n) = apndn(n)
else
apeffn(n) = c0
endif
if (apndn(n) < puny) apeffn(n) = c0
endif
enddo ! ncat
endif ! tr_pond_topo
! Initialize for safety
do n = 1, ncat
alvdrn(n) = c0
alidrn(n) = c0
alvdfn(n) = c0
alidfn(n) = c0
fswsfcn(n) = c0
fswintn(n) = c0
fswthrun(n) = c0
enddo ! ncat
Iswabsn(:,:) = c0
Sswabsn(:,:) = c0
endif ! calc_Tsfc
if (present(fswthrun_vdr)) fswthrun_vdr = l_fswthrun_vdr
if (present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf
if (present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr
if (present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf
deallocate(l_fswthrun_vdr)
deallocate(l_fswthrun_vdf)
deallocate(l_fswthrun_idr)
deallocate(l_fswthrun_idf)
end subroutine icepack_step_radiation
! Delta-Eddington solution expressions
!=======================================================================
real(kind=dbl_kind) function alpha(w,uu,gg,e)
real(kind=dbl_kind), intent(in) :: w, uu, gg, e
alpha = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
end function alpha
!=======================================================================
real(kind=dbl_kind) function agamm(w,uu,gg,e)
real(kind=dbl_kind), intent(in) :: w, uu, gg, e
agamm = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
end function agamm
!=======================================================================
real(kind=dbl_kind) function n(uu,et)
real(kind=dbl_kind), intent(in) :: uu, et
n = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
end function n
!=======================================================================
real(kind=dbl_kind) function u(w,gg,e)
real(kind=dbl_kind), intent(in) :: w, gg, e
u = c1p5*(c1 - w*gg)/e
end function u
!=======================================================================
real(kind=dbl_kind) function el(w,gg)
real(kind=dbl_kind), intent(in) :: w, gg
el = sqrt(c3*(c1-w)*(c1 - w*gg))
end function el
!=======================================================================
real(kind=dbl_kind) function taus(w,f,t)
real(kind=dbl_kind), intent(in) :: w, f, t
taus = (c1 - w*f)*t
end function taus
!=======================================================================
real(kind=dbl_kind) function omgs(w,f)
real(kind=dbl_kind), intent(in) :: w, f
omgs = (c1 - f)*w/(c1 - w*f)
end function omgs
!=======================================================================
real(kind=dbl_kind) function asys(gg,f)
real(kind=dbl_kind), intent(in) :: gg, f
asys = (gg - f)/(c1 - f)
end function asys
!=======================================================================
end module icepack_shortwave
!=======================================================================