diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index e116f04da..dcd337868 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 @@ -366,7 +366,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, real(kind=kind_phys), dimension(:), intent(out) :: ushfsfci ! 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 @@ -543,7 +543,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 @@ -562,7 +562,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 @@ -606,7 +606,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 diff --git a/physics/GFS_PBL_generic.meta b/physics/GFS_PBL_generic.meta index bbc45c15d..d31d82441 100644 --- a/physics/GFS_PBL_generic.meta +++ b/physics/GFS_PBL_generic.meta @@ -1317,15 +1317,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 diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 39cee585f..8a293fcef 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -1188,7 +1188,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 ) @@ -1231,7 +1230,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 ) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index f3b87d1f5..6429f3f90 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -419,20 +419,22 @@ 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_lnd, weasd_ice, snowd, 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, 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, & @@ -441,13 +443,15 @@ 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_lnd, weasd_ice, & 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, 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 @@ -468,6 +472,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 = '' @@ -490,6 +498,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)) then if (islmsk(i) == 2) then @@ -573,8 +583,32 @@ 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 +! +! 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 +! + hflxq(i) = hflx(i) + hffac(i) = 1.0 + if (lheatstrg) 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 +! + 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 diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 08f66fe2f..17400a982 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -1899,6 +1899,77 @@ 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 +[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 + 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 diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index a4bf92543..b6dd30cfe 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -204,20 +204,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, cplchm, cplwav, lssav, icy, wet, dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,& + subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, 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, & - runoff, srunoff, runof, drain, lheatstrg, z0fac, e0fac, zorl, hflx, evap, hflxq, evapq, hffac, hefac, errmsg, errflg) + runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, errmsg, errflg) implicit none integer, intent(in) :: im logical, intent(in) :: cplflx, cplchm, 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, & @@ -235,11 +236,11 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy, ! For canopy heat storage logical, intent(in) :: lheatstrg - real(kind=kind_phys), intent(in) :: z0fac, e0fac - real(kind=kind_phys), dimension(:), intent(in) :: zorl + real(kind=kind_phys), intent(in) :: h0facu, h0facs + 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, 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) :: hffac ! CCPP error handling variables character(len=*), intent(out) :: errmsg @@ -248,13 +249,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy, ! 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.2, z0max=1.0 - real(kind=kind_phys), parameter :: u10min=2.5, u10max=7.5 - 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 = '' @@ -359,32 +355,28 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy, 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 in a factor of hffac given as a function of surface roughness & +! green vegetation fraction (zvfun) ! do i=1,im hflxq(i) = hflx(i) - evapq(i) = evap(i) - hffac(i) = one - hefac(i) = one + hffac(i) = 1.0 enddo if (lheatstrg) then do i=1,im - tem = 0.01_kind_phys * zorl(i) ! change unit from cm to m - tem1 = (tem - z0min) / (z0max - z0min) - hffac(i) = z0fac * min(max(tem1, zero), one) - tem = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i)) - tem1 = (tem - u10min) / (u10max - u10min) - tem2 = one - min(max(tem1, zero), one) - hffac(i) = tem2 * hffac(i) - hefac(i) = one + e0fac * hffac(i) - hffac(i) = one + hffac(i) - hflxq(i) = hflx(i) / hffac(i) - evapq(i) = evap(i) / hefac(i) + if (dry(i)) 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 diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 893d07dd4..3d7021b03 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -476,6 +476,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 @@ -1229,33 +1237,24 @@ 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 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 @@ -1275,22 +1274,22 @@ 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 - intent = out + intent = in optional = F [hffac] standard_name = surface_upward_sensible_heat_flux_reduction_factor @@ -1301,15 +1300,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 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 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 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) 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)) 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 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 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) diff --git a/physics/moninshoc.meta b/physics/moninshoc.meta index 5d6bebcee..175a2be51 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) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index dd4b93a31..95e99e7be 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -10,9 +10,9 @@ 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 +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 @@ -80,10 +80,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 +99,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 +108,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 +129,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 +138,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 +156,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 +171,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 +198,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 +228,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 +253,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 +265,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 +320,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 @@ -432,16 +436,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 ! @@ -465,6 +470,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 ! @@ -473,12 +479,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 @@ -515,6 +523,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. @@ -599,14 +608,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) @@ -648,8 +675,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)) + @@ -693,14 +720,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 @@ -754,13 +773,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)) @@ -768,6 +886,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 ! @@ -807,13 +927,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 @@ -835,7 +967,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 @@ -882,6 +1015,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 @@ -912,7 +1046,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 @@ -1079,14 +1213,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)) @@ -1162,13 +1296,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)) @@ -1203,21 +1336,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 @@ -1652,60 +1786,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) @@ -1745,7 +1853,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 @@ -1794,7 +1902,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 @@ -1943,7 +2051,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 @@ -2255,7 +2363,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 @@ -2280,7 +2388,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 @@ -2432,40 +2540,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 : @@ -2505,9 +2614,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 !! @@ -2721,10 +2831,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) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index ada8a5cf3..ca7cad4df 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 diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index b0339f8b3..c314809cc 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -14,7 +14,7 @@ subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, & integer, intent(in) :: imfshalcnv integer, intent(in) :: imfshalcnv_samf - + ! CCPP error handling character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -25,7 +25,7 @@ subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, & & ' shallow convection is different from SAMF' errflg = 1 return - end if + end if end subroutine samfshalcnv_init subroutine samfshalcnv_finalize() @@ -61,7 +61,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 +87,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 +108,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 +120,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 +176,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 +197,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 +207,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 +260,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 +268,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)) xmb(i) = 0. scaldfunc(i)=-1.0 ! wang initialized @@ -280,6 +286,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 @@ -287,10 +294,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)) xmb(i) = 0. enddo @@ -344,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 @@ -439,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. @@ -512,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) @@ -563,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)) + @@ -668,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 @@ -737,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 @@ -1008,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 @@ -1351,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? @@ -1523,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 @@ -1558,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 ! @@ -1724,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) diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index f37ec354d..8a346e75b 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 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) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 4bbfe61cc..206c64456 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), @@ -195,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, @@ -212,50 +211,50 @@ 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=3000.) + 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 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) ! @@ -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)*(1.+fv*max(q1(i,1,1),qmin))) 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 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) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 445eb0dc4..0941b1144 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -62,7 +62,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) @@ -72,9 +72,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & thsfc_loc, & !intent(in) & tskin_wat, tskin_lnd, tskin_ice, & !intent(in) & tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in) - & snwdph_lnd,snwdph_ice, & !intent(in) - & z0rl_wat, z0rl_lnd, z0rl_ice, & !intent(inout) - & z0rl_wav, & !intent(inout) + & z0rl_wat, z0rl_lnd, z0rl_ice, & !intent(inout) + & z0rl_wav, & !intent(inout) & ustar_wat, ustar_lnd, ustar_ice, & !intent(inout) & cm_wat, cm_lnd, cm_ice, & !intent(inout) & ch_wat, ch_lnd, ch_ice, & !intent(inout) @@ -85,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 @@ -103,13 +103,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_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) :: & @@ -124,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 @@ -132,17 +132,17 @@ 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) :: tv1 - - real(kind=kind_phys) :: tvs, z0, z0max, snwdph_wat + real(kind=kind_phys) :: tvs, z0, z0max, ztmax, 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) @@ -174,8 +174,6 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! write(0,*)'in sfc_diff, sfc_z0_type=',sfc_z0_type - snwdph_wat = zero - do i=1,im if(flag_iter(i)) then @@ -193,6 +191,9 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) thv1 = t1(i) / prslk1(i) * virtfac endif + zvfun(i) = zero + gdx = sqrt(garea(i)) + ! compute stability dependent exchange coefficients ! this portion of the code is presently suppressed ! @@ -251,24 +252,37 @@ 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_lnd(i) = 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_lnd(i) = 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 = max(min(czilc, 0.8_kp), 0.08_kp) + tem1 = 1.0_kp - sigmaf(i) + czilc = czilc * tem1 * tem1 + ztmax_lnd(i) = 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_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(i) = sqrt(tem1 * tem2) ! call stability ! --- inputs: - & (z1(i), snwdph_lnd(i), thv1, wind(i), - & z0max, ztmax_lnd(i), tvs, grav, tv1, thsfc_loc, + & (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), & cm_lnd(i), ch_lnd(i), stress_lnd(i), ustar_lnd(i)) @@ -276,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 @@ -298,19 +314,27 @@ 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 - - tem1 = 1.0_kp - sigmaf(i) - ztmax_ice(i) = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) +! czilc = 0.8_kp +! +! tem1 = 1.0_kp - sigmaf(i) +! ztmax_ice(i) = z0max*exp( - tem1*tem1 +! & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) +! + czilc = 10.0_kp ** (- 4.0_kp * z0max) + czilc = max(min(czilc, 0.8_kp), 0.08_kp) + tem1 = 1.0_kp - sigmaf(i) + czilc = czilc * tem1 * tem1 + ztmax_ice(i) = z0max * exp( - czilc * ca + & * 258.2_kp * sqrt(ustar_ice(i)*z0max) ) +! ztmax_ice(i) = max(ztmax_ice(i), 1.0e-6) ! call stability ! --- inputs: - & (z1(i), snwdph_ice(i), thv1, wind(i), - & z0max, ztmax_ice(i), tvs, grav, tv1, thsfc_loc, + & (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), & cm_ice(i), ch_ice(i), stress_ice(i), ustar_ice(i)) @@ -320,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 @@ -330,7 +356,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) 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 @@ -360,8 +386,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! call stability ! --- inputs: - & (z1(i), snwdph_wat, thv1, wind(i), - & z0max, ztmax_wat(i), tvs, grav, tv1, thsfc_loc, + & (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), & cm_wat(i), ch_wat(i), stress_wat(i), ustar_wat(i)) @@ -370,7 +396,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) @@ -398,7 +427,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) @@ -420,8 +451,8 @@ end subroutine sfc_diff_run !>\ingroup GFS_diff_main subroutine stability & ! --- inputs: - & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, & - & tv1, thsfc_loc, & + & ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & + & thsfc_loc, & ! --- outputs: & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) !----- @@ -429,8 +460,7 @@ subroutine stability & integer, parameter :: kp = kind_phys ! --- inputs: real(kind=kind_phys), intent(in) :: & - & z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav - real(kind=kind_phys), intent(in) :: tv1 + & z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav logical, intent(in) :: thsfc_loc ! --- outputs: @@ -440,27 +470,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) @@ -483,7 +527,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 ! @@ -502,7 +546,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 @@ -517,11 +561,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) ) @@ -533,7 +575,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 @@ -543,10 +585,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 @@ -556,11 +596,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 diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta index e7551cf99..f079b4357 100644 --- a/physics/sfc_diff.meta +++ b/physics/sfc_diff.meta @@ -87,6 +87,24 @@ 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 = out + optional = F [wind] standard_name = wind_speed_at_lowest_model_layer long_name = wind speed at lowest model level @@ -312,24 +330,6 @@ 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) diff --git a/physics/sfc_nst.f b/physics/sfc_nst.f index 529fa7828..9258b5256 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 ! @@ -190,12 +196,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 @@ -248,6 +257,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 @@ -642,7 +663,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 diff --git a/physics/sfc_nst.meta b/physics/sfc_nst.meta index 685f6f59b..92884faef 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 diff --git a/physics/sfc_ocean.F b/physics/sfc_ocean.F index 67a6df04f..79a9eb295 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 ! @@ -93,15 +98,19 @@ 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) :: 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 +122,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 diff --git a/physics/sfc_ocean.meta b/physics/sfc_ocean.meta index 844eaed88..ace1fcf70 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 diff --git a/physics/shalcnv.meta b/physics/shalcnv.meta index 38436c8bd..68880b773 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) 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) 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)