diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index fe63c1cea..ed26b795f 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -747,9 +747,9 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%det_thl ', Diag%det_thl) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%det_sqv ', Diag%det_sqv) end if - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%nupdraft ', Diag%nupdraft) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%maxwidth ', Diag%maxwidth) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%maxMF ', Diag%maxMF) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%ktop_plume ', Diag%ktop_plume) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%ztop_plume ', Diag%ztop_plume) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%exch_h ', Diag%exch_h) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%exch_m ', Diag%exch_m) end if diff --git a/physics/lsm_ruc.F90 b/physics/lsm_ruc.F90 index 665fe6d14..ba1b1b4e9 100644 --- a/physics/lsm_ruc.F90 +++ b/physics/lsm_ruc.F90 @@ -359,6 +359,8 @@ subroutine lsm_ruc_run & ! inputs & qsurf_ice, gflux_ice, evap_ice, ep1d_ice, hflx_ice, & & cm_ice, ch_ice, snowfallac_ice, acsnow_ice, snowmt_ice, & & albdvis_ice, albdnir_ice, albivis_ice, albinir_ice, & + & add_fire_heat_flux, fire_heat_flux_out, & + & frac_grid_burned_out, & ! --- out & rhosnf, sbsno, & & cmm_lnd, chh_lnd, cmm_ice, chh_ice, & @@ -381,7 +383,7 @@ subroutine lsm_ruc_run & ! inputs real (kind_phys), dimension(:), intent(in) :: oro, sigma real (kind_phys), dimension(:), intent(in) :: & - & t1, sigmaf, laixy, dlwflx, dswsfc, tg3, & + & t1, sigmaf, dlwflx, dswsfc, tg3, & & coszen, prsl1, wind, shdmin, shdmax, & & sfalb_lnd_bck, snoalb, zf, qc, q1, & ! for land @@ -417,7 +419,7 @@ subroutine lsm_ruc_run & ! inputs real (kind_phys), dimension(:), intent(in) :: zs real (kind_phys), dimension(:), intent(in) :: srflag real (kind_phys), dimension(:), intent(inout) :: & - & canopy, trans, smcwlt2, smcref2, & + & canopy, trans, smcwlt2, smcref2, laixy, & ! for land & weasd_lnd, snwdph_lnd, tskin_lnd, & & tsurf_lnd, z0rl_lnd, tsnow_lnd, & @@ -430,6 +432,9 @@ subroutine lsm_ruc_run & ! inputs ! --- in real (kind_phys), dimension(:), intent(in) :: & & rainnc, rainc, ice, snow, graupel, rhonewsn1 + real (kind_phys), dimension(:), intent(in) :: fire_heat_flux_out, & + frac_grid_burned_out + logical, intent(in) :: add_fire_heat_flux ! --- in/out: ! --- on RUC levels real (kind_phys), dimension(:,:), intent(inout) :: & @@ -505,12 +510,13 @@ subroutine lsm_ruc_run & ! inputs & solnet_lnd, sfcexc, & & runoff1, runoff2, acrunoff, semis_bck, & & sfcems_lnd, hfx_lnd, shdfac, shdmin1d, shdmax1d, & + & fire_heat_flux1d, & & sneqv_lnd, snoalb1d_lnd, snowh_lnd, snoh_lnd, tsnav_lnd, & & snomlt_lnd, sncovr_lnd, soilw, soilm, ssoil_lnd, & & soilt_lnd, tbot, & & xlai, swdn, z0_lnd, znt_lnd, rhosnfr, infiltr, & - & precipfr, snfallac_lnd, acsn_lnd, & - & qsfc_lnd, qsg_lnd, qvg_lnd, qcg_lnd, soilt1_lnd, chklowq + & precipfr, snfallac_lnd, acsn_lnd, soilt1_lnd, chklowq, & + & qsfc_lnd, qsg_lnd, qvg_lnd, qcg_lnd, smcwlt, smcref ! ice real (kind_phys),dimension (im,1) :: & & albbck_ice, alb_ice, chs_ice, flhc_ice, flqc_ice, & @@ -540,7 +546,7 @@ subroutine lsm_ruc_run & ! inputs integer :: l, k, i, j, fractional_seaice, ilst real (kind_phys) :: dm, cimin(im) logical :: flag(im), flag_ice(im), flag_ice_uncoupled(im) - logical :: rdlai2d, myj, frpcpn + logical :: myj, frpcpn logical :: debug_print !-- diagnostic point @@ -645,15 +651,27 @@ subroutine lsm_ruc_run & ! inputs nsoil = lsoil_ruc do i = 1, im ! i - horizontal loop - ! reassign smcref2 and smcwlt2 to RUC values if(.not. land(i)) then !water and sea ice - smcref2 (i) = one - smcwlt2 (i) = zero + smcref (i,1) = one + smcwlt (i,1) = zero + xlai (i,1) = zero + elseif (kdt == 1) then + !land + ! reassign smcref2 and smcwlt2 to RUC values at kdt=1 + smcref (i,1) = REFSMC(stype(i)) + smcwlt (i,1) = WLTSMC(stype(i)) + !-- rdlai is .true. when the LAI data is available in the INPUT/sfc_data.nc on cold-start + if(rdlai) then + xlai(i,1) = laixy(i) + else + xlai(i,1) = LAITBL(vtype(i)) + endif else - !land - smcref2 (i) = REFSMC(stype(i)) - smcwlt2 (i) = WLTSMC(stype(i)) + !-- land and kdt > 1, parameters has sub-grid heterogeneity + smcref (i,1) = smcref2 (i) + smcwlt (i,1) = smcwlt2 (i) + xlai (i,1) = laixy (i) endif enddo @@ -813,10 +831,6 @@ subroutine lsm_ruc_run & ! inputs ffrozp(i,j) = real(nint(srflag(i)),kind_phys) endif - !-- rdlai is .false. when the LAI data is not available in the - ! - INPUT/sfc_data.nc - - rdlai2d = rdlai conflx2(i,1,j) = zf(i) * 2._kind_phys ! factor 2. is needed to get the height of ! atm. forcing inside RUC LSM (inherited @@ -843,14 +857,15 @@ subroutine lsm_ruc_run & ! inputs !!\n \a graupelncv - time-step graupel (\f$kg m^{-2} \f$) !!\n \a snowncv - time-step snow (\f$kg m^{-2} \f$) !!\n \a precipfr - time-step precipitation in solid form (\f$kg m^{-2} \f$) -!!\n \a shdfac - areal fractional coverage of green vegetation (0.0-1.0) -!!\n \a shdmin - minimum areal fractional coverage of green vegetation -> !shdmin1d -!!\n \a shdmax - maximum areal fractional coverage of green vegetation -> !shdmax1d +!!\n \a shdfac - areal fractional coverage of green vegetation (0.0-100.%) +!!\n \a shdmin - minimum areal fractional coverage of green vegetation in % -> !shdmin1d +!!\n \a shdmax - maximum areal fractional coverage of green vegetation in % -> !shdmax1d !!\n \a tbot - bottom soil temperature (local yearly-mean sfc air temp) lwdn(i,j) = dlwflx(i) !..downward lw flux at sfc in w/m2 swdn(i,j) = dswsfc(i) !..downward sw flux at sfc in w/m2 + ! all precip input to RUC LSM is in [mm] !prcp(i,j) = rhoh2o * tprcp(i) ! tprcp in [m] - convective plus explicit !raincv(i,j) = rhoh2o * rainc(i) ! total time-step convective precip @@ -918,17 +933,12 @@ subroutine lsm_ruc_run & ! inputs write (0,*)'MODIS landuse is not available' endif - if(rdlai2d) then - xlai(i,j) = laixy(i) - else - xlai(i,j) = zero - endif - semis_bck(i,j) = semisbase(i) ! --- units % shdfac(i,j) = sigmaf(i)*100._kind_phys shdmin1d(i,j) = shdmin(i)*100._kind_phys shdmax1d(i,j) = shdmax(i)*100._kind_phys + fire_heat_flux1d(i,j) = fire_heat_flux_out(i) ! JLS if (land(i)) then ! at least some land in the grid cell @@ -976,7 +986,6 @@ subroutine lsm_ruc_run & ! inputs snoalb1d_lnd(i,j) = snoalb(i) albbck_lnd(i,j) = min(0.9_kind_phys,albbcksol(i)) !sfalb_lnd_bck(i) - !-- spp_lsm if (spp_lsm == 1) then !-- spp for LSM is dimentioned as (1:lsoil_ruc) @@ -999,6 +1008,19 @@ subroutine lsm_ruc_run & ! inputs alb_lnd(i,j) = albbck_lnd(i,j) * (one-sncovr_lnd(i,j)) + snoalb(i) * sncovr_lnd(i,j) ! sfalb_lnd(i) solnet_lnd(i,j) = dswsfc(i)*(one-alb_lnd(i,j)) !..net sw rad flx (dn-up) at sfc in w/m2 + IF ( add_fire_heat_flux .and. fire_heat_flux_out(i) > 0) then ! JLS + if (debug_print) then + print *,'alb_lnd before fire, xlat/xlon ', alb_lnd(i,j), xlat_d(i),xlon_d(i) + print *,'fire_heat_flux_out, frac_grid_burned_out, xlat/xlon ', & + fire_heat_flux_out(i),frac_grid_burned_out(i),xlat_d(i),xlon_d(i) + endif + ! limit albedo in the areas affected by the fire + alb_lnd(i,j) = alb_lnd(i,j) * (one-frac_grid_burned_out(i)) + 0.08_kind_phys*frac_grid_burned_out(i) + if (debug_print) then + print *,'alb_lnd after fire, xlat/xlon ', alb_lnd(i,j), xlat_d(i),xlon_d(i) + endif + ENDIF + cmc(i,j) = canopy(i) ! [mm] soilt_lnd(i,j) = tsurf_lnd(i) ! sanity check for snow temperature tsnow @@ -1163,7 +1185,7 @@ subroutine lsm_ruc_run & ! inputs & wet(i,j), cmc(i,j), shdfac(i,j), alb_lnd(i,j), znt_lnd(i,j), & & z0_lnd(i,j), snoalb1d_lnd(i,j), albbck_lnd(i,j), & & xlai(i,j), landusef(i,:,j), nlcat, & - & soilctop(i,:,j), nscat, & + & soilctop(i,:,j), nscat, smcwlt(i,j), smcref(i,j), & & qsfc_lnd(i,j), qsg_lnd(i,j), qvg_lnd(i,j), qcg_lnd(i,j), & & dew_lnd(i,j), soilt1_lnd(i,j), & & tsnav_lnd(i,j), tbot(i,j), vtype_lnd(i,j), stype_lnd(i,j), & @@ -1178,8 +1200,9 @@ subroutine lsm_ruc_run & ! inputs & infiltr(i,j), runoff1(i,j), runoff2(i,j), acrunoff(i,j), & & sfcexc(i,j), acceta(i,j), ssoil_lnd(i,j), & & snfallac_lnd(i,j), acsn_lnd(i,j), snomlt_lnd(i,j), & - & smfrsoil(i,:,j),keepfrsoil(i,:,j), .false., & - & shdmin1d(i,j), shdmax1d(i,j), rdlai2d, & + & smfrsoil(i,:,j),keepfrsoil(i,:,j), & + & add_fire_heat_flux,fire_heat_flux1d(i,j), .false., & + & shdmin1d(i,j), shdmax1d(i,j), rdlai, & & ims,ime, jms,jme, kms,kme, & & its,ite, jts,jte, kts,kte, errmsg, errflg ) if(debug_print) then @@ -1218,7 +1241,7 @@ subroutine lsm_ruc_run & ! inputs 'ssoil(i,j) =',ssoil_lnd(i,j), & 'snfallac(i,j) =',snfallac_lnd(i,j), & 'acsn_lnd(i,j) =',acsn_lnd(i,j), & - 'snomlt(i,j) =',snomlt_lnd(i,j) + 'snomlt(i,j) =',snomlt_lnd(i,j),'xlai(i,j) =',xlai(i,j) endif endif @@ -1289,6 +1312,10 @@ subroutine lsm_ruc_run & ! inputs ! --- ... unit conversion (from m to mm) snwdph_lnd(i) = snowh_lnd(i,j) * rhoh2o + laixy(i) = xlai(i,j) + smcwlt2(i) = smcwlt(i,j) + smcref2(i) = smcref(i,j) + canopy(i) = cmc(i,j) ! mm weasd_lnd(i) = sneqv_lnd(i,j) ! mm sncovr1_lnd(i) = sncovr_lnd(i,j) @@ -1318,6 +1345,7 @@ subroutine lsm_ruc_run & ! inputs write (0,*)'LAND -i,j,stype_lnd,vtype_lnd',i,j,stype_lnd(i,j),vtype_lnd(i,j) write (0,*)'i,j,tsurf_lnd(i)',i,j,tsurf_lnd(i) write (0,*)'kdt,iter,stsoil(i,:,j)',kdt,iter,stsoil(i,:,j) + write (0,*)'laixy(i)',laixy(i) endif endif ! end of land @@ -1449,7 +1477,7 @@ subroutine lsm_ruc_run & ! inputs & wet_ice(i,j), cmc(i,j), shdfac(i,j), alb_ice(i,j), & & znt_ice(i,j), z0_ice(i,j), snoalb1d_ice(i,j), & & albbck_ice(i,j), xlai(i,j),landusef(i,:,j), nlcat, & - & soilctop(i,:,j), nscat, & + & soilctop(i,:,j), nscat, smcwlt(i,j), smcref(i,j), & & qsfc_ice(i,j), qsg_ice(i,j), qvg_ice(i,j), qcg_ice(i,j), & & dew_ice(i,j), soilt1_ice(i,j), & & tsnav_ice(i,j), tbot(i,j), vtype_ice(i,j), stype_ice(i,j), & @@ -1464,8 +1492,9 @@ subroutine lsm_ruc_run & ! inputs & infiltr(i,j), runoff1(i,j), runoff2(i,j), acrunoff(i,j), & & sfcexc(i,j), acceta(i,j), ssoil_ice(i,j), & & snfallac_ice(i,j), acsn_ice(i,j), snomlt_ice(i,j), & - & smfrice(i,:,j),keepfrice(i,:,j), .false., & - & shdmin1d(i,j), shdmax1d(i,j), rdlai2d, & + & smfrice(i,:,j),keepfrice(i,:,j), & + & add_fire_heat_flux,fire_heat_flux1d(i,j), .false., & + & shdmin1d(i,j), shdmax1d(i,j), rdlai, & & ims,ime, jms,jme, kms,kme, & & its,ite, jts,jte, kts,kte, & & errmsg, errflg) @@ -1502,6 +1531,10 @@ subroutine lsm_ruc_run & ! inputs albivis_ice(i) = sfalb_ice(i) albinir_ice(i) = sfalb_ice(i) + laixy(i) = zero + smcwlt2(i) = zero + smcref2(i) = one + stm(i) = 3.e3_kind_phys ! kg m-2 do k = 1, lsoil_ruc tsice(i,k) = stsice(i,k,j) @@ -1517,6 +1550,7 @@ subroutine lsm_ruc_run & ! inputs write (0,*)'ICE - i,j,stype_ice,vtype_ice)',i,j,stype_ice(i,j),vtype_ice(i,j) write (0,*)'i,j,tsurf_ice(i)',i,j,tsurf_ice(i) write (0,*)'kdt,iter,stsice(i,:,j)',kdt,iter,stsice(i,:,j) + write (0,*)'laixy(i)',laixy(i) endif endif ! ice @@ -1762,6 +1796,8 @@ subroutine rucinit (lsm_cold_start, im, lsoil_ruc, lsoil, & ! in tbot(i,j) = tg3(i) ivgtyp(i,j) = vtype(i) isltyp(i,j) = stype(i) + if(isltyp(i,j)==0) isltyp(i,j)=14 + if(ivgtyp(i,j)==0) ivgtyp(i,j)=17 if (landfrac(i) > zero .or. fice(i) > zero) then !-- land or ice tsk(i,j) = tskin_lnd(i) diff --git a/physics/lsm_ruc.meta b/physics/lsm_ruc.meta index 34a5b8a8b..9bc7fa10a 100644 --- a/physics/lsm_ruc.meta +++ b/physics/lsm_ruc.meta @@ -813,7 +813,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = in + intent = inout [dlwflx] standard_name = surface_downwelling_longwave_flux long_name = surface downwelling longwave flux at current time @@ -1747,6 +1747,29 @@ dimensions = () type = logical intent = in +[add_fire_heat_flux] + standard_name = flag_for_fire_heat_flux + long_name = flag to add fire heat flux to LSM + units = flag + dimensions = () + type = logical + intent = in +[fire_heat_flux_out] + standard_name = surface_fire_heat_flux + long_name = heat flux of fire at the surface + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[frac_grid_burned_out] + standard_name = fraction_of_grid_cell_burning + long_name = ration of the burnt area to the grid cell area + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index ec6b5700d..6840f80bf 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -232,6 +232,18 @@ ! bl_mynn_cloudpdf = 2 (Chab-Becht). ! Removed WRF_CHEM dependencies. ! Many miscellaneous tweaks. +! v4.5.2 / CCPP +! Some code optimization. Removed many conditions from loops. Redesigned the mass- +! flux scheme to use 8 plumes instead of a variable n plumes. This results in +! the removal of the output variable "nudprafts" and adds maxwidth and ztop_plume. +! Revision option bl_mynn_cloudpdf = 2, which now ensures cloud fractions for all +! optically relevant mixing ratios (tip from Greg Thompson). Also, added flexibility +! for tuning near-surface cloud fractions to remove excess fog/low ceilings. +! Now outputs all SGS cloud mixing ratios as grid-mean values, not in-cloud. This +! results in a change in the pre-radiation code to no longer multiply mixing ratios +! by cloud fractions. +! Lots of code cleanup: removal of test code, comments, changing text case, etc. +! Many misc tuning/tweaks. ! ! Many of these changes are now documented in references listed above. !==================================================================== @@ -256,11 +268,11 @@ MODULE module_bl_mynn !=================================================================== ! From here on, these are MYNN-specific parameters: ! The parameters below depend on stability functions of module_sf_mynn. - real(kind_phys), PARAMETER :: cphm_st=5.0, cphm_unst=16.0, & + real(kind_phys), parameter :: cphm_st=5.0, cphm_unst=16.0, & cphh_st=5.0, cphh_unst=16.0 ! Closure constants - real(kind_phys), PARAMETER :: & + real(kind_phys), parameter :: & &pr = 0.74, & &g1 = 0.235, & ! NN2009 = 0.235 &b1 = 24.0, & @@ -275,7 +287,7 @@ MODULE module_bl_mynn &a2 = a1*( g1-c1 )/( g1*pr ), & &g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 ) - real(kind_phys), PARAMETER :: & + real(kind_phys), parameter :: & &cc2 = 1.0-c2, & &cc3 = 1.0-c3, & &e1c = 3.0*a2*b2*cc3, & @@ -286,15 +298,15 @@ MODULE module_bl_mynn ! Constants for min tke in elt integration (qmin), max z/L in els (zmax), ! and factor for eddy viscosity for TKE (Kq = Sqfac*Km): - real(kind_phys), PARAMETER :: qmin=0.0, zmax=1.0, Sqfac=3.0 + real(kind_phys), parameter :: qmin=0.0, zmax=1.0, Sqfac=3.0 ! Note that the following mixing-length constants are now specified in mym_length ! &cns=3.5, alp1=0.23, alp2=0.3, alp3=3.0, alp4=10.0, alp5=0.2 - real(kind_phys), PARAMETER :: gpw=5./3., qcgmin=1.e-8, qkemin=1.e-12 - real(kind_phys), PARAMETER :: tliq = 269. !all hydrometeors are liquid when T > tliq + real(kind_phys), parameter :: gpw=5./3., qcgmin=1.e-8, qkemin=1.e-12 + real(kind_phys), parameter :: tliq = 269. !all hydrometeors are liquid when T > tliq ! Constants for cloud PDF (mym_condensation) - real(kind_phys), PARAMETER :: rr2=0.7071068, rrp=0.3989423 + real(kind_phys), parameter :: rr2=0.7071068, rrp=0.3989423 !>Use Canuto/Kitamura mod (remove Ric and negative TKE) (1:yes, 0:no) !!For more info, see Canuto et al. (2008 JAS) and Kitamura (Journal of the @@ -304,35 +316,35 @@ MODULE module_bl_mynn !!(above) back to NN2009 values (see commented out lines next to the !!parameters above). This only removes the negative TKE problem !!but does not necessarily improve performance - neutral impact. - real(kind_phys), PARAMETER :: CKmod=1. + real(kind_phys), parameter :: CKmod=1. !>Use Ito et al. (2015, BLM) scale-aware (0: no, 1: yes). Note that this also has impacts !!on the cloud PDF and mass-flux scheme, using Honnert et al. (2011) similarity function !!for TKE in the upper PBL/cloud layer. - real(kind_phys), PARAMETER :: scaleaware=1. + real(kind_phys), parameter :: scaleaware=1. !>Of the following the options, use one OR the other, not both. !>Adding top-down diffusion driven by cloud-top radiative cooling - INTEGER, PARAMETER :: bl_mynn_topdown = 0 + integer, parameter :: bl_mynn_topdown = 0 !>Option to activate downdrafts, from Elynn Wu (0: deactive, 1: active) - INTEGER, PARAMETER :: bl_mynn_edmf_dd = 0 + integer, parameter :: bl_mynn_edmf_dd = 0 !>Option to activate heating due to dissipation of TKE (to activate, set to 1.0) - INTEGER, PARAMETER :: dheat_opt = 1 + integer, parameter :: dheat_opt = 1 !Option to activate environmental subsidence in mass-flux scheme - LOGICAL, PARAMETER :: env_subs = .false. + logical, parameter :: env_subs = .false. !Option to switch flux-profile relationship for surface (from Puhales et al. 2020) !0: use original Dyer-Hicks, 1: use Cheng-Brustaert and Blended COARE - INTEGER, PARAMETER :: bl_mynn_stfunc = 1 + integer, parameter :: bl_mynn_stfunc = 1 !option to print out more stuff for debugging purposes - LOGICAL, PARAMETER :: debug_code = .false. - INTEGER, PARAMETER :: idbg = 23 !specific i-point to write out + logical, parameter :: debug_code = .false. + integer, parameter :: idbg = 23 !specific i-point to write out ! Used in WRF-ARW module_physics_init.F - INTEGER :: mynn_level + integer :: mynn_level CONTAINS @@ -388,7 +400,8 @@ SUBROUTINE mynn_bl_driver( & &edmf_thl,edmf_ent,edmf_qc, & &sub_thl3D,sub_sqv3D, & &det_thl3D,det_sqv3D, & - &nupdraft,maxMF,ktop_plume, & + &maxwidth,maxMF,ztop_plume, & + &ktop_plume, & &spp_pbl,pattern_spp_pbl, & &rthraten, & &FLAG_QC,FLAG_QI,FLAG_QNC, & @@ -401,30 +414,30 @@ SUBROUTINE mynn_bl_driver( & !------------------------------------------------------------------- - INTEGER, INTENT(in) :: initflag + integer, intent(in) :: initflag !INPUT NAMELIST OPTIONS: - LOGICAL, INTENT(in) :: restart,cycling - INTEGER, INTENT(in) :: tke_budget - INTEGER, INTENT(in) :: bl_mynn_cloudpdf - INTEGER, INTENT(in) :: bl_mynn_mixlength - INTEGER, INTENT(in) :: bl_mynn_edmf - LOGICAL, INTENT(in) :: bl_mynn_tkeadvect - INTEGER, INTENT(in) :: bl_mynn_edmf_mom - INTEGER, INTENT(in) :: bl_mynn_edmf_tke - INTEGER, INTENT(in) :: bl_mynn_mixscalars - INTEGER, INTENT(in) :: bl_mynn_output - INTEGER, INTENT(in) :: bl_mynn_cloudmix - INTEGER, INTENT(in) :: bl_mynn_mixqt - INTEGER, INTENT(in) :: icloud_bl - real(kind_phys), INTENT(in) :: closure - - LOGICAL, INTENT(in) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,& + logical, intent(in) :: restart,cycling + integer, intent(in) :: tke_budget + integer, intent(in) :: bl_mynn_cloudpdf + integer, intent(in) :: bl_mynn_mixlength + integer, intent(in) :: bl_mynn_edmf + logical, intent(in) :: bl_mynn_tkeadvect + integer, intent(in) :: bl_mynn_edmf_mom + integer, intent(in) :: bl_mynn_edmf_tke + integer, intent(in) :: bl_mynn_mixscalars + integer, intent(in) :: bl_mynn_output + integer, intent(in) :: bl_mynn_cloudmix + integer, intent(in) :: bl_mynn_mixqt + integer, intent(in) :: icloud_bl + real(kind_phys), intent(in) :: closure + + logical, intent(in) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QNC,& FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, & FLAG_OZONE,FLAG_QS - LOGICAL, INTENT(IN) :: mix_chem,enh_mix,rrfs_sd,smoke_dbg + logical, intent(in) :: mix_chem,enh_mix,rrfs_sd,smoke_dbg - INTEGER, INTENT(in) :: & + integer, intent(in) :: & & IDS,IDE,JDS,JDE,KDS,KDE & &,IMS,IME,JMS,JME,KMS,KME & &,ITS,ITE,JTS,JTE,KTS,KTE @@ -444,81 +457,82 @@ SUBROUTINE mynn_bl_driver( & ! to prevent a crash on Cheyenne. Do not change it back without testing if the code runs ! on Cheyenne with the GNU compiler. - real(kind_phys), INTENT(in) :: delt - real(kind_phys), DIMENSION(:), INTENT(in) :: dx - real(kind_phys), DIMENSION(:,:), INTENT(in) :: dz, & + real(kind_phys), intent(in) :: delt + real(kind_phys), dimension(:), intent(in) :: dx + real(kind_phys), dimension(:,:), intent(in) :: dz, & &u,v,w,th,sqv3D,p,exner,rho,T3D - real(kind_phys), DIMENSION(:,:), INTENT(in) :: & + real(kind_phys), dimension(:,:), intent(in) :: & &sqc3D,sqi3D,sqs3D,qni,qnc,qnwfa,qnifa,qnbca - real(kind_phys), DIMENSION(:,:), INTENT(in):: ozone - real(kind_phys), DIMENSION(:), INTENT(in):: ust, & + real(kind_phys), dimension(:,:), intent(in):: ozone + real(kind_phys), dimension(:), intent(in):: ust, & &ch,qsfc,ps,wspd - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: & + real(kind_phys), dimension(:,:), intent(inout) :: & &Qke,Tsq,Qsq,Cov,qke_adv - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: & + real(kind_phys), dimension(:,:), intent(inout) :: & &rublten,rvblten,rthblten,rqvblten,rqcblten, & &rqiblten,rqsblten,rqniblten,rqncblten, & &rqnwfablten,rqnifablten,rqnbcablten - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: dozone - real(kind_phys), DIMENSION(:,:), INTENT(in) :: rthraten + real(kind_phys), dimension(:,:), intent(inout) :: dozone + real(kind_phys), dimension(:,:), intent(in) :: rthraten - real(kind_phys), DIMENSION(:,:), INTENT(out) :: exch_h,exch_m - real(kind_phys), DIMENSION(:), INTENT(in) :: xland, & + real(kind_phys), dimension(:,:), intent(out) :: exch_h,exch_m + real(kind_phys), dimension(:), intent(in) :: xland, & &ts,znt,hfx,qfx,uoce,voce !These 10 arrays are only allocated when bl_mynn_output > 0 - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: & + real(kind_phys), dimension(:,:), intent(inout) :: & & edmf_a,edmf_w,edmf_qt,edmf_thl,edmf_ent,edmf_qc, & & sub_thl3D,sub_sqv3D,det_thl3D,det_sqv3D -! real, DIMENSION(IMS:IME,KMS:KME) :: & +! real, dimension(ims:ime,kms:kme) :: & ! & edmf_a_dd,edmf_w_dd,edmf_qt_dd,edmf_thl_dd,edmf_ent_dd,edmf_qc_dd - real(kind_phys), DIMENSION(:), INTENT(inout) :: Pblh - real(kind_phys), DIMENSION(:), INTENT(inout) :: rmol + real(kind_phys), dimension(:), intent(inout) :: Pblh + real(kind_phys), dimension(:), intent(inout) :: rmol - real(kind_phys), DIMENSION(IMS:IME) :: psig_bl,psig_shcu + real(kind_phys), dimension(ims:ime) :: psig_bl,psig_shcu - INTEGER,DIMENSION(:),INTENT(INOUT) :: & - &KPBL,nupdraft,ktop_plume + integer,dimension(:),intent(INOUT) :: & + &KPBL,ktop_plume - real(kind_phys), DIMENSION(:), INTENT(out) :: maxmf + real(kind_phys), dimension(:), intent(out) :: & + &maxmf,maxwidth,ztop_plume - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: el_pbl + real(kind_phys), dimension(:,:), intent(inout) :: el_pbl - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: & + real(kind_phys), dimension(:,:), intent(inout) :: & &qWT,qSHEAR,qBUOY,qDISS,dqke ! 3D budget arrays are not allocated when tke_budget == 0 ! 1D (local) budget arrays are used for passing between subroutines. - real(kind_phys), DIMENSION(kts:kte) :: & + real(kind_phys), dimension(kts:kte) :: & &qwt1,qshear1,qbuoy1,qdiss1,dqke1,diss_heat - real(kind_phys), DIMENSION(:,:), intent(out) :: Sh3D,Sm3D + real(kind_phys), dimension(:,:), intent(out) :: Sh3D,Sm3D - real(kind_phys), DIMENSION(:,:), INTENT(inout) :: & + real(kind_phys), dimension(:,:), intent(inout) :: & &qc_bl,qi_bl,cldfra_bl - real(kind_phys), DIMENSION(KTS:KTE) :: qc_bl1D,qi_bl1D, & + real(kind_phys), dimension(kts:kte) :: qc_bl1D,qi_bl1D, & &cldfra_bl1D,qc_bl1D_old,qi_bl1D_old,cldfra_bl1D_old ! smoke/chemical arrays - INTEGER, INTENT(IN ) :: nchem, kdvel, ndvel - real(kind_phys), DIMENSION(:,:,:), INTENT(INOUT) :: chem3d - real(kind_phys), DIMENSION(:,:), INTENT(IN) :: vdep - real(kind_phys), DIMENSION(:), INTENT(IN) :: frp,EMIS_ANT_NO + integer, intent(IN ) :: nchem, kdvel, ndvel + real(kind_phys), dimension(:,:,:), intent(INOUT) :: chem3d + real(kind_phys), dimension(:,:), intent(IN) :: vdep + real(kind_phys), dimension(:), intent(IN) :: frp,EMIS_ANT_NO !local - real(kind_phys), DIMENSION(kts:kte ,nchem) :: chem1 - real(kind_phys), DIMENSION(kts:kte+1,nchem) :: s_awchem1 - real(kind_phys), DIMENSION(ndvel) :: vd1 - INTEGER :: ic + real(kind_phys), dimension(kts:kte ,nchem) :: chem1 + real(kind_phys), dimension(kts:kte+1,nchem) :: s_awchem1 + real(kind_phys), dimension(ndvel) :: vd1 + integer :: ic !local vars - INTEGER :: ITF,JTF,KTF, IMD,JMD - INTEGER :: i,j,k,kproblem - real(kind_phys), DIMENSION(KTS:KTE) :: & + integer :: ITF,JTF,KTF, IMD,JMD + integer :: i,j,k,kproblem + real(kind_phys), dimension(kts:kte) :: & &thl,tl,qv1,qc1,qi1,qs1,sqw, & &el, dfm, dfh, dfq, tcd, qcd, pdk, pdt, pdq, pdc, & - &vt, vq, sgm - real(kind_phys), DIMENSION(KTS:KTE) :: & + &vt, vq, sgm, kzero + real(kind_phys), dimension(kts:kte) :: & &thetav,sh,sm,u1,v1,w1,p1, & &ex1,dz1,th1,tk1,rho1,qke1,tsq1,qsq1,cov1, & &sqv,sqi,sqc,sqs, & @@ -527,45 +541,45 @@ SUBROUTINE mynn_bl_driver( & &qnbca1,dqnwfa1,dqnifa1,dqnbca1,dozone1 !mass-flux variables - real(kind_phys), DIMENSION(KTS:KTE) :: & + real(kind_phys), dimension(kts:kte) :: & &dth1mf,dqv1mf,dqc1mf,du1mf,dv1mf - real(kind_phys), DIMENSION(KTS:KTE) :: & + real(kind_phys), dimension(kts:kte) :: & &edmf_a1,edmf_w1,edmf_qt1,edmf_thl1, & &edmf_ent1,edmf_qc1 - real(kind_phys), DIMENSION(KTS:KTE) :: & + real(kind_phys), dimension(kts:kte) :: & &edmf_a_dd1,edmf_w_dd1,edmf_qt_dd1,edmf_thl_dd1, & &edmf_ent_dd1,edmf_qc_dd1 - real(kind_phys), DIMENSION(KTS:KTE) :: & + real(kind_phys), dimension(kts:kte) :: & &sub_thl,sub_sqv,sub_u,sub_v, & &det_thl,det_sqv,det_sqc,det_u,det_v - real(kind_phys), DIMENSION(KTS:KTE+1) :: & + real(kind_phys), dimension(kts:kte+1) :: & &s_aw1,s_awthl1,s_awqt1, & &s_awqv1,s_awqc1,s_awu1,s_awv1,s_awqke1, & &s_awqnc1,s_awqni1,s_awqnwfa1,s_awqnifa1, & &s_awqnbca1 - real(kind_phys), DIMENSION(KTS:KTE+1) :: & + real(kind_phys), dimension(kts:kte+1) :: & &sd_aw1,sd_awthl1,sd_awqt1, & &sd_awqv1,sd_awqc1,sd_awu1,sd_awv1,sd_awqke1 - real(kind_phys), DIMENSION(KTS:KTE+1) :: zw + real(kind_phys), dimension(kts:kte+1) :: zw real(kind_phys) :: cpm,sqcg,flt,fltv,flq,flqv,flqc, & &pmz,phh,exnerg,zet,phi_m, & &afk,abk,ts_decay, qc_bl2, qi_bl2, & - &th_sfc,ztop_plume,wsp + &th_sfc,wsp !top-down diffusion - real(kind_phys), DIMENSION(ITS:ITE) :: maxKHtopdown - real(kind_phys), DIMENSION(KTS:KTE) :: KHtopdown,TKEprodTD + real(kind_phys), dimension(ITS:ITE) :: maxKHtopdown + real(kind_phys), dimension(kts:kte) :: KHtopdown,TKEprodTD - LOGICAL :: INITIALIZE_QKE,problem + logical :: INITIALIZE_QKE,problem ! Stochastic fields - INTEGER, INTENT(IN) :: spp_pbl - real(kind_phys), DIMENSION(:,:), INTENT(IN) :: pattern_spp_pbl - real(kind_phys), DIMENSION(KTS:KTE) :: rstoch_col + integer, intent(IN) :: spp_pbl + real(kind_phys), dimension(:,:), intent(IN) :: pattern_spp_pbl + real(kind_phys), dimension(KTS:KTE) :: rstoch_col ! Substepping TKE - INTEGER :: nsub + integer :: nsub real(kind_phys) :: delt2 @@ -629,9 +643,11 @@ SUBROUTINE mynn_bl_driver( & !edmf_qc_dd(its:ite,kts:kte)=0. ENDIF ktop_plume(its:ite)=0 !int - nupdraft(its:ite)=0 !int + ztop_plume(its:ite)=0. + maxwidth(its:ite)=0. maxmf(its:ite)=0. maxKHtopdown(its:ite)=0. + kzero(kts:kte)=0. ! DH* CHECK HOW MUCH OF THIS INIT IF-BLOCK IS ACTUALLY NEEDED FOR RESTARTS !> - Within the MYNN-EDMF, there is a dependecy check for the first time step, @@ -740,7 +756,7 @@ SUBROUTINE mynn_bl_driver( & !keep snow out for now - increases ceiling bias sqw(k)=sqv(k)+sqc(k)+sqi(k)!+sqs(k) thl(k)=th1(k) - xlvcp/ex1(k)*sqc(k) & - & - xlscp/ex1(k)*(sqi(k)+sqs(k)) + & - xlscp/ex1(k)*(sqi(k))!+sqs(k)) !Use form from Tripoli and Cotton (1981) with their !suggested min temperature to improve accuracy. !thl(k)=th(i,k)*(1.- xlvcp/MAX(tk1(k),TKmin)*sqc(k) & @@ -987,10 +1003,10 @@ SUBROUTINE mynn_bl_driver( & else zw(k)=zw(k-1)+dz(i,k-1) endif - !keep snow out for now - increases ceiling bias + !keep snow out for now - increases ceiling bias sqw(k)= sqv(k)+sqc(k)+sqi(k)!+sqs(k) thl(k)= th1(k) - xlvcp/ex1(k)*sqc(k) & - & - xlscp/ex1(k)*(sqi(k)+sqs(k)) + & - xlscp/ex1(k)*(sqi(k))!+sqs(k)) !Use form from Tripoli and Cotton (1981) with their !suggested min temperature to improve accuracy. !thl(k)=th(i,k)*(1.- xlvcp/MAX(tk1(k),TKmin)*sqc(k) & @@ -1021,7 +1037,7 @@ SUBROUTINE mynn_bl_driver( & endif s_awchem1 = 0.0 -!> - Call get_pblh() to calculate the hybrid \f$\theta_{vli}-TKE\f$ +!> - Call get_pblh() to calculate the hybrid \f$\theta_{v}-TKE\f$ !! PBL height diagnostic. CALL GET_PBLH(KTS,KTE,PBLH(i),thetav,& & Qke1,zw,dz1,xland(i),KPBL(i)) @@ -1147,8 +1163,8 @@ SUBROUTINE mynn_bl_driver( & &FLAG_QNC,FLAG_QNI, & &FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA, & &Psig_shcu(i), & - &nupdraft(i),ktop_plume(i), & - &maxmf(i),ztop_plume, & + &maxwidth(i),ktop_plume(i), & + &maxmf(i),ztop_plume(i), & &spp_pbl,rstoch_col ) endif @@ -1220,9 +1236,9 @@ SUBROUTINE mynn_bl_driver( & call mynn_tendencies(kts,kte,i, & &delt, dz1, rho1, & &u1, v1, th1, tk1, qv1, & - &qc1, qi1, qs1, qnc1, qni1, & + &qc1, qi1, kzero, qnc1, qni1, & !kzero replaces qs1 - not mixing snow &ps(i), p1, ex1, thl, & - &sqv, sqc, sqi, sqs, sqw, & + &sqv, sqc, sqi, kzero, sqw, & !kzero replaces sqs - not mixing snow &qnwfa1, qnifa1, qnbca1, ozone1, & &ust(i),flt,flq,flqv,flqc, & &wspd(i),uoce(i),voce(i), & @@ -1295,27 +1311,27 @@ SUBROUTINE mynn_bl_driver( & &dfm, dfh, dz1, K_m1, K_h1 ) !UPDATE 3D ARRAYS - exch_m(i,:) =k_m1(:) - exch_h(i,:) =k_h1(:) - rublten(i,:) =du1(:) - rvblten(i,:) =dv1(:) - rthblten(i,:)=dth1(:) - rqvblten(i,:)=dqv1(:) + exch_m(i,kts:kte) =k_m1(kts:kte) + exch_h(i,kts:kte) =k_h1(kts:kte) + rublten(i,kts:kte) =du1(kts:kte) + rvblten(i,kts:kte) =dv1(kts:kte) + rthblten(i,kts:kte)=dth1(kts:kte) + rqvblten(i,kts:kte)=dqv1(kts:kte) if (bl_mynn_cloudmix > 0) then - if (flag_qc) rqcblten(i,:)=dqc1(:) - if (flag_qi) rqiblten(i,:)=dqi1(:) - if (flag_qs) rqsblten(i,:)=dqs1(:) + if (flag_qc) rqcblten(i,kts:kte)=dqc1(kts:kte) + if (flag_qi) rqiblten(i,kts:kte)=dqi1(kts:kte) + if (flag_qs) rqsblten(i,kts:kte)=dqs1(kts:kte) else if (flag_qc) rqcblten(i,:)=0. if (flag_qi) rqiblten(i,:)=0. if (flag_qs) rqsblten(i,:)=0. endif if (bl_mynn_cloudmix > 0 .and. bl_mynn_mixscalars > 0) then - if (flag_qnc) rqncblten(i,:) =dqnc1(:) - if (flag_qni) rqniblten(i,:) =dqni1(:) - if (flag_qnwfa) rqnwfablten(i,:)=dqnwfa1(:) - if (flag_qnifa) rqnifablten(i,:)=dqnifa1(:) - if (flag_qnbca) rqnbcablten(i,:)=dqnbca1(:) + if (flag_qnc) rqncblten(i,kts:kte) =dqnc1(kts:kte) + if (flag_qni) rqniblten(i,kts:kte) =dqni1(kts:kte) + if (flag_qnwfa) rqnwfablten(i,kts:kte)=dqnwfa1(kts:kte) + if (flag_qnifa) rqnifablten(i,kts:kte)=dqnifa1(kts:kte) + if (flag_qnbca) rqnbcablten(i,kts:kte)=dqnbca1(kts:kte) else if (flag_qnc) rqncblten(i,:) =0. if (flag_qni) rqniblten(i,:) =0. @@ -1323,19 +1339,19 @@ SUBROUTINE mynn_bl_driver( & if (flag_qnifa) rqnifablten(i,:)=0. if (flag_qnbca) rqnbcablten(i,:)=0. endif - dozone(i,:)=dozone1(:) + dozone(i,kts:kte)=dozone1(kts:kte) if (icloud_bl > 0) then - qc_bl(i,:) =qc_bl1D(:) - qi_bl(i,:) =qi_bl1D(:) - cldfra_bl(i,:)=cldfra_bl1D(:) + qc_bl(i,kts:kte) =qc_bl1D(kts:kte) + qi_bl(i,kts:kte) =qi_bl1D(kts:kte) + cldfra_bl(i,kts:kte)=cldfra_bl1D(kts:kte) endif - el_pbl(i,:)=el(:) - qke(i,:) =qke1(:) - tsq(i,:) =tsq1(:) - qsq(i,:) =qsq1(:) - cov(i,:) =cov1(:) - sh3d(i,:) =sh(:) - sm3d(i,:) =sm(:) + el_pbl(i,kts:kte)=el(kts:kte) + qke(i,kts:kte) =qke1(kts:kte) + tsq(i,kts:kte) =tsq1(kts:kte) + qsq(i,kts:kte) =qsq1(kts:kte) + cov(i,kts:kte) =cov1(kts:kte) + sh3d(i,kts:kte) =sh(kts:kte) + sm3d(i,kts:kte) =sm(kts:kte) if (tke_budget .eq. 1) then !! TKE budget is now given in m**2/s**-3 (Puhales, 2020) @@ -1363,24 +1379,24 @@ SUBROUTINE mynn_bl_driver( & !update updraft/downdraft properties if (bl_mynn_output > 0) then !research mode == 1 if (bl_mynn_edmf > 0) then - edmf_a(i,:) =edmf_a1(:) - edmf_w(i,:) =edmf_w1(:) - edmf_qt(i,:) =edmf_qt1(:) - edmf_thl(i,:) =edmf_thl1(:) - edmf_ent(i,:) =edmf_ent1(:) - edmf_qc(i,:) =edmf_qc1(:) - sub_thl3D(i,:)=sub_thl(:) - sub_sqv3D(i,:)=sub_sqv(:) - det_thl3D(i,:)=det_thl(:) - det_sqv3D(i,:)=det_sqv(:) + edmf_a(i,kts:kte) =edmf_a1(kts:kte) + edmf_w(i,kts:kte) =edmf_w1(kts:kte) + edmf_qt(i,kts:kte) =edmf_qt1(kts:kte) + edmf_thl(i,kts:kte) =edmf_thl1(kts:kte) + edmf_ent(i,kts:kte) =edmf_ent1(kts:kte) + edmf_qc(i,kts:kte) =edmf_qc1(kts:kte) + sub_thl3D(i,kts:kte)=sub_thl(kts:kte) + sub_sqv3D(i,kts:kte)=sub_sqv(kts:kte) + det_thl3D(i,kts:kte)=det_thl(kts:kte) + det_sqv3D(i,kts:kte)=det_sqv(kts:kte) endif !if (bl_mynn_edmf_dd > 0) THEN - ! edmf_a_dd(i,:) =edmf_a_dd1(:) - ! edmf_w_dd(i,:) =edmf_w_dd1(:) - ! edmf_qt_dd(i,:) =edmf_qt_dd1(:) - ! edmf_thl_dd(i,:)=edmf_thl_dd1(:) - ! edmf_ent_dd(i,:)=edmf_ent_dd1(:) - ! edmf_qc_dd(i,:) =edmf_qc_dd1(:) + ! edmf_a_dd(i,kts:kte) =edmf_a_dd1(kts:kte) + ! edmf_w_dd(i,kts:kte) =edmf_w_dd1(kts:kte) + ! edmf_qt_dd(i,kts:kte) =edmf_qt_dd1(kts:kte) + ! edmf_thl_dd(i,kts:kte)=edmf_thl_dd1(kts:kte) + ! edmf_ent_dd(i,kts:kte)=edmf_ent_dd1(kts:kte) + ! edmf_qc_dd(i,kts:kte) =edmf_qc_dd1(kts:kte) !endif endif @@ -1509,27 +1525,27 @@ SUBROUTINE mym_initialize ( & ! !------------------------------------------------------------------- - integer, INTENT(IN) :: kts,kte - integer, INTENT(IN) :: bl_mynn_mixlength - logical, INTENT(IN) :: INITIALIZE_QKE -! real(kind_phys), INTENT(IN) :: ust, rmo, pmz, phh, flt, flq - real(kind_phys), INTENT(IN) :: rmo, Psig_bl, xland - real(kind_phys), INTENT(IN) :: dx, ust, zi - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: dz - real(kind_phys), DIMENSION(kts:kte+1), INTENT(in) :: zw - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: u,v,thl,& + integer, intent(in) :: kts,kte + integer, intent(in) :: bl_mynn_mixlength + logical, intent(in) :: INITIALIZE_QKE +! real(kind_phys), intent(in) :: ust, rmo, pmz, phh, flt, flq + real(kind_phys), intent(in) :: rmo, Psig_bl, xland + real(kind_phys), intent(in) :: dx, ust, zi + real(kind_phys), dimension(kts:kte), intent(in) :: dz + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw + real(kind_phys), dimension(kts:kte), intent(in) :: u,v,thl,& &qw,cldfra_bl1D,edmf_w1,edmf_a1 - real(kind_phys), DIMENSION(kts:kte), INTENT(out) :: tsq,qsq,cov - real(kind_phys), DIMENSION(kts:kte), INTENT(inout) :: el,qke - real(kind_phys), DIMENSION(kts:kte) :: & + real(kind_phys), dimension(kts:kte), intent(out) :: tsq,qsq,cov + real(kind_phys), dimension(kts:kte), intent(inout) :: el,qke + real(kind_phys), dimension(kts:kte) :: & &ql,pdk,pdt,pdq,pdc,dtl,dqw,dtv, & &gm,gh,sm,sh,qkw,vt,vq - INTEGER :: k,l,lmax + integer :: k,l,lmax real(kind_phys):: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1., & &flt=0.,fltv=0.,flq=0.,tmpq - real(kind_phys), DIMENSION(kts:kte) :: theta,thetav - real(kind_phys), DIMENSION(kts:kte) :: rstoch_col - INTEGER ::spp_pbl + real(kind_phys), dimension(kts:kte) :: theta,thetav + real(kind_phys), dimension(kts:kte) :: rstoch_col + integer ::spp_pbl !> - At first ql, vt and vq are set to zero. DO k = kts,kte @@ -1706,17 +1722,17 @@ SUBROUTINE mym_level2 (kts,kte, & ! !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte + integer, intent(in) :: kts,kte #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: dz - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: u,v, & + real(kind_phys), dimension(kts:kte), intent(in) :: dz + real(kind_phys), dimension(kts:kte), intent(in) :: u,v, & &thl,qw,ql,vt,vq,thetav - real(kind_phys), DIMENSION(kts:kte), INTENT(out) :: & + real(kind_phys), dimension(kts:kte), intent(out) :: & &dtl,dqw,dtv,gm,gh,sm,sh integer :: k @@ -1844,25 +1860,25 @@ SUBROUTINE mym_length ( & !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte + integer, intent(in) :: kts,kte #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - INTEGER, INTENT(IN) :: bl_mynn_mixlength - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: dz - real(kind_phys), DIMENSION(kts:kte+1), INTENT(in) :: zw - real(kind_phys), INTENT(in) :: rmo,flt,fltv,flq,Psig_bl,xland - real(kind_phys), INTENT(IN) :: dx,zi - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: u1,v1, & + integer, intent(in) :: bl_mynn_mixlength + real(kind_phys), dimension(kts:kte), intent(in) :: dz + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw + real(kind_phys), intent(in) :: rmo,flt,fltv,flq,Psig_bl,xland + real(kind_phys), intent(in) :: dx,zi + real(kind_phys), dimension(kts:kte), intent(in) :: u1,v1, & &qke,vt,vq,cldfra_bl1D,edmf_w1,edmf_a1 - real(kind_phys), DIMENSION(kts:kte), INTENT(out) :: qkw, el - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: dtv + real(kind_phys), dimension(kts:kte), intent(out) :: qkw, el + real(kind_phys), dimension(kts:kte), intent(in) :: dtv real(kind_phys):: elt,vsc - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: theta - real(kind_phys), DIMENSION(kts:kte) :: qtke,elBLmin,elBLavg,thetaw + real(kind_phys), dimension(kts:kte), intent(in) :: theta + real(kind_phys), dimension(kts:kte) :: qtke,elBLmin,elBLavg,thetaw real(kind_phys):: wt,wt2,zi2,h1,h2,hs,elBLmin0,elBLavg0,cldavg ! THE FOLLOWING CONSTANTS ARE IMPORTANT FOR REGULATING THE @@ -1879,22 +1895,22 @@ SUBROUTINE mym_length ( & !THEY ONLY IMPOSE LIMITS ON THE CALCULATION OF THE MIXING LENGTH !SCALES SO THAT THE BOULAC MIXING LENGTH (IN FREE ATMOS) DOES !NOT ENCROACH UPON THE BOUNDARY LAYER MIXING LENGTH (els, elb & elt). - real(kind_phys), PARAMETER :: minzi = 300. !< min mixed-layer height - real(kind_phys), PARAMETER :: maxdz = 750. !< max (half) transition layer depth + real(kind_phys), parameter :: minzi = 300. !< min mixed-layer height + real(kind_phys), parameter :: maxdz = 750. !< max (half) transition layer depth !! =0.3*2500 m PBLH, so the transition !! layer stops growing for PBLHs > 2.5 km. - real(kind_phys), PARAMETER :: mindz = 300. !< 300 !min (half) transition layer depth + real(kind_phys), parameter :: mindz = 300. !< 300 !min (half) transition layer depth !SURFACE LAYER LENGTH SCALE MODS TO REDUCE IMPACT IN UPPER BOUNDARY LAYER - real(kind_phys), PARAMETER :: ZSLH = 100. !< Max height correlated to surface conditions (m) - real(kind_phys), PARAMETER :: CSL = 2. !< CSL = constant of proportionality to L O(1) + real(kind_phys), parameter :: ZSLH = 100. !< Max height correlated to surface conditions (m) + real(kind_phys), parameter :: CSL = 2. !< CSL = constant of proportionality to L O(1) - INTEGER :: i,j,k + integer :: i,j,k real(kind_phys):: afk,abk,zwk,zwk1,dzk,qdz,vflx,bv,tau_cloud, & & wstar,elb,els,elf,el_stab,el_mf,el_stab_mf,elb_mf, & & PBLH_PLUS_ENT,Uonset,Ugrid,wt_u,el_les - real(kind_phys), PARAMETER :: ctau = 1000. !constant for tau_cloud + real(kind_phys), parameter :: ctau = 1000. !constant for tau_cloud ! tv0 = 0.61*tref ! gtr = 9.81/tref @@ -2254,13 +2270,13 @@ SUBROUTINE boulac_length0(k,kts,kte,zw,dz,qtke,theta,lb1,lb2) ! lb2 = the average of the length up and length down !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: k,kts,kte - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: qtke,dz,theta - real(kind_phys), INTENT(OUT) :: lb1,lb2 - real(kind_phys), DIMENSION(kts:kte+1), INTENT(IN) :: zw + integer, intent(in) :: k,kts,kte + real(kind_phys), dimension(kts:kte), intent(in) :: qtke,dz,theta + real(kind_phys), intent(out) :: lb1,lb2 + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw !LOCAL VARS - INTEGER :: izz, found + integer :: izz, found real(kind_phys):: dlu,dld real(kind_phys):: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz @@ -2404,15 +2420,15 @@ SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2) ! lb2 = the average of the length up and length down !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: qtke,dz,theta - real(kind_phys), DIMENSION(kts:kte), INTENT(OUT):: lb1,lb2 - real(kind_phys), DIMENSION(kts:kte+1), INTENT(IN) :: zw + integer, intent(in) :: kts,kte + real(kind_phys), dimension(kts:kte), intent(in) :: qtke,dz,theta + real(kind_phys), dimension(kts:kte), intent(out):: lb1,lb2 + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw !LOCAL VARS - INTEGER :: iz, izz, found - real(kind_phys), DIMENSION(kts:kte) :: dlu,dld - real(kind_phys), PARAMETER :: Lmax=2000. !soft limit + integer :: iz, izz, found + real(kind_phys), dimension(kts:kte) :: dlu,dld + real(kind_phys), parameter :: Lmax=2000. !soft limit real(kind_phys):: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz !print*,"IN MYNN-BouLac",kts, kte @@ -2618,40 +2634,40 @@ SUBROUTINE mym_turbulence ( & !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte + integer, intent(in) :: kts,kte #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - INTEGER, INTENT(IN) :: bl_mynn_mixlength,tke_budget - real(kind_phys), INTENT(IN) :: closure - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: dz - real(kind_phys), DIMENSION(kts:kte+1), INTENT(in) :: zw - real(kind_phys), INTENT(in) :: rmo,flt,fltv,flq, & + integer, intent(in) :: bl_mynn_mixlength,tke_budget + real(kind_phys), intent(in) :: closure + real(kind_phys), dimension(kts:kte), intent(in) :: dz + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw + real(kind_phys), intent(in) :: rmo,flt,fltv,flq, & &Psig_bl,Psig_shcu,xland,dx,zi - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: u,v,thl,thetav,qw, & + real(kind_phys), dimension(kts:kte), intent(in) :: u,v,thl,thetav,qw, & &ql,vt,vq,qke,tsq,qsq,cov,cldfra_bl1D,edmf_w1,edmf_a1, & &TKEprodTD - real(kind_phys), DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq, & + real(kind_phys), dimension(kts:kte), intent(out) :: dfm,dfh,dfq, & &pdk,pdt,pdq,pdc,tcd,qcd,el - real(kind_phys), DIMENSION(kts:kte), INTENT(inout) :: & + real(kind_phys), dimension(kts:kte), intent(inout) :: & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D real(kind_phys):: q3sq_old,dlsq1,qWTP_old,qWTP_new real(kind_phys):: dudz,dvdz,dTdz,upwp,vpwp,Tpwp - real(kind_phys), DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh + real(kind_phys), dimension(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh - INTEGER :: k + integer :: k ! real(kind_phys):: cc2,cc3,e1c,e2c,e3c,e4c,e5c real(kind_phys):: e6c,dzk,afk,abk,vtt,vqq, & &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh real(kind_phys):: cldavg - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: theta + real(kind_phys), dimension(kts:kte), intent(in) :: theta real(kind_phys):: a2fac, duz, ri !JOE-Canuto/Kitamura mod @@ -2664,10 +2680,10 @@ SUBROUTINE mym_turbulence ( & DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden ! Stochastic - INTEGER, INTENT(IN) :: spp_pbl - real(kind_phys), DIMENSION(KTS:KTE) :: rstoch_col + integer, intent(in) :: spp_pbl + real(kind_phys), dimension(kts:kte) :: rstoch_col real(kind_phys):: Prnum, shb - real(kind_phys), PARAMETER :: Prlimit = 5.0 + real(kind_phys), parameter :: Prlimit = 5.0 ! ! tv0 = 0.61*tref @@ -3042,7 +3058,7 @@ SUBROUTINE mym_turbulence ( & ! q-variance (pdq), and covariance (pdc) pdk(k) = elq*( sm(k)*gm(k) & & +sh(k)*gh(k)+gamv ) + & - & TKEprodTD(k) + & 0.5*TKEprodTD(k) ! xmchen pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k) pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k) pdc(k) = elh*( sh(k)*dtl(k)+gamt ) & @@ -3086,9 +3102,9 @@ SUBROUTINE mym_turbulence ( & !qBUOY1D(k) = elq*(sh(k)*(-dTdz*grav/thl(k)) + gamv) !! ORIGINAL CODE !! Buoyncy term takes the TKEprodTD(k) production now - qBUOY1D(k) = elq*(sh(k)*gh(k)+gamv)+TKEprodTD(k) !staggered + qBUOY1D(k) = elq*(sh(k)*gh(k)+gamv)+0.5*TKEprodTD(k) ! xmchen - !!!Dissipation Term (now it evaluated on mym_predict) + !!!Dissipation Term (now it evaluated in mym_predict) !qDISS1D(k) = (q3sq**(3./2.))/(b1*MAX(el(k),1.)) !! ORIGINAL CODE !! >> EOB @@ -3113,8 +3129,6 @@ SUBROUTINE mym_turbulence ( & qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk ) END DO ! - - if (spp_pbl==1) then DO k = kts,kte dfm(k)= dfm(k) + dfm(k)* rstoch_col(k) * 1.5 * MAX(exp(-MAX(zw(k)-8000.,0.0)/2000.),0.001) @@ -3181,43 +3195,43 @@ SUBROUTINE mym_predict (kts,kte, & & delt, & & dz, & & ust, flt, flq, pmz, phh, & - & el, dfq, rho, & + & el, dfq, rho, & & pdk, pdt, pdq, pdc, & & qke, tsq, qsq, cov, & & s_aw,s_awqke,bl_mynn_edmf_tke, & & qWT1D, qDISS1D,tke_budget) !! TKE budget (Puhales, 2020) !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte + integer, intent(in) :: kts,kte #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - real(kind_phys), INTENT(IN) :: closure - INTEGER, INTENT(IN) :: bl_mynn_edmf_tke,tke_budget - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: dz, dfq, el, rho - real(kind_phys), DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc - real(kind_phys), INTENT(IN) :: flt, flq, pmz, phh - real(kind_phys), INTENT(IN) :: ust, delt - real(kind_phys), DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov + real(kind_phys), intent(in) :: closure + integer, intent(in) :: bl_mynn_edmf_tke,tke_budget + real(kind_phys), dimension(kts:kte), intent(in) :: dz, dfq, el, rho + real(kind_phys), dimension(kts:kte), intent(inout) :: pdk, pdt, pdq, pdc + real(kind_phys), intent(in) :: flt, flq, pmz, phh + real(kind_phys), intent(in) :: ust, delt + real(kind_phys), dimension(kts:kte), intent(inout) :: qke,tsq, qsq, cov ! WA 8/3/15 - real(kind_phys), DIMENSION(kts:kte+1), INTENT(INOUT) :: s_awqke,s_aw + real(kind_phys), dimension(kts:kte+1), intent(inout) :: s_awqke,s_aw !! TKE budget (Puhales, 2020, WRF 4.2.1) << EOB - real(kind_phys), DIMENSION(kts:kte), INTENT(OUT) :: qWT1D, qDISS1D - real(kind_phys), DIMENSION(kts:kte) :: tke_up,dzinv + real(kind_phys), dimension(kts:kte), intent(out) :: qWT1D, qDISS1D + real(kind_phys), dimension(kts:kte) :: tke_up,dzinv !! >> EOB - INTEGER :: k - real(kind_phys), DIMENSION(kts:kte) :: qkw, bp, rp, df3q + integer :: k + real(kind_phys), dimension(kts:kte) :: qkw, bp, rp, df3q real(kind_phys):: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l,onoff - real(kind_phys), DIMENSION(kts:kte) :: dtz - real(kind_phys), DIMENSION(kts:kte) :: a,b,c,d,x + real(kind_phys), dimension(kts:kte) :: dtz + real(kind_phys), dimension(kts:kte) :: a,b,c,d,x - real(kind_phys), DIMENSION(kts:kte) :: rhoinv - real(kind_phys), DIMENSION(kts:kte+1) :: rhoz,kqdz,kmdz + real(kind_phys), dimension(kts:kte) :: rhoinv + real(kind_phys), dimension(kts:kte+1) :: rhoz,kqdz,kmdz ! REGULATE THE MOMENTUM MIXING FROM THE MASS-FLUX SCHEME (on or off) IF (bl_mynn_edmf_tke == 0) THEN @@ -3263,7 +3277,7 @@ SUBROUTINE mym_predict (kts,kte, & kmdz(k) = MAX(kmdz(k), 0.5* s_aw(k)) kmdz(k) = MAX(kmdz(k), -0.5*(s_aw(k)-s_aw(k+1))) ENDDO -!JOE-end conservation mods + !end conservation mods pdk1 = 2.0*ust**3*pmz/( vkz ) phm = 2.0/ust *phh/( vkz ) @@ -3271,8 +3285,8 @@ SUBROUTINE mym_predict (kts,kte, & pdq1 = phm*flq**2 pdc1 = phm*flt*flq ! -! ** pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1. ** - pdk(kts) = pdk1 -pdk(kts+1) +! ** pdk(1)+pdk(2) corresponds to pdk1. ** + pdk(kts) = pdk1 - pdk(kts+1) !! pdt(kts) = pdt1 -pdt(kts+1) !! pdq(kts) = pdq1 -pdq(kts+1) @@ -3367,7 +3381,7 @@ SUBROUTINE mym_predict (kts,kte, & ENDDO k=kte qWT1D(k)=dzinv(k)*(-kqdz(k)*(tke_up(k)-tke_up(k-1)) & - & + 0.5*rhoinv(k)*(-s_aw(k)*tke_up(k)-s_aw(k)*tke_up(k-1)+s_awqke(k))*onoff) !unstaggared + & + 0.5*rhoinv(k)*(-s_aw(k)*tke_up(k)-s_aw(k)*tke_up(k-1)+s_awqke(k))*onoff) !unstaggered !! >> EOBvt qDISS1D=bp*tke_up !! TKE dissipation rate !unstaggered END IF @@ -3596,39 +3610,43 @@ SUBROUTINE mym_condensation (kts,kte, & !------------------------------------------------------------------- - INTEGER, INTENT(IN) :: kts,kte, bl_mynn_cloudpdf + integer, intent(in) :: kts,kte, bl_mynn_cloudpdf #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - real(kind_phys), INTENT(IN) :: HFX1,rmo,xland - real(kind_phys), INTENT(IN) :: dx,pblh1 - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: dz - real(kind_phys), DIMENSION(kts:kte+1), INTENT(IN) :: zw - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: p,exner,thl,qw, & + real(kind_phys), intent(in) :: HFX1,rmo,xland + real(kind_phys), intent(in) :: dx,pblh1 + real(kind_phys), dimension(kts:kte), intent(in) :: dz + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw + real(kind_phys), dimension(kts:kte), intent(in) :: p,exner,thl,qw, & &qv,qc,qi,qs,tsq,qsq,cov,th - real(kind_phys), DIMENSION(kts:kte), INTENT(INOUT) :: vt,vq,sgm + real(kind_phys), dimension(kts:kte), intent(inout) :: vt,vq,sgm - real(kind_phys), DIMENSION(kts:kte) :: alp,a,bet,b,ql,q1,RH - real(kind_phys), DIMENSION(kts:kte), INTENT(OUT) :: qc_bl1D,qi_bl1D, & + real(kind_phys), dimension(kts:kte) :: alp,a,bet,b,ql,q1,RH + real(kind_phys), dimension(kts:kte), intent(out) :: qc_bl1D,qi_bl1D, & &cldfra_bl1D DOUBLE PRECISION :: t3sq, r3sq, c3sq real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, & &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, & - &ls,wt,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & - &qmq,qsat_tk,q1_rh,rh_hack - real(kind_phys), PARAMETER :: rhcrit=0.83 !for hom pdf min sigma - INTEGER :: i,j,k + &ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & + &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc + real(kind_phys), parameter :: qpct_sfc=0.025 + real(kind_phys), parameter :: qpct_pbl=0.030 + real(kind_phys), parameter :: qpct_trp=0.040 + real(kind_phys), parameter :: rhcrit =0.83 !for cloudpdf = 2 + real(kind_phys), parameter :: rhmax =1.01 !for cloudpdf = 2 + integer :: i,j,k real(kind_phys):: erf !VARIABLES FOR ALTERNATIVE SIGMA real:: dth,dtl,dqw,dzk,els - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: Sh,el + real(kind_phys), dimension(kts:kte), intent(in) :: Sh,el !variables for SGS BL clouds real(kind_phys) :: zagl,damp,PBLH2 @@ -3636,11 +3654,11 @@ SUBROUTINE mym_condensation (kts,kte, & !JAYMES: variables for tropopause-height estimation real(kind_phys) :: theta1, theta2, ht1, ht2 - INTEGER :: k_tropo + integer :: k_tropo ! Stochastic - INTEGER, INTENT(IN) :: spp_pbl - real(kind_phys), DIMENSION(KTS:KTE) :: rstoch_col + integer, intent(in) :: spp_pbl + real(kind_phys), dimension(kts:kte) :: rstoch_col real(kind_phys) :: qw_pert ! First, obtain an estimate for the tropopause height (k), using the method employed in the @@ -3794,29 +3812,31 @@ SUBROUTINE mym_condensation (kts,kte, & !Diagnostic statistical scheme of Chaboureau and Bechtold (2002), JAS !but with use of higher-order moments to estimate sigma - PBLH2=MAX(10.,PBLH1) + pblh2=MAX(10._kind_phys,pblh1) zagl = 0. + dzm1 = 0. DO k = kts,kte-1 - zagl = zagl + dz(k) - t = th(k)*exner(k) + zagl = zagl + 0.5*(dz(k) + dzm1) + dzm1 = dz(k) - xl = xl_blend(t) ! obtain latent heat - qsat_tk = qsat_blend(t, p(k)) ! saturation water vapor mixing ratio at tk and p - rh(k)=MAX(MIN(1.00,qw(k)/MAX(1.E-10,qsat_tk)),0.001) + t = th(k)*exner(k) + xl = xl_blend(t) ! obtain latent heat + qsat_tk= qsat_blend(t, p(k)) ! saturation water vapor mixing ratio at tk and p + rh(k) = MAX(MIN(rhmax, qw(k)/MAX(1.E-10,qsat_tk)),0.001_kind_phys) !dqw/dT: Clausius-Clapeyron - dqsl = qsat_tk*ep_2*xlv/( r_d*t**2 ) + dqsl = qsat_tk*ep_2*xlv/( r_d*t**2 ) alp(k) = 1.0/( 1.0+dqsl*xlvcp ) bet(k) = dqsl*exner(k) - rsl = xl*qsat_tk / (r_v*t**2) ! slope of C-C curve at t (=abs temperature) + rsl = xl*qsat_tk / (r_v*t**2) ! slope of C-C curve at t (=abs temperature) ! CB02, Eqn. 4 - cpm = cp + qw(k)*cpv ! CB02, sec. 2, para. 1 - a(k) = 1./(1. + xl*rsl/cpm) ! CB02 variable "a" - b(k) = a(k)*rsl ! CB02 variable "b" + cpm = cp + qw(k)*cpv ! CB02, sec. 2, para. 1 + a(k) = 1./(1. + xl*rsl/cpm) ! CB02 variable "a" + b(k) = a(k)*rsl ! CB02 variable "b" !SPP - qw_pert = qw(k) + qw(k)*0.5*rstoch_col(k)*real(spp_pbl) + qw_pert= qw(k) + qw(k)*0.5*rstoch_col(k)*real(spp_pbl) !This form of qmq (the numerator of Q1) no longer uses the a(k) factor qmq = qw_pert - qsat_tk ! saturation deficit/excess; @@ -3826,28 +3846,46 @@ SUBROUTINE mym_condensation (kts,kte, & r3sq = max( qsq(k), 0.0 ) !Calculate sigma using higher-order moments: sgm(k) = SQRT( r3sq ) - !Set limits on sigma relative to saturation water vapor + !Set constraints on sigma relative to saturation water vapor sgm(k) = min( sgm(k), qsat_tk*0.666 ) - sgm(k) = max( sgm(k), qsat_tk*0.035 ) + !sgm(k) = max( sgm(k), qsat_tk*0.035 ) + + !introduce vertical grid spacing dependence on min sgm + wt = max(500. - max(dz(k)-100.,0.0), 0.0_kind_phys)/500. !=1 for dz < 100 m, =0 for dz > 600 m + sgm(k) = sgm(k) + sgm(k)*0.2*(1.0-wt) !inflate sgm for coarse dz + + !allow min sgm to vary with dz and z. + qpct = qpct_pbl*wt + qpct_trp*(1.0-wt) + qpct = min(qpct, max(qpct_sfc, qpct_pbl*zagl/500.) ) + sgm(k) = max( sgm(k), qsat_tk*qpct ) + q1(k) = qmq / sgm(k) ! Q1, the normalized saturation !Add condition for falling/settling into low-RH layers, so at least - !some cloud fraction is applied for all qc and qi. - rh_hack = rh(k) - !ensure adequate RH & q1 when qi is at least 1e-9 - if (qi(k)>1.e-9) then - rh_hack =min(1.0, rhcrit + 0.06*(9.0 + log10(qi(k)))) + !some cloud fraction is applied for all qc, qs, and qi. + rh_hack= rh(k) + !ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH) + if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then + rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 - q1_rh =-3. + 3.*(rh_hack-rhcrit)/(1.-rhcrit) + q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) q1(k) =max(q1_rh, q1(k) ) endif !ensure adequate RH & q1 when qc is at least 1e-6 if (qc(k)>1.e-6) then - rh_hack =min(1.0, rhcrit + 0.09*(6.0 + log10(qc(k)))) + rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k)))) + rh(k) =max(rh(k), rh_hack) + !add rh-based q1 + q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) + q1(k) =max(q1_rh, q1(k) ) + endif + !ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH) + if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then + rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 - q1_rh =-3. + 3.*(rh_hack-rhcrit)/(1.-rhcrit) + q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) q1(k) =max(q1_rh, q1(k) ) endif @@ -3864,20 +3902,17 @@ SUBROUTINE mym_condensation (kts,kte, & ! Specify hydrometeors ! JAYMES- this option added 8 May 2015 ! The cloud water formulations are taken from CB02, Eq. 8. - IF (q1k < 0.) THEN !unsaturated -#ifdef SINGLE_PREC - ql_water = sgm(k)*EXP(1.2*q1k-1.) -#else - ql_water = sgm(k)*EXP(1.2*q1k-1.) -#endif - ql_ice = sgm(k)*EXP(1.2*q1k-1.) - ELSE IF (q1k > 2.) THEN !supersaturated - ql_water = sgm(k)*q1k - ql_ice = sgm(k)*q1k - ELSE !slightly saturated (0 > q1 < 2) - ql_water = sgm(k)*(EXP(-1.) + 0.66*q1k + 0.086*q1k**2) - ql_ice = sgm(k)*(EXP(-1.) + 0.66*q1k + 0.086*q1k**2) - ENDIF + maxqc = max(qw(k) - qsat_tk, 0.0) + if (q1k < 0.) then !unsaturated + ql_water = sgm(k)*exp(1.2*q1k-1.) + ql_ice = sgm(k)*exp(1.2*q1k-1.) + elseif (q1k > 2.) then !supersaturated + ql_water = min(sgm(k)*q1k, maxqc) + ql_ice = sgm(k)*q1k + else !slightly saturated (0 > q1 < 2) + ql_water = min(sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2), maxqc) + ql_ice = sgm(k)*(exp(-1.) + 0.66*q1k + 0.086*q1k**2) + endif !In saturated grid cells, use average of SGS and resolved values !if ( qc(k) > 1.e-6 ) ql_water = 0.5 * ( ql_water + qc(k) ) @@ -3922,17 +3957,22 @@ SUBROUTINE mym_condensation (kts,kte, & ! Fng = 1.-1.5*q1k !ENDIF ! Use the form of "Fng" from Bechtold and Siebesma (1998, JAS) - IF (q1k .GE. 1.0) THEN + if (q1k .ge. 1.0) then Fng = 1.0 - ELSEIF (q1k .GE. -1.7 .AND. q1k .LT. 1.0) THEN - Fng = EXP(-0.4*(q1k-1.0)) - ELSEIF (q1k .GE. -2.5 .AND. q1k .LT. -1.7) THEN - Fng = 3.0 + EXP(-3.8*(q1k+1.7)) - ELSE - Fng = MIN(23.9 + EXP(-1.6*(q1k+2.5)), 60.) - ENDIF + elseif (q1k .ge. -1.7 .and. q1k .lt. 1.0) then + Fng = exp(-0.4*(q1k-1.0)) + elseif (q1k .ge. -2.5 .and. q1k .lt. -1.7) then + Fng = 3.0 + exp(-3.8*(q1k+1.7)) + else + Fng = min(23.9 + exp(-1.6*(q1k+2.5)), 60._kind_phys) + endif + + cfmax = min(cldfra_bl1D(k), 0.6_kind_phys) + !Further limit the cf going into vt & vq near the surface + zsl = min(max(25., 0.1*pblh2), 100.) + wt = min(zagl/zsl, 1.0) !=0 at z=0 m, =1 above ekman layer + cfmax = cfmax*wt - cfmax= min(cldfra_bl1D(k), 0.6) bb = b(k)*t/th(k) ! bb is "b" in BCMT95. Their "b" differs from ! "b" in CB02 (i.e., b(k) above) by a factor ! of T/theta. Strictly, b(k) above is formulated in @@ -4023,17 +4063,17 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & &bl_mynn_mixscalars ) !------------------------------------------------------------------- - INTEGER, INTENT(in) :: kts,kte,i + integer, intent(in) :: kts,kte,i #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - INTEGER, INTENT(in) :: bl_mynn_cloudmix,bl_mynn_mixqt, & + integer, intent(in) :: bl_mynn_cloudmix,bl_mynn_mixqt, & bl_mynn_edmf,bl_mynn_edmf_mom, & bl_mynn_mixscalars - LOGICAL, INTENT(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QS, & + logical, intent(IN) :: FLAG_QI,FLAG_QNI,FLAG_QC,FLAG_QS, & &FLAG_QNC,FLAG_QNWFA,FLAG_QNIFA,FLAG_QNBCA ! thl - liquid water potential temperature @@ -4043,47 +4083,47 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & ! flq - surface flux of qw ! mass-flux plumes - real(kind_phys), DIMENSION(kts:kte+1), INTENT(in) :: s_aw, & + real(kind_phys), dimension(kts:kte+1), intent(in) :: s_aw, & &s_awthl,s_awqt,s_awqnc,s_awqni,s_awqv,s_awqc,s_awu,s_awv, & &s_awqnwfa,s_awqnifa,s_awqnbca, & &sd_aw,sd_awthl,sd_awqt,sd_awqv,sd_awqc,sd_awu,sd_awv ! tendencies from mass-flux environmental subsidence and detrainment - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: sub_thl,sub_sqv, & + real(kind_phys), dimension(kts:kte), intent(in) :: sub_thl,sub_sqv, & &sub_u,sub_v,det_thl,det_sqv,det_sqc,det_u,det_v - real(kind_phys), DIMENSION(kts:kte), INTENT(in) :: u,v,th,tk,qv,qc,qi,& + real(kind_phys), dimension(kts:kte), intent(in) :: u,v,th,tk,qv,qc,qi,& &qs,qni,qnc,rho,p,exner,dfq,dz,tsq,qsq,cov,tcd,qcd, & &cldfra_bl1d,diss_heat - real(kind_phys), DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc,& + real(kind_phys), dimension(kts:kte), intent(inout) :: thl,sqw,sqv,sqc,& &sqi,sqs,qnwfa,qnifa,qnbca,ozone,dfm,dfh - real(kind_phys), DIMENSION(kts:kte), INTENT(inout) :: du,dv,dth,dqv, & + real(kind_phys), dimension(kts:kte), intent(inout) :: du,dv,dth,dqv, & &dqc,dqi,dqs,dqni,dqnc,dqnwfa,dqnifa,dqnbca,dozone - real(kind_phys), INTENT(IN) :: flt,flq,flqv,flqc,uoce,voce - real(kind_phys), INTENT(IN) :: ust,delt,psfc,wspd + real(kind_phys), intent(in) :: flt,flq,flqv,flqc,uoce,voce + real(kind_phys), intent(in) :: ust,delt,psfc,wspd !debugging real(kind_phys):: wsp,wsp2,tk2,th2 - LOGICAL :: problem + logical :: problem integer :: kproblem -! real(kind_phys), INTENT(IN) :: gradu_top,gradv_top,gradth_top,gradqv_top +! real(kind_phys), intent(in) :: gradu_top,gradv_top,gradth_top,gradqv_top !local vars - real(kind_phys), DIMENSION(kts:kte) :: dtz,dfhc,dfmc,delp - real(kind_phys), DIMENSION(kts:kte) :: sqv2,sqc2,sqi2,sqs2,sqw2, & + real(kind_phys), dimension(kts:kte) :: dtz,dfhc,dfmc,delp + real(kind_phys), dimension(kts:kte) :: sqv2,sqc2,sqi2,sqs2,sqw2, & &qni2,qnc2,qnwfa2,qnifa2,qnbca2,ozone2 - real(kind_phys), DIMENSION(kts:kte) :: zfac,plumeKh,rhoinv - real(kind_phys), DIMENSION(kts:kte) :: a,b,c,d,x - real(kind_phys), DIMENSION(kts:kte+1) :: rhoz, & !rho on model interface + real(kind_phys), dimension(kts:kte) :: zfac,plumeKh,rhoinv + real(kind_phys), dimension(kts:kte) :: a,b,c,d,x + real(kind_phys), dimension(kts:kte+1) :: rhoz, & !rho on model interface &khdz,kmdz real(kind_phys):: rhs,gfluxm,gfluxp,dztop,maxdfh,mindfh,maxcf,maxKh,zw real(kind_phys):: t,esat,qsl,onoff,kh,km,dzk,rhosfc real(kind_phys):: ustdrag,ustdiff,qvflux real(kind_phys):: th_new,portion_qc,portion_qi,condensate,qsat - INTEGER :: k,kk + integer :: k,kk !Activate nonlocal mixing from the mass-flux scheme for !number concentrations and aerosols (0.0 = no; 1.0 = yes) - real(kind_phys), PARAMETER :: nonloc = 1.0 + real(kind_phys), parameter :: nonloc = 1.0 dztop=.5*(dz(kte)+dz(kte-1)) @@ -4586,7 +4626,8 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & !============================================ ! MIX SNOW ( sqs ) !============================================ -IF (bl_mynn_cloudmix > 0 .AND. FLAG_QS) THEN +!hard-code to not mix snow +IF (bl_mynn_cloudmix > 0 .AND. .false.) THEN k=kts !rho-weighted: @@ -4813,8 +4854,8 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & d(kte)=qnbca(kte) ! CALL tridiag(kte,a,b,c,d) -! CALL tridiag2(kte,a,b,c,d,x) - CALL tridiag3(kte,a,b,c,d,x) + CALL tridiag2(kte,a,b,c,d,x) +! CALL tridiag3(kte,a,b,c,d,x) DO k=kts,kte !qnbca2(k)=d(k-kts+1) @@ -4891,9 +4932,6 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & sqi2(k) = 0.0 ! if sqw2 > qsat sqc2(k) = 0.0 ENDIF - !dqv(k) = (sqv2(k) - sqv(k))/delt - !dqc(k) = (sqc2(k) - sqc(k))/delt - !dqi(k) = (sqi2(k) - sqi(k))/delt ENDDO ENDIF @@ -4902,7 +4940,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & ! WATER VAPOR TENDENCY !===================== DO k=kts,kte - Dqv(k)=(sqv2(k)/(1.-sqv2(k)) - qv(k))/delt + Dqv(k)=(sqv2(k) - sqv(k))/delt !if (sqv2(k) < 0.0)print*,"neg qv:",sqv2(k),k ENDDO @@ -4913,7 +4951,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & !print*,"FLAG_QC:",FLAG_QC IF (FLAG_QC) THEN DO k=kts,kte - Dqc(k)=(sqc2(k)/(1.-sqv2(k)) - qc(k))/delt + Dqc(k)=(sqc2(k) - sqc(k))/delt !if (sqc2(k) < 0.0)print*,"neg qc:",sqc2(k),k ENDDO ELSE @@ -4941,7 +4979,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & !=================== IF (FLAG_QI) THEN DO k=kts,kte - Dqi(k)=(sqi2(k)/(1.-sqv2(k)) - qi(k))/delt + Dqi(k)=(sqi2(k) - sqi(k))/delt !if (sqi2(k) < 0.0)print*,"neg qi:",sqi2(k),k ENDDO ELSE @@ -4953,9 +4991,9 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & !=================== ! CLOUD SNOW TENDENCY !=================== - IF (FLAG_QS) THEN + IF (.false.) THEN !disabled DO k=kts,kte - Dqs(k)=(sqs2(k)/(1.-sqs2(k)) - qs(k))/delt + Dqs(k)=(sqs2(k) - sqs(k))/delt ENDDO ELSE DO k=kts,kte @@ -4979,10 +5017,11 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & ELSE !-MIX CLOUD SPECIES? !CLOUDS ARE NOT NIXED (when bl_mynn_cloudmix == 0) DO k=kts,kte - Dqc(k)=0. + Dqc(k) =0. Dqnc(k)=0. - Dqi(k)=0. + Dqi(k) =0. Dqni(k)=0. + Dqs(k) =0. ENDDO ENDIF @@ -5207,36 +5246,36 @@ SUBROUTINE mynn_mix_chem(kts,kte,i, & enh_mix, smoke_dbg ) !------------------------------------------------------------------- - INTEGER, INTENT(in) :: kts,kte,i - real(kind_phys), DIMENSION(kts:kte), INTENT(IN) :: dfh,dz,tcd,qcd - real(kind_phys), DIMENSION(kts:kte), INTENT(INOUT) :: rho - real(kind_phys), INTENT(IN) :: flt - real(kind_phys), INTENT(IN) :: delt,pblh - INTEGER, INTENT(IN) :: nchem, kdvel, ndvel - real(kind_phys), DIMENSION( kts:kte+1), INTENT(IN) :: s_aw - real(kind_phys), DIMENSION( kts:kte, nchem ), INTENT(INOUT) :: chem1 - real(kind_phys), DIMENSION( kts:kte+1,nchem), INTENT(IN) :: s_awchem - real(kind_phys), DIMENSION( ndvel ), INTENT(IN) :: vd1 - real(kind_phys), INTENT(IN) :: emis_ant_no,frp - LOGICAL, INTENT(IN) :: rrfs_sd,enh_mix,smoke_dbg + integer, intent(in) :: kts,kte,i + real(kind_phys), dimension(kts:kte), intent(in) :: dfh,dz,tcd,qcd + real(kind_phys), dimension(kts:kte), intent(inout) :: rho + real(kind_phys), intent(in) :: flt + real(kind_phys), intent(in) :: delt,pblh + integer, intent(in) :: nchem, kdvel, ndvel + real(kind_phys), dimension( kts:kte+1), intent(in) :: s_aw + real(kind_phys), dimension( kts:kte, nchem ), intent(inout) :: chem1 + real(kind_phys), dimension( kts:kte+1,nchem), intent(in) :: s_awchem + real(kind_phys), dimension( ndvel ), intent(in) :: vd1 + real(kind_phys), intent(in) :: emis_ant_no,frp + logical, intent(in) :: rrfs_sd,enh_mix,smoke_dbg !local vars - real(kind_phys), DIMENSION(kts:kte) :: dtz - real(kind_phys), DIMENSION(kts:kte) :: a,b,c,d,x + real(kind_phys), dimension(kts:kte) :: dtz + real(kind_phys), dimension(kts:kte) :: a,b,c,d,x real(kind_phys):: rhs,dztop real(kind_phys):: t,dzk real(kind_phys):: hght real(kind_phys):: khdz_old, khdz_back - INTEGER :: k,kk,kmaxfire ! JLS 12/21/21 - INTEGER :: ic ! Chemical array loop index + integer :: k,kk,kmaxfire ! JLS 12/21/21 + integer :: ic ! Chemical array loop index - INTEGER, SAVE :: icall + integer, SAVE :: icall - real(kind_phys), DIMENSION(kts:kte) :: rhoinv - real(kind_phys), DIMENSION(kts:kte+1) :: rhoz,khdz - real(kind_phys), PARAMETER :: NO_threshold = 10.0 ! For anthropogenic sources - real(kind_phys), PARAMETER :: frp_threshold = 10.0 ! RAR 02/11/22: I increased the frp threshold to enhance mixing over big fires - real(kind_phys), PARAMETER :: pblh_threshold = 100.0 + real(kind_phys), dimension(kts:kte) :: rhoinv + real(kind_phys), dimension(kts:kte+1) :: rhoz,khdz + real(kind_phys), parameter :: NO_threshold = 10.0 ! For anthropogenic sources + real(kind_phys), parameter :: frp_threshold = 10.0 ! RAR 02/11/22: I increased the frp threshold to enhance mixing over big fires + real(kind_phys), parameter :: pblh_threshold = 100.0 dztop=.5*(dz(kte)+dz(kte-1)) @@ -5335,14 +5374,14 @@ SUBROUTINE retrieve_exchange_coeffs(kts,kte,& !------------------------------------------------------------------- - INTEGER , INTENT(in) :: kts,kte + integer , intent(in) :: kts,kte - real(kind_phys), DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh + real(kind_phys), dimension(KtS:KtE), intent(in) :: dz,dfm,dfh - real(kind_phys), DIMENSION(KtS:KtE), INTENT(out) :: K_m, K_h + real(kind_phys), dimension(KtS:KtE), intent(out) :: K_m, K_h - INTEGER :: k + integer :: k real(kind_phys):: dzk K_m(kts)=0. @@ -5368,13 +5407,13 @@ SUBROUTINE tridiag(n,a,b,c,d) !------------------------------------------------------------------- - INTEGER, INTENT(in):: n - real(kind_phys), DIMENSION(n), INTENT(in) :: a,b - real(kind_phys), DIMENSION(n), INTENT(inout) :: c,d + integer, intent(in):: n + real(kind_phys), dimension(n), intent(in) :: a,b + real(kind_phys), dimension(n), intent(inout) :: c,d - INTEGER :: i + integer :: i real(kind_phys):: p - real(kind_phys), DIMENSION(n) :: q + real(kind_phys), dimension(n) :: q c(n)=0. q(1)=-c(1)/b(1) @@ -5508,23 +5547,23 @@ SUBROUTINE GET_PBLH(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea,kzi) !value could be found to work best in all conditions. !--------------------------------------------------------------- - INTEGER,INTENT(IN) :: KTS,KTE + integer,intent(in) :: KTS,KTE #ifdef HARDCODE_VERTICAL # define kts 1 # define kte HARDCODE_VERTICAL #endif - real(kind_phys), INTENT(OUT) :: zi - real(kind_phys), INTENT(IN) :: landsea - real(kind_phys), DIMENSION(KTS:KTE), INTENT(IN) :: thetav1D, qke1D, dz1D - real(kind_phys), DIMENSION(KTS:KTE+1), INTENT(IN) :: zw1D + real(kind_phys), intent(out) :: zi + real(kind_phys), intent(in) :: landsea + real(kind_phys), dimension(kts:kte), intent(in) :: thetav1D, qke1D, dz1D + real(kind_phys), dimension(kts:kte+1), intent(in) :: zw1D !LOCAL VARS real(kind_phys):: PBLH_TKE,qtke,qtkem1,wt,maxqke,TKEeps,minthv real(kind_phys):: delt_thv !delta theta-v; dependent on land/sea point - real(kind_phys), PARAMETER :: sbl_lim = 200. !upper limit of stable BL height (m). - real(kind_phys), PARAMETER :: sbl_damp = 400. !transition length for blending (m). - INTEGER :: I,J,K,kthv,ktke,kzi + real(kind_phys), parameter :: sbl_lim = 200. !upper limit of stable BL height (m). + real(kind_phys), parameter :: sbl_damp = 400. !transition length for blending (m). + integer :: I,J,K,kthv,ktke,kzi !Initialize KPBL (kzi) kzi = 2 @@ -5689,12 +5728,12 @@ SUBROUTINE DMP_mf( & & F_QNWFA,F_QNIFA,F_QNBCA, & & Psig_shcu, & ! output info - & nup2,ktop,maxmf,ztop, & + & maxwidth,ktop,maxmf,ztop, & ! inputs for stochastic perturbations & spp_pbl,rstoch_col ) ! inputs: - INTEGER, INTENT(IN) :: KTS,KTE,KPBL,momentum_opt,tke_opt,scalar_opt + integer, intent(in) :: KTS,KTE,KPBL,momentum_opt,tke_opt,scalar_opt #ifdef HARDCODE_VERTICAL # define kts 1 @@ -5702,133 +5741,137 @@ SUBROUTINE DMP_mf( & #endif ! Stochastic - INTEGER, INTENT(IN) :: spp_pbl - real(kind_phys), DIMENSION(KTS:KTE) :: rstoch_col + integer, intent(in) :: spp_pbl + real(kind_phys), dimension(kts:kte) :: rstoch_col - real(kind_phys),DIMENSION(KTS:KTE), INTENT(IN) :: & + real(kind_phys),dimension(kts:kte), intent(in) :: & &U,V,W,TH,THL,TK,QT,QV,QC, & &exner,dz,THV,P,rho,qke,qnc,qni,qnwfa,qnifa,qnbca - real(kind_phys),DIMENSION(KTS:KTE+1), INTENT(IN) :: zw !height at full-sigma - real(kind_phys), INTENT(IN) :: flt,fltv,flq,flqv,Psig_shcu, & + real(kind_phys),dimension(kts:kte+1), intent(in) :: zw !height at full-sigma + real(kind_phys), intent(in) :: flt,fltv,flq,flqv,Psig_shcu, & &landsea,ts,dx,dt,ust,pblh - LOGICAL, OPTIONAL :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA,F_QNBCA + logical, optional :: F_QC,F_QI,F_QNC,F_QNI,F_QNWFA,F_QNIFA,F_QNBCA ! outputs - updraft properties - real(kind_phys),DIMENSION(KTS:KTE), INTENT(OUT) :: edmf_a,edmf_w, & + real(kind_phys),dimension(kts:kte), intent(out) :: edmf_a,edmf_w, & & edmf_qt,edmf_thl,edmf_ent,edmf_qc !add one local edmf variable: - real(kind_phys),DIMENSION(KTS:KTE) :: edmf_th + real(kind_phys),dimension(kts:kte) :: edmf_th ! output - INTEGER, INTENT(OUT) :: nup2,ktop - real(kind_phys), INTENT(OUT) :: maxmf - real(kind_phys), INTENT(OUT) :: ztop + integer, intent(out) :: ktop + real(kind_phys), intent(out) :: maxmf,ztop,maxwidth ! outputs - variables needed for solver - real(kind_phys),DIMENSION(KTS:KTE+1) :: s_aw, & !sum ai*rho*wis_awphi + real(kind_phys),dimension(kts:kte+1) :: s_aw, & !sum ai*rho*wis_awphi &s_awthl,s_awqt,s_awqv,s_awqc,s_awqnc,s_awqni, & &s_awqnwfa,s_awqnifa,s_awqnbca,s_awu,s_awv, & &s_awqke,s_aw2 - real(kind_phys),DIMENSION(KTS:KTE), INTENT(INOUT) :: & + real(kind_phys),dimension(kts:kte), intent(inout) :: & &qc_bl1d,cldfra_bl1d,qc_bl1d_old,cldfra_bl1d_old - INTEGER, PARAMETER :: nup=10, debug_mf=0 + integer, parameter :: nup=8, debug_mf=0 + real(kind_phys) :: nup2 !------------- local variables ------------------- ! updraft properties defined on interfaces (k=1 is the top of the ! first model layer - real(kind_phys),DIMENSION(KTS:KTE+1,1:NUP) :: & + real(kind_phys),dimension(kts:kte+1,1:NUP) :: & &UPW,UPTHL,UPQT,UPQC,UPQV, & &UPA,UPU,UPV,UPTHV,UPQKE,UPQNC, & &UPQNI,UPQNWFA,UPQNIFA,UPQNBCA ! entrainment variables - real(kind_phys),DIMENSION(KTS:KTE,1:NUP) :: ENT,ENTf - INTEGER,DIMENSION(KTS:KTE,1:NUP) :: ENTi + real(kind_phys),dimension(kts:kte,1:NUP) :: ENT,ENTf + integer,dimension(kts:kte,1:NUP) :: ENTi ! internal variables - INTEGER :: K,I,k50 + integer :: K,I,k50 real(kind_phys):: fltv2,wstar,qstar,thstar,sigmaW,sigmaQT, & &sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,Psig_w,maxw,maxqc,wpbl real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,QNCn,QNIn, & - QNWFAn,QNIFAn,QNBCAn, & - Wn2,Wn,EntEXP,EntEXM,EntW,BCOEFF,THVkm1,THVk,Pk,rho_int + & QNWFAn,QNIFAn,QNBCAn, & + & Wn2,Wn,EntEXP,EntEXM,EntW,BCOEFF,THVkm1,THVk,Pk,rho_int ! w parameters - real(kind_phys), PARAMETER :: & + real(kind_phys), parameter :: & &Wa=2./3., & &Wb=0.002, & &Wc=1.5 ! Lateral entrainment parameters ( L0=100 and ENT0=0.1) were taken from ! Suselj et al (2013, jas). Note that Suselj et al (2014,waf) use L0=200 and ENT0=0.2. - real(kind_phys),PARAMETER :: & - & L0=100., & - & ENT0=0.1 - - ! Implement ideas from Neggers (2016, JAMES): - real(kind_phys), PARAMETER :: Atot = 0.10 ! Maximum total fractional area of all updrafts - real(kind_phys), PARAMETER :: lmax = 1000.! diameter of largest plume - real(kind_phys), PARAMETER :: dl = 100. ! diff size of each plume - the differential multiplied by the integrand - real(kind_phys), PARAMETER :: dcut = 1.2 ! max diameter of plume to parameterize relative to dx (km) - real(kind_phys):: d != -2.3 to -1.7 ;=-1.9 in Neggers paper; power law exponent for number density (N=Cl^d). + real(kind_phys),parameter :: & + & L0=100., & + & ENT0=0.1 + + ! Parameters/variables for regulating plumes: + real(kind_phys), parameter :: Atot = 0.10 ! Maximum total fractional area of all updrafts + real(kind_phys), parameter :: lmax = 1000.! diameter of largest plume (absolute maximum, can be smaller) + real(kind_phys), parameter :: lmin = 300. ! diameter of smallest plume (absolute minimum, can be larger) + real(kind_phys), parameter :: dlmin = 0. ! delta increase in the diameter of smallest plume (large fltv) + real(kind_phys) :: minwidth ! actual width of smallest plume + real(kind_phys) :: dl ! variable increment of plume size + real(kind_phys), parameter :: dcut = 1.2 ! max diameter of plume to parameterize relative to dx (km) + real(kind_phys):: d != -2.3 to -1.7 ;=-1.9 in Neggers paper; power law exponent for number density (N=Cl^d). ! Note that changing d to -2.0 makes each size plume equally contribute to the total coverage of all plumes. ! Note that changing d to -1.7 doubles the area coverage of the largest plumes relative to the smallest plumes. - real(kind_phys):: cn,c,l,n,an2,hux,maxwidth,wspd_pbl,cloud_base,width_flx + real(kind_phys):: cn,c,l,n,an2,hux,wspd_pbl,cloud_base,width_flx ! chem/smoke - INTEGER, INTENT(IN) :: nchem - real(kind_phys),DIMENSION(:, :) :: chem1 - real(kind_phys),DIMENSION(kts:kte+1, nchem) :: s_awchem - real(kind_phys),DIMENSION(nchem) :: chemn - real(kind_phys),DIMENSION(KTS:KTE+1,1:NUP, nchem) :: UPCHEM - INTEGER :: ic - real(kind_phys),DIMENSION(KTS:KTE+1, nchem) :: edmf_chem - LOGICAL, INTENT(IN) :: mix_chem + integer, intent(in) :: nchem + real(kind_phys),dimension(:, :) :: chem1 + real(kind_phys),dimension(kts:kte+1, nchem) :: s_awchem + real(kind_phys),dimension(nchem) :: chemn + real(kind_phys),dimension(kts:kte+1,1:NUP, nchem) :: UPCHEM + integer :: ic + real(kind_phys),dimension(kts:kte+1, nchem) :: edmf_chem + logical, intent(in) :: mix_chem !JOE: add declaration of ERF real(kind_phys):: ERF - LOGICAL :: superadiabatic + logical :: superadiabatic ! VARIABLES FOR CHABOUREAU-BECHTOLD CLOUD FRACTION - real(kind_phys),DIMENSION(KTS:KTE), INTENT(INOUT) :: vt, vq, sgm + real(kind_phys),dimension(kts:kte), intent(inout) :: vt, vq, sgm real(kind_phys):: sigq,xl,rsl,cpm,a,qmq,mf_cf,Aup,Q1,diffqt,qsat_tk,& Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid, & Ac_mf,Ac_strat,qc_mf - real(kind_phys), PARAMETER :: cf_thresh = 0.5 ! only overwrite stratus CF less than this value + real(kind_phys), parameter :: cf_thresh = 0.5 ! only overwrite stratus CF less than this value ! Variables for plume interpolation/saturation check - real(kind_phys),DIMENSION(KTS:KTE) :: exneri,dzi + real(kind_phys),dimension(kts:kte) :: exneri,dzi,rhoz real(kind_phys):: THp, QTp, QCp, QCs, esat, qsl - real(kind_phys):: csigma,acfac,ac_wsp,ac_cld + real(kind_phys):: csigma,acfac,ac_wsp !plume overshoot - INTEGER :: overshoot + integer :: overshoot real(kind_phys):: bvf, Frz, dzp !Flux limiter: not let mass-flux of heat between k=1&2 exceed (fluxportion)*(surface heat flux). !This limiter makes adjustments to the entire column. real(kind_phys):: adjustment, flx1 - real(kind_phys), PARAMETER :: fluxportion=0.75 ! set liberally, so has minimal impact. 0.5 starts to have a noticeable impact + real(kind_phys), parameter :: fluxportion=0.75 ! set liberally, so has minimal impact. Note that + ! 0.5 starts to have a noticeable impact ! over land (decrease maxMF by 10-20%), but no impact over water. !Subsidence - real(kind_phys),DIMENSION(KTS:KTE) :: sub_thl,sub_sqv,sub_u,sub_v, & !tendencies due to subsidence + real(kind_phys),dimension(kts:kte) :: sub_thl,sub_sqv,sub_u,sub_v, & !tendencies due to subsidence det_thl,det_sqv,det_sqc,det_u,det_v, & !tendencied due to detrainment envm_a,envm_w,envm_thl,envm_sqv,envm_sqc, & envm_u,envm_v !environmental variables defined at middle of layer - real(kind_phys),DIMENSION(KTS:KTE+1) :: envi_a,envi_w !environmental variables defined at model interface + real(kind_phys),dimension(kts:kte+1) :: envi_a,envi_w !environmental variables defined at model interface real(kind_phys):: temp,sublim,qc_ent,qv_ent,qt_ent,thl_ent,detrate, & detrateUV,oow,exc_fac,aratio,detturb,qc_grid,qc_sgs, & - qc_plume,exc_heat,exc_moist,tk_int - real(kind_phys), PARAMETER :: Cdet = 1./45. - real(kind_phys), PARAMETER :: dzpmax = 300. !limit dz used in detrainment - can be excessing in thick layers + qc_plume,exc_heat,exc_moist,tk_int,tvs + real(kind_phys), parameter :: Cdet = 1./45. + real(kind_phys), parameter :: dzpmax = 300. !limit dz used in detrainment - can be excessing in thick layers !parameter "Csub" determines the propotion of upward vertical velocity that contributes to !environmenatal subsidence. Some portion is expected to be compensated by downdrafts instead of !gentle environmental subsidence. 1.0 assumes all upward vertical velocity in the mass-flux scheme !is compensated by "gentle" environmental subsidence. - real(kind_phys), PARAMETER :: Csub=0.25 + real(kind_phys), parameter :: Csub=0.25 !Factor for the pressure gradient effects on momentum transport - real(kind_phys), PARAMETER :: pgfac = 0.00 ! Zhang and Wu showed 0.4 is more appropriate for lower troposphere + real(kind_phys), parameter :: pgfac = 0.00 ! Zhang and Wu showed 0.4 is more appropriate for lower troposphere real(kind_phys):: Uk,Ukm1,Vk,Vkm1,dxsa ! check the inputs @@ -5859,9 +5902,9 @@ SUBROUTINE DMP_mf( & UPQNWFA=0. UPQNIFA=0. UPQNBCA=0. - IF ( mix_chem ) THEN - UPCHEM(KTS:KTE+1,1:NUP,1:nchem)=0.0 - ENDIF + if ( mix_chem ) then + UPCHEM(kts:kte+1,1:NUP,1:nchem)=0.0 + endif ENT=0.001 ! Initialize mean updraft properties @@ -5871,9 +5914,9 @@ SUBROUTINE DMP_mf( & edmf_thl=0. edmf_ent=0. edmf_qc =0. - IF ( mix_chem ) THEN + if ( mix_chem ) then edmf_chem(kts:kte+1,1:nchem) = 0.0 - ENDIF + endif ! Initialize the variables needed for implicit solver s_aw=0. @@ -5889,153 +5932,163 @@ SUBROUTINE DMP_mf( & s_awqnwfa=0. s_awqnifa=0. s_awqnbca=0. - IF ( mix_chem ) THEN + if ( mix_chem ) then s_awchem(kts:kte+1,1:nchem) = 0.0 - ENDIF + endif ! Initialize explicit tendencies for subsidence & detrainment sub_thl = 0. sub_sqv = 0. - sub_u = 0. - sub_v = 0. + sub_u = 0. + sub_v = 0. det_thl = 0. det_sqv = 0. det_sqc = 0. - det_u = 0. - det_v = 0. + det_u = 0. + det_v = 0. + nup2 = nup !start with nup, but set to zero if activation criteria fails ! Taper off MF scheme when significant resolved-scale motions ! are present This function needs to be asymetric... - k = 1 - maxw = 0.0 + maxw = 0.0 cloud_base = 9000.0 -! DO WHILE (ZW(k) < pblh + 500.) - DO k=1,kte-1 - IF(zw(k) > pblh + 500.) exit + do k=1,kte-1 + if (zw(k) > pblh + 500.) exit wpbl = w(k) - IF(w(k) < 0.)wpbl = 2.*w(k) - maxw = MAX(maxw,ABS(wpbl)) + if (w(k) < 0.)wpbl = 2.*w(k) + maxw = max(maxw,abs(wpbl)) !Find highest k-level below 50m AGL - IF(ZW(k)<=50.)k50=k + if (ZW(k)<=50.)k50=k !Search for cloud base - qc_sgs = MAX(qc(k), qc_bl1d(k)*cldfra_bl1d(k)) - IF(qc_sgs> 1E-5 .AND. cloud_base == 9000.0)THEN + qc_sgs = max(qc(k), qc_bl1d(k)) + if (qc_sgs> 1E-5 .and. (cldfra_bl1d(k) .ge. 0.5) .and. cloud_base == 9000.0) then cloud_base = 0.5*(ZW(k)+ZW(k+1)) - ENDIF + endif + enddo - !k = k + 1 - ENDDO - !print*," maxw before manipulation=", maxw - maxw = MAX(0.,maxw - 1.0) ! do nothing for small w (< 1 m/s), but - Psig_w = MAX(0.0, 1.0 - maxw) ! linearly taper off for w > 1.0 m/s - Psig_w = MIN(Psig_w, Psig_shcu) - !print*," maxw=", maxw," Psig_w=",Psig_w," Psig_shcu=",Psig_shcu + !do nothing for small w (< 1 m/s), but linearly taper off for w > 1.0 m/s + maxw = max(0.,maxw - 1.0) + Psig_w = max(0.0, 1.0 - maxw) + Psig_w = min(Psig_w, Psig_shcu) !Completely shut off MF scheme for strong resolved-scale vertical velocities. fltv2 = fltv - IF(Psig_w == 0.0 .and. fltv > 0.0) fltv2 = -1.*fltv + if(Psig_w == 0.0 .and. fltv > 0.0) fltv2 = -1.*fltv ! If surface buoyancy is positive we do integration, otherwise no. ! Also, ensure that it is at least slightly superadiabatic up through 50 m superadiabatic = .false. - IF((landsea-1.5).GE.0)THEN + if ((landsea-1.5).ge.0) then hux = -0.001 ! WATER ! dT/dz must be < - 0.1 K per 100 m. - ELSE + else hux = -0.005 ! LAND ! dT/dz must be < - 0.5 K per 100 m. - ENDIF - DO k=1,MAX(1,k50-1) !use "-1" because k50 used interface heights (zw). - IF (k == 1) then - IF ((th(k)-ts)/(0.5*dz(k)) < hux) THEN + endif + tvs = ts*(1.0+p608*qv(kts)) + do k=1,max(1,k50-1) !use "-1" because k50 used interface heights (zw). + if (k == 1) then + if ((thv(k)-tvs)/(0.5*dz(k)) < hux) then superadiabatic = .true. - ELSE + else superadiabatic = .false. exit - ENDIF - ELSE - IF ((th(k)-th(k-1))/(0.5*(dz(k)+dz(k-1))) < hux) THEN + endif + else + if ((thv(k)-thv(k-1))/(0.5*(dz(k)+dz(k-1))) < hux) then superadiabatic = .true. - ELSE + else superadiabatic = .false. exit - ENDIF - ENDIF - ENDDO + endif + endif + enddo ! Determine the numer of updrafts/plumes in the grid column: ! Some of these criteria may be a little redundant but useful for bullet-proofing. - ! (1) largest plume = 1.0 * dx. - ! (2) Apply a scale-break, assuming no plumes with diameter larger than PBLH can exist. + ! (1) largest plume = 1.2 * dx. + ! (2) Apply a scale-break, assuming no plumes with diameter larger than 1.1*PBLH can exist. ! (3) max plume size beneath clouds deck approx = 0.5 * cloud_base. ! (4) add wspd-dependent limit, when plume model breaks down. (hurricanes) ! (5) limit to reduce max plume sizes in weakly forced conditions. This is only ! meant to "soften" the activation of the mass-flux scheme. ! Criteria (1) - NUP2 = max(1,min(NUP,INT(dx*dcut/dl))) + maxwidth = min(dx*dcut, lmax) !Criteria (2) - maxwidth = 1.1*PBLH + maxwidth = min(maxwidth, 1.1_kind_phys*PBLH) ! Criteria (3) - maxwidth = MIN(maxwidth,0.5*cloud_base) + if ((landsea-1.5) .lt. 0) then !land + maxwidth = MIN(maxwidth, 0.5_kind_phys*cloud_base) + else !water + maxwidth = MIN(maxwidth, 0.9_kind_phys*cloud_base) + endif ! Criteria (4) - wspd_pbl=SQRT(MAX(u(kts)**2 + v(kts)**2, 0.01)) + wspd_pbl=SQRT(MAX(u(kts)**2 + v(kts)**2, 0.01_kind_phys)) !Note: area fraction (acfac) is modified below ! Criteria (5) - only a function of flt (not fltv) if ((landsea-1.5).LT.0) then !land - !width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.050)/0.03) + .5),1000.), 0.) - width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.040)/0.03) + .5),1000.), 0.) + width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.040)/0.04) + .5),1000._kind_phys), 0._kind_phys) else !water - width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.003)/0.01) + .5),1000.), 0.) + width_flx = MAX(MIN(1000.*(0.6*tanh((fltv - 0.007)/0.02) + .5),1000._kind_phys), 0._kind_phys) + endif + maxwidth = MIN(maxwidth, width_flx) + minwidth = lmin + !allow min plume size to increase in large flux conditions (eddy diffusivity should be + !large enough to handle the representation of small plumes). + if (maxwidth .ge. (lmax - 1.0) .and. fltv .gt. 0.2)minwidth = lmin + dlmin*min((fltv-0.2)/0.3, 1._kind_phys) + + if (maxwidth .le. minwidth) then ! deactivate MF component + nup2 = 0 + maxwidth = 0.0 endif - maxwidth = MIN(maxwidth,width_flx) - ! Convert maxwidth to number of plumes - NUP2 = MIN(MAX(INT((maxwidth - MOD(maxwidth,100.))/100), 0), NUP2) - !Initialize values for 2d output fields: - ktop = 0 - ztop = 0.0 - maxmf= 0.0 + ! Initialize values for 2d output fields: + ktop = 0 + ztop = 0.0 + maxmf= 0.0 - IF ( fltv2 > 0.002 .AND. NUP2 .GE. 1 .AND. superadiabatic) then - !PRINT*," Conditions met to run mass-flux scheme",fltv2,pblh +!Begin plume processing if passes criteria +if ( fltv2 > 0.002 .AND. (maxwidth > minwidth) .AND. superadiabatic) then ! Find coef C for number size density N cn = 0. - d=-1.9 !set d to value suggested by Neggers 2015 (JAMES). - !d=-1.9 + .2*tanh((fltv2 - 0.05)/0.15) - do I=1,NUP !NUP2 - IF(I > NUP2) exit - l = dl*I ! diameter of plume + d =-1.9 !set d to value suggested by Neggers 2015 (JAMES). + dl = (maxwidth - minwidth)/real(nup-1,kind=kind_phys) + do i=1,NUP + ! diameter of plume + l = minwidth + dl*real(i-1) cn = cn + l**d * (l*l)/(dx*dx) * dl ! sum fractional area of each plume enddo C = Atot/cn !Normalize C according to the defined total fraction (Atot) ! Make updraft area (UPA) a function of the buoyancy flux if ((landsea-1.5).LT.0) then !land - !acfac = .5*tanh((fltv2 - 0.03)/0.09) + .5 - !acfac = .5*tanh((fltv2 - 0.02)/0.09) + .5 acfac = .5*tanh((fltv2 - 0.02)/0.05) + .5 else !water acfac = .5*tanh((fltv2 - 0.01)/0.03) + .5 endif !add a windspeed-dependent adjustment to acfac that tapers off - !the mass-flux scheme linearly above sfc wind speeds of 20 m/s: - ac_wsp = 1.0 - min(max(wspd_pbl - 20.0, 0.0), 10.0)/10.0 - !reduce area fraction beneath cloud bases < 1200 m AGL - ac_cld = min(cloud_base/1200., 1.0) - acfac = acfac * min(ac_wsp, ac_cld) + !the mass-flux scheme linearly above sfc wind speeds of 10 m/s. + !Note: this effect may be better represented by an increase in + !entrainment rate for high wind consitions (more ambient turbulence). + if (wspd_pbl .le. 10.) then + ac_wsp = 1.0 + else + ac_wsp = 1.0 - min((wspd_pbl - 10.0)/15., 1.0) + endif + acfac = acfac * ac_wsp ! Find the portion of the total fraction (Atot) of each plume size: An2 = 0. - do I=1,NUP !NUP2 - IF(I > NUP2) exit - l = dl*I ! diameter of plume + do i=1,NUP + ! diameter of plume + l = minwidth + dl*real(i-1) N = C*l**d ! number density of plume n - UPA(1,I) = N*l*l/(dx*dx) * dl ! fractional area of plume n + UPA(1,i) = N*l*l/(dx*dx) * dl ! fractional area of plume n - UPA(1,I) = UPA(1,I)*acfac - An2 = An2 + UPA(1,I) ! total fractional area of all plumes + UPA(1,i) = UPA(1,i)*acfac + An2 = An2 + UPA(1,i) ! total fractional area of all plumes !print*," plume size=",l,"; area=",UPA(1,I),"; total=",An2 end do @@ -6048,23 +6101,25 @@ SUBROUTINE DMP_mf( & qstar=max(flq,1.0E-5)/wstar thstar=flt/wstar - IF((landsea-1.5).GE.0)THEN + if ((landsea-1.5) .ge. 0) then csigma = 1.34 ! WATER - ELSE + else csigma = 1.34 ! LAND - ENDIF + endif if (env_subs) then exc_fac = 0.0 else if ((landsea-1.5).GE.0) then !water: increase factor to compensate for decreased pwmin/pwmax - exc_fac = 0.58*4.0*min(cloud_base/1000., 1.0) + exc_fac = 0.58*4.0 else !land: no need to increase factor - already sufficiently large superadiabatic layers exc_fac = 0.58 endif endif + !decrease excess for large wind speeds + exc_fac = exc_fac * ac_wsp !Note: sigmaW is typically about 0.5*wstar sigmaW =csigma*wstar*(z0/pblh)**(onethird)*(1 - 0.8*z0/pblh) @@ -6077,14 +6132,11 @@ SUBROUTINE DMP_mf( & wmax=MIN(sigmaW*pwmax,0.5) !SPECIFY SURFACE UPDRAFT PROPERTIES AT MODEL INTERFACE BETWEEN K = 1 & 2 - DO I=1,NUP !NUP2 - IF(I > NUP2) exit + do i=1,NUP wlv=wmin+(wmax-wmin)/NUP2*(i-1) !SURFACE UPDRAFT VERTICAL VELOCITY UPW(1,I)=wmin + real(i)/real(NUP)*(wmax-wmin) - !IF (UPW(1,I) > 0.5*ZW(2)/dt) UPW(1,I) = 0.5*ZW(2)/dt - UPU(1,I)=(U(KTS)*DZ(KTS+1)+U(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPV(1,I)=(V(KTS)*DZ(KTS+1)+V(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPQC(1,I)=0.0 @@ -6093,21 +6145,11 @@ SUBROUTINE DMP_mf( & exc_heat = exc_fac*UPW(1,I)*sigmaTH/sigmaW UPTHV(1,I)=(THV(KTS)*DZ(KTS+1)+THV(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) & & + exc_heat -!was UPTHL(1,I)= UPTHV(1,I)/(1.+svp1*UPQT(1,I)) !assume no saturated parcel at surface UPTHL(1,I)=(THL(KTS)*DZ(KTS+1)+THL(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) & & + exc_heat !calculate exc_moist by use of surface fluxes exc_moist=exc_fac*UPW(1,I)*sigmaQT/sigmaW - !calculate exc_moist by conserving rh: -! tk_int =(tk(kts)*dz(kts+1)+tk(kts+1)*dz(kts))/(dz(kts+1)+dz(kts)) -! pk =(p(kts)*dz(kts+1)+p(kts+1)*dz(kts))/(dz(kts+1)+dz(kts)) -! qtk =(qt(kts)*dz(kts+1)+qt(kts+1)*dz(kts))/(dz(kts)+dz(kts+1)) -! qsat_tk = qsat_blend(tk_int, pk) ! saturation water vapor mixing ratio at tk and p -! rhgrid =MAX(MIN(1.0,qtk/MAX(1.E-8,qsat_tk)),0.001) -! tk_int = tk_int + exc_heat -! qsat_tk = qsat_blend(tk_int, pk) -! exc_moist= max(rhgrid*qsat_tk - qtk, 0.0) UPQT(1,I)=(QT(KTS)*DZ(KTS+1)+QT(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1))& & +exc_moist @@ -6117,36 +6159,36 @@ SUBROUTINE DMP_mf( & UPQNWFA(1,I)=(QNWFA(KTS)*DZ(KTS+1)+QNWFA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPQNIFA(1,I)=(QNIFA(KTS)*DZ(KTS+1)+QNIFA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) UPQNBCA(1,I)=(QNBCA(KTS)*DZ(KTS+1)+QNBCA(KTS+1)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) - ENDDO + enddo - IF ( mix_chem ) THEN - DO I=1,NUP !NUP2 - IF(I > NUP2) exit + if ( mix_chem ) then + do i=1,NUP do ic = 1,nchem - UPCHEM(1,I,ic)=(chem1(KTS,ic)*DZ(KTS+1)+chem1(KTS+1,ic)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) + UPCHEM(1,i,ic)=(chem1(KTS,ic)*DZ(KTS+1)+chem1(KTS+1,ic)*DZ(KTS))/(DZ(KTS)+DZ(KTS+1)) enddo - ENDDO - ENDIF + enddo + endif !Initialize environmental variables which can be modified by detrainment - DO k=kts,kte - envm_thl(k)=THL(k) - envm_sqv(k)=QV(k) - envm_sqc(k)=QC(k) - envm_u(k)=U(k) - envm_v(k)=V(k) - ENDDO + envm_thl(kts:kte)=THL(kts:kte) + envm_sqv(kts:kte)=QV(kts:kte) + envm_sqc(kts:kte)=QC(kts:kte) + envm_u(kts:kte)=U(kts:kte) + envm_v(kts:kte)=V(kts:kte) + do k=kts,kte-1 + rhoz(k) = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k)) + enddo + rhoz(kte) = rho(kte) !dxsa is scale-adaptive factor governing the pressure-gradient term of the momentum transport dxsa = 1. - MIN(MAX((12000.0-dx)/(12000.0-3000.0), 0.), 1.) ! do integration updraft - DO I=1,NUP !NUP2 - IF(I > NUP2) exit + do i=1,NUP QCn = 0. overshoot = 0 - l = dl*I ! diameter of plume - DO k=KTS+1,KTE-1 + l = minwidth + dl*real(i-1) ! diameter of plume + do k=kts+1,kte-1 !Entrainment from Tian and Kuang (2016) !ENT(k,i) = 0.35/(MIN(MAX(UPW(K-1,I),0.75),1.9)*l) wmin = 0.3 + l*0.0005 !* MAX(pblh-ZW(k+1), 0.0)/pblh @@ -6161,7 +6203,7 @@ SUBROUTINE DMP_mf( & ENT(k,i) = max(ENT(k,i),0.0003) !ENT(k,i) = max(ENT(k,i),0.05/ZW(k)) !not needed for Tian and Kuang - !JOE - increase entrainment for plumes extending very high. + !increase entrainment for plumes extending very high. IF(ZW(k) >= MIN(pblh+1500., 4000.))THEN ENT(k,i)=ENT(k,i) + (ZW(k)-MIN(pblh+1500.,4000.))*5.0E-6 ENDIF @@ -6339,6 +6381,7 @@ SUBROUTINE DMP_mf( & exit !exit k-loop END IF ENDDO + IF (debug_mf == 1) THEN IF (MAXVAL(UPW(:,I)) > 10.0 .OR. MINVAL(UPA(:,I)) < 0.0 .OR. & MAXVAL(UPA(:,I)) > Atot .OR. NUP2 > 10) THEN @@ -6358,30 +6401,26 @@ SUBROUTINE DMP_mf( & ENDIF ENDIF ENDDO - ELSE +ELSE !At least one of the conditions was not met for activating the MF scheme. NUP2=0. - END IF !end criteria for mass-flux scheme +END IF !end criteria check for mass-flux scheme - ktop=MIN(ktop,KTE-1) ! Just to be safe... - IF (ktop == 0) THEN - ztop = 0.0 - ELSE - ztop=zw(ktop) - ENDIF - - IF(nup2 > 0) THEN +ktop=MIN(ktop,KTE-1) +IF (ktop == 0) THEN + ztop = 0.0 +ELSE + ztop=zw(ktop) +ENDIF - !Calculate the fluxes for each variable - !All s_aw* variable are == 0 at k=1 - DO i=1,NUP !NUP2 - IF(I > NUP2) exit +IF (nup2 > 0) THEN + !Calculate the fluxes for each variable + !All s_aw* variable are == 0 at k=1 + DO i=1,NUP DO k=KTS,KTE-1 - IF(k > ktop) exit - rho_int = (rho(k)*DZ(k+1)+rho(k+1)*DZ(k))/(DZ(k+1)+DZ(k)) - s_aw(k+1) = s_aw(k+1) + rho_int*UPA(K,i)*UPW(K,i)*Psig_w - s_awthl(k+1)= s_awthl(k+1) + rho_int*UPA(K,i)*UPW(K,i)*UPTHL(K,i)*Psig_w - s_awqt(k+1) = s_awqt(k+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQT(K,i)*Psig_w + s_aw(k+1) = s_aw(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*Psig_w + s_awthl(k+1)= s_awthl(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPTHL(K,i)*Psig_w + s_awqt(k+1) = s_awqt(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQT(K,i)*Psig_w !to conform to grid mean properties, move qc to qv in grid mean !saturated layers, so total water fluxes are preserved but !negative qc fluxes in unsaturated layers is reduced. @@ -6390,72 +6429,76 @@ SUBROUTINE DMP_mf( & ! else ! qc_plume = 0.0 ! endif - s_awqc(k+1) = s_awqc(k+1) + rho_int*UPA(K,i)*UPW(K,i)*qc_plume*Psig_w - IF (momentum_opt > 0) THEN - s_awu(k+1) = s_awu(k+1) + rho_int*UPA(K,i)*UPW(K,i)*UPU(K,i)*Psig_w - s_awv(k+1) = s_awv(k+1) + rho_int*UPA(K,i)*UPW(K,i)*UPV(K,i)*Psig_w - ENDIF - IF (tke_opt > 0) THEN - s_awqke(k+1)= s_awqke(k+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQKE(K,i)*Psig_w - ENDIF + s_awqc(k+1) = s_awqc(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*qc_plume*Psig_w s_awqv(k+1) = s_awqt(k+1) - s_awqc(k+1) ENDDO - ENDDO - - IF ( mix_chem ) THEN - DO k=KTS,KTE - IF(k > KTOP) exit - rho_int = (rho(k)*DZ(k+1)+rho(k+1)*DZ(k))/(DZ(k+1)+DZ(k)) - DO i=1,NUP !NUP2 - IF(I > NUP2) exit - do ic = 1,nchem - s_awchem(k+1,ic) = s_awchem(k+1,ic) + rho_int*UPA(K,i)*UPW(K,i)*UPCHEM(K,i,ic)*Psig_w - enddo - ENDDO - ENDDO - ENDIF - - IF (scalar_opt > 0) THEN - DO k=KTS,KTE - IF(k > KTOP) exit - rho_int = (rho(k)*DZ(k+1)+rho(k+1)*DZ(k))/(DZ(k+1)+DZ(k)) - DO I=1,NUP !NUP2 - IF (I > NUP2) exit - s_awqnc(k+1)= s_awqnc(K+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQNC(K,i)*Psig_w - s_awqni(k+1)= s_awqni(K+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQNI(K,i)*Psig_w - s_awqnwfa(k+1)= s_awqnwfa(K+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQNWFA(K,i)*Psig_w - s_awqnifa(k+1)= s_awqnifa(K+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQNIFA(K,i)*Psig_w - s_awqnbca(k+1)= s_awqnbca(K+1) + rho_int*UPA(K,i)*UPW(K,i)*UPQNBCA(K,i)*Psig_w - ENDDO - ENDDO - ENDIF + ENDDO + !momentum + if (momentum_opt > 0) then + do i=1,nup + do k=kts,kte-1 + s_awu(k+1) = s_awu(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPU(K,i)*Psig_w + s_awv(k+1) = s_awv(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPV(K,i)*Psig_w + enddo + enddo + endif + !tke + if (tke_opt > 0) then + do i=1,nup + do k=kts,kte-1 + s_awqke(k+1)= s_awqke(k+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQKE(K,i)*Psig_w + enddo + enddo + endif + !chem + if ( mix_chem ) then + do k=kts,kte + do i=1,nup + do ic = 1,nchem + s_awchem(k+1,ic) = s_awchem(k+1,ic) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPCHEM(K,i,ic)*Psig_w + enddo + enddo + enddo + endif + + if (scalar_opt > 0) then + do k=kts,kte + do I=1,nup + s_awqnc(k+1) = s_awqnc(K+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQNC(K,i)*Psig_w + s_awqni(k+1) = s_awqni(K+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQNI(K,i)*Psig_w + s_awqnwfa(k+1)= s_awqnwfa(K+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQNWFA(K,i)*Psig_w + s_awqnifa(k+1)= s_awqnifa(K+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQNIFA(K,i)*Psig_w + s_awqnbca(k+1)= s_awqnbca(K+1) + rhoz(k)*UPA(K,i)*UPW(K,i)*UPQNBCA(K,i)*Psig_w + enddo + enddo + endif - !Flux limiter: Check ratio of heat flux at top of first model layer - !and at the surface. Make sure estimated flux out of the top of the - !layer is < fluxportion*surface_heat_flux - IF (s_aw(kts+1) /= 0.) THEN + !Flux limiter: Check ratio of heat flux at top of first model layer + !and at the surface. Make sure estimated flux out of the top of the + !layer is < fluxportion*surface_heat_flux + IF (s_aw(kts+1) /= 0.) THEN dzi(kts) = 0.5*(DZ(kts)+DZ(kts+1)) !dz centered at model interface flx1 = MAX(s_aw(kts+1)*(TH(kts)-TH(kts+1))/dzi(kts),1.0e-5) - ELSE + ELSE flx1 = 0.0 !print*,"ERROR: s_aw(kts+1) == 0, NUP=",NUP," NUP2=",NUP2,& ! " superadiabatic=",superadiabatic," KTOP=",KTOP - ENDIF - adjustment=1.0 - !Print*,"Flux limiter in MYNN-EDMF, adjustment=",fluxportion*flt/dz(kts)/flx1 - !Print*,"flt/dz=",flt/dz(kts)," flx1=",flx1," s_aw(kts+1)=",s_aw(kts+1) - IF (flx1 > fluxportion*flt/dz(kts) .AND. flx1>0.0) THEN + ENDIF + adjustment=1.0 + !Print*,"Flux limiter in MYNN-EDMF, adjustment=",fluxportion*flt/dz(kts)/flx1 + !Print*,"flt/dz=",flt/dz(kts)," flx1=",flx1," s_aw(kts+1)=",s_aw(kts+1) + IF (flx1 > fluxportion*flt/dz(kts) .AND. flx1>0.0) THEN adjustment= fluxportion*flt/dz(kts)/flx1 - s_aw = s_aw*adjustment - s_awthl= s_awthl*adjustment - s_awqt = s_awqt*adjustment - s_awqc = s_awqc*adjustment - s_awqv = s_awqv*adjustment - s_awqnc= s_awqnc*adjustment - s_awqni= s_awqni*adjustment - s_awqnwfa= s_awqnwfa*adjustment - s_awqnifa= s_awqnifa*adjustment - s_awqnbca= s_awqnbca*adjustment + s_aw = s_aw*adjustment + s_awthl = s_awthl*adjustment + s_awqt = s_awqt*adjustment + s_awqc = s_awqc*adjustment + s_awqv = s_awqv*adjustment + s_awqnc = s_awqnc*adjustment + s_awqni = s_awqni*adjustment + s_awqnwfa = s_awqnwfa*adjustment + s_awqnifa = s_awqnifa*adjustment + s_awqnbca = s_awqnbca*adjustment IF (momentum_opt > 0) THEN s_awu = s_awu*adjustment s_awv = s_awv*adjustment @@ -6467,62 +6510,57 @@ SUBROUTINE DMP_mf( & s_awchem = s_awchem*adjustment ENDIF UPA = UPA*adjustment - ENDIF - !Print*,"adjustment=",adjustment," fluxportion=",fluxportion," flt=",flt - - !Calculate mean updraft properties for output: - !all edmf_* variables at k=1 correspond to the interface at top of first model layer - DO k=KTS,KTE-1 - IF(k > KTOP) exit - rho_int = (rho(k)*DZ(k+1)+rho(k+1)*DZ(k))/(DZ(k+1)+DZ(k)) - DO I=1,NUP !NUP2 - IF(I > NUP2) exit - edmf_a(K) =edmf_a(K) +UPA(K,i) - edmf_w(K) =edmf_w(K) +rho_int*UPA(K,i)*UPW(K,i) - edmf_qt(K) =edmf_qt(K) +rho_int*UPA(K,i)*UPQT(K,i) - edmf_thl(K)=edmf_thl(K)+rho_int*UPA(K,i)*UPTHL(K,i) - edmf_ent(K)=edmf_ent(K)+rho_int*UPA(K,i)*ENT(K,i) - edmf_qc(K) =edmf_qc(K) +rho_int*UPA(K,i)*UPQC(K,i) - ENDDO - + ENDIF + !Print*,"adjustment=",adjustment," fluxportion=",fluxportion," flt=",flt + + !Calculate mean updraft properties for output: + !all edmf_* variables at k=1 correspond to the interface at top of first model layer + do k=kts,kte-1 + do I=1,nup + edmf_a(K) =edmf_a(K) +UPA(K,i) + edmf_w(K) =edmf_w(K) +rhoz(k)*UPA(K,i)*UPW(K,i) + edmf_qt(K) =edmf_qt(K) +rhoz(k)*UPA(K,i)*UPQT(K,i) + edmf_thl(K)=edmf_thl(K)+rhoz(k)*UPA(K,i)*UPTHL(K,i) + edmf_ent(K)=edmf_ent(K)+rhoz(k)*UPA(K,i)*ENT(K,i) + edmf_qc(K) =edmf_qc(K) +rhoz(k)*UPA(K,i)*UPQC(K,i) + enddo + enddo + do k=kts,kte-1 !Note that only edmf_a is multiplied by Psig_w. This takes care of the !scale-awareness of the subsidence below: - IF (edmf_a(k)>0.) THEN - edmf_w(k)=edmf_w(k)/edmf_a(k) - edmf_qt(k)=edmf_qt(k)/edmf_a(k) - edmf_thl(k)=edmf_thl(k)/edmf_a(k) - edmf_ent(k)=edmf_ent(k)/edmf_a(k) - edmf_qc(k)=edmf_qc(k)/edmf_a(k) - edmf_a(k)=edmf_a(k)*Psig_w - - !FIND MAXIMUM MASS-FLUX IN THE COLUMN: - IF(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k) - ENDIF - ENDDO ! end k - - !smoke/chem - IF ( mix_chem ) THEN - DO k=kts,kte-1 - IF(k > KTOP) exit - rho_int = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k)) - DO I=1,NUP !NUP2 - IF(I > NUP2) exit + if (edmf_a(k)>0.) then + edmf_w(k)=edmf_w(k)/edmf_a(k) + edmf_qt(k)=edmf_qt(k)/edmf_a(k) + edmf_thl(k)=edmf_thl(k)/edmf_a(k) + edmf_ent(k)=edmf_ent(k)/edmf_a(k) + edmf_qc(k)=edmf_qc(k)/edmf_a(k) + edmf_a(k)=edmf_a(k)*Psig_w + !FIND MAXIMUM MASS-FLUX IN THE COLUMN: + if(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k) + endif + enddo ! end k + + !smoke/chem + if ( mix_chem ) then + do k=kts,kte-1 + do I=1,nup do ic = 1,nchem - edmf_chem(k,ic) = edmf_chem(k,ic) + rho_int*UPA(K,I)*UPCHEM(k,i,ic) + edmf_chem(k,ic) = edmf_chem(k,ic) + rhoz(k)*UPA(K,I)*UPCHEM(k,i,ic) enddo - ENDDO - - IF (edmf_a(k)>0.) THEN + enddo + enddo + do k=kts,kte-1 + if (edmf_a(k)>0.) then do ic = 1,nchem edmf_chem(k,ic) = edmf_chem(k,ic)/edmf_a(k) enddo - ENDIF - ENDDO ! end k - ENDIF + endif + enddo ! end k + endif - !Calculate the effects environmental subsidence. - !All envi_*variables are valid at the interfaces, like the edmf_* variables - IF (env_subs) THEN + !Calculate the effects environmental subsidence. + !All envi_*variables are valid at the interfaces, like the edmf_* variables + IF (env_subs) THEN DO k=kts+1,kte-1 !First, smooth the profiles of w & a, since sharp vertical gradients !in plume variables are not likely extended to env variables @@ -6557,18 +6595,16 @@ SUBROUTINE DMP_mf( & !calculate tendencies from subsidence and detrainment valid at the middle of !each model layer. The lowest model layer uses an assumes w=0 at the surface. dzi(kts) = 0.5*(dz(kts)+dz(kts+1)) - rho_int = (rho(kts)*dz(kts+1)+rho(kts+1)*dz(kts))/(dz(kts+1)+dz(kts)) sub_thl(kts)= 0.5*envi_w(kts)*envi_a(kts)* & - (rho(kts+1)*thl(kts+1)-rho(kts)*thl(kts))/dzi(kts)/rho_int + (rho(kts+1)*thl(kts+1)-rho(kts)*thl(kts))/dzi(kts)/rhoz(k) sub_sqv(kts)= 0.5*envi_w(kts)*envi_a(kts)* & - (rho(kts+1)*qv(kts+1)-rho(kts)*qv(kts))/dzi(kts)/rho_int + (rho(kts+1)*qv(kts+1)-rho(kts)*qv(kts))/dzi(kts)/rhoz(k) DO k=kts+1,kte-1 dzi(k) = 0.5*(dz(k)+dz(k+1)) - rho_int = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k)) sub_thl(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * & - (rho(k+1)*thl(k+1)-rho(k)*thl(k))/dzi(k)/rho_int + (rho(k+1)*thl(k+1)-rho(k)*thl(k))/dzi(k)/rhoz(k) sub_sqv(k)= 0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * & - (rho(k+1)*qv(k+1)-rho(k)*qv(k))/dzi(k)/rho_int + (rho(k+1)*qv(k+1)-rho(k)*qv(k))/dzi(k)/rhoz(k) ENDDO DO k=KTS,KTE-1 @@ -6578,17 +6614,15 @@ SUBROUTINE DMP_mf( & ENDDO IF (momentum_opt > 0) THEN - rho_int = (rho(kts)*dz(kts+1)+rho(kts+1)*dz(kts))/(dz(kts+1)+dz(kts)) sub_u(kts)=0.5*envi_w(kts)*envi_a(kts)* & - (rho(kts+1)*u(kts+1)-rho(kts)*u(kts))/dzi(kts)/rho_int + (rho(kts+1)*u(kts+1)-rho(kts)*u(kts))/dzi(kts)/rhoz(k) sub_v(kts)=0.5*envi_w(kts)*envi_a(kts)* & - (rho(kts+1)*v(kts+1)-rho(kts)*v(kts))/dzi(kts)/rho_int + (rho(kts+1)*v(kts+1)-rho(kts)*v(kts))/dzi(kts)/rhoz(k) DO k=kts+1,kte-1 - rho_int = (rho(k)*dz(k+1)+rho(k+1)*dz(k))/(dz(k+1)+dz(k)) sub_u(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * & - (rho(k+1)*u(k+1)-rho(k)*u(k))/dzi(k)/rho_int + (rho(k+1)*u(k+1)-rho(k)*u(k))/dzi(k)/rhoz(k) sub_v(k)=0.5*(envi_w(k)+envi_w(k-1))*0.5*(envi_a(k)+envi_a(k-1)) * & - (rho(k+1)*v(k+1)-rho(k)*v(k))/dzi(k)/rho_int + (rho(k+1)*v(k+1)-rho(k)*v(k))/dzi(k)/rhoz(k) ENDDO DO k=KTS,KTE-1 @@ -6596,23 +6630,23 @@ SUBROUTINE DMP_mf( & det_v(k) = Cdet*(envm_v(k)-v(k))*envi_a(k)*Psig_w ENDDO ENDIF - ENDIF !end subsidence/env detranment + ENDIF !end subsidence/env detranment - !First, compute exner, plume theta, and dz centered at interface - !Here, k=1 is the top of the first model layer. These values do not - !need to be defined at k=kte (unused level). - DO K=KTS,KTE-1 - exneri(k) = (exner(k)*DZ(k+1)+exner(k+1)*DZ(k))/(DZ(k+1)+DZ(k)) + !First, compute exner, plume theta, and dz centered at interface + !Here, k=1 is the top of the first model layer. These values do not + !need to be defined at k=kte (unused level). + DO K=KTS,KTE-1 + exneri(k) = (exner(k)*dz(k+1)+exner(k+1)*dz(k))/(dz(k+1)+dz(k)) edmf_th(k)= edmf_thl(k) + xlvcp/exneri(k)*edmf_qc(K) - dzi(k) = 0.5*(DZ(k)+DZ(k+1)) - ENDDO + dzi(k) = 0.5*(dz(k)+dz(k+1)) + ENDDO !JOE: ADD CLDFRA_bl1d, qc_bl1d. Note that they have already been defined in ! mym_condensation. Here, a shallow-cu component is added, but no cumulus ! clouds can be added at k=1 (start loop at k=2). - do k=kts+1,kte-2 - IF(k > KTOP) exit - IF(0.5*(edmf_qc(k)+edmf_qc(k-1))>0.0 .and. (cldfra_bl1d(k) < cf_thresh))THEN + do k=kts+1,kte-2 + if (k > KTOP) exit + if(0.5*(edmf_qc(k)+edmf_qc(k-1))>0.0 .and. (cldfra_bl1d(k) < cf_thresh))THEN !interpolate plume quantities to mass levels Aup = (edmf_a(k)*dzi(k-1)+edmf_a(k-1)*dzi(k))/(dzi(k-1)+dzi(k)) THp = (edmf_th(k)*dzi(k-1)+edmf_th(k-1)*dzi(k))/(dzi(k-1)+dzi(k)) @@ -6686,8 +6720,8 @@ SUBROUTINE DMP_mf( & !mf_cf = min(max(0.5 + 0.36 * atan(1.20*(Q1+0.4)),0.01),0.6) !Original CB mf_cf = min(max(0.5 + 0.36 * atan(1.55*Q1),0.01),0.6) - mf_cf = max(mf_cf, 1.75 * Aup) - mf_cf = min(mf_cf, 5.0 * Aup) + mf_cf = max(mf_cf, 1.8 * Aup) + mf_cf = min(mf_cf, 5.0 * Aup) endif !IF ( debug_code ) THEN @@ -6705,10 +6739,7 @@ SUBROUTINE DMP_mf( & if (QCp * Aup > 5e-5) then qc_bl1d(k) = 1.86 * (QCp * Aup) - 2.2e-5 else - qc_bl1d(k) = 1.18 * (QCp * Aup) - endif - if (mf_cf .ge. Aup) then - qc_bl1d(k) = qc_bl1d(k) / mf_cf + qc_bl1d(k) = 1.18 * (QCp * Aup) endif cldfra_bl1d(k) = mf_cf Ac_mf = mf_cf @@ -6718,9 +6749,6 @@ SUBROUTINE DMP_mf( & else qc_bl1d(k) = 1.18 * (QCp * Aup) endif - if (mf_cf .ge. Aup) then - qc_bl1d(k) = qc_bl1d(k) / mf_cf - endif cldfra_bl1d(k) = mf_cf Ac_mf = mf_cf endif @@ -6752,13 +6780,13 @@ SUBROUTINE DMP_mf( & endif !check for (qc in plume) .and. (cldfra_bl < threshold) enddo !k-loop - ENDIF !end nup2 > 0 +ENDIF !end nup2 > 0 - !modify output (negative: dry plume, positive: moist plume) - if (ktop > 0) then - maxqc = maxval(edmf_qc(1:ktop)) - if ( maxqc < 1.E-8) maxmf = -1.0*maxmf - endif +!modify output (negative: dry plume, positive: moist plume) +if (ktop > 0) then + maxqc = maxval(edmf_qc(1:ktop)) + if ( maxqc < 1.E-8) maxmf = -1.0*maxmf +endif ! ! debugging @@ -6927,62 +6955,68 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & &qc_bl1d,cldfra_bl1d, & &rthraten ) - INTEGER, INTENT(IN) :: KTS,KTE,KPBL - real(kind_phys),DIMENSION(KTS:KTE), INTENT(IN) :: U,V,TH,THL,TK,QT,QV,QC,& - THV,P,rho,exner,dz - real(kind_phys),DIMENSION(KTS:KTE), INTENT(IN) :: rthraten + integer, intent(in) :: KTS,KTE,KPBL + real(kind_phys), dimension(kts:kte), intent(in) :: & + U,V,TH,THL,TK,QT,QV,QC,THV,P,rho,exner,dz + real(kind_phys), dimension(kts:kte), intent(in) :: rthraten ! zw .. heights of the downdraft levels (edges of boxes) - real(kind_phys),DIMENSION(KTS:KTE+1), INTENT(IN) :: ZW - real(kind_phys), INTENT(IN) :: WTHL,WQT - real(kind_phys), INTENT(IN) :: dt,ust,pblh + real(kind_phys), dimension(kts:kte+1), intent(in) :: ZW + real(kind_phys), intent(in) :: WTHL,WQT + real(kind_phys), intent(in) :: dt,ust,pblh ! outputs - downdraft properties - real(kind_phys),DIMENSION(KTS:KTE), INTENT(OUT) :: edmf_a_dd,edmf_w_dd, & - & edmf_qt_dd,edmf_thl_dd, edmf_ent_dd,edmf_qc_dd + real(kind_phys), dimension(kts:kte), intent(out) :: & + edmf_a_dd,edmf_w_dd, & + edmf_qt_dd,edmf_thl_dd, edmf_ent_dd,edmf_qc_dd ! outputs - variables needed for solver (sd_aw - sum ai*wi, sd_awphi - sum ai*wi*phii) - real(kind_phys),DIMENSION(KTS:KTE+1) :: sd_aw, sd_awthl, sd_awqt, sd_awu, & - sd_awv, sd_awqc, sd_awqv, sd_awqke, sd_aw2 + real(kind_phys), dimension(kts:kte+1) :: & + sd_aw, sd_awthl, sd_awqt, sd_awu, & + sd_awv, sd_awqc, sd_awqv, sd_awqke, sd_aw2 - real(kind_phys),DIMENSION(KTS:KTE), INTENT(IN) :: qc_bl1d, cldfra_bl1d + real(kind_phys), dimension(kts:kte), intent(in) :: & + qc_bl1d, cldfra_bl1d - INTEGER, PARAMETER :: NDOWN=5, debug_mf=0 !fixing number of plumes to 5 + integer, parameter:: ndown = 5 ! draw downdraft starting height randomly between cloud base and cloud top - INTEGER, DIMENSION(1:NDOWN) :: DD_initK - real(kind_phys) , DIMENSION(1:NDOWN) :: randNum + integer, dimension(1:NDOWN) :: DD_initK + real(kind_phys), dimension(1:NDOWN) :: randNum ! downdraft properties - real(kind_phys),DIMENSION(KTS:KTE+1,1:NDOWN) :: DOWNW,DOWNTHL,DOWNQT,& - DOWNQC,DOWNA,DOWNU,DOWNV,DOWNTHV + real(kind_phys), dimension(kts:kte+1,1:NDOWN) :: & + DOWNW,DOWNTHL,DOWNQT,DOWNQC,DOWNA,DOWNU,DOWNV,DOWNTHV ! entrainment variables - Real(Kind_phys),DIMENSION(KTS+1:KTE+1,1:NDOWN) :: ENT,ENTf - INTEGER,DIMENSION(KTS+1:KTE+1,1:NDOWN) :: ENTi + real(kind_phys), dimension(KTS+1:KTE+1,1:NDOWN) :: ENT,ENTf + integer, dimension(KTS+1:KTE+1,1:NDOWN) :: ENTi ! internal variables - INTEGER :: K,I,ki, kminrad, qlTop, p700_ind, qlBase - real(kind_phys):: wthv,wstar,qstar,thstar,sigmaW,sigmaQT,sigmaTH,z0, & - pwmin,pwmax,wmin,wmax,wlv,wtv,went,mindownw - real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,Wn,THVk,Pk, & - EntEXP,EntW, Beta_dm, EntExp_M, rho_int - real(kind_phys):: jump_thetav, jump_qt, jump_thetal, & + integer :: K,I,ki, kminrad, qlTop, p700_ind, qlBase + real(kind_phys):: wthv,wstar,qstar,thstar,sigmaW,sigmaQT, & + sigmaTH,z0,pwmin,pwmax,wmin,wmax,wlv,wtv,went,mindownw + real(kind_phys):: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,Wn, & + THVk,Pk,EntEXP,EntW,beta_dm,EntExp_M,rho_int + real(kind_phys):: jump_thetav, jump_qt, jump_thetal, & refTHL, refTHV, refQT ! DD specific internal variables real(kind_phys):: minrad,zminrad, radflux, F0, wst_rad, wst_dd logical :: cloudflg - - real(kind_phys):: sigq,xl,rsl,cpm,a,mf_cf,diffqt,& + real(kind_phys):: sigq,xl,rsl,cpm,a,mf_cf,diffqt, & Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid ! w parameters - real(kind_phys),PARAMETER :: & - &Wa=1., & - &Wb=1.5,& - &Z00=100.,& - &BCOEFF=0.2 + real(kind_phys),parameter :: & + &Wa=1., Wb=1.5, Z00=100., BCOEFF=0.2 ! entrainment parameters - real(kind_phys),PARAMETER :: & - & L0=80,& - & ENT0=0.2 - + real(kind_phys),parameter :: & + &L0=80, ENT0=0.2 + !downdraft properties + real(kind_phys):: & + & dp, & !diameter of plume + & dl, & !diameter increment + & Adn !total area of downdrafts + !additional printouts for debugging + integer, parameter :: debug_mf=0 + + dl = (1000.-500.)/real(ndown) pwmin=-3. ! drawing from the negative tail -3sigma to -1sigma pwmax=-1. @@ -7052,6 +7086,14 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & if ( radflux < 0.0 ) F0 = abs(radflux) + F0 enddo F0 = max(F0, 1.0) + + !Allow the total fractional area of the downdrafts to be proportional + !to the radiative forcing: + !for 50 W/m2, Adn = 0.10 + !for 100 W/m2, Adn = 0.15 + !for 150 W/m2, Adn = 0.20 + Adn = min( 0.05 + F0*0.001, 0.3) + !found Sc cloud and cloud not at surface, trigger downdraft if (cloudflg) then @@ -7066,14 +7108,14 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & ! call Poisson(1,NDOWN,kts+1,kte,ENTf,ENTi) - ! entrainent: Ent=Ent0/dz*P(dz/L0) - do i=1,NDOWN - do k=kts+1,kte -! ENT(k,i)=real(ENTi(k,i))*Ent0/(ZW(k+1)-ZW(k)) - ENT(k,i) = 0.002 - ENT(k,i) = min(ENT(k,i),0.9/(ZW(k+1)-ZW(k))) - enddo - enddo +! ! entrainent: Ent=Ent0/dz*P(dz/L0) +! do i=1,NDOWN +! do k=kts+1,kte +!! ENT(k,i)=real(ENTi(k,i))*Ent0/(ZW(k+1)-ZW(k)) +! ENT(k,i) = 0.002 +! ENT(k,i) = min(ENT(k,i),0.9/(ZW(k+1)-ZW(k))) +! enddo +! enddo !!![EW: INVJUMP] find 700mb height then subtract trpospheric lapse rate!!! p700_ind = MINLOC(ABS(p-70000),1)!p1D is 70000 @@ -7116,8 +7158,10 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & !DOWNW(ki,I)=0.5*(wlv+wtv) DOWNW(ki,I)=wlv + !multiply downa by cloud fraction, so it's impact will diminish if + !clouds are mixed away over the course of the longer radiation time step !DOWNA(ki,I)=0.5*ERF(wtv/(sqrt(2.)*sigmaW))-0.5*ERF(wlv/(sqrt(2.)*sigmaW)) - DOWNA(ki,I)=.1/real(NDOWN) + DOWNA(ki,I)=Adn/real(NDOWN) DOWNU(ki,I)=(u(ki-1)*DZ(ki) + u(ki)*DZ(ki-1)) /(DZ(ki)+DZ(ki-1)) DOWNV(ki,I)=(v(ki-1)*DZ(ki) + v(ki)*DZ(ki-1)) /(DZ(ki)+DZ(ki-1)) @@ -7144,16 +7188,21 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & enddo - !print*, " Begin integration of downdrafts:" DO I=1,NDOWN + dp = 500. + dl*real(I) ! diameter of plume (meters) !print *, "Plume # =", I,"=======================" DO k=DD_initK(I)-1,KTS+1,-1 + + !Entrainment from Tian and Kuang (2016), with constraints + wmin = 0.3 + dp*0.0005 + ENT(k,i) = 0.33/(MIN(MAX(-1.0*DOWNW(k+1,I),wmin),0.9)*dp) + !starting at the first interface level below cloud top !EntExp=exp(-ENT(K,I)*dz(k)) !EntExp_M=exp(-ENT(K,I)/3.*dz(k)) - EntExp =ENT(K,I)*dz(k) - EntExp_M=ENT(K,I)*0.333*dz(k) + EntExp =ENT(K,I)*dz(k) !for all scalars + EntExp_M=ENT(K,I)*0.333*dz(k) !test for momentum QTn =DOWNQT(k+1,I) *(1.-EntExp) + QT(k)*EntExp THLn=DOWNTHL(k+1,I)*(1.-EntExp) + THL(k)*EntExp @@ -7187,11 +7236,11 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & BCOEFF*B/mindownw)*MIN(dz(k), 250.) !Do not allow a parcel to accelerate more than 1.25 m/s over 200 m. - !Add max increase of 2.0 m/s for coarse vertical resolution. - IF (Wn < DOWNW(K+1,I) - MIN(1.25*dz(k)/200., 2.0))THEN - Wn = DOWNW(K+1,I) - MIN(1.25*dz(k)/200., 2.0) + !Add max acceleration of -2.0 m/s for coarse vertical resolution. + IF (Wn < DOWNW(K+1,I) - MIN(1.25*dz(k)/200., -2.0))THEN + Wn = DOWNW(K+1,I) - MIN(1.25*dz(k)/200., -2.0) ENDIF - !Add symmetrical max decrease in w + !Add symmetrical max decrease in velocity (less negative) IF (Wn > DOWNW(K+1,I) + MIN(1.25*dz(k)/200., 2.0))THEN Wn = DOWNW(K+1,I) + MIN(1.25*dz(k)/200., 2.0) ENDIF @@ -7237,7 +7286,6 @@ SUBROUTINE DDMF_JPL(kts,kte,dt,zw,dz,p, & ! Even though downdraft starts at different height, average all up to qlTop DO k=qlTop,KTS,-1 DO I=1,NDOWN - IF (I > NDOWN) exit edmf_a_dd(K) =edmf_a_dd(K) +DOWNA(K-1,I) edmf_w_dd(K) =edmf_w_dd(K) +DOWNA(K-1,I)*DOWNW(K-1,I) edmf_qt_dd(K) =edmf_qt_dd(K) +DOWNA(K-1,I)*DOWNQT(K-1,I) @@ -7287,8 +7335,8 @@ SUBROUTINE SCALE_AWARE(dx,PBL1,Psig_bl,Psig_shcu) ! Psig_bl tapers local mixing ! Psig_shcu tapers nonlocal mixing - real(kind_phys), INTENT(IN) :: dx,pbl1 - real(kind_phys), INTENT(OUT) :: Psig_bl,Psig_shcu + real(kind_phys), intent(in) :: dx,pbl1 + real(kind_phys), intent(out) :: Psig_bl,Psig_shcu real(kind_phys) :: dxdh Psig_bl=1.0 @@ -7361,28 +7409,28 @@ FUNCTION esat_blend(t) IMPLICIT NONE - real(kind_phys), INTENT(IN):: t + real(kind_phys), intent(in):: t real(kind_phys):: esat_blend,XC,ESL,ESI,chi !liquid - real(kind_phys), PARAMETER:: J0= .611583699E03 - real(kind_phys), PARAMETER:: J1= .444606896E02 - real(kind_phys), PARAMETER:: J2= .143177157E01 - real(kind_phys), PARAMETER:: J3= .264224321E-1 - real(kind_phys), PARAMETER:: J4= .299291081E-3 - real(kind_phys), PARAMETER:: J5= .203154182E-5 - real(kind_phys), PARAMETER:: J6= .702620698E-8 - real(kind_phys), PARAMETER:: J7= .379534310E-11 - real(kind_phys), PARAMETER:: J8=-.321582393E-13 + real(kind_phys), parameter:: J0= .611583699E03 + real(kind_phys), parameter:: J1= .444606896E02 + real(kind_phys), parameter:: J2= .143177157E01 + real(kind_phys), parameter:: J3= .264224321E-1 + real(kind_phys), parameter:: J4= .299291081E-3 + real(kind_phys), parameter:: J5= .203154182E-5 + real(kind_phys), parameter:: J6= .702620698E-8 + real(kind_phys), parameter:: J7= .379534310E-11 + real(kind_phys), parameter:: J8=-.321582393E-13 !ice - real(kind_phys), PARAMETER:: K0= .609868993E03 - real(kind_phys), PARAMETER:: K1= .499320233E02 - real(kind_phys), PARAMETER:: K2= .184672631E01 - real(kind_phys), PARAMETER:: K3= .402737184E-1 - real(kind_phys), PARAMETER:: K4= .565392987E-3 - real(kind_phys), PARAMETER:: K5= .521693933E-5 - real(kind_phys), PARAMETER:: K6= .307839583E-7 - real(kind_phys), PARAMETER:: K7= .105785160E-9 - real(kind_phys), PARAMETER:: K8= .161444444E-12 + real(kind_phys), parameter:: K0= .609868993E03 + real(kind_phys), parameter:: K1= .499320233E02 + real(kind_phys), parameter:: K2= .184672631E01 + real(kind_phys), parameter:: K3= .402737184E-1 + real(kind_phys), parameter:: K4= .565392987E-3 + real(kind_phys), parameter:: K5= .521693933E-5 + real(kind_phys), parameter:: K6= .307839583E-7 + real(kind_phys), parameter:: K7= .105785160E-9 + real(kind_phys), parameter:: K8= .161444444E-12 XC=MAX(-80.,t - t0c) !note t0c = 273.15, tice is set in module mynn_common to 240 @@ -7412,28 +7460,28 @@ FUNCTION qsat_blend(t, P) IMPLICIT NONE - real(kind_phys), INTENT(IN):: t, P + real(kind_phys), intent(in):: t, P real(kind_phys):: qsat_blend,XC,ESL,ESI,RSLF,RSIF,chi !liquid - real(kind_phys), PARAMETER:: J0= .611583699E03 - real(kind_phys), PARAMETER:: J1= .444606896E02 - real(kind_phys), PARAMETER:: J2= .143177157E01 - real(kind_phys), PARAMETER:: J3= .264224321E-1 - real(kind_phys), PARAMETER:: J4= .299291081E-3 - real(kind_phys), PARAMETER:: J5= .203154182E-5 - real(kind_phys), PARAMETER:: J6= .702620698E-8 - real(kind_phys), PARAMETER:: J7= .379534310E-11 - real(kind_phys), PARAMETER:: J8=-.321582393E-13 + real(kind_phys), parameter:: J0= .611583699E03 + real(kind_phys), parameter:: J1= .444606896E02 + real(kind_phys), parameter:: J2= .143177157E01 + real(kind_phys), parameter:: J3= .264224321E-1 + real(kind_phys), parameter:: J4= .299291081E-3 + real(kind_phys), parameter:: J5= .203154182E-5 + real(kind_phys), parameter:: J6= .702620698E-8 + real(kind_phys), parameter:: J7= .379534310E-11 + real(kind_phys), parameter:: J8=-.321582393E-13 !ice - real(kind_phys), PARAMETER:: K0= .609868993E03 - real(kind_phys), PARAMETER:: K1= .499320233E02 - real(kind_phys), PARAMETER:: K2= .184672631E01 - real(kind_phys), PARAMETER:: K3= .402737184E-1 - real(kind_phys), PARAMETER:: K4= .565392987E-3 - real(kind_phys), PARAMETER:: K5= .521693933E-5 - real(kind_phys), PARAMETER:: K6= .307839583E-7 - real(kind_phys), PARAMETER:: K7= .105785160E-9 - real(kind_phys), PARAMETER:: K8= .161444444E-12 + real(kind_phys), parameter:: K0= .609868993E03 + real(kind_phys), parameter:: K1= .499320233E02 + real(kind_phys), parameter:: K2= .184672631E01 + real(kind_phys), parameter:: K3= .402737184E-1 + real(kind_phys), parameter:: K4= .565392987E-3 + real(kind_phys), parameter:: K5= .521693933E-5 + real(kind_phys), parameter:: K6= .307839583E-7 + real(kind_phys), parameter:: K7= .105785160E-9 + real(kind_phys), parameter:: K8= .161444444E-12 XC=MAX(-80.,t - t0c) @@ -7470,7 +7518,7 @@ FUNCTION xl_blend(t) IMPLICIT NONE - real(kind_phys), INTENT(IN):: t + real(kind_phys), intent(in):: t real(kind_phys):: xl_blend,xlvt,xlst,chi !note: t0c = 273.15, tice is set in mynn_common @@ -7499,11 +7547,11 @@ FUNCTION phim(zet) ! stable conditions [z/L ~ O(10)]. IMPLICIT NONE - real(kind_phys), INTENT(IN):: zet + real(kind_phys), intent(in):: zet real(kind_phys):: dummy_0,dummy_1,dummy_11,dummy_2,dummy_22,dummy_3,dummy_33,dummy_4,dummy_44,dummy_psi - real(kind_phys), PARAMETER :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st - real(kind_phys), PARAMETER :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st - real(kind_phys), PARAMETER :: am_unst=10., ah_unst=34. + real(kind_phys), parameter :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st + real(kind_phys), parameter :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st + real(kind_phys), parameter :: am_unst=10., ah_unst=34. real(kind_phys):: phi_m,phim if ( zet >= 0.0 ) then @@ -7551,11 +7599,11 @@ FUNCTION phih(zet) ! stable conditions [z/L ~ O(10)]. IMPLICIT NONE - real(kind_phys), INTENT(IN):: zet + real(kind_phys), intent(in):: zet real(kind_phys):: dummy_0,dummy_1,dummy_11,dummy_2,dummy_22,dummy_3,dummy_33,dummy_4,dummy_44,dummy_psi - real(kind_phys), PARAMETER :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st - real(kind_phys), PARAMETER :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st - real(kind_phys), PARAMETER :: am_unst=10., ah_unst=34. + real(kind_phys), parameter :: am_st=6.1, bm_st=2.5, rbm_st=1./bm_st + real(kind_phys), parameter :: ah_st=5.3, bh_st=1.1, rbh_st=1./bh_st + real(kind_phys), parameter :: am_unst=10., ah_unst=34. real(kind_phys):: phh,phih if ( zet >= 0.0 ) then diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 271db11d0..44e552160 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -993,6 +993,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rainprod, evapprod, & #endif refl_10cm, diagflag, do_radar_ref, & + max_hail_diam_sfc, & vt_dbz_wt, first_time_step, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & @@ -1062,6 +1063,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & GRAUPELNC, GRAUPELNCV REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & refl_10cm + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & + max_hail_diam_sfc REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & vt_dbz_wt LOGICAL, INTENT(IN) :: first_time_step @@ -1416,6 +1419,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & qcten1(k) = 0. endif initialize_extended_diagnostics enddo + lsml = lsm(i,j) if (is_aerosol_aware .or. merra2_aerosol_aware) then do k = kts, kte nc1d(k) = nc(i,k,j) @@ -1423,7 +1427,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nifa1d(k) = nifa(i,k,j) enddo else - lsml = lsm(i,j) do k = kts, kte if(lsml == 1) then nc1d(k) = Nt_c_l/rho(k) @@ -1679,6 +1682,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & (nsteps>1 .and. istep==nsteps) .or. & (nsteps==1 .and. ndt==1)) THEN + max_hail_diam_sfc(i,j) = hail_mass_99th_percentile(kts, kte, qg1d, t1d, p1d, qv1d) + !> - Call calc_refl10cm() diagflag_present: IF ( PRESENT (diagflag) ) THEN @@ -2464,17 +2469,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo - + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -3541,17 +3536,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo - + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -3589,7 +3574,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION if (clap .gt. eps) then if (is_aerosol_aware .or. merra2_aerosol_aware) then - xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k))) + xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k), lsml)) else if(lsml == 1) then xnc = Nt_c_l @@ -5366,14 +5351,15 @@ end subroutine table_ccnAct ! TO_DO ITEM: For radiation cooling producing fog, in which case the !.. updraft velocity could easily be negative, we could use the temp !.. and its tendency to diagnose a pretend postive updraft velocity. - real function activ_ncloud(Tt, Ww, NCCN) + real function activ_ncloud(Tt, Ww, NCCN, lsm_in) implicit none REAL, INTENT(IN):: Tt, Ww, NCCN + INTEGER, INTENT(IN):: lsm_in REAL:: n_local, w_local INTEGER:: i, j, k, l, m, n REAL:: A, B, C, D, t, u, x1, x2, y1, y2, nx, wy, fraction - + REAL:: lower_lim_nuc_frac ! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc ! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw @@ -5420,6 +5406,14 @@ real function activ_ncloud(Tt, Ww, NCCN) l = 3 m = 2 + if (lsm_in .eq. 1) then ! land + lower_lim_nuc_frac = 0. + else if (lsm_in .eq. 0) then ! water + lower_lim_nuc_frac = 0.15 + else + lower_lim_nuc_frac = 0.15 ! catch-all for anything else + endif + A = tnccn_act(i-1,j-1,k,l,m) B = tnccn_act(i,j-1,k,l,m) C = tnccn_act(i,j,k,l,m) @@ -5434,7 +5428,8 @@ real function activ_ncloud(Tt, Ww, NCCN) ! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1)) fraction = (1.0-t)*(1.0-u)*A + t*(1.0-u)*B + t*u*C + (1.0-t)*u*D - + fraction = MAX(fraction, lower_lim_nuc_frac) + ! if (NCCN*fraction .gt. 0.75*Nt_c_max) then ! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k ! endif @@ -6085,16 +6080,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then - do k = kte, kts, -1 - ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 - N0_exp = 10.**(zans1) - N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 - lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg - ilamg(k) = 1./lamg - N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) - enddo + call graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) endif !+---+-----------------------------------------------------------------+ @@ -6471,6 +6457,88 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) END SUBROUTINE semi_lagrange_sedim +!>\ingroup aathompson +!! @brief Calculates graupel size distribution parameters +!! +!! Calculates graupel intercept and slope parameters for +!! for a vertical column +!! +!! @param[in] kts integer start index for vertical column +!! @param[in] kte integer end index for vertical column +!! @param[in] rand1 real random number for stochastic physics +!! @param[in] rg real array, size(kts:kte) for graupel mass concentration [kg m^3] +!! @param[out] ilamg double array, size(kts:kte) for inverse graupel slope parameter [m] +!! @param[out] N0_g double array, size(kts:kte) for graupel intercept paramter [m-4] +subroutine graupel_psd_parameters(kts, kte, rand1, rg, ilamg, N0_g) + + implicit none + + integer, intent(in) :: kts, kte + real, intent(in) :: rand1 + real, intent(in) :: rg(:) + double precision, intent(out) :: ilamg(:), N0_g(:) + + integer :: k + real :: ygra1, zans1 + double precision :: N0_exp, lam_exp, lamg + + do k = kte, kts, -1 + ygra1 = alog10(max(1.e-9, rg(k))) + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 + N0_exp = 10.**(zans1) + N0_exp = max(dble(gonv_min), min(N0_exp, dble(gonv_max))) + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ilamg(k) = 1./lamg + N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + enddo + +end subroutine graupel_psd_parameters + +!>\ingroup aathompson +!! @brief Calculates graupel/hail maximum diameter +!! +!! Calculates graupel/hail maximum diameter (currently the 99th percentile of mass distribtuion) +!! for a vertical column +!! +!! @param[in] kts integer start index for vertical column +!! @param[in] kte integer end index for vertical column +!! @param[in] qg real array, size(kts:kte) for graupel mass mixing ratio [kg kg^-1] +!! @param[in] temperature double array, size(kts:kte) temperature [K] +!! @param[in] pressure double array, size(kts:kte) pressure [Pa] +!! @param[in] qv real array, size(kts:kte) water vapor mixing ratio [kg kg^-1] +!! @param[out] max_hail_diam real maximum hail diameter [m] +function hail_mass_99th_percentile(kts, kte, qg, temperature, pressure, qv) result(max_hail_diam) + + implicit none + + integer, intent(in) :: kts, kte + real, intent(in) :: qg(:), temperature(:), pressure(:), qv(:) + real :: max_hail_diam + + integer :: k + real :: rho(kts:kte), rg(kts:kte), max_hail_column(kts:kte) + double precision :: ilamg(kts:kte), N0_g(kts:kte) + real, parameter :: random_number = 0. + + max_hail_column = 0. + rg = 0. + do k = kts, kte + rho(k) = 0.622*pressure(k)/(R*temperature(k)*(max(1.e-10, qv(k))+0.622)) + if (qg(k) .gt. R1) then + rg(k) = qg(k)*rho(k) + else + rg(k) = R1 + endif + enddo + + call graupel_psd_parameters(kts, kte, random_number, rg, ilamg, N0_g) + + where(rg .gt. 1.e-9) max_hail_column = 10.05 * ilamg + max_hail_diam = max_hail_column(kts) + +end function hail_mass_99th_percentile + !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index eecc5493c..3d847348d 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -26,7 +26,7 @@ MODULE module_sf_mynn ! ! LAND only: ! "iz0tlnd" namelist option is used to select the following momentum options: -! (default) =0: Zilitinkevich (1995); Czil now set to 0.085 +! (default) =0: Zilitinkevich (1995); Czil now set to 0.095 ! =1: Czil_new (modified according to Chen & Zhang 2008) ! =2: Modified Yang et al (2002, 2008) - generalized for all landuse ! =3: constant zt = z0/7.4 (original form; Garratt 1992) @@ -225,7 +225,7 @@ SUBROUTINE SFCLAY_mynn( & ! (water =1: z0 from Davis et al (2008), zt & zq from COARE3.0/3.5 ! only) =2: z0 from Davis et al (2008), zt & zq from Garratt (1992) ! =3: z0 from Taylor and Yelland (2004), zt and zq from COARE 3.0/3.5 -!-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.085, +!-- iz0tlnd =0: Zilitinkevich (1995) with Czil=0.095, ! (land =1: Czil_new (modified according to Chen & Zhang 2008) ! only) =2: Modified Yang et al (2002, 2008) - generalized for all landuse ! =3: constant zt = z0/7.4 (Garratt 1992) @@ -947,7 +947,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE: THSK_lnd(I) = TSK_lnd(I)*THCON(I) !(K) THVSK_lnd(I) = THSK_lnd(I)*(1.+EP1*qsfc_lnd(I)) - if(THVSK_lnd(I) < 170. .or. THVSK_lnd(I) > 360.) & + if(THVSK_lnd(I) < 160. .or. THVSK_lnd(I) > 390.) & print *,'THVSK_lnd(I)',itimestep,i,THVSK_lnd(I),THSK_lnd(i),tsurf_lnd(i),tskin_lnd(i),qsfc_lnd(i) endif if(icy(i)) then @@ -956,7 +956,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE: THSK_ice(I) = TSK_ice(I)*THCON(I) !(K) THVSK_ice(I) = THSK_ice(I)*(1.+EP1*qsfc_ice(I)) !(K) - if(THVSK_ice(I) < 170. .or. THVSK_ice(I) > 360.) & + if(THVSK_ice(I) < 160. .or. THVSK_ice(I) > 390.) & print *,'THVSK_ice(I)',itimestep,i,THVSK_ice(I),THSK_ice(i),tsurf_ice(i),tskin_ice(i),qsfc_ice(i) endif if(wet(i)) then @@ -965,7 +965,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & ! CONVERT SKIN TEMPERATURES TO POTENTIAL TEMPERATURE: THSK_wat(I) = TSK_wat(I)*THCON(I) !(K) THVSK_wat(I) = THSK_wat(I)*(1.+EP1*QVSH(I)) !(K) - if(THVSK_wat(I) < 170. .or. THVSK_wat(I) > 360.) & + if(THVSK_wat(I) < 160. .or. THVSK_wat(I) > 390.) & print *,'THVSK_wat(I)',i,THVSK_wat(I),THSK_wat(i),tsurf_wat(i),tskin_wat(i),qsfc_wat(i) endif endif ! flag_iter @@ -1380,6 +1380,8 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, & else ZNTstoch_lnd(I) = ZNT_lnd(I) endif + !add limit to prevent ridiculous values of z0 (more than dz/15) + ZNTstoch_lnd(I) = min(ZNTstoch_lnd(I), dz8w1d(i)*0.0666_kind_phys) !-------------------------------------- ! LAND @@ -2604,7 +2606,7 @@ SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,& IF ( IZ0TLND2 .EQ. 1 ) THEN CZIL = 10.0 ** ( -0.40 * ( Z_0 / 0.07 ) ) ELSE - CZIL = 0.085 !0.075 !0.10 + CZIL = 0.095 !0.075 !0.10 END IF Zt = Z_0*EXP(-KARMAN*CZIL*SQRT(restar)) diff --git a/physics/module_sf_noahmp_glacier.F90 b/physics/module_sf_noahmp_glacier.F90 index 6e34c43af..fcbe40a70 100644 --- a/physics/module_sf_noahmp_glacier.F90 +++ b/physics/module_sf_noahmp_glacier.F90 @@ -2652,7 +2652,7 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in !to obtain equilibrium state of snow in glacier region - if(sneqv > mwd) then ! 100 mm -> maximum water depth + if(sneqv > mwd .and. isnow /= 0) then ! 100 mm -> maximum water depth bdsnow = snice(0) / dzsnso(0) snoflow = (sneqv - mwd) snice(0) = snice(0) - snoflow diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index 160127e43..52fbc8123 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -97,6 +97,7 @@ SUBROUTINE LSMRUC(xlat,xlon, & MAVAIL,CANWAT,VEGFRA, & ALB,ZNT,Z0,SNOALB,ALBBCK,LAI, & landusef, nlcat, soilctop, nscat, & + smcwlt, smcref, & QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, & TBOT,IVGTYP,ISLTYP,XLAND, & ISWATER,ISICE,XICE,XICE_THRESHOLD, & @@ -107,6 +108,7 @@ SUBROUTINE LSMRUC(xlat,xlon, & RUNOFF1,RUNOFF2,ACRUNOFF,SFCEXC, & SFCEVP,GRDFLX,SNOWFALLAC,ACSNOW,SNOM, & SMFR3D,KEEPFR3DFLAG, & + add_fire_heat_flux,fire_heat_flux, & myj,shdmin,shdmax,rdlai2d, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & @@ -239,6 +241,8 @@ SUBROUTINE LSMRUC(xlat,xlon, & real (kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMIN real (kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: hgt real (kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: stdev + LOGICAL, intent(in) :: add_fire_heat_flux + real (kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: fire_heat_flux LOGICAL, intent(in) :: rdlai2d real (kind_phys), DIMENSION( 1:nsl), INTENT(IN ) :: ZS @@ -252,6 +256,8 @@ SUBROUTINE LSMRUC(xlat,xlon, & SNOALB, & ALB, & LAI, & + SMCWLT, & + SMCREF, & EMISS, & EMISBCK, & MAVAIL, & @@ -757,6 +763,8 @@ SUBROUTINE LSMRUC(xlat,xlon, & !-- update background emissivity for land points, can have vegetation mosaic effect EMISBCK(I,J) = EMISSL(I,J) + smcwlt(i,j) = wilt + smcref(i,j) = ref IF (debug_print ) THEN if(init)then @@ -961,6 +969,7 @@ SUBROUTINE LSMRUC(xlat,xlon, & snoalb(i,j),albbck(i,j),lai(i,j), & hgt(i,j),stdev(i,j), & !new myj,seaice(i,j),isice, & + add_fire_heat_flux,fire_heat_flux(i,j), & !--- soil fixed fields QWRTZ, & rhocs,dqm,qmin,ref, & @@ -1212,6 +1221,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia QKMS,TKMS,PC,MAVAIL,CST,VEGFRA,ALB,ZNT, & ALB_SNOW,ALB_SNOW_FREE,lai,hgt,stdev, & MYJ,SEAICE,ISICE, & + add_fire_heat_flux,fire_heat_flux, & QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & !--- soil fixed fields sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, & !--- constants @@ -1256,7 +1266,9 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia SEAICE, & RHO, & QKMS, & - TKMS + TKMS, & + fire_heat_flux + LOGICAL, INTENT(IN ) :: add_fire_heat_flux INTEGER, INTENT(IN ) :: IVGTYP, ISLTYP !--- 2-D variables @@ -1509,7 +1521,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia ENDIF if(snhei.gt.0.0081_kind_phys*rhowater/rhosn) then -!*** Update snow density for current temperature (Koren et al. 1999) +!*** Update snow density for current temperature (Koren et al 1999,doi:10.1029/1999JD900232.) BSN=delt/3600._kind_phys*c1sn*exp(0.08_kind_phys*min(zero,tsnav)-c2sn*rhosn*1.e-3_kind_phys) if(bsn*snwe*100._kind_phys.lt.1.e-4_kind_phys) goto 777 XSN=rhosn*(exp(bsn*snwe*100._kind_phys)-one)/(bsn*snwe*100._kind_phys) @@ -1813,6 +1825,13 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia UPFLUX = T3 *SOILT XINET = EMISS_snowfree*(GLW-UPFLUX) RNET = GSWnew + XINET + IF ( add_fire_heat_flux .and. fire_heat_flux >0 ) then ! JLS + IF (debug_print ) THEN + print *,'RNET snow-free, fire_heat_flux, xlat/xlon',RNET, fire_heat_flux,xlat,xlon + ENDIF + RNET = RNET + fire_heat_flux + ENDIF + IF (debug_print ) THEN print *,'Fractional snow - snowfrac=',snowfrac print *,'Snowfrac<1 GSWin,GSWnew -',GSWin,GSWnew,'SOILT, RNET',soilt,rnet @@ -1837,7 +1856,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia ilands = ivgtyp - CALL SOIL(debug_print,xlat,xlon, & + CALL SOIL(debug_print,xlat, xlon, testptlat, testptlon,& !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,gswin, & @@ -1933,6 +1952,12 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia if (SEAICE .LT. 0.5_kind_phys) then ! LAND + IF ( add_fire_heat_flux .and. fire_heat_flux>0 ) then ! JLS + IF (debug_print ) THEN + print *,'RNET snow, fire_heat_flux, xlat/xlon',RNET, fire_heat_flux,xlat,xlon + ENDIF + RNET = RNET + fire_heat_flux + ENDIF if(snow_mosaic==one)then snfr=one else @@ -2223,7 +2248,14 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia if(SEAICE .LT. 0.5_kind_phys) then ! LAND - CALL SOIL(debug_print,xlat,xlon, & + IF ( add_fire_heat_flux .and. fire_heat_flux>0) then ! JLS + IF (debug_print ) THEN + print *,'RNET no snow, fire_heat_flux, xlat/xlon',RNET, fire_heat_flux,xlat,xlon + endif + RNET = RNET + fire_heat_flux + ENDIF + + CALL SOIL(debug_print,xlat, xlon, testptlat, testptlon,& !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,GSWin, & @@ -2316,7 +2348,7 @@ END FUNCTION QSN !>\ingroup lsm_ruc_group !> This subroutine calculates energy and moisture budget for vegetated surfaces !! without snow, heat diffusion and Richards eqns in soil. - SUBROUTINE SOIL (debug_print,xlat,xlon, & + SUBROUTINE SOIL (debug_print,xlat,xlon,testptlat,testptlon,& i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,& !--- input variables PRCPMS,RAINF,PATM,QVATM,QCATM, & GLW,GSW,GSWin,EMISS,RNET, & @@ -2398,7 +2430,8 @@ SUBROUTINE SOIL (debug_print,xlat,xlon, & INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,iland,isoil - real (kind_phys), INTENT(IN ) :: DELT,CONFLX,xlat,xlon + real (kind_phys), INTENT(IN ) :: DELT,CONFLX + real (kind_phys), INTENT(IN ) :: xlat,xlon,testptlat,testptlon LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables real (kind_phys), & @@ -2622,6 +2655,7 @@ SUBROUTINE SOIL (debug_print,xlat,xlon, & ! hydraulic condeuctivities !****************************************************************** CALL SOILPROP( debug_print, & + xlat, xlon, testptlat, testptlon, & !--- input variables nzs,fwsat,lwsat,tav,keepfr, & soilmois,soiliqw,soilice, & @@ -2657,6 +2691,7 @@ SUBROUTINE SOIL (debug_print,xlat,xlon, & ! TRANSF computes transpiration function !************************************************************** CALL TRANSF(debug_print, & + xlat, xlon, testptlat, testptlon, & !--- input variables nzs,nroot,soiliqw,tabs,lai,gswin, & !--- soil fixed fields @@ -2714,7 +2749,7 @@ SUBROUTINE SOIL (debug_print,xlat,xlon, & ! SOILTEMP soilves heat budget and diffusion eqn. in soil !************************************************************** - CALL SOILTEMP(debug_print,xlat,xlon, & + CALL SOILTEMP(debug_print,xlat,xlon,testptlat,testptlon,& !--- input variables i,j,iland,isoil, & delt,ktau,conflx,nzs,nddzs,nroot, & @@ -2784,6 +2819,7 @@ SUBROUTINE SOIL (debug_print,xlat,xlon, & ! and Richards eqn. !************************************************************************* CALL SOILMOIST (debug_print, & + xlat, xlon, testptlat, testptlon, & !-- input delt,nzs,nddzs,DTDZS,DTDZS2,RIW, & zsmain,zshalf,diffu,hydro, & @@ -3578,6 +3614,7 @@ SUBROUTINE SNOWSOIL ( debug_print,xlat,xlon, & ! hydraulic condeuctivities !****************************************************************** CALL SOILPROP(debug_print, & + xlat, xlon, testptlat, testptlon, & !--- input variables nzs,fwsat,lwsat,tav,keepfr, & soilmois,soiliqw,soilice, & @@ -3628,6 +3665,7 @@ SUBROUTINE SNOWSOIL ( debug_print,xlat,xlon, & ! TRANSF computes transpiration function !************************************************************** CALL TRANSF(debug_print, & + xlat, xlon, testptlat, testptlon, & !--- input variables nzs,nroot,soiliqw,tabs,lai,gswin, & !--- soil fixed fields @@ -3723,7 +3761,7 @@ SUBROUTINE SNOWSOIL ( debug_print,xlat,xlon, & !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28) ! AND TSO,ETA PROFILES !************************************************************************* - CALL SOILMOIST (debug_print, & + CALL SOILMOIST (debug_print,xlat,xlon,testptlat,testptlon,& !-- input delt,nzs,nddzs,DTDZS,DTDZS2,RIW, & zsmain,zshalf,diffu,hydro, & @@ -4046,35 +4084,25 @@ SUBROUTINE SNOWSEAICE( debug_print,xlat,xlon, & RHOnewCSN=sheatsn * RHOnewSN if(isncond_opt == 1) then - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then - !-- some areas with large snow depth have unrealistically - !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. - !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys - else - !-- old version thdifsn = 0.265/RHOCSN - THDIFSN = 0.265_kind_phys/RHOCSN - endif + !-- old version thdifsn = 0.265/RHOCSN + THDIFSN = 0.265_kind_phys/RHOCSN else !-- 07Jun19 - thermal conductivity (K_eff) from Sturm et al.(1997) !-- keff = 10. ** (2.650 * RHOSN*1.e-3 - 1.652) fact = one if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys)) then keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys - !-- fact is added by tgs based on 4 Jan 2017 testing - fact = 5._kind_phys else keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys - fact = 2._kind_phys endif - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then + if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys) then !-- some areas with large snow depth have unrealistically !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. + !-- Based on Sturm et al. keff=0.452 typical for hard snow slabs + !-- with rhosn=488 kg/m^3. Thdifsn = 0.452/(2090*488)=4.431718e-7 !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys + thdifsn = 4.431718e-7_kind_phys else thdifsn = keff/rhocsn * fact endif @@ -4510,35 +4538,25 @@ SUBROUTINE SNOWSEAICE( debug_print,xlat,xlon, & RHOCSN=sheatsn* RHOSN if(isncond_opt == 1) then - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then - !-- some areas with large snow depth have unrealistically - !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. - !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys - else - !-- old version thdifsn = 0.265/RHOCSN - THDIFSN = 0.265_kind_phys/RHOCSN - endif + !-- old version thdifsn = 0.265/RHOCSN + THDIFSN = 0.265_kind_phys/RHOCSN else !-- 07Jun19 - thermal conductivity (K_eff) from Sturm et al.(1997) !-- keff = 10. ** (2.650 * RHOSN*1.e-3 - 1.652) fact = one if(rhosn < 156._kind_phys .or. (newsn > zero .and. rhonewsn < 156._kind_phys)) then keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys - !-- fact is added by tgs based on 4 Jan 2017 testing - fact = 5._kind_phys else keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys - fact = 2._kind_phys endif - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then + if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys) then !-- some areas with large snow depth have unrealistically !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. + !-- Based on Sturm et al. keff=0.452 typical for hard snow slabs + !-- with rhosn=488 kg/m^3. Thdifsn = 0.452/(2090*488)=4.431718e-7 !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys + thdifsn = 4.431718e-7_kind_phys else thdifsn = keff/rhocsn * fact endif @@ -4679,7 +4697,7 @@ END SUBROUTINE SNOWSEAICE !>\ingroup lsm_ruc_group !> This subroutine solves energy budget equation and heat diffusion !! equation. - SUBROUTINE SOILTEMP( debug_print,xlat,xlon, & + SUBROUTINE SOILTEMP( debug_print,xlat,xlon,testptlat,testptlon,& i,j,iland,isoil, & !--- input variables delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, & @@ -4749,7 +4767,8 @@ SUBROUTINE SOILTEMP( debug_print,xlat,xlon, & INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,iland,isoil - real (kind_phys), INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF,xlat,xlon + real (kind_phys), INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF + real (kind_phys), INTENT(IN ) :: xlat, xlon, testptlat, testptlon real (kind_phys), INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM !--- 3-D Atmospheric variables real (kind_phys), & @@ -5193,27 +5212,16 @@ SUBROUTINE SNOWTEMP( debug_print,xlat,xlon, & RHOCSN=sheatsn* RHOSN RHOnewCSN=sheatsn* RHOnewSN if(isncond_opt == 1) then - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then - !-- some areas with large snow depth have unrealistically - !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. - !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys - else - !-- old version thdifsn = 0.265/RHOCSN - THDIFSN = 0.265_kind_phys/RHOCSN - endif + !-- old version thdifsn = 0.265/RHOCSN + THDIFSN = 0.265_kind_phys/RHOCSN else !-- 07Jun19 - thermal conductivity (K_eff) from Sturm et al.(1997) !-- keff = 10. ** (2.650 * RHOSN*1.e-3 - 1.652) fact = one if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys)) then keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys - !-- fact is added by tgs based on 4 Jan 2017 testing - fact = 5._kind_phys else keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys - fact = 2._kind_phys if(debug_print) then print *,'SnowTemp xlat,xlon,rhosn,keff', xlat,xlon,rhosn,keff,keff/rhocsn*fact print *,'SNOWTEMP - 0.265/rhocsn',0.265_kind_phys/rhocsn @@ -5223,12 +5231,13 @@ SUBROUTINE SNOWTEMP( debug_print,xlat,xlon, & print *,'SNOWTEMP - xlat,xlon,newsnow,rhonewsn,rhosn,fact,keff',xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff endif - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then + if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys) then !-- some areas with large snow depth have unrealistically !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. + !-- Based on Sturm et al. keff=0.452 typical for hard snow slabs + !-- with rhosn=488 kg/m^3. Thdifsn = 0.452/(2090*488)=4.431718e-7 !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys + thdifsn = 4.431718e-7_kind_phys else thdifsn = keff/rhocsn * fact endif @@ -5776,27 +5785,16 @@ SUBROUTINE SNOWTEMP( debug_print,xlat,xlon, & RHOCSN=sheatsn* RHOSN if(isncond_opt == 1) then - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then - !-- some areas with large snow depth have unrealistically - !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. - !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys - else - !-- old version thdifsn = 0.265/RHOCSN - THDIFSN = 0.265_kind_phys/RHOCSN - endif + !-- old version thdifsn = 0.265/RHOCSN + THDIFSN = 0.265_kind_phys/RHOCSN else !-- 07Jun19 - thermal conductivity (K_eff) from Sturm et al.(1997) !-- keff = 10. ** (2.650 * RHOSN*1.e-3 - 1.652) fact = one if(rhosn < 156._kind_phys .or. (newsnow > zero .and. rhonewsn < 156._kind_phys)) then keff = 0.023_kind_phys + 0.234_kind_phys * rhosn * 1.e-3_kind_phys - !-- fact is added by tgs based on 4 Jan 2017 testing - fact = 5._kind_phys else keff = 0.138_kind_phys - 1.01_kind_phys * rhosn*1.e-3_kind_phys + 3.233_kind_phys * rhosn**2 * 1.e-6_kind_phys - fact = 2._kind_phys if(debug_print) then print *,'End SNOWTEMP - xlat,xlon,rhosn,keff',xlat,xlon,rhosn,keff print *,'End SNOWTEMP - 0.265/rhocsn',0.265/rhocsn @@ -5807,12 +5805,13 @@ SUBROUTINE SNOWTEMP( debug_print,xlat,xlon, & xlat,xlon,newsnow, rhonewsn,rhosn,fact,keff,keff/rhocsn*fact endif - if(newsnow <= zero .and. snhei > 3.0_kind_phys*SNHEI_crit .and. rhosn > 250._kind_phys) then + if(newsnow <= zero .and. snhei > one .and. rhosn > 250._kind_phys) then !-- some areas with large snow depth have unrealistically !-- low snow density (in the Rockie's with snow depth > 1 m). - !-- Based on Sturm et al. the 2.5e-6 is typical for hard snow slabs. + !-- Based on Sturm et al. keff=0.452 typical for hard snow slabs + !-- with rhosn=488 kg/m^3. Thdifsn = 0.452/(2090*488)=4.431718e-7 !-- In future a better compaction scheme is needed for these areas. - thdifsn = 2.5e-6_kind_phys + thdifsn = 4.431718e-7_kind_phys else thdifsn = keff/rhocsn * fact endif @@ -5959,6 +5958,7 @@ END SUBROUTINE SNOWTEMP !! This subroutine solves moisture budget and computes soil moisture !! and surface and sub-surface runoffs. SUBROUTINE SOILMOIST ( debug_print, & + xlat, xlon, testptlat, testptlon, & DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, & !--- input parameters ZSMAIN,ZSHALF,DIFFU,HYDRO, & QSG,QVG,QCG,QCATM,QVATM,PRCP, & @@ -6012,6 +6012,7 @@ SUBROUTINE SOILMOIST ( debug_print, & !--- input variables LOGICAL, INTENT(IN ) :: debug_print real (kind_phys), INTENT(IN ) :: DELT + real (kind_phys), INTENT(IN ) :: xlat, xlon, testptlat, testptlon INTEGER, INTENT(IN ) :: NZS,NDDZS ! input variables @@ -6099,8 +6100,12 @@ SUBROUTINE SOILMOIST ( debug_print, & DENOM=one+X2+X4-Q2*COSMC(K) COSMC(K+1)=Q4/DENOM IF (debug_print ) THEN - print *,'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' & - ,q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' & + ,q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k + endif ENDIF RHSMC(K+1)=(SOILMOIS(KN)+Q2*RHSMC(K) & +TRANSP(KN) & @@ -6131,8 +6136,12 @@ SUBROUTINE SOILMOIST ( debug_print, & TOTLIQ=PRCP-DRIP/DELT-(one-VEGFRAC)*DEW*RAS-SMELT IF (debug_print ) THEN -print *,'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', & - UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', & + UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT + endif ENDIF FLX=TOTLIQ @@ -6175,7 +6184,7 @@ SUBROUTINE SOILMOIST ( debug_print, & INFMAX1 = zero ENDIF IF (debug_print ) THEN - print *,'INFMAX1 before frozen part',INFMAX1 + print *,'INFMAX1 before frozen part',INFMAX1 ENDIF ! ----------- FROZEN GROUND VERSION -------------------------- @@ -6209,8 +6218,8 @@ SUBROUTINE SOILMOIST ( debug_print, & INFMAX = MAX(INFMAX1,HYDRO(1)*SOILMOIS(1)) INFMAX = MIN(INFMAX, -TOTLIQ) IF (debug_print ) THEN -print *,'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', & - INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ + print *,'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', & + INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ ENDIF !---- IF (-TOTLIQ.GT.INFMAX)THEN @@ -6260,8 +6269,12 @@ SUBROUTINE SOILMOIST ( debug_print, & END IF IF (debug_print ) THEN - print *,'SOILMOIS,SOILIQW, soilice',SOILMOIS,SOILIQW,soilice*riw - print *,'COSMC,RHSMC',COSMC,RHSMC + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'SOILMOIS,SOILIQW, soilice',SOILMOIS,SOILIQW,soilice*riw + print *,'COSMC,RHSMC',COSMC,RHSMC + endif ENDIF !--- FINAL SOLUTION FOR SOILMOIS ! DO K=2,NZS1 @@ -6287,7 +6300,11 @@ SUBROUTINE SOILMOIST ( debug_print, & END IF END DO IF (debug_print ) THEN - print *,'END soilmois,soiliqw,soilice',soilmois,SOILIQW,soilice*riw + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'END soilmois,soiliqw,soilice',soilmois,SOILIQW,soilice*riw + endif ENDIF MAVAIL=max(.00001_kind_phys,min(one,(SOILMOIS(1)/(REF-QMIN)*(one-snowfrac)+one*snowfrac))) @@ -6299,6 +6316,7 @@ END SUBROUTINE SOILMOIST !! This subroutine computes thermal diffusivity, and diffusional and !! hydraulic condeuctivities in soil. SUBROUTINE SOILPROP( debug_print, & + xlat, xlon, testptlat, testptlon, & nzs,fwsat,lwsat,tav,keepfr, & !--- input variables soilmois,soiliqw,soilice, & soilmoism,soiliqwm,soilicem, & @@ -6332,6 +6350,8 @@ SUBROUTINE SOILPROP( debug_print, & !--- soil properties LOGICAL, INTENT(IN ) :: debug_print INTEGER, INTENT(IN ) :: NZS + real (kind_phys), INTENT(IN ) :: xlat, xlon, testptlat, testptlon + real (kind_phys) , & INTENT(IN ) :: RHOCS, & BCLH, & @@ -6508,6 +6528,7 @@ END SUBROUTINE SOILPROP !> This subroutine solves the transpiration function (EQs. 18,19 in !! Smirnova et al.(1997) \cite Smirnova_1997) SUBROUTINE TRANSF( debug_print, & + xlat,xlon,testptlat,testptlon, & nzs,nroot,soiliqw,tabs,lai,gswin, & !--- input variables dqm,qmin,ref,wilt,zshalf,pc,iland, & !--- soil fixed fields tranf,transum) !--- output variables @@ -6528,6 +6549,7 @@ SUBROUTINE TRANSF( debug_print, & LOGICAL, INTENT(IN ) :: debug_print INTEGER, INTENT(IN ) :: nroot,nzs,iland + real (kind_phys), INTENT(IN ) :: xlat,xlon,testptlat,testptlon real (kind_phys) , & INTENT(IN ) :: GSWin, TABS, lai @@ -6574,7 +6596,7 @@ SUBROUTINE TRANSF( debug_print, & ap4=59.656_kind_phys gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4 if(totliq.ge.ref) gx=one - if(totliq.le.zero) gx=zero + if(totliq.le.wilt) gx=zero if(gx.gt.one) gx=one if(gx.lt.zero) gx=zero DID=zshalf(2) @@ -6587,7 +6609,7 @@ SUBROUTINE TRANSF( debug_print, & TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID ENDIF !-- uncomment next line for non-linear root distribution - TRANF(1)=part(1) + !TRANF(1)=part(1) DO K=2,NROOT totliq=soiliqw(k)+qmin @@ -6597,7 +6619,7 @@ SUBROUTINE TRANSF( debug_print, & sm4=sm3*sm1 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4 if(totliq.ge.ref) gx=one - if(totliq.le.zero) gx=zero + if(totliq.le.wilt) gx=zero if(gx.gt.one) gx=one if(gx.lt.zero) gx=zero DID=zshalf(K+1)-zshalf(K) @@ -6611,8 +6633,16 @@ SUBROUTINE TRANSF( debug_print, & /(REF-WILT)*DID ENDIF !-- uncomment next line for non-linear root distribution -! TRANF(k)=part(k) + !TRANF(k)=part(k) END DO + IF (debug_print ) THEN + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'soiliqw =',soiliqw,'wilt=',wilt,'qmin= ',qmin + print *,'tranf = ',tranf + endif + ENDIF ! For LAI> 3 => transpiration at potential rate (F.Tardieu, 2013) if(lai > 4._kind_phys) then @@ -6624,7 +6654,11 @@ SUBROUTINE TRANSF( debug_print, & ! pctot=min(0.8,max(pc,pc*lai)) endif IF ( debug_print ) THEN - print *,'pctot,lai,pc',pctot,lai,pc + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'pctot,lai,pc',pctot,lai,pc + endif ENDIF !--- !--- air temperature function @@ -6634,9 +6668,6 @@ SUBROUTINE TRANSF( debug_print, & ELSE FTEM = one / (one + EXP(0.5_kind_phys * (TABS - 314.0_kind_phys))) ENDIF - IF ( debug_print ) THEN - print *,'tabs,ftem',tabs,ftem - ENDIF !--- incoming solar function cmin = one/rsmax_data cmax = one/rstbl(iland) @@ -6659,27 +6690,33 @@ SUBROUTINE TRANSF( debug_print, & else fsol = one endif - IF ( debug_print ) THEN - print *,'GSWin,lai,f1,fsol',gswin,lai,f1,fsol - ENDIF !--- total conductance totcnd =(cmin + (cmax - cmin)*pctot*ftem*fsol)/cmax IF ( debug_print ) THEN - print *,'iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd' & - ,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'GSWin,Tabs,lai,f1,cmax,cmin,pc,pctot,ftem,fsol',GSWin,Tabs,lai,f1,cmax,cmin,pc,pctot,ftem,fsol + print *,'iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd' & + ,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd + endif ENDIF !-- TRANSUM - total for the rooting zone transum=zero DO K=1,NROOT ! linear root distribution - TRANF(k)=max(cmin,TRANF(k)*totcnd) + TRANF(k)=max(zero,TRANF(k)*totcnd) transum=transum+tranf(k) END DO IF ( debug_print ) THEN - print *,'transum,TRANF',transum,tranf - endif + if (abs(xlat-testptlat).lt.0.05 .and. & + abs(xlon-testptlon).lt.0.05)then + print *,'xlat,xlon=',xlat,xlon + print *,'transum,TRANF',transum,tranf + endif + ENDIF !----------------------------------------------------------------- END SUBROUTINE TRANSF diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index c456e87cd..7b5b83b37 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -329,6 +329,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step, istep, nsteps, & prcp, rain, graupel, ice, snow, sr, & refl_10cm, fullradar_diag, & + max_hail_diam_sfc, & do_radar_ref, aerfld, & mpicomm, mpirank, mpiroot, blkno, & ext_diag, diag3d, reset_diag3d, & @@ -387,6 +388,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & real(kind_phys), intent( out) :: sr(:) ! Radar reflectivity real(kind_phys), intent(inout) :: refl_10cm(:,:) + real(kind_phys), intent(inout) :: max_hail_diam_sfc(:) logical, intent(in ) :: do_radar_ref logical, intent(in) :: sedi_semi integer, intent(in) :: decfl @@ -698,6 +700,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & kme_stoch=kme_stoch, & @@ -738,6 +741,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + max_hail_diam_sfc=max_hail_diam_sfc, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 5918e4dd9..ae1072d39 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -610,6 +610,14 @@ type = real kind = kind_phys intent = out +[max_hail_diam_sfc] + standard_name = max_hail_diameter_sfc + long_name = instantaneous maximum hail diameter at lowest model level + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout [fullradar_diag] standard_name = do_full_radar_reflectivity long_name = flag for computing full radar reflectivity diff --git a/physics/mynnedmf_wrapper.F90 b/physics/mynnedmf_wrapper.F90 index 3c7de235f..487753027 100644 --- a/physics/mynnedmf_wrapper.F90 +++ b/physics/mynnedmf_wrapper.F90 @@ -131,7 +131,8 @@ SUBROUTINE mynnedmf_wrapper_run( & & edmf_a,edmf_w,edmf_qt, & & edmf_thl,edmf_ent,edmf_qc, & & sub_thl,sub_sqv,det_thl,det_sqv,& - & nupdraft,maxMF,ktop_plume, & + & maxwidth,maxMF,ztop_plume, & + & ktop_plume, & & dudt, dvdt, dtdt, & & dqdt_water_vapor, dqdt_liquid_cloud, & ! <=== ntqv, ntcw & dqdt_ice, dqdt_snow, & ! <=== ntiw, ntsw @@ -310,9 +311,9 @@ SUBROUTINE mynnedmf_wrapper_run( & real(kind_phys), dimension(:), intent(out) :: & & ch,dtsfc1,dqsfc1,dusfc1,dvsfc1, & & dtsfci_diag,dqsfci_diag,dusfci_diag,dvsfci_diag, & - & maxMF + & maxMF,maxwidth,ztop_plume integer, dimension(:), intent(inout) :: & - & kpbl,nupdraft,ktop_plume + & kpbl,ktop_plume real(kind_phys), dimension(:), intent(inout) :: & & dusfc_cpl,dvsfc_cpl,dtsfc_cpl,dqsfc_cpl @@ -325,6 +326,7 @@ SUBROUTINE mynnedmf_wrapper_run( & integer :: idtend real(kind_phys), dimension(im) :: dusfci1,dvsfci1,dtsfci1,dqsfci1 real(kind_phys), allocatable :: save_qke_adv(:,:) + real(kind_phys), dimension(levs) :: kzero ! Initialize CCPP error handling variables errmsg = '' @@ -355,6 +357,7 @@ SUBROUTINE mynnedmf_wrapper_run( & !print*,"in MYNN, initflag=",initflag endif + kzero = zero !generic zero array !initialize arrays for test EMIS_ANT_NO = 0. @@ -391,7 +394,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QNI= .true. FLAG_QC = .true. FLAG_QNC= .true. - FLAG_QS = .false. !.true. + FLAG_QS = .true. FLAG_QNWFA= nssl_ccn_on ! ERM: Perhaps could use this field for CCN field? FLAG_QNIFA= .false. FLAG_QNBCA= .false. @@ -400,7 +403,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = 0.0 !qgrs_snow(i,k) + sqs(i,k) = qgrs_snow(i,k) ozone(i,k) = qgrs_ozone(i,k) qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) @@ -418,7 +421,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .false. + FLAG_QS = .true. !pipe it in, but do not mix FLAG_QNC= .true. FLAG_QNWFA= .true. FLAG_QNIFA= .true. @@ -428,7 +431,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = 0. !qgrs_snow(i,k) + sqs(i,k) = qgrs_snow(i,k) qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -441,7 +444,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .false. + FLAG_QS = .true. FLAG_QNC= .true. FLAG_QNWFA= .false. FLAG_QNIFA= .false. @@ -451,7 +454,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = 0. !qgrs_snow(i,k) + sqs(i,k) = qgrs_snow(i,k) qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -464,7 +467,7 @@ SUBROUTINE mynnedmf_wrapper_run( & FLAG_QI = .true. FLAG_QNI= .true. FLAG_QC = .true. - FLAG_QS = .false. + FLAG_QS = .true. FLAG_QNC= .false. FLAG_QNWFA= .false. FLAG_QNIFA= .false. @@ -474,7 +477,7 @@ SUBROUTINE mynnedmf_wrapper_run( & sqv(i,k) = qgrs_water_vapor(i,k) sqc(i,k) = qgrs_liquid_cloud(i,k) sqi(i,k) = qgrs_ice(i,k) - sqs(i,k) = 0. !qgrs_snow(i,k) + sqs(i,k) = qgrs_snow(i,k) qnc(i,k) = 0. qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) @@ -565,7 +568,7 @@ SUBROUTINE mynnedmf_wrapper_run( & call moisture_check2(levs, delt, & delp(i,:), exner(i,:), & sqv(i,:), sqc(i,:), & - sqi(i,:), sqs(i,:), & + sqi(i,:), kzero(:), & t3d(i,:) ) enddo @@ -748,7 +751,7 @@ SUBROUTINE mynnedmf_wrapper_run( & & edmf_thl=edmf_thl,edmf_ent=edmf_ent,edmf_qc=edmf_qc,& !output & sub_thl3D=sub_thl,sub_sqv3D=sub_sqv, & & det_thl3D=det_thl,det_sqv3D=det_sqv, & - & nupdraft=nupdraft,maxMF=maxMF, & !output + & maxwidth=maxwidth,maxMF=maxMF,ztop_plume=ztop_plume,& !output & ktop_plume=ktop_plume, & !output & spp_pbl=spp_pbl,pattern_spp_pbl=spp_wts_pbl, & !input & RTHRATEN=htrlw, & !input @@ -834,7 +837,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) !dqdt_ozone(i,k) = 0.0 dqdt_water_aer_num_conc(i,k) = RQNWFABLTEN(i,k) dqdt_ice_aer_num_conc(i,k) = RQNIFABLTEN(i,k) @@ -869,7 +872,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) enddo enddo if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then @@ -887,7 +890,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_liquid_cloud(i,k) = RQCBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) !dqdt_ozone(i,k) = 0.0 enddo enddo @@ -917,7 +920,7 @@ SUBROUTINE mynnedmf_wrapper_run( & dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k) dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k)) dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k) - !dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) + dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k)) IF ( nssl_ccn_on ) THEN ! dqdt_cccn(i,k) = RQNWFABLTEN(i,k) ENDIF @@ -1005,8 +1008,8 @@ SUBROUTINE mynnedmf_wrapper_run( & print*,"dudt:",dudt(1,1),dudt(1,2),dudt(1,levs) print*,"dvdt:",dvdt(1,1),dvdt(1,2),dvdt(1,levs) print*,"dqdt:",dqdt_water_vapor(1,1),dqdt_water_vapor(1,2),dqdt_water_vapor(1,levs) - print*,"ktop_plume:",ktop_plume(1)," maxmf:",maxmf(1) - print*,"nup:",nupdraft(1) + print*,"ztop_plume:",ztop_plume(1)," maxmf:",maxmf(1) + print*,"maxwidth:",maxwidth(1) print* endif diff --git a/physics/mynnedmf_wrapper.meta b/physics/mynnedmf_wrapper.meta index ec4706aba..8614d3ba2 100644 --- a/physics/mynnedmf_wrapper.meta +++ b/physics/mynnedmf_wrapper.meta @@ -964,13 +964,14 @@ type = real kind = kind_phys intent = inout -[nupdraft] - standard_name = number_of_plumes - long_name = number of plumes per grid column - units = count +[maxwidth] + standard_name = maximum_width_of_plumes + long_name = maximum width of plumes per grid column + units = m dimensions = (horizontal_loop_extent) - type = integer - intent = inout + type = real + kind = kind_phys + intent = out [maxMF] standard_name = maximum_mass_flux long_name = maximum mass flux within a column @@ -979,6 +980,14 @@ type = real kind = kind_phys intent = out +[ztop_plume] + standard_name = height_of_tallest_plume_in_a_column + long_name = height of tallest plume in a column + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [ktop_plume] standard_name = k_level_of_highest_plume long_name = k-level of highest plume diff --git a/physics/sgscloud_radpre.F90 b/physics/sgscloud_radpre.F90 index 44ab87bcc..936393d5b 100644 --- a/physics/sgscloud_radpre.F90 +++ b/physics/sgscloud_radpre.F90 @@ -216,10 +216,10 @@ subroutine sgscloud_radpre_run( & qi(i,k) = ice_frac*qi_bl(i,k) !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b) - if(qi(i,k)>1.E-8)clouds5(i,k)=max(173.45 + 2.14*Tc, 20.) + clouds5(i,k)=max(173.45 + 2.14*Tc, 20.) !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 8b) !iwc = qi(i,k)*1.0e6*rho(i,k) - !IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(139.7 + 1.76*Tc + 13.49*LOG(iwc), 20.) + !clouds5(i,k)=MAX(139.7 + 1.76*Tc + 13.49*LOG(iwc), 20.) !calculate the ice water path using additional BL clouds clouds4(i,k) = max(0.0, qi(i,k) * gfac * delp(i,k)) @@ -229,7 +229,7 @@ subroutine sgscloud_radpre_run( & qs(i,k) = snow_frac*qi_bl(i,k) !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b) - if(qs(i,k)>1.E-8)clouds9(i,k)=max(2.*(173.45 + 2.14*Tc), 50.) + clouds9(i,k)=max(2.*(173.45 + 2.14*Tc), 50.) !calculate the snow water path using additional BL clouds clouds8(i,k) = max(0.0, qs(i,k) * gfac * delp(i,k)) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 7b69fc9e3..c9a6344b8 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -54,7 +54,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, dust_alpha_in, dust_gamma_in, fire_in, & seas_opt_in, dust_opt_in, drydep_opt_in, coarsepm_settling_in, & do_plumerise_in, plumerisefire_frq_in, addsmoke_flag_in, & - wetdep_ls_opt_in,wetdep_ls_alpha_in, & + wetdep_ls_opt_in,wetdep_ls_alpha_in, fire_heat_flux_out, & + frac_grid_burned_out, & smoke_forecast_in, aero_ind_fdb_in,dbg_opt_in,errmsg,errflg) implicit none @@ -88,6 +89,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:), intent(inout) :: coef_bb, fhist real(kind_phys), dimension(:,:), intent(inout) :: ebu_smoke real(kind_phys), dimension(:,:), intent(inout) :: fire_in + real(kind_phys), dimension(:), intent(out) :: fire_heat_flux_out + real(kind_phys), dimension(:), intent(out) :: frac_grid_burned_out real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume real(kind_phys), dimension(:), intent( out) :: hwp real(kind_phys), dimension(:,:), intent(out) :: smoke_ext, dust_ext @@ -339,7 +342,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ! the plumerise is controlled by the namelist option of plumerise_flag if (call_fire) then call ebu_driver ( & - flam_frac,ebu_in,ebu, & + flam_frac,ebu_in,ebu, & t_phy,moist(:,:,:,p_qv), & rho_phy,vvel,u_phy,v_phy,p_phy, & z_at_w,zmid,g,con_cp,con_rd, & @@ -348,8 +351,16 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, errmsg, errflg ) if(errflg/=0) return + do i = its,ite + if ( plume_frp(i,1,p_frp_hr) .ge. 1.E7 ) then + fire_heat_flux_out(i) = min(max(0.,0.88*plume_frp(i,1,p_frp_hr)/0.55/dxy(i,1)) ,50000.) ! JLS - W m-2 [0 - 10,000] + frac_grid_burned_out(i) = min(max(0., 1.3*0.0006*plume_frp(i,1,p_frp_hr)/dxy(i,1) ),1.) + else + fire_heat_flux_out(i) = 0.0 + frac_grid_burned_out(i) = 0.0 + endif + enddo end if - ! -- add biomass burning emissions at every timestep if (addsmoke_flag == 1) then call add_emis_burn(dt,dz8w,rho_phy,rel_hum,chem, & diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index cddc20fbc..7b22b9799 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -740,6 +740,22 @@ dimensions = () type = logical intent = in +[fire_heat_flux_out] + standard_name = surface_fire_heat_flux + long_name = heat flux of fire at the surface + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[frac_grid_burned_out] + standard_name = fraction_of_grid_cell_burning + long_name = ration of the burnt area to the grid cell area + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [dbg_opt_in] standard_name = do_smoke_debug long_name = flag for rrfs smoke plumerise debug