#include <misc.h>
#include <preproc.h>
module SNICARMod 8,4
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: SNICARMod
!
! !DESCRIPTION:
! Calculate albedo of snow containing impurities
! and the evolution of snow effective radius
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_sys_mod , only : shr_sys_flush
use clm_varctl , only : iulog
use shr_const_mod , only : SHR_CONST_RHOICE
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
public :: SNICAR_RT ! Snow albedo and vertically-resolved solar absorption
public :: SnowAge_grain ! Snow effective grain size evolution
public :: SnowAge_init ! Initial read in of snow-aging file
public :: SnowOptics_init ! Initial read in of snow-optics file
!
! !PUBLIC DATA MEMBERS:
real(r8), public, parameter :: snw_rds_min = 54.526_r8 ! minimum allowed snow effective radius (also "fresh snow" value) [microns]
integer, public, parameter :: sno_nbr_aer = 8 ! number of aerosol species in snowpack (indices described above) [nbr]
logical, public, parameter :: DO_SNO_OC = .false. ! parameter to include organic carbon (OC) in snowpack radiative calculations
logical, public, parameter :: DO_SNO_AER = .true. ! parameter to include aerosols in snowpack radiative calculations
real(r8), public, parameter :: scvng_fct_mlt_bcphi = 0.20_r8 ! scavenging factor for hydrophillic BC inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_bcpho = 0.03_r8 ! scavenging factor for hydrophobic BC inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_ocphi = 0.20_r8 ! scavenging factor for hydrophillic OC inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_ocpho = 0.03_r8 ! scavenging factor for hydrophobic OC inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_dst1 = 0.02_r8 ! scavenging factor for dust species 1 inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_dst2 = 0.02_r8 ! scavenging factor for dust species 2 inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_dst3 = 0.01_r8 ! scavenging factor for dust species 3 inclusion in meltwater [frc]
real(r8), public, parameter :: scvng_fct_mlt_dst4 = 0.01_r8 ! scavenging factor for dust species 4 inclusion in meltwater [frc]
! !PRIVATE MEMBER FUNCTIONS:
!
! !PRIVATE DATA MEMBERS:
! Aerosol species indices:
! 1= hydrophillic black carbon
! 2= hydrophobic black carbon
! 3= hydrophilic organic carbon
! 4= hydrophobic organic carbon
! 5= dust species 1
! 6= dust species 2
! 7= dust species 3
! 8= dust species 4
integer, parameter :: numrad_snw = 5 ! number of spectral bands used in snow model [nbr]
integer, parameter :: nir_bnd_bgn = 2 ! first band index in near-IR spectrum [idx]
integer, parameter :: nir_bnd_end = 5 ! ending near-IR band index [idx]
integer, parameter :: idx_Mie_snw_mx = 1471 ! number of effective radius indices used in Mie lookup table [idx]
integer, parameter :: idx_T_max = 11 ! maxiumum temperature index used in aging lookup table [idx]
integer, parameter :: idx_T_min = 1 ! minimum temperature index used in aging lookup table [idx]
integer, parameter :: idx_Tgrd_max = 31 ! maxiumum temperature gradient index used in aging lookup table [idx]
integer, parameter :: idx_Tgrd_min = 1 ! minimum temperature gradient index used in aging lookup table [idx]
integer, parameter :: idx_rhos_max = 8 ! maxiumum snow density index used in aging lookup table [idx]
integer, parameter :: idx_rhos_min = 1 ! minimum snow density index used in aging lookup table [idx]
integer, parameter :: snw_rds_max_tbl = 1500 ! maximum effective radius defined in Mie lookup table [microns]
integer, parameter :: snw_rds_min_tbl = 30 ! minimium effective radius defined in Mie lookup table [microns]
real(r8), parameter :: snw_rds_max = 1500._r8 ! maximum allowed snow effective radius [microns]
real(r8), parameter :: snw_rds_refrz = 1000._r8 ! effective radius of re-frozen snow [microns]
real(r8), parameter :: min_snw = 1.0E-30_r8 ! minimum snow mass required for SNICAR RT calculation [kg m-2]
!real(r8), parameter :: C1_liq_Brun89 = 1.28E-17_r8 ! constant for liquid water grain growth [m3 s-1], from Brun89
real(r8), parameter :: C1_liq_Brun89 = 0._r8 ! constant for liquid water grain growth [m3 s-1], from Brun89: zeroed to accomodate dry snow aging
real(r8), parameter :: C2_liq_Brun89 = 4.22E-13_r8 ! constant for liquid water grain growth [m3 s-1], from Brun89: corrected for LWC in units of percent
real(r8), parameter :: tim_cns_bc_rmv = 2.2E-8_r8 ! time constant for removal of BC in snow on sea-ice [s-1] (50% mass removal/year)
real(r8), parameter :: tim_cns_oc_rmv = 2.2E-8_r8 ! time constant for removal of OC in snow on sea-ice [s-1] (50% mass removal/year)
real(r8), parameter :: tim_cns_dst_rmv = 2.2E-8_r8 ! time constant for removal of dust in snow on sea-ice [s-1] (50% mass removal/year)
! scaling of the snow aging rate (tuning option):
logical :: flg_snoage_scl = .false. ! flag for scaling the snow aging rate by some arbitrary factor
real(r8), parameter :: xdrdt = 1.0_r8 ! arbitrary factor applied to snow aging rate
! snow and aerosol Mie parameters:
! (arrays declared here, but are set in iniTimeConst)
! (idx_Mie_snw_mx is number of snow radii with defined parameters (i.e. from 30um to 1500um))
! direct-beam weighted ice optical properties
real(r8) :: ss_alb_snw_drc(idx_Mie_snw_mx,numrad_snw)
real(r8) :: asm_prm_snw_drc(idx_Mie_snw_mx,numrad_snw)
real(r8) :: ext_cff_mss_snw_drc(idx_Mie_snw_mx,numrad_snw)
! diffuse radiation weighted ice optical properties
real(r8) :: ss_alb_snw_dfs(idx_Mie_snw_mx,numrad_snw)
real(r8) :: asm_prm_snw_dfs(idx_Mie_snw_mx,numrad_snw)
real(r8) :: ext_cff_mss_snw_dfs(idx_Mie_snw_mx,numrad_snw)
! hydrophiliic BC
real(r8) :: ss_alb_bc1(1,numrad_snw)
real(r8) :: asm_prm_bc1(1,numrad_snw)
real(r8) :: ext_cff_mss_bc1(1,numrad_snw)
! hydrophobic BC
real(r8) :: ss_alb_bc2(1,numrad_snw)
real(r8) :: asm_prm_bc2(1,numrad_snw)
real(r8) :: ext_cff_mss_bc2(1,numrad_snw)
! hydrophobic OC
real(r8) :: ss_alb_oc1(1,numrad_snw)
real(r8) :: asm_prm_oc1(1,numrad_snw)
real(r8) :: ext_cff_mss_oc1(1,numrad_snw)
! hydrophilic OC
real(r8) :: ss_alb_oc2(1,numrad_snw)
real(r8) :: asm_prm_oc2(1,numrad_snw)
real(r8) :: ext_cff_mss_oc2(1,numrad_snw)
! dust species 1:
real(r8) :: ss_alb_dst1(1,numrad_snw)
real(r8) :: asm_prm_dst1(1,numrad_snw)
real(r8) :: ext_cff_mss_dst1(1,numrad_snw)
! dust species 2:
real(r8) :: ss_alb_dst2(1,numrad_snw)
real(r8) :: asm_prm_dst2(1,numrad_snw)
real(r8) :: ext_cff_mss_dst2(1,numrad_snw)
! dust species 3:
real(r8) :: ss_alb_dst3(1,numrad_snw)
real(r8) :: asm_prm_dst3(1,numrad_snw)
real(r8) :: ext_cff_mss_dst3(1,numrad_snw)
! dust species 4:
real(r8) :: ss_alb_dst4(1,numrad_snw)
real(r8) :: asm_prm_dst4(1,numrad_snw)
real(r8) :: ext_cff_mss_dst4(1,numrad_snw)
! best-fit parameters for snow aging defined over:
! 11 temperatures from 225 to 273 K
! 31 temperature gradients from 0 to 300 K/m
! 8 snow densities from 0 to 350 kg/m3
! (arrays declared here, but are set in iniTimeConst)
real(r8) :: snowage_tau(idx_rhos_max,idx_Tgrd_max,idx_T_max)
real(r8) :: snowage_kappa(idx_rhos_max,idx_Tgrd_max,idx_T_max)
real(r8) :: snowage_drdt0(idx_rhos_max,idx_Tgrd_max,idx_T_max)
!
! !REVISION HISTORY:
! Created by Mark Flanner
!
!EOP
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SNICAR_RT
!
!
! !CALLED FROM:
! subroutine SurfaceAlbedo in module SurfaceAlbedoMod (CLM)
! subroutine albice (CSIM)
!
! !REVISION HISTORY:
! Author: Mark Flanner
!
! !INTERFACE:
subroutine SNICAR_RT (flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, & 10,18
coszen, flg_slr_in, h2osno_liq, h2osno_ice, snw_rds, &
mss_cnc_aer_in, albsfc, albout, flx_abs)
!
! !DESCRIPTION:
! Determine reflectance of, and vertically-resolved solar absorption in,
! snow with impurities.
!
! Original references on physical models of snow reflectance include:
! Wiscombe and Warren [1980] and Warren and Wiscombe [1980],
! Journal of Atmospheric Sciences, 37,
!
! The multi-layer solution for multiple-scattering used here is from:
! Toon et al. [1989], Rapid calculation of radiative heating rates
! and photodissociation rates in inhomogeneous multiple scattering atmospheres,
! J. Geophys. Res., 94, D13, 16287-16301
!
! The implementation of the SNICAR model in CLM/CSIM is described in:
! Flanner, M., C. Zender, J. Randerson, and P. Rasch [2007],
! Present-day climate forcing and response from black carbon in snow,
! J. Geophys. Res., 112, D11202, doi: 10.1029/2006JD008003
! !USES:
use clmtype
use clm_varpar , only : nlevsno, numrad
use clm_time_manager , only : get_nstep
use abortutils , only : endrun
use shr_const_mod , only : SHR_CONST_PI
!
! !ARGUMENTS:
implicit none
integer , intent(in) :: flg_snw_ice ! flag: =1 when called from CLM, =2 when called from CSIM
integer , intent(in) :: lbc, ubc ! column index bounds [unitless]
integer , intent(in) :: num_nourbanc ! number of columns in non-urban filter
integer , intent(in) :: filter_nourbanc(ubc-lbc+1) ! column filter for non-urban points
real(r8), intent(in) :: coszen(lbc:ubc) ! cosine of solar zenith angle for next time step (col) [unitless]
integer , intent(in) :: flg_slr_in ! flag: =1 for direct-beam incident flux, =2 for diffuse incident flux
real(r8), intent(in) :: h2osno_liq(lbc:ubc,-nlevsno+1:0) ! liquid water content (col,lyr) [kg/m2]
real(r8), intent(in) :: h2osno_ice(lbc:ubc,-nlevsno+1:0) ! ice content (col,lyr) [kg/m2]
integer, intent(in) :: snw_rds(lbc:ubc,-nlevsno+1:0) ! snow effective radius (col,lyr) [microns, m^-6]
real(r8), intent(in) :: mss_cnc_aer_in(lbc:ubc,-nlevsno+1:0,sno_nbr_aer) ! mass concentration of all aerosol species (col,lyr,aer) [kg/kg]
real(r8), intent(in) :: albsfc(lbc:ubc,numrad) ! albedo of surface underlying snow (col,bnd) [frc]
real(r8), intent(out) :: albout(lbc:ubc,numrad) ! snow albedo, averaged into 2 bands (=0 if no sun or no snow) (col,bnd) [frc]
real(r8), intent(out) :: flx_abs(lbc:ubc,-nlevsno+1:1,numrad) ! absorbed flux in each layer per unit flux incident on top of snowpack (col,lyr,bnd) [frc]
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
integer, pointer :: snl(:) ! negative number of snow layers (col) [nbr]
real(r8), pointer :: h2osno(:) ! snow liquid water equivalent (col) [kg/m2]
integer, pointer :: clandunit(:) ! corresponding landunit of column (col) [idx] (debugging only)
integer, pointer :: cgridcell(:) ! columns's gridcell index (col) [idx] (debugging only)
integer, pointer :: ltype(:) ! landunit type (lnd) (debugging only)
real(r8), pointer :: londeg(:) ! longitude (degrees) (debugging only)
real(r8), pointer :: latdeg(:) ! latitude (degrees) (debugging only)
!
! !OTHER LOCAL VARIABLES:
!EOP
!-----------------------------------------------------------------------
!
! variables for snow radiative transfer calculations
! Local variables representing single-column values of arrays:
integer :: snl_lcl ! negative number of snow layers [nbr]
integer :: snw_rds_lcl(-nlevsno+1:0) ! snow effective radius [m^-6]
real(r8):: flx_slrd_lcl(1:numrad_snw) ! direct beam incident irradiance [W/m2] (set to 1)
real(r8):: flx_slri_lcl(1:numrad_snw) ! diffuse incident irradiance [W/m2] (set to 1)
real(r8):: mss_cnc_aer_lcl(-nlevsno+1:0,1:sno_nbr_aer) ! aerosol mass concentration (lyr,aer_nbr) [kg/kg]
real(r8):: h2osno_lcl ! total column snow mass [kg/m2]
real(r8):: h2osno_liq_lcl(-nlevsno+1:0) ! liquid water mass [kg/m2]
real(r8):: h2osno_ice_lcl(-nlevsno+1:0) ! ice mass [kg/m2]
real(r8):: albsfc_lcl(1:numrad_snw) ! albedo of underlying surface [frc]
real(r8):: ss_alb_snw_lcl(-nlevsno+1:0) ! single-scatter albedo of ice grains (lyr) [frc]
real(r8):: asm_prm_snw_lcl(-nlevsno+1:0) ! asymmetry parameter of ice grains (lyr) [frc]
real(r8):: ext_cff_mss_snw_lcl(-nlevsno+1:0) ! mass extinction coefficient of ice grains (lyr) [m2/kg]
real(r8):: ss_alb_aer_lcl(sno_nbr_aer) ! single-scatter albedo of aerosol species (aer_nbr) [frc]
real(r8):: asm_prm_aer_lcl(sno_nbr_aer) ! asymmetry parameter of aerosol species (aer_nbr) [frc]
real(r8):: ext_cff_mss_aer_lcl(sno_nbr_aer) ! mass extinction coefficient of aerosol species (aer_nbr) [m2/kg]
! Other local variables
integer :: APRX_TYP ! two-stream approximation type (1=Eddington, 2=Quadrature, 3=Hemispheric Mean) [nbr]
integer :: DELTA ! flag to use Delta approximation (Joseph, 1976) (1= use, 0= don't use)
real(r8):: flx_wgt(1:numrad_snw) ! weights applied to spectral bands, specific to direct and diffuse cases (bnd) [frc]
integer :: flg_nosnl ! flag: =1 if there is snow, but zero snow layers, =0 if at least 1 snow layer [flg]
integer :: trip ! flag: =1 to redo RT calculation if result is unrealistic
integer :: flg_dover ! defines conditions for RT redo (explained below)
real(r8):: albedo ! temporary snow albedo [frc]
real(r8):: flx_sum ! temporary summation variable for NIR weighting
real(r8):: albout_lcl(numrad_snw) ! snow albedo by band [frc]
real(r8):: flx_abs_lcl(-nlevsno+1:1,numrad_snw)! absorbed flux per unit incident flux at top of snowpack (lyr,bnd) [frc]
real(r8):: L_snw(-nlevsno+1:0) ! h2o mass (liquid+solid) in snow layer (lyr) [kg/m2]
real(r8):: tau_snw(-nlevsno+1:0) ! snow optical depth (lyr) [unitless]
real(r8):: L_aer(-nlevsno+1:0,sno_nbr_aer) ! aerosol mass in snow layer (lyr,nbr_aer) [kg/m2]
real(r8):: tau_aer(-nlevsno+1:0,sno_nbr_aer) ! aerosol optical depth (lyr,nbr_aer) [unitless]
real(r8):: tau_sum ! cumulative (snow+aerosol) optical depth [unitless]
real(r8):: tau_clm(-nlevsno+1:0) ! column optical depth from layer bottom to snowpack top (lyr) [unitless]
real(r8):: omega_sum ! temporary summation of single-scatter albedo of all aerosols [frc]
real(r8):: g_sum ! temporary summation of asymmetry parameter of all aerosols [frc]
real(r8):: tau(-nlevsno+1:0) ! weighted optical depth of snow+aerosol layer (lyr) [unitless]
real(r8):: omega(-nlevsno+1:0) ! weighted single-scatter albedo of snow+aerosol layer (lyr) [frc]
real(r8):: g(-nlevsno+1:0) ! weighted asymmetry parameter of snow+aerosol layer (lyr) [frc]
real(r8):: tau_star(-nlevsno+1:0) ! transformed (i.e. Delta-Eddington) optical depth of snow+aerosol layer (lyr) [unitless]
real(r8):: omega_star(-nlevsno+1:0) ! transformed (i.e. Delta-Eddington) SSA of snow+aerosol layer (lyr) [frc]
real(r8):: g_star(-nlevsno+1:0) ! transformed (i.e. Delta-Eddington) asymmetry paramater of snow+aerosol layer (lyr) [frc]
integer :: nstep ! current timestep [nbr] (debugging only)
integer :: g_idx, c_idx, l_idx ! gridcell, column, and landunit indices [idx]
integer :: bnd_idx ! spectral band index (1 <= bnd_idx <= numrad_snw) [idx]
integer :: rds_idx ! snow effective radius index for retrieving Mie parameters from lookup table [idx]
integer :: snl_btm ! index of bottom snow layer (0) [idx]
integer :: snl_top ! index of top snow layer (-4 to 0) [idx]
integer :: fc ! column filter index
integer :: i ! layer index [idx]
integer :: j ! aerosol number index [idx]
integer :: n ! tridiagonal matrix index [idx]
integer :: m ! secondary layer index [idx]
real(r8):: F_direct(-nlevsno+1:0) ! direct-beam radiation at bottom of layer interface (lyr) [W/m^2]
real(r8):: F_net(-nlevsno+1:0) ! net radiative flux at bottom of layer interface (lyr) [W/m^2]
real(r8):: F_abs(-nlevsno+1:0) ! net absorbed radiative energy (lyr) [W/m^2]
real(r8):: F_abs_sum ! total absorbed energy in column [W/m^2]
real(r8):: F_sfc_pls ! upward radiative flux at snowpack top [W/m^2]
real(r8):: F_btm_net ! net flux at bottom of snowpack [W/m^2]
real(r8):: F_sfc_net ! net flux at top of snowpack [W/m^2]
real(r8):: energy_sum ! sum of all energy terms; should be 0.0 [W/m^2]
real(r8):: F_direct_btm ! direct-beam radiation at bottom of snowpack [W/m^2]
real(r8):: mu_not ! cosine of solar zenith angle (used locally) [frc]
integer :: err_idx ! counter for number of times through error loop [nbr]
real(r8):: lat_coord ! gridcell latitude (debugging only)
real(r8):: lon_coord ! gridcell longitude (debugging only)
integer :: sfctype ! underlying surface type (debugging only)
real(r8):: pi ! 3.1415...
! intermediate variables for radiative transfer approximation:
real(r8):: gamma1(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: gamma2(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: gamma3(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: gamma4(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: lambda(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: GAMMA(-nlevsno+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: mu_one ! two-stream coefficient from Toon et al. (lyr) [unitless]
real(r8):: e1(-nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (lyr)
real(r8):: e2(-nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (lyr)
real(r8):: e3(-nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (lyr)
real(r8):: e4(-nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (lyr)
real(r8):: C_pls_btm(-nlevsno+1:0) ! intermediate variable: upward flux at bottom interface (lyr) [W/m2]
real(r8):: C_mns_btm(-nlevsno+1:0) ! intermediate variable: downward flux at bottom interface (lyr) [W/m2]
real(r8):: C_pls_top(-nlevsno+1:0) ! intermediate variable: upward flux at top interface (lyr) [W/m2]
real(r8):: C_mns_top(-nlevsno+1:0) ! intermediate variable: downward flux at top interface (lyr) [W/m2]
real(r8):: A(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: B(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: D(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: E(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: AS(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: DS(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: X(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
real(r8):: Y(-2*nlevsno+1:0) ! tri-diag intermediate variable from Toon et al. (2*lyr)
! Assign local pointers to derived subtypes components (column-level)
! (CLM-specific)
if (flg_snw_ice == 1) then
snl => clm3%g%l%c%cps%snl
h2osno => clm3%g%l%c%cws%h2osno
clandunit => clm3%g%l%c%landunit ! (debug only)
cgridcell => clm3%g%l%c%gridcell ! (debug only)
ltype => clm3%g%l%itype ! (debug only)
londeg => clm3%g%londeg ! (debug only)
latdeg => clm3%g%latdeg ! (debug only)
endif
! Define constants
pi = SHR_CONST_PI
! always use Delta approximation for snow
DELTA = 1
! Get current timestep
nstep = get_nstep()
! Loop over all non-urban columns
! (when called from CSIM, there is only one column)
do fc = 1,num_nourbanc
c_idx = filter_nourbanc(fc)
! Zero absorbed radiative fluxes:
do i=-nlevsno+1,1,1
flx_abs_lcl(:,:) = 0._r8
flx_abs(c_idx,i,:) = 0._r8
enddo
! set snow/ice mass to be used for RT:
if (flg_snw_ice == 1) then
h2osno_lcl = h2osno(c_idx)
else
h2osno_lcl = h2osno_ice(c_idx,0)
endif
! Qualifier for computing snow RT:
! 1) sunlight from atmosphere model
! 2) minimum amount of snow on ground.
! Otherwise, set snow albedo to zero
if ((coszen(c_idx) > 0._r8) .and. (h2osno_lcl > min_snw)) then
! Set variables specific to CLM
if (flg_snw_ice == 1) then
! Assign local (single-column) variables to global values
! If there is snow, but zero snow layers, we must create a layer locally.
! This layer is presumed to have the fresh snow effective radius.
if (snl(c_idx) > -1) then
flg_nosnl = 1
snl_lcl = -1
h2osno_ice_lcl(0) = h2osno_lcl
h2osno_liq_lcl(0) = 0._r8
snw_rds_lcl(0) = nint(snw_rds_min)
else
flg_nosnl = 0
snl_lcl = snl(c_idx)
h2osno_liq_lcl(:) = h2osno_liq(c_idx,:)
h2osno_ice_lcl(:) = h2osno_ice(c_idx,:)
snw_rds_lcl(:) = snw_rds(c_idx,:)
endif
snl_btm = 0
snl_top = snl_lcl+1
! for debugging only
l_idx = clandunit(c_idx)
g_idx = cgridcell(c_idx)
sfctype = ltype(l_idx)
lat_coord = latdeg(g_idx)
lon_coord = londeg(g_idx)
! Set variables specific to CSIM
else
flg_nosnl = 0
snl_lcl = -1
h2osno_liq_lcl(:) = h2osno_liq(c_idx,:)
h2osno_ice_lcl(:) = h2osno_ice(c_idx,:)
snw_rds_lcl(:) = snw_rds(c_idx,:)
snl_btm = 0
snl_top = 0
sfctype = -1
lat_coord = -90
lon_coord = 0
endif
! Set local aerosol array
do j=1,sno_nbr_aer
mss_cnc_aer_lcl(:,j) = mss_cnc_aer_in(c_idx,:,j)
enddo
! Set spectral underlying surface albedos to their corresponding VIS or NIR albedos
albsfc_lcl(1) = albsfc(c_idx,1)
albsfc_lcl(nir_bnd_bgn:nir_bnd_end) = albsfc(c_idx,2)
! Error check for snow grain size:
do i=snl_top,snl_btm,1
if ((snw_rds_lcl(i) < snw_rds_min_tbl) .or. (snw_rds_lcl(i) > snw_rds_max_tbl)) then
write (iulog,*) "SNICAR ERROR: snow grain radius of ", snw_rds_lcl(i), " out of bounds."
write (iulog,*) "NSTEP= ", nstep
write (iulog,*) "flg_snw_ice= ", flg_snw_ice
write (iulog,*) "column: ", c_idx, " level: ", i, " snl(c)= ", snl_lcl
write (iulog,*) "lat= ", lat_coord, " lon= ", lon_coord
write (iulog,*) "h2osno(c)= ", h2osno_lcl
call endrun()
endif
enddo
! Incident flux weighting parameters
! - sum of all VIS bands must equal 1
! - sum of all NIR bands must equal 1
!
! Spectral bands (5-band case)
! Band 1: 0.3-0.7um (VIS)
! Band 2: 0.7-1.0um (NIR)
! Band 3: 1.0-1.2um (NIR)
! Band 4: 1.2-1.5um (NIR)
! Band 5: 1.5-5.0um (NIR)
!
! The following weights are appropriate for surface-incident flux in a mid-latitude winter atmosphere
!
! 3-band weights
if (numrad_snw==3) then
! Direct:
if (flg_slr_in == 1) then
flx_wgt(1) = 1._r8
flx_wgt(2) = 0.66628670195247_r8
flx_wgt(3) = 0.33371329804753_r8
! Diffuse:
elseif (flg_slr_in == 2) then
flx_wgt(1) = 1._r8
flx_wgt(2) = 0.77887652162877_r8
flx_wgt(3) = 0.22112347837123_r8
endif
! 5-band weights
elseif(numrad_snw==5) then
! Direct:
if (flg_slr_in == 1) then
flx_wgt(1) = 1._r8
flx_wgt(2) = 0.49352158521175_r8
flx_wgt(3) = 0.18099494230665_r8
flx_wgt(4) = 0.12094898498813_r8
flx_wgt(5) = 0.20453448749347_r8
! Diffuse:
elseif (flg_slr_in == 2) then
flx_wgt(1) = 1._r8
flx_wgt(2) = 0.58581507618433_r8
flx_wgt(3) = 0.20156903770812_r8
flx_wgt(4) = 0.10917889346386_r8
flx_wgt(5) = 0.10343699264369_r8
endif
endif
! Loop over snow spectral bands
do bnd_idx = 1,numrad_snw
mu_not = coszen(c_idx) ! must set here, because of error handling
flg_dover = 1 ! default is to redo
err_idx = 0 ! number of times through loop
do while (flg_dover > 0)
! DEFAULT APPROXIMATIONS:
! VIS: Delta-Eddington
! NIR (all): Delta-Hemispheric Mean
! WARNING: DO NOT USE DELTA-EDDINGTON FOR NIR DIFFUSE - this sometimes results in negative albedo
!
! ERROR CONDITIONS:
! Conditions which cause "trip", resulting in redo of RT approximation:
! 1. negative absorbed flux
! 2. total absorbed flux greater than incident flux
! 3. negative albedo
! NOTE: These errors have only been encountered in spectral bands 4 and 5
!
! ERROR HANDLING
! 1st error (flg_dover=2): switch approximation (Edd->HM or HM->Edd)
! 2nd error (flg_dover=3): change zenith angle by 0.02 (this happens about 1 in 10^6 cases)
! 3rd error (flg_dover=4): switch approximation with new zenith
! Subsequent errors: repeatedly change zenith and approximations...
if (bnd_idx == 1) then
if (flg_dover == 2) then
APRX_TYP = 3
elseif (flg_dover == 3) then
APRX_TYP = 1
if (coszen(c_idx) > 0.5_r8) then
mu_not = mu_not - 0.02_r8
else
mu_not = mu_not + 0.02_r8
endif
elseif (flg_dover == 4) then
APRX_TYP = 3
else
APRX_TYP = 1
endif
else
if (flg_dover == 2) then
APRX_TYP = 1
elseif (flg_dover == 3) then
APRX_TYP = 3
if (coszen(c_idx) > 0.5_r8) then
mu_not = mu_not - 0.02_r8
else
mu_not = mu_not + 0.02_r8
endif
elseif (flg_dover == 4) then
APRX_TYP = 1
else
APRX_TYP = 3
endif
endif
! Set direct or diffuse incident irradiance to 1
! (This has to be within the bnd loop because mu_not is adjusted in rare cases)
if (flg_slr_in == 1) then
flx_slrd_lcl(bnd_idx) = 1._r8/(mu_not*pi) ! this corresponds to incident irradiance of 1.0
flx_slri_lcl(bnd_idx) = 0._r8
else
flx_slrd_lcl(bnd_idx) = 0._r8
flx_slri_lcl(bnd_idx) = 1._r8
endif
! Pre-emptive error handling: aerosols can reap havoc on these absorptive bands.
! Since extremely high soot concentrations have a negligible effect on these bands, zero them.
if ( (numrad_snw == 5).and.((bnd_idx == 5).or.(bnd_idx == 4)) ) then
mss_cnc_aer_lcl(:,:) = 0._r8
endif
if ( (numrad_snw == 3).and.(bnd_idx == 3) ) then
mss_cnc_aer_lcl(:,:) = 0._r8
endif
! Define local Mie parameters based on snow grain size and aerosol species,
! retrieved from a lookup table.
if (flg_slr_in == 1) then
do i=snl_top,snl_btm,1
rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
! snow optical properties (direct radiation)
ss_alb_snw_lcl(i) = ss_alb_snw_drc(rds_idx,bnd_idx)
asm_prm_snw_lcl(i) = asm_prm_snw_drc(rds_idx,bnd_idx)
ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_drc(rds_idx,bnd_idx)
enddo
elseif (flg_slr_in == 2) then
do i=snl_top,snl_btm,1
rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
! snow optical properties (diffuse radiation)
ss_alb_snw_lcl(i) = ss_alb_snw_dfs(rds_idx,bnd_idx)
asm_prm_snw_lcl(i) = asm_prm_snw_dfs(rds_idx,bnd_idx)
ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_dfs(rds_idx,bnd_idx)
enddo
endif
! aerosol species 1 optical properties
ss_alb_aer_lcl(1) = ss_alb_bc1(1,bnd_idx)
asm_prm_aer_lcl(1) = asm_prm_bc1(1,bnd_idx)
ext_cff_mss_aer_lcl(1) = ext_cff_mss_bc1(1,bnd_idx)
! aerosol species 2 optical properties
ss_alb_aer_lcl(2) = ss_alb_bc2(1,bnd_idx)
asm_prm_aer_lcl(2) = asm_prm_bc2(1,bnd_idx)
ext_cff_mss_aer_lcl(2) = ext_cff_mss_bc2(1,bnd_idx)
! aerosol species 3 optical properties
ss_alb_aer_lcl(3) = ss_alb_oc1(1,bnd_idx)
asm_prm_aer_lcl(3) = asm_prm_oc1(1,bnd_idx)
ext_cff_mss_aer_lcl(3) = ext_cff_mss_oc1(1,bnd_idx)
! aerosol species 4 optical properties
ss_alb_aer_lcl(4) = ss_alb_oc2(1,bnd_idx)
asm_prm_aer_lcl(4) = asm_prm_oc2(1,bnd_idx)
ext_cff_mss_aer_lcl(4) = ext_cff_mss_oc2(1,bnd_idx)
! aerosol species 5 optical properties
ss_alb_aer_lcl(5) = ss_alb_dst1(1,bnd_idx)
asm_prm_aer_lcl(5) = asm_prm_dst1(1,bnd_idx)
ext_cff_mss_aer_lcl(5) = ext_cff_mss_dst1(1,bnd_idx)
! aerosol species 6 optical properties
ss_alb_aer_lcl(6) = ss_alb_dst2(1,bnd_idx)
asm_prm_aer_lcl(6) = asm_prm_dst2(1,bnd_idx)
ext_cff_mss_aer_lcl(6) = ext_cff_mss_dst2(1,bnd_idx)
! aerosol species 7 optical properties
ss_alb_aer_lcl(7) = ss_alb_dst3(1,bnd_idx)
asm_prm_aer_lcl(7) = asm_prm_dst3(1,bnd_idx)
ext_cff_mss_aer_lcl(7) = ext_cff_mss_dst3(1,bnd_idx)
! aerosol species 8 optical properties
ss_alb_aer_lcl(8) = ss_alb_dst4(1,bnd_idx)
asm_prm_aer_lcl(8) = asm_prm_dst4(1,bnd_idx)
ext_cff_mss_aer_lcl(8) = ext_cff_mss_dst4(1,bnd_idx)
! 1. snow and aerosol layer column mass (L_snw, L_aer [kg/m^2])
! 2. optical Depths (tau_snw, tau_aer)
! 3. weighted Mie properties (tau, omega, g)
! Weighted Mie parameters of each layer
do i=snl_top,snl_btm,1
L_snw(i) = h2osno_ice_lcl(i)+h2osno_liq_lcl(i)
tau_snw(i) = L_snw(i)*ext_cff_mss_snw_lcl(i)
do j=1,sno_nbr_aer
L_aer(i,j) = L_snw(i)*mss_cnc_aer_lcl(i,j)
tau_aer(i,j) = L_aer(i,j)*ext_cff_mss_aer_lcl(j)
enddo
tau_sum = 0._r8
omega_sum = 0._r8
g_sum = 0._r8
do j=1,sno_nbr_aer
tau_sum = tau_sum + tau_aer(i,j)
omega_sum = omega_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j))
g_sum = g_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j)*asm_prm_aer_lcl(j))
enddo
tau(i) = tau_sum + tau_snw(i)
omega(i) = (1/tau(i))*(omega_sum+(ss_alb_snw_lcl(i)*tau_snw(i)))
g(i) = (1/(tau(i)*omega(i)))*(g_sum+ (asm_prm_snw_lcl(i)*ss_alb_snw_lcl(i)*tau_snw(i)))
enddo
! DELTA transformations, if requested
if (DELTA == 1) then
do i=snl_top,snl_btm,1
g_star(i) = g(i)/(1+g(i))
omega_star(i) = ((1-(g(i)**2))*omega(i)) / (1-(omega(i)*(g(i)**2)))
tau_star(i) = (1-(omega(i)*(g(i)**2)))*tau(i)
enddo
else
do i=snl_top,snl_btm,1
g_star(i) = g(i)
omega_star(i) = omega(i)
tau_star(i) = tau(i)
enddo
endif
! Total column optical depth:
! tau_clm(i) = total optical depth above the bottom of layer i
tau_clm(snl_top) = 0._r8
do i=snl_top+1,snl_btm,1
tau_clm(i) = tau_clm(i-1)+tau_star(i-1)
enddo
! Direct radiation at bottom of snowpack:
F_direct_btm = albsfc_lcl(bnd_idx)*mu_not*exp(-(tau_clm(snl_btm)+tau_star(snl_btm))/mu_not)*pi*flx_slrd_lcl(bnd_idx)
! Intermediates
! Gamma values are approximation-specific.
! Eddington
if (APRX_TYP==1) then
do i=snl_top,snl_btm,1
gamma1(i) = (7-(omega_star(i)*(4+(3*g_star(i)))))/4
gamma2(i) = -(1-(omega_star(i)*(4-(3*g_star(i)))))/4
gamma3(i) = (2-(3*g_star(i)*mu_not))/4
gamma4(i) = 1-gamma3(i)
mu_one = 0.5
enddo
! Quadrature
elseif (APRX_TYP==2) then
do i=snl_top,snl_btm,1
gamma1(i) = (3**0.5)*(2-(omega_star(i)*(1+g_star(i))))/2
gamma2(i) = omega_star(i)*(3**0.5)*(1-g_star(i))/2
gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
gamma4(i) = 1-gamma3(i)
mu_one = 1/(3**0.5)
enddo
! Hemispheric Mean
elseif (APRX_TYP==3) then
do i=snl_top,snl_btm,1
gamma1(i) = 2 - (omega_star(i)*(1+g_star(i)))
gamma2(i) = omega_star(i)*(1-g_star(i))
gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
gamma4(i) = 1-gamma3(i)
mu_one = 0.5
enddo
endif
! Intermediates for tri-diagonal solution
do i=snl_top,snl_btm,1
lambda(i) = sqrt(abs((gamma1(i)**2) - (gamma2(i)**2)))
GAMMA(i) = gamma2(i)/(gamma1(i)+lambda(i))
e1(i) = 1+(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
e2(i) = 1-(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
e3(i) = GAMMA(i) + exp(-lambda(i)*tau_star(i))
e4(i) = GAMMA(i) - exp(-lambda(i)*tau_star(i))
enddo !enddo over snow layers
! Intermediates for tri-diagonal solution
do i=snl_top,snl_btm,1
if (flg_slr_in == 1) then
C_pls_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-(tau_clm(i)+tau_star(i))/mu_not)* &
(((gamma1(i)-(1/mu_not))*gamma3(i))+ &
(gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_mns_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-(tau_clm(i)+tau_star(i))/mu_not)* &
(((gamma1(i)+(1/mu_not))*gamma4(i))+ &
(gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_pls_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-tau_clm(i)/mu_not)*(((gamma1(i)-(1/mu_not))* &
gamma3(i))+(gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_mns_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-tau_clm(i)/mu_not)*(((gamma1(i)+(1/mu_not))* &
gamma4(i))+(gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))
else
C_pls_btm(i) = 0._r8
C_mns_btm(i) = 0._r8
C_pls_top(i) = 0._r8
C_mns_top(i) = 0._r8
endif
enddo
! Coefficients for tridiaganol matrix solution
do i=2*snl_lcl+1,0,1
!Boundary values for i=1 and i=2*snl_lcl, specifics for i=odd and i=even
if (i==(2*snl_lcl+1)) then
A(i) = 0
B(i) = e1(snl_top)
D(i) = -e2(snl_top)
E(i) = flx_slri_lcl(bnd_idx)-C_mns_top(snl_top)
elseif(i==0) then
A(i) = e1(snl_btm)-(albsfc_lcl(bnd_idx)*e3(snl_btm))
B(i) = e2(snl_btm)-(albsfc_lcl(bnd_idx)*e4(snl_btm))
D(i) = 0
E(i) = F_direct_btm-C_pls_btm(snl_btm)+(albsfc_lcl(bnd_idx)*C_mns_btm(snl_btm))
elseif(mod(i,2)==-1) then ! If odd and i>=3 (n=1 for i=3)
n=floor(i/2.0)
A(i) = (e2(n)*e3(n))-(e4(n)*e1(n))
B(i) = (e1(n)*e1(n+1))-(e3(n)*e3(n+1))
D(i) = (e3(n)*e4(n+1))-(e1(n)*e2(n+1))
E(i) = (e3(n)*(C_pls_top(n+1)-C_pls_btm(n)))+(e1(n)*(C_mns_btm(n)-C_mns_top(n+1)))
elseif(mod(i,2)==0) then ! If even and i<=2*snl_lcl
n=(i/2)
A(i) = (e2(n+1)*e1(n))-(e3(n)*e4(n+1))
B(i) = (e2(n)*e2(n+1))-(e4(n)*e4(n+1))
D(i) = (e1(n+1)*e4(n+1))-(e2(n+1)*e3(n+1))
E(i) = (e2(n+1)*(C_pls_top(n+1)-C_pls_btm(n)))+(e4(n+1)*(C_mns_top(n+1)-C_mns_btm(n)))
endif
enddo
AS(0) = A(0)/B(0)
DS(0) = E(0)/B(0)
do i=-1,(2*snl_lcl+1),-1
X(i) = 1/(B(i)-(D(i)*AS(i+1)))
AS(i) = A(i)*X(i)
DS(i) = (E(i)-(D(i)*DS(i+1)))*X(i)
enddo
Y(2*snl_lcl+1) = DS(2*snl_lcl+1)
do i=(2*snl_lcl+2),0,1
Y(i) = DS(i)-(AS(i)*Y(i-1))
enddo
! Downward direct-beam and net flux (F_net) at the base of each layer:
do i=snl_top,snl_btm,1
F_direct(i) = mu_not*pi*flx_slrd_lcl(bnd_idx)*exp(-(tau_clm(i)+tau_star(i))/mu_not)
F_net(i) = (Y(2*i-1)*(e1(i)-e3(i))) + (Y(2*i)*(e2(i)-e4(i))) + &
C_pls_btm(i) - C_mns_btm(i) - F_direct(i)
enddo
! Upward flux at snowpack top:
F_sfc_pls = (Y(2*snl_lcl+1)*(exp(-lambda(snl_top)*tau_star(snl_top))+ &
GAMMA(snl_top))) + (Y(2*snl_lcl+2)*(exp(-lambda(snl_top)* &
tau_star(snl_top))-GAMMA(snl_top))) + C_pls_top(snl_top)
! Net flux at bottom = absorbed radiation by underlying surface:
F_btm_net = -F_net(snl_btm)
! Bulk column albedo and surface net flux
albedo = F_sfc_pls/((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))
F_sfc_net = F_sfc_pls - ((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))
trip = 0
! Absorbed flux in each layer
do i=snl_top,snl_btm,1
if(i==snl_top) then
F_abs(i) = F_net(i)-F_sfc_net
else
F_abs(i) = F_net(i)-F_net(i-1)
endif
flx_abs_lcl(i,bnd_idx) = F_abs(i)
! ERROR check: negative absorption
if (flx_abs_lcl(i,bnd_idx) < -0.00001) then
trip = 1
endif
enddo
flx_abs_lcl(1,bnd_idx) = F_btm_net
if (flg_nosnl == 1) then
! If there are no snow layers (but still snow), all absorbed energy must be in top soil layer
!flx_abs_lcl(:,bnd_idx) = 0._r8
!flx_abs_lcl(1,bnd_idx) = F_abs(0) + F_btm_net
! changed on 20070408:
! OK to put absorbed energy in the fictitous snow layer because routine SurfaceRadiation
! handles the case of no snow layers. Then, if a snow layer is addded between now and
! SurfaceRadiation (called in Hydrology1), absorbed energy will be properly distributed.
flx_abs_lcl(0,bnd_idx) = F_abs(0)
flx_abs_lcl(1,bnd_idx) = F_btm_net
endif
!Underflow check (we've already tripped the error condition above)
do i=snl_top,1,1
if (flx_abs_lcl(i,bnd_idx) < 0._r8) then
flx_abs_lcl(i,bnd_idx) = 0._r8
endif
enddo
F_abs_sum = 0._r8
do i=snl_top,snl_btm,1
F_abs_sum = F_abs_sum + F_abs(i)
enddo
!ERROR check: absorption greater than incident flux
! (should make condition more generic than "1._r8")
if (F_abs_sum > 1._r8) then
trip = 1
endif
!ERROR check:
if ((albedo < 0._r8).and.(trip==0)) then
trip = 1
endif
! Set conditions for redoing RT calculation
if ((trip == 1).and.(flg_dover == 1)) then
flg_dover = 2
elseif ((trip == 1).and.(flg_dover == 2)) then
flg_dover = 3
elseif ((trip == 1).and.(flg_dover == 3)) then
flg_dover = 4
elseif((trip == 1).and.(flg_dover == 4).and.(err_idx < 20)) then
flg_dover = 3
err_idx = err_idx + 1
elseif((trip == 1).and.(flg_dover == 4).and.(err_idx >= 20)) then
flg_dover = 0
write(iulog,*) "SNICAR ERROR: FOUND A WORMHOLE. STUCK IN INFINITE LOOP! Called from: ", flg_snw_ice
write(iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(c_idx,0)
write(iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
write(iulog,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
write(iulog,*) "SNICAR STATS: soot1(0)= ", mss_cnc_aer_lcl(0,1)
write(iulog,*) "SNICAR STATS: soot2(0)= ", mss_cnc_aer_lcl(0,2)
write(iulog,*) "SNICAR STATS: dust1(0)= ", mss_cnc_aer_lcl(0,3)
write(iulog,*) "SNICAR STATS: dust2(0)= ", mss_cnc_aer_lcl(0,4)
write(iulog,*) "SNICAR STATS: dust3(0)= ", mss_cnc_aer_lcl(0,5)
write(iulog,*) "SNICAR STATS: dust4(0)= ", mss_cnc_aer_lcl(0,6)
call endrun()
else
flg_dover = 0
endif
enddo !enddo while (flg_dover > 0)
! Energy conservation check:
! Incident direct+diffuse radiation equals (absorbed+bulk_transmitted+bulk_reflected)
energy_sum = (mu_not*pi*flx_slrd_lcl(bnd_idx)) + flx_slri_lcl(bnd_idx) - (F_abs_sum + F_btm_net + F_sfc_pls)
if (abs(energy_sum) > 0.00001_r8) then
write (iulog,"(a,e12.6,a,i6,a,i6)") "SNICAR ERROR: Energy conservation error of : ", energy_sum, &
" at timestep: ", nstep, " at column: ", c_idx
call endrun()
endif
albout_lcl(bnd_idx) = albedo
! Check that albedo is less than 1
if (albout_lcl(bnd_idx) > 1.0) then
write (iulog,*) "SNICAR ERROR: Albedo > 1.0 at c: ", c_idx, " NSTEP= ",nstep
write (iulog,*) "SNICAR STATS: bnd_idx= ",bnd_idx
write (iulog,*) "SNICAR STATS: albout_lcl(bnd)= ",albout_lcl(bnd_idx), " albsfc_lcl(bnd_idx)= ",albsfc_lcl(bnd_idx)
write (iulog,*) "SNICAR STATS: landtype= ", sfctype
write (iulog,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
write (iulog,*) "SNICAR STATS: coszen= ", coszen(c_idx), " flg_slr= ", flg_slr_in
write (iulog,*) "SNICAR STATS: soot(-4)= ", mss_cnc_aer_lcl(-4,1)
write (iulog,*) "SNICAR STATS: soot(-3)= ", mss_cnc_aer_lcl(-3,1)
write (iulog,*) "SNICAR STATS: soot(-2)= ", mss_cnc_aer_lcl(-2,1)
write (iulog,*) "SNICAR STATS: soot(-1)= ", mss_cnc_aer_lcl(-1,1)
write (iulog,*) "SNICAR STATS: soot(0)= ", mss_cnc_aer_lcl(0,1)
write (iulog,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
write (iulog,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
write (iulog,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
write (iulog,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
write (iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
write (iulog,*) "SNICAR STATS: snw_rds(-4)= ", snw_rds(c_idx,-4)
write (iulog,*) "SNICAR STATS: snw_rds(-3)= ", snw_rds(c_idx,-3)
write (iulog,*) "SNICAR STATS: snw_rds(-2)= ", snw_rds(c_idx,-2)
write (iulog,*) "SNICAR STATS: snw_rds(-1)= ", snw_rds(c_idx,-1)
write (iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(c_idx,0)
call endrun()
endif
enddo ! loop over wvl bands
! Weight output NIR albedo appropriately
albout(c_idx,1) = albout_lcl(1)
flx_sum = 0._r8
do bnd_idx= nir_bnd_bgn,nir_bnd_end
flx_sum = flx_sum + flx_wgt(bnd_idx)*albout_lcl(bnd_idx)
end do
albout(c_idx,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
! Weight output NIR absorbed layer fluxes (flx_abs) appropriately
flx_abs(c_idx,:,1) = flx_abs_lcl(:,1)
do i=snl_top,1,1
flx_sum = 0._r8
do bnd_idx= nir_bnd_bgn,nir_bnd_end
flx_sum = flx_sum + flx_wgt(bnd_idx)*flx_abs_lcl(i,bnd_idx)
enddo
flx_abs(c_idx,i,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
end do
! Write diagnostics, if desired. (default is to not compile this)
#if 0
write(iulog,*) "SNICAR STATS: NSTEP= ", nstep
write(iulog,*) "SNICAR STATS: Col: ", c_idx
write(iulog,*) "SNICAR STATS: snl(c)= ",snl_lcl
write(iulog,*) "SNICAR STATS: cosine zenith= ", coszen(c_idx)
write(iulog,*) "SNICAR STATS: h2osno(c): ", h2osno_lcl
write(iulog,*) "SNICAR STATS: albout_lcl(1): ", albout_lcl(1)
write(iulog,*) "SNICAR STATS: albout_lcl(2): ", albout_lcl(2)
write(iulog,*) "SNICAR STATS: albout_lcl(3): ", albout_lcl(3)
write(iulog,*) "SNICAR STATS: albout_lcl(4): ", albout_lcl(4)
write(iulog,*) "SNICAR STATS: albout_lcl(5): ", albout_lcl(5)
write(iulog,*) "SNICAR STATS: albout(1): ", albout(c_idx,1)
write(iulog,*) "SNICAR STATS: albout(2): ", albout(c_idx,2)
write(iulog,*) "SNICAR STATS: NIR flx_abs(-4): ", flx_abs(c_idx,-4,2)
write(iulog,*) "SNICAR STATS: NIR flx_abs(-3): ", flx_abs(c_idx,-3,2)
write(iulog,*) "SNICAR STATS: NIR flx_abs(-2): ", flx_abs(c_idx,-2,2)
write(iulog,*) "SNICAR STATS: NIR flx_abs(-1): ", flx_abs(c_idx,-1,2)
write(iulog,*) "SNICAR STATS: NIR flx_abs(0): ", flx_abs(c_idx,0,2)
write(iulog,*) "SNICAR STATS: TOPLYR ABS, BND 1= ", flx_abs_lcl(snl_top,1)
write(iulog,*) "SNICAR STATS: TOPLYR ABS, BND 2= ", flx_abs_lcl(snl_top,2)
write(iulog,*) "SNICAR STATS: TOPLYR ABS, BND 3= ", flx_abs_lcl(snl_top,3)
write(iulog,*) "SNICAR STATS: TOPLYR ABS, BND 4= ", flx_abs_lcl(snl_top,4)
write(iulog,*) "SNICAR STATS: TOPLYR ABS, BND 5= ", flx_abs_lcl(snl_top,5)
write (iulog,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
write (iulog,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
write (iulog,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
write (iulog,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
write (iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
write (iulog,*) "SNICAR STATS: snw_rds(-4)= ", snw_rds(c_idx,-4)
write (iulog,*) "SNICAR STATS: snw_rds(-3)= ", snw_rds(c_idx,-3)
write (iulog,*) "SNICAR STATS: snw_rds(-2)= ", snw_rds(c_idx,-2)
write (iulog,*) "SNICAR STATS: snw_rds(-1)= ", snw_rds(c_idx,-1)
write (iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(c_idx,0)
#endif
! If snow < minimum_snow, but > 0, and there is sun, set albedo to underlying surface albedo
elseif ( (coszen(c_idx) > 0._r8) .and. (h2osno_lcl < min_snw) .and. (h2osno_lcl > 0._r8) ) then
albout(c_idx,1) = albsfc(c_idx,1)
albout(c_idx,2) = albsfc(c_idx,2)
! There is either zero snow, or no sun
else
albout(c_idx,1) = 0._r8
albout(c_idx,2) = 0._r8
endif ! if column has snow and coszen > 0
enddo ! loop over all columns
end subroutine SNICAR_RT
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: SnowAge_grain
!
! !INTERFACE:
subroutine SnowAge_grain(lbc, ubc, num_snowc, filter_snowc, num_nosnowc, filter_nosnowc) 1,7
!
! !DESCRIPTION:
! Updates the snow effective grain size (radius).
! Contributions to grain size evolution are from:
! 1. vapor redistribution (dry snow)
! 2. liquid water redistribution (wet snow)
! 3. re-freezing of liquid water
!
! Vapor redistribution: Method is to retrieve 3 best-bit parameters that
! depend on snow temperature, temperature gradient, and density,
! that are derived from the microphysical model described in:
! Flanner and Zender (2006), Linking snowpack microphysics and albedo
! evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834.
! The parametric equation has the form:
! dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), where:
! r is the effective radius,
! tau and kappa are best-fit parameters,
! drdt_0 is the initial rate of change of effective radius, and
! dr_fresh is the difference between the current and fresh snow states
! (r_current - r_fresh).
!
! Liquid water redistribution: Apply the grain growth function from:
! Brun, E. (1989), Investigation of wet-snow metamorphism in respect of
! liquid-water content, Annals of Glaciology, 13, 22-26.
! There are two parameters that describe the grain growth rate as
! a function of snow liquid water content (LWC). The "LWC=0" parameter
! is zeroed here because we are accounting for dry snowing with a
! different representation
!
! Re-freezing of liquid water: Assume that re-frozen liquid water clumps
! into an arbitrarily large effective grain size (snw_rds_refrz).
! The phenomenon is observed (Grenfell), but so far unquantified, as far as
! I am aware.
!
!
! !USES:
use clmtype
use clm_time_manager , only : get_step_size, get_nstep
use clm_varpar , only : nlevsno
use clm_varcon , only : spval
use abortutils , only : endrun
use shr_const_mod , only : SHR_CONST_RHOICE, SHR_CONST_PI
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbc, ubc ! column bounds
integer, intent(in) :: num_snowc ! number of column snow points in column filter
integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
integer, intent(in) :: num_nosnowc ! number of column non-snow points in column filter
integer, intent(in) :: filter_nosnowc(ubc-lbc+1) ! column filter for non-snow points
!
!
! !CALLED FROM: clm_driver1
!
! !LOCAL VARIABLES:
!
! local pointers to implicit arguments
!
real(r8), pointer :: t_soisno(:,:) ! soil and snow temperature (col,lyr) [K]
integer, pointer :: snl(:) ! negative number of snow layers (col) [nbr]
real(r8), pointer :: t_grnd(:) ! ground temperature (col) [K]
real(r8), pointer :: dz(:,:) ! layer thickness (col,lyr) [m]
real(r8), pointer :: h2osno(:) ! snow water (col) [mm H2O]
real(r8), pointer :: snw_rds(:,:) ! effective grain radius (col,lyr) [microns, m-6]
real(r8), pointer :: snw_rds_top(:) ! effective grain radius, top layer (col) [microns, m-6]
real(r8), pointer :: sno_liq_top(:) ! liquid water fraction (mass) in top snow layer (col) [frc]
real(r8), pointer :: h2osoi_liq(:,:) ! liquid water content (col,lyr) [kg m-2]
real(r8), pointer :: h2osoi_ice(:,:) ! ice content (col,lyr) [kg m-2]
real(r8), pointer :: snot_top(:) ! snow temperature in top layer (col) [K]
real(r8), pointer :: dTdz_top(:) ! temperature gradient in top layer (col) [K m-1]
real(r8), pointer :: qflx_snow_grnd_col(:) ! snow on ground after interception (col) [kg m-2 s-1]
real(r8), pointer :: qflx_snwcp_ice(:) ! excess precipitation due to snow capping [kg m-2 s-1]
real(r8), pointer :: qflx_snofrz_lyr(:,:) ! snow freezing rate (col,lyr) [kg m-2 s-1]
logical , pointer :: do_capsnow(:) ! true => do snow capping
!
! !OTHER LOCAL VARIABLES:
!
integer :: snl_top ! top snow layer index [idx]
integer :: snl_btm ! bottom snow layer index [idx]
integer :: i ! layer index [idx]
integer :: c_idx ! column index [idx]
integer :: fc ! snow column filter index [idx]
integer :: T_idx ! snow aging lookup table temperature index [idx]
integer :: Tgrd_idx ! snow aging lookup table temperature gradient index [idx]
integer :: rhos_idx ! snow aging lookup table snow density index [idx]
real(r8) :: t_snotop ! temperature at upper layer boundary [K]
real(r8) :: t_snobtm ! temperature at lower layer boundary [K]
real(r8) :: dTdz(lbc:ubc,-nlevsno:0) ! snow temperature gradient (col,lyr) [K m-1]
real(r8) :: bst_tau ! snow aging parameter retrieved from lookup table [hour]
real(r8) :: bst_kappa ! snow aging parameter retrieved from lookup table [unitless]
real(r8) :: bst_drdt0 ! snow aging parameter retrieved from lookup table [um hr-1]
real(r8) :: dr ! incremental change in snow effective radius [um]
real(r8) :: dr_wet ! incremental change in snow effective radius from wet growth [um]
real(r8) :: dr_fresh ! difference between fresh snow r_e and current r_e [um]
real(r8) :: newsnow ! fresh snowfall [kg m-2]
real(r8) :: refrzsnow ! re-frozen snow [kg m-2]
real(r8) :: frc_newsnow ! fraction of layer mass that is new snow [frc]
real(r8) :: frc_oldsnow ! fraction of layer mass that is old snow [frc]
real(r8) :: frc_refrz ! fraction of layer mass that is re-frozen snow [frc]
real(r8) :: frc_liq ! fraction of layer mass that is liquid water[frc]
real(r8) :: dtime ! land model time step [sec]
real(r8) :: rhos ! snow density [kg m-3]
real(r8) :: h2osno_lyr ! liquid + solid H2O in snow layer [kg m-2]
! Assign local pointers to derived subtypes components (column-level)
t_soisno => clm3%g%l%c%ces%t_soisno
snl => clm3%g%l%c%cps%snl
t_grnd => clm3%g%l%c%ces%t_grnd
dz => clm3%g%l%c%cps%dz
h2osno => clm3%g%l%c%cws%h2osno
snw_rds => clm3%g%l%c%cps%snw_rds
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
snot_top => clm3%g%l%c%cps%snot_top
dTdz_top => clm3%g%l%c%cps%dTdz_top
snw_rds_top => clm3%g%l%c%cps%snw_rds_top
sno_liq_top => clm3%g%l%c%cps%sno_liq_top
qflx_snow_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_snow_grnd
qflx_snwcp_ice => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_ice
qflx_snofrz_lyr => clm3%g%l%c%cwf%qflx_snofrz_lyr
do_capsnow => clm3%g%l%c%cps%do_capsnow
! set timestep and step interval
dtime = get_step_size()
! loop over columns that have at least one snow layer
do fc = 1, num_snowc
c_idx = filter_snowc(fc)
snl_btm = 0
snl_top = snl(c_idx) + 1
! loop over snow layers
do i=snl_top,snl_btm,1
!
!********** 1. DRY SNOW AGING ***********
!
h2osno_lyr = h2osoi_liq(c_idx,i) + h2osoi_ice(c_idx,i)
! temperature gradient
if (i == snl_top) then
! top layer
t_snotop = t_grnd(c_idx)
t_snobtm = (t_soisno(c_idx,i+1)*dz(c_idx,i) + t_soisno(c_idx,i)*dz(c_idx,i+1)) / (dz(c_idx,i)+dz(c_idx,i+1))
else
t_snotop = (t_soisno(c_idx,i-1)*dz(c_idx,i) + t_soisno(c_idx,i)*dz(c_idx,i-1)) / (dz(c_idx,i)+dz(c_idx,i-1))
t_snobtm = (t_soisno(c_idx,i+1)*dz(c_idx,i) + t_soisno(c_idx,i)*dz(c_idx,i+1)) / (dz(c_idx,i)+dz(c_idx,i+1))
endif
dTdz(c_idx,i) = abs((t_snotop - t_snobtm) / dz(c_idx,i))
! snow density
rhos = (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i)) / dz(c_idx,i)
! best-fit table indecies
T_idx = nint((t_soisno(c_idx,i)-223) / 5) + 1
Tgrd_idx = nint(dTdz(c_idx,i) / 10) + 1
rhos_idx = nint((rhos-50) / 50) + 1
! boundary check:
if (T_idx < idx_T_min) then
T_idx = idx_T_min
endif
if (T_idx > idx_T_max) then
T_idx = idx_T_max
endif
if (Tgrd_idx < idx_Tgrd_min) then
Tgrd_idx = idx_Tgrd_min
endif
if (Tgrd_idx > idx_Tgrd_max) then
Tgrd_idx = idx_Tgrd_max
endif
if (rhos_idx < idx_rhos_min) then
rhos_idx = idx_rhos_min
endif
if (rhos_idx > idx_rhos_max) then
rhos_idx = idx_rhos_max
endif
! best-fit parameters
bst_tau = snowage_tau(rhos_idx,Tgrd_idx,T_idx)
bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx)
bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx)
! change in snow effective radius, using best-fit parameters
dr_fresh = snw_rds(c_idx,i)-snw_rds_min
dr = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1/bst_kappa)) * (dtime/3600)
!
!********** 2. WET SNOW AGING ***********
!
! We are assuming wet and dry evolution occur simultaneously, and
! the contributions from both can be summed.
! This is justified by setting the linear offset constant C1_liq_Brun89 to zero [Brun, 1989]
! liquid water faction
frc_liq = min(0.1_r8, (h2osoi_liq(c_idx,i) / (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i))))
!dr_wet = 1E6_r8*(dtime*(C1_liq_Brun89 + C2_liq_Brun89*(frc_liq**(3))) / (4*SHR_CONST_PI*(snw_rds(c_idx,i)/1E6)**(2)))
!simplified, units of microns:
dr_wet = 1E18_r8*(dtime*(C2_liq_Brun89*(frc_liq**(3))) / (4*SHR_CONST_PI*snw_rds(c_idx,i)**(2)))
dr = dr + dr_wet
!
!********** 3. SNOWAGE SCALING (TURNED OFF BY DEFAULT) *************
!
! Multiply rate of change of effective radius by some constant, xdrdt
if (flg_snoage_scl) then
dr = dr*xdrdt
endif
!
!********** 4. INCREMENT EFFECTIVE RADIUS, ACCOUNTING FOR: ***********
! DRY AGING
! WET AGING
! FRESH SNOW
! RE-FREEZING
!
! new snowfall [kg/m2]
if (do_capsnow(c_idx)) then
newsnow = max(0._r8, (qflx_snwcp_ice(c_idx)*dtime))
else
newsnow = max(0._r8, (qflx_snow_grnd_col(c_idx)*dtime))
endif
! snow that has re-frozen [kg/m2]
refrzsnow = max(0._r8, (qflx_snofrz_lyr(c_idx,i)*dtime))
! fraction of layer mass that is re-frozen
frc_refrz = refrzsnow / h2osno_lyr
! fraction of layer mass that is new snow
if (i == snl_top) then
frc_newsnow = newsnow / h2osno_lyr
else
frc_newsnow = 0._r8
endif
if ((frc_refrz + frc_newsnow) > 1._r8) then
frc_refrz = frc_refrz / (frc_refrz + frc_newsnow)
frc_newsnow = 1._r8 - frc_refrz
frc_oldsnow = 0._r8
else
frc_oldsnow = 1._r8 - frc_refrz - frc_newsnow
endif
! mass-weighted mean of fresh snow, old snow, and re-frozen snow effective radius
snw_rds(c_idx,i) = (snw_rds(c_idx,i)+dr)*frc_oldsnow + snw_rds_min*frc_newsnow + snw_rds_refrz*frc_refrz
!
!********** 5. CHECK BOUNDARIES ***********
!
! boundary check
if (snw_rds(c_idx,i) < snw_rds_min) then
snw_rds(c_idx,i) = snw_rds_min
endif
if (snw_rds(c_idx,i) > snw_rds_max) then
snw_rds(c_idx,i) = snw_rds_max
end if
! set top layer variables for history files
if (i == snl_top) then
snot_top(c_idx) = t_soisno(c_idx,i)
dTdz_top(c_idx) = dTdz(c_idx,i)
snw_rds_top(c_idx) = snw_rds(c_idx,i)
sno_liq_top(c_idx) = h2osoi_liq(c_idx,i) / (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i))
endif
enddo
enddo
! Special case: snow on ground, but not enough to have defined a snow layer:
! set snw_rds to fresh snow grain size:
do fc = 1, num_nosnowc
c_idx = filter_nosnowc(fc)
if (h2osno(c_idx) > 0._r8) then
snw_rds(c_idx,0) = snw_rds_min
endif
enddo
end subroutine SnowAge_grain
subroutine SnowOptics_init( ) 1,67
use fileutils , only : getfil
use CLM_varctl , only : fsnowoptics
use spmdMod , only : mpicom, MPI_REAL8, masterproc
use ncdio , only : nf_close, nf_inq_varid, nf_open, check_ret, &
nf_get_var_double
integer :: ncid ! netCDF file id
character(len=256) :: locfn ! local filename
character(len= 32) :: subname = 'SnowOptics_init' ! subroutine name
integer :: varid ! netCDF id's
integer :: ier ! error status
!
! Open optics file:
if (masterproc) then
write(iulog,*) 'Attempting to read snow optical properties .....'
call getfil (fsnowoptics, locfn, 0)
call check_ret(nf_open(locfn, 0, ncid), subname)
write(iulog,*) subname,trim(fsnowoptics)
! direct-beam snow Mie parameters:
call check_ret(nf_inq_varid(ncid, 'ss_alb_ice_drc', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_snw_drc), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_ice_drc', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_snw_drc), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_ice_drc', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_snw_drc), subname)
! diffuse snow Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_ice_dfs', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_snw_dfs), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_ice_dfs', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_snw_dfs), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_ice_dfs', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_snw_dfs), subname)
! BC species 1 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_bcphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_bc1), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_bcphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_bc1), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_bcphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_bc1), subname)
! BC species 2 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_bcphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_bc2), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_bcphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_bc2), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_bcphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_bc2), subname)
! OC species 1 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_ocphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_oc1), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_ocphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_oc1), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_ocphil', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_oc1), subname)
! OC species 2 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_ocphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_oc2), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_ocphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_oc2), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_ocphob', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_oc2), subname)
! dust species 1 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_dust01', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_dst1), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_dust01', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_dst1), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_dust01', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_dst1), subname)
! dust species 2 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_dust02', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_dst2), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_dust02', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_dst2), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_dust02', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_dst2), subname)
! dust species 3 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_dust03', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_dst3), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_dust03', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_dst3), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_dust03', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_dst3), subname)
! dust species 4 Mie parameters
call check_ret(nf_inq_varid(ncid, 'ss_alb_dust04', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ss_alb_dst4), subname)
call check_ret(nf_inq_varid(ncid, 'asm_prm_dust04', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, asm_prm_dst4), subname)
call check_ret(nf_inq_varid(ncid, 'ext_cff_mss_dust04', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, ext_cff_mss_dst4), subname)
endif
call mpi_bcast (ss_alb_snw_dfs, size(ss_alb_snw_dfs), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_snw_dfs, size(asm_prm_snw_dfs), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_snw_dfs, size(ext_cff_mss_snw_dfs), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_snw_drc, size(ss_alb_snw_drc), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_snw_drc, size(asm_prm_snw_drc), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_snw_drc, size(ext_cff_mss_snw_drc), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_bc1, size(ss_alb_bc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_bc2, size(ss_alb_bc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_oc1, size(ss_alb_oc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_oc2, size(ss_alb_oc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_bc1, size(asm_prm_bc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_bc2, size(asm_prm_bc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_oc1, size(asm_prm_oc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_oc2, size(asm_prm_oc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_dst1, size(ss_alb_dst1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_dst2, size(ss_alb_dst2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_dst3, size(ss_alb_dst3), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ss_alb_dst4, size(ss_alb_dst4), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_dst1, size(asm_prm_dst1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_dst2, size(asm_prm_dst2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_dst3, size(asm_prm_dst3), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (asm_prm_dst4, size(asm_prm_dst4), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_bc1, size(ext_cff_mss_bc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_bc2, size(ext_cff_mss_bc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_oc1, size(ext_cff_mss_oc1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_oc2, size(ext_cff_mss_oc2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_dst1, size(ext_cff_mss_dst1), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_dst2, size(ext_cff_mss_dst2), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_dst3, size(ext_cff_mss_dst3), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (ext_cff_mss_dst4, size(ext_cff_mss_dst4), MPI_REAL8 , 0, mpicom, ier)
if (masterproc) then
call check_ret(nf_close(ncid), subname)
write(iulog,*) 'Successfully read snow optical properties'
! print some diagnostics:
write (iulog,*) 'SNICAR: Mie single scatter albedos for direct-beam ice, rds=100um: ', &
ss_alb_snw_drc(71,1), ss_alb_snw_drc(71,2), ss_alb_snw_drc(71,3), &
ss_alb_snw_drc(71,4), ss_alb_snw_drc(71,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for diffuse ice, rds=100um: ', &
ss_alb_snw_dfs(71,1), ss_alb_snw_dfs(71,2), ss_alb_snw_dfs(71,3), &
ss_alb_snw_dfs(71,4), ss_alb_snw_dfs(71,5)
if (DO_SNO_OC) then
write (iulog,*) 'SNICAR: Including OC aerosols from snow radiative transfer calculations'
else
write (iulog,*) 'SNICAR: Excluding OC aerosols from snow radiative transfer calculations'
endif
write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophillic BC: ', &
ss_alb_bc1(1,1), ss_alb_bc1(1,2), ss_alb_bc1(1,3), ss_alb_bc1(1,4), ss_alb_bc1(1,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophobic BC: ', &
ss_alb_bc2(1,1), ss_alb_bc2(1,2), ss_alb_bc2(1,3), ss_alb_bc2(1,4), ss_alb_bc2(1,5)
if (DO_SNO_OC) then
write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophillic OC: ', &
ss_alb_oc1(1,1), ss_alb_oc1(1,2), ss_alb_oc1(1,3), ss_alb_oc1(1,4), ss_alb_oc1(1,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophobic OC: ', &
ss_alb_oc2(1,1), ss_alb_oc2(1,2), ss_alb_oc2(1,3), ss_alb_oc2(1,4), ss_alb_oc2(1,5)
endif
write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 1: ', &
ss_alb_dst1(1,1), ss_alb_dst1(1,2), ss_alb_dst1(1,3), ss_alb_dst1(1,4), ss_alb_dst1(1,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 2: ', &
ss_alb_dst2(1,1), ss_alb_dst2(1,2), ss_alb_dst2(1,3), ss_alb_dst2(1,4), ss_alb_dst2(1,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 3: ', &
ss_alb_dst3(1,1), ss_alb_dst3(1,2), ss_alb_dst3(1,3), ss_alb_dst3(1,4), ss_alb_dst3(1,5)
write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 4: ', &
ss_alb_dst4(1,1), ss_alb_dst4(1,2), ss_alb_dst4(1,3), ss_alb_dst4(1,4), ss_alb_dst4(1,5)
write(iulog,*)
end if
end subroutine SnowOptics_init
subroutine SnowAge_init( ) 1,13
use CLM_varctl , only : fsnowaging
use fileutils , only : getfil
use spmdMod , only : mpicom, MPI_REAL8, masterproc
use ncdio , only : nf_close, nf_inq_varid, nf_open, check_ret, &
nf_get_var_double
integer :: ncid ! netCDF file id
character(len=256) :: locfn ! local filename
character(len= 32) :: subname = 'SnowOptics_init' ! subroutine name
integer :: varid ! netCDF id's
integer :: ier ! error status
! Open snow aging (effective radius evolution) file:
if (masterproc) then
write(iulog,*) 'Attempting to read snow aging parameters .....'
call getfil (fsnowaging, locfn, 0)
call check_ret(nf_open(locfn, 0, ncid), subname)
write(iulog,*) subname,trim(fsnowaging)
! snow aging parameters
call check_ret(nf_inq_varid(ncid, 'tau', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, snowage_tau), subname)
call check_ret(nf_inq_varid(ncid, 'kappa', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, snowage_kappa), subname)
call check_ret(nf_inq_varid(ncid, 'drdsdt0', varid), subname)
call check_ret(nf_get_var_double(ncid, varid, snowage_drdt0), subname)
endif
call mpi_bcast (snowage_tau, size(snowage_tau), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (snowage_kappa, size(snowage_kappa), MPI_REAL8 , 0, mpicom, ier)
call mpi_bcast (snowage_drdt0, size(snowage_drdt0), MPI_REAL8 , 0, mpicom, ier)
if (masterproc) then
call check_ret(nf_close(ncid), subname)
write(iulog,*) 'Successfully read snow aging properties'
! print some diagnostics:
write (iulog,*) 'SNICAR: snowage tau for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_tau(3,11,9)
write (iulog,*) 'SNICAR: snowage kappa for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_kappa(3,11,9)
write (iulog,*) 'SNICAR: snowage dr/dt_0 for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_drdt0(3,11,9)
endif
end subroutine SnowAge_init
end module SNICARMod