From c4e0873b447287cd69b56bfafded2d9bb42388aa Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Fri, 16 Apr 2021 18:49:23 +0000 Subject: [PATCH 1/3] Add initialization of water vapor mixing ratio at the surface to lsm_ruc_init. This is needed for MYNN surface layer scheme at the first time step. --- physics/sfc_drv_ruc.F90 | 28 +++++++++++++--- physics/sfc_drv_ruc.meta | 72 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 4 deletions(-) diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90 index 8586737c9..517581c56 100644 --- a/physics/sfc_drv_ruc.F90 +++ b/physics/sfc_drv_ruc.F90 @@ -31,17 +31,18 @@ module lsm_ruc !! \htmlinclude lsm_ruc_init.html !! subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & - flag_restart, flag_init, & + flag_restart, flag_init, con_fvirt, con_rd, & im, lsoil_ruc, lsoil, kice, nlev, & ! in lsm_ruc, lsm, slmsk, stype, vtype, & ! in - tsfc_lnd, tsfc_wat, & ! in + t1, q1, prsl1, tsfc_lnd, tsfc_ice, tsfc_wat, & ! in tg3, smc, slc, stc, fice, min_seaice, & ! in sncovr_lnd, sncovr_ice, snoalb, & ! in facsf, facwf, alvsf, alvwf, alnsf, alnwf, & ! in + sfcqv_lnd, sfcqv_ice, & ! out sfalb_lnd_bck, & ! out + semisbase, semis_lnd, semis_ice, & ! out albdvis_lnd,albdnir_lnd,albivis_lnd,albinir_lnd, & ! out albdvis_ice,albdnir_ice,albivis_ice,albinir_ice, & ! out - semisbase, semis_lnd, semis_ice, & ! out zs, sh2o, smfrkeep, tslb, smois, wetness, & ! out tsice, pores, resid, errmsg, errflg) @@ -56,12 +57,18 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & integer, intent(in) :: kice integer, intent(in) :: nlev integer, intent(in) :: lsm_ruc, lsm + real (kind=kind_phys),intent(in) :: con_fvirt + real (kind=kind_phys),intent(in) :: con_rd real (kind=kind_phys), dimension(im), intent(in) :: slmsk real (kind=kind_phys), dimension(im), intent(in) :: stype real (kind=kind_phys), dimension(im), intent(in) :: vtype + real (kind=kind_phys), dimension(im), intent(in) :: t1 + real (kind=kind_phys), dimension(im), intent(in) :: q1 + real (kind=kind_phys), dimension(im), intent(in) :: prsl1 real (kind=kind_phys), dimension(im), intent(in) :: tsfc_lnd + real (kind=kind_phys), dimension(im), intent(in) :: tsfc_ice real (kind=kind_phys), dimension(im), intent(in) :: tsfc_wat real (kind=kind_phys), dimension(im), intent(in) :: tg3 real (kind=kind_phys), dimension(im), intent(in) :: sncovr_lnd @@ -87,7 +94,8 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & real (kind=kind_phys), dimension(im), intent(inout) :: semis_ice real (kind=kind_phys), dimension(im), intent(inout) :: & albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, & - albdvis_ice, albdnir_ice, albivis_ice, albinir_ice + albdvis_ice, albdnir_ice, albivis_ice, albinir_ice, & + sfcqv_lnd, sfcqv_ice ! --- out real (kind=kind_phys), dimension(:), intent(out) :: zs @@ -102,6 +110,7 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & ! --- local real (kind=kind_phys), dimension(lsoil_ruc) :: dzs real (kind=kind_phys) :: alb_lnd, alb_ice + real (kind=kind_phys) :: q0, qs1, rho integer :: ipr, i, k logical :: debug_print integer, dimension(im) :: soiltyp, vegtype @@ -193,6 +202,17 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, & albivis_ice(i) = alb_ice albinir_ice(i) = alb_ice + if (.not.flag_restart) then + !-- initialize QV mixing ratio at the surface from atm. 1st level + q0 = max(q1(i)/(1.-q1(i)), 1.e-8) ! q1=specific humidity at level 1 (kg/kg) + rho = prsl1(i) / (con_rd*t1(i)*(1.0+con_fvirt*q0)) + qs1 = rslf(prsl1(i),tsfc_lnd(i)) !* qs1=sat. mixing ratio at level 1 (kg/kg) + q0 = min(qs1, q0) + sfcqv_lnd(i) = q0 + qs1 = rslf(prsl1(i),tsfc_ice(i)) + sfcqv_ice(i) = qs1 + endif + enddo ! i call init_soil_depth_3 ( zs , dzs , lsoil_ruc ) diff --git a/physics/sfc_drv_ruc.meta b/physics/sfc_drv_ruc.meta index 8198a3c99..e622b0372 100644 --- a/physics/sfc_drv_ruc.meta +++ b/physics/sfc_drv_ruc.meta @@ -63,6 +63,24 @@ type = logical intent = in optional = F +[con_fvirt] + standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one + long_name = rv/rd - 1 (rv = ideal gas constant for water vapor) + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[con_rd] + standard_name = gas_constant_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent @@ -146,6 +164,33 @@ kind = kind_phys intent = inout optional = F +[t1] + standard_name = air_temperature_at_lowest_model_layer + long_name = mean temperature at lowest model layer + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[q1] + standard_name = water_vapor_specific_humidity_at_lowest_model_layer + long_name = water vapor specific humidity at lowest model layer + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[prsl1] + standard_name = air_pressure_at_lowest_model_layer + long_name = mean pressure at lowest model layer + units = Pa + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [tsfc_lnd] standard_name = surface_skin_temperature long_name = surface skin temperature @@ -155,6 +200,15 @@ kind = kind_phys intent = inout optional = F +[tsfc_ice] + standard_name = surface_skin_temperature_over_ice_interstitial + long_name = surface skin temperature over ice (temporary use as interstitial) + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F [tsfc_wat] standard_name = sea_surface_temperature long_name = sea surface temperature @@ -299,6 +353,24 @@ kind = kind_phys intent = inout optional = F +[sfcqv_lnd] + standard_name = water_vapor_mixing_ratio_at_surface_over_land + long_name = water vapor mixing ratio at surface over land + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out + optional = F +[sfcqv_ice] + standard_name = water_vapor_mixing_ratio_at_surface_over_ice + long_name = water vapor mixing ratio at surface over ice + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out + optional = F [sfalb_lnd_bck] standard_name =surface_snow_free_albedo_over_land long_name = surface snow-free albedo over ice From cbbafc58e888a42642720c6e1de62cc7085e6040 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Fri, 16 Apr 2021 20:58:53 +0000 Subject: [PATCH 2/3] Bug fix in the composite for ialb=2 option. --- physics/radiation_surface.f | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/physics/radiation_surface.f b/physics/radiation_surface.f index 41d647796..8e098b37d 100644 --- a/physics/radiation_surface.f +++ b/physics/radiation_surface.f @@ -647,14 +647,14 @@ subroutine setalb & !-- Composite mean surface albedo from land, open water and !-- ice fractions - sfcalb(i,1) = min(0.99,max(0.01,lsmalbdnir(i)))*fracl(i) & - & + asenb_wat*fraco(i) + asenb_ice*fraci(i) - sfcalb(i,2) = min(0.99,max(0.01,lsmalbinir(i)))*fracl(i) & + sfcalb(i,1) = min(0.99,max(0.01,lsmalbdnir(i)))*fracl(i) & ! direct beam NIR + & + asenb_wat*fraco(i) + asenb_ice*fraci(i) + sfcalb(i,2) = min(0.99,max(0.01,lsmalbinir(i)))*fracl(i) & ! diffuse NIR & + asend_wat*fraco(i) + asend_ice*fraci(i) - sfcalb(i,3) = min(0.99,max(0.01,lsmalbdvis(i)))*fracl(i) & - & + asevb_wat*fraco(i) + asenb_ice*fraci(i) - sfcalb(i,4) = min(0.99,max(0.01,lsmalbivis(i)))*fracl(i) & - & + asevd_wat*fraco(i) + asend_ice*fraci(i) + sfcalb(i,3) = min(0.99,max(0.01,lsmalbdvis(i)))*fracl(i) & ! direct beam visible + & + asevb_wat*fraco(i) + asevb_ice*fraci(i) + sfcalb(i,4) = min(0.99,max(0.01,lsmalbivis(i)))*fracl(i) & ! diffuse visible + & + asevd_wat*fraco(i) + asevd_ice*fraci(i) enddo ! end_do_i_loop From d83c1a154c12094c61758e27cfe802b6b73369bc Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Mon, 19 Apr 2021 18:19:45 +0000 Subject: [PATCH 3/3] Add fractional code to ialb=1 option used with the Noah LSM. --- physics/radiation_surface.f | 205 ++++++++++++++++++------------------ 1 file changed, 100 insertions(+), 105 deletions(-) diff --git a/physics/radiation_surface.f b/physics/radiation_surface.f index 8e098b37d..66911c71c 100644 --- a/physics/radiation_surface.f +++ b/physics/radiation_surface.f @@ -449,118 +449,113 @@ subroutine setalb & do i = 1, IMAX -!> - Calculate snow cover input directly for land model, no -!! conversion needed. + !-- water albedo + asevd_wat = 0.06 + asend_wat = 0.06 + asevb_wat = asevd_wat + asenb_wat = asevd_wat + + ! direct albedo CZA dependence over water + if (fraco(i) > f_zero .and. coszf(i) > 0.0001) then + if (tsknf(i) >= con_t0c) then + asevb_wat = max (asevd_wat, 0.026/(coszf(i)**1.7 + 0.065) & + & + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) & + & * (coszf(i)-f_one)) + asenb_wat = asevb_wat + endif - fsno0 = sncovr(i) ! snow fraction on land - - if (nint(slmsk(i))==0 .and. tsknf(i)>con_tice) fsno0 = f_zero - - if (nint(slmsk(i)) == 2) then - if(lsm == lsm_ruc) then - !-- use RUC LSM's snow-cover fraction for ice - fsno0 = sncovr_ice(i) ! snow fraction on ice - else - asnow = 0.02*snowf(i) - argh = min(0.50, max(.025, 0.01*zorlf(i))) - hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) ) - fsno0 = asnow / (argh + asnow) * hrgh - endif - endif - - fsno1 = f_one - fsno0 ! snow-free fraction (land or ice), 1-sea - flnd0 = min(f_one, facsf(i)+facwf(i)) ! 1-land, 0-sea/ice - fsea0 = max(f_zero, f_one-flnd0)! ! 1-sea/ice, 0-land - fsno = fsno0 ! snow cover, >0 - land/ice - fsea = fsea0 * fsno1 ! 1-sea/ice, 0-land - flnd = flnd0 * fsno1 ! <=1-land,0-sea/ice - -!> - Calculate diffused sea surface albedo. - - if (tsknf(i) >= 271.5) then - asevd = 0.06 - asend = 0.06 - elseif (tsknf(i) < 271.1) then - asevd = 0.70 - asend = 0.65 - else - a1 = (tsknf(i) - 271.1)**2 - asevd = 0.7 - 4.0*a1 - asend = 0.65 - 3.6875*a1 - endif - -!> - Calculate diffused snow albedo, land area use input max snow -!! albedo. - - if (nint(slmsk(i)) == 2) then - ffw = f_one - fice(i) - if (ffw < f_one) then - dtgd = max(f_zero, min(5.0, (con_ttp-tisfc(i)) )) - b1 = 0.03 * dtgd + if (icy(i)) then + !-- Computation of ice albedo + asnow = 0.02*snowf(i) + argh = min(0.50, max(.025, 0.01*zorlf(i))) + hrgh = min(f_one,max(0.20,1.0577-1.1538e-3*hprif(i))) + fsno0 = asnow / (argh + asnow) * hrgh + ! diffused + if (tsknf(i) < 271.1) then + asevd_ice = 0.70 + asend_ice = 0.65 else - b1 = f_zero + a1 = (tsknf(i) - 271.1)**2 + asevd_ice = 0.7 - 4.0*a1 + asend_ice = 0.65 - 3.6875*a1 endif + ! direct + asevb_ice = asevd_ice + asenb_ice = asend_ice + + if (fsno0 > f_zero) then + ! Snow on ice + dtgd = max(f_zero, min(5.0, (con_ttp-tisfc(i)) )) + b1 = 0.03 * dtgd + asnvd = (asevd_ice + b1) ! diffused snow albedo + asnnd = (asend_ice + b1) + if (coszf(i) > 0.0001 .and. coszf(i) < 0.5) then ! direct snow albedo + csnow = 0.5 * (3.0 / (f_one+4.0*coszf(i)) - f_one) + asnvb = min( 0.98, asnvd+(f_one-asnvd)*csnow ) + asnnb = min( 0.98, asnnd+(f_one-asnnd)*csnow ) + else + asnvb = asnvd + asnnb = asnnd + endif - b3 = 0.06 * ffw - asnvd = (0.70 + b1) * fice(i) + b3 - asnnd = (0.60 + b1) * fice(i) + b3 - asevd = 0.70 * fice(i) + b3 - asend = 0.60 * fice(i) + b3 - else - asnvd = snoalb(i) - asnnd = snoalb(i) - endif - -!> - Calculate direct snow albedo. - - if (nint(slmsk(i)) == 2) then - if (coszf(i) < 0.5) then - csnow = 0.5 * (3.0 / (f_one+4.0*coszf(i)) - f_one) - asnvb = min( 0.98, asnvd+(f_one-asnvd)*csnow ) - asnnb = min( 0.98, asnnd+(f_one-asnnd)*csnow ) - else - asnvb = asnvd - asnnb = asnnd - endif - else - asnvb = snoalb(i) - asnnb = snoalb(i) - endif - -!> - Calculate direct sea surface albedo, use fanglin's zenith angle -!! treatment. - - if (coszf(i) > 0.0001) then - -! rfcs = 1.89 - 3.34*coszf(i) + 4.13*coszf(i)*coszf(i) & -! & - 2.02*coszf(i)*coszf(i)*coszf(i) - rfcs = 1.775/(1.0+1.55*coszf(i)) + ! composite ice and snow albedos + asevd_ice = asevd_ice * (1. - fsno0) + asnvd * fsno0 + asend_ice = asend_ice * (1. - fsno0) + asnnd * fsno0 + asevb_ice = asevb_ice * (1. - fsno0) + asnvb * fsno0 + asenb_ice = asenb_ice * (1. - fsno0) + asnnb * fsno0 + endif ! snow + else + ! icy = false, fill in values + asevd_ice = 0.70 + asend_ice = 0.65 + asevb_ice = 0.70 + asenb_ice = 0.65 + endif ! end icy - if (tsknf(i) >= con_t0c) then - !- sea - asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) & - & + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) & - & * (coszf(i)-f_one)) - asenb = asevb + if (fracl(i) > f_zero) then +!> - Calculate snow cover input directly for land model, no +!! conversion needed. + + fsno0 = sncovr(i) ! snow fraction on land + + fsno1 = f_one - fsno0 + flnd0 = min(f_one, facsf(i)+facwf(i)) + flnd = flnd0 * fsno1 ! snow-free fraction + fsno = f_one - flnd ! snow-covered fraction + + !> - use Fanglin's zenith angle treatment. + if (coszf(i) > 0.0001) then + rfcs = 1.775/(1.0+1.55*coszf(i)) else - !- ice - asevb = asevd - asenb = asend + !- no sun + rfcs = f_one endif - else - !- no sun - rfcs = f_one - asevb = asevd - asenb = asend - endif - - !- zenith dependence is applied only to direct beam albedo - ab1bm = min(0.99, alnsf(i)*rfcs) - ab2bm = min(0.99, alvsf(i)*rfcs) - sfcalb(i,1) = ab1bm *flnd + asenb*fsea + asnnb*fsno - sfcalb(i,2) = alnwf(i)*flnd + asend*fsea + asnnd*fsno - sfcalb(i,3) = ab2bm *flnd + asevb*fsea + asnvb*fsno - sfcalb(i,4) = alvwf(i)*flnd + asevd*fsea + asnvd*fsno + !- zenith dependence is applied only to direct beam albedo + ab1bm = min(0.99, alnsf(i)*rfcs) + ab2bm = min(0.99, alvsf(i)*rfcs) + + alndnb = ab1bm *flnd + snoalb(i) * fsno + alndnd = alnwf(i)*flnd + snoalb(i) * fsno + alndvb = ab2bm *flnd + snoalb(i) * fsno + alndvd = alvwf(i)*flnd + snoalb(i) * fsno + else + !-- fill in values of land albedo + alndnb = 0. + alndnd = 0. + alndvb = 0. + alndvd = 0. + endif ! end land + + !-- Composite mean surface albedo from land, open water and + !-- ice fractions + sfcalb(i,1) = min(0.99,max(0.01,alndnb))*fracl(i) & ! direct beam NIR + & + asenb_wat*fraco(i) + asenb_ice*fraci(i) + sfcalb(i,2) = min(0.99,max(0.01,alndnd))*fracl(i) & ! diffuse NIR + & + asend_wat*fraco(i) + asend_ice*fraci(i) + sfcalb(i,3) = min(0.99,max(0.01,alndvb))*fracl & ! direct beam visible + & + asevb_wat*fraco(i) + asevb_ice*fraci(i) + sfcalb(i,4) = min(0.99,max(0.01,alndvd))*fracl & ! diffuse visible + & + asevd_wat*fraco(i) + asevd_ice*fraci(i) enddo ! end_do_i_loop