From 3307baa8827bad76efa7597f910e835110015e7d Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 15:19:14 +0000 Subject: [PATCH 01/52] updated hfx2 & qfx2 --- physics/cu_gf_driver.meta | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/cu_gf_driver.meta b/physics/cu_gf_driver.meta index 73ce19754..6d1315381 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -286,8 +286,8 @@ intent = in optional = F [hfx2] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -295,8 +295,8 @@ intent = in optional = F [qfx2] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_latent_heat_flux + long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) type = real From 89cd97f87881a614f8d5d71162a66bbbf08667e2 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 15:24:46 +0000 Subject: [PATCH 02/52] updated evap & hfx --- physics/cu_ntiedtke.meta | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/cu_ntiedtke.meta b/physics/cu_ntiedtke.meta index 4d4c6597a..235168f83 100644 --- a/physics/cu_ntiedtke.meta +++ b/physics/cu_ntiedtke.meta @@ -204,8 +204,8 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_latent_heat_flux + long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) type = real @@ -213,8 +213,8 @@ intent = in optional = F [hfx] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real From ccac8315af699a31f7462be7606baac1e24a5648 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 16:34:40 +0000 Subject: [PATCH 03/52] updated hflx & evap --- physics/gcm_shoc.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/gcm_shoc.meta b/physics/gcm_shoc.meta index b021fa306..f44560890 100644 --- a/physics/gcm_shoc.meta +++ b/physics/gcm_shoc.meta @@ -288,8 +288,8 @@ intent = in optional = F [hflx] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -297,7 +297,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From 51f2842b98ee6b9cbd06bdcbede1bcf4058cad4a Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 16:35:41 +0000 Subject: [PATCH 04/52] updated GFS debug --- physics/GFS_debug.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 99f36f077..57bc05ede 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -1091,7 +1091,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_ice ', Interstitial%ep1d_ice ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_land ', Interstitial%ep1d_land ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_water ', Interstitial%ep1d_water ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evapq ', Interstitial%evapq ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_ice ', Interstitial%evap_ice ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_land ', Interstitial%evap_land ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_water ', Interstitial%evap_water ) @@ -1134,7 +1133,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gflx_water ', Interstitial%gflx_water ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gwdcu ', Interstitial%gwdcu ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gwdcv ', Interstitial%gwdcv ) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hefac ', Interstitial%hefac ) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%zvfun ', Interstitial%zvfun ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hffac ', Interstitial%hffac ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hflxq ', Interstitial%hflxq ) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hflx_ice ', Interstitial%hflx_ice ) From 662bb1f4821ebc5e389ca371af126e9349008081 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 16:40:55 +0000 Subject: [PATCH 05/52] updated GFS PBL generic --- physics/GFS_PBL_generic.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index 5183836ec..403a68b35 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -317,7 +317,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqsfc_cpl, dusfci_cpl, dvsfci_cpl, dtsfci_cpl, dqsfci_cpl, dusfc_diag, dvsfc_diag, dtsfc_diag, dqsfc_diag, & dusfci_diag, dvsfci_diag, dtsfci_diag, dqsfci_diag, dt3dt, du3dt_PBL, du3dt_OGWD, dv3dt_PBL, dv3dt_OGWD, dq3dt, & dq3dt_ozone, rd, cp, fvirt, hvap, t1, q1, prsl, hflx, ushfsfci, oceanfrac, kdt, dusfc_cice, dvsfc_cice, & - dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, dkt_cpl, dkt, hffac, hefac, & + dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, dkt_cpl, dkt, hffac, & ugrs, vgrs, tgrs, qgrs, save_u, save_v, save_t, save_q, errmsg, errflg) use machine, only : kind_phys @@ -370,7 +370,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, real(kind=kind_phys), dimension(:,:), intent(in) :: dkt ! From canopy heat storage - reduction factors in latent/sensible heat flux due to surface roughness - real(kind=kind_phys), dimension(:), intent(in) :: hffac, hefac + real(kind=kind_phys), dimension(:), intent(in) :: hffac character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -556,7 +556,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, else !use PBL fluxes when CICE fluxes is unavailable dusfci_cpl(i) = dusfc1(i) dvsfci_cpl(i) = dvsfc1(i) - dtsfci_cpl(i) = dtsfc1(i) + dtsfci_cpl(i) = dtsfc1(i)*hffac(i) dqsfci_cpl(i) = dqsfc1(i) end if elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point @@ -575,7 +575,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dusfci_cpl(i) = dusfc1(i) dvsfci_cpl(i) = dvsfc1(i) dtsfci_cpl(i) = dtsfc1(i)*hffac(i) - dqsfci_cpl(i) = dqsfc1(i)*hefac(i) + dqsfci_cpl(i) = dqsfc1(i) endif ! dusfc_cpl (i) = dusfc_cpl(i) + dusfci_cpl(i) * dtf @@ -601,7 +601,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dusfci_diag(i) = dusfc1(i) dvsfci_diag(i) = dvsfc1(i) dtsfci_diag(i) = dtsfc1(i)*hffac(i) - dqsfci_diag(i) = dqsfc1(i)*hefac(i) + dqsfci_diag(i) = dqsfc1(i) dtsfc_diag (i) = dtsfc_diag(i) + dtsfci_diag(i) * dtf dqsfc_diag (i) = dqsfc_diag(i) + dqsfci_diag(i) * dtf enddo From 0dcde4af5652493772e9f619a926566db3ee0457 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 16:45:25 +0000 Subject: [PATCH 06/52] updated GFS PBL generic meta --- physics/GFS_PBL_generic.meta | 9 --------- 1 file changed, 9 deletions(-) diff --git a/physics/GFS_PBL_generic.meta b/physics/GFS_PBL_generic.meta index 01f5ec3ce..9ca84695c 100644 --- a/physics/GFS_PBL_generic.meta +++ b/physics/GFS_PBL_generic.meta @@ -1343,15 +1343,6 @@ kind = kind_phys intent = in optional = F -[hefac] - standard_name = surface_upward_latent_heat_flux_reduction_factor - long_name = surface upward latent heat flux reduction factor from canopy heat storage - units = none - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [ugrs] standard_name = x_wind long_name = zonal wind From 989701b0a2d4e7ca273bb7b09adc113d63715bf4 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 16:52:11 +0000 Subject: [PATCH 07/52] updated GFS surface generic --- physics/GFS_surface_generic.F90 | 63 +++++++++++++++++++-------------- 1 file changed, 37 insertions(+), 26 deletions(-) diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 70a5b3541..5c52c6c04 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -216,12 +216,13 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl, & dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl, & v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl, & - nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, & - runoff, srunoff, runof, drain, lheatstrg, z0fac, e0fac, zorl, hflx, evap, hflxq, evapq, hffac, hefac, errmsg, errflg) + nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, islmsk, sigmaf, & + runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zorl, hflx, evap, hflxq, zvfun, hffac, errmsg, errflg) implicit none integer, intent(in) :: im + integer, dimension(im), intent(in) :: islmsk logical, intent(in) :: cplflx, cplwav, lssav logical, dimension(:), intent(in) :: icy, wet real(kind=kind_phys), intent(in) :: dtf @@ -237,15 +238,15 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt evcwa, transa, sbsnoa, snowca, snohfa, ep real(kind=kind_phys), dimension(:), intent(inout) :: runoff, srunoff - real(kind=kind_phys), dimension(:), intent(in) :: drain, runof + real(kind=kind_phys), dimension(:), intent(in) :: drain, runof, sigmaf ! For canopy heat storage logical, intent(in) :: lheatstrg - real(kind=kind_phys), intent(in) :: z0fac, e0fac + real(kind=kind_phys), intent(in) :: h0facu, h0facs real(kind=kind_phys), dimension(:), intent(in) :: zorl real(kind=kind_phys), dimension(:), intent(in) :: hflx, evap - real(kind=kind_phys), dimension(:), intent(out) :: hflxq, evapq - real(kind=kind_phys), dimension(:), intent(out) :: hffac, hefac + real(kind=kind_phys), dimension(:), intent(out) :: hflxq + real(kind=kind_phys), dimension(:), intent(out) :: zvfun, hffac ! CCPP error handling variables character(len=*), intent(out) :: errmsg @@ -255,8 +256,7 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt real(kind=kind_phys), parameter :: albdf = 0.06_kind_phys ! Parameters for canopy heat storage parametrization - real(kind=kind_phys), parameter :: z0min=0.2, z0max=1.0 - real(kind=kind_phys), parameter :: u10min=2.5, u10max=7.5 + real(kind=kind_phys), parameter :: z0min=0.1, z0max=1.0 integer :: i real(kind=kind_phys) :: xcosz_loc, ocalnirdf_cpl, ocalnirbm_cpl, ocalvisdf_cpl, ocalvisbm_cpl @@ -361,32 +361,43 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt enddo endif -! --- ... Boundary Layer and Free atmospheic turbulence parameterization ! -! in order to achieve heat storage within canopy layer, in the canopy heat -! storage parameterization the kinematic sensible and latent heat fluxes -! (hflx & evap) as surface boundary forcings to the pbl scheme are -! reduced as a function of surface roughness +! in order to achieve heat storage within canopy layer, in the canopy +! heat torage parameterization the kinematic sensible heat flux +! (hflx) as surface boundary forcing to the pbl scheme is +! reduced as a function of surface roughness & green vegetation +! fraction ! +! background diffusivity & background mixing length are also given by +! a function of surface roughness & green vegetation fraction +! + do i=1,im + if(islmsk(i) == 1) then + tem = 0.01 * zorl(i) ! change unit from cm to m + tem1 = (tem - z0min) / (z0max - z0min) + tem1 = min(max(tem1, 0.0), 1.0) + tem2 = max(sigmaf(i), 0.1) +! tem2 = sigmaf(i) + zvfun(i) = sqrt(tem1 * tem2) + else + zvfun(i) = 0. + endif + enddo do i=1,im hflxq(i) = hflx(i) - evapq(i) = evap(i) hffac(i) = 1.0 - hefac(i) = 1.0 enddo if (lheatstrg) then do i=1,im - tem = 0.01 * zorl(i) ! change unit from cm to m - tem1 = (tem - z0min) / (z0max - z0min) - hffac(i) = z0fac * min(max(tem1, 0.0), 1.0) - tem = sqrt(u10m(i)**2+v10m(i)**2) - tem1 = (tem - u10min) / (u10max - u10min) - tem2 = 1.0 - min(max(tem1, 0.0), 1.0) - hffac(i) = tem2 * hffac(i) - hefac(i) = 1. + e0fac * hffac(i) - hffac(i) = 1. + hffac(i) - hflxq(i) = hflx(i) / hffac(i) - evapq(i) = evap(i) / hefac(i) + if(islmsk(i) == 1) then + if(hflx(i) > 0.) then + hffac(i) = h0facu * zvfun(i) + else + hffac(i) = h0facs * zvfun(i) + endif + hffac(i) = 1. + hffac(i) + hflxq(i) = hflx(i) / hffac(i) + endif enddo endif From b542460cda828d83f4c6544b6a93c8523299fec8 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:03:52 +0000 Subject: [PATCH 08/52] updated GFS surface generic meta --- physics/GFS_surface_generic.meta | 50 ++++++++++++++++++-------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index e174de153..963fde747 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -1222,6 +1222,23 @@ kind = kind_phys intent = inout optional = F +[islmsk] + standard_name = sea_land_ice_mask + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F +[sigmaf] + standard_name = bounded_vegetation_area_fraction + long_name = areal fractional cover of green vegetation bounded on the bottom + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [runoff] standard_name = total_runoff long_name = total water runoff @@ -1266,18 +1283,18 @@ type = logical intent = in optional = F -[z0fac] - standard_name = surface_roughness_fraction_factor - long_name = surface roughness fraction factor for canopy heat storage parameterization +[h0facu] + standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_unstable_surface_layer + long_name = canopy heat storage factor for sensible heat flux in unstable surface layer units = none dimensions = () type = real kind = kind_phys intent = in optional = F -[e0fac] - standard_name = latent_heat_flux_fraction_factor_relative_to_sensible_heat_flux - long_name = latent heat flux fraction factor relative to sensible heat flux for canopy heat storage parameterization +[h0facs] + standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_stable_surface_layer + long_name = canopy heat storage factor for sensible heat flux in stable surface layer units = none dimensions = () type = real @@ -1312,18 +1329,18 @@ intent = in optional = F [hflxq] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = out optional = F -[evapq] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness - units = kg kg-1 m s-1 +[zvfun] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction + long_name = function of surface roughness length and green vegetation fraction + units = none dimensions = (horizontal_loop_extent) type = real kind = kind_phys @@ -1338,15 +1355,6 @@ kind = kind_phys intent = out optional = F -[hefac] - standard_name = surface_upward_latent_heat_flux_reduction_factor - long_name = surface upward latent heat flux reduction factor from canopy heat storage - units = none - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out - optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From b8a3a7ebdc6d9a39fe6837571a40b9891b3cc882 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:05:36 +0000 Subject: [PATCH 09/52] updated mfpblq --- physics/mfpbltq.f | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index a6fc22cef..b906052cd 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -40,12 +40,13 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, c local variables and arrays ! integer i, j, k, n, ndc + integer kpblx(im), kpbly(im) ! real(kind=kind_phys) dt2, dz, ce0, cm, & factor, gocp, & g, b1, f1, & bb1, bb2, - & a1, pgcon, + & alp, vpertmax,a1, pgcon, & qmin, qlmin, xmmx, rbint, & tem, tem1, tem2, & ptem, ptem1, ptem2 @@ -54,7 +55,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, & tlu, gamma, qlu, & thup, thvu, dq ! - real(kind=kind_phys) rbdn(im), rbup(im), xlamuem(im,km-1) + real(kind=kind_phys) rbdn(im), rbup(im), hpblx(im), + & xlamuem(im,km-1) real(kind=kind_phys) delz(im), xlamax(im) ! real(kind=kind_phys) wu2(im,km), thlu(im,km), @@ -71,7 +73,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) parameter(ce0=0.4,cm=1.0) parameter(qmin=1.e-8,qlmin=1.e-12) - parameter(pgcon=0.55) + parameter(alp=1.5,vpertmax=3.0,pgcon=0.55) parameter(b1=0.5,f1=0.15) ! !************************************************************************ @@ -99,9 +101,11 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! do i=1,im if(cnvflg(i)) then - thlu(i,1)= thlx(i,1) + vpert(i) + ptem = alp * vpert(i) + ptem = min(ptem, vpertmax) + thlu(i,1)= thlx(i,1) + ptem qtu(i,1) = qtx(i,1) - buo(i,1) = g * vpert(i) / thvx(i,1) + buo(i,1) = g * ptem / thvx(i,1) endif enddo ! @@ -213,6 +217,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! do i=1,im flg(i) = .true. + kpblx(i) = 1 + kpbly(i) = kpbl(i) if(cnvflg(i)) then flg(i) = .false. rbup(i) = wu2(i,1) @@ -223,14 +229,14 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, if(.not.flg(i)) then rbdn(i) = rbup(i) rbup(i) = wu2(i,k) - kpbl(i)= k + kpblx(i)= k flg(i) = rbup(i).le.0. endif enddo enddo do i = 1,im if(cnvflg(i)) then - k = kpbl(i) + k = kpblx(i) if(rbdn(i) <= 0.) then rbint = 0. elseif(rbup(i) >= 0.) then @@ -238,7 +244,17 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, else rbint = rbdn(i)/(rbdn(i)-rbup(i)) endif - hpbl(i) = zm(i,k-1) + rbint*(zm(i,k)-zm(i,k-1)) + hpblx(i) = zm(i,k-1) + rbint*(zm(i,k)-zm(i,k-1)) + endif + enddo +! + do i = 1,im + if(cnvflg(i)) then + if(kpblx(i) < kpbl(i)) then + kpbl(i) = kpblx(i) + hpbl(i) = hpblx(i) + endif + if(kpbl(i) <= 1) cnvflg(i)=.false. endif enddo ! @@ -255,7 +271,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! do k = 1, kmpbl do i=1,im - if(cnvflg(i)) then + if(cnvflg(i) .and. kpblx(i) < kpbly(i)) then +! if(cnvflg(i)) then if(k < kpbl(i)) then ptem = 1./(zm(i,k)+delz(i)) tem = max((hpbl(i)-zm(i,k)+delz(i)) ,delz(i)) From 4d301b16e6825393b3f0bd703248d43c67c0437a Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:11:33 +0000 Subject: [PATCH 10/52] updated evap & hflx --- physics/module_MYJPBL_wrapper.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_MYJPBL_wrapper.meta b/physics/module_MYJPBL_wrapper.meta index 9d70397e7..a064cbd85 100644 --- a/physics/module_MYJPBL_wrapper.meta +++ b/physics/module_MYJPBL_wrapper.meta @@ -474,7 +474,7 @@ intent = inout optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) @@ -483,8 +483,8 @@ intent = in optional = F [hflx] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real From 3bdbf1cfcfc3c3d71ca12546a12fe0808d5d4c8c Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:15:27 +0000 Subject: [PATCH 11/52] updated hflx & qflx --- physics/module_MYNNPBL_wrapper.meta | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/module_MYNNPBL_wrapper.meta b/physics/module_MYNNPBL_wrapper.meta index de24fcbef..6d87959f9 100644 --- a/physics/module_MYNNPBL_wrapper.meta +++ b/physics/module_MYNNPBL_wrapper.meta @@ -343,8 +343,8 @@ intent = out optional = F [hflx] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -352,8 +352,8 @@ intent = in optional = F [qflx] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_latent_heat_flux + long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) type = real From a034b146aab50861ff07eb6e9f2ef71fa502a845 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:22:10 +0000 Subject: [PATCH 12/52] updated heat & evap --- physics/moninedmf.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/moninedmf.meta b/physics/moninedmf.meta index 24096cbe6..9862efe9f 100644 --- a/physics/moninedmf.meta +++ b/physics/moninedmf.meta @@ -250,8 +250,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -259,7 +259,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From b567312333a6b0f32589ecce4075571622a96b20 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:25:51 +0000 Subject: [PATCH 13/52] updated heat & evap --- physics/moninshoc.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/moninshoc.meta b/physics/moninshoc.meta index 51f2c4536..18c047331 100644 --- a/physics/moninshoc.meta +++ b/physics/moninshoc.meta @@ -248,8 +248,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -257,7 +257,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From 26838633e1e516323a8ec1cf408534317da6e3bc Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:32:14 +0000 Subject: [PATCH 14/52] updated deep convection --- physics/samfdeepcnv.f | 389 +++++++++++++++++++++++++++--------------- 1 file changed, 250 insertions(+), 139 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 425aa92a9..a8905696c 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -10,9 +10,10 @@ module samfdeepcnv contains - subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, & + subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, +& & errmsg, errflg) - + integer, intent(in) :: imfdeepcnv integer, intent(in) :: imfdeepcnv_samf character(len=*), intent(out) :: errmsg @@ -21,7 +22,8 @@ subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, & ! Consistency checks if (imfdeepcnv/=imfdeepcnv_samf) then - write(errmsg,'(*(a))') 'Logic error: namelist choice of', & + write(errmsg,'(*(a))') 'Logic error: namelist choice of', +& & ' deep convection is different from SAMF scheme' errflg = 1 return @@ -80,10 +82,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & t0c,delt,ntk,ntr,delp, & & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav,hwrf_samfdeep, & & cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, & - & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & + & dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& - & clam,c0s,c1,betal,betas,evfact,evfactl,pgcon,asolfac, & + & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & & rainevap, & & errmsg,errflg) @@ -99,7 +101,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & fv, grav, hvap, rd, rv, t0c real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & - & prslp(:,:), garea(:), dot(:,:), phil(:,:) + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: hwrf_samfdeep real(kind=kind_phys), intent(in) :: nthresh @@ -108,6 +110,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger integer, intent(inout) :: kcnv(:) + ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH real(kind=kind_phys), intent(inout) :: qtr(:,:,:), & & q1(:,:), t1(:,:), u1(:,:), v1(:,:), & & cnvw(:,:), cnvc(:,:) @@ -128,7 +131,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: clam, c0s, c1, & & betal, betas, asolfac, & - & evfact, evfactl, pgcon + & evef, pgcon character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! @@ -137,12 +140,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! integer latd,lond ! real(kind=kind_phys) clamd, tkemx, tkemn, dtke, - & beta, dbeta, betamx, betamn, + & beta, clamca, & cxlame, cxlamd, - & cxlamu, & xlamde, xlamdd, - & crtlamd, - & crtlame + & crtlame, crtlamd ! ! real(kind=kind_phys) detad real(kind=kind_phys) adw, aup, aafac, d0, @@ -157,7 +158,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & edtmaxl, edtmaxs, el2orc, elocp, & es, etah, & cthk, dthk, - & evef, fact1, fact2, factor, +! & evfact, evfactl, + & fact1, fact2, factor, & gamma, pprime, cm, & qlk, qrch, qs, & rain, rfact, shear, tfac, @@ -171,15 +173,16 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & xqrch, tem, tem1, tem2, & ptem, ptem1, ptem2 ! - integer kb(im), kbcon(im), kbcon1(im), + integer kb(im), kb1(im), kbcon(im), kbcon1(im), & ktcon(im), ktcon1(im), ktconn(im), & jmin(im), lmin(im), kbmax(im), - & kbm(im), kmax(im) + & kbm(im), kmax(im), kd94(im) ! ! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), real(kind=kind_phys) aa1(im), tkemean(im),clamt(im), & ps(im), del(im,km), prsl(im,km), - & umean(im), tauadv(im), gdx(im), +! & umean(im), tauadv(im), gdx(im), + & gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), & deltv(im), dtconv(im), edt(im), @@ -197,10 +200,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! & xpwev(im), delebar(im,ntr), & delubar(im), delvbar(im) ! - real(kind=kind_phys) c0(im) + real(kind=kind_phys) c0(im), sfcpbl(im) cj real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, - & cinacr, cinacrmx, cinacrmn + & cinacr, cinacrmx, cinacrmn, + & sfclfac, rhcrt cj ! ! parameters for updraft velocity calculation @@ -226,9 +230,9 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & parameter(cm=1.0) ! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) parameter(clamd=0.03,tkemx=0.65,tkemn=0.05) + parameter(clamca=0.03) parameter(dtke=tkemx-tkemn) - parameter(dbeta=0.1) - parameter(cthk=150.,dthk=25.) + parameter(cthk=200.,dthk=25.,sfclfac=0.2,rhcrt=0.75) parameter(cinpcrmx=180.,cinpcrmn=120.) ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) @@ -251,7 +255,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km), & dbyo(im,km), zo(im,km), & xlamue(im,km), xlamud(im,km), - & fent1(im,km), fent2(im,km), frh(im,km), + & fent1(im,km), fent2(im,km), + & rh(im,km), frh(im,km), & heo(im,km), heso(im,km), & qrcd(im,km), dellah(im,km), dellaq(im,km), & dellae(im,km,ntr), @@ -262,7 +267,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & qrcko(im,km), qrcdo(im,km), & pwo(im,km), pwdo(im,km), c0t(im,km), & tx1(im), sumx(im), cnvwt(im,km) -! &, rhbar(im) + &, rhbar(im) ! logical do_aerosols, totflg, cnvflg(im), asqecflg(im), flg(im) ! @@ -317,6 +322,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & c do i=1,im cnvflg(i) = .true. + sfcpbl(i) = sfclfac * hpbl(i) rn(i)=0. mbdt(i)=10. kbot(i)=km+1 @@ -424,16 +430,17 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! model tunable parameters are all here edtmaxl = .3 edtmaxs = .3 +! evfact = 0.3 +! evfactl = 0.3 + aafac = .1 if (hwrf_samfdeep) then - aafac = .1 - cxlamu = 1.0e-3 + cxlame = 1.0e-3 else - aafac = .05 cxlame = 1.0e-4 endif - crtlamd = 1.0e-4 + cxlamd = 0.75e-4 crtlame = 1.0e-4 - cxlamd = 1.0e-4 + crtlamd = 1.0e-4 xlamde = 1.0e-4 xlamdd = 1.0e-4 ! @@ -457,6 +464,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & kbmax(i) = km kbm(i) = km kmax(i) = km + kd94(i) = km tx1(i) = 1.0 / ps(i) enddo ! @@ -465,12 +473,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if (prsl(i,k)*tx1(i) > 0.04) kmax(i) = k + 1 if (prsl(i,k)*tx1(i) > 0.45) kbmax(i) = k + 1 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1 + if (prsl(i,k)*tx1(i) > 0.94) kd94(i) = k + 1 enddo enddo do i=1,im kmax(i) = min(km,kmax(i)) kbmax(i) = min(kbmax(i),kmax(i)) kbm(i) = min(kbm(i),kmax(i)) + kd94(i) = min(kd94(i),kmax(i)) enddo c c hydrostatic height assume zero terr and initially assume @@ -507,6 +517,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & eta(i,k) = 1. fent1(i,k)= 1. fent2(i,k)= 1. + rh(i,k) = 0. frh(i,k) = 0. hcko(i,k) = 0. qcko(i,k) = 0. @@ -591,14 +602,32 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & c this is the level where updraft starts c !> ## Perform calculations related to the updraft of the entraining/detraining cloud model ("static control"). -!> - Search below index "kbm" for the level of maximum moist static energy. +!> - Find the index for a level of sfclfac*hpbl which is initial guess for the parcel starting level. + do i=1,im + flg(i) = .true. + kb1(i) = 1 + enddo + do k = 2, km1 + do i=1,im + if (flg(i) .and. zo(i,k) <= sfcpbl(i)) then + kb1(i) = k + else + flg(i) = .false. + endif + enddo + enddo do i=1,im - hmax(i) = heo(i,1) - kb(i) = 1 + kb1(i) = min(kb1(i),kbm(i)) + enddo +c +!> - Search below index "kbm" and above kb1 for the level of maximum moist static energy. + do i=1,im + hmax(i) = heo(i,kb1(i)) + kb(i) = kb1(i) enddo do k = 2, km do i=1,im - if (k <= kbm(i)) then + if (k > kb1(i) .and. k <= kbm(i)) then if(heo(i,k) > hmax(i)) then kb(i) = k hmax(i) = heo(i,k) @@ -640,8 +669,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & val2 = 1.e-10 qo(i,k) = max(qo(i,k), val2 ) ! qo(i,k) = min(qo(i,k),qeso(i,k)) - tem = min(qo(i,k)/qeso(i,k), 1.) - frh(i,k) = 1. - tem + rh(i,k) = min(qo(i,k)/qeso(i,k), 1.) + frh(i,k) = 1. - rh(i,k) heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) + & cp * to(i,k) + hvap * qo(i,k) heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) + @@ -685,14 +714,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i=1,im if(kbcon(i) == kmax(i)) cnvflg(i) = .false. enddo -!! - if(do_ca .and. ca_trigger)then - do i=1,im - if(ca_deep(i) > nthresh) then - cnvflg(i) = .true. - endif - enddo - endif !! totflg = .true. do i=1,im @@ -746,13 +767,112 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo !! if(do_ca .and. ca_trigger)then + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo + endif + + totflg = .true. do i=1,im - if(ca_deep(i) > nthresh) then - cnvflg(i) = .true. - endif + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! +! +! re-define kb & kbcon +! + do i=1,im + if (cnvflg(i)) then + hmax(i) = heo(i,1) + kb(i) = 1 + endif + enddo + do k = 2, km + do i=1,im + if (cnvflg(i) .and. k <= kbm(i)) then + if(heo(i,k) > hmax(i)) then + kb(i) = k + hmax(i) = heo(i,k) + endif + endif + enddo + enddo +! + do i=1,im + flg(i) = cnvflg(i) + if(flg(i)) kbcon(i) = kmax(i) enddo + do k = 1, km1 + do i=1,im + if (flg(i) .and. k <= kbmax(i)) then + if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then + kbcon(i) = k + flg(i) = .false. + endif + endif + enddo + enddo +! + do i=1,im + if(cnvflg(i) .and. kbcon(i) == kmax(i)) then + cnvflg(i) = .false. + endif + enddo +!! + if(do_ca .and. ca_trigger)then + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo endif + + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! + do i=1,im + if(cnvflg(i)) then +! pdot(i) = 10.* dot(i,kbcon(i)) + pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s + endif + enddo +! +!> - if the mean relative humidity in the subcloud layers is less than a threshold value (rhcrt), convection is not triggered. +! + do i = 1, im + rhbar(i) = 0. + sumx(i) = 0. + enddo + do k = 1, km1 + do i = 1, im + if (cnvflg(i)) then + if(k >= kb(i) .and. k < kbcon(i)) then + dz = zo(i,k+1) - zo(i,k) + rhbar(i) = rhbar(i) + rh(i,k) * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo + do i= 1, im + if(cnvflg(i)) then + rhbar(i) = rhbar(i) / sumx(i) + if(rhbar(i) < rhcrt) then + cnvflg(i) = .false. + endif + endif + enddo !! + if(do_ca .and. ca_trigger)then + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo + endif + totflg = .true. do i=1,im totflg = totflg .and. (.not. cnvflg(i)) @@ -760,6 +880,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(totflg) return !! ! +!Lisa: at this point only trigger criteria is set + ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE ! @@ -799,13 +921,25 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo ! + if(do_ca .and. ca_entr)then + do i=1,im + if(cnvflg(i)) then + if(ca_deep(i) > nthresh)then + clamt(i) = clam - clamca + else + clamt(i) = clam + endif + endif + enddo + endif + else ! if(do_ca .and. ca_entr)then do i=1,im if(cnvflg(i)) then if(ca_deep(i) > nthresh)then - clamt(i) = clam - clamd + clamt(i) = clam - clamca else clamt(i) = clam endif @@ -827,7 +961,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im if(cnvflg(i)) then - xlamue(i,k) = clamt(i) / zi(i,k) + dz =zo(i,k+1) - zo(i,k) + xlamue(i,k) = clamt(i) / (zi(i,k) + dz) xlamue(i,k) = max(xlamue(i,k), crtlame) endif enddo @@ -874,6 +1009,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im if(cnvflg(i) .and. k < kmax(i)) then +! xlamud(i,k) = crtlamd xlamud(i,k) = 0.001 * clamt(i) endif enddo @@ -904,7 +1040,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tem = cxlamu * frh(i,k) * fent2(i,k) + tem = cxlame * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem endif enddo @@ -1071,14 +1207,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo !! + if(do_ca .and. ca_trigger)then - do i=1,im - if(ca_deep(i) > nthresh) then - cnvflg(i) = .true. - endif - enddo + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo endif -!! + totflg = .true. do i = 1, im totflg = totflg .and. (.not. cnvflg(i)) @@ -1154,13 +1290,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif !hwrf_samfdeep !! if(do_ca .and. ca_trigger)then - do i=1,im - if(ca_deep(i) > nthresh) then - cnvflg(i) = .true. - endif - enddo + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo endif -!! + totflg = .true. do i=1,im totflg = totflg .and. (.not. cnvflg(i)) @@ -1195,21 +1330,22 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(tem < cthk) cnvflg(i) = .false. endif enddo -!! + + if(do_ca .and. ca_trigger)then - do i=1,im - if(ca_deep(i) > nthresh) then - cnvflg(i) = .true. - endif - enddo + do i=1,im + if(ca_deep(i) > nthresh) cnvflg(i) = .true. + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + enddo endif -!! + totflg = .true. - do i = 1, im + do i=1,im totflg = totflg .and. (.not. cnvflg(i)) enddo if(totflg) return !! + c c search for downdraft originating level above theta-e minimum c @@ -1644,60 +1780,34 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i = 1, im if(cnvflg(i)) then - if(k >= 1 .and. k < kbcon(i)) then + if(k >= 1 .and. k < kd94(i)) then dz = zi(i,k+1) - zi(i,k) sumx(i) = sumx(i) + dz endif endif enddo enddo - - if (hwrf_samfdeep) then - do i = 1, im + do i = 1, im beta = betas if(islimsk(i) == 1) beta = betal if(cnvflg(i)) then - dz = (sumx(i)+zi(i,1))/float(kbcon(i)) - tem = 1./float(kbcon(i)) + dz = (sumx(i)+zi(i,1))/float(kd94(i)) + tem = 1./float(kd94(i)) xlamd(i) = (1.-beta**tem)/dz endif - enddo - else - do i = 1, im - if(cnvflg(i)) then - betamn = betas - if(islimsk(i) == 1) betamn = betal - if(ntk > 0) then - betamx = betamn + dbeta - if(tkemean(i) > tkemx) then - beta = betamn - else if(tkemean(i) < tkemn) then - beta = betamx - else - tem = (betamx - betamn) * (tkemean(i) - tkemn) - beta = betamx - tem / dtke - endif - else - beta = betamn - endif - dz = (sumx(i)+zi(i,1))/float(kbcon(i)) - tem = 1./float(kbcon(i)) - xlamd(i) = (1.-beta**tem)/dz - endif - enddo - endif + enddo c c determine downdraft mass flux c -!> - Calculate the normalized downdraft mass flux from equation 1 of Pan and Wu (1995) \cite pan_and_wu_1995 . Downdraft entrainment and detrainment rates are constants from the downdraft origination to the LFC. +!> - Calculate the normalized downdraft mass flux from equation 1 of Pan and Wu (1995) \cite pan_and_wu_1995 . Downdraft entrainment and detrainment rates are constants from the downdraft origination to the level of 60mb above the ground surface (kd94). do k = km1, 1, -1 do i = 1, im if (cnvflg(i) .and. k <= kmax(i)-1) then - if(k < jmin(i) .and. k >= kbcon(i)) then + if(k < jmin(i) .and. k >= kd94(i)) then dz = zi(i,k+1) - zi(i,k) ptem = xlamdd - xlamde etad(i,k) = etad(i,k+1) * (1. - ptem * dz) - else if(k < kbcon(i)) then + else if(k < kd94(i)) then dz = zi(i,k+1) - zi(i,k) ptem = xlamd(i) + xlamdd - xlamde etad(i,k) = etad(i,k+1) * (1. - ptem * dz) @@ -1737,7 +1847,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (cnvflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) - if(k >= kbcon(i)) then + if(k >= kd94(i)) then tem = xlamde * dz tem1 = 0.5 * xlamdd * dz else @@ -1786,7 +1896,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! detad = etad(i,k+1) - etad(i,k) cj dz = zi(i,k+1) - zi(i,k) - if(k >= kbcon(i)) then + if(k >= kd94(i)) then tem = xlamde * dz tem1 = 0.5 * xlamdd * dz else @@ -1935,7 +2045,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1)) c - if(k <= kbcon(i)) then + if(k <= kd94(i)) then ptem = xlamde ptem1 = xlamd(i)+xlamdd else @@ -2247,7 +2357,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (asqecflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) - if(k >= kbcon(i)) then + if(k >= kd94(i)) then tem = xlamde * dz tem1 = 0.5 * xlamdd * dz else @@ -2272,7 +2382,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! detad = etad(i,k+1) - etad(i,k) cj dz = zi(i,k+1) - zi(i,k) - if(k >= kbcon(i)) then + if(k >= kd94(i)) then tem = xlamde * dz tem1 = 0.5 * xlamdd * dz else @@ -2424,40 +2534,41 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. - do i= 1, im - if(cnvflg(i)) then - sumx(i) = 0. - umean(i) = 0. - endif - enddo - do k = 2, km1 - do i = 1, im - if(cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon1(i)) then - dz = zi(i,k) - zi(i,k-1) - tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) - umean(i) = umean(i) + tem * dz - sumx(i) = sumx(i) + dz - endif - endif - enddo - enddo - do i= 1, im - if(cnvflg(i)) then - umean(i) = umean(i) / sumx(i) - umean(i) = max(umean(i), 1.) - tauadv(i) = gdx(i) / umean(i) - endif - enddo +! do i= 1, im +! if(cnvflg(i)) then +! sumx(i) = 0. +! umean(i) = 0. +! endif +! enddo +! do k = 2, km1 +! do i = 1, im +! if(cnvflg(i)) then +! if(k >= kbcon1(i) .and. k < ktcon1(i)) then +! dz = zi(i,k) - zi(i,k-1) +! tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) +! umean(i) = umean(i) + tem * dz +! sumx(i) = sumx(i) + dz +! endif +! endif +! enddo +! enddo +! do i= 1, im +! if(cnvflg(i)) then +! umean(i) = umean(i) / sumx(i) +! umean(i) = max(umean(i), 1.) +! tauadv(i) = gdx(i) / umean(i) +! endif +! enddo !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - tfac = tauadv(i) / dtconv(i) - tfac = min(tfac, 1.) - xmb(i) = tfac*betaw*rho*wc(i) +! tfac = tauadv(i) / dtconv(i) +! tfac = min(tfac, 1.) +! xmb(i) = tfac*betaw*rho*wc(i) + xmb(i) = betaw*rho*wc(i) endif enddo !> - For the cases where the quasi-equilibrium assumption of Arakawa-Schubert is valid, first calculate the large scale destabilization as in equation 5 of Pan and Wu (1995) \cite pan_and_wu_1995 : @@ -2497,9 +2608,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & !! !! Again when dtconv is larger than tauadv, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. if(asqecflg(i)) then - tfac = tauadv(i) / dtconv(i) - tfac = min(tfac, 1.) - xmb(i) = -tfac * fld(i) / xk(i) +! tfac = tauadv(i) / dtconv(i) +! tfac = min(tfac, 1.) +! xmb(i) = -tfac * fld(i) / xk(i) + xmb(i) = -fld(i) / xk(i) endif enddo !! @@ -2713,10 +2825,9 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & rn(i) = rn(i) + rain * xmb(i) * .001 * dt2 endif if(flg(i) .and. k < ktcon(i)) then - evef = edt(i) * evfact - if(islimsk(i) == 1) evef=edt(i) * evfactl +! evef = edt(i) * evfact +! if(islimsk(i) == 1) evef=edt(i) * evfactl ! if(islimsk(i) == 1) evef=.07 -c if(islimsk(i) == 1) evef = 0. qcond(i) = evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) From f44d15a20f5d9b543ce08be58eb2e9c96262bef9 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:35:26 +0000 Subject: [PATCH 15/52] updated samfdeepcnv.meta --- physics/samfdeepcnv.meta | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index ff3c0d115..8a2a16fe7 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -375,6 +375,15 @@ type = integer intent = in optional = F +[hpbl] + standard_name = atmosphere_boundary_layer_thickness + long_name = PBL top height + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [ud_mf] standard_name = instantaneous_atmosphere_updraft_convective_mass_flux long_name = (updraft mass flux) * delt @@ -571,18 +580,9 @@ kind = kind_phys intent = in optional = F -[evfact] - standard_name = rain_evaporation_coefficient_deep_convection - long_name = convective rain evaporation coefficient for deep conv. - units = frac - dimensions = () - type = real - kind = kind_phys - intent = in - optional = F -[evfactl] - standard_name = rain_evaporation_coefficient_over_land_deep_convection - long_name = convective rain evaporation coefficient over land for deep conv. +[evef] + standard_name = rain_evaporation_coefficient_convection + long_name = convective rain evaporation coefficient for convection units = frac dimensions = () type = real From 64b3d2d183f742a1a5ecb5a45bedb85e4d074722 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:44:21 +0000 Subject: [PATCH 16/52] updated shallow convection --- physics/samfshalcnv.f | 315 ++++++++++++++++++++++++++++-------------- 1 file changed, 209 insertions(+), 106 deletions(-) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 1697cfe35..cfc7654f8 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -9,23 +9,25 @@ module samfshalcnv contains - subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, & + subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, +& & errmsg, errflg) integer, intent(in) :: imfshalcnv integer, intent(in) :: imfshalcnv_samf - + ! CCPP error handling character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Consistency checks if (imfshalcnv/=imfshalcnv_samf) then - write(errmsg,'(*(a))') 'Logic error: namelist choice of', & + write(errmsg,'(*(a))') 'Logic error: namelist choice of', +& & ' shallow convection is different from SAMF' errflg = 1 return - end if + end if end subroutine samfshalcnv_init subroutine samfshalcnv_finalize() @@ -61,7 +63,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,pgcon,asolfac,hwrf_samfshal,errmsg,errflg) + & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -87,7 +89,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:) ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & - & asolfac, pgcon + & asolfac, evef, pgcon logical, intent(in) :: hwrf_samfshal character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -108,8 +110,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & dz, dz1, e1, & el2orc, elocp, aafac, cm, & es, etah, h1, - & evef, evfact, evfactl, fact1, - & fact2, factor, dthk, +! & evfact, evfactl, + & fact1, fact2, factor, dthk, & gamma, pprime, betaw, & qlk, qrch, qs, & rfact, shear, tfac, @@ -120,30 +122,34 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & rho, tem, tem1, tem2, & ptem, ptem1 ! - integer kb(im), kbcon(im), kbcon1(im), + integer kb(im), kb1(im), kbcon(im), kbcon1(im), & ktcon(im), ktcon1(im), ktconn(im), & kbm(im), kmax(im) ! real(kind=kind_phys) aa1(im), cina(im), & tkemean(im), clamt(im), & ps(im), del(im,km), prsl(im,km), - & umean(im), tauadv(im), gdx(im), +! & umean(im), tauadv(im), gdx(im), + & gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), - & deltv(im), dtconv(im), edt(im), +! & deltv(im), dtconv(im), edt(im), + & deltv(im), dtconv(im), & pdot(im), po(im,km), & qcond(im), qevap(im), hmax(im), - & rntot(im), vshear(im), +! & rntot(im), vshear(im), + & rntot(im), & xlamud(im), xmb(im), xmbmax(im), & delebar(im,ntr), & delubar(im), delvbar(im) ! - real(kind=kind_phys) c0(im) + real(kind=kind_phys) c0(im), sfcpbl(im) c - real(kind=kind_phys) crtlamd + real(kind=kind_phys) crtlame, crtlamd ! real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, - & cinacr, cinacrmx, cinacrmn + & cinacr, cinacrmx, cinacrmn, + & sfclfac, rhcrt ! ! parameters for updraft velocity calculation real(kind=kind_phys) bet1, cd1, f1, gam1, @@ -172,10 +178,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) parameter(clamd=0.1,tkemx=0.65,tkemn=0.05) parameter(dtke=tkemx-tkemn) - parameter(dthk=25.) + parameter(dthk=25.,sfclfac=0.2,rhcrt=0.75) parameter(cinpcrmx=180.,cinpcrmn=120.) parameter(cinacrmx=-120.) - parameter(crtlamd=3.e-4) parameter(dtmax=10800.,dtmin=600.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrt=15.e3) @@ -194,6 +199,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km), real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), & dbyo(im,km), zo(im,km), xlamue(im,km), + & rh(im,km), & heo(im,km), heso(im,km), & dellah(im,km), dellaq(im,km), & dellae(im,km,ntr), @@ -203,6 +209,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eta(im,km), & zi(im,km), pwo(im,km), c0t(im,km), & sumx(im), tx1(im), cnvwt(im,km) + &, rhbar(im) ! logical do_aerosols, totflg, cnvflg(im), flg(im) ! @@ -255,6 +262,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kbot(i)=km+1 ktop(i)=0 endif + sfcpbl(i) = sfclfac * hpbl(i) rn(i)=0. kbcon(i)=km ktcon(i)=1 @@ -262,10 +270,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kb(i)=km pdot(i) = 0. qlko_ktcon(i) = 0. - edt(i) = 0. +! edt(i) = 0. aa1(i) = 0. cina(i) = 0. - vshear(i) = 0. +! vshear(i) = 0. gdx(i) = sqrt(garea(i)) scaldfunc(i)=-1.0 ! wang initialized sigmagfm(i)=-1.0 @@ -279,6 +287,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kbot(i)=km+1 ktop(i)=0 endif + sfcpbl(i) = sfclfac * hpbl(i) rn(i)=0. kbcon(i)=km ktcon(i)=1 @@ -286,10 +295,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kb(i)=km pdot(i) = 0. qlko_ktcon(i) = 0. - edt(i) = 0. +! edt(i) = 0.0 aa1(i) = 0. cina(i) = 0. - vshear(i) = 0. +! vshear(i) = 0. gdx(i) = sqrt(garea(i)) enddo endif @@ -342,14 +351,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & dt2 = delt ! c model tunable parameters are all here - if (hwrf_samfshal) then - aafac = .1 - else - aafac = .05 - endif -c evef = 0.07 - evfact = 0.3 - evfactl = 0.3 + aafac = .1 +! evfact = 0.3 +! evfactl = 0.3 +! + crtlame = 1.0e-4 + crtlamd = 3.0e-4 ! w1l = -8.e-3 w2l = -4.e-2 @@ -437,6 +444,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (cnvflg(i) .and. k <= kmax(i)) then pfld(i,k) = prsl(i,k) * 10.0 eta(i,k) = 1. + rh(i,k) = 0. hcko(i,k) = 0. qcko(i,k) = 0. qrcko(i,k)= 0. @@ -510,16 +518,34 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c this is the level where updraft starts c !> ## Perform calculations related to the updraft of the entraining/detraining cloud model ("static control"). +!> - Find the index for a level of sfclfac*hpbl which is initial guess for the parcel starting level. + do i=1,im + flg(i) = cnvflg(i) + kb1(i) = 1 + enddo + do k = 1, km1 + do i=1,im + if (flg(i) .and. zo(i,k) <= sfcpbl(i)) then + kb1(i) = k + else + flg(i) = .false. + endif + enddo + enddo + do i=1,im + kb1(i) = min(kb1(i),kpbl(i)) + enddo +c !> - Search in the PBL for the level of maximum moist static energy to start the ascending parcel. do i=1,im if (cnvflg(i)) then - hmax(i) = heo(i,1) - kb(i) = 1 + hmax(i) = heo(i,kb1(i)) + kb(i) = kb1(i) endif enddo do k = 2, km do i=1,im - if (cnvflg(i) .and. k <= kpbl(i)) then + if(cnvflg(i) .and. (k > kb1(i) .and. k <= kpbl(i))) then if(heo(i,k) > hmax(i)) then kb(i) = k hmax(i) = heo(i,k) @@ -561,6 +587,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & val2 = 1.e-10 qo(i,k) = max(qo(i,k), val2 ) ! qo(i,k) = min(qo(i,k),qeso(i,k)) + rh(i,k) = min(qo(i,k)/qeso(i,k), 1.) heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) + & cp * to(i,k) + hvap * qo(i,k) heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) + @@ -666,11 +693,95 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo if(totflg) return ! +! re-define kb & kbcon +! + do i=1,im + if (cnvflg(i)) then + hmax(i) = heo(i,1) + kb(i) = 1 + endif + enddo + do k = 2, km + do i=1,im + if (cnvflg(i) .and. k <= kpbl(i)) then + if(heo(i,k) > hmax(i)) then + kb(i) = k + hmax(i) = heo(i,k) + endif + endif + enddo + enddo +! + do i=1,im + flg(i) = cnvflg(i) + if(flg(i)) kbcon(i) = kmax(i) + enddo + do k = 2, km1 + do i=1,im + if (flg(i) .and. k < kbm(i)) then + if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then + kbcon(i) = k + flg(i) = .false. + endif + endif + enddo + enddo +! + do i=1,im + if(cnvflg(i)) then + if(kbcon(i) == kmax(i)) cnvflg(i) = .false. + endif + enddo +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! + do i=1,im + if(cnvflg(i)) then +! pdot(i) = 10.* dot(i,kbcon(i)) + pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s + endif + enddo +! +!> - if the mean relative humidity in the subcloud layers is less than a threshold value (rhcrt), convection is not triggered. +! + do i = 1, im + rhbar(i) = 0. + sumx(i) = 0. + enddo + do k = 1, km1 + do i = 1, im + if (cnvflg(i)) then + if(k >= kb(i) .and. k < kbcon(i)) then + dz = zo(i,k+1) - zo(i,k) + rhbar(i) = rhbar(i) + rh(i,k) * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo + do i= 1, im + if(cnvflg(i)) then + rhbar(i) = rhbar(i) / sumx(i) + if(rhbar(i) < rhcrt) then + cnvflg(i) = .false. + endif + endif + enddo +!! + totflg = .true. + do i=1,im + totflg = totflg .and. (.not. cnvflg(i)) + enddo + if(totflg) return +!! +! ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE ! -! - !c !c specify the detrainment rate for the updrafts !c @@ -735,7 +846,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im if(cnvflg(i)) then - xlamue(i,k) = clamt(i) / zi(i,k) + dz = zo(i,k+1) - zo(i,k) + xlamue(i,k) = clamt(i) / (zi(i,k) + dz) + xlamue(i,k) = max(xlamue(i,k), crtlame) endif enddo enddo @@ -1006,23 +1119,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c specify upper limit of mass flux at cloud base c !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. - if(hwrf_samfshal) then - do i = 1, im + do i = 1, im if(cnvflg(i)) then k = kbcon(i) dp = 1000. * del(i,k) xmbmax(i) = dp / (grav * dt2) endif - enddo - else - do i = 1, im - if(cnvflg(i)) then - k = kbcon(i) - dp = 1000. * del(i,k) - xmbmax(i) = dp / (2. * grav * dt2) - endif - enddo - endif + enddo c c compute cloud moisture property and precipitation c @@ -1349,34 +1452,34 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & !! E = 1.591 - 0.639\frac{\Delta V}{\Delta z} + 0.0953\left(\frac{\Delta V}{\Delta z}\right)^2 - 0.00496\left(\frac{\Delta V}{\Delta z}\right)^3 !! \f] !! where \f$\Delta V\f$ is the integrated horizontal shear over the cloud depth, \f$\Delta z\f$, (the ratio is converted to units of \f$10^{-3} s^{-1}\f$). The variable "edt" is \f$1-E\f$ and is constrained to the range \f$[0,0.9]\f$. - do i = 1, im - if(cnvflg(i)) then - vshear(i) = 0. - endif - enddo - do k = 2, km - do i = 1, im - if (cnvflg(i)) then - if(k > kb(i) .and. k <= ktcon(i)) then - shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 - & + (vo(i,k)-vo(i,k-1)) ** 2) - vshear(i) = vshear(i) + shear - endif - endif - enddo - enddo - do i = 1, im - if(cnvflg(i)) then - vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i))) - e1=1.591-.639*vshear(i) - & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3) - edt(i)=1.-e1 - val = .9 - edt(i) = min(edt(i),val) - val = .0 - edt(i) = max(edt(i),val) - endif - enddo +! do i = 1, im +! if(cnvflg(i)) then +! vshear(i) = 0. +! endif +! enddo +! do k = 2, km +! do i = 1, im +! if (cnvflg(i)) then +! if(k > kb(i) .and. k <= ktcon(i)) then +! shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 +! & + (vo(i,k)-vo(i,k-1)) ** 2) +! vshear(i) = vshear(i) + shear +! endif +! endif +! enddo +! enddo +! do i = 1, im +! if(cnvflg(i)) then +! vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i))) +! e1=1.591-.639*vshear(i) +! & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3) +! edt(i)=1.-e1 +! val = .9 +! edt(i) = min(edt(i),val) +! val = .0 +! edt(i) = max(edt(i),val) +! endif +! enddo c c--- what would the change be, that a cloud with unit mass c--- will do to the environment? @@ -1521,31 +1624,31 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. - do i= 1, im - if(cnvflg(i)) then - sumx(i) = 0. - umean(i) = 0. - endif - enddo - do k = 2, km1 - do i = 1, im - if(cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon1(i)) then - dz = zi(i,k) - zi(i,k-1) - tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) - umean(i) = umean(i) + tem * dz - sumx(i) = sumx(i) + dz - endif - endif - enddo - enddo - do i= 1, im - if(cnvflg(i)) then - umean(i) = umean(i) / sumx(i) - umean(i) = max(umean(i), 1.) - tauadv(i) = gdx(i) / umean(i) - endif - enddo +! do i= 1, im +! if(cnvflg(i)) then +! sumx(i) = 0. +! umean(i) = 0. +! endif +! enddo +! do k = 2, km1 +! do i = 1, im +! if(cnvflg(i)) then +! if(k >= kbcon1(i) .and. k < ktcon1(i)) then +! dz = zi(i,k) - zi(i,k-1) +! tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) +! umean(i) = umean(i) + tem * dz +! sumx(i) = sumx(i) + dz +! endif +! endif +! enddo +! enddo +! do i= 1, im +! if(cnvflg(i)) then +! umean(i) = umean(i) / sumx(i) +! umean(i) = max(umean(i), 1.) +! tauadv(i) = gdx(i) / umean(i) +! endif +! enddo c c compute cloud base mass flux as a function of the mean c updraft velcoity @@ -1556,9 +1659,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - tfac = tauadv(i) / dtconv(i) - tfac = min(tfac, 1.) - xmb(i) = tfac*betaw*rho*wc(i) +! tfac = tauadv(i) / dtconv(i) +! tfac = min(tfac, 1.) +! xmb(i) = tfac*betaw*rho*wc(i) + xmb(i) = betaw*rho*wc(i) endif enddo ! @@ -1722,10 +1826,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif endif if(flg(i) .and. k < ktcon(i)) then - evef = edt(i) * evfact - if(islimsk(i) == 1) evef=edt(i) * evfactl +! evef = edt(i) * evfact +! if(islimsk(i) == 1) evef=edt(i) * evfactl ! if(islimsk(i) == 1) evef=.07 -c if(islimsk(i) == 1) evef = 0. qcond(i) = evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) From b5fca2756388f8e2b4b73e93a5ad6e47ca50c87c Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:48:18 +0000 Subject: [PATCH 17/52] updated samfshalcnv.meta --- physics/samfshalcnv.meta | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index a454da3e7..686f3ae7c 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -430,6 +430,15 @@ kind = kind_phys intent = in optional = F +[evef] + standard_name = rain_evaporation_coefficient_convection + long_name = convective rain evaporation coefficient for convection + units = frac + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [pgcon] standard_name = momentum_transport_reduction_factor_pgf_shallow_convection long_name = reduction factor in momentum transport due to shal conv. induced pressure gradient force From 3e155b4264fd76df2bccc9bc7a1fe4854b211468 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:52:59 +0000 Subject: [PATCH 18/52] updated heat & evap --- physics/satmedmfvdif.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/satmedmfvdif.meta b/physics/satmedmfvdif.meta index baea94ad5..14b717f88 100644 --- a/physics/satmedmfvdif.meta +++ b/physics/satmedmfvdif.meta @@ -365,8 +365,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -374,7 +374,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From ed61f02c012a993c6dc77dc857ff6bd4181b2c6f Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 17:57:49 +0000 Subject: [PATCH 19/52] updated TKE-EDMF --- physics/satmedmfvdifq.F | 252 ++++++++++++++++++++++++---------------- 1 file changed, 154 insertions(+), 98 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 4bbfe61cc..c1ad33b08 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -35,7 +35,7 @@ subroutine satmedmfvdifq_init (satmedmf, ! Consistency checks if (.not. satmedmf) then - write(errmsg,fmt='(*(a))') 'Logic error: satmedmf = .false.' + write(errmsg,fmt='(*(a))') 'Logic error: satmedmf = .false.' errflg = 1 return end if @@ -69,8 +69,8 @@ end subroutine satmedmfvdifq_finalize !! @{ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & - & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,islimsk, & - & snwdph_lnd,psk,rbsoil,zorl,u10m,v10m,fm,fh, & + & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & + & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku, & @@ -86,7 +86,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !---------------------------------------------------------------------- integer, intent(in) :: im, km, ntrac, ntcw, ntiw, ntke, ntoz integer, intent(in) :: kinver(:) - integer, intent(in) :: islimsk(:) integer, intent(out) :: kpbl(:) logical, intent(in) :: gen_tend,ldiag3d,qdiag3d ! @@ -101,7 +100,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & t1(:,:), q1(:,:,:), & & swh(:,:), hlw(:,:), & & xmu(:), garea(:), & - & snwdph_lnd(:), & + & zvfun(:), & & psk(:), rbsoil(:), & & zorl(:), tsea(:), & & u10m(:), v10m(:), & @@ -116,8 +115,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & dt3dt(:,:), dq3dt(:,:), & & do3dt(:,:) real(kind=kind_phys), intent(out) :: & - & dusfc(:), dvsfc(:), & - & dtsfc(:), dqsfc(:), & + & dusfc(:), dvsfc(:), & + & dtsfc(:), dqsfc(:), & & hpbl(:) real(kind=kind_phys), intent(out) :: & & dkt(:,:), dku(:,:) @@ -151,7 +150,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & phih(im), phim(im), prn(im,km-1), & rbdn(im), rbup(im), thermal(im), & ustar(im), wstar(im), hpblx(im), - & ust3(im), wst3(im), + & ust3(im), wst3(im), rho_a(im), & z0(im), crb(im), & hgamt(im), hgamq(im), & wscale(im),vpert(im), @@ -168,7 +167,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & f1(im,km), f2(im,km*(ntrac-1)) ! real(kind=kind_phys) elm(im,km), ele(im,km), - & ckz(im,km), chz(im,km), frik(im), + & ckz(im,km), chz(im,km), & diss(im,km-1),prod(im,km-1), & bf(im,km-1), shr2(im,km-1), & xlamue(im,km-1), xlamde(im,km-1), @@ -212,41 +211,41 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & epsi, beta, chx, cqx, & rdt, rdz, qmin, qlmin, & rimin, rbcr, rbint, tdzmin, - & rlmn, rlmn1, rlmn2, + & rlmn, rlmn0, rlmn1, rlmn2, & rlmx, elmx, & ttend, utend, vtend, qtend, & zfac, zfmin, vk, spdk2, - & tkmin, tkminx, xkzinv, xkgdx, - & zlup, zldn, bsum, - & tem, tem1, tem2, + & tkmin, tkbmx, xkgdx, + & xkinv1, xkinv2, + & zlup, zldn, bsum, cs0, + & tem, tem1, tem2, tem3, & ptem, ptem0, ptem1, ptem2 -! - real(kind=kind_phys) xkzm_mp, xkzm_hp ! real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck ! - real(kind=kind_phys) qlcr, zstblmax + real(kind=kind_phys) qlcr, zstblmax, hcrinv ! real(kind=kind_phys) h1 !! - parameter(wfac=7.0,cfac=3.0) + parameter(wfac=7.0,cfac=4.5) parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) parameter(vk=0.4,rimin=-100.) parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) - parameter(rlmn=30.,rlmn1=5.,rlmn2=10.) + parameter(rlmn=30.,rlmn0=5.,rlmn1=5.,rlmn2=10.) parameter(rlmx=300.,elmx=300.) parameter(prmin=0.25,prmax=4.0) parameter(pr0=1.0,prtke=1.0,prscu=0.67) parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) - parameter(tkmin=1.e-9,tkminx=0.2,dspmax=10.0) + parameter(tkmin=1.e-9,tkbmx=0.2,dspmax=10.0) parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) parameter(aphi5=5.,aphi16=16.) parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=5000.) - parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.1) - parameter(h1=0.33333333) + parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.) + parameter(qlcr=3.5e-5,zstblmax=2500.) + parameter(xkinv1=0.15,xkinv2=0.3) + parameter(h1=0.33333333,hcrinv=250.) parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15) - parameter(ce0=0.4) + parameter(ce0=0.4,cs0=0.2) parameter(rchck=1.5,ndt=20) gravi=1.0/grav @@ -287,12 +286,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & buod(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 - rlmnz(i,k) = rlmn + rlmnz(i,k) = rlmn0 enddo enddo - do i=1,im - frik(i) = 1.0 - enddo do i=1,im zi(i,km+1) = phii(i,km+1) * gravi enddo @@ -331,41 +327,22 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo) -!> - set background diffusivities as a function of -!! horizontal grid size with xkzm_h & xkzm_m for gdx >= 25km -!! and 0.01 for gdx=5m, i.e., -!! \n xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) -!! \n xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) - - do i=1,im - xkzm_mp = xkzm_m - xkzm_hp = xkzm_h -! - if( islimsk(i) == 1 .and. snwdph_lnd(i) > 10.0 ) then ! over land - if (rbsoil(i) > 0. .and. rbsoil(i) <= 0.25) then - xkzm_mp = xkzm_m * (1.0 - rbsoil(i)/0.25)**2 + - & 0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2) - xkzm_hp = xkzm_h * (1.0 - rbsoil(i)/0.25)**2 + - & 0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2) - else if (rbsoil(i) > 0.25) then - xkzm_mp = 0.1 - xkzm_hp = 0.1 - endif - endif +!> - set background diffusivities with xkzm_h & xkzm_m for gdx >= xkgdx and +!! as a function of horizontal grid size for gdx < xkgdx +!! \n xkzm_hx = xkzm_h * (gdx / xkgdx) +!! \n xkzm_mx = xkzm_m * (gdx / xkgdx) ! + do i=1,im kx1(i) = 1 tx1(i) = 1.0 / prsi(i,1) tx2(i) = tx1(i) if(gdx(i) >= xkgdx) then - xkzm_hx(i) = xkzm_hp - xkzm_mx(i) = xkzm_mp + xkzm_hx(i) = xkzm_h + xkzm_mx(i) = xkzm_m else - tem = 1. / (xkgdx - 5.) - tem1 = (xkzm_hp - 0.01) * tem - tem2 = (xkzm_mp - 0.01) * tem - ptem = gdx(i) - 5. - xkzm_hx(i) = 0.01 + tem1 * ptem - xkzm_mx(i) = 0.01 + tem2 * ptem + tem = gdx(i) / xkgdx + xkzm_hx(i) = xkzm_h * tem + xkzm_mx(i) = xkzm_m * tem endif enddo do k = 1,km1 @@ -374,19 +351,18 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & xkzmo(i,k) = 0.0 if (k < kinver(i)) then ! minimum turbulent mixing length - ptem = prsl(i,k) * tx1(i) + ptem = prsi(i,k+1) * tx1(i) tem1 = 1.0 - ptem tem2 = tem1 * tem1 * 2.5 tem2 = min(1.0, exp(-tem2)) rlmnz(i,k)= rlmn * tem2 - rlmnz(i,k)= max(rlmnz(i,k), rlmn1) + rlmnz(i,k)= max(rlmnz(i,k), rlmn0) ! vertical background diffusivity - ptem = prsi(i,k+1) * tx1(i) - tem1 = 1.0 - ptem tem2 = tem1 * tem1 * 10.0 tem2 = min(1.0, exp(-tem2)) xkzo(i,k) = xkzm_hx(i) * tem2 -! vertical background diffusivity for momentum +! vertical background diffusivity for +! momentum if (ptem >= xkzm_s) then xkzmo(i,k) = xkzm_mx(i) kx1(i) = k + 1 @@ -399,10 +375,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & endif enddo enddo - +! !> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) + rho_a(i) = prsl(i,1) / (rd * t1(i,1)) dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. @@ -710,10 +687,54 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & hgamq(i) = evap(i)/wscale(i) vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1) vpert(i) = max(vpert(i),0.) - vpert(i) = min(cfac*vpert(i),gamcrt) + tem = min(cfac*vpert(i),gamcrt) + thermal(i)= thermal(i) + tem endif enddo ! +! enhance the pbl height by considering the thermal excess +! (overshoot pbl top) +! + do i=1,im + flg(i) = .true. + if(pcnvflg(i)) then + flg(i) = .false. + rbup(i) = rbsoil(i) + endif + enddo + do k = 2, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.) + rbup(i) = (thlvx(i,k)-thermal(i))* + & (g*zl(i,k)/thlvx(i,1))/spdk2 + kpbl(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + do i = 1,im + if(pcnvflg(i)) then + k = kpbl(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpbl(i) < zi(i,kpbl(i))) then + kpbl(i) = kpbl(i) - 1 + endif + if(kpbl(i) <= 1) then + pcnvflg(i) = .false. + pblflg(i) = .false. + endif + endif + enddo +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! look for stratocumulus !> ## Determine whether stratocumulus layers exist and compute quantities @@ -848,38 +869,43 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! -! background diffusivity decreasing with increasing surface layer stability -! -! do i = 1, im -! if(.not.sfcflg(i)) then -! tem = (1. + 5. * rbsoil(i))**2. -!! tem = (1. + 5. * zol(i))**2. -! frik(i) = 0.1 + 0.9 / tem -! endif -! enddo +! Above a threshold height (hcrinv), the background vertical +! diffusivities & mixing length +! in the inversion layers are set to much smaller values (xkinv1 & +! rlmn1) ! -! do k = 1,km1 -! do i=1,im -! xkzo(i,k) = frik(i) * xkzo(i,k) -! xkzmo(i,k)= frik(i) * xkzmo(i,k) -! enddo -! enddo -! -!> ## The background vertical diffusivities in the inversion layers are limited -!! to be less than or equal to xkzinv +! Below the threshold height (hcrinv), the background vertical +! diffusivities & mixing length +! in the inversion layers are increased with increasing roughness +! length & vegetation fraction ! do k = 1,km1 do i=1,im -! tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k) -! if(tem1 > 1.e-5) then - tem1 = tvx(i,k+1)-tvx(i,k) - if(tem1 > 0. .and. islimsk(i) /= 1) then - xkzo(i,k) = min(xkzo(i,k), xkzinv) - xkzmo(i,k) = min(xkzmo(i,k), xkzinv) - rlmnz(i,k) = min(rlmnz(i,k), rlmn2) + if(zi(i,k+1) > hcrinv) then + tem1 = tvx(i,k+1)-tvx(i,k) + if(tem1 >= 0.) then + xkzo(i,k) = min(xkzo(i,k), xkinv1) + xkzmo(i,k) = min(xkzmo(i,k), xkinv1) + rlmnz(i,k) = min(rlmnz(i,k), rlmn1) + endif + else + tem1 = tvx(i,k+1)-tvx(i,k) + if(tem1 > 0.) then + ptem = xkzo(i,k) * zvfun(i) + xkzo(i,k) = min(max(ptem, xkinv2), xkzo(i,k)) + ptem = xkzmo(i,k) * zvfun(i) + xkzmo(i,k) = min(max(ptem, xkinv2), xkzmo(i,k)) + ptem = rlmnz(i,k) * zvfun(i) + rlmnz(i,k) = min(max(ptem, rlmn2), rlmnz(i,k)) + endif endif enddo enddo + do k = 2,km1 + do i=1,im + rlmnz(i,k) = 0.5 * (rlmnz(i,k-1) + rlmnz(i,k)) + enddo + enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> ## Compute an asymtotic mixing length @@ -892,8 +918,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & do n = k, km1 if(mlenflg) then dz = zl(i,n+1) - zl(i,n) - ptem = gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))*dz -! ptem = gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k))*dz +! tem1 = 0.5 * (thvx(i,n) + thvx(i,n+1)) +!! tem1 = 0.5 * (thlvx(i,n) + thlvx(i,n+1)) + tem3=((u1(i,n+1)-u1(i,n))/dz)**2 + tem3=tem3+((v1(i,n+1)-v1(i,n))/dz)**2 + tem3=cs0*sqrt(tem3)*sqrt(tke(i,k)) + ptem = (gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))+tem3)*dz +! ptem = (gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k)+tem3)*dz +! ptem = (gotvx(i,n)*(tem1-thvx(i,k))+tem3)*dz +!! ptem = (gotvx(i,n)*(tem1-thlvx(i,k)+tem3)*dz bsum = bsum + ptem zlup = zlup + dz if(bsum >= tke(i,k)) then @@ -917,13 +950,23 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & if(n == 1) then dz = zl(i,1) tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) +! tem1 = 0.5 * (tem1 + thvx(i,n)) +!! tem1 = 0.5 * (tem1 + thlvx(i,n)) + tem3 = (u1(i,1)/dz)**2 + tem3 = tem3+(v1(i,1)/dz)**2 + tem3 = cs0*sqrt(tem3)*sqrt(tke(i,1)) else dz = zl(i,n) - zl(i,n-1) tem1 = thvx(i,n-1) ! tem1 = thlvx(i,n-1) +! tem1 = 0.5 * (thvx(i,n-1) + thvx(i,n)) +!! tem1 = 0.5 * (thlvx(i,n-1) + thlvx(i,n)) + tem3 = ((u1(i,n)-u1(i,n-1))/dz)**2 + tem3 = tem3+((v1(i,n)-v1(i,n-1))/dz)**2 + tem3 = cs0*sqrt(tem3)*sqrt(tke(i,k)) endif - ptem = gotvx(i,n)*(thvx(i,k)-tem1)*dz -! ptem = gotvx(i,n)*(thlvx(i,k)-tem1)*dz + ptem = (gotvx(i,n)*(thvx(i,k)-tem1)+tem3)*dz +! ptem = (gotvx(i,n)*(thlvx(i,k)-tem1)+tem3)*dz bsum = bsum + ptem zldn = zldn + dz if(bsum >= tke(i,k)) then @@ -954,6 +997,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel !! having an initial TKE can travel upward and downward before being stopped !! by buoyancy effects. +! +! Following Rodier et. al (2017), environmental wind shear effect on +! mixing length was included. +! ptem2 = min(zlup,zldn) rlam(i,k) = elmfac * ptem2 rlam(i,k) = max(rlam(i,k), tem1) @@ -1063,7 +1110,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & endif ptem = tem1 / (tem * elm(i,k)) tkmnz(i,k) = ptem * ptem - tkmnz(i,k) = min(tkmnz(i,k), tkminx) + tkmnz(i,k) = min(tkmnz(i,k), tkbmx) tkmnz(i,k) = max(tkmnz(i,k), tkmin) enddo enddo @@ -1432,10 +1479,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend - dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend - dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend +! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend +! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo + do i = 1,im + dtsfc(i) = rho_a(i) * cp * heat(i) + dqsfc(i) = rho_a(i) * hvap * evap(i) + enddo +! if(ldiag3d .and. .not. gen_tend) then do k = 1,km do i = 1,im @@ -1553,7 +1605,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & ! enddo enddo - c !> - Call tridi2() to solve tridiagonal problem for momentum c @@ -1567,10 +1618,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & vtend = (f2(i,k)-v1(i,k))*rdt du(i,k) = du(i,k)+utend dv(i,k) = dv(i,k)+vtend - dusfc(i) = dusfc(i)+conw*del(i,k)*utend - dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend +! dusfc(i) = dusfc(i)+conw*del(i,k)*utend +! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend enddo enddo + do i = 1,im + dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i) + dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i) + enddo +! if(ldiag3d .and. .not. gen_tend) then do k = 1,km do i = 1,im From fd2031051aebd898648c947e825e0d2b354e9b75 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:03:01 +0000 Subject: [PATCH 20/52] updated satmedmfvdifq.meta --- physics/satmedmfvdifq.meta | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index fd2dbe887..45a6fa5a1 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -290,18 +290,10 @@ kind = kind_phys intent = in optional = F -[islimsk] - standard_name = sea_land_ice_mask - long_name = sea/land/ice mask (=0/1/2) - units = flag - dimensions = (horizontal_loop_extent) - type = integer - intent = in - optional = F -[snwdph_lnd] - standard_name = surface_snow_thickness_water_equivalent_over_land - long_name = water equivalent snow depth over land - units = mm +[zvfun] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction + long_name = function of surface roughness length and green vegetation fraction + units = none dimensions = (horizontal_loop_extent) type = real kind = kind_phys @@ -380,8 +372,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -389,7 +381,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From c4f53940cd8b164ab023e08937be44ca43b58914 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:05:57 +0000 Subject: [PATCH 21/52] updated TKE-EDMF --- physics/satmedmfvdifq.F | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index c1ad33b08..a891fa387 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -150,7 +150,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & phih(im), phim(im), prn(im,km-1), & rbdn(im), rbup(im), thermal(im), & ustar(im), wstar(im), hpblx(im), - & ust3(im), wst3(im), rho_a(im), +! & ust3(im), wst3(im), rho_a(im), + & ust3(im), wst3(im), & z0(im), crb(im), & hgamt(im), hgamq(im), & wscale(im),vpert(im), @@ -379,7 +380,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) - rho_a(i) = prsl(i,1) / (rd * t1(i,1)) +! rho_a(i) = prsl(i,1) / (rd * t1(i,1)) dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. @@ -1479,14 +1480,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend -! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend -! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend + dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend + dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo - do i = 1,im - dtsfc(i) = rho_a(i) * cp * heat(i) - dqsfc(i) = rho_a(i) * hvap * evap(i) - enddo +! do i = 1,im +! dtsfc(i) = rho_a(i) * cp * heat(i) +! dqsfc(i) = rho_a(i) * hvap * evap(i) +! enddo ! if(ldiag3d .and. .not. gen_tend) then do k = 1,km @@ -1618,14 +1619,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & vtend = (f2(i,k)-v1(i,k))*rdt du(i,k) = du(i,k)+utend dv(i,k) = dv(i,k)+vtend -! dusfc(i) = dusfc(i)+conw*del(i,k)*utend -! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend + dusfc(i) = dusfc(i)+conw*del(i,k)*utend + dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend enddo enddo - do i = 1,im - dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i) - dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i) - enddo +! do i = 1,im +! dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i) +! dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i) +! enddo ! if(ldiag3d .and. .not. gen_tend) then do k = 1,km From 3175790051d3042b867b58e30d0155012a5d36bb Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:09:14 +0000 Subject: [PATCH 22/52] updated TKE-EDMF --- physics/satmedmfvdifq.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index a891fa387..ea61924aa 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -241,7 +241,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) parameter(aphi5=5.,aphi16=16.) parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.) + parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=3000.) parameter(qlcr=3.5e-5,zstblmax=2500.) parameter(xkinv1=0.15,xkinv2=0.3) parameter(h1=0.33333333,hcrinv=250.) From a3c4d52d586a25d5dbf13c629207c4e9c8439429 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:10:10 +0000 Subject: [PATCH 23/52] updated GFS surface layer scheme --- physics/sfc_diff.f | 129 ++++++++++++++++++++++++++++----------------- 1 file changed, 82 insertions(+), 47 deletions(-) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 93102e467..d5695319d 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -61,7 +61,7 @@ end subroutine sfc_diff_finalize !! - Calculate the exchange coefficients:\f$cm\f$, \f$ch\f$, and \f$stress\f$ as inputs of other \a sfc schemes. !! subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) - & ps,t1,q1,z1,wind, & !intent(in) + & ps,t1,q1,z1,garea,wind, & !intent(in) & prsl1,prslki,prsik1,prslk1, & !intent(in) & sigmaf,vegtype,shdmax,ivegsrc, & !intent(in) & z0pert,ztpert, & ! mg, sfc-perts !intent(in) @@ -70,7 +70,6 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & wet,dry,icy, & !intent(in) & tskin_wat, tskin_lnd, tskin_ice, & !intent(in) & tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in) - & snwdph_wat,snwdph_lnd,snwdph_ice, & !intent(in) & z0rl_wat, z0rl_lnd, z0rl_ice, & !intent(inout) & z0rl_wav, & !intent(inout) & ustar_wat, ustar_lnd, ustar_ice, & !intent(inout) @@ -98,13 +97,12 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) real(kind=kind_phys), dimension(:), intent(in) :: u10m,v10m real(kind=kind_phys), intent(in) :: rvrdm1, eps, epsm1, grav real(kind=kind_phys), dimension(:), intent(in) :: & - & ps,t1,q1,z1,prsl1,prslki,prsik1,prslk1, & + & ps,t1,q1,z1,garea,prsl1,prslki,prsik1,prslk1, & & wind,sigmaf,shdmax, & & z0pert,ztpert ! mg, sfc-perts real(kind=kind_phys), dimension(:), intent(in) :: & & tskin_wat, tskin_lnd, tskin_ice, & - & tsurf_wat, tsurf_lnd, tsurf_ice, & - & snwdph_wat,snwdph_lnd,snwdph_ice + & tsurf_wat, tsurf_lnd, tsurf_ice real(kind=kind_phys), dimension(:), intent(in) :: z0rl_wav real(kind=kind_phys), dimension(:), intent(inout) :: & @@ -125,14 +123,16 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! integer i ! - real(kind=kind_phys) :: rat, thv1, restar, wind10m, + real(kind=kind_phys) :: rat, tv1, thv1, restar, wind10m, & czilc, tem1, tem2, virtfac - real(kind=kind_phys) :: tvs, z0, z0max, ztmax + real(kind=kind_phys) :: tvs, z0, z0max, ztmax, zvfun, gdx +! + real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 ! real(kind=kind_phys), parameter :: & one=1.0_kp, zero=0.0_kp, half=0.5_kp, qmin=1.0e-8_kp - &, charnock=.014_kp, z0s_max=.317e-2_kp &! a limiting value at high winds over sea + &, charnock=.018_kp, z0s_max=.317e-2_kp &! a limiting value at high winds over sea &, zmin=1.0e-6_kp & &, vis=1.4e-5_kp, rnu=1.51e-5_kp, visi=one/vis & &, log01=log(0.01_kp), log05=log(0.05_kp), log07=log(0.07_kp) @@ -167,7 +167,10 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) do i=1,im if(flag_iter(i)) then virtfac = one + rvrdm1 * max(q1(i),qmin) - thv1 = t1(i) * prslki(i) * virtfac + tv1 = t1(i) * virtfac + thv1 = tv1 * prslki(i) + zvfun = zero + gdx = sqrt(garea(i)) ! compute stability dependent exchange coefficients ! this portion of the code is presently suppressed @@ -225,23 +228,34 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) z0max = max(z0max, zmin) -! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil - czilc = 0.8_kp - - tem1 = 1.0_kp - sigmaf(i) - ztmax = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) - +!! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil +! czilc = 0.8_kp +! tem1 = 1.0_kp - sigmaf(i) +! ztmax = z0max*exp( - tem1*tem1 +! & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) +! + czilc = 10.0_kp ** (- 4.0_kp * z0max) ! Trier et al. (2011,WAF) + czilc = min(czilc, 0.8_kp) + tem1 = 1.0_kp - sigmaf(i) + czilc = czilc * tem1 * tem1 + ztmax = z0max * exp( - czilc * ca + & * 258.2_kp * sqrt(ustar_lnd(i)*z0max) ) +! ! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land if (ztpert(i) /= zero) then ztmax = ztmax * (10.0_kp**ztpert(i)) endif ztmax = max(ztmax, zmin) +! + tem1 = (z0max - z0lo) / (z0up - z0lo) + tem1 = min(max(tem1, zero), 1.0_kp) + tem2 = max(sigmaf(i), 0.1_kp) + zvfun = sqrt(tem1 * tem2) ! call stability ! --- inputs: - & (z1(i), snwdph_lnd(i), thv1, wind(i), + & (z1(i), zvfun, gdx, tv1, thv1, wind(i), & z0max, ztmax, tvs, grav, ! --- outputs: & rb_lnd(i), fm_lnd(i), fh_lnd(i), fm10_lnd(i), fh2_lnd(i), @@ -265,18 +279,26 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) z0max = max(z0max, zmin) -! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height +!! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height ! dependance of czil - czilc = 0.8_kp +! czilc = 0.8_kp - tem1 = 1.0_kp - sigmaf(i) - ztmax = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) +! tem1 = 1.0_kp - sigmaf(i) +! ztmax = z0max*exp( - tem1*tem1 +! & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) +! + czilc = 10.0_kp ** (- 4.0_kp * z0max) + czilc = min(czilc, 0.8_kp) + tem1 = 1.0_kp - sigmaf(i) + czilc = czilc * tem1 * tem1 + ztmax = z0max * exp( - czilc * ca + & * 258.2_kp * sqrt(ustar_ice(i)*z0max) ) +! ztmax = max(ztmax, 1.0e-6) ! call stability ! --- inputs: - & (z1(i), snwdph_ice(i), thv1, wind(i), + & (z1(i), zvfun, gdx, tv1, thv1, wind(i), & z0max, ztmax, tvs, grav, ! --- outputs: & rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i), @@ -290,7 +312,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac z0 = 0.01_kp * z0rl_wat(i) z0max = max(zmin, min(z0,z1(i))) - ustar_wat(i) = sqrt(grav * z0 / charnock) +! ustar_wat(i) = sqrt(grav * z0 / charnock) wind10m = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i)) !** test xubin's new z0 @@ -320,7 +342,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! call stability ! --- inputs: - & (z1(i), snwdph_wat(i), thv1, wind(i), + & (z1(i), zvfun, gdx, tv1, thv1, wind(i), & z0max, ztmax, tvs, grav, ! --- outputs: & rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i), @@ -330,7 +352,10 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! if (sfc_z0_type >= 0) then if (sfc_z0_type == 0) then - z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) +! z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) + tem1 = 0.11 * vis / ustar_wat(i) + z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i) + ! mbek -- toga-coare flux algorithm ! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i) @@ -358,7 +383,9 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) endif elseif (z0rl_wav(i) <= 1.0e-7_kp) then - z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) +! z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) + tem1 = 0.11 * vis / ustar_wat(i) + z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i) if (redrag) then z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),1.0e-7_kp) @@ -380,7 +407,7 @@ end subroutine sfc_diff_run !>\ingroup GFS_diff_main subroutine stability & ! --- inputs: - & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, & + & ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & ! --- outputs: & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) !----- @@ -388,7 +415,7 @@ subroutine stability & integer, parameter :: kp = kind_phys ! --- inputs: real(kind=kind_phys), intent(in) :: & - & z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav + & z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav ! --- outputs: real(kind=kind_phys), intent(out) :: & @@ -397,27 +424,41 @@ subroutine stability & ! --- locals: real(kind=kind_phys), parameter :: alpha=5.0_kp, a0=-3.975_kp & &, a1=12.32_kp, alpha4=4.0_kp*alpha & - &, b1=-7.755_kp, b2=6.041_kp, alpha2=alpha+alpha & - &, beta=1.0_kp & + &, b1=-7.755_kp, b2=6.041_kp & + &, xkrefsqr=0.3_kp, xkmin=0.05_kp & + &, xkgdx=3000.0_kp & &, a0p=-7.941_kp, a1p=24.75_kp, b1p=-8.705_kp, b2p=7.899_kp& - &, ztmin1=-999.0_kp, zero=0.0_kp, one=1.0_kp + &, zolmin=-10.0_kp, zero=0.0_kp, one=1.0_kp real(kind=kind_phys) aa, aa0, bb, bb0, dtv, adtv, & hl1, hl12, pm, ph, pm10, ph2, & z1i, & fms, fhs, hl0, hl0inf, hlinf, & hl110, hlt, hltinf, olinf, - & tem1, tem2, ztmax1 + & tem1, tem2, zolmax + + real(kind=kind_phys) xkzo z1i = one / z1 - tem1 = z0max/z1 - if (abs(one-tem1) > 1.0e-6_kp) then - ztmax1 = - beta*log(tem1)/(alpha2*(one-tem1)) +! +! set background diffusivities with one for gdx >= xkgdx and +! as a function of horizontal grid size for gdx < xkgdx +! (i.e., gdx/xkgdx for gdx < xkgdx) +! + if(gdx >= xkgdx) then + xkzo = one else - ztmax1 = 99.0_kp + xkzo = gdx / xkgdx endif - if( z0max < 0.05_kp .and. snwdph < 10.0_kp ) ztmax1 = 99.0_kp + + tem1 = tv1 - tvs + if(tem1 > zero) then + tem2 = xkzo * zvfun + xkzo = min(max(tem2, xkmin), xkzo) + endif + + zolmax = xkrefsqr / sqrt(xkzo) ! compute stability indices (rb and hlinf) @@ -438,7 +479,7 @@ subroutine stability & fm10 = log((z0max+10.0_kp) * tem1) fh2 = log((ztmax+2.0_kp) * tem2) hlinf = rb * fm * fm / fh - hlinf = min(max(hlinf,ztmin1),ztmax1) + hlinf = min(max(hlinf,zolmin),zolmax) ! ! stable case ! @@ -457,7 +498,7 @@ subroutine stability & fms = fm - pm fhs = fh - ph hl1 = fms * fms * rb / fhs - hl1 = min(max(hl1, ztmin1), ztmax1) + hl1 = min(hl1, zolmax) endif ! ! second iteration @@ -472,11 +513,9 @@ subroutine stability & pm = aa0 - aa + log( (one+aa)/(one+aa0) ) ph = bb0 - bb + log( (one+bb)/(one+bb0) ) hl110 = hl1 * 10.0_kp * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) aa = sqrt(one + alpha4 * hl110) pm10 = aa0 - aa + log( (one+aa)/(one+aa0) ) hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12,ztmin1),ztmax1) ! aa = sqrt(one + alpha4 * hl12) bb = sqrt(one + alpha4 * hl12) ph2 = bb0 - bb + log( (one+bb)/(one+bb0) ) @@ -488,7 +527,7 @@ subroutine stability & tem1 = 50.0_kp * z0max if(abs(olinf) <= tem1) then hlinf = -z1 / tem1 - hlinf = min(max(hlinf,ztmin1),ztmax1) + hlinf = max(hlinf, zolmin) endif ! ! get pm and ph @@ -498,10 +537,8 @@ subroutine stability & pm = (a0 + a1*hl1) * hl1 / (one+ (b1+b2*hl1) *hl1) ph = (a0p + a1p*hl1) * hl1 / (one+ (b1p+b2p*hl1)*hl1) hl110 = hl1 * 10.0_kp * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) pm10 = (a0 + a1*hl110) * hl110/(one+(b1+b2*hl110)*hl110) hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = (a0p + a1p*hl12) * hl12/(one+(b1p+b2p*hl12)*hl12) else ! hlinf < 0.05 hl1 = -hlinf @@ -511,11 +548,9 @@ subroutine stability & ! pm = log(hl1) + 2.0 * hl1 ** (-.25) - .8776 ! ph = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386 hl110 = hl1 * 10.0_kp * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) pm10 = log(hl110) + 2.0_kp/sqrt(sqrt(hl110)) - 0.8776_kp ! pm10 = log(hl110) + 2. * hl110 ** (-.25) - .8776 hl12 = (hl1+hl1) * z1i - hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = log(hl12) + 0.5_kp / sqrt(hl12) + 1.386_kp ! ph2 = log(hl12) + .5 * hl12 ** (-.5) + 1.386 endif From f16f1af56d14c2dec4a883f8583ac8a7bd48a0aa Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:13:27 +0000 Subject: [PATCH 24/52] updated sfc_diff.meta --- physics/sfc_diff.meta | 36 +++++++++--------------------------- 1 file changed, 9 insertions(+), 27 deletions(-) diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta index 342eaeea5..0ec39aa5c 100644 --- a/physics/sfc_diff.meta +++ b/physics/sfc_diff.meta @@ -87,6 +87,15 @@ kind = kind_phys intent = in optional = F +[garea] + standard_name = cell_area + long_name = area of the grid cell + units = m2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [wind] standard_name = wind_speed_at_lowest_model_layer long_name = wind speed at lowest model level @@ -304,33 +313,6 @@ kind = kind_phys intent = in optional = F -[snwdph_wat] - standard_name = surface_snow_thickness_water_equivalent_over_water - long_name = water equivalent snow depth over water - units = mm - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F -[snwdph_lnd] - standard_name = surface_snow_thickness_water_equivalent_over_land - long_name = water equivalent snow depth over land - units = mm - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F -[snwdph_ice] - standard_name = surface_snow_thickness_water_equivalent_over_ice - long_name = water equivalent snow depth over ice - units = mm - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [z0rl_wat] standard_name = surface_roughness_length_over_water long_name = surface roughness length over water (temporary use as interstitial) From 2a518b28c75f666396f326b958ffcfae09b270cf Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:17:41 +0000 Subject: [PATCH 25/52] add sea spray effect parameterization --- physics/sfc_nst.f | 53 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/physics/sfc_nst.f b/physics/sfc_nst.f index f03e725f3..663276224 100644 --- a/physics/sfc_nst.f +++ b/physics/sfc_nst.f @@ -28,6 +28,7 @@ end subroutine sfc_nst_finalize subroutine sfc_nst_run & & ( im, hvap, cp, hfus, jcal, eps, epsm1, rvrdm1, rd, rhw0, & ! --- inputs: & pi, tgice, sbc, ps, u1, v1, t1, q1, tref, cm, ch, & + & lseaspray, fm, fm10, & & prsl1, prslki, prsik1, prslk1, wet, use_flake, xlon, & & sinlat, stress, & & sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, & @@ -47,6 +48,7 @@ subroutine sfc_nst_run & ! call sfc_nst ! ! inputs: ! ! ( im, ps, u1, v1, t1, q1, tref, cm, ch, ! +! lseaspray, fm, fm10, ! ! prsl1, prslki, wet, use_flake, xlon, sinlat, stress, ! ! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, ! ! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, ! @@ -89,6 +91,10 @@ subroutine sfc_nst_run & ! tref - real, reference/foundation temperature ( k ) im ! ! cm - real, surface exchange coeff for momentum (m/s) im ! ! ch - real, surface exchange coeff heat & moisture(m/s) im ! +! lseaspray- logical, .t. for parameterization for sea spray 1 ! +! fm - real, a stability profile function for momentum im ! +! fm10 - real, a stability profile function for momentum im ! +! at 10m ! ! prsl1 - real, surface layer mean pressure (pa) im ! ! prslki - real, im ! ! prsik1 - real, im ! @@ -189,12 +195,15 @@ subroutine sfc_nst_run & real (kind=kind_phys), intent(in) :: hvap, cp, hfus, jcal, eps, & & epsm1, rvrdm1, rd, rhw0, sbc, pi, tgice real (kind=kind_phys), dimension(:), intent(in) :: ps, u1, v1, & - & t1, q1, tref, cm, ch, prsl1, prslki, prsik1, prslk1, & - & xlon,xcosz, & + & t1, q1, tref, cm, ch, fm, fm10, & + & prsl1, prslki, prsik1, prslk1, xlon, xcosz, & & sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, wind real (kind=kind_phys), intent(in) :: timestep real (kind=kind_phys), intent(in) :: solhr +! For sea spray effect + logical, intent(in) :: lseaspray +! logical, dimension(:), intent(in) :: flag_iter, flag_guess, wet, & & use_flake ! &, icy @@ -245,6 +254,18 @@ subroutine sfc_nst_run & ! external functions called: iw3jdn integer :: iw3jdn +! +! parameters for sea spray effect +! + real (kind=kind_phys) :: f10m, u10m, v10m, ws10, ru10, qss1, + & bb1, hflxs, evaps, ptem +! +! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.1, +! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.0, +! real (kind=kind_phys), parameter :: alps=1.0, bets=1.0, gams=0.2, + real (kind=kind_phys), parameter :: alps=0.75,bets=0.75,gams=0.15, + & ws10cr=30., conlf=7.2e-9, consf=6.4e-8 +! !====================================================================================================== cc ! Initialize CCPP error handling variables @@ -629,7 +650,33 @@ subroutine sfc_nst_run & endif enddo endif ! if ( nstf_name1 > 1 ) then - +! +! include sea spray effects +! + do i=1,im + if(lseaspray .and. flag(i)) then + f10m = fm10(i) / fm(i) + u10m = f10m * u1(i) + v10m = f10m * v1(i) + ws10 = sqrt(u10m*u10m + v10m*v10m) + ws10 = max(ws10,1.) + ws10 = min(ws10,ws10cr) + tem = .015 * ws10 * ws10 + ru10 = 1. - .087 * log(10./tem) + qss1 = fpvs(t1(i)) + qss1 = eps * qss1 / (prsl1(i) + epsm1 * qss1) + tem = rd * cp * t1(i) * t1(i) + tem = 1. + eps * hvap * hvap * qss1 / tem + bb1 = 1. / tem + evaps = conlf * (ws10**5.4) * ru10 * bb1 + evaps = evaps * rho_a(i) * hvap * (qss1 - q0(i)) + evap(i) = evap(i) + alps * evaps + hflxs = consf * (ws10**3.4) * ru10 + hflxs = hflxs * rho_a(i) * cp * (tskin(i) - t1(i)) + ptem = alps - gams + hflx(i) = hflx(i) + bets * hflxs - ptem * evaps + endif + enddo ! do i=1,im if ( flag(i) ) then From 20189cbc29f374cf5f596d79c5aa2f9baa1b4207 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:19:26 +0000 Subject: [PATCH 26/52] updated sfc_nst.meta --- physics/sfc_nst.meta | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/physics/sfc_nst.meta b/physics/sfc_nst.meta index a29f10f90..08dde266f 100644 --- a/physics/sfc_nst.meta +++ b/physics/sfc_nst.meta @@ -195,6 +195,32 @@ kind = kind_phys intent = in optional = F +[lseaspray] + standard_name = flag_for_sea_spray + long_name = flag for sea spray parameterization + units = flag + dimensions = () + type = logical + intent = in + optional = F +[fm] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_water + long_name = Monin-Obukhov similarity function for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[fm10] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_water + long_name = Monin-Obukhov similarity parameter for momentum at 10m over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [prsl1] standard_name = air_pressure_at_lowest_model_layer long_name = surface layer mean pressure From 3d1e7f3b0d2fedac9a0cc1afe7466e1239777c12 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:21:45 +0000 Subject: [PATCH 27/52] add sea spray effect parameterization --- physics/sfc_ocean.F | 97 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 81 insertions(+), 16 deletions(-) diff --git a/physics/sfc_ocean.F b/physics/sfc_ocean.F index 67a6df04f..3f0fd23bc 100644 --- a/physics/sfc_ocean.F +++ b/physics/sfc_ocean.F @@ -26,8 +26,9 @@ end subroutine sfc_ocean_finalize subroutine sfc_ocean_run & !................................... ! --- inputs: - & ( im, rd, eps, epsm1, rvrdm1, ps, t1, q1, & - & tskin, cm, ch, prsl1, prslki, wet, use_flake, wind, &, ! --- inputs + & ( im, hvap, cp, rd, eps, epsm1, rvrdm1, ps, u1, v1, t1, q1, & + & tskin, cm, ch, lseaspray, fm, fm10, & + & prsl1, prslki, wet, use_flake, wind, &, ! --- inputs & flag_iter, & & qsurf, cmm, chh, gflux, evap, hflx, ep, & ! --- outputs & errmsg, errflg & @@ -40,8 +41,7 @@ subroutine sfc_ocean_run & ! ! ! call sfc_ocean ! ! inputs: ! -! ( im, ps, t1, q1, tskin, cm, ch, ! -!! ( im, ps, u1, v1, t1, q1, tskin, cm, ch, ! +! ( im, ps, u1, v1, t1, q1, tskin, cm, ch, lseaspray, fm, fm10, ! ! prsl1, prslki, wet, use_flake, wind, flag_iter, ! ! outputs: ! ! qsurf, cmm, chh, gflux, evap, hflx, ep ) ! @@ -65,11 +65,16 @@ subroutine sfc_ocean_run & ! inputs: size ! ! im - integer, horizontal dimension 1 ! ! ps - real, surface pressure im ! +! u1, v1 - real, u/v component of surface layer wind (m/s) im ! ! t1 - real, surface layer mean temperature ( k ) im ! ! q1 - real, surface layer mean specific humidity im ! ! tskin - real, ground surface skin temperature ( k ) im ! ! cm - real, surface exchange coeff for momentum (m/s) im ! ! ch - real, surface exchange coeff heat & moisture(m/s) im ! +! lseaspray- logical, .t. for parameterization for sea spray 1 ! +! fm - real, a stability profile function for momentum im ! +! fm10 - real, a stability profile function for momentum im ! +! at 10m ! ! prsl1 - real, surface layer mean pressure im ! ! prslki - real, im ! ! wet - logical, =T if any ocean/lak, =F otherwise im ! @@ -97,11 +102,14 @@ subroutine sfc_ocean_run & &, qmin = 1.0e-8_kind_phys ! --- inputs: integer, intent(in) :: im - real (kind=kind_phys), intent(in) :: rd, eps, epsm1, rvrdm1 + real (kind=kind_phys), intent(in) :: hvap, cp, rd, eps, epsm1, rvrdm1 - real (kind=kind_phys), dimension(:), intent(in) :: ps, & - & t1, q1, tskin, cm, ch, prsl1, prslki, wind + real (kind=kind_phys), dimension(:), intent(in) :: ps, u1, v1, & + & t1, q1, tskin, cm, ch, fm, fm10, prsl1, prslki, wind +! For sea spray effect + logical, intent(in) :: lseaspray +! logical, dimension(:), intent(in) :: flag_iter, wet, use_flake ! --- outputs: @@ -113,48 +121,105 @@ subroutine sfc_ocean_run & ! --- locals: - real (kind=kind_phys) :: q0, qss, rch, rho, tem + real (kind=kind_phys) :: qss, rch, tem, + & elocp, cpinv, hvapi + real (kind=kind_phys), dimension(im) :: rho, q0 integer :: i + logical :: flag(im) +! +! parameters for sea spray effect +! + real (kind=kind_phys) :: f10m, u10m, v10m, ws10, ru10, qss1, + & bb1, hflxs, evaps, ptem +! +! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.1, +! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.0, +! real (kind=kind_phys), parameter :: alps=1.0, bets=1.0, gams=0.2, + real (kind=kind_phys), parameter :: alps=0.75,bets=0.75,gams=0.15, + & ws10cr=30., conlf=7.2e-9, consf=6.4e-8 +! +!====================================================================================================== !===> ... begin here ! ! -- ... initialize CCPP error handling variables errmsg = '' errflg = 0 + + cpinv = one/cp + hvapi = one/hvap + elocp = hvap/cp ! ! --- ... flag for open water do i = 1, im - + flag(i) = (wet(i) .and. flag_iter(i) .and. .not. use_flake(i)) ! --- ... initialize variables. all units are supposedly m.k.s. unless specified ! ps is in pascals, wind is wind speed, ! rho is density, qss is sat. hum. at surface - if (wet(i) .and. flag_iter(i) .and. .not. use_flake(i)) then - q0 = max( q1(i), qmin ) - rho = prsl1(i) / (rd*t1(i)*(one + rvrdm1*q0)) + if ( flag(i) ) then + q0(i) = max( q1(i), qmin ) + rho(i) = prsl1(i) / (rd*t1(i)*(one + rvrdm1*q0(i))) qss = fpvs( tskin(i) ) qss = eps*qss / (ps(i) + epsm1*qss) ! --- ... rcp = rho cp ch v + rch = rho(i) * cp * ch(i) * wind(i) tem = ch(i) * wind(i) cmm(i) = cm(i) * wind(i) - chh(i) = rho * tem + chh(i) = rho(i) * tem ! --- ... sensible and latent heat flux over open water - hflx(i) = tem * (tskin(i) - t1(i) * prslki(i)) + hflx(i) = rch * (tskin(i) - t1(i) * prslki(i)) - evap(i) = tem * (qss - q0) + evap(i) = elocp * rch * (qss - q0(i)) - ep(i) = evap(i) qsurf(i) = qss gflux(i) = zero endif enddo ! +! include sea spray effects +! + do i=1,im + if(lseaspray .and. flag(i)) then + f10m = fm10(i) / fm(i) + u10m = f10m * u1(i) + v10m = f10m * v1(i) + ws10 = sqrt(u10m*u10m + v10m*v10m) + ws10 = max(ws10,1.) + ws10 = min(ws10,ws10cr) + tem = .015 * ws10 * ws10 + ru10 = 1. - .087 * log(10./tem) + qss1 = fpvs(t1(i)) + qss1 = eps * qss1 / (prsl1(i) + epsm1 * qss1) + tem = rd * cp * t1(i) * t1(i) + tem = 1. + eps * hvap * hvap * qss1 / tem + bb1 = 1. / tem + evaps = conlf * (ws10**5.4) * ru10 * bb1 + evaps = evaps * rho(i) * hvap * (qss1 - q0(i)) + evap(i) = evap(i) + alps * evaps + hflxs = consf * (ws10**3.4) * ru10 + hflxs = hflxs * rho(i) * cp * (tskin(i) - t1(i)) + ptem = alps - gams + hflx(i) = hflx(i) + bets * hflxs - ptem * evaps + endif + enddo +! + do i = 1, im + if ( flag(i) ) then + tem = 1.0 / rho(i) + hflx(i) = hflx(i) * tem * cpinv + evap(i) = evap(i) * tem * hvapi + ep(i) = evap(i) + endif + enddo +! + return !................................... end subroutine sfc_ocean_run From f4bdbf5ce883c98a6764737c89c3139020774d32 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:24:07 +0000 Subject: [PATCH 28/52] updated sfc_ocean.meta --- physics/sfc_ocean.meta | 62 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/physics/sfc_ocean.meta b/physics/sfc_ocean.meta index f27c2207d..0fcc4e646 100644 --- a/physics/sfc_ocean.meta +++ b/physics/sfc_ocean.meta @@ -15,6 +15,24 @@ type = integer intent = in optional = F +[hvap] + standard_name = latent_heat_of_vaporization_of_water_at_0C + long_name = latent heat of evaporation/sublimation + units = J kg-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[cp] + standard_name = specific_heat_of_dry_air_at_constant_pressure + long_name = specific heat of dry air at constant pressure + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [rd] standard_name = gas_constant_dry_air long_name = ideal gas constant for dry air @@ -60,6 +78,24 @@ kind = kind_phys intent = in optional = F +[u1] + standard_name = x_wind_at_lowest_model_layer + long_name = x component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[v1] + standard_name = y_wind_at_lowest_model_layer + long_name = y component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [t1] standard_name = air_temperature_at_lowest_model_layer long_name = surface layer mean temperature @@ -105,6 +141,32 @@ kind = kind_phys intent = in optional = F +[lseaspray] + standard_name = flag_for_sea_spray + long_name = flag for sea spray parameterization + units = flag + dimensions = () + type = logical + intent = in + optional = F +[fm] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_water + long_name = Monin-Obukhov similarity function for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[fm10] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_water + long_name = Monin-Obukhov similarity parameter for momentum at 10m over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [prsl1] standard_name = air_pressure_at_lowest_model_layer long_name = surface layer mean pressure From f20791979152db0d4a2e41f3f8cfbad12037ce4b Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:26:21 +0000 Subject: [PATCH 29/52] updated heat & evap --- physics/shalcnv.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/shalcnv.meta b/physics/shalcnv.meta index 7986d28f8..b1a5332ed 100644 --- a/physics/shalcnv.meta +++ b/physics/shalcnv.meta @@ -351,8 +351,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -360,7 +360,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From 9b889271cc704aa0a3a4dc398933cfbbfd088a25 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:29:15 +0000 Subject: [PATCH 30/52] updated heat & evap --- physics/shinhongvdif.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/shinhongvdif.meta b/physics/shinhongvdif.meta index 6b12b64f5..e7d1baccc 100644 --- a/physics/shinhongvdif.meta +++ b/physics/shinhongvdif.meta @@ -265,8 +265,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -274,7 +274,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From baaaf6b4a4e551d6feefff033c97e39401f8c460 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 21 May 2021 18:31:39 +0000 Subject: [PATCH 31/52] updated heat & evap --- physics/ysuvdif.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/ysuvdif.meta b/physics/ysuvdif.meta index 1ee952d45..ed8ada624 100644 --- a/physics/ysuvdif.meta +++ b/physics/ysuvdif.meta @@ -292,8 +292,8 @@ intent = in optional = F [heat] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real @@ -301,7 +301,7 @@ intent = in optional = F [evap] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness + standard_name = kinematic_surface_upward_latent_heat_flux long_name = kinematic surface upward latent heat flux units = kg kg-1 m s-1 dimensions = (horizontal_loop_extent) From 16a71292056fd84ee2597203836b391fc8c343bd Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Mon, 24 May 2021 02:31:53 +0000 Subject: [PATCH 32/52] updated TKE-EDMF --- physics/satmedmfvdifq.F | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index ea61924aa..f6ab554ca 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -150,8 +150,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & phih(im), phim(im), prn(im,km-1), & rbdn(im), rbup(im), thermal(im), & ustar(im), wstar(im), hpblx(im), -! & ust3(im), wst3(im), rho_a(im), - & ust3(im), wst3(im), + & ust3(im), wst3(im), rho_a(im), & z0(im), crb(im), & hgamt(im), hgamq(im), & wscale(im),vpert(im), @@ -380,7 +379,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) -! rho_a(i) = prsl(i,1) / (rd * t1(i,1)) + rho_a(i) = prsl(i,1) / (rd * t1(i,1)) dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. @@ -1480,14 +1479,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend - dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend - dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend +! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend +! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo -! do i = 1,im -! dtsfc(i) = rho_a(i) * cp * heat(i) -! dqsfc(i) = rho_a(i) * hvap * evap(i) -! enddo + do i = 1,im + dtsfc(i) = rho_a(i) * cp * heat(i) + dqsfc(i) = rho_a(i) * hvap * evap(i) + enddo ! if(ldiag3d .and. .not. gen_tend) then do k = 1,km @@ -1619,14 +1618,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & vtend = (f2(i,k)-v1(i,k))*rdt du(i,k) = du(i,k)+utend dv(i,k) = dv(i,k)+vtend - dusfc(i) = dusfc(i)+conw*del(i,k)*utend - dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend +! dusfc(i) = dusfc(i)+conw*del(i,k)*utend +! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend enddo enddo -! do i = 1,im -! dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i) -! dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i) -! enddo + do i = 1,im + dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i) + dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i) + enddo ! if(ldiag3d .and. .not. gen_tend) then do k = 1,km From 108a82f2b8aa0a896abd601d36c05d18c1b8b325 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Thu, 27 May 2021 12:29:41 +0000 Subject: [PATCH 33/52] Fix compilation failure --- physics/sfc_ocean.F | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/physics/sfc_ocean.F b/physics/sfc_ocean.F index 3f0fd23bc..79a9eb295 100644 --- a/physics/sfc_ocean.F +++ b/physics/sfc_ocean.F @@ -98,11 +98,12 @@ subroutine sfc_ocean_run & implicit none ! --- constant parameters: - real (kind=kind_phys), parameter :: one = 1.0_kind_phys, zero = 0.0_kind_phys & - &, qmin = 1.0e-8_kind_phys + real (kind=kind_phys), parameter :: one = 1.0_kind_phys, & + & zero = 0.0_kind_phys, qmin = 1.0e-8_kind_phys ! --- inputs: integer, intent(in) :: im - real (kind=kind_phys), intent(in) :: hvap, cp, rd, eps, epsm1, rvrdm1 + real (kind=kind_phys), intent(in) :: hvap, cp, rd, & + & eps, epsm1, rvrdm1 real (kind=kind_phys), dimension(:), intent(in) :: ps, u1, v1, & & t1, q1, tskin, cm, ch, fm, fm10, prsl1, prslki, wind From bc605882a92c47ff880688bc1f0dc3309531e408 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Thu, 27 May 2021 12:40:11 +0000 Subject: [PATCH 34/52] Fix line continuation --- physics/samfdeepcnv.f | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index a8905696c..edcf55d66 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -10,8 +10,7 @@ module samfdeepcnv contains - subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, -& + subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, & & errmsg, errflg) integer, intent(in) :: imfdeepcnv @@ -22,8 +21,7 @@ subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, ! Consistency checks if (imfdeepcnv/=imfdeepcnv_samf) then - write(errmsg,'(*(a))') 'Logic error: namelist choice of', -& + write(errmsg,'(*(a))') 'Logic error: namelist choice of', & & ' deep convection is different from SAMF scheme' errflg = 1 return From 5a3c4e0dacb25cb3c1eee7048bee884bd2c4cc1e Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Thu, 27 May 2021 12:45:40 +0000 Subject: [PATCH 35/52] Fix line continuation problem --- physics/samfshalcnv.f | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index cfc7654f8..7fd60f5f8 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -9,8 +9,7 @@ module samfshalcnv contains - subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, -& + subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, & & errmsg, errflg) integer, intent(in) :: imfshalcnv @@ -22,8 +21,7 @@ subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, ! Consistency checks if (imfshalcnv/=imfshalcnv_samf) then - write(errmsg,'(*(a))') 'Logic error: namelist choice of', -& + write(errmsg,'(*(a))') 'Logic error: namelist choice of', & & ' shallow convection is different from SAMF' errflg = 1 return From 916a5cd8de49c953d529b4462ddc0c994e507ce5 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Thu, 27 May 2021 13:45:00 +0000 Subject: [PATCH 36/52] Updated TKE-EDMF --- physics/satmedmfvdifq.F | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index f6ab554ca..a59c43e53 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -179,6 +179,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & ucko(im,km), vcko(im,km), & buou(im,km), xmf(im,km) ! +! temporary local variables for instantaneous surface fluxes +! (they will be removed later) +! + real(kind=kind_phys) dusfc1(im),dvsfc1(im),dtsfc1(im),dqsfc1(im) +! ! variables for stratocumulus-top induced downdrafts ! real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac), @@ -379,11 +384,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) - rho_a(i) = prsl(i,1) / (rd * t1(i,1)) - dusfc(i) = 0. - dvsfc(i) = 0. - dtsfc(i) = 0. - dqsfc(i) = 0. + rho_a(i) = prsl(i,1)/(rd*t1(i,1)*(1.+fv*max(q1(i,1,1),qmin))) + dusfc1(i) = 0. + dvsfc1(i) = 0. + dtsfc1(i) = 0. + dqsfc1(i) = 0. kpbl(i) = 1 hpbl(i) = 0. kpblx(i) = 1 @@ -1479,8 +1484,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend -! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend -! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend + dtsfc1(i) = dtsfc1(i)+cont*del(i,k)*ttend + dqsfc1(i) = dqsfc1(i)+conq*del(i,k)*qtend enddo enddo do i = 1,im @@ -1618,8 +1623,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & vtend = (f2(i,k)-v1(i,k))*rdt du(i,k) = du(i,k)+utend dv(i,k) = dv(i,k)+vtend -! dusfc(i) = dusfc(i)+conw*del(i,k)*utend -! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend + dusfc1(i) = dusfc1(i)+conw*del(i,k)*utend + dvsfc1(i) = dvsfc1(i)+conw*del(i,k)*vtend enddo enddo do i = 1,im From c7fa4b940f8aa643d9a75f833beb4159138877cb Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Tue, 1 Jun 2021 12:38:23 +0000 Subject: [PATCH 37/52] Updated TKE-EDMF for surface flux output --- physics/satmedmfvdifq.F | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index a59c43e53..206c64456 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -179,11 +179,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & ucko(im,km), vcko(im,km), & buou(im,km), xmf(im,km) ! -! temporary local variables for instantaneous surface fluxes -! (they will be removed later) -! - real(kind=kind_phys) dusfc1(im),dvsfc1(im),dtsfc1(im),dqsfc1(im) -! ! variables for stratocumulus-top induced downdrafts ! real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac), @@ -199,7 +194,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & real(kind=kind_phys) aphi16, aphi5, & wfac, cfac, & gamcrt, gamcrq, sfcfrac, - & conq, cont, conw, +! & conq, cont, conw, & dsdz2, dsdzt, dkmax, & dsig, dt2, dtodsd, & dtodsu, g, factor, dz, @@ -256,10 +251,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & gravi=1.0/grav g=grav gocp=g/cp - cont=cp/g - conq=hvap/g - conw=1.0/g ! for del in pa -! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa +! cont=cp/g +! conq=hvap/g +! conw=1.0/g ! for del in pa +!! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa elocp=hvap/cp el2orc=hvap*hvap/(rv*cp) ! @@ -385,10 +380,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & do i = 1,im z0(i) = 0.01 * zorl(i) rho_a(i) = prsl(i,1)/(rd*t1(i,1)*(1.+fv*max(q1(i,1,1),qmin))) - dusfc1(i) = 0. - dvsfc1(i) = 0. - dtsfc1(i) = 0. - dqsfc1(i) = 0. + dusfc(i) = 0. + dvsfc(i) = 0. + dtsfc(i) = 0. + dqsfc(i) = 0. kpbl(i) = 1 hpbl(i) = 0. kpblx(i) = 1 @@ -1484,8 +1479,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend - dtsfc1(i) = dtsfc1(i)+cont*del(i,k)*ttend - dqsfc1(i) = dqsfc1(i)+conq*del(i,k)*qtend +! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend +! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo do i = 1,im @@ -1623,8 +1618,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & vtend = (f2(i,k)-v1(i,k))*rdt du(i,k) = du(i,k)+utend dv(i,k) = dv(i,k)+vtend - dusfc1(i) = dusfc1(i)+conw*del(i,k)*utend - dvsfc1(i) = dvsfc1(i)+conw*del(i,k)*vtend +! dusfc(i) = dusfc(i)+conw*del(i,k)*utend +! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend enddo enddo do i = 1,im From cd548906c661a156bb6a98ee53ac534e24b38451 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 16:37:48 +0000 Subject: [PATCH 38/52] add zvfun as intent out --- physics/sfc_diff.meta | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta index c0d448655..f079b4357 100644 --- a/physics/sfc_diff.meta +++ b/physics/sfc_diff.meta @@ -96,6 +96,15 @@ kind = kind_phys intent = in optional = F +[zvfun] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction + long_name = function of surface roughness length and green vegetation fraction + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out + optional = F [wind] standard_name = wind_speed_at_lowest_model_layer long_name = wind speed at lowest model level From fce8641a5bcc71c5db6c73fda2a21f3a73ddf5ce Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 16:38:38 +0000 Subject: [PATCH 39/52] define zvfun as intent out --- physics/GFS_surface_generic.meta | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 053613913..145c9f82c 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -1335,7 +1335,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = in optional = F [hffac] standard_name = surface_upward_sensible_heat_flux_reduction_factor From bfcacf6f8a3099daf362b2b1293a8b77c803fde7 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 17:08:02 +0000 Subject: [PATCH 40/52] add zvfun as intent out --- physics/sfc_diff.f | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 9cd6050c4..d797e2176 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -84,6 +84,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & fm10_wat, fm10_lnd, fm10_ice, & !intent(inout) & fh2_wat, fh2_lnd, fh2_ice, & !intent(inout) & ztmax_wat, ztmax_lnd, ztmax_ice, & !intent(inout) + & zvfun, & !intent(out) & errmsg, errflg) !intent(out) ! implicit none @@ -122,6 +123,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & fm10_wat, fm10_lnd, fm10_ice, & & fh2_wat, fh2_lnd, fh2_ice, & & ztmax_wat, ztmax_lnd, ztmax_ice + real(kind=kind_phys), dimension(:), intent(out) :: zvfun ! character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -134,7 +136,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & czilc, tem1, tem2, virtfac ! - real(kind=kind_phys) :: tvs, z0, z0max, ztmax, zvfun, gdx + real(kind=kind_phys) :: tvs, z0, z0max, ztmax, gdx ! real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 ! @@ -189,7 +191,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) thv1 = t1(i) / prslk1(i) * virtfac endif - zvfun = zero + zvfun(i) = zero gdx = sqrt(garea(i)) ! compute stability dependent exchange coefficients @@ -269,15 +271,17 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ztmax_lnd(i) = ztmax_lnd(i) * (10.0_kp**ztpert(i)) endif ztmax_lnd(i) = max(ztmax_lnd(i), zmin) +! +! compute a function of surface roughness & green vegetation fraction (zvfun) ! tem1 = (z0max - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, zero), 1.0_kp) tem2 = max(sigmaf(i), 0.1_kp) - zvfun = sqrt(tem1 * tem2) + zvfun(i) = sqrt(tem1 * tem2) ! call stability ! --- inputs: - & (z1(i), zvfun, gdx, tv1, thv1, wind(i), + & (z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & z0max, ztmax_lnd(i), tvs, grav, thsfc_loc, ! --- outputs: & rb_lnd(i), fm_lnd(i), fh_lnd(i), fm10_lnd(i), fh2_lnd(i), @@ -327,7 +331,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! call stability ! --- inputs: - & (z1(i), zvfun, gdx, tv1, thv1, wind(i), + & (z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & z0max, ztmax_ice(i), tvs, grav, thsfc_loc, ! --- outputs: & rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i), @@ -378,7 +382,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! call stability ! --- inputs: - & (z1(i), zvfun, gdx, tv1, thv1, wind(i), + & (z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & z0max, ztmax_wat(i), tvs, grav, thsfc_loc, ! --- outputs: & rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i), From 70c1e55368f753928c454daccce6c3de232e8947 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 17:20:00 +0000 Subject: [PATCH 41/52] update of GFS_surface_generic_post --- physics/GFS_surface_generic.F90 | 33 +++++++-------------------------- 1 file changed, 7 insertions(+), 26 deletions(-) diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 6b62b2865..b7e3843fb 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -215,8 +215,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl, & dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl, & v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl, & - nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, islmsk, sigmaf, & - runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zorl, hflx, evap, hflxq, zvfun, hffac, errmsg, errflg) + nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, islmsk, & + runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, errmsg, errflg) implicit none @@ -237,15 +237,15 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt evcwa, transa, sbsnoa, snowca, snohfa, ep real(kind=kind_phys), dimension(:), intent(inout) :: runoff, srunoff - real(kind=kind_phys), dimension(:), intent(in) :: drain, runof, sigmaf + real(kind=kind_phys), dimension(:), intent(in) :: drain, runof ! For canopy heat storage logical, intent(in) :: lheatstrg real(kind=kind_phys), intent(in) :: h0facu, h0facs - real(kind=kind_phys), dimension(:), intent(in) :: zorl + real(kind=kind_phys), dimension(:), intent(in) :: zvfun real(kind=kind_phys), dimension(:), intent(in) :: hflx, evap real(kind=kind_phys), dimension(:), intent(out) :: hflxq - real(kind=kind_phys), dimension(:), intent(out) :: zvfun, hffac + real(kind=kind_phys), dimension(:), intent(out) :: hffac ! CCPP error handling variables character(len=*), intent(out) :: errmsg @@ -254,12 +254,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt ! Local variables real(kind=kind_phys), parameter :: albdf = 0.06_kind_phys - ! Parameters for canopy heat storage parametrization - real(kind=kind_phys), parameter :: z0min=0.1, z0max=1.0 - integer :: i real(kind=kind_phys) :: xcosz_loc, ocalnirdf_cpl, ocalnirbm_cpl, ocalvisdf_cpl, ocalvisbm_cpl - real(kind=kind_phys) :: tem, tem1, tem2 ! Initialize CCPP error handling variables errmsg = '' @@ -364,24 +360,9 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt ! in order to achieve heat storage within canopy layer, in the canopy ! heat torage parameterization the kinematic sensible heat flux ! (hflx) as surface boundary forcing to the pbl scheme is -! reduced as a function of surface roughness & green vegetation -! fraction -! -! background diffusivity & background mixing length are also given by -! a function of surface roughness & green vegetation fraction +! reduced in a factor of hffac given as a function of surface roughness & +! green vegetation fraction (zvfun) ! - do i=1,im - if(islmsk(i) == 1) then - tem = 0.01 * zorl(i) ! change unit from cm to m - tem1 = (tem - z0min) / (z0max - z0min) - tem1 = min(max(tem1, 0.0), 1.0) - tem2 = max(sigmaf(i), 0.1) -! tem2 = sigmaf(i) - zvfun(i) = sqrt(tem1 * tem2) - else - zvfun(i) = 0. - endif - enddo do i=1,im hflxq(i) = hflx(i) hffac(i) = 1.0 From d111179b8c0e317353048cd03df03b17920c8119 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 17:22:31 +0000 Subject: [PATCH 42/52] update of GFS_surface_generic.meta --- physics/GFS_surface_generic.meta | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 145c9f82c..e6393f228 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -1221,15 +1221,6 @@ type = integer intent = in optional = F -[sigmaf] - standard_name = bounded_vegetation_area_fraction - long_name = areal fractional cover of green vegetation bounded on the bottom - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [runoff] standard_name = total_runoff long_name = total water runoff @@ -1292,15 +1283,6 @@ kind = kind_phys intent = in optional = F -[zorl] - standard_name = surface_roughness_length - long_name = surface roughness length - units = cm - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [hflx] standard_name = kinematic_surface_upward_sensible_heat_flux long_name = kinematic surface upward sensible heat flux From 44a6088c6e9eb1964e851c7c4c16554c3b118eb4 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 18:17:43 +0000 Subject: [PATCH 43/52] update of GFS_surface_composites.F90 with zvfun --- physics/GFS_surface_composites.F90 | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index ee99e0f85..31becefa4 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -375,14 +375,14 @@ end subroutine GFS_surface_composites_post_finalize !! subroutine GFS_surface_composites_post_run ( & im, kice, km, rd, rvrdm1, cplflx, cplwav2atm, frac_grid, flag_cice, thsfc_loc, islmsk, dry, wet, icy, wind, t1, q1, prsl1, & - landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, & + landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, garea, & cd, cd_wat, cd_lnd, cd_ice, cdq, cdq_wat, cdq_lnd, cdq_ice, rb, rb_wat, rb_lnd, rb_ice, stress, stress_wat, stress_lnd, & stress_ice, ffmm, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh, ffhh_wat, ffhh_lnd, ffhh_ice, uustar, uustar_wat, uustar_lnd, & uustar_ice, fm10, fm10_wat, fm10_lnd, fm10_ice, fh2, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, & cmm, cmm_wat, cmm_lnd, cmm_ice, chh, chh_wat, chh_lnd, chh_ice, gflx, gflx_wat, gflx_lnd, gflx_ice, ep1d, ep1d_wat, & ep1d_lnd, ep1d_ice, weasd, weasd_wat, weasd_lnd, weasd_ice, snowd, snowd_wat, snowd_lnd, snowd_ice, tprcp, tprcp_wat, & tprcp_lnd, tprcp_ice, evap, evap_wat, evap_lnd, evap_ice, hflx, hflx_wat, hflx_lnd, hflx_ice, qss, qss_wat, qss_lnd, & - qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice, stc, & + qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice, sigmaf, zvfun, stc, & grav, prsik1, prslk1, prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice, errmsg, errflg) implicit none @@ -397,13 +397,14 @@ subroutine GFS_surface_composites_post_run ( fm10_wat, fm10_lnd, fm10_ice, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, cmm_wat, cmm_lnd, cmm_ice, & chh_wat, chh_lnd, chh_ice, gflx_wat, gflx_lnd, gflx_ice, ep1d_wat, ep1d_lnd, ep1d_ice, weasd_wat, weasd_lnd, weasd_ice, & snowd_wat, snowd_lnd, snowd_ice,tprcp_wat, tprcp_lnd, tprcp_ice, evap_wat, evap_lnd, evap_ice, hflx_wat, hflx_lnd, & - hflx_ice, qss_wat, qss_lnd, qss_ice, tsfc_wat, tsfc_lnd, tsfc_ice, zorlo, zorll, zorli + hflx_ice, qss_wat, qss_lnd, qss_ice, tsfc_wat, tsfc_lnd, tsfc_ice, zorlo, zorll, zorli, garea real(kind=kind_phys), dimension(:), intent(inout) :: zorl, cd, cdq, rb, stress, ffmm, ffhh, uustar, fm10, & fh2, cmm, chh, gflx, ep1d, weasd, snowd, tprcp, evap, hflx, qss, tsfc, tsfco, tsfcl, tisfc real(kind=kind_phys), dimension(:), intent(in ) :: tice ! interstitial sea ice temperature real(kind=kind_phys), dimension(:), intent(inout) :: hice, cice + real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun real(kind=kind_phys), intent(in ) :: min_seaice real(kind=kind_phys), intent(in ) :: rd, rvrdm1 @@ -424,6 +425,10 @@ subroutine GFS_surface_composites_post_run ( real(kind=kind_phys) :: txl, txi, txo, wfrac, q0, rho ! For calling "stability" real(kind=kind_phys) :: tsurf, virtfac, tv1, thv1, tvs, z0max, ztmax +! + real(kind=kind_phys) :: tem1, tem2, gdx + real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0 +! ! Initialize CCPP error handling variables errmsg = '' @@ -448,6 +453,8 @@ subroutine GFS_surface_composites_post_run ( weasd(i) = txl*weasd_lnd(i) + txi*weasd_ice(i) snowd(i) = txl*snowd_lnd(i) + txi*snowd_ice(i) !tprcp(i) = txl*tprcp_lnd(i) + txi*tprcp_ice(i) + txo*tprcp_wat(i) +! + sigmaf(i) = txl*sigmaf(i) if (.not. flag_cice(i) .and. islmsk(i) == 2) then evap(i) = txl*evap_lnd(i) + wfrac*evap_ice(i) @@ -524,8 +531,17 @@ subroutine GFS_surface_composites_post_run ( stress(i) = stress_ice(i) uustar(i) = uustar_ice(i) else ! Mix of multiple surface types (land, water, and/or ice) - call stability(z1(i), snowd(i), thv1, wind(i), z0max, ztmax, tvs, grav, & ! inputs - tv1, thsfc_loc, & ! inputs +! +! compute zvfun with composite surface roughness & green vegetation fraction +! + tem1 = (z0max - z0lo) / (z0up - z0lo) + tem1 = min(max(tem1, zero), one) + tem2 = max(sigmaf(i), 0.1) + zvfun(i) = sqrt(tem1 * tem2) + gdx = sqrt(garea(i)) +! + call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & ! inputs + z0max, ztmax, tvs, grav, thsfc_loc, & ! inputs rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs stress(i), uustar(i)) endif ! Checking to see if point is one or multiple surface types From 847cf62a8233bfd8f0ee6bc8766b96b63b63483e Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 4 Jun 2021 18:20:55 +0000 Subject: [PATCH 44/52] update of GFS_surface_composites.meta with zvfun, sigmaf, & garea --- physics/GFS_surface_composites.meta | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 95f2c6e4e..88dae1ae4 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -1919,6 +1919,33 @@ kind = kind_phys intent = in optional = F +[garea] + standard_name = cell_area + long_name = area of the grid cell + units = m2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[zvfun] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction + long_name = function of surface roughness length and green vegetation fraction + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F +[sigmaf] + standard_name = bounded_vegetation_area_fraction + long_name = areal fractional cover of green vegetation bounded on the bottom + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F [ztmax_wat] standard_name = bounded_surface_roughness_length_for_heat_over_water long_name = bounded surface roughness length for heat over water From bdfefddd73ee580c7c79d7ded4be81a8e1c7c27d Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Wed, 9 Jun 2021 19:30:57 +0000 Subject: [PATCH 45/52] update sfc_diff with zvfun --- physics/sfc_diff.f | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index d797e2176..0cd373b87 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -290,6 +290,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) if (icy(i)) then ! Some ice + zvfun(i) = zero + if(thsfc_loc) then ! Use local potential temperature tvs = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac else ! Use potential temperature referenced to 1000 hPa @@ -342,6 +344,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! the stuff now put into "stability" if (wet(i)) then ! Some open ocean + + zvfun(i) = zero if(thsfc_loc) then ! Use local potential temperature tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac From 445e3edf0495c63e9f6f8d65bd07889787a73273 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Wed, 9 Jun 2021 19:49:26 +0000 Subject: [PATCH 46/52] update GFS_surface_composites with canopy heat storage variables in the fractional grid --- physics/GFS_surface_composites.F90 | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 31becefa4..973286d93 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -382,13 +382,15 @@ subroutine GFS_surface_composites_post_run ( cmm, cmm_wat, cmm_lnd, cmm_ice, chh, chh_wat, chh_lnd, chh_ice, gflx, gflx_wat, gflx_lnd, gflx_ice, ep1d, ep1d_wat, & ep1d_lnd, ep1d_ice, weasd, weasd_wat, weasd_lnd, weasd_ice, snowd, snowd_wat, snowd_lnd, snowd_ice, tprcp, tprcp_wat, & tprcp_lnd, tprcp_ice, evap, evap_wat, evap_lnd, evap_ice, hflx, hflx_wat, hflx_lnd, hflx_ice, qss, qss_wat, qss_lnd, & - qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice, sigmaf, zvfun, stc, & + qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice, & + sigmaf, zvfun, lheatstrg, h0facu, h0facs, hflxq, hffac, stc, & grav, prsik1, prslk1, prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice, errmsg, errflg) implicit none integer, intent(in) :: im, kice, km logical, intent(in) :: cplflx, frac_grid, cplwav2atm + logical, intent(in) :: lheatstrg logical, dimension(:), intent(in) :: flag_cice, dry, wet, icy integer, dimension(:), intent(in) :: islmsk real(kind=kind_phys), dimension(:), intent(in) :: wind, t1, q1, prsl1, landfrac, lakefrac, oceanfrac, & @@ -404,7 +406,8 @@ subroutine GFS_surface_composites_post_run ( real(kind=kind_phys), dimension(:), intent(in ) :: tice ! interstitial sea ice temperature real(kind=kind_phys), dimension(:), intent(inout) :: hice, cice - real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun + real(kind=kind_phys), dimension(:), intent(inout) :: sigmaf, zvfun, hflxq, hffac + real(kind=kind_phys), intent(in ) :: h0facu, h0facs real(kind=kind_phys), intent(in ) :: min_seaice real(kind=kind_phys), intent(in ) :: rd, rvrdm1 @@ -532,13 +535,33 @@ subroutine GFS_surface_composites_post_run ( uustar(i) = uustar_ice(i) else ! Mix of multiple surface types (land, water, and/or ice) ! -! compute zvfun with composite surface roughness & green vegetation fraction +! re-compute zvfun with composite surface roughness & green vegetation fraction ! tem1 = (z0max - z0lo) / (z0up - z0lo) tem1 = min(max(tem1, zero), one) tem2 = max(sigmaf(i), 0.1) zvfun(i) = sqrt(tem1 * tem2) gdx = sqrt(garea(i)) +! +! re-compute variables for canopy heat storage parameterization with the updated zvfun +! in the fractional grid +! +! + do i=1,im + hflxq(i) = hflx(i) + hffac(i) = 1.0 + enddo + if (lheatstrg) then + do i=1,im + if(hflx(i) > 0.) then + hffac(i) = h0facu * zvfun(i) + else + hffac(i) = h0facs * zvfun(i) + endif + hffac(i) = 1. + hffac(i) + hflxq(i) = hflx(i) / hffac(i) + enddo + endif ! call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & ! inputs z0max, ztmax, tvs, grav, thsfc_loc, & ! inputs From 9e62687333a589621cf3f3a3e46dcb7e18ad54df Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Wed, 9 Jun 2021 19:53:54 +0000 Subject: [PATCH 47/52] update GFS_surface_composites with canopy heat storage variables in the fractional grid --- physics/GFS_surface_composites.meta | 44 +++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 88dae1ae4..73ca42982 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -1937,6 +1937,50 @@ kind = kind_phys intent = inout optional = F +[lheatstrg] + standard_name = flag_for_canopy_heat_storage + long_name = flag for canopy heat storage parameterization + units = flag + dimensions = () + type = logical + intent = in + optional = F +[h0facu] + standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_unstable_surface_layer + long_name = canopy heat storage factor for sensible heat flux in unstable surface layer + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[h0facs] + standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_stable_surface_layer + long_name = canopy heat storage factor for sensible heat flux in stable surface layer + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[hflxq] + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation + units = K m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F +[hffac] + standard_name = surface_upward_sensible_heat_flux_reduction_factor + long_name = surface upward sensible heat flux reduction factor from canopy heat storage + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F [sigmaf] standard_name = bounded_vegetation_area_fraction long_name = areal fractional cover of green vegetation bounded on the bottom From 964fe63b264733b6bccccad1756120c03814d255 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Wed, 9 Jun 2021 21:07:34 +0000 Subject: [PATCH 48/52] replace islmsk=1 with dry --- physics/GFS_surface_generic.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index b7e3843fb..353540b70 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -209,21 +209,21 @@ end subroutine GFS_surface_generic_post_finalize !> \section arg_table_GFS_surface_generic_post_run Argument Table !! \htmlinclude GFS_surface_generic_post_run.html !! - subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,& + subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, dry, icy, wet, & + dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1, & adjsfcdlw, adjsfcdsw, adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfculw, adjsfculw_wat, adjnirbmu, adjnirdfu, & adjvisbmu, adjvisdfu,t2m, q2m, u10m, v10m, tsfc, tsfc_wat, pgr, xcosz, evbs, evcw, trans, sbsno, snowc, snohf, & epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl, & dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl, & v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl, & - nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, islmsk, & + nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep, & runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, errmsg, errflg) implicit none integer, intent(in) :: im - integer, dimension(im), intent(in) :: islmsk logical, intent(in) :: cplflx, cplwav, lssav - logical, dimension(:), intent(in) :: icy, wet + logical, dimension(:), intent(in) :: dry, icy, wet real(kind=kind_phys), intent(in) :: dtf real(kind=kind_phys), dimension(:), intent(in) :: ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1, adjsfcdlw, adjsfcdsw, & @@ -369,7 +369,7 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt enddo if (lheatstrg) then do i=1,im - if(islmsk(i) == 1) then + if (dry(i)) then if(hflx(i) > 0.) then hffac(i) = h0facu * zvfun(i) else From 9da3165bccdbf75a024905688772ce51d3bdab8b Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Wed, 9 Jun 2021 21:11:12 +0000 Subject: [PATCH 49/52] add dry --- physics/GFS_surface_generic.meta | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index e6393f228..fbdb59fcc 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -504,6 +504,14 @@ type = logical intent = in optional = F +[dry] + standard_name = flag_nonzero_land_surface_fraction + long_name = flag indicating presence of some land surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in + optional = F [icy] standard_name = flag_nonzero_sea_ice_surface_fraction long_name = flag indicating presence of some sea ice surface area fraction @@ -1213,14 +1221,6 @@ kind = kind_phys intent = inout optional = F -[islmsk] - standard_name = sea_land_ice_mask - long_name = landmask: sea/land/ice=0/1/2 - units = flag - dimensions = (horizontal_loop_extent) - type = integer - intent = in - optional = F [runoff] standard_name = total_runoff long_name = total water runoff From a7995e39f58775d678112044f4d1af4496e1cfa6 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Thu, 10 Jun 2021 02:27:10 +0000 Subject: [PATCH 50/52] fix a bug --- physics/GFS_surface_composites.F90 | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 973286d93..9feb5f25e 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -545,22 +545,17 @@ subroutine GFS_surface_composites_post_run ( ! ! re-compute variables for canopy heat storage parameterization with the updated zvfun ! in the fractional grid -! ! - do i=1,im - hflxq(i) = hflx(i) - hffac(i) = 1.0 - enddo + hflxq(i) = hflx(i) + hffac(i) = 1.0 if (lheatstrg) then - do i=1,im - if(hflx(i) > 0.) then - hffac(i) = h0facu * zvfun(i) - else - hffac(i) = h0facs * zvfun(i) - endif - hffac(i) = 1. + hffac(i) - hflxq(i) = hflx(i) / hffac(i) - enddo + if(hflx(i) > 0.) then + hffac(i) = h0facu * zvfun(i) + else + hffac(i) = h0facs * zvfun(i) + endif + hffac(i) = 1. + hffac(i) + hflxq(i) = hflx(i) / hffac(i) endif ! call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i), & ! inputs From 686e82a566ecf4b9d69f1ff298ec27b75429305b Mon Sep 17 00:00:00 2001 From: grantfirl Date: Mon, 21 Jun 2021 09:32:38 -0600 Subject: [PATCH 51/52] update scm_sfc_flux_spec.F90/meta to reflect changes in surface fluxes expected by PBL schemes from ccpp-physics PR #665 --- physics/scm_sfc_flux_spec.F90 | 7 +++---- physics/scm_sfc_flux_spec.meta | 13 ++----------- 2 files changed, 5 insertions(+), 15 deletions(-) diff --git a/physics/scm_sfc_flux_spec.F90 b/physics/scm_sfc_flux_spec.F90 index 398bff5d6..e4f425eb2 100644 --- a/physics/scm_sfc_flux_spec.F90 +++ b/physics/scm_sfc_flux_spec.F90 @@ -53,7 +53,7 @@ end subroutine scm_sfc_flux_spec_finalize !! -# Calculate the surface drag coefficient for heat and moisture. !! -# Calculate the u and v wind at 10m. subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec_sh_flux, spec_lh_flux, & - exner_inverse, T_surf, cp, grav, hvap, rd, fvirt, vonKarman, sh_flux, lh_flux, sh_flux_chs, lh_flux_chs, u_star, sfc_stress, cm, ch, & + exner_inverse, T_surf, cp, grav, hvap, rd, fvirt, vonKarman, sh_flux, lh_flux, sh_flux_chs, u_star, sfc_stress, cm, ch, & fm, fh, rb, u10m, v10m, wind1, qss, t2m, q2m, errmsg, errflg) use machine, only: kind_phys @@ -63,7 +63,7 @@ subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec real(kind=kind_phys), intent(in) :: cp, grav, hvap, rd, fvirt, vonKarman real(kind=kind_phys), intent(out) :: sh_flux(:), lh_flux(:), u_star(:), sfc_stress(:), & cm(:), ch(:), fm(:), fh(:), rb(:), u10m(:), v10m(:), wind1(:), qss(:), t2m(:), q2m(:), & - sh_flux_chs(:), lh_flux_chs(:) + sh_flux_chs(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -83,8 +83,7 @@ subroutine scm_sfc_flux_spec_run (u1, v1, z1, t1, q1, p1, roughness_length, spec sh_flux(i) = spec_sh_flux(i) lh_flux(i) = spec_lh_flux(i) sh_flux_chs(i) = sh_flux(i) - lh_flux_chs(i) = lh_flux(i) - + roughness_length_m = 0.01*roughness_length(i) wind1(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) diff --git a/physics/scm_sfc_flux_spec.meta b/physics/scm_sfc_flux_spec.meta index 7edc8c9e5..13187392e 100644 --- a/physics/scm_sfc_flux_spec.meta +++ b/physics/scm_sfc_flux_spec.meta @@ -209,23 +209,14 @@ intent = out optional = F [sh_flux_chs] - standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward sensible heat flux reduced by surface roughness + standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation + long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation units = K m s-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = out optional = F -[lh_flux_chs] - standard_name = kinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughness - long_name = kinematic surface upward latent heat flux reduced by surface roughness - units = kg kg-1 m s-1 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out - optional = F [u_star] standard_name = surface_friction_velocity long_name = boundary layer parameter From 41cbdc25b4f343f7a40176ce4d8c2e1606401c49 Mon Sep 17 00:00:00 2001 From: "Jongil.Han" Date: Fri, 2 Jul 2021 12:29:28 +0000 Subject: [PATCH 52/52] fix the conflict --- physics/GFS_PBL_generic.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index cfa9b2df5..c8c514757 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -317,7 +317,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqsfc_cpl, dusfci_cpl, dvsfci_cpl, dtsfci_cpl, dqsfci_cpl, dusfc_diag, dvsfc_diag, dtsfc_diag, dqsfc_diag, & dusfci_diag, dvsfci_diag, dtsfci_diag, dqsfci_diag, dt3dt, du3dt_PBL, du3dt_OGWD, dv3dt_PBL, dv3dt_OGWD, dq3dt, & dq3dt_ozone, rd, cp, fvirt, hvap, t1, q1, prsl, hflx, ushfsfci, oceanfrac, kdt, dusfc_cice, dvsfc_cice, & - dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, hffac, hefac, & + dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, hffac, & ugrs, vgrs, tgrs, qgrs, save_u, save_v, save_t, save_q, errmsg, errflg) use machine, only : kind_phys