From 4cbea925e49b885a8e7b63aa300cc13402f28a52 Mon Sep 17 00:00:00 2001 From: wx20hw Date: Sun, 30 Jan 2022 05:10:14 +0000 Subject: [PATCH 01/21] set up option for thermal roughness --- physics/module_sf_noahmp_glacier.f90 | 24 +++++++--- physics/module_sf_noahmplsm.f90 | 66 ++++++++++++++++++++++++++-- physics/noahmp_tables.f90 | 2 +- physics/sfc_noahmp_drv.F90 | 10 +++-- physics/sfc_noahmp_drv.meta | 7 +++ 5 files changed, 96 insertions(+), 13 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90 index 4c3a53c88..26dd810c7 100644 --- a/physics/module_sf_noahmp_glacier.f90 +++ b/physics/module_sf_noahmp_glacier.f90 @@ -62,6 +62,7 @@ module noahmp_glacier_globals INTEGER :: OPT_GLA != 1 !(suggested 1) INTEGER :: OPT_SFC != 1 !(suggested 1) + INTEGER :: OPT_TRS != 1 !(suggested 2) ! adjustable parameters for snow processes @@ -1129,8 +1130,10 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso real (kind=kind_phys) :: b !< temporary calculation real (kind=kind_phys) :: t, tdc !< kelvin to degree celsius with limit -50 to +50 real (kind=kind_phys), dimension( 1:nsoil) :: sice !< soil ice + real (kind=kind_phys) :: czil !< calculate roughness length of heat tdc(t) = min( 50., max(-50.,(t-tfrz)) ) + czil=0.5 ! ----------------------------------------------------------------- ! initialization variables that do not depend on stability iteration @@ -1155,10 +1158,18 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso fv = ur*vkc/log(zlvli/z0m) reyni = fv*z0m/(1.5e-05) !introduction of fv dependent z0h for the iter - if (reyni .gt. 2.0) then - z0h = z0m/exp(2.46*(reyni)**0.25 - log(7.4)) !Brutsaert 1982 - else - z0h = z0m/exp(-log(0.397)) !Brusaert 1982, table 4 + if (opt_trs == 1) then + z0h = z0m + elseif (opt_trs == 2) then + z0h = z0m*exp(-czil*0.4*258.2*sqrt(fv*z0m)) + elseif (opt_trs == 3) then + z0h = z0m*0.0001 + elseif (opt_trs == 4) then + if (reyni .gt. 2.0) then + z0h = z0m/exp(2.46*(reyni)**0.25 - log(7.4)) !Brutsaert 1982 + else + z0h = z0m/exp(-log(0.397)) !Brusaert 1982, table 4 + endif endif z0h_total = z0h @@ -3328,7 +3339,8 @@ end subroutine error_glacier ! ================================================================================================== !>\ingroup NoahMP_LSM - subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iopt_gla, iopt_sfc) + subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iopt_gla,& + iopt_sfc, iopt_trs) implicit none @@ -3339,6 +3351,7 @@ subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iop !! 1 -> semi-implicit; 2 -> full implicit (original noah) integer, intent(in) :: iopt_gla !< glacier option (1->phase change; 2->simple) integer, intent(in) :: iopt_sfc !< sfc scheme option + integer, intent(in) :: iopt_trs !< thermal roughness option ! ------------------------------------------------------------------------------------------------- @@ -3348,6 +3361,7 @@ subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iop opt_stc = iopt_stc opt_gla = iopt_gla opt_sfc = iopt_sfc + opt_trs = iopt_trs end subroutine noahmp_options_glacier diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 944446085..b602a683e 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -159,6 +159,11 @@ module module_sf_noahmplsm ! **0 -> no crop model, will run default dynamic vegetation ! 1 -> liu, et al. 2016 + integer :: opt_trs !< options for thermal roughness scheme + ! **1 -> z0h=z0 + ! 2 -> czil + ! 3 -> ec style + ! 4 -> kb inversed !------------------------------------------------------------------------------------------! ! physical constants: ! !------------------------------------------------------------------------------------------! @@ -2241,6 +2246,21 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in q1 = fveg * (eah*0.622/(sfcprs - 0.378*eah)) + (1.0 - fveg)*qsfc q2e = fveg * q2v + (1.0 - fveg) * q2b + if (opt_trs == 1) then + z0wrf = fveg * z0m + (1.0 - fveg) * z0mg + z0hwrf = z0wrf + elseif (opt_trs == 2) then + z0wrf = fveg * z0m + (1.0 - fveg) * z0mg + z0hwrf = fveg * z0m*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0m)) & + +(1.0 - fveg) * z0mg*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0mg)) + elseif (opt_trs == 3) then + z0wrf = fveg * z0m + (1.0 - fveg) * z0mg + if (vegtyp.le.5) then + z0hwrf = fveg * z0m + (1.0 - fveg) * z0mg*0.1 + else + z0hwrf = fveg * z0m*0.01 + (1.0 - fveg) * z0mg*0.1 + endif + elseif (opt_trs == 4) then coeffa = (csigmaf0 - csigmaf1)/(1.0 - exp(-1.0*aone)) coeffb = csigmaf0 - coeffa csigmafveg = coeffa * exp(-1.0*aone*fveg) + coeffb @@ -2259,6 +2279,9 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in kbsigmafveg = csigmafveg/log((zlvl-ezpd)/z0wrf) - log((zlvl-ezpd)/z0wrf) z0hwrf = z0wrf/exp(kbsigmafveg) +! place holder doe other roughness scheme +! elseif (opt_trs == x) then + endif else taux = tauxb @@ -2283,7 +2306,19 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in chv = chb z0wrf = z0mg + if (opt_trs == 1) then + z0hwrf = z0wrf + elseif (opt_trs == 2) then + z0hwrf = z0wrf*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0wrf)) + elseif (opt_trs == 3) then + if (vegtyp.le.5) then + z0hwrf = z0wrf + else + z0hwrf = z0wrf*0.01 + endif + elseif (opt_trs == 4) then z0hwrf =z0wrf/exp( csigmaf0/log((zlvl-ezpd)/z0wrf) - log((zlvl-ezpd)/z0wrf) ) + endif end if @@ -3965,11 +4000,22 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & cir = (2.-emv*(1.-emg))*emv*sb ! --------------------------------------------------------------------------------------------- + if (opt_trs == 1) then + z0h = z0m + elseif (opt_trs == 2) then + z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) + elseif (opt_trs == 3) then + if (vegtyp.le.5) then + z0h = z0m + else + z0h = z0m*0.01 + endif + elseif (opt_trs == 4) then sigmaa = 1.0 - (0.5/(0.5+vaie))*exp(-vaie**2/8.0) kbsigmaf1 = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(dlf*ur/log((zlvl-zpd)/z0m)) z0h = z0m/exp(kbsigmaf1) csigmaf1 = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m)+kbsigmaf1) ! for output for interpolation - + endif ! -- tem1 = (z0m - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) @@ -4582,7 +4628,19 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & csigmaf0 = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m) + kbsigmaf0) - z0h = max(z0m/exp(kbsigmaf0),1.0e-6) + if (opt_trs == 1) then + z0h = z0m + elseif (opt_trs == 2) then + z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) + elseif (opt_trs == 3) then + if (vegtyp.le.5) then + z0h = z0m + else + z0h = z0m*0.01 + endif + elseif (opt_trs == 4) then + z0h = max(z0m/exp(kbsigmaf0),1.0e-6) + endif ! ! for sfcdiff3; maybe should move to inside the option ! @@ -9782,7 +9840,7 @@ end subroutine psn_crop !>\ingroup NoahMP_LSM subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, & - iopt_rsf , iopt_soil, iopt_pedo, iopt_crop ) + iopt_rsf , iopt_soil, iopt_pedo, iopt_crop ,iopt_trs ) implicit none @@ -9804,6 +9862,7 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc integer, intent(in) :: iopt_soil !soil parameters set-up option integer, intent(in) :: iopt_pedo !pedo-transfer function (1->saxton and rawls) integer, intent(in) :: iopt_crop !crop model option (0->none; 1->liu et al.) + integer, intent(in) :: iopt_trs !thermal roughness scheme option (1->z0h=z0; 2->rb reversed) ! ------------------------------------------------------------------------------------------------- @@ -9824,6 +9883,7 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc opt_soil = iopt_soil opt_pedo = iopt_pedo opt_crop = iopt_crop + opt_trs = iopt_trs end subroutine noahmp_options diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90 index 9cb25b3f3..5f6246a0f 100644 --- a/physics/noahmp_tables.f90 +++ b/physics/noahmp_tables.f90 @@ -735,7 +735,7 @@ module noahmp_tables real :: refkdt_table = 3.0 !< parameter in the surface runoff parameterization real :: frzk_table =0.15 !< frozen ground parameter real :: zbot_table = -8.0 !< depth [m] of lower boundary soil temperature - real :: czil_table = 0.1 !< parameter used in the calculation of the roughness length for heat + real :: czil_table = 0.5 !< parameter used in the calculation of the roughness length for heat ! mptable.tbl radiation parameters diff --git a/physics/sfc_noahmp_drv.F90 b/physics/sfc_noahmp_drv.F90 index 1fd9773ff..397a09674 100644 --- a/physics/sfc_noahmp_drv.F90 +++ b/physics/sfc_noahmp_drv.F90 @@ -111,7 +111,7 @@ subroutine noahmpdrv_run & shdmin, shdmax, snoalb, sfalb, flag_iter,con_g, & idveg, iopt_crs, iopt_btr, iopt_run, iopt_sfc, iopt_frz, & iopt_inf, iopt_rad, iopt_alb, iopt_snf, iopt_tbot, & - iopt_stc, xlatin, xcoszin, iyrlen, julian, garea, & + iopt_stc, iopt_trs,xlatin, xcoszin, iyrlen, julian, garea, & rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, & con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1, & con_fvirt, con_rd, con_hfus, thsfc_loc, & @@ -213,6 +213,7 @@ subroutine noahmpdrv_run & integer , intent(in) :: iopt_snf ! option for partitioning precipitation into rainfall & snowfall integer , intent(in) :: iopt_tbot ! option for lower boundary condition of soil temperature integer , intent(in) :: iopt_stc ! option for snow/soil temperature time scheme (only layer 1) + integer , intent(in) :: iopt_trs ! option for thermal roughness scheme real(kind=kind_phys), dimension(:) , intent(in) :: xlatin ! latitude real(kind=kind_phys), dimension(:) , intent(in) :: xcoszin ! cosine of zenith angle integer , intent(in) :: iyrlen ! year length [days] @@ -700,8 +701,8 @@ subroutine noahmpdrv_run & call noahmp_options(idveg ,iopt_crs, iopt_btr , iopt_run, iopt_sfc, & iopt_frz, iopt_inf , iopt_rad, iopt_alb, & - iopt_snf, iopt_tbot, iopt_stc, & - iopt_rsf, iopt_soil, iopt_pedo, iopt_crop ) + iopt_snf, iopt_tbot, iopt_stc, iopt_rsf, & + iopt_soil,iopt_pedo, iopt_crop,iopt_trs ) if ( vegetation_category == isice_table ) then @@ -714,7 +715,8 @@ subroutine noahmpdrv_run & ice_flag = -1 temperature_soil_bot = min(temperature_soil_bot,263.15) - call noahmp_options_glacier(iopt_alb, iopt_snf, iopt_tbot, iopt_stc, iopt_gla, iopt_sfc ) + call noahmp_options_glacier(iopt_alb, iopt_snf, iopt_tbot, iopt_stc, iopt_gla, & + iopt_sfc ,iopt_trs) call noahmp_glacier ( & i_location ,1 ,cosine_zenith ,nsnow , & diff --git a/physics/sfc_noahmp_drv.meta b/physics/sfc_noahmp_drv.meta index e37036c32..712a457a6 100644 --- a/physics/sfc_noahmp_drv.meta +++ b/physics/sfc_noahmp_drv.meta @@ -424,6 +424,13 @@ dimensions = () type = integer intent = in +[iopt_trs] + standard_name = control_for_land_surface_scheme_surface_thermal_roughness + long_name = choice for surface thermal roughness option (see noahmp module for definition) + units = index + dimensions = () + type = integer + intent = in [xlatin] standard_name = latitude long_name = latitude From 22de66b8306b687585c5521367ab1a706e04209d Mon Sep 17 00:00:00 2001 From: wx20hw Date: Mon, 31 Jan 2022 17:08:50 +0000 Subject: [PATCH 02/21] change czil --- physics/module_sf_noahmp_glacier.f90 | 4 ++-- physics/noahmp_tables.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90 index 26dd810c7..c4c03aaf8 100644 --- a/physics/module_sf_noahmp_glacier.f90 +++ b/physics/module_sf_noahmp_glacier.f90 @@ -1133,7 +1133,7 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso real (kind=kind_phys) :: czil !< calculate roughness length of heat tdc(t) = min( 50., max(-50.,(t-tfrz)) ) - czil=0.5 + czil=0.1 ! ----------------------------------------------------------------- ! initialization variables that do not depend on stability iteration @@ -1163,7 +1163,7 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso elseif (opt_trs == 2) then z0h = z0m*exp(-czil*0.4*258.2*sqrt(fv*z0m)) elseif (opt_trs == 3) then - z0h = z0m*0.0001 + z0h = z0m*0.1 elseif (opt_trs == 4) then if (reyni .gt. 2.0) then z0h = z0m/exp(2.46*(reyni)**0.25 - log(7.4)) !Brutsaert 1982 diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90 index 5f6246a0f..9cb25b3f3 100644 --- a/physics/noahmp_tables.f90 +++ b/physics/noahmp_tables.f90 @@ -735,7 +735,7 @@ module noahmp_tables real :: refkdt_table = 3.0 !< parameter in the surface runoff parameterization real :: frzk_table =0.15 !< frozen ground parameter real :: zbot_table = -8.0 !< depth [m] of lower boundary soil temperature - real :: czil_table = 0.5 !< parameter used in the calculation of the roughness length for heat + real :: czil_table = 0.1 !< parameter used in the calculation of the roughness length for heat ! mptable.tbl radiation parameters From 10fa17e895ecd21db0d24d1cef7b12523cabce40 Mon Sep 17 00:00:00 2001 From: wx20hw Date: Wed, 16 Feb 2022 20:03:11 +0000 Subject: [PATCH 03/21] canopy height dependant czil --- physics/module_sf_noahmplsm.f90 | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index b602a683e..0fc4e8948 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -1895,6 +1895,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) :: csigmaf0 real (kind=kind_phys) :: csigmaf1 real (kind=kind_phys) :: csigmafveg + real (kind=kind_phys) :: czil1 real (kind=kind_phys) :: cdmnv real (kind=kind_phys) :: ezpdv @@ -2251,8 +2252,11 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in z0hwrf = z0wrf elseif (opt_trs == 2) then z0wrf = fveg * z0m + (1.0 - fveg) * z0mg - z0hwrf = fveg * z0m*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0m)) & - +(1.0 - fveg) * z0mg*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0mg)) +! z0hwrf = fveg * z0m*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0m)) & +! +(1.0 - fveg) * z0mg*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0mg)) + czil1=10.0 ** (- (0.40/0.07) * parameters%hvt) + z0hwrf = fveg * z0m*exp(-czil1*0.4*258.2*sqrt(ustarx*z0m)) & + +(1.0 - fveg) * z0mg*exp(-czil1*0.4*258.2*sqrt(ustarx*z0mg)) elseif (opt_trs == 3) then z0wrf = fveg * z0m + (1.0 - fveg) * z0mg if (vegtyp.le.5) then @@ -2309,7 +2313,9 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in if (opt_trs == 1) then z0hwrf = z0wrf elseif (opt_trs == 2) then - z0hwrf = z0wrf*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0wrf)) +! z0hwrf = z0wrf*exp(-parameters%czil*0.4*258.2*sqrt(ustarx*z0wrf)) + czil1=10.0 ** (- (0.40/0.07) * parameters%hvt) + z0hwrf = z0wrf*exp(-czil1*0.4*258.2*sqrt(ustarx*z0wrf)) elseif (opt_trs == 3) then if (vegtyp.le.5) then z0hwrf = z0wrf @@ -3866,7 +3872,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: laisune !sunlit leaf area index, one-sided (m2/m2),effective real (kind=kind_phys) :: laishae !shaded leaf area index, one-sided (m2/m2),effective - real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx + real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx,czil1 real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 integer :: k !index @@ -4003,7 +4009,9 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if (opt_trs == 1) then z0h = z0m elseif (opt_trs == 2) then - z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) +! z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) + czil1= 10.0 ** (- (0.40/0.07) * hcan) + z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) elseif (opt_trs == 3) then if (vegtyp.le.5) then z0h = z0m @@ -4581,7 +4589,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys) :: fh2 !monin-obukhov heat adjustment at 2m real (kind=kind_phys) :: ch2 !surface exchange at 2m - real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx + real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx,czil1 real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 integer :: iter !iteration index @@ -4631,7 +4639,9 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if (opt_trs == 1) then z0h = z0m elseif (opt_trs == 2) then - z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) +! z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) + czil1= 10.0 ** (- (0.40/0.07) * parameters%hvt) + z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) elseif (opt_trs == 3) then if (vegtyp.le.5) then z0h = z0m From b90d4e2c7b2770295406d5344806fb52c1c1d41b Mon Sep 17 00:00:00 2001 From: helin wei Date: Mon, 7 Mar 2022 16:08:57 +0000 Subject: [PATCH 04/21] add canopy heat storage and gvf impact on thermal conductivity --- physics/module_sf_noahmp_glacier.f90 | 3 +- physics/module_sf_noahmplsm.f90 | 88 ++++++++++++++++++---------- physics/sfc_diag_post.F90 | 12 +++- physics/sfc_diag_post.meta | 16 +++++ physics/sfc_noahmp_drv.F90 | 4 +- 5 files changed, 88 insertions(+), 35 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90 index c4c03aaf8..1ea4a45b8 100644 --- a/physics/module_sf_noahmp_glacier.f90 +++ b/physics/module_sf_noahmp_glacier.f90 @@ -1152,7 +1152,8 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ! the following only applies to opt_sfc =3, opt_sfc = 1 still done its old way snwd = snowh*1000.0 - zlvli = zlvl - zpd +! zlvli = zlvl - zpd + zlvli = zlvl ! fv = ustarx ! the input maybe too high for glacial fv = ur*vkc/log(zlvli/z0m) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 0fc4e8948..0913531f8 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -678,18 +678,21 @@ subroutine noahmp_sflx (parameters, & real (kind=kind_phys) :: latheag !< latent heat vap./sublimation (j/kg) logical :: frozen_ground !< used to define latent heat pathway logical :: frozen_canopy !< used to define latent heat pathway - LOGICAL :: dveg_active !< flag to run dynamic vegetation - LOGICAL :: crop_active !< flag to run crop model + logical :: dveg_active !< flag to run dynamic vegetation + logical :: crop_active !< flag to run crop model +! add canopy heat storage (C.He added based on GY Niu's communication) + real :: canhs ! canopy heat storage change w/m2 ! intent (out) variables need to be assigned a value. these normally get assigned values ! only if dveg == 2. nee = 0.0 npp = 0.0 gpp = 0.0 - pahv = 0. - pahg = 0. - pahb = 0. - pah = 0. + pahv = 0. + pahg = 0. + pahb = 0. + pah = 0. + canhs = 0. ! -------------------------------------------------------------------------------------------------- ! re-process atmospheric forcing @@ -774,7 +777,7 @@ subroutine noahmp_sflx (parameters, & co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in elai ,esai ,fwet ,foln , & !in - fveg ,pahv ,pahg ,pahb , & !in + fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in z0wrf ,z0hwrf , & !out @@ -797,7 +800,7 @@ subroutine noahmp_sflx (parameters, & t2mv ,t2mb ,fsrv , & fsrg ,rssun ,rssha ,albd ,albi ,albsnd,albsni, bgap ,wgap, tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb , & !out - emissi ,pah , & + emissi ,pah ,canhs, & shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out qsfc = q1 ! @@ -868,9 +871,9 @@ subroutine noahmp_sflx (parameters, & nsnow ,ist ,errwat ,iloc , jloc ,fveg , & sav ,sag ,fsrv ,fsrg ,zwt ,pah , & #ifdef CCPP - pahv ,pahg ,pahb ,errmsg, errflg) !in ( except errwat [out] and errmsg, errflg [inout] ) + pahv ,pahg ,pahb ,canhs,errmsg, errflg) !in ( except errwat [out] and errmsg, errflg [inout] ) #else - pahv ,pahg ,pahb ) !in ( except errwat, which is out ) + pahv ,pahg ,pahb, canhs ) !in ( except errwat, which is out ) #endif #ifdef CCPP @@ -1405,9 +1408,9 @@ subroutine error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & nsnow ,ist ,errwat, iloc ,jloc ,fveg , & sav ,sag ,fsrv ,fsrg ,zwt ,pah , & #ifdef CCPP - pahv ,pahg ,pahb ,errmsg, errflg) + pahv ,pahg ,pahb ,canhs,errmsg, errflg) #else - pahv ,pahg ,pahb ) + pahv ,pahg ,pahb ,canhs) #endif ! -------------------------------------------------------------------------------------------------- ! check surface energy balance and water balance @@ -1456,6 +1459,7 @@ subroutine error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & real (kind=kind_phys), intent(in) :: pahv !precipitation advected heat - total (w/m2) real (kind=kind_phys), intent(in) :: pahg !precipitation advected heat - total (w/m2) real (kind=kind_phys), intent(in) :: pahb !precipitation advected heat - total (w/m2) + real (kind=kind_phys), intent(in) :: canhs !canopy heat storage change (w/m2) C.He added based on GY Niu's communication #ifdef CCPP character(len=*) , intent(inout) :: errmsg @@ -1501,7 +1505,7 @@ subroutine error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & #endif end if - erreng = sav+sag-(fira+fsh+fcev+fgev+fctr+ssoil) +pah + erreng = sav+sag-(fira+fsh+fcev+fgev+fctr+ssoil+canhs) +pah ! erreng = fveg*sav+sag-(fira+fsh+fcev+fgev+fctr+ssoil) if(abs(erreng) > 0.01) then write(message,*) 'erreng =',erreng,' at i,j: ',iloc,jloc @@ -1551,6 +1555,12 @@ subroutine error (parameters,swdown ,fsa ,fsr ,fira ,fsh ,fcev , & errmsg = trim(errmsg)//NEW_LINE('A')//trim(message) #else call wrf_message(trim(message)) +#endif + write(message,'(a17,f10.4)') "canopy heat storage: ",canhs +#ifdef CCPP + errmsg = trim(errmsg)//NEW_LINE('A')//trim(message) +#else + call wrf_message(trim(message)) #endif write(message,'(a17,4f10.4)') "precip advected: ",pah,pahv,pahg,pahb #ifdef CCPP @@ -1605,7 +1615,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in elai ,esai ,fwet ,foln , & !in - fveg ,pahv ,pahg ,pahb , & !in + fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in z0wrf ,z0hwrf , & !out @@ -1627,7 +1637,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in qc ,qsfc ,psfc , & !in t2mv ,t2mb ,fsrv , & fsrg ,rssun ,rssha ,albd ,albi,albsnd ,albsni,bgap ,wgap,tgv,tgb,& - q1 ,q2v ,q2b ,q2e ,chv ,chb, emissi,pah ,& + q1 ,q2v ,q2b ,q2e ,chv ,chb, emissi,pah,canhs,& shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out !jref:end @@ -1701,6 +1711,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: esai !lai adjusted for burying by snow real (kind=kind_phys) , intent(in) :: fwet !fraction of canopy that is wet [-] real (kind=kind_phys) , intent(in) :: fveg !greeness vegetation fraction (-) + real (kind=kind_phys) , intent(in) :: shdfac !< green vegetation fraction [0.0-1.0] real (kind=kind_phys) , intent(in) :: lat !latitude (radians) real (kind=kind_phys) , intent(in) :: canliq !canopy-intercepted liquid water (mm) real (kind=kind_phys) , intent(in) :: canice !canopy-intercepted ice mass (mm) @@ -1774,6 +1785,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(out) :: t2mb !2-m air temperature over bare ground part [k] real (kind=kind_phys) , intent(out) :: bgap real (kind=kind_phys) , intent(out) :: wgap + real (kind=kind_phys) , intent(out) :: canhs !canopy heat storage change (w/m2) real (kind=kind_phys), dimension(1:2) , intent(out) :: albd !albedo (direct) real (kind=kind_phys), dimension(1:2) , intent(out) :: albi !albedo (diffuse) real (kind=kind_phys), dimension(1:2) , intent(out) :: albsnd !snow albedo (direct) @@ -2032,7 +2044,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in call thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , & !in + lat ,z0m ,zlvl ,vegtyp , fveg, & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out @@ -2157,7 +2169,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag , & !in fwet ,laisun ,laisha ,cwp ,dzsnso , & !in - zlvl ,zpd ,z0m ,fveg , & !in + zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in rsurf ,latheav ,latheag ,parsun ,parsha ,igs , & !in @@ -2172,7 +2184,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in #endif tauxv ,tauyv ,irg ,irc ,shg , & !out shc ,evg ,evc ,tr ,ghv , & !out - t2mv ,psnsun ,psnsha ,csigmaf1, & !out + t2mv ,psnsun ,psnsha ,csigmaf1,canhs, & !out !jref:start qc ,qsfc ,psfc , & !in q2v ,chv2, chleaf, chuc) !inout @@ -2196,7 +2208,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in dzsnso ,zlvl ,zpdg ,z0mg ,fsno, & !in emg ,stc ,df ,rsurf ,latheag , & !in gammag ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in - thsfc_loc, prslkix,prsik1x,prslk1x,fveg,garea1, & !in + thsfc_loc, prslkix,prsik1x,prslk1x,fveg,shdfac,garea1, & !in #ifdef CCPP tgb ,cmb ,chb, ustarx,errmsg ,errflg , & !inout #else @@ -2415,7 +2427,7 @@ end subroutine energy subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , & !in + lat ,z0m ,zlvl ,vegtyp , fveg,& !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out ! ------------------------------------------------------------------------------------------------- @@ -2441,6 +2453,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), intent(in) :: z0m !roughness length (m) real (kind=kind_phys), intent(in) :: zlvl !reference height (m) integer , intent(in) :: vegtyp !vegtyp type + real (kind=kind_phys), intent(in) :: fveg !green vegetation fraction [0.0-1.0] ! outputs real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: df !thermal conductivity [w/m/k] @@ -2456,6 +2469,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), dimension(-nsnow+1: 0) :: cvsno !volumetric specific heat (j/m3/k) real (kind=kind_phys), dimension(-nsnow+1: 0) :: tksno !snow thermal conductivity (j/m3/k) real (kind=kind_phys), dimension( 1:nsoil) :: sice !soil ice content + real (kind=kind_phys), parameter :: sbeta = -2.0 ! -------------------------------------------------------------------------------------------------- ! compute snow thermal conductivity and heat capacity @@ -2488,6 +2502,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) + df = df * exp (sbeta * fveg) ! compute lake thermal properties ! (no consideration of turbulent mixing for this version) @@ -3634,7 +3649,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag, & !in fwet ,laisun ,laisha ,cwp ,dzsnso , & !in - zlvl ,zpd ,z0m ,fveg, & !in + zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in rsurf ,latheav ,latheag ,parsun ,parsha ,igs , & !in @@ -3649,7 +3664,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & #endif tauxv ,tauyv ,irg ,irc ,shg , & !out shc ,evg ,evc ,tr ,gh , & !out - t2mv ,psnsun ,psnsha ,csigmaf1, & !out + t2mv ,psnsun ,psnsha ,csigmaf1,canhs, & !out qc ,qsfc ,psfc , & !in q2v ,cah2 ,chleaf ,chuc ) !inout @@ -3658,7 +3673,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & ! ground (tg) temperatures that balance the surface energy budgets ! vegetated: -! -sav + irc[tv] + shc[tv] + evc[tv] + tr[tv] = 0 +! -sav + irc[tv] + shc[tv] + evc[tv] + tr[tv] + canhs(tv) = 0 ! -sag + irg[tg] + shg[tg] + evg[tg] + gh[tg] = 0 ! -------------------------------------------------------------------------------------------------- implicit none @@ -3673,6 +3688,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & integer, intent(in) :: isnow !actual no. of snow layers integer, intent(in) :: vegtyp !vegetation physiology type real (kind=kind_phys), intent(in) :: fveg !greeness vegetation fraction (-) + real (kind=kind_phys), intent(in) :: shdfac !greeness vegetation fraction (-) real (kind=kind_phys), intent(in) :: sav !solar rad absorbed by veg (w/m2) real (kind=kind_phys), intent(in) :: sag !solar rad absorbed by ground (w/m2) real (kind=kind_phys), intent(in) :: lwdn !atmospheric longwave radiation (w/m2) @@ -3753,7 +3769,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & #endif ! output -! -fsa + fira + fsh + (fcev + fctr + fgev) + fcst + ssoil = 0 +! -fsa + fira + fsh + (fcev + fctr + fgev) + fcst + ssoil + canhs = 0 real (kind=kind_phys), intent(out) :: tauxv !wind stress: e-w (n/m2) real (kind=kind_phys), intent(out) :: tauyv !wind stress: n-s (n/m2) real (kind=kind_phys), intent(out) :: irc !net longwave radiation (w/m2) [+= to atm] @@ -3770,6 +3786,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(out) :: csigmaf1 real (kind=kind_phys), intent(out) :: chleaf !leaf exchange coefficient real (kind=kind_phys), intent(out) :: chuc !under canopy exchange coefficient + real (kind=kind_phys), intent(out) :: canhs !canopy heat storage change (w/m2) real (kind=kind_phys), intent(out) :: q2v real (kind=kind_phys) :: cah !sensible heat conductance, canopy air to zlvl air (m/s) @@ -3864,8 +3881,9 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: ch2v !exchange coefficient for 2m over vegetation. real (kind=kind_phys) :: cq2v !exchange coefficient for 2m over vegetation. real (kind=kind_phys) :: eah2 !2m vapor pressure over canopy - real (kind=kind_phys) :: qfx !moisture flux + real (kind=kind_phys) :: qfx !moisture flux real (kind=kind_phys) :: e1 + real (kind=kind_phys) :: hcv !canopy heat capacity j/m2/k, C.He added real (kind=kind_phys) :: vaie !total leaf area index + stem area index,effective @@ -3929,7 +3947,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & ! for sfcdiff3 snwd = snowh*1000.0 - zlvlv = zlvl - zpd +! zlvlv = zlvl - zpd + zlvlv = zlvl virtfacv = 1.0 + 0.61 * max(qair, 1.e-8) tv1v = sfctmp * virtfacv @@ -4027,7 +4046,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & ! -- tem1 = (z0m - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(fveg, 0.1_kind_phys) + tem2 = max(shdfac, 0.1_kind_phys) zvfun1= sqrt(tem1 * tem2) gdx=sqrt(garea1) if(opt_sfc == 1 .or. opt_sfc == 2) then @@ -4156,14 +4175,19 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & evc = min(canice*latheav/dt,evc) end if +! canopy heat capacity + hcv = 0.02*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice !j/m2/k + b = sav-irc-shc-evc-tr+pahv !additional w/m2 - a = fveg*(4.*cir*tv**3 + csh + (cev+ctr)*destv) !volumetric heat capacity +! a = fveg*(4.*cir*tv**3 + csh + (cev+ctr)*destv) !volumetric heat capacity + a = fveg*(4.*cir*tv**3 + csh + (cev+ctr)*destv) + hcv/dt !volumetric heat capacity dtv = b/a irc = irc + fveg*4.*cir*tv**3*dtv shc = shc + fveg*csh*dtv evc = evc + fveg*cev*destv*dtv tr = tr + fveg*ctr*destv*dtv + canhs = dtv*hcv/dt ! update vegetation surface temperature tv = tv + dtv @@ -4413,7 +4437,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & dzsnso ,zlvl ,zpd ,z0m ,fsno , & !in emg ,stc ,df ,rsurf ,lathea , & !in gamma ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in - thsfc_loc, prslkix,prsik1x,prslk1x,fveg,garea1, & !in + thsfc_loc, prslkix,prsik1x,prslk1x,fveg,shdfac,garea1, & !in #ifdef CCPP tgb ,cm ,ch,ustarx,errmsg ,errflg , & !inout #else @@ -4470,6 +4494,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys) , intent(in) :: prsik1x ! in exner function real (kind=kind_phys) , intent(in) :: prslk1x ! in exner function real (kind=kind_phys) , intent(in) :: fveg + real (kind=kind_phys) , intent(in) :: shdfac real (kind=kind_phys) , intent(in) :: garea1 !jref:start; in @@ -4655,7 +4680,8 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & ! for sfcdiff3; maybe should move to inside the option ! snwd = snowh*1000.0 - zlvlb = zlvl - zpd +! zlvlb = zlvl - zpd + zlvlb = zlvl virtfacb = 1.0 + 0.61 * max(qair, 1.e-8) tv1b = sfctmp * virtfacb @@ -4672,7 +4698,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & ! ----------------------------------------------------------------- tem1 = (z0m - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(fveg, 0.1_kind_phys) + tem2 = max(shdfac, 0.1_kind_phys) zvfun1= sqrt(tem1 * tem2) gdx=sqrt(garea1) diff --git a/physics/sfc_diag_post.F90 b/physics/sfc_diag_post.F90 index 6f14fe93d..36541b0fc 100644 --- a/physics/sfc_diag_post.F90 +++ b/physics/sfc_diag_post.F90 @@ -16,7 +16,7 @@ end subroutine sfc_diag_post_finalize !! #endif subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con_epsm1, pgr,& - t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax, & + t2mmp,q2mp, t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax, & wind10mmax, u10mmax, v10mmax, dpt2m, errmsg, errflg) use machine, only: kind_phys @@ -29,6 +29,7 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con logical , dimension(:), intent(in) :: dry real(kind=kind_phys), dimension(:), intent(in) :: pgr, u10m, v10m real(kind=kind_phys), dimension(:), intent(inout) :: t2m, q2m, tmpmin, tmpmax, spfhmin, spfhmax + real(kind=kind_phys), dimension(:), intent(inout) :: t2mmp, q2mp real(kind=kind_phys), dimension(:), intent(inout) :: wind10mmax, u10mmax, v10mmax, dpt2m character(len=*), intent(out) :: errmsg @@ -41,6 +42,15 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con errmsg = '' errflg = 0 +! if (lsm == lsm_noahmp) then +! do i=1,im +! if(dry(i)) then +! t2m(i) = t2mmp(i) +! q2m(i) = q2mp(i) +! endif +! enddo +! endif + if (lssav) then do i=1,im tmpmax(i) = max(tmpmax(i),t2m(i)) diff --git a/physics/sfc_diag_post.meta b/physics/sfc_diag_post.meta index 21d76a147..95e8d8428 100644 --- a/physics/sfc_diag_post.meta +++ b/physics/sfc_diag_post.meta @@ -74,6 +74,22 @@ type = real kind = kind_phys intent = in +[t2mmp] + standard_name = temperature_at_2m_from_noahmp + long_name = 2 meter temperature from noahmp + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[q2mp] + standard_name = specific_humidity_at_2m_from_noahmp + long_name = 2 meter specific humidity from noahmp + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [t2m] standard_name = air_temperature_at_2m long_name = 2 meter temperature diff --git a/physics/sfc_noahmp_drv.F90 b/physics/sfc_noahmp_drv.F90 index 397a09674..0ebcbd615 100644 --- a/physics/sfc_noahmp_drv.F90 +++ b/physics/sfc_noahmp_drv.F90 @@ -923,7 +923,7 @@ subroutine noahmpdrv_run & snowc (i) = snow_cover_fraction sncovr1 (i) = snow_cover_fraction -! qsurf (i) = spec_humidity_surface + qsurf (i) = spec_humidity_surface tsurf (i) = tskin(i) tvxy (i) = temperature_leaf @@ -998,7 +998,7 @@ subroutine noahmpdrv_run & cmm (i) = cmxy(i) * wind(i) snwdph (i) = snow_depth * 1000.0 ! convert from m to mm; wait after the stability call - qsurf (i) = q1(i) + evap(i)/(con_hvap*density*ch(i)*wind(i)) +! qsurf (i) = q1(i) + evap(i)/(con_hvap*density*ch(i)*wind(i)) ! ! --- change units for output From ae7ac42b3679be983aac2dbd5b9f959c4f6db86f Mon Sep 17 00:00:00 2001 From: barlage Date: Mon, 7 Mar 2022 14:37:38 -0700 Subject: [PATCH 05/21] add sfcdif3 as a separate subroutine --- physics/module_sf_noahmplsm.f90 | 549 ++++++++++---------------------- 1 file changed, 167 insertions(+), 382 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 0913531f8..0248a116b 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2184,7 +2184,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in #endif tauxv ,tauyv ,irg ,irc ,shg , & !out shc ,evg ,evc ,tr ,ghv , & !out - t2mv ,psnsun ,psnsha ,csigmaf1,canhs, & !out + t2mv ,psnsun ,psnsha ,canhs , & !out + csigmaf1, & !out !jref:start qc ,qsfc ,psfc , & !in q2v ,chv2, chleaf, chuc) !inout @@ -3664,7 +3665,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & #endif tauxv ,tauyv ,irg ,irc ,shg , & !out shc ,evg ,evc ,tr ,gh , & !out - t2mv ,psnsun ,psnsha ,csigmaf1,canhs, & !out + t2mv ,psnsun ,psnsha ,canhs , & !out + csigmaf1, & !out qc ,qsfc ,psfc , & !in q2v ,cah2 ,chleaf ,chuc ) !inout @@ -3673,7 +3675,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & ! ground (tg) temperatures that balance the surface energy budgets ! vegetated: -! -sav + irc[tv] + shc[tv] + evc[tv] + tr[tv] + canhs(tv) = 0 +! -sav + irc[tv] + shc[tv] + evc[tv] + tr[tv] + canhs[tv] = 0 ! -sag + irg[tg] + shg[tg] + evg[tg] + gh[tg] = 0 ! -------------------------------------------------------------------------------------------------- implicit none @@ -3688,7 +3690,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & integer, intent(in) :: isnow !actual no. of snow layers integer, intent(in) :: vegtyp !vegetation physiology type real (kind=kind_phys), intent(in) :: fveg !greeness vegetation fraction (-) - real (kind=kind_phys), intent(in) :: shdfac !greeness vegetation fraction (-) real (kind=kind_phys), intent(in) :: sav !solar rad absorbed by veg (w/m2) real (kind=kind_phys), intent(in) :: sag !solar rad absorbed by ground (w/m2) real (kind=kind_phys), intent(in) :: lwdn !atmospheric longwave radiation (w/m2) @@ -3703,12 +3704,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: dt !time step (s) real (kind=kind_phys), intent(in) :: fsno !snow fraction - logical , intent(in) :: thsfc_loc - real (kind=kind_phys) , intent(in) :: prslkix ! in exner function - real (kind=kind_phys) , intent(in) :: prsik1x ! in exner function - real (kind=kind_phys) , intent(in) :: prslk1x ! in exner function - real (kind=kind_phys) , intent(in) :: garea1 ! - real (kind=kind_phys), intent(in) :: snowh !actual snow depth [m] real (kind=kind_phys), intent(in) :: fwet !wetted fraction of canopy real (kind=kind_phys), intent(in) :: cwp !canopy wind parameter @@ -3761,7 +3756,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(inout) :: tg !ground temperature (k) real (kind=kind_phys), intent(inout) :: cm !momentum drag coefficient real (kind=kind_phys), intent(inout) :: ch !sensible heat exchange coefficient - real (kind=kind_phys), intent(inout) :: ustarx !< friction velocity #ifdef CCPP character(len=*), intent(inout) :: errmsg @@ -3783,11 +3777,9 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(out) :: t2mv !2 m height air temperature (k) real (kind=kind_phys), intent(out) :: psnsun !sunlit leaf photosynthesis (umolco2/m2/s) real (kind=kind_phys), intent(out) :: psnsha !shaded leaf photosynthesis (umolco2/m2/s) - real (kind=kind_phys), intent(out) :: csigmaf1 real (kind=kind_phys), intent(out) :: chleaf !leaf exchange coefficient real (kind=kind_phys), intent(out) :: chuc !under canopy exchange coefficient real (kind=kind_phys), intent(out) :: canhs !canopy heat storage change (w/m2) - real (kind=kind_phys), intent(out) :: q2v real (kind=kind_phys) :: cah !sensible heat conductance, canopy air to zlvl air (m/s) real (kind=kind_phys) :: u10v !10 m wind speed in eastward dir (m/s) @@ -3857,22 +3849,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: ch2 !surface exchange at 2m real (kind=kind_phys) :: thstar !surface exchange at 2m - real (kind=kind_phys) :: dlf ! leaf dimension - real(kind=kind_phys) :: sigmaa ! momentum partition parameter - real(kind=kind_phys) :: kbsigmaf1 ! kb^-1 for fully convered by vegetation - real(kind=kind_phys) :: kbsigmafc ! kb^-1 under canopy ground - - real (kind=kind_phys) :: fm10 !monin-obukhov momentum adjustment at 10m - real (kind=kind_phys) :: rb1v !Bulk Richardson # over vegetation - real (kind=kind_phys) :: stress1v !Stress over vegetation - real (kind=kind_phys) :: snwd - real (kind=kind_phys) :: virtfacv - real (kind=kind_phys) :: thv1v - real (kind=kind_phys) :: tvsv - real (kind=kind_phys) :: tv1v - real (kind=kind_phys) :: zlvlv - - real (kind=kind_phys) :: thvair real (kind=kind_phys) :: thah real (kind=kind_phys) :: rahc2 !aerodynamic resistance for sensible heat (s/m) @@ -3885,14 +3861,10 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: e1 real (kind=kind_phys) :: hcv !canopy heat capacity j/m2/k, C.He added - real (kind=kind_phys) :: vaie !total leaf area index + stem area index,effective real (kind=kind_phys) :: laisune !sunlit leaf area index, one-sided (m2/m2),effective real (kind=kind_phys) :: laishae !shaded leaf area index, one-sided (m2/m2),effective - real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx,czil1 - real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 - integer :: k !index integer :: iter !iteration index @@ -3905,8 +3877,16 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & integer :: liter !last iteration - integer :: niter !for sfcdiff3 +! New variables for sfcdif3 + logical , intent(in ) :: thsfc_loc + real (kind=kind_phys), intent(in ) :: prslkix ! in exner function + real (kind=kind_phys), intent(in ) :: prsik1x ! in exner function + real (kind=kind_phys), intent(in ) :: prslk1x ! in exner function + real (kind=kind_phys), intent(in ) :: garea1 + real (kind=kind_phys), intent(in ) :: shdfac ! greeness vegetation fraction (-) + real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity + real (kind=kind_phys), intent( out) :: csigmaf1 ! real (kind=kind_phys) :: t, tdc !kelvin to degree celsius with limit -50 to +50 @@ -3918,11 +3898,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & mpe = 1e-6 liter = 0 - fv = ustarx - - niter = 1 - if (ur < 2.0) niter = 2 - ! --------------------------------------------------------------------------------------------- ! initialization variables that do not depend on stability iteration ! --------------------------------------------------------------------------------------------- @@ -3936,31 +3911,12 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & h = 0. qfx = 0. - csigmaf1 = 0. - ! limit lai vaie = min(6.,vai ) laisune = min(6.,laisun) laishae = min(6.,laisha) -! for sfcdiff3 - - snwd = snowh*1000.0 -! zlvlv = zlvl - zpd - zlvlv = zlvl - - virtfacv = 1.0 + 0.61 * max(qair, 1.e-8) - tv1v = sfctmp * virtfacv - - if(thsfc_loc) then ! Use local potential temperature - thv1v = sfctmp * prslkix * virtfacv - else ! Use potential temperature reference to 1000 hPa - thv1v = sfctmp / prslk1x * virtfacv - endif -! - - ! saturation vapor pressure at ground temperature t = tdc(tg) @@ -3975,8 +3931,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & qsfc = 0.622*eair/(psfc-0.378*eair) - dlf = parameters%dleaf !leaf dimension - ! canopy height hcan = parameters%hvt @@ -4024,37 +3978,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & air = -emv*(1.+(1.-emv)*(1.-emg))*lwdn - emv*emg*sb*tg**4 cir = (2.-emv*(1.-emg))*emv*sb ! --------------------------------------------------------------------------------------------- - - if (opt_trs == 1) then - z0h = z0m - elseif (opt_trs == 2) then -! z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) - czil1= 10.0 ** (- (0.40/0.07) * hcan) - z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) - elseif (opt_trs == 3) then - if (vegtyp.le.5) then - z0h = z0m - else - z0h = z0m*0.01 - endif - elseif (opt_trs == 4) then - sigmaa = 1.0 - (0.5/(0.5+vaie))*exp(-vaie**2/8.0) - kbsigmaf1 = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(dlf*ur/log((zlvl-zpd)/z0m)) - z0h = z0m/exp(kbsigmaf1) - csigmaf1 = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m)+kbsigmaf1) ! for output for interpolation - endif -! -- - tem1 = (z0m - z0lo) / (z0up - z0lo) - tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(shdfac, 0.1_kind_phys) - zvfun1= sqrt(tem1 * tem2) - gdx=sqrt(garea1) - if(opt_sfc == 1 .or. opt_sfc == 2) then - loop1: do iter = 1, niterc ! begin stability iteration -! use newly derived z0m/z0h - if(iter == 1) then z0hg = z0mg else @@ -4089,6 +4014,15 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & cm = cm / ur endif + if(opt_sfc == 3) then + call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in + zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in + zpd ,snowh ,fveg ,garea1 ,.true. ,vaie , & !in + ustarx ,fm ,fh ,fm2 ,fh2 , & !inout + z0h ,fv ,csigmaf1,cm ,ch ) !out + + endif + ramc = max(1.,1./(cm*ur)) rahc = max(1.,1./(ch*ur)) rawc = rahc @@ -4209,135 +4143,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & end do loop1 ! end stability iteration - endif !opt_sfc 1 or 2 -! -! sfcdiff3 -! - if (opt_sfc == 3) then - - z0hg = z0mg - - do iter = 1, niter !1 or 2; depending on ur - - if(thsfc_loc) then ! Use local potential temperature - tvsv = tah * virtfacv - else ! Use potential temperature referenced to 1000 hPa - tvsv = tah/prsik1x * virtfacv - endif - - call stability & - (zlvlv, zvfun1, gdx,tv1v,thv1v, ur, z0m, z0h, tvsv, grav,thsfc_loc, & - rb1v, fm,fh,fm10,fh2,cm,ch,stress1v,fv) - - ramc = max(1.,1./(cm*ur)) - rahc = max(1.,1./(ch*ur)) - rawc = rahc - -! aerodyn resistance between heights z0g and d+z0v, rag, and leaf -! boundary layer resistance, rb - - call ragrb(parameters,iter ,vaie ,rhoair ,hg ,tah , & !in - zpd ,z0mg ,z0hg ,hcan ,uc , & !in - z0h ,fv ,cwp ,vegtyp ,mpe , & !in - tv ,mozg ,fhg ,iloc ,jloc , & !inout - ramg ,rahg ,rawg ,rb ) !out - -! es and d(es)/dt evaluated at tv - - t = tdc(tv) - call esat(t, esatw, esati, dsatw, dsati) - if (t .gt. 0.) then - estv = esatw - destv = dsatw - else - estv = esati - destv = dsati - end if - -! stomatal resistance - - if(iter == 1) then - if (opt_crs == 1) then ! ball-berry - call stomata (parameters,vegtyp,mpe ,parsun ,foln ,iloc , jloc , & !in - tv ,estv ,eah ,sfctmp,sfcprs, & !in - o2air ,co2air,igs ,btran ,rb , & !in - rssun ,psnsun) !out - - call stomata (parameters,vegtyp,mpe ,parsha ,foln ,iloc , jloc , & !in - tv ,estv ,eah ,sfctmp,sfcprs, & !in - o2air ,co2air,igs ,btran ,rb , & !in - rssha ,psnsha) !out - end if - - if (opt_crs == 2) then ! jarvis - call canres (parameters,parsun,tv ,btran ,eah ,sfcprs, & !in - rssun ,psnsun,iloc ,jloc ) !out - - call canres (parameters,parsha,tv ,btran ,eah ,sfcprs, & !in - rssha ,psnsha,iloc ,jloc ) !out - end if - end if - -! prepare for sensible heat flux above veg. - - cah = 1./rahc - cvh = 2.*vaie/rb - cgh = 1./rahg - cond = cah + cvh + cgh - ata = (sfctmp*cah + tg*cgh) / cond - bta = cvh/cond - csh = (1.-bta)*rhoair*cpair*cvh - -! prepare for latent heat flux above veg. - - caw = 1./rawc - cew = fwet*vaie/rb - ctw = (1.-fwet)*(laisune/(rb+rssun) + laishae/(rb+rssha)) - cgw = 1./(rawg+rsurf) - cond = caw + cew + ctw + cgw - aea = (eair*caw + estg*cgw) / cond - bea = (cew+ctw)/cond - cev = (1.-bea)*cew*rhoair*cpair/gammav ! barlage: change to vegetation v3.6 - ctr = (1.-bea)*ctw*rhoair*cpair/gammav - -! evaluate surface fluxes with current temperature and solve for dts - - tah = ata + bta*tv ! canopy air t. - eah = aea + bea*estv ! canopy air e - - irc = fveg*(air + cir*tv**4) - shc = fveg*rhoair*cpair*cvh * ( tv-tah) - evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 - tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav - if (tv > tfrz) then - evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 - else - evc = min(canice*latheav/dt,evc) - end if - - b = sav-irc-shc-evc-tr+pahv !additional w/m2 - a = fveg*(4.*cir*tv**3 + csh + (cev+ctr)*destv) !volumetric heat capacity - dtv = b/a - - irc = irc + fveg*4.*cir*tv**3*dtv - shc = shc + fveg*csh*dtv - evc = evc + fveg*cev*destv*dtv - tr = tr + fveg*ctr*destv*dtv - -! update vegetation surface temperature - tv = tv + dtv -! tah = ata + bta*tv ! canopy air t; update here for consistency - -! for computing m-o length in the next iteration - h = rhoair*cpair*(tah - sfctmp) /rahc - hg = rhoair*cpair*(tg - tah) /rahg - -! consistent specific humidity from canopy air vapor pressure - qsfc = (0.622*eah)/(sfcprs-0.378*eah) - - enddo ! iteration - endif ! sfcdiff3 - ! under-canopy fluxes and tg air = - emg*(1.-emv)*lwdn - emg*emv*sb*tv**4 @@ -4443,7 +4248,8 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & #else tgb ,cm ,ch,ustarx, & !inout #endif - tauxb ,tauyb ,irb ,shb ,evb,csigmaf0,& !out + tauxb ,tauyb ,irb ,shb ,evb , & !out + csigmaf0, & !out ghb ,t2mb ,dx ,dz8w ,ivgtyp , & !out qc ,qsfc ,psfc , & !in sfcprs ,q2b ,ehb2 ) !in @@ -4489,14 +4295,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(in) :: rhsur !raltive humidity in surface soil/snow air space (-) real (kind=kind_phys), intent(in) :: fsno !snow fraction - logical , intent(in) :: thsfc_loc - real (kind=kind_phys) , intent(in) :: prslkix ! in exner function - real (kind=kind_phys) , intent(in) :: prsik1x ! in exner function - real (kind=kind_phys) , intent(in) :: prslk1x ! in exner function - real (kind=kind_phys) , intent(in) :: fveg - real (kind=kind_phys) , intent(in) :: shdfac - real (kind=kind_phys) , intent(in) :: garea1 - !jref:start; in integer , intent(in) :: ivgtyp real (kind=kind_phys) , intent(in) :: qc !cloud water mixing ratio @@ -4513,7 +4311,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(inout) :: tgb !ground temperature (k) real (kind=kind_phys), intent(inout) :: cm !momentum drag coefficient real (kind=kind_phys), intent(inout) :: ch !sensible heat exchange coefficient - real (kind=kind_phys), intent(inout) :: ustarx !friction velocity #ifdef CCPP character(len=*), intent(inout) :: errmsg integer, intent(inout) :: errflg @@ -4529,7 +4326,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(out) :: evb !latent heat flux (w/m2) [+ to atm] real (kind=kind_phys), intent(out) :: ghb !ground heat flux (w/m2) [+ to soil] real (kind=kind_phys), intent(out) :: t2mb !2 m height air temperature (k) - real (kind=kind_phys), intent(out) :: csigmaf0 ! !jref:start real (kind=kind_phys), intent(out) :: q2b !bare ground heat conductance real (kind=kind_phys) :: ehb !bare ground heat conductance @@ -4540,17 +4336,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & ! local variables - real (kind=kind_phys) :: rb1b !Bulk Richardson # over bare soil - real (kind=kind_phys) :: stress1b !Stress over bare soil - real (kind=kind_phys) :: snwd - real (kind=kind_phys) :: virtfacb - real (kind=kind_phys) :: thv1b - real (kind=kind_phys) :: tvsb - real (kind=kind_phys) :: tv1b - real (kind=kind_phys) :: zlvlb - - real (kind=kind_phys) :: fm10 - real (kind=kind_phys) :: taux !wind stress: e-w (n/m2) real (kind=kind_phys) :: tauy !wind stress: n-s (n/m2) real (kind=kind_phys) :: fira !total net longwave rad (w/m2) [+ to atm] @@ -4577,9 +4362,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys) :: cev !coefficients for ev as function of esat[ts] real (kind=kind_phys) :: cgh !coefficients for st as function of ts - real(kind=kind_phys) :: kbsigmaf0 - real(kind=kind_phys) :: reynb - !jref:start real (kind=kind_phys) :: rahb2 !aerodynamic resistance for sensible heat 2m (s/m) real (kind=kind_phys) :: rawb2 !aerodynamic resistance for water vapor 2m (s/m) @@ -4614,18 +4396,26 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys) :: fh2 !monin-obukhov heat adjustment at 2m real (kind=kind_phys) :: ch2 !surface exchange at 2m - real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx,czil1 - real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 - integer :: iter !iteration index integer :: niterb !number of iterations for surface temperature - integer :: niter - real (kind=kind_phys) :: mpe !prevents overflow error if division by zero !jref:start ! data niterb /3/ data niterb /5/ save niterb + +! New variables for sfcdif3 + + logical , intent(in ) :: thsfc_loc + real (kind=kind_phys), intent(in ) :: prslkix ! in exner function + real (kind=kind_phys), intent(in ) :: prsik1x ! in exner function + real (kind=kind_phys), intent(in ) :: prslk1x ! in exner function + real (kind=kind_phys), intent(in ) :: fveg + real (kind=kind_phys), intent(in ) :: shdfac + real (kind=kind_phys), intent(in ) :: garea1 + real (kind=kind_phys), intent(inout) :: ustarx !friction velocity + real (kind=kind_phys), intent( out) :: csigmaf0 ! + real (kind=kind_phys) :: t, tdc !kelvin to degree celsius with limit -50 to +50 tdc(t) = min( 50., max(-50.,(t-tfrz)) ) @@ -4641,69 +4431,10 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & h = 0. qfx = 0. - csigmaf0 = 0. - kbsigmaf0 = 0. - - niter = 1 - if (ur < 2.0) niter = 2 - - fv = ustarx - -! fv = ur*vkc/log((zlvl-zpd)/z0m) - - reynb = fv*z0m/(1.5e-05) - - if (reynb .gt. 2.0) then - kbsigmaf0 = 2.46*reynb**0.25 - log(7.4) - else - kbsigmaf0 = - log(0.397) - endif - - csigmaf0 = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m) + kbsigmaf0) - - if (opt_trs == 1) then - z0h = z0m - elseif (opt_trs == 2) then -! z0h = z0m*exp(-parameters%czil*0.4*258.2*sqrt(fv*z0m)) - czil1= 10.0 ** (- (0.40/0.07) * parameters%hvt) - z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) - elseif (opt_trs == 3) then - if (vegtyp.le.5) then - z0h = z0m - else - z0h = z0m*0.01 - endif - elseif (opt_trs == 4) then - z0h = max(z0m/exp(kbsigmaf0),1.0e-6) - endif -! -! for sfcdiff3; maybe should move to inside the option -! - snwd = snowh*1000.0 -! zlvlb = zlvl - zpd - zlvlb = zlvl - - virtfacb = 1.0 + 0.61 * max(qair, 1.e-8) - tv1b = sfctmp * virtfacb - - if(thsfc_loc) then ! Use local potential temperature - thv1b = sfctmp * prslkix * virtfacb - else ! Use potential temperature reference to 1000 hPa - thv1b = sfctmp / prslk1x * virtfacb - endif - cir = emg*sb cgh = 2.*df(isnow+1)/dzsnso(isnow+1) ! ----------------------------------------------------------------- - tem1 = (z0m - z0lo) / (z0up - z0lo) - tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(shdfac, 0.1_kind_phys) - zvfun1= sqrt(tem1 * tem2) - gdx=sqrt(garea1) - - if (opt_sfc == 1 .or. opt_sfc == 2) then - loop3: do iter = 1, niterb ! begin stability iteration ! if(iter == 1) then @@ -4743,6 +4474,15 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & endif + if(opt_sfc == 3) then + call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in + zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in + zpd ,snowh ,shdfac ,garea1 ,.false. ,0.0 , & !in + ustarx ,fm ,fh ,fm2 ,fh2 , & !inout + z0h ,fv ,csigmaf0,cm ,ch ) !out + + endif + ramb = max(1.,1./(cm*ur)) rahb = max(1.,1./(ch*ur)) rawb = rahb @@ -4800,83 +4540,8 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & qfx = (qsfc-qair)*cev*gamma/cpair end do loop3 ! end stability iteration - endif ! opt_sfc 1/2 ! ----------------------------------------------------------------- - if (opt_sfc == 3) then - - do iter = 1, niter !1 or 2; depending on ur - - if(thsfc_loc) then ! Use local potential temperature - tvsb = tgb * virtfacb - else ! Use potential temperature referenced to 1000 hPa - tvsb = tgb/prsik1x * virtfacb - endif - - call stability & - (zlvlb, zvfun1, gdx,tv1b,thv1b, ur, z0m, z0h, tvsb, grav,thsfc_loc, & - rb1b, fm,fh,fm10,fh2,cm,ch,stress1b,fv) - - - ramb = max(1.,1./(cm*ur)) - rahb = max(1.,1./(ch*ur)) - rawb = rahb - -!jref - variables for diagnostics - emb = 1./ramb - ehb = 1./rahb - -! es and d(es)/dt evaluated at tg - - t = tdc(tgb) - call esat(t, esatw, esati, dsatw, dsati) - if (t .gt. 0.) then - estg = esatw - destg = dsatw - else - estg = esati - destg = dsati - end if - - csh = rhoair*cpair/rahb - cev = rhoair*cpair/gamma/(rsurf+rawb) - -! surface fluxes and dtg - - irb = cir * tgb**4 - emg*lwdn - shb = csh * (tgb - sfctmp ) - evb = cev * (estg*rhsur - eair ) - ghb = cgh * (tgb - stc(isnow+1)) - - b = sag-irb-shb-evb-ghb+pahb - a = 4.*cir*tgb**3 + csh + cev*destg + cgh - dtg = b/a - - irb = irb + 4.*cir*tgb**3*dtg - shb = shb + csh*dtg - evb = evb + cev*destg*dtg - ghb = ghb + cgh*dtg - -! update ground surface temperature - tgb = tgb + dtg - -! for m-o length -! h = csh * (tgb - sfctmp) - - t = tdc(tgb) - call esat(t, esatw, esati, dsatw, dsati) - if (t .gt. 0.) then - estg = esatw - else - estg = esati - end if - qsfc = 0.622*(estg*rhsur)/(psfc-0.378*(estg*rhsur)) - - qfx = (qsfc-qair)*cev*gamma/cpair - - end do ! end stability iteration - endif ! sfcdiff3 - ! if snow on ground and tg > tfrz: reset tg = tfrz. reevaluate ground fluxes. if(opt_stc == 1 .or. opt_stc == 3) then @@ -5409,6 +5074,126 @@ subroutine sfcdif2(parameters,iter ,z0 ,thz0 ,thlm ,sfcspd , & !in ! ---------------------------------------------------------------------- end subroutine sfcdif2 +!== begin sfcdif3 ================================================================================== + +!>\ingroup NoahMP_LSM +!! compute surface drag coefficient cm for momentum and ch for heat. + subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in + zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in + zpd ,snowh ,fveg ,garea1 ,vegetated,vaie , & !in + ustarx ,fm ,fh ,fm2 ,fh2 , & !inout + z0h ,fv ,csigmaf ,cm ,ch ) !out + +! ------------------------------------------------------------------------------------------------- +! computing surface drag coefficient cm for momentum and ch for heat +! ------------------------------------------------------------------------------------------------- + implicit none +! ------------------------------------------------------------------------------------------------- +! inputs + + type (noahmp_parameters), intent(in) :: parameters + integer, intent(in ) :: iloc ! grid index + integer, intent(in ) :: jloc ! grid index + integer, intent(in ) :: iter ! iteration index + real (kind=kind_phys), intent(in ) :: sfctmp ! temperature at reference height [K] + real (kind=kind_phys), intent(in ) :: qair ! specific humidity at reference height [kg/kg] + real (kind=kind_phys), intent(in ) :: ur ! wind speed [m/s] + real (kind=kind_phys), intent(in ) :: zlvl ! reference height [m] + real (kind=kind_phys), intent(in ) :: tgb ! ground temperature [K] + logical, intent(in ) :: thsfc_loc ! flag for using sfc-based theta + real (kind=kind_phys), intent(in ) :: prslkix ! in exner function + real (kind=kind_phys), intent(in ) :: prsik1x ! in exner function + real (kind=kind_phys), intent(in ) :: prslk1x ! in exner function + real (kind=kind_phys), intent(in ) :: z0m ! roughness length, momentum, ground [m] + real (kind=kind_phys), intent(in ) :: zpd ! zero plane displacement [m] + real (kind=kind_phys), intent(in ) :: snowh ! snow depth [m] + real (kind=kind_phys), intent(in ) :: fveg ! fractional vegetation cover + real (kind=kind_phys), intent(in ) :: garea1 ! grid area [km2] + logical, intent(in ) :: vegetated ! .true. if vegetated + real (kind=kind_phys), intent(in ) :: vaie ! vegetation area index [m2/m2] + real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity [m/s] + real (kind=kind_phys), intent(inout) :: fm ! momentum stability correction, weighted by prior iters + real (kind=kind_phys), intent(inout) :: fh ! sen heat stability correction, weighted by prior iters + real (kind=kind_phys), intent(inout) :: fm2 ! sen heat stability correction, weighted by prior iters + real (kind=kind_phys), intent(inout) :: fh2 ! sen heat stability correction, weighted by prior iters + real (kind=kind_phys), intent( out) :: z0h ! roughness length, sensible heat, ground [m] + real (kind=kind_phys), intent( out) :: fv ! friction velocity (m/s) + real (kind=kind_phys), intent( out) :: csigmaf ! + real (kind=kind_phys), intent( out) :: cm ! drag coefficient for momentum + real (kind=kind_phys), intent( out) :: ch ! drag coefficient for heat + + real (kind=kind_phys) :: reyn ! reynolds number + real (kind=kind_phys) :: kbsigmaf ! kb factor + real (kind=kind_phys) :: snwd ! snow depth [mm] + real (kind=kind_phys) :: zlvlb ! reference height - zpd [m] + real (kind=kind_phys) :: virtfac ! virtual temperature factor [-] + real (kind=kind_phys) :: tv1 ! virtual temperature at reference [K] + real (kind=kind_phys) :: thv1 ! virtual theta at reference [K] + real (kind=kind_phys) :: tvs ! virtural surface temperature [K] + real (kind=kind_phys) :: rb1 ! bulk Richardson - stability output + real (kind=kind_phys) :: stress1 ! stress - stability output + real (kind=kind_phys) :: fm10 ! 10-m stability adjustment - stability output + real (kind=kind_phys) :: dlf ! leaf dimension + real (kind=kind_phys) :: sigmaa ! momentum partition parameter + real (kind=kind_phys) :: tem1,tem2,zvfun1,gdx + real (kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 + +! ------------------------------------------------------------------------------------------------- + + fv = ustarx +! fv = ur*vkc/log((zlvl-zpd)/z0m) + + if(vegetated) then + + dlf = parameters%dleaf + sigmaa = 1.0 - (0.5/(0.5+vaie)) * exp(-vaie**2/8.0) + kbsigmaf = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(dlf*ur/log((zlvl-zpd)/z0m)) + z0h = z0m/exp(kbsigmaf) + csigmaf = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m) + kbsigmaf) ! for output for interpolation + + else + + reyn = fv*z0m/(1.5e-05) + if (reyn .gt. 2.0) then + kbsigmaf = 2.46*reyn**0.25 - log(7.4) + else + kbsigmaf = - log(0.397) + endif + + z0h = max(z0m/exp(kbsigmaf),1.0e-6) + csigmaf = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m) + kbsigmaf) + + end if + + snwd = snowh*1000.0 + zlvlb = zlvl! - zpd + + virtfac = 1.0 + 0.61 * max(qair, 1.0e-8) + tv1 = sfctmp * virtfac + + if(thsfc_loc) then ! Use local potential temperature + thv1 = sfctmp * prslkix * virtfac + else ! Use potential temperature reference to 1000 hPa + thv1 = sfctmp / prslk1x * virtfac + endif + + tem1 = (z0m - z0lo) / (z0up - z0lo) + tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) + tem2 = max(fveg, 0.1_kind_phys) + zvfun1 = sqrt(tem1 * tem2) + gdx = sqrt(garea1) + + if(thsfc_loc) then ! Use local potential temperature + tvs = tgb * virtfac + else ! Use potential temperature referenced to 1000 hPa + tvs = tgb/prsik1x * virtfac + endif + + call stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & + rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) + + end subroutine sfcdif3 + !== begin esat ===================================================================================== !>\ingroup NoahMP_LSM From c50f50a07c78e9eed773f8c936d8360b61a1d5c9 Mon Sep 17 00:00:00 2001 From: barlage Date: Mon, 7 Mar 2022 14:42:58 -0700 Subject: [PATCH 06/21] change fveg to shdfac in sfcdif3 vege call --- physics/module_sf_noahmplsm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 0248a116b..5964e4575 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -4017,7 +4017,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,fveg ,garea1 ,.true. ,vaie , & !in + zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf1,cm ,ch ) !out From 53c0c7acfe64609c021e551f1edf67376bcd1387 Mon Sep 17 00:00:00 2001 From: barlage Date: Mon, 7 Mar 2022 15:07:35 -0700 Subject: [PATCH 07/21] move trs options to sfcdif3 --- physics/module_sf_noahmplsm.f90 | 44 ++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 5964e4575..42a213fed 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -3986,6 +3986,19 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & z0hg = z0mg !* exp(-czil*0.4*258.2*sqrt(fv*z0mg)) end if + if (opt_trs == 1) then + z0h = z0m + elseif (opt_trs == 2) then + czil1= 10.0 ** (- (0.40/0.07) * parameters%hvt) + z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) + elseif (opt_trs == 3) then + if (vegtyp.le.5) then + z0h = z0m + else + z0h = z0m*0.01 + endif + endif + ! aerodyn resistances between heights zlvl and d+z0v if(opt_sfc == 1) then @@ -4017,7 +4030,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie , & !in + zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,vegtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf1,cm ,ch ) !out @@ -4477,7 +4490,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.false. ,0.0 , & !in + zpd ,snowh ,shdfac ,garea1 ,.false. ,0.0 ,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out @@ -5080,7 +5093,7 @@ end subroutine sfcdif2 !! compute surface drag coefficient cm for momentum and ch for heat. subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,fveg ,garea1 ,vegetated,vaie , & !in + zpd ,snowh ,fveg ,garea1 ,vegetated,vaie ,vegtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf ,cm ,ch ) !out @@ -5111,6 +5124,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: garea1 ! grid area [km2] logical, intent(in ) :: vegetated ! .true. if vegetated real (kind=kind_phys), intent(in ) :: vaie ! vegetation area index [m2/m2] + integer , intent(in ) :: vegtyp ! vegetation category real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm ! momentum stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fh ! sen heat stability correction, weighted by prior iters @@ -5132,8 +5146,8 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys) :: tvs ! virtural surface temperature [K] real (kind=kind_phys) :: rb1 ! bulk Richardson - stability output real (kind=kind_phys) :: stress1 ! stress - stability output + real (kind=kind_phys) :: czil1 ! canopy based czil real (kind=kind_phys) :: fm10 ! 10-m stability adjustment - stability output - real (kind=kind_phys) :: dlf ! leaf dimension real (kind=kind_phys) :: sigmaa ! momentum partition parameter real (kind=kind_phys) :: tem1,tem2,zvfun1,gdx real (kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 @@ -5145,11 +5159,23 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur if(vegetated) then - dlf = parameters%dleaf - sigmaa = 1.0 - (0.5/(0.5+vaie)) * exp(-vaie**2/8.0) - kbsigmaf = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(dlf*ur/log((zlvl-zpd)/z0m)) - z0h = z0m/exp(kbsigmaf) - csigmaf = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m) + kbsigmaf) ! for output for interpolation + if (opt_trs == 1) then + z0h = z0m + elseif (opt_trs == 2) then + czil1= 10.0 ** (- (0.40/0.07) * parameters%hvt) + z0h = z0m*exp(-czil1*0.4*258.2*sqrt(fv*z0m)) + elseif (opt_trs == 3) then + if (vegtyp.le.5) then + z0h = z0m + else + z0h = z0m*0.01 + endif + elseif (opt_trs == 4) then + sigmaa = 1.0 - (0.5/(0.5+vaie))*exp(-vaie**2/8.0) + kbsigmaf = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(parameters%dleaf*ur/log((zlvl-zpd)/z0m)) + z0h = z0m/exp(kbsigmaf) + csigmaf = log((zlvl-zpd)/z0m)*(log((zlvl-zpd)/z0m)+kbsigmaf) ! for output for interpolation + endif else From f093f77d40f4d7e0c5097caa20191050964ef5d5 Mon Sep 17 00:00:00 2001 From: barlage Date: Mon, 7 Mar 2022 15:10:50 -0700 Subject: [PATCH 08/21] fix missing czil1 in vege_flux --- physics/module_sf_noahmplsm.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 42a213fed..4a296debb 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -3887,6 +3887,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in ) :: shdfac ! greeness vegetation fraction (-) real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity real (kind=kind_phys), intent( out) :: csigmaf1 ! + real (kind=kind_phys) :: czil1 ! canopy based czil real (kind=kind_phys) :: t, tdc !kelvin to degree celsius with limit -50 to +50 From 13a1d1480410c3c7b453e252f2d3c159e48cb04f Mon Sep 17 00:00:00 2001 From: barlage Date: Tue, 8 Mar 2022 06:55:33 -0700 Subject: [PATCH 09/21] add some clean up to energy --- physics/module_sf_noahmplsm.f90 | 65 +++++++++++++++++---------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 4a296debb..4ff484dfb 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -1952,26 +1952,21 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in chv2 = 0. rb = 0. -! - cdmnv = 0. - ezpdv = 0. - - cdmng = 0. - ezpdg = 0. - - cdmn = 0. - ezpd = 0. - - gsigma = 0. - - z0hwrf = 0. - csigmaf1 = 0. - csigmaf0 = 0. - csigmafveg= 0. - kbsigmafveg = 0. - aone = 0. - coeffa = 0. - coeffb = 0. + cdmnv = 0.0 + ezpdv = 0.0 + cdmng = 0.0 + ezpdg = 0.0 + cdmn = 0.0 + ezpd = 0.0 + gsigma = 0.0 + z0hwrf = 0.0 + csigmaf1 = 0.0 + csigmaf0 = 0.0 + csigmafveg= 0.0 + kbsigmafveg = 0.0 + aone = 0.0 + coeffa = 0.0 + coeffb = 0.0 ! @@ -2190,9 +2185,11 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in qc ,qsfc ,psfc , & !in q2v ,chv2, chleaf, chuc) !inout - cdmnv = 0.4*0.4/log((zlvl-zpd)/z0m)**2 - aone = 2.6*(10.0*parameters%hvt/(zlvl-zpd))**0.355 - ezpdv = zpd*fveg !for the grid +! new coupling code + + cdmnv = 0.4*0.4/log((zlvl-zpd)/z0m)**2 + aone = 2.6*(10.0*parameters%hvt/(zlvl-zpd))**0.355 + ezpdv = zpd*fveg !for the grid !jref:end #ifdef CCPP @@ -2221,18 +2218,20 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in qc ,qsfc ,psfc , & !in sfcprs ,q2b, chb2) !in - cdmng = 0.4*0.4/log((zlvl-zpdg)/z0mg)**2 - ezpdg = zpdg +! new coupling code + + cdmng = 0.4*0.4/log((zlvl-zpdg)/z0mg)**2 + ezpdg = zpdg ! ! vegetation is optional; use the larger one ! - if (ezpdv .ge. ezpdg ) then - ezpd = ezpdv - elseif (ezpdv .gt. 0.0 .and. ezpdv .lt. ezpdg) then - ezpd = (1.0 -fveg)*ezpdg - else - ezpd = ezpdg - endif + if (ezpdv .ge. ezpdg ) then + ezpd = ezpdv + elseif (ezpdv .gt. 0.0 .and. ezpdv .lt. ezpdg) then + ezpd = (1.0 -fveg)*ezpdg + else + ezpd = ezpdg + endif !jref:end #ifdef CCPP @@ -2260,6 +2259,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in q1 = fveg * (eah*0.622/(sfcprs - 0.378*eah)) + (1.0 - fveg)*qsfc q2e = fveg * q2v + (1.0 - fveg) * q2b +! new coupling code + if (opt_trs == 1) then z0wrf = fveg * z0m + (1.0 - fveg) * z0mg z0hwrf = z0wrf From ebb4fa16d3d2494850431a97b3772e60875d8975 Mon Sep 17 00:00:00 2001 From: barlage Date: Tue, 8 Mar 2022 07:01:12 -0700 Subject: [PATCH 10/21] add some groundwater mods from ncar code --- physics/module_sf_noahmplsm.f90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 4ff484dfb..445034741 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -7621,8 +7621,10 @@ subroutine soilwater (parameters,nsoil ,nsnow ,dt ,zsoil ,dzsnso , & !in if ( parameters%urban_flag ) fcr(1)= 0.95 if(opt_run == 1) then - fff = 6.0 - fsat = parameters%fsatmx*exp(-0.5*fff*(zwt-2.0)) +! fff = 6.0 + fff = parameters%bexp(1) / 3.0 ! calibratable, c.he changed based on gy niu's update +! fsat = parameters%fsatmx*exp(-0.5*fff*(zwt-2.0)) + fsat = parameters%fsatmx*exp(-0.5*fff*zwt) ! c.he changed based on gy niu's update if(qinsur > 0.) then runsrf = qinsur * ( (1.0-fcr(1))*fsat + fcr(1) ) pddum = qinsur - runsrf ! m/s @@ -8337,8 +8339,9 @@ subroutine groundwater(parameters,nsnow ,nsoil ,dt ,sice ,zsoil , & !in real (kind=kind_phys) :: watmin!minimum soil vol soil moisture [m3/m3] real (kind=kind_phys) :: xs !excessive water above saturation [mm] real (kind=kind_phys), parameter :: rous = 0.2 !specific yield [-] - real (kind=kind_phys), parameter :: cmic = 0.20 !microprore content (0.0-1.0) +! real (kind=kind_phys), parameter :: cmic = 0.20 !microprore content (0.0-1.0) !0.0-close to free drainage + real (kind=kind_phys), parameter :: cmic = 0.80 ! calibratable, c.he changed based on gy niu's update ! ------------------------------------------------------------- qdis = 0.0 qin = 0.0 @@ -8380,8 +8383,10 @@ subroutine groundwater(parameters,nsnow ,nsoil ,dt ,sice ,zsoil , & !in ! groundwater discharge [mm/s] - fff = 6.0 - rsbmx = 5.0 +! fff = 6.0 +! rsbmx = 5.0 + fff = parameters%bexp(iwt) / 3.0 ! calibratable, c.he changed based on gy niu's update + rsbmx = hk(iwt) * 1.0e3 * exp(3.0) ! mm/s, calibratable, c.he changed based on gy niu's update qdis = (1.0-fcrmax)*rsbmx*exp(-parameters%timean)*exp(-fff*(zwt-2.0)) From 41cf4ecb44a6a983a8132c2986b91a0c2964595b Mon Sep 17 00:00:00 2001 From: helin wei Date: Tue, 8 Mar 2022 14:31:07 +0000 Subject: [PATCH 11/21] gvf impact on thermal conductivity limited to the first soil layer --- physics/module_sf_noahmplsm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 4a296debb..0601e98f1 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2503,7 +2503,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) - df = df * exp (sbeta * fveg) + df(1) = df(1) * exp (sbeta * fveg) ! compute lake thermal properties ! (no consideration of turbulent mixing for this version) From c1d813e21bd5238c65c95974264965e2f01540f6 Mon Sep 17 00:00:00 2001 From: helin wei Date: Tue, 8 Mar 2022 19:19:40 +0000 Subject: [PATCH 12/21] correct the reference height --- physics/module_sf_noahmp_glacier.f90 | 3 +-- physics/module_sf_noahmplsm.f90 | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90 index 1ea4a45b8..c4c03aaf8 100644 --- a/physics/module_sf_noahmp_glacier.f90 +++ b/physics/module_sf_noahmp_glacier.f90 @@ -1152,8 +1152,7 @@ subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ! the following only applies to opt_sfc =3, opt_sfc = 1 still done its old way snwd = snowh*1000.0 -! zlvli = zlvl - zpd - zlvli = zlvl + zlvli = zlvl - zpd ! fv = ustarx ! the input maybe too high for glacial fv = ur*vkc/log(zlvli/z0m) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index a93284475..919d81507 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -5194,7 +5194,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur end if snwd = snowh*1000.0 - zlvlb = zlvl! - zpd + zlvlb = zlvl - zpd virtfac = 1.0 + 0.61 * max(qair, 1.0e-8) tv1 = sfctmp * virtfac From 11b50ca1f939042faf1ec3d10707fb73c2af5f4b Mon Sep 17 00:00:00 2001 From: helin wei Date: Wed, 9 Mar 2022 18:21:04 +0000 Subject: [PATCH 13/21] to read new hig-res ice climatology data --- physics/sfcsub.F | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/sfcsub.F b/physics/sfcsub.F index e8b61f083..cdc91cca9 100644 --- a/physics/sfcsub.F +++ b/physics/sfcsub.F @@ -34,7 +34,8 @@ module sfccyc_module integer, parameter :: kpdalf(2)=(/214,217/) ! real (kind=kind_io8), parameter :: ten=10.0, one=1.0, zero=0.0 - integer, parameter :: xdata=5000, ydata=2500, mdata=xdata*ydata +! integer, parameter :: xdata=5000, ydata=2500, mdata=xdata*ydata + integer, parameter :: xdata=7200, ydata=3600, mdata=xdata*ydata integer :: veg_type_landice integer :: soil_type_landice integer :: num_threads From 4ed3982e48694e9abe52ae2dbf48d79068484c43 Mon Sep 17 00:00:00 2001 From: helin wei Date: Thu, 10 Mar 2022 17:49:49 +0000 Subject: [PATCH 14/21] replace fveg by lai/laimax to be used for dependent --- physics/module_sf_noahmplsm.f90 | 60 ++++++++++++++++++++++++--------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 919d81507..7e17f511d 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -682,6 +682,10 @@ subroutine noahmp_sflx (parameters, & logical :: crop_active !< flag to run crop model ! add canopy heat storage (C.He added based on GY Niu's communication) real :: canhs ! canopy heat storage change w/m2 +! maximum lai/sai used for some parameterizations based on plant growthi + real (kind=kind_phys) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided + ! intent (out) variables need to be assigned a value. these normally get assigned values ! only if dveg == 2. @@ -732,7 +736,7 @@ subroutine noahmp_sflx (parameters, & ! vegetation phenology call phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yearlen , julian , & !in - lai , sai , troot , elai , esai ,igs, pgs) + lai , sai , laimax, saimax, troot , elai , esai ,igs, pgs) !input gvf should be consistent with lai if(dveg == 1 .or. dveg == 6 .or. dveg == 7) then @@ -776,7 +780,7 @@ subroutine noahmp_sflx (parameters, & sfctmp ,thair ,lwdn ,uu ,vv ,zlvl , & !in co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in - elai ,esai ,fwet ,foln , & !in + elai ,esai ,laimax, saimax, fwet ,foln , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in @@ -1055,7 +1059,7 @@ end subroutine atm !!vegetation phenology considering vegetation canopy being buried by snow and !!evolution in time. subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yearlen , julian , & !in - lai , sai , troot , elai , esai , igs, pgs) + lai , sai , laimax, saimax, troot , elai , esai , igs, pgs) ! -------------------------------------------------------------------------------------------------- ! vegetation phenology considering vegeation canopy being buries by snow and evolution in time @@ -1076,6 +1080,8 @@ subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yea real (kind=kind_phys) , intent(inout) :: sai !sai, unadjusted for burying by snow ! outputs + real (kind=kind_phys) , intent(out ) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys) , intent(out ) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys) , intent(out ) :: elai !leaf area index, after burying by snow real (kind=kind_phys) , intent(out ) :: esai !stem area index, after burying by snow real (kind=kind_phys) , intent(out ) :: igs !growing season index (0=off, 1=on) @@ -1095,6 +1101,23 @@ subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yea real (kind=kind_phys) :: t !current month (1.00, ..., 12.00) ! -------------------------------------------------------------------------------------------------- +! derive monthly maximum lai and sai from monthly lai + + laimax=parameters%laim(1) + saimax=parameters%saim(1) + + do k=1,12 + + if(parameters%laim(k).ge.laimax)then + laimax=parameters%laim(k) + endif + + if(parameters%saim(k).ge.saimax)then + saimax=parameters%saim(k) + endif + + enddo + if (croptype == 0) then if ( dveg == 1 .or. dveg == 3 .or. dveg == 4 ) then @@ -1614,7 +1637,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in sfctmp ,thair ,lwdn ,uu ,vv ,zref , & !in co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in - elai ,esai ,fwet ,foln , & !in + elai ,esai ,laimax, saimax, fwet ,foln , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in @@ -1709,6 +1732,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: cosz !cosine solar zenith angle (0-1) real (kind=kind_phys) , intent(in) :: elai !lai adjusted for burying by snow real (kind=kind_phys) , intent(in) :: esai !lai adjusted for burying by snow + real (kind=kind_phys) , intent(in) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys) , intent(in) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys) , intent(in) :: fwet !fraction of canopy that is wet [-] real (kind=kind_phys) , intent(in) :: fveg !greeness vegetation fraction (-) real (kind=kind_phys) , intent(in) :: shdfac !< green vegetation fraction [0.0-1.0] @@ -2039,7 +2064,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in call thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , fveg, & !in + lat ,z0m ,zlvl ,vegtyp , elai,laimax, & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out @@ -2163,7 +2188,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in dt ,sav ,sag ,lwdn ,ur , & !in uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag , & !in - fwet ,laisun ,laisha ,cwp ,dzsnso , & !in + laimax, saimax,fwet ,laisun ,laisha ,cwp ,dzsnso , & !in zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in @@ -2429,7 +2454,7 @@ end subroutine energy subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , fveg,& !in + lat ,z0m ,zlvl ,vegtyp , elai, laimax,& !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out ! ------------------------------------------------------------------------------------------------- @@ -2454,8 +2479,9 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), intent(in) :: lat !latitude (radians) real (kind=kind_phys), intent(in) :: z0m !roughness length (m) real (kind=kind_phys), intent(in) :: zlvl !reference height (m) - integer , intent(in) :: vegtyp !vegtyp type - real (kind=kind_phys), intent(in) :: fveg !green vegetation fraction [0.0-1.0] + real (kind=kind_phys), intent(in) :: elai !lai adjusted for burying by snow + real (kind=kind_phys), intent(in) :: laimax !< monthly maximum leaf area index, one-sided + integer , intent(in) :: vegtyp !vegtyp type ! outputs real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: df !thermal conductivity [w/m/k] @@ -2504,7 +2530,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) - df(1) = df(1) * exp (sbeta * fveg) + df(1) = df(1) * exp (sbeta * elai/laimax) ! compute lake thermal properties ! (no consideration of turbulent mixing for this version) @@ -3650,7 +3676,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & dt ,sav ,sag ,lwdn ,ur , & !in uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag, & !in - fwet ,laisun ,laisha ,cwp ,dzsnso , & !in + laimax, saimax,fwet ,laisun ,laisha ,cwp ,dzsnso , & !in zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in @@ -3706,6 +3732,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: fsno !snow fraction real (kind=kind_phys), intent(in) :: snowh !actual snow depth [m] + real (kind=kind_phys), intent(in) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys), intent(in) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys), intent(in) :: fwet !wetted fraction of canopy real (kind=kind_phys), intent(in) :: cwp !canopy wind parameter @@ -4032,7 +4060,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,vegtyp , & !in + zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,laimax,saimax,vegtyp, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf1,cm ,ch ) !out @@ -4492,7 +4520,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.false. ,0.0 ,ivgtyp , & !in + zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,4.5,1.4,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out @@ -5095,7 +5123,7 @@ end subroutine sfcdif2 !! compute surface drag coefficient cm for momentum and ch for heat. subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,fveg ,garea1 ,vegetated,vaie ,vegtyp , & !in + zpd ,snowh ,fveg ,garea1 ,vegetated,vaie,laimax,saimax,vegtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf ,cm ,ch ) !out @@ -5126,6 +5154,8 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: garea1 ! grid area [km2] logical, intent(in ) :: vegetated ! .true. if vegetated real (kind=kind_phys), intent(in ) :: vaie ! vegetation area index [m2/m2] + real (kind=kind_phys), intent(in ) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys), intent(in ) :: laimax !< monthly maximum leaf area index, one-sided integer , intent(in ) :: vegtyp ! vegetation category real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm ! momentum stability correction, weighted by prior iters @@ -5207,7 +5237,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur tem1 = (z0m - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(fveg, 0.1_kind_phys) + tem2 = max(vaie/(laimax+saimax), 0.1_kind_phys) zvfun1 = sqrt(tem1 * tem2) gdx = sqrt(garea1) From 8e1b316051e6039101e8c3173a5af9dd2df63590 Mon Sep 17 00:00:00 2001 From: helin wei Date: Thu, 10 Mar 2022 19:11:01 +0000 Subject: [PATCH 15/21] simplify the code with internal function maxval --- physics/module_sf_noahmplsm.f90 | 56 ++++++++++----------------------- 1 file changed, 17 insertions(+), 39 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 7e17f511d..c945e66ff 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -683,8 +683,6 @@ subroutine noahmp_sflx (parameters, & ! add canopy heat storage (C.He added based on GY Niu's communication) real :: canhs ! canopy heat storage change w/m2 ! maximum lai/sai used for some parameterizations based on plant growthi - real (kind=kind_phys) :: saimax !< monthly maximum stem area index, one-sided - real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided ! intent (out) variables need to be assigned a value. these normally get assigned values @@ -736,7 +734,7 @@ subroutine noahmp_sflx (parameters, & ! vegetation phenology call phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yearlen , julian , & !in - lai , sai , laimax, saimax, troot , elai , esai ,igs, pgs) + lai , sai , troot , elai , esai ,igs, pgs) !input gvf should be consistent with lai if(dveg == 1 .or. dveg == 6 .or. dveg == 7) then @@ -780,7 +778,7 @@ subroutine noahmp_sflx (parameters, & sfctmp ,thair ,lwdn ,uu ,vv ,zlvl , & !in co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in - elai ,esai ,laimax, saimax, fwet ,foln , & !in + elai ,esai ,fwet ,foln , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in @@ -1059,7 +1057,7 @@ end subroutine atm !!vegetation phenology considering vegetation canopy being buried by snow and !!evolution in time. subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yearlen , julian , & !in - lai , sai , laimax, saimax, troot , elai , esai , igs, pgs) + lai , sai , troot , elai , esai , igs, pgs) ! -------------------------------------------------------------------------------------------------- ! vegetation phenology considering vegeation canopy being buries by snow and evolution in time @@ -1080,8 +1078,6 @@ subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yea real (kind=kind_phys) , intent(inout) :: sai !sai, unadjusted for burying by snow ! outputs - real (kind=kind_phys) , intent(out ) :: saimax !< monthly maximum stem area index, one-sided - real (kind=kind_phys) , intent(out ) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys) , intent(out ) :: elai !leaf area index, after burying by snow real (kind=kind_phys) , intent(out ) :: esai !stem area index, after burying by snow real (kind=kind_phys) , intent(out ) :: igs !growing season index (0=off, 1=on) @@ -1101,23 +1097,6 @@ subroutine phenology (parameters,vegtyp ,croptype, snowh , tv , lat , yea real (kind=kind_phys) :: t !current month (1.00, ..., 12.00) ! -------------------------------------------------------------------------------------------------- -! derive monthly maximum lai and sai from monthly lai - - laimax=parameters%laim(1) - saimax=parameters%saim(1) - - do k=1,12 - - if(parameters%laim(k).ge.laimax)then - laimax=parameters%laim(k) - endif - - if(parameters%saim(k).ge.saimax)then - saimax=parameters%saim(k) - endif - - enddo - if (croptype == 0) then if ( dveg == 1 .or. dveg == 3 .or. dveg == 4 ) then @@ -1637,7 +1616,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in sfctmp ,thair ,lwdn ,uu ,vv ,zref , & !in co2air ,o2air ,solad ,solai ,cosz ,igs , & !in eair ,tbot ,zsnso ,zsoil , & !in - elai ,esai ,laimax, saimax, fwet ,foln , & !in + elai ,esai ,fwet ,foln , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in @@ -1732,8 +1711,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: cosz !cosine solar zenith angle (0-1) real (kind=kind_phys) , intent(in) :: elai !lai adjusted for burying by snow real (kind=kind_phys) , intent(in) :: esai !lai adjusted for burying by snow - real (kind=kind_phys) , intent(in) :: saimax !< monthly maximum stem area index, one-sided - real (kind=kind_phys) , intent(in) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys) , intent(in) :: fwet !fraction of canopy that is wet [-] real (kind=kind_phys) , intent(in) :: fveg !greeness vegetation fraction (-) real (kind=kind_phys) , intent(in) :: shdfac !< green vegetation fraction [0.0-1.0] @@ -2064,7 +2041,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in call thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , elai,laimax, & !in + lat ,z0m ,zlvl ,vegtyp , elai, & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out @@ -2188,7 +2165,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in dt ,sav ,sag ,lwdn ,ur , & !in uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag , & !in - laimax, saimax,fwet ,laisun ,laisha ,cwp ,dzsnso , & !in + fwet ,laisun ,laisha ,cwp ,dzsnso , & !in zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in @@ -2454,7 +2431,7 @@ end subroutine energy subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , elai, laimax,& !in + lat ,z0m ,zlvl ,vegtyp , elai, & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out ! ------------------------------------------------------------------------------------------------- @@ -2480,7 +2457,6 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), intent(in) :: z0m !roughness length (m) real (kind=kind_phys), intent(in) :: zlvl !reference height (m) real (kind=kind_phys), intent(in) :: elai !lai adjusted for burying by snow - real (kind=kind_phys), intent(in) :: laimax !< monthly maximum leaf area index, one-sided integer , intent(in) :: vegtyp !vegtyp type ! outputs @@ -2498,6 +2474,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), dimension(-nsnow+1: 0) :: tksno !snow thermal conductivity (j/m3/k) real (kind=kind_phys), dimension( 1:nsoil) :: sice !soil ice content real (kind=kind_phys), parameter :: sbeta = -2.0 + real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided ! -------------------------------------------------------------------------------------------------- ! compute snow thermal conductivity and heat capacity @@ -2530,6 +2507,7 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) + laimax = maxval(parameters%laim) df(1) = df(1) * exp (sbeta * elai/laimax) ! compute lake thermal properties @@ -3676,7 +3654,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & dt ,sav ,sag ,lwdn ,ur , & !in uu ,vv ,sfctmp ,thair ,qair , & !in eair ,rhoair ,snowh ,vai ,gammav ,gammag, & !in - laimax, saimax,fwet ,laisun ,laisha ,cwp ,dzsnso , & !in + fwet ,laisun ,laisha ,cwp ,dzsnso , & !in zlvl ,zpd ,z0m ,fveg ,shdfac, & !in z0mg ,emv ,emg ,canliq ,fsno, & !in canice ,stc ,df ,rssun ,rssha , & !in @@ -3732,8 +3710,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: fsno !snow fraction real (kind=kind_phys), intent(in) :: snowh !actual snow depth [m] - real (kind=kind_phys), intent(in) :: saimax !< monthly maximum stem area index, one-sided - real (kind=kind_phys), intent(in) :: laimax !< monthly maximum leaf area index, one-sided real (kind=kind_phys), intent(in) :: fwet !wetted fraction of canopy real (kind=kind_phys), intent(in) :: cwp !canopy wind parameter @@ -4060,7 +4036,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,laimax,saimax,vegtyp, & !in + zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,vegtyp, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf1,cm ,ch ) !out @@ -4520,7 +4496,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,4.5,1.4,ivgtyp , & !in + zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out @@ -5123,7 +5099,7 @@ end subroutine sfcdif2 !! compute surface drag coefficient cm for momentum and ch for heat. subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,fveg ,garea1 ,vegetated,vaie,laimax,saimax,vegtyp , & !in + zpd ,snowh ,fveg ,garea1 ,vegetated,vaie,vegtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf ,cm ,ch ) !out @@ -5154,8 +5130,6 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: garea1 ! grid area [km2] logical, intent(in ) :: vegetated ! .true. if vegetated real (kind=kind_phys), intent(in ) :: vaie ! vegetation area index [m2/m2] - real (kind=kind_phys), intent(in ) :: saimax !< monthly maximum stem area index, one-sided - real (kind=kind_phys), intent(in ) :: laimax !< monthly maximum leaf area index, one-sided integer , intent(in ) :: vegtyp ! vegetation category real (kind=kind_phys), intent(inout) :: ustarx ! friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm ! momentum stability correction, weighted by prior iters @@ -5183,10 +5157,14 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys) :: sigmaa ! momentum partition parameter real (kind=kind_phys) :: tem1,tem2,zvfun1,gdx real (kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 + real (kind=kind_phys) :: saimax !< monthly maximum stem area index, one-sided + real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided ! ------------------------------------------------------------------------------------------------- fv = ustarx + laimax = maxval(parameters%laim) + saimax = maxval(parameters%saim) ! fv = ur*vkc/log((zlvl-zpd)/z0m) if(vegetated) then From 70507a0cdca6274e23943f9e96ac06750d7bf410 Mon Sep 17 00:00:00 2001 From: helin wei Date: Thu, 10 Mar 2022 22:23:34 +0000 Subject: [PATCH 16/21] to avoid exception floating point --- physics/module_sf_noahmplsm.f90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index c945e66ff..99b0cde7f 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2507,8 +2507,12 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) - laimax = maxval(parameters%laim) - df(1) = df(1) * exp (sbeta * elai/laimax) + if(elai.gt.0.) then + laimax = maxval(parameters%laim) + laimax = min(laimax, 0.1) + + df(1) = df(1) * exp (sbeta * elai/laimax) + endif ! compute lake thermal properties ! (no consideration of turbulent mixing for this version) @@ -5165,6 +5169,9 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur fv = ustarx laimax = maxval(parameters%laim) saimax = maxval(parameters%saim) + laimax = min(laimax, 0.1) + saimax = min(saimax, 0.1) + ! fv = ur*vkc/log((zlvl-zpd)/z0m) if(vegetated) then From d33598bb0079f56aeec3af97689fb24cb04049ba Mon Sep 17 00:00:00 2001 From: helin wei Date: Fri, 11 Mar 2022 16:28:53 +0000 Subject: [PATCH 17/21] revert the df1 change due to some negative impact on surface temperature --- physics/module_sf_noahmplsm.f90 | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 99b0cde7f..f9024c321 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2041,7 +2041,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in call thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , elai, & !in + lat ,z0m ,zlvl ,vegtyp , & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out @@ -2431,7 +2431,7 @@ end subroutine energy subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , & !in dt ,snowh ,snice ,snliq , & !in smc ,sh2o ,tg ,stc ,ur , & !in - lat ,z0m ,zlvl ,vegtyp , elai, & !in + lat ,z0m ,zlvl ,vegtyp , & !in df ,hcpct ,snicev ,snliqv ,epore , & !out fact ) !out ! ------------------------------------------------------------------------------------------------- @@ -2456,7 +2456,6 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), intent(in) :: lat !latitude (radians) real (kind=kind_phys), intent(in) :: z0m !roughness length (m) real (kind=kind_phys), intent(in) :: zlvl !reference height (m) - real (kind=kind_phys), intent(in) :: elai !lai adjusted for burying by snow integer , intent(in) :: vegtyp !vegtyp type ! outputs @@ -2474,7 +2473,6 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , real (kind=kind_phys), dimension(-nsnow+1: 0) :: tksno !snow thermal conductivity (j/m3/k) real (kind=kind_phys), dimension( 1:nsoil) :: sice !soil ice content real (kind=kind_phys), parameter :: sbeta = -2.0 - real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided ! -------------------------------------------------------------------------------------------------- ! compute snow thermal conductivity and heat capacity @@ -2507,12 +2505,6 @@ subroutine thermoprop (parameters,nsoil ,nsnow ,isnow ,ist ,dzsnso , ! not in use because of the separation of the canopy layer from the ground. ! but this may represent the effects of leaf litter (niu comments) ! df1 = df1 * exp (sbeta * shdfac) - if(elai.gt.0.) then - laimax = maxval(parameters%laim) - laimax = min(laimax, 0.1) - - df(1) = df(1) * exp (sbeta * elai/laimax) - endif ! compute lake thermal properties ! (no consideration of turbulent mixing for this version) From 3095d719239fbc804d632eeca711e7d5ed2680fd Mon Sep 17 00:00:00 2001 From: helin wei Date: Mon, 14 Mar 2022 21:06:18 +0000 Subject: [PATCH 18/21] correct the condition to avoid a divide by zero exception --- physics/module_sf_noahmplsm.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index f9024c321..1460e61f4 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -5151,7 +5151,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys) :: czil1 ! canopy based czil real (kind=kind_phys) :: fm10 ! 10-m stability adjustment - stability output real (kind=kind_phys) :: sigmaa ! momentum partition parameter - real (kind=kind_phys) :: tem1,tem2,zvfun1,gdx + real (kind=kind_phys) :: tem1,tem2,zvfun1,gdx,slaifrac real (kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 real (kind=kind_phys) :: saimax !< monthly maximum stem area index, one-sided real (kind=kind_phys) :: laimax !< monthly maximum leaf area index, one-sided @@ -5161,8 +5161,12 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur fv = ustarx laimax = maxval(parameters%laim) saimax = maxval(parameters%saim) - laimax = min(laimax, 0.1) - saimax = min(saimax, 0.1) + + if(laimax+saimax .gt. 0) then + slaifrac=vaie/(laimax+saimax) + else + slaifrac=0.1_kind_phys + endif ! fv = ur*vkc/log((zlvl-zpd)/z0m) @@ -5214,7 +5218,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur tem1 = (z0m - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys) - tem2 = max(vaie/(laimax+saimax), 0.1_kind_phys) + tem2 = max(slaifrac, 0.1_kind_phys) zvfun1 = sqrt(tem1 * tem2) gdx = sqrt(garea1) From 27ea849d8e88c70f6e3a1d014a0c85d0dd6ef2b9 Mon Sep 17 00:00:00 2001 From: helin wei Date: Tue, 15 Mar 2022 13:25:51 +0000 Subject: [PATCH 19/21] further refinement of the impact of vegetation on zvfun --- physics/module_sf_noahmplsm.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 1460e61f4..360536ec3 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -5164,6 +5164,8 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur if(laimax+saimax .gt. 0) then slaifrac=vaie/(laimax+saimax) + slaifrac=min(slaifrac,1.) + slaifrac=fveg*slaifrac else slaifrac=0.1_kind_phys endif From c722905e5240250ac7986624af0688f12737d8fb Mon Sep 17 00:00:00 2001 From: helin wei Date: Wed, 16 Mar 2022 03:20:38 +0000 Subject: [PATCH 20/21] replace shdfac by fveg for zvfun --- physics/module_sf_noahmplsm.f90 | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 360536ec3..ef022b4ee 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -4032,7 +4032,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh ,shdfac ,garea1 ,.true. ,vaie ,vegtyp, & !in + zpd ,snowh ,fveg ,garea1 ,.true. ,vaie ,vegtyp, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf1,cm ,ch ) !out @@ -4492,7 +4492,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,ivgtyp , & !in + zpd ,snowh,fveg ,garea1 ,.false. ,0.0,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out @@ -5161,14 +5161,17 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur fv = ustarx laimax = maxval(parameters%laim) saimax = maxval(parameters%saim) - - if(laimax+saimax .gt. 0) then + if(dveg.eq.4 .or. dveg.eq.5) then + if(laimax+saimax .gt. 0 .and. fveg .gt. 0) then slaifrac=vaie/(laimax+saimax) slaifrac=min(slaifrac,1.) slaifrac=fveg*slaifrac else slaifrac=0.1_kind_phys endif + else + slaifrac=fveg + endif ! fv = ur*vkc/log((zlvl-zpd)/z0m) From 4284846e2110f9c6e6781de2c002ae35964f03a1 Mon Sep 17 00:00:00 2001 From: weizhong zheng Date: Fri, 18 Mar 2022 14:45:38 +0000 Subject: [PATCH 21/21] modify the eddy diffusivity for heat at the top of the canopy --- physics/module_sf_noahmplsm.f90 | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index ef022b4ee..6e59407bb 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -3828,6 +3828,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: fm !momentum stability correction, weighted by prior iters real (kind=kind_phys) :: fh !sen heat stability correction, weighted by prior iters real (kind=kind_phys) :: fhg !sen heat stability correction, ground + real (kind=kind_phys) :: fhgh !sen heat stability correction, canopy real (kind=kind_phys) :: hcan !canopy height (m) [note: hcan >= z0mg] real (kind=kind_phys) :: a !temporary calculation @@ -4048,7 +4049,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & call ragrb(parameters,iter ,vaie ,rhoair ,hg ,tah , & !in zpd ,z0mg ,z0hg ,hcan ,uc , & !in z0h ,fv ,cwp ,vegtyp ,mpe , & !in - tv ,mozg ,fhg ,iloc ,jloc , & !inout + tv ,mozg ,fhg ,fhgh ,iloc ,jloc , & !inout ramg ,rahg ,rawg ,rb ) !out ! es and d(es)/dt evaluated at tv @@ -4604,7 +4605,7 @@ end subroutine bare_flux subroutine ragrb(parameters,iter ,vai ,rhoair ,hg ,tah , & !in zpd ,z0mg ,z0hg ,hcan ,uc , & !in z0h ,fv ,cwp ,vegtyp ,mpe , & !in - tv ,mozg ,fhg ,iloc ,jloc , & !inout + tv ,mozg ,fhg ,fhgh ,iloc ,jloc , & !inout ramg ,rahg ,rawg ,rb ) !out ! -------------------------------------------------------------------------------------------------- ! compute under-canopy aerodynamic resistance rag and leaf boundary layer @@ -4638,6 +4639,7 @@ subroutine ragrb(parameters,iter ,vai ,rhoair ,hg ,tah , & !in real (kind=kind_phys), intent(inout) :: mozg !monin-obukhov stability parameter real (kind=kind_phys), intent(inout) :: fhg !stability correction + real (kind=kind_phys), intent(inout) :: fhgh !stability correction, canopy ! outputs real (kind=kind_phys) :: ramg !aerodynamic resistance for momentum (s/m) @@ -4652,29 +4654,36 @@ subroutine ragrb(parameters,iter ,vai ,rhoair ,hg ,tah , & !in real (kind=kind_phys) :: tmprah2 !temporary calculation for aerodynamic resistances real (kind=kind_phys) :: tmprb !temporary calculation for rb real (kind=kind_phys) :: molg,fhgnew,cwpc + real (kind=kind_phys) :: mozgh, fhgnewh ! -------------------------------------------------------------------------------------------------- ! stability correction to below canopy resistance mozg = 0. molg = 0. + mozgh = 0. if(iter > 1) then tmp1 = vkc * (grav/tah) * hg/(rhoair*cpair) if (abs(tmp1) .le. mpe) tmp1 = mpe molg = -1. * fv**3 / tmp1 mozg = min( (zpd-z0mg)/molg, 1.) + mozgh = min( (hcan - zpd)/molg, 1.) end if if (mozg < 0.) then fhgnew = (1. - 15.*mozg)**(-0.25) + fhgnewh = 0.74 * (1. - 9.*mozg)**(-0.5) ! PHIh else fhgnew = 1.+ 4.7*mozg + fhgnewh = 0.74 + 4.7*mozgh ! PHIh endif if (iter == 1) then fhg = fhgnew + fhgh = fhgnewh else fhg = 0.5 * (fhg+fhgnew) + fhgh = 0.5 * (fhgh+fhgnewh) endif cwpc = (cwp * vai * hcan * fhg)**0.5 @@ -4686,7 +4695,7 @@ subroutine ragrb(parameters,iter ,vai ,rhoair ,hg ,tah , & !in ! aerodynamic resistances raw and rah between heights zpd+z0h and z0hg. - kh = max ( vkc*fv*(hcan-zpd), mpe ) + kh = max ( vkc*fv*(hcan-zpd)/fhgh, mpe ) ramg = 0. rahg = tmprah2 / kh rawg = rahg