diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml index 4fe7f81874e4..7c9480979f19 100644 --- a/components/elm/bld/namelist_files/namelist_defaults.xml +++ b/components/elm/bld/namelist_files/namelist_defaults.xml @@ -556,6 +556,10 @@ this mask will have smb calculated over the entire global land surface lnd/clm2/snicardata/snicar_optics_5bnd_mam_c160322.nc lnd/clm2/snicardata/snicar_drdt_bst_fit_60_c070416.nc .true. +1 +.false. +.false. +0 2000 diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml index 872feb7d347b..58a8f12a78f3 100644 --- a/components/elm/bld/namelist_files/namelist_definition.xml +++ b/components/elm/bld/namelist_files/namelist_definition.xml @@ -511,7 +511,7 @@ Per tape series history write frequency. +group="elm_inparm" valid_values=""> + +snow shape: 1=sphere (original assumption in SNICAR) + 2=spheroid; 3=hexagonal plate; 4=koch snowflake + + + +is_dust_internal_mixing + + + +is_BC_internal_mixing + + + +atmospheric types for snicar + + Toggle to use 16 snow layers instead of 5. This option enables more realistic diff --git a/components/elm/src/biogeophys/AerosolType.F90 b/components/elm/src/biogeophys/AerosolType.F90 index 14e811ad4082..4f868894b778 100644 --- a/components/elm/src/biogeophys/AerosolType.F90 +++ b/components/elm/src/biogeophys/AerosolType.F90 @@ -231,6 +231,39 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='mass of dust in top snow layer', & ptr_col=this%mss_dst_top_col, set_urb=spval) + ! add ouput of mass of aerosols + this%mss_cnc_bcphi_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_bcphi_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_bcphi_col', & + ptr_col=this%mss_cnc_bcphi_col, no_snow_behavior=no_snow_normal, default='inactive') + + this%mss_cnc_bcpho_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_bcpho_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_bcpho_col', & + ptr_col=this%mss_cnc_bcpho_col, no_snow_behavior=no_snow_normal, default='inactive') + + this%mss_cnc_dst1_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_dst1_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_dst1_col', & + ptr_col=this%mss_cnc_dst1_col, no_snow_behavior=no_snow_normal, default='inactive') + + this%mss_cnc_dst2_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_dst2_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_dst2_col', & + ptr_col=this%mss_cnc_dst2_col, no_snow_behavior=no_snow_normal, default='inactive') + + this%mss_cnc_dst3_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_dst3_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_dst3_col', & + ptr_col=this%mss_cnc_dst3_col, no_snow_behavior=no_snow_normal, default='inactive') + + this%mss_cnc_dst4_col(begc:endc,:) = spval + call hist_addfld2d (fname='mss_cnc_dst4_col', units='kg/kg', type2d='levsno', & + avgflag='A', long_name='mss_cnc_dst4_col', & + ptr_col=this%mss_cnc_dst4_col, no_snow_behavior=no_snow_normal, default='inactive') + + + end subroutine InitHistory !----------------------------------------------------------------------- diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90 index 3c7b569e2f66..4508451593c0 100644 --- a/components/elm/src/biogeophys/SnowSnicarMod.F90 +++ b/components/elm/src/biogeophys/SnowSnicarMod.F90 @@ -32,6 +32,7 @@ module SnowSnicarMod 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: integer, public, parameter :: sno_nbr_aer = 8 ! number of aerosol species in snowpack @@ -122,13 +123,19 @@ module SnowSnicarMod 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); + !!! direct & diffuse + real(r8) :: flx_wgt_dir (6,90,numrad_snw); ! 6 atmospheric types, 0-89 SZA + real(r8) :: flx_wgt_dif (6, numrad_snw); ! 6 atmospheric types + !$acc declare create(ss_alb_snw_drc ) !$acc declare create(asm_prm_snw_drc ) !$acc declare create(ext_cff_mss_snw_drc) !$acc declare create(ss_alb_snw_dfs ) !$acc declare create(asm_prm_snw_dfs ) !$acc declare create(ext_cff_mss_snw_dfs) - + !$acc declare create(flx_wgt_dir ) + !$acc declare create(flx_wgt_dif ) + #ifdef MODAL_AER !mgf++ ! Size-dependent BC optical properties. Currently a fixed BC size is @@ -1482,7 +1489,7 @@ end subroutine SnowAge_grain subroutine SnowOptics_init( ) use fileutils , only : getfil - use elm_varctl , only : fsnowoptics + use elm_varctl , only : fsnowoptics, snicar_atm_type use spmdMod , only : masterproc use ncdio_pio , only : file_desc_t, ncd_io, ncd_pio_openfile, ncd_pio_closefile use ncdio_pio , only : ncd_pio_openfile, ncd_inqfdims, ncd_pio_closefile, ncd_inqdid, ncd_inqdlen @@ -1512,6 +1519,26 @@ subroutine SnowOptics_init( ) call ncd_io( 'ss_alb_ice_dfs', ss_alb_snw_dfs, 'read', ncid, posNOTonfile=.true.) call ncd_io( 'asm_prm_ice_dfs', asm_prm_snw_dfs, 'read', ncid, posNOTonfile=.true.) call ncd_io( 'ext_cff_mss_ice_dfs', ext_cff_mss_snw_dfs, 'read', ncid, posNOTonfile=.true.) + + !!! direct and diffuse flux under different atmospheric conditions + if (snicar_atm_type > 0)then + ! direct-beam incident spectral flux: + call ncd_io( 'flx_wgt_dir', flx_wgt_dir, 'read', ncid, posNOTonfile=.true.) + ! diffuse incident spectral flux: + call ncd_io( 'flx_wgt_dif', flx_wgt_dif, 'read', ncid, posNOTonfile=.true.) + + ! Test by Dalei Hao + !write (iulog,*) 'flx_wgt_dir', flx_wgt_dir(1,10,2), & + ! flx_wgt_dir(1,10,3), flx_wgt_dir(1,10,4), flx_wgt_dir(1,10,5), & + ! flx_wgt_dir(1,20,2), & + ! flx_wgt_dir(1,20,3), flx_wgt_dir(1,20,4), flx_wgt_dir(1,20,5) + + !write (iulog,*) 'flx_wgt_dif', flx_wgt_dif(1,2), & + ! flx_wgt_dif(1,3), flx_wgt_dif(1,4), flx_wgt_dif(1,5), & + ! flx_wgt_dif(2,2), & + ! flx_wgt_dif(2,3), flx_wgt_dif(2,4), flx_wgt_dif(2,5) + endif + !$acc update device( & !$acc ss_alb_snw_drc ,& !$acc asm_prm_snw_drc ,& @@ -1753,6 +1780,7 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & use elm_varpar , only : nlevsno, numrad use clm_time_manager , only : get_nstep use shr_const_mod , only : SHR_CONST_PI + use elm_varctl , only: snow_shape_defined,is_dust_internal_mixing,is_BC_internal_mixing, snicar_atm_type ! ! !ARGUMENTS: integer , intent(in) :: flg_snw_ice ! flag: =1 when called from CLM, =2 when called from CSIM @@ -1857,6 +1885,37 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & integer :: sfctype ! underlying surface type (debugging only) real(r8):: pi ! 3.1415... + !!!!!!!!!!!!!!!!!!!!!!!! + ! New variales for non-spherical snow shape ! He et al., 2017 + integer :: snw_shp_lcl(-nlevsno+1:0) ! Snow grain shape option: + ! 1=sphere; 2=spheroid; 3=hexagonal plate; 4=koch snowflake + integer :: snw_fs_lcl(-nlevsno+1:0) ! Shape factor: ratio of nonspherical grain effective radii to that of equal-volume sphere + ! 0=use recommended default value (He et al. 2017); + ! others(0 1 (i.e. nonspherical) + integer :: snw_ar_lcl(-nlevsno+1:0) ! % Aspect ratio: ratio of grain width to length + ! 0=use recommended default value (He et al. 2017); + ! others(0.1 1 (i.e. nonspherical) + real(r8):: & + diam_ice , & ! + fs_sphd , & ! + fs_hex0 , & ! + fs_hex , & ! + fs_koch , & ! + AR_tmp , & ! + g_ice_Cg_tmp(7) , & ! + gg_ice_F07_tmp(7) , & ! + g_ice_F07 , & ! + g_ice , & ! + gg_F07_intp , & ! + g_Cg_intp, & ! + R_1_omega_tmp , & ! !!! BC internal mixing + C_BC_total , & ! !!! BC concentrtion + C_dust_total !! dust concentration + integer :: slr_zen + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! SNICAR_AD new variables, follow sea-ice shortwave conventions real(r8):: & trndir(-nlevsno+1:1) , & ! solar beam down transmission from top @@ -2002,6 +2061,30 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & !mgf-- #endif + ! Constants for aspherical ice particles %%% + ! g_snw asymmetry factor parameterization coefficients (6 bands) from + ! Table 3 & Eqs. 6-7 in He et al. (2017) + ! assume same values for 4-5 um band, which leads to very small biases (<3%) + + real(r8) :: g_b2(7) + real(r8) :: g_b1(7) + real(r8) :: g_b0(7) + real(r8) :: g_F07_c2(7) + real(r8) :: g_F07_c1(7) + real(r8) :: g_F07_c0(7) + real(r8) :: g_F07_p2(7) + real(r8) :: g_F07_p1(7) + real(r8) :: g_F07_p0(7) + real(r8) :: BC_d0(3) + real(r8) :: BC_d1(3) + real(r8) :: BC_d2(3) + real(r8) :: dust_clear_d0(3) + real(r8) :: dust_clear_d1(3) + real(r8) :: dust_clear_d2(3) + real(r8) :: dust_cloudy_d0(3) + real(r8) :: dust_cloudy_d1(3) + real(r8) :: dust_cloudy_d2(3) + ! Enforce expected array sizes associate(& @@ -2031,8 +2114,42 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & 0.0951585_r8, 0.1246290_r8, & 0.1495960_r8, 0.1691565_r8, & 0.1826034_r8, 0.1894506_r8/) - - + !!!!!!!!!!!! snow albedo improvement by Dalei Hao 2021 + !!!!!!!!! snow shape + ! define snow shape + snw_shp_lcl(:) = snow_shape_defined + snw_fs_lcl(:) = 0._r8 + snw_ar_lcl(:) = 0._r8 + !data g_wvl(:) /0.25,0.70,1.41,1.90,2.50,3.50,4.00,5.00/ ! wavelength (um) division point + !g_wvl_center = g_wvl(2:8)/2 + g_wvl(1:7)/2 ; ! center point for wavelength band + data g_b0(:) /9.76029E-01,9.67798E-01,1.00111E+00,1.00224E+00,9.64295E-01,9.97475E-01,9.97475E-01/ + data g_b1(:) /5.21042E-01,4.96181E-01,1.83711E-01,1.37082E-01,5.50598E-02,8.48743E-02,8.48743E-02/ + data g_b2(:) /-2.66792E-04,1.14088E-03,2.37011E-04,-2.35905E-04,8.40449E-04,-4.71484E-04,-4.71484E-04/ + ! Tables 1 & 2 and Eqs. 3.1-3.4 from Fu, 2007 + data g_F07_c2(:) /1.349959E-1,1.115697E-1,9.853958E-2,5.557793E-2,-1.233493E-1,0.0,0.0/ + data g_F07_c1(:) /-3.987320E-1,-3.723287E-1,-3.924784E-1,-3.259404E-1,4.429054E-2,-1.726586E-1,-1.726586E-1/ + data g_F07_c0(:) /7.938904E-1,8.030084E-1,8.513932E-1,8.692241E-1,7.085850E-1,6.412701E-1,6.412701E-1/ + data g_F07_p2(:) /3.165543E-3,2.014810E-3,1.780838E-3,6.987734E-4,-1.882932E-2,-2.277872E-2,-2.277872E-2/ + data g_F07_p1(:) /1.140557E-1,1.143152E-1,1.143814E-1,1.071238E-1,1.353873E-1,1.914431E-1,1.914431E-1/ + data g_F07_p0(:) /5.292852E-1,5.425909E-1,5.601598E-1,6.023407E-1,6.473899E-1,4.634944E-1,4.634944E-1/ + + !!! BC internal mixing + data BC_d0(:) /3.50098E-02,6.51688E-03,7.96544E-01/ + data BC_d1(:) /9.91050E-01,7.36315E-01,4.36649E-02/ + data BC_d2(:) /3.00370E+01,9.52134E+02,2.57288E+02/ + + !!! dust internal mixing + data dust_clear_d0(:) /1.0413E+00,1.0168E+00,1.0189E+00/ + data dust_clear_d1(:) /1.0016E+00,1.0070E+00,1.0840E+00/ + data dust_clear_d2(:) /2.4208E-01,1.5300E-03,1.1230E-04/ + + data dust_cloudy_d0(:) /1.0388E+00,1.0167E+00,1.0189E+00/ + data dust_cloudy_d1(:) /1.0015E+00,1.0061E+00,1.0823E+00/ + data dust_cloudy_d2(:) /2.5973E-01,1.6200E-03,1.1721E-04/ +!!!!!!!!!!!!! + + + ! Loop over all non-urban columns ! (when called from CSIM, there is only one column) do fc = 1,num_nourbanc @@ -2109,8 +2226,8 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & ! 40nm) assumed for freshly-emitted BC in MAM. Future ! implementations may prognose the BC effective radius in ! snow. - rds_bcint_lcl(:) = 100._r8 - rds_bcext_lcl(:) = 100._r8 + rds_bcint_lcl(:) = 100._r8 ! default 100 + rds_bcext_lcl(:) = 100._r8 ! default 100 !mgf-- #endif @@ -2169,18 +2286,46 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & 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 + if (snicar_atm_type == 0) 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 + else + slr_zen = nint(acosd(coszen(c_idx))) + if (slr_zen>89) then + slr_zen = 89 + endif + flx_wgt(1) = 1._r8 + flx_wgt(2) = flx_wgt_dir(snicar_atm_type, slr_zen+1, 2) + flx_wgt(3) = flx_wgt_dir(snicar_atm_type, slr_zen+1, 3) + flx_wgt(4) = flx_wgt_dir(snicar_atm_type, slr_zen+1, 4) + flx_wgt(5) = flx_wgt_dir(snicar_atm_type, slr_zen+1, 5) + + ! write(iulog,*) "SNICAR_AD STATS: coszen(c_idx) (0)= ", coszen(c_idx) ! add by Dalei check + ! write(iulog,*) "SNICAR_AD STATS: slr_zen (0)= ", slr_zen ! add by Dalei check + ! write(iulog,*) "SNICAR_AD STATS: flx_wgt (2)= ", flx_wgt(2) ! add by Dalei check + ! write(iulog,*) "SNICAR_AD STATS: flx_wgt (4)= ", flx_wgt(4) ! add by Dalei check + endif ! 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 + if (snicar_atm_type == 0) 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 + else + flx_wgt(1) = 1._r8 + flx_wgt(2) = flx_wgt_dif(snicar_atm_type, 2) + flx_wgt(3) = flx_wgt_dif(snicar_atm_type, 3) + flx_wgt(4) = flx_wgt_dif(snicar_atm_type, 4) + flx_wgt(5) = flx_wgt_dif(snicar_atm_type, 5) + + !write(iulog,*) "SNICAR_AD STATS: flx_wgt (0)= ", flx_wgt(2) ! add by Dalei check + !write(iulog,*) "SNICAR_AD STATS: flx_wgt (0)= ", flx_wgt(4) ! add by Dalei check + endif endif endif ! end if numrad_snw @@ -2239,6 +2384,102 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & enddo endif + !!! Dalei Hao + ! shape-dependent asymetry factors (He et al., 2017) + do i=snl_top,snl_btm,1 + if(snw_shp_lcl(i) == 2) then ! spheroid + + diam_ice = 2._r8*snw_rds_lcl(i) + if(snw_fs_lcl(i) == 0) then + fs_sphd = 0.929_r8 + else + fs_sphd = snw_fs_lcl(i) + endif + fs_hex = 0.788_r8 + if(snw_ar_lcl(i) == 0) then + AR_tmp = 0.5_r8 + else + AR_tmp = snw_ar_lcl(i) + endif + g_ice_Cg_tmp = g_b0 * ((fs_sphd/fs_hex)**g_b1) * (diam_ice**g_b2) ! Eq.7, He et al. (2017) + gg_ice_F07_tmp = g_F07_c0 + g_F07_c1 * AR_tmp + g_F07_c2 * (AR_tmp**2) ! Eqn. 3.1 in Fu (2007) + + elseif(snw_shp_lcl(i) == 3) then ! hexagonal plate + diam_ice = 2._r8*snw_rds_lcl(i) + if(snw_fs_lcl(i) == 0) then + fs_hex0 = 0.788_r8 + else + fs_hex0 = snw_fs_lcl(i) + endif + fs_hex = 0.788_r8 + if(snw_ar_lcl(i) == 0) then + AR_tmp = 2.5_r8 + else + AR_tmp = snw_ar_lcl(i) + endif + g_ice_Cg_tmp = g_b0 * ((fs_hex0/fs_hex)**g_b1) * (diam_ice**g_b2) ! Eq.7, He et al. (2017) + gg_ice_F07_tmp = g_F07_p0 + g_F07_p1 * log(AR_tmp) + g_F07_p2 * ((log(AR_tmp))**2) ! Eqn. 3.3 in Fu (2007) + + elseif(snw_shp_lcl(i) == 4) then ! koch snowflake + diam_ice = 2._r8 * snw_rds_lcl(i) /0.544_r8 + if(snw_fs_lcl(i) == 0) then + fs_koch = 0.712_r8 + else + fs_koch = snw_fs_lcl(i) + endif + fs_hex = 0.788_r8 + if(snw_ar_lcl(i) == 0) then + AR_tmp = 2.5_r8 + else + AR_tmp = snw_ar_lcl(i) + endif + + g_ice_Cg_tmp = g_b0 * ((fs_koch/fs_hex)**g_b1) * (diam_ice**g_b2) ! Eq.7, He et al. (2017) + gg_ice_F07_tmp = g_F07_p0 + g_F07_p1 * log(AR_tmp) + g_F07_p2 * ((log(AR_tmp))**2) ! Eqn. 3.3 in Fu (2007) + + endif + + ! 6 wavelength bands for g_ice to be interpolated into 480-bands of SNICAR + ! shape-preserving piecewise interpolation into 480-bands + if(snw_shp_lcl(i) > 1) then + !g_Cg_intp = pchip(g_wvl_center,g_ice_Cg_tmp,wvl) ; + !gg_F07_intp = pchip(g_wvl_center,gg_ice_F07_tmp,wvl) ; + !data g_wvl(:) /0.25,0.70,1.41,1.90,2.50,3.50,4.00,5.00/ ! wavelength (um) division point + !g_wvl_center = g_wvl(2:8)/2 + g_wvl(1:7)/2 ; ! center point for wavelength band + ! elm wavelength/ /0.3,0.7,1.0,1.2,1.5,5/ + !wvl_5bd = [0.5 0.85 1.1 1.35 3.25]; + ! He /0.475 1.055 1.655 2.2 3 3.75 4.5 + ! linear interpolation to get the Cg and G_f80 for band_idx. + if(bnd_idx == 1) then + g_Cg_intp = (g_ice_Cg_tmp(2)-g_ice_Cg_tmp(1))/(1.055_r8-0.475_r8)*(0.5_r8-0.475_r8)+g_ice_Cg_tmp(1); + gg_F07_intp = (gg_ice_F07_tmp(2)-gg_ice_F07_tmp(1))/(1.055_r8-0.475_r8)*(0.5_r8-0.475_r8)+gg_ice_F07_tmp(1); + elseif(bnd_idx == 2) then + g_Cg_intp = (g_ice_Cg_tmp(2)-g_ice_Cg_tmp(1))/(1.055_r8-0.475_r8)*(0.85_r8-0.475_r8)+g_ice_Cg_tmp(1); + gg_F07_intp = (gg_ice_F07_tmp(2)-gg_ice_F07_tmp(1))/(1.055_r8-0.475_r8)*(0.85_r8-0.475_r8)+gg_ice_F07_tmp(1); + + elseif(bnd_idx == 3) then + g_Cg_intp = (g_ice_Cg_tmp(3)-g_ice_Cg_tmp(2))/(1.655_r8-1.055_r8)*(1.1_r8-1.055_r8)+g_ice_Cg_tmp(2); + gg_F07_intp = (gg_ice_F07_tmp(3)-gg_ice_F07_tmp(2))/(1.655_r8-1.055_r8)*(1.1_r8-1.055_r8)+gg_ice_F07_tmp(2); + elseif(bnd_idx == 4) then + g_Cg_intp = (g_ice_Cg_tmp(3)-g_ice_Cg_tmp(2))/(1.655_r8-1.055_r8)*(1.35_r8-1.055_r8)+g_ice_Cg_tmp(2); + gg_F07_intp = (gg_ice_F07_tmp(3)-gg_ice_F07_tmp(2))/(1.655_r8-1.055_r8)*(1.35_r8-1.055_r8)+gg_ice_F07_tmp(2); + elseif(bnd_idx == 5) then + g_Cg_intp = (g_ice_Cg_tmp(6)-g_ice_Cg_tmp(5))/(3.75_r8-3.0_r8)*(3.25_r8-3.0_r8)+g_ice_Cg_tmp(5); + gg_F07_intp = (gg_ice_F07_tmp(6)-gg_ice_F07_tmp(5))/(3.75_r8-3.0_r8)*(3.25_r8-3.0_r8)+gg_ice_F07_tmp(5); + endif + + g_ice_F07 = gg_F07_intp + (1._r8 - gg_F07_intp) / ss_alb_snw_lcl(i) / 2._r8 ! Eq.2.2 in Fu (2007) + g_ice = g_ice_F07 * g_Cg_intp ! Eq.6, He et al. (2017) + asm_prm_snw_lcl(i) = g_ice; + endif + + if(asm_prm_snw_lcl(i) > 0.99_r8) then + asm_prm_snw_lcl(i) = 0.99_r8 + endif + + enddo + !!!-end + !H. Wang ! aerosol species 1 optical properties ! ss_alb_aer_lcl(1) = ss_alb_bc1(bnd_idx) @@ -2323,16 +2564,27 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & ! retrieve absorption enhancement factor for within-ice BC enh_fct = bcenh(bnd_idx,idx_bcint_nclrds,idx_bcint_icerds) + if (is_BC_internal_mixing) then ! get BC optical properties (moved from above) ! aerosol species 1 optical properties (within-ice BC) - ss_alb_aer_lcl(1) = ss_alb_bc1(bnd_idx,idx_bcint_nclrds) - asm_prm_aer_lcl(1) = asm_prm_bc1(bnd_idx,idx_bcint_nclrds) - ext_cff_mss_aer_lcl(1) = ext_cff_mss_bc1(bnd_idx,idx_bcint_nclrds)*enh_fct + ss_alb_aer_lcl(1) = ss_alb_bc1(bnd_idx,idx_bcint_nclrds) + asm_prm_aer_lcl(1) = asm_prm_bc1(bnd_idx,idx_bcint_nclrds) + ext_cff_mss_aer_lcl(1) = ext_cff_mss_bc1(bnd_idx,idx_bcint_nclrds)*enh_fct ! aerosol species 2 optical properties (external BC) - ss_alb_aer_lcl(2) = ss_alb_bc2(bnd_idx,idx_bcext_nclrds) - asm_prm_aer_lcl(2) = asm_prm_bc2(bnd_idx,idx_bcext_nclrds) - ext_cff_mss_aer_lcl(2) = ext_cff_mss_bc2(bnd_idx,idx_bcext_nclrds) + ss_alb_aer_lcl(2) = ss_alb_bc2(bnd_idx,idx_bcext_nclrds) + asm_prm_aer_lcl(2) = asm_prm_bc2(bnd_idx,idx_bcext_nclrds) + ext_cff_mss_aer_lcl(2) = ext_cff_mss_bc2(bnd_idx,idx_bcext_nclrds) + else + ss_alb_aer_lcl(1) = ss_alb_bc1(bnd_idx,idx_bcint_nclrds) + asm_prm_aer_lcl(1) = asm_prm_bc1(bnd_idx,idx_bcint_nclrds) + ext_cff_mss_aer_lcl(1) = ext_cff_mss_bc1(bnd_idx,idx_bcint_nclrds) + + ! aerosol species 2 optical properties (external BC) + ss_alb_aer_lcl(2) = ss_alb_bc2(bnd_idx,idx_bcext_nclrds) + asm_prm_aer_lcl(2) = asm_prm_bc2(bnd_idx,idx_bcext_nclrds) + ext_cff_mss_aer_lcl(2) = ext_cff_mss_bc2(bnd_idx,idx_bcext_nclrds) + endif #else ! bulk aerosol treatment (BC optical properties independent @@ -2346,14 +2598,83 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & ss_alb_aer_lcl(2) = ss_alb_bc2(bnd_idx) asm_prm_aer_lcl(2) = asm_prm_bc2(bnd_idx) ext_cff_mss_aer_lcl(2) = ext_cff_mss_bc2(bnd_idx) + + + if (is_BC_internal_mixing) then + + if(bnd_idx < 4) then + ! R_1_omega: BC-induced enhancement in snow single-scattering coalbedo + ! R_1_omega from Eq.8b in He et al.(2017,JC) is based on BC Re=0.1um & + ! MAC=6.81 m2/g (@550 nm) & BC density=1.7g/cm3. + ! To be consistent with SNICAR default (BC MAC=7.5 m2/g @550nm), we + ! made adjustments on BC size & density as follows to get MAC=7.5m2/g. + ! (1) We use BC Re=0.045um [geometric mean diameter=0.06um (Dentener et al.2006, + ! Yu and Luo,2009) & geometric std=1.5 (Flanner et al.2007;Aoki et al., 2011)] + ! (2) We tune BC density from 1.7 to 1.49 g/cm3 (Aoki et al., 2011) to match BC MAC=7.5 m2/g @550 nm. 100 + C_BC_total = mss_cnc_aer_lcl(i,1) * 1.7_r8/1.49_r8 * 1.0E+09; ! kg/kg to ng/g + + if (C_BC_total > 0) then + R_1_omega_tmp = BC_d0(bnd_idx) * ((C_BC_total + BC_d2(bnd_idx))**BC_d1(bnd_idx)) ! Eq. 8b in He et al.2017,JC + ! Adust R_1_omega_tmp due to BC Re from 0.1 to 0.045um based on + ! Eq. 1 & Table S1 in He et al.2018 (GRL) + if(bnd_idx == 1) then + R_1_omega_tmp = (R_1_omega_tmp / ((0.1_r8/0.05_r8)**(-0.1866_r8)))**((0.1_r8/0.05_r8)**0.1918_r8) ! visible + R_1_omega_tmp = ((0.045_r8/0.05_r8)**(-0.1866_r8)) * (R_1_omega_tmp ** ((0.045_r8/0.05_r8)**(-0.1918_r8))) ! visible + else + R_1_omega_tmp = (R_1_omega_tmp / ((0.1_r8/0.05_r8)**(-0.0046_r8)))** ((0.1_r8/0.05_r8)**0.5177_r8) ! NIR + R_1_omega_tmp = ((0.045_r8/0.05_r8)**(-0.0046_r8)) * (R_1_omega_tmp ** ((0.045_r8/0.05_r8)**(-0.5177_r8))) ! NIR + endif + ! new omega for entire BC-snow internal mixture + ss_alb_snw_lcl(i) = 1_r8 - (1.0_r8 - ss_alb_snw_lcl(i))*R_1_omega_tmp + + endif + endif + + ss_alb_aer_lcl(1) = 0._r8 + asm_prm_aer_lcl(1) = 0._r8 + ext_cff_mss_aer_lcl(1) = 0._r8 + !mss_cnc_aer_lcl(i,1) = 0._r8 + endif #endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! dust internal mixing + if (is_dust_internal_mixing) then + if (bnd_idx < 4) then + C_dust_total = mss_cnc_aer_lcl(i,5) + mss_cnc_aer_lcl(i,6) + mss_cnc_aer_lcl(i,7) + mss_cnc_aer_lcl(i,8) + C_dust_total = C_dust_total * 1.0E+06 ! kg/kg to ug/g + if(C_dust_total > 0) then + ! Direct: + if (flg_slr_in == 1) then + R_1_omega_tmp = dust_clear_d0(bnd_idx) + dust_clear_d2(bnd_idx)*(C_dust_total**dust_clear_d1(bnd_idx)) ! Eq. 1 in He et al.2019,JAMES + else + R_1_omega_tmp = dust_cloudy_d0(bnd_idx) + dust_cloudy_d2(bnd_idx)*(C_dust_total**dust_cloudy_d1(bnd_idx)) ! Eq. 1 in He et al.2019,JAMES + endif + + + ! new omega for entire BC-snow internal mixture + ss_alb_snw_lcl(i) = 1.0_r8 - (1.0_r8 - ss_alb_snw_lcl(i)) *R_1_omega_tmp + + endif + endif + do j = 5,8,1 + ss_alb_aer_lcl(j) = 0._r8 + asm_prm_aer_lcl(j) = 0._r8 + ext_cff_mss_aer_lcl(j) = 0._r8 + !mss_cnc_aer_lcl(i,j) = 0._r8 + enddo + endif + !mgf-- 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) + if (is_dust_internal_mixing .and. (j >= 5)) then + L_aer(i,j) = 0._r8 + else + L_aer(i,j) = L_snw(i)*mss_cnc_aer_lcl(i,j) + endif tau_aer(i,j) = L_aer(i,j)*ext_cff_mss_aer_lcl(j) enddo @@ -2639,6 +2960,29 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & write(iulog,*) "SNICAR_AD STATS: dust2(0)= ", mss_cnc_aer_lcl(0,4) write(iulog,*) "SNICAR_AD STATS: dust3(0)= ", mss_cnc_aer_lcl(0,5) write(iulog,*) "SNICAR_AD STATS: dust4(0)= ", mss_cnc_aer_lcl(0,6) + write(iulog,*) "SNICAR_AD STATS: ss_alb_snw_lcl (0)= ", ss_alb_snw_lcl(i) ! add by Dalei check + write(iulog,*) "SNICAR_AD STATS: asm_prm_snw_lcl (0)= ", asm_prm_snw_lcl(i) ! add by Dalei check + write(iulog,*) "SNICAR_AD STATS: ext_cff_mss_snw_lcl (0)= ", ext_cff_mss_snw_lcl(i) ! add by Dalei check + write(iulog,*) "SNICAR_AD STATS: g_star (0)= ", g_star(i) ! add by Dalei check + write(iulog,*) "SNICAR_AD STATS: omega_star (0)= ", omega_star(i) ! add by Dalei check + write(iulog,*) "SNICAR_AD STATS: tau_star (0)= ", tau_star(i) ! add by Dalei check + write(iulog,*) "g_Cg_intp ", g_Cg_intp ! add by Dalei check + write(iulog,*) "gg_F07_intp ", gg_F07_intp ! add by Dalei check + write(iulog,*) "g_ice_F07 ", g_ice_F07 ! add by Dalei check + write(iulog,*) "bnd_idx ", bnd_idx ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(1) ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(2) ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(3) ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(4) ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(5) ! add by Dalei check + write(iulog,*) "gg_ice_F07_tmp ", gg_ice_F07_tmp(6) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(1) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(2) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(3) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(4) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(5) ! add by Dalei check + write(iulog,*) "g_ice_Cg_tmp ", g_ice_Cg_tmp(6) ! add by Dalei check + call endrun(decomp_index=c_idx, elmlevel=namec, msg=errmsg(__FILE__, __LINE__)) endif enddo diff --git a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 index 0ed3cebb10ad..20154a162a2b 100644 --- a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 +++ b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 @@ -94,6 +94,7 @@ subroutine SurfaceAlbedo(bounds, & ! !USES: !$acc routine seq use elm_varctl , only : iulog, subgridflag, use_snicar_frc, use_fates, use_snicar_ad + use elm_varcon , only: spval ! added by Dalei Hao use shr_orb_mod ! @@ -232,7 +233,15 @@ subroutine SurfaceAlbedo(bounds, & fabd_sun_z => surfalb_vars%fabd_sun_z_patch , & ! Output: [real(r8) (:,:) ] absorbed sunlit leaf direct PAR (per unit lai+sai) for each canopy layer fabd_sha_z => surfalb_vars%fabd_sha_z_patch , & ! Output: [real(r8) (:,:) ] absorbed shaded leaf direct PAR (per unit lai+sai) for each canopy layer fabi_sun_z => surfalb_vars%fabi_sun_z_patch , & ! Output: [real(r8) (:,:) ] absorbed sunlit leaf diffuse PAR (per unit lai+sai) for each canopy layer - fabi_sha_z => surfalb_vars%fabi_sha_z_patch & ! Output: [real(r8) (:,:) ] absorbed shaded leaf diffuse PAR (per unit lai+sai) for each canopy layer + fabi_sha_z => surfalb_vars%fabi_sha_z_patch , & ! Output: [real(r8) (:,:) ] absorbed shaded leaf diffuse PAR (per unit lai+sai) for each canopy layer + albsnd_pur_hst => surfalb_vars%albsnd_pur_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo direct added by Dalei Hao + albsni_pur_hst => surfalb_vars%albsni_pur_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo diffuse added by Dalei Hao + albsnd_nodust_hst => surfalb_vars%albsnd_nodust_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo direct added by Dalei Hao + albsni_nodust_hst => surfalb_vars%albsni_nodust_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo diffuse added by Dalei Hao + albsnd_nobc_hst => surfalb_vars%albsnd_nobc_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo direct added by Dalei Hao + albsni_nobc_hst => surfalb_vars%albsni_nobc_hst_col , & ! Output: [real(r8) (:,:) ] pure snow albedo diffuse added by Dalei Hao + albsnd_hst2 => surfalb_vars%albsnd_hst_col2 , & ! Output: [real(r8) (:,:) ] snow albedo, direct, for history files (col,bnd) [frc] + albsni_hst2 => surfalb_vars%albsni_hst_col2 & ! Output: [real(r8) (:,:) ] snow ground albedo, diffuse, for history files (col,bnd) [frc] ) ! Cosine solar zenith angle for next time step @@ -679,6 +688,7 @@ subroutine SurfaceAlbedo(bounds, & ! pure snow albedo for all-aerosol radiative forcing albgrd_pur(c,ib) = albsod(c,ib)*(1.-frac_sno(c)) + albsnd_pur(c,ib)*frac_sno(c) albgri_pur(c,ib) = albsoi(c,ib)*(1.-frac_sno(c)) + albsni_pur(c,ib)*frac_sno(c) + end if ! also in this loop (but optionally in a different loop for vectorized code) @@ -713,7 +723,7 @@ subroutine SurfaceAlbedo(bounds, & ! For diagnostics, set snow albedo to spval over non-snow non-urban points ! so that it is not averaged in history buffer (OPTIONAL) - ! TODO - this is set to 0 not spval - seems wrong since it will be averaged in + ! TODO - this is set to 0 not spval - seems wrong since it will be averaged in % check do ib = 1, nband do fc = 1,num_nourbanc @@ -721,9 +731,28 @@ subroutine SurfaceAlbedo(bounds, & if ((coszen_col(c) > 0._r8) .and. (h2osno(c) > 0._r8)) then albsnd_hst(c,ib) = albsnd(c,ib) albsni_hst(c,ib) = albsni(c,ib) + ! output pure snow albedo added by Dalei Hao + albsnd_hst2(c,ib) = albsnd(c,ib) + albsni_hst2(c,ib) = albsni(c,ib) + albsnd_pur_hst(c,ib) = albsnd_pur(c,ib) + albsni_pur_hst(c,ib) = albsni_pur(c,ib) + albsnd_nobc_hst(c,ib) = albsnd_bc(c,ib) + albsni_nobc_hst(c,ib) = albsni_bc(c,ib) + albsnd_nodust_hst(c,ib) = albsnd_dst(c,ib) + albsni_nodust_hst(c,ib) = albsni_dst(c,ib) else albsnd_hst(c,ib) = 0._r8 albsni_hst(c,ib) = 0._r8 + ! output pure snow albedo added by Dalei Hao + albsnd_hst2(c,ib) = spval + albsni_hst2(c,ib) = spval + albsnd_pur_hst(c,ib) = spval + albsni_pur_hst(c,ib) = spval + albsnd_nobc_hst(c,ib) = spval + albsni_nobc_hst(c,ib) = spval + albsnd_nodust_hst(c,ib) = spval + albsni_nodust_hst(c,ib) = spval + endif enddo enddo @@ -966,6 +995,7 @@ subroutine SurfaceAlbedo(bounds, & ftii(p,ib) = 1._r8 albd(p,ib) = albgrd(c,ib) albi(p,ib) = albgri(c,ib) + end do end do diff --git a/components/elm/src/biogeophys/SurfaceAlbedoType.F90 b/components/elm/src/biogeophys/SurfaceAlbedoType.F90 index 34d9e87216a8..cea7a9aa4a7a 100644 --- a/components/elm/src/biogeophys/SurfaceAlbedoType.F90 +++ b/components/elm/src/biogeophys/SurfaceAlbedoType.F90 @@ -81,6 +81,22 @@ module SurfaceAlbedoType real(r8), pointer :: fi_top_adjust (:,:) => null() !adjustment factor for diffuse flux (numrad) !!=== End + !!! total broadband albedo added by Dalei Hao + real(r8), pointer :: albsn_hst_col (:) => null() ! col snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsn_pur_hst_col (:) => null() ! col pure snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + real(r8), pointer :: alb_hst_patch (:) => null() ! col patch surface albedo, total, for history files (col) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsn_nodust_hst_col (:) => null() ! col snow albedo (no dust), total, for history files (col) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsn_nobc_hst_col (:) => null() ! col snow albedo (no BC), total, for history files (col) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsnd_pur_hst_col (:,:) => null() ! col pure snow albedo, direct, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsni_pur_hst_col (:,:) => null() ! col pure snow albedo, diffuse, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsnd_nodust_hst_col (:,:) => null() ! col no-dust snow albedo, direct, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsni_nodust_hst_col (:,:) => null() ! col no-dust snow albedo, diffuse, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsnd_nobc_hst_col (:,:) => null() ! col no-bc snow albedo, direct, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsni_nobc_hst_col (:,:) => null() ! col no-bc snow albedo, diffuse, for history files (col,bnd) [frc] ! Added by Dalei Hao + real(r8), pointer :: albsnd_hst_col2 (:,:) => null() ! col snow albedo, direct , for history files (col,bnd) [frc] + real(r8), pointer :: albsni_hst_col2 (:,:) => null() ! col snow albedo, diffuse, for history files (col,bnd) [frc] + + real(r8), pointer :: ftdd_patch (:,:) => null() ! patch down direct flux below canopy per unit direct flx (numrad) real(r8), pointer :: ftid_patch (:,:) => null() ! patch down diffuse flux below canopy per unit direct flx (numrad) real(r8), pointer :: ftii_patch (:,:) => null() ! patch down diffuse flux below canopy per unit diffuse flx (numrad) @@ -301,6 +317,21 @@ subroutine InitAllocate(this, bounds) allocate(this%vcmaxcintsun_patch (begp:endp)) ; this%vcmaxcintsun_patch (:) =spval allocate(this%vcmaxcintsha_patch (begp:endp)) ; this%vcmaxcintsha_patch (:) =spval + ! output snow albedo braodband added by Daleihao + allocate(this%albsn_hst_col (begc:endc)) ; this%albsn_hst_col (:) = spval + allocate(this%albsn_pur_hst_col (begc:endc)) ; this%albsn_pur_hst_col (:) = spval + allocate(this%albsn_nodust_hst_col (begc:endc)) ; this%albsn_nodust_hst_col (:) = spval + allocate(this%albsn_nobc_hst_col (begc:endc)) ; this%albsn_nobc_hst_col (:) = spval + allocate(this%alb_hst_patch (begp:endp)) ; this%alb_hst_patch (:) = spval + allocate(this%albsnd_pur_hst_col (begc:endc,numrad)) ; this%albsnd_pur_hst_col (:,:) = spval + allocate(this%albsni_pur_hst_col (begc:endc,numrad)) ; this%albsni_pur_hst_col (:,:) = spval + allocate(this%albsnd_nodust_hst_col (begc:endc,numrad)) ; this%albsnd_nodust_hst_col (:,:) = spval + allocate(this%albsni_nodust_hst_col (begc:endc,numrad)) ; this%albsni_nodust_hst_col (:,:) = spval + allocate(this%albsnd_nobc_hst_col (begc:endc,numrad)) ; this%albsnd_nobc_hst_col (:,:) = spval + allocate(this%albsni_nobc_hst_col (begc:endc,numrad)) ; this%albsni_nobc_hst_col (:,:) = spval + allocate(this%albsnd_hst_col2 (begc:endc,numrad)) ; this%albsnd_hst_col2 (:,:) = spval + allocate(this%albsni_hst_col2 (begc:endc,numrad)) ; this%albsni_hst_col2 (:,:) = spval + end subroutine InitAllocate !----------------------------------------------------------------------- @@ -362,6 +393,74 @@ subroutine InitHistory(this, bounds) ptr_patch=this%fi_top_adjust) !!=== End + + ! output broadband albedo added by Dalei Hao + this%albsn_hst_col(begc:endc) = spval + call hist_addfld1d (fname='ALBSN', units='1', & + avgflag='A', long_name='broadband snow albedo (total)', & + ptr_col=this%albsn_hst_col, default='inactive') + + this%albsn_pur_hst_col(begc:endc) = spval + call hist_addfld1d (fname='ALBSN_PUR', units='1', & + avgflag='A', long_name='broadband pure snow albedo (total)', & + ptr_col=this%albsn_pur_hst_col, default='inactive') + + this%albsn_nodust_hst_col(begc:endc) = spval + call hist_addfld1d (fname='ALBSN_NODUST', units='1', & + avgflag='A', long_name='broadband no-dust snow albedo (total)', & + ptr_col=this%albsn_nodust_hst_col, default='inactive') + + this%albsn_nobc_hst_col(begc:endc) = spval + call hist_addfld1d (fname='ALBSN_NOBC', units='1', & + avgflag='A', long_name='broadband no-BC snow albedo (total)', & + ptr_col=this%albsn_nobc_hst_col, default='inactive') + + this%alb_hst_patch(begp:endp) = spval + call hist_addfld1d (fname='ALB', units='1', & + avgflag='A', long_name='broadband surface albedo (total)', & + ptr_patch=this%alb_hst_patch, default='inactive') + + this%albsni_hst_col2(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSNI', units='proportion', type2d='numrad', & + avgflag='A', long_name='snow albedo (indirect)', & + ptr_col=this%albsni_hst_col2, default='inactive') + + this%albsni_pur_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSNI_PUR', units='proportion', type2d='numrad', & + avgflag='A', long_name='pure snow albedo (indirect)', & + ptr_col=this%albsni_pur_hst_col, default='inactive') + + this%albsni_nodust_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSNI_NODUST', units='proportion', type2d='numrad', & + avgflag='A', long_name='no-dust snow albedo (indirect)', & + ptr_col=this%albsni_nodust_hst_col, default='inactive') + + this%albsni_nobc_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSNI_NOBC', units='proportion', type2d='numrad', & + avgflag='A', long_name='no-bc snow albedo (indirect)', & + ptr_col=this%albsni_nobc_hst_col, default='inactive') + + this%albsnd_hst_col2(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSND', units='proportion', type2d='numrad', & + avgflag='A', long_name='snow albedo (direct)', & + ptr_col=this%albsnd_hst_col2, default='inactive') + + this%albsnd_pur_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSND_PUR', units='proportion', type2d='numrad', & + avgflag='A', long_name='pure snow albedo (direct)', & + ptr_col=this%albsnd_pur_hst_col, default='inactive') + + this%albsnd_nodust_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSND_NODUST', units='proportion', type2d='numrad', & + avgflag='A', long_name='no-dust snow albedo (direct)', & + ptr_col=this%albsnd_nodust_hst_col, default='inactive') + + this%albsnd_nobc_hst_col(begc:endc,:) = spval + call hist_addfld2d (fname='ALBSND_NOBC', units='proportion', type2d='numrad', & + avgflag='A', long_name='no-bc snow albedo (direct)', & + ptr_col=this%albsnd_nobc_hst_col, default='inactive') + + end subroutine InitHistory !----------------------------------------------------------------------- @@ -409,10 +508,12 @@ subroutine InitCold(this, bounds) this%ftid_patch (begp:endp, :) = 0.0_r8 this%ftii_patch (begp:endp, :) = 1.0_r8 + !!=== TOP solar radiation parameterization ===== this%fd_top_adjust (begp:endp, :) = 1.0_r8 this%fi_top_adjust (begp:endp, :) = 1.0_r8 !!=== End + end subroutine InitCold @@ -501,6 +602,18 @@ subroutine Restart(this, bounds, ncid, flag, & dim1name='pft', dim2name='levcan', switchdim=.true., & long_name='tlai increment for canopy layer', units='', & interpinic_flag='interp', readvar=readvar, data=this%tlai_z_patch) + + !!! added by Dalei Hao + call restartvar(ncid=ncid, flag=flag, varname='albsn_hst_col', xtype=ncd_double, & + dim1name='column', & + long_name='snow albedo total', units='1', & + interpinic_flag='interp', readvar=readvar, data=this%albsn_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='alb_hst_patch', xtype=ncd_double, & + dim1name='pft', & + long_name='surface albedo total', units='1', & + interpinic_flag='interp', readvar=readvar, data=this%alb_hst_patch) + if (flag=='read' .and. .not. readvar) then if (masterproc) then write(iulog,*) "can't find tlai_z in restart (or initial) file..." @@ -654,7 +767,56 @@ subroutine Restart(this, bounds, ncid, flag, & if (masterproc) write(iulog,*) "Initialize albgri_dst to albgri" this%albgri_dst_col(begc:endc,:) = this%albgri_col(begc:endc,:) end if - + + !!! added by Dalei Hao + call restartvar(ncid=ncid, flag=flag, varname='albsn_pur_hst_col', xtype=ncd_double, & + dim1name='column', & + long_name='pure snow albedo total', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsn_pur_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsn_nodust_hst_col', xtype=ncd_double, & + dim1name='column', & + long_name=' snow albedo no dust total', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsn_nodust_hst_col) + + + call restartvar(ncid=ncid, flag=flag, varname='albsn_nobc_hst_col', xtype=ncd_double, & + dim1name='column', & + long_name=' snow albedo no BC total', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsn_nobc_hst_col) + + + call restartvar(ncid=ncid, flag=flag, varname='albsnd_pur_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='pure snow albedo direct (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsnd_pur_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsni_pur_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='pure snow albedo diffuse (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsni_pur_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsnd_nodust_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='no-dust snow albedo direct (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsnd_nodust_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsni_nodust_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='no-dust snow albedo diffuse (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsni_nodust_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsnd_nobc_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='no-BC snow albedo direct (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsnd_nobc_hst_col) + + call restartvar(ncid=ncid, flag=flag, varname='albsni_nobc_hst_col', xtype=ncd_double, & + dim1name='column', dim2name='numrad', switchdim=.true., & + long_name='no-BC snow albedo diffuse (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albsni_nobc_hst_col) + + end if ! end of if-use_snicar_frc ! patch type physical state variable - fabd diff --git a/components/elm/src/biogeophys/SurfaceRadiationMod.F90 b/components/elm/src/biogeophys/SurfaceRadiationMod.F90 index 47fe244c530e..3ff26c74a00e 100644 --- a/components/elm/src/biogeophys/SurfaceRadiationMod.F90 +++ b/components/elm/src/biogeophys/SurfaceRadiationMod.F90 @@ -470,7 +470,20 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsds_sno_vd => surfrad_vars%fsds_sno_vd_patch , & ! Output: [real(r8) (:) ] incident visible, direct radiation on snow (for history files) (pft) [W/m2] fsds_sno_nd => surfrad_vars%fsds_sno_nd_patch , & ! Output: [real(r8) (:) ] incident near-IR, direct radiation on snow (for history files) (pft) [W/m2] fsds_sno_vi => surfrad_vars%fsds_sno_vi_patch , & ! Output: [real(r8) (:) ] incident visible, diffuse radiation on snow (for history files) (pft) [W/m2] - fsds_sno_ni => surfrad_vars%fsds_sno_ni_patch & ! Output: [real(r8) (:) ] incident near-IR, diffuse radiation on snow (for history files) (pft) [W/m2] + fsds_sno_ni => surfrad_vars%fsds_sno_ni_patch , & ! Output: [real(r8) (:) ] incident near-IR, diffuse radiation on snow (for history files) (pft) [W/m2] + !!! total broadband albedo added by Dalei Hao + albsn_hst => surfalb_vars%albsn_hst_col , & ! col snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + albsn_pur_hst => surfalb_vars%albsn_pur_hst_col , & ! col pure snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + albsn_nodust_hst => surfalb_vars%albsn_nodust_hst_col , & ! col no-dust snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + albsn_nobc_hst => surfalb_vars%albsn_nobc_hst_col , & ! col no-BC snow albedo, total , for history files (col) [frc] ! Added by Dalei Hao + alb_hst => surfalb_vars%alb_hst_patch , & ! col patch surface albedo, total, for history files (col) [frc] ! Added by Dalei Hao + albsnd_pur_hst => surfalb_vars%albsnd_pur_hst_col , & ! col pure snow albedo, direct , for history files (col) [frc] ! Added by Dalei Hao + albsni_pur_hst => surfalb_vars%albsni_pur_hst_col , & ! col pure snow albedo, diffuse , for history files (col) [frc] ! Added by Dalei Hao + albsnd_nodust_hst => surfalb_vars%albsnd_nodust_hst_col , & ! col no-dust snow albedo, direct , for history files (col) [frc] ! Added by Dalei Hao + albsni_nodust_hst => surfalb_vars%albsni_nodust_hst_col , & ! col no-dust snow albedo, diffuse , for history files (col) [frc] ! Added by Dalei Hao + albsnd_nobc_hst => surfalb_vars%albsnd_nobc_hst_col , & ! col no-BC snow albedo, direct , for history files (col) [frc] ! Added by Dalei Hao + albsni_nobc_hst => surfalb_vars%albsni_nobc_hst_col & ! col no-BC snow albedo, diffuse , for history files (col) [frc] ! Added by Dalei Hao + ) dtime = dtime_mod @@ -567,6 +580,7 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & ! Solar radiation absorbed by ground surface without any aerosols absrad_pur = trd(p,ib)*(1._r8-albgrd_pur(c,ib)) + tri(p,ib)*(1._r8-albgri_pur(c,ib)) sabg_pur(p) = sabg_pur(p) + absrad_pur + end if end do ! end of pft loop @@ -727,6 +741,9 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & rnir = albd(p,2)*forc_solad(t,2) + albi(p,2)*forc_solai(t,2) fsr(p) = rvis + rnir + !! total broadband albedo added by Dalei Hao + alb_hst(p) = fsr(p) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + fsds_vis_d(p) = forc_solad(t,1) fsds_nir_d(p) = forc_solad(t,2) fsds_vis_i(p) = forc_solai(t,1) @@ -767,6 +784,20 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_sno_nd(p) = fsds_nir_d(p)*albsnd_hst(c,2) fsr_sno_vi(p) = fsds_vis_i(p)*albsni_hst(c,1) fsr_sno_ni(p) = fsds_nir_i(p)*albsni_hst(c,2) + + !! total broadband snow albedo added by Dalei Hao + albsn_hst(c) = (albsnd_hst(c,1) * forc_solad(t,1) + albsni_hst(c,1) * forc_solai(t,1) + albsnd_hst(c,2) *forc_solad(t,2) + albsni_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + + if (use_snicar_frc) then + !! total broadband pure snow albedo added by Dalei Hao + albsn_pur_hst(c) = (albsnd_pur_hst(c,1) * forc_solad(t,1) + albsni_pur_hst(c,1) * forc_solai(t,1) + albsnd_pur_hst(c,2) *forc_solad(t,2) + albsni_pur_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + !! total broadband no-dust snow albedo added by Dalei Hao + albsn_nodust_hst(c) = (albsnd_nodust_hst(c,1) * forc_solad(t,1) + albsni_nodust_hst(c,1) * forc_solai(t,1) + albsnd_nodust_hst(c,2) *forc_solad(t,2) + albsni_nodust_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + !! total broadband no-BC snow albedo added by Dalei Hao + albsn_nobc_hst(c) = (albsnd_nobc_hst(c,1) * forc_solad(t,1) + albsni_nobc_hst(c,1) * forc_solai(t,1) + albsnd_nobc_hst(c,2) *forc_solad(t,2) + albsni_nobc_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + + endif + else fsds_sno_vd(p) = spval fsds_sno_nd(p) = spval @@ -777,13 +808,21 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_sno_nd(p) = spval fsr_sno_vi(p) = spval fsr_sno_ni(p) = spval + albsn_hst(c) = spval ! Dalei Hao + + if (use_snicar_frc) then + albsn_pur_hst(c) = spval ! Dalei Hao + albsn_nodust_hst(c) = spval ! Dalei Hao + albsn_nobc_hst(c) = spval ! Dalei Hao + endif endif end do do fp = 1,num_urbanp p = filter_urbanp(fp) t = veg_pp%topounit(p) g = veg_pp%gridcell(p) - + c = veg_pp%column(p) ! added by Dalei Hao + local_secp1 = secs + nint((grc_pp%londeg(g)/degpsec)/dtime)*dtime local_secp1 = mod(local_secp1,isecspday) @@ -829,6 +868,21 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_nir_d_ln(p) = spval endif fsr(p) = fsr_vis_d(p) + fsr_nir_d(p) + fsr_vis_i(p) + fsr_nir_i(p) + + !! total broadband albedo added by Dalei Hao + alb_hst(p) = fsr(p) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + !! total broadband snow albedo added by Dalei Hao + albsn_hst(c) = (albsnd_hst(c,1) * forc_solad(t,1) + albsni_hst(c,1) * forc_solai(t,1) + albsnd_hst(c,2) *forc_solad(t,2) + albsni_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + + if (use_snicar_frc) then + !! total broadband pure snow albedo added by Dalei Hao + albsn_pur_hst(c) = (albsnd_pur_hst(c,1) * forc_solad(t,1) + albsni_pur_hst(c,1) * forc_solai(t,1) + albsnd_pur_hst(c,2) *forc_solad(t,2) + albsni_pur_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + !! total broadband no-dust snow albedo added by Dalei Hao + albsn_nodust_hst(c) = (albsnd_nodust_hst(c,1) * forc_solad(t,1) + albsni_nodust_hst(c,1) * forc_solai(t,1) + albsnd_nodust_hst(c,2) *forc_solad(t,2) + albsni_nodust_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + !! total broadband no-BC snow albedo added by Dalei Hao + albsn_nobc_hst(c) = (albsnd_nobc_hst(c,1) * forc_solad(t,1) + albsni_nobc_hst(c,1) * forc_solai(t,1) + albsnd_nobc_hst(c,2) *forc_solad(t,2) + albsni_nobc_hst(c,2) * forc_solai(t,2)) / (forc_solad(t,1) + forc_solai(t,1) + forc_solad(t,2) + forc_solai(t,2)) + + endif end do end associate diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90 index 897d05295864..faef710e88e1 100755 --- a/components/elm/src/main/controlMod.F90 +++ b/components/elm/src/main/controlMod.F90 @@ -53,6 +53,8 @@ module controlMod use elm_varctl , only: const_climate_hist use elm_varctl , only: use_top_solar_rad ! TOP solar radiation parameterization ! + use elm_varctl , only: snow_shape_defined,is_dust_internal_mixing,is_BC_internal_mixing,snicar_atm_type + ! ! !PUBLIC TYPES: implicit none save @@ -309,6 +311,10 @@ subroutine control_init( ) use_top_solar_rad ! End + ! snow shape + namelist /elm_inparm/ & + snow_shape_defined,is_dust_internal_mixing,is_BC_internal_mixing, snicar_atm_type + ! ---------------------------------------------------------------------- ! Default values ! ---------------------------------------------------------------------- @@ -806,8 +812,14 @@ subroutine control_spmd() call mpi_bcast (albice, 2, MPI_REAL8,0, mpicom, ier) call mpi_bcast (more_vertlayers,1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (const_climate_hist, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (use_top_solar_rad, 1, MPI_LOGICAL, 0, mpicom, ier) ! TOP solar radiation parameterization + call mpi_bcast (snow_shape_defined, 1, MPI_LOGICAL, 0, mpicom, ier) ! for snow shape + call mpi_bcast (is_dust_internal_mixing, 1, MPI_LOGICAL, 0, mpicom, ier) ! for snow-dust internal mixing + call mpi_bcast (is_BC_internal_mixing, 1, MPI_LOGICAL, 0, mpicom, ier) ! for snow-BC internal mixing + call mpi_bcast (snicar_atm_type, 1, MPI_LOGICAL, 0, mpicom, ier) ! for atmospheric condition + ! glacier_mec variables call mpi_bcast (create_glacier_mec_landunit, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (maxpatch_glcmec, 1, MPI_INTEGER, 0, mpicom, ier) @@ -942,6 +954,10 @@ subroutine control_print () write(iulog,*) ' two-way irrigation = ', tw_irr write(iulog,*) ' use_snicar_frc = ', use_snicar_frc write(iulog,*) ' use_snicar_ad = ', use_snicar_ad + write(iulog,*) ' snow_shape_defined = ', snow_shape_defined ! for snow shape + write(iulog,*) ' is_dust_internal_mixing = ', is_dust_internal_mixing ! for snow shape + write(iulog,*) ' is_BC_internal_mixing = ', is_BC_internal_mixing ! for snow shape + write(iulog,*) ' snicar_atm_type = ', snicar_atm_type ! for snow shape write(iulog,*) ' use_vancouver = ', use_vancouver write(iulog,*) ' use_mexicocity = ', use_mexicocity write(iulog,*) ' use_noio = ', use_noio diff --git a/components/elm/src/main/elm_varcon.F90 b/components/elm/src/main/elm_varcon.F90 index ae29c8633cfb..bcfd598de7f5 100644 --- a/components/elm/src/main/elm_varcon.F90 +++ b/components/elm/src/main/elm_varcon.F90 @@ -131,7 +131,7 @@ module elm_varcon real(r8) :: ac_wasteheat_factor = 0.0_r8 !wasteheat factor for urban air conditioning (-) real(r8) :: wasteheat_limit = 100._r8 !limit on wasteheat (W/m2) - real(r8) :: h2osno_max = 1000._r8 ! max allowed snow thickness (mm H2O) + real(r8) :: h2osno_max = 2000._r8 ! max allowed snow thickness (mm H2O) real(r8), parameter :: lapse_glcmec = 0.006_r8 ! surface temperature lapse rate (deg m-1) ! Pritchard et al. (GRL, 35, 2008) use 0.006 real(r8), parameter :: glcmec_rain_snow_threshold = SHR_CONST_TKFRZ ! temperature dividing rain & snow in downscaling (K) diff --git a/components/elm/src/main/elm_varctl.F90 b/components/elm/src/main/elm_varctl.F90 index 2a47478f1c0a..14fd201e7665 100644 --- a/components/elm/src/main/elm_varctl.F90 +++ b/components/elm/src/main/elm_varctl.F90 @@ -11,7 +11,7 @@ module elm_varctl ! ! !PUBLIC MEMBER FUNCTIONS: implicit none - public :: elm_varctl_set ! Set variables + public :: elm_varctl_set ! Set variables public :: cnallocate_carbon_only_set public :: cnallocate_carbon_only public :: cnallocate_carbonnitrogen_only_set @@ -498,6 +498,26 @@ module elm_varctl integer, public :: budget_ann = 1 integer, public :: budget_ltann = 1 integer, public :: budget_ltend = 0 + + !---------------------------------------------------------- + ! SNICAR + !---------------------------------------------------------- + integer, public :: snow_shape_defined = 1 + logical, public :: is_dust_internal_mixing = .false. + logical, public :: is_BC_internal_mixing = .false. + integer, public :: snicar_atm_type = 0 !Atmospheric profile used to obtain surface-incident spectral flux distribution and subsequent broadband albedo: ! 0 = default + ! 1 = mid-latitude winter + ! 2 = mid-latitude summer + ! 3 = sub-Arctic winter + ! 4 = sub-Arctic summer + ! 5 = Summit,Greenland (sub-Arctic summer, surface pressure of 796hPa) + ! 6 = High Mountain (summer, surface pressure of 556 hPa) + + !$acc declare copyin(snow_shape_defined) + !$acc declare copyin(is_dust_internal_mixing) + !$acc declare copyin(is_BC_internal_mixing) + !$acc declare copyin(snicar_atm_type) + contains !--------------------------------------------------------------------------- diff --git a/share/streams/shr_dmodel_mod.F90 b/share/streams/shr_dmodel_mod.F90 index 8c2c2c9ce6e9..e2e2532b49f6 100644 --- a/share/streams/shr_dmodel_mod.F90 +++ b/share/streams/shr_dmodel_mod.F90 @@ -672,7 +672,7 @@ subroutine shr_dmodel_readLBUB(stream,pio_subsystem,pio_iotype,pio_iodesc_r8, pi rDateUB = real(mDateUB,R8) + real(mSecUB,R8)/spd call t_stopf(trim(lstr)//'_setup') - if (rDateM < rDateLB .or. rDateM > rDateUB) then + if (rDateM < rDateLB .or. rDateM >= rDateUB) then call t_startf(trim(lstr)//'_fbound') if (my_task == master_task) then ! call shr_stream_findBounds(stream,mDate,mSec, & diff --git a/share/streams/shr_tInterp_mod.F90 b/share/streams/shr_tInterp_mod.F90 index 04e1cb5639d6..139aa2b2740c 100644 --- a/share/streams/shr_tInterp_mod.F90 +++ b/share/streams/shr_tInterp_mod.F90 @@ -305,6 +305,11 @@ subroutine shr_tInterp_getAvgCosz(tavCosz,lonr,latr,ymd1,tod1,ymd2,tod2,eccen,mv tod = tod1 do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB] + !--- get next cosz value for t-avg --- + call shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar) + n = n + ldt + tavCosz = tavCosz + cosz*real(ldt,SHR_KIND_R8) ! add to partial sum + !--- advance to next time in [LB,UB] --- ymd0 = ymd tod0 = tod @@ -312,19 +317,6 @@ subroutine shr_tInterp_getAvgCosz(tavCosz,lonr,latr,ymd1,tod1,ymd2,tod2,eccen,mv call shr_cal_advDateInt(ldt,'seconds',ymd0,tod0,ymd,tod,calendar) call shr_cal_timeSet(reday,ymd,tod,calendar) - if (reday > reday2) then - ymd = ymd2 - tod = tod2 - timeint = reday2-reday0 - call ESMF_TimeIntervalGet(timeint,s_i8=dtsec) - ldt = dtsec - endif - - !--- get next cosz value for t-avg --- - call shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr,calendar) - n = n + ldt - tavCosz = tavCosz + cosz*real(ldt,SHR_KIND_R8) ! add to partial sum - end do tavCosz = tavCosz/real(n,SHR_KIND_R8) ! form t-avg