diff --git a/CODEOWNERS b/CODEOWNERS index b6c597371..0d5230f89 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -3,7 +3,7 @@ # These owners will be the default owners for everything in the repo. #* @defunkt -* @DomHeinzeller +* @climbfuji @llpcarson @grantfirl @JulieSchramm # Order is important. The last matching pattern has the most precedence. # So if a pull request only touches javascript files, only these owners diff --git a/physics/GFS_DCNV_generic.F90 b/physics/GFS_DCNV_generic.F90 index bfe97bc70..cbab52377 100644 --- a/physics/GFS_DCNV_generic.F90 +++ b/physics/GFS_DCNV_generic.F90 @@ -169,10 +169,10 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, cs dt3dt(i,k) = dt3dt(i,k) + (gt0(i,k)-save_t(i,k)) * frain du3dt(i,k) = du3dt(i,k) + (gu0(i,k)-save_u(i,k)) * frain dv3dt(i,k) = dv3dt(i,k) + (gv0(i,k)-save_v(i,k)) * frain - -! upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * (con_g*frain) -! dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * (con_g*frain) -! det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * (con_g*frain) + ! convective mass fluxes + upd_mf(i,k) = upd_mf(i,k) + ud_mf(i,k) * (con_g*frain) + dwn_mf(i,k) = dwn_mf(i,k) + dd_mf(i,k) * (con_g*frain) + det_mf(i,k) = det_mf(i,k) + dt_mf(i,k) * (con_g*frain) enddo enddo if(qdiag3d) then diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index 435a80509..672be5ba6 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -89,7 +89,7 @@ subroutine GFS_MP_generic_post_run(im, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_qv, rain0, ice0, snow0, & graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice, & totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, dt3dt, dq3dt, rain_cpl, rainc_cpl, snow_cpl, pwat, & - do_sppt, ca_global, dtdtr, dtdtc, drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & + drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, errmsg, errflg) ! use machine, only: kind_phys @@ -119,9 +119,6 @@ subroutine GFS_MP_generic_post_run(im, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac, real(kind=kind_phys), dimension(:,:), intent(inout) :: dq3dt ! only if ldiag3d and qdiag3d ! Stochastic physics / surface perturbations - logical, intent(in) :: do_sppt, ca_global - real(kind=kind_phys), dimension(im,levs), intent(inout) :: dtdtr - real(kind=kind_phys), dimension(im,levs), intent(in) :: dtdtc real(kind=kind_phys), dimension(im), intent(inout) :: drain_cpl real(kind=kind_phys), dimension(im), intent(inout) :: dsnow_cpl @@ -391,11 +388,6 @@ subroutine GFS_MP_generic_post_run(im, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac, pwat(i) = pwat(i) * onebg enddo - ! Stochastic physics / surface perturbations - if (do_sppt .or. ca_global) then -!--- radiation heating rate - dtdtr(1:im,:) = dtdtr(1:im,:) + dtdtc(1:im,:)*dtf - endif end subroutine GFS_MP_generic_post_run !> @} diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 981f5478d..c42b12c31 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -728,40 +728,6 @@ kind = kind_phys intent = inout optional = F -[do_sppt] - standard_name = flag_for_stochastic_physics_perturbations - long_name = flag for stochastic physics perturbations - units = flag - dimensions = () - type = logical - intent = in - optional = F -[ca_global] - standard_name = flag_for_global_cellular_automata - long_name = switch for global ca - units = flag - dimensions = () - type = logical - intent = in - optional = F -[dtdtr] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_on_physics_time_step - long_name = temp. change due to radiative heating per time step - units = K - dimensions = (horizontal_loop_extent,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[dtdtc] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_assuming_clear_sky - long_name = clear sky radiative (shortwave + longwave) heating rate at current time - units = K s-1 - dimensions = (horizontal_loop_extent,vertical_dimension) - type = real - kind = kind_phys - intent = in - optional = F [drain_cpl] standard_name = tendency_of_lwe_thickness_of_precipitation_amount_for_coupling long_name = change in rain_cpl (coupling_type) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 3e8e987c7..af3f4e147 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -535,7 +535,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%acvt' , Tbd%acvt) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%hpbl' , Tbd%hpbl) if (Model%do_sppt) then - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dtdtr' , Tbd%dtdtr) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dtdtnp' , Tbd%dtdtnp) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dtotprcp' , Tbd%dtotprcp) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dcnvprcp' , Tbd%dcnvprcp) end if diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index 94fc5e36b..0cc6a66b8 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -22,7 +22,7 @@ module GFS_phys_time_vary use h2ointerp, only : read_h2odata, setindxh2o, h2ointerpol use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm - use aerinterp, only : read_aerdata, setindxaer, aerinterpol + use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf use iccn_def, only : ciplin, ccnin, ci_pres use iccninterp, only : read_cidata, setindxci, ciinterpol @@ -166,7 +166,7 @@ subroutine GFS_phys_time_vary_init ( integer, intent(out) :: errflg ! Local variables - integer :: i, j, ix, vegtyp + integer :: i, j, ix, vegtyp, iamin, iamax, jamin, jamax real(kind_phys) :: rsnow !--- Noah MP @@ -182,12 +182,17 @@ subroutine GFS_phys_time_vary_init ( errflg = 0 if (is_initialized) return + iamin=999 + iamax=-999 + jamin=999 + jamax=-999 !$OMP parallel num_threads(nthrds) default(none) & !$OMP shared (me,master,ntoz,h2o_phys,im,nx,ny,idate) & !$OMP shared (xlat_d,xlon_d,imap,jmap,errmsg,errflg) & !$OMP shared (levozp,oz_coeff,oz_pres,ozpl) & !$OMP shared (levh2o,h2o_coeff,h2o_pres,h2opl) & +!$OMP shared (iamin, iamax, jamin, jamax) & !$OMP shared (iaerclm,ntrcaer,aer_nm,iflip,iccn) & !$OMP shared (jindx1_o3,jindx2_o3,ddy_o3,jindx1_h,jindx2_h,ddy_h) & !$OMP shared (jindx1_aer,jindx2_aer,ddy_aer,iindx1_aer,iindx2_aer,ddx_aer) & @@ -300,16 +305,20 @@ subroutine GFS_phys_time_vary_init ( call setindxh2o (im, xlat_d, jindx1_h, jindx2_h, ddy_h) endif -!$OMP section !> - Call setindxaer() to initialize aerosols data +!$OMP section if (iaerclm) then call setindxaer (im, xlat_d, jindx1_aer, & jindx2_aer, ddy_aer, xlon_d, & iindx1_aer, iindx2_aer, ddx_aer, & me, master) + iamin=min(minval(iindx1_aer), iamin) + iamax=max(maxval(iindx2_aer), iamax) + jamin=min(minval(jindx1_aer), jamin) + jamax=max(maxval(jindx2_aer), jamax) endif - !$OMP section + !> - Call setindxci() to initialize IN and CCN data if (iccn == 1) then call setindxci (im, xlat_d, jindx1_ci, & @@ -367,6 +376,10 @@ subroutine GFS_phys_time_vary_init ( !$OMP end sections !$OMP end parallel + if (iaerclm) then + call read_aerdataf (iamin, iamax, jamin, jamax, me,master,iflip, & + idate,errmsg,errflg) + endif if (lsm == lsm_noahmp) then if (all(tvxy < zero)) then @@ -433,9 +446,6 @@ subroutine GFS_phys_time_vary_init ( eahxy(ix) = 2000.0_kind_phys -! eahxy = psfc*qv/(0.622+qv); qv is mixing ratio, converted from sepcific -! humidity specific humidity /(1.0 - specific humidity) - cmxy(ix) = zero chxy(ix) = zero fwetxy(ix) = zero @@ -804,10 +814,10 @@ subroutine GFS_phys_time_vary_timestep_init ( !> - Call ciinterpol() to make IN and CCN data interpolation if (iccn == 1) then - call ciinterpol (me, im, idate, fhour, & - jindx1_ci, jindx2_ci, & - ddy_ci, iindx1_ci, & - iindx2_ci, ddx_ci, & + call ciinterpol (me, im, idate, fhour, & + jindx1_ci, jindx2_ci, & + ddy_ci, iindx1_ci, & + iindx2_ci, ddx_ci, & levs, prsl, in_nm, ccn_nm) endif diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index a1acc3fa0..fb46de2bd 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -2,8 +2,8 @@ !! Contains code related to GFS physics suite setup (physics part of time_vary_step) !>\defgroup mod_GFS_phys_time_vary GFS Physics Time Update -!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor, -!! aerosol, IN&CCN and surface properties updates. +!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor, +!! aerosol, IN&CCN and surface properties updates. !> @{ module GFS_phys_time_vary @@ -23,10 +23,19 @@ module GFS_phys_time_vary use iccn_def, only : ciplin, ccnin, ci_pres use iccninterp, only : read_cidata, setindxci, ciinterpol -#if 0 + use cires_tauamf_data, only: cires_indx_ugwp, read_tau_amf, tau_amf_interp + use cires_tauamf_data, only: tau_limb, days_limb, ugwp_taulat + !--- variables needed for calculating 'sncovr' use namelist_soilveg, only: salp_data, snupx -#endif + use set_soilveg_mod, only: set_soilveg + + ! --- needed for Noah MP init + use noahmp_tables, only: laim_table,saim_table,sla_table, & + bexp_table,smcmax_table,smcwlt_table, & + dwsat_table,dksat_table,psisat_table, & + isurban_table,isbarren_table, & + isice_table,iswater_table implicit none @@ -36,9 +45,13 @@ module GFS_phys_time_vary logical :: is_initialized = .false. - real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_phys), parameter :: con_99 = 99.0_kind_phys - real(kind=kind_phys), parameter :: con_100 = 100.0_kind_phys + real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys + real(kind=kind_phys), parameter :: con_99 = 99.0_kind_phys + real(kind=kind_phys), parameter :: con_100 = 100.0_kind_phys + real(kind=kind_phys), parameter :: missing_value = 9.99e20_kind_phys + real(kind=kind_phys), parameter :: drythresh = 1.e-4_kind_phys + real(kind=kind_phys), parameter :: zero = 0.0_kind_phys + real(kind=kind_phys), parameter :: one = 1.0_kind_phys contains @@ -52,7 +65,14 @@ subroutine GFS_phys_time_vary_init ( jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap, & - errmsg, errflg) + do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau, & + isot, ivegsrc, nlunit, sncovr, sncovr_ice, lsm, lsm_noahmp, lsm_ruc, min_seaice, & + fice, landfrac, vtype, weasd, lsoil, zs, dzs, lsnow_lsm_lbound, lsnow_lsm_ubound, & + tvxy, tgxy, tahxy, canicexy, canliqxy, eahxy, cmxy, chxy, fwetxy, sneqvoxy, alboldxy,& + qsnowxy, wslakexy, albdvis, albdnir, albivis, albinir, emiss, taussxy, waxy, wtxy, & + zwtxy, xlaixy, xsaixy, lfmassxy, stmassxy, rtmassxy, woodxy, stblcpxy, fastcpxy, & + smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy, & + slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, nthrds, errmsg, errflg) implicit none @@ -71,12 +91,85 @@ subroutine GFS_phys_time_vary_init ( integer, intent(inout) :: jindx1_ci(:), jindx2_ci(:), iindx1_ci(:), iindx2_ci(:) real(kind_phys), intent(inout) :: ddy_ci(:), ddx_ci(:) integer, intent(inout) :: imap(:), jmap(:) - + logical, intent(in) :: do_ugwp_v1 + real(kind_phys), intent(inout) :: ddy_j1tau(:), ddy_j2tau(:) + integer, intent(inout) :: jindx1_tau(:), jindx2_tau(:) + + integer, intent(in) :: isot, ivegsrc, nlunit + real(kind_phys), intent(inout) :: sncovr(:), sncovr_ice(:) + integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc + real(kind_phys), intent(in) :: min_seaice, fice(:) + real(kind_phys), intent(in) :: landfrac(:), vtype(:) + real(kind_phys), intent(inout) :: weasd(:) + + ! NoahMP - only allocated when NoahMP is used + integer, intent(in) :: lsoil, lsnow_lsm_lbound, lsnow_lsm_ubound + real(kind_phys), intent(in) :: zs(:) + real(kind_phys), intent(in) :: dzs(:) + real(kind_phys), intent(inout) :: tvxy(:) + real(kind_phys), intent(inout) :: tgxy(:) + real(kind_phys), intent(inout) :: tahxy(:) + real(kind_phys), intent(inout) :: canicexy(:) + real(kind_phys), intent(inout) :: canliqxy(:) + real(kind_phys), intent(inout) :: eahxy(:) + real(kind_phys), intent(inout) :: cmxy(:) + real(kind_phys), intent(inout) :: chxy(:) + real(kind_phys), intent(inout) :: fwetxy(:) + real(kind_phys), intent(inout) :: sneqvoxy(:) + real(kind_phys), intent(inout) :: alboldxy(:) + real(kind_phys), intent(inout) :: qsnowxy(:) + real(kind_phys), intent(inout) :: wslakexy(:) + real(kind_phys), intent(inout) :: albdvis(:) + real(kind_phys), intent(inout) :: albdnir(:) + real(kind_phys), intent(inout) :: albivis(:) + real(kind_phys), intent(inout) :: albinir(:) + real(kind_phys), intent(inout) :: emiss(:) + real(kind_phys), intent(inout) :: taussxy(:) + real(kind_phys), intent(inout) :: waxy(:) + real(kind_phys), intent(inout) :: wtxy(:) + real(kind_phys), intent(inout) :: zwtxy(:) + real(kind_phys), intent(inout) :: xlaixy(:) + real(kind_phys), intent(inout) :: xsaixy(:) + real(kind_phys), intent(inout) :: lfmassxy(:) + real(kind_phys), intent(inout) :: stmassxy(:) + real(kind_phys), intent(inout) :: rtmassxy(:) + real(kind_phys), intent(inout) :: woodxy(:) + real(kind_phys), intent(inout) :: stblcpxy(:) + real(kind_phys), intent(inout) :: fastcpxy(:) + real(kind_phys), intent(inout) :: smcwtdxy(:) + real(kind_phys), intent(inout) :: deeprechxy(:) + real(kind_phys), intent(inout) :: rechxy(:) + real(kind_phys), intent(inout) :: snowxy(:) + real(kind_phys), intent(inout) :: snicexy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: snliqxy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: tsnoxy (:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: smoiseq(:,:) + real(kind_phys), intent(inout) :: zsnsoxy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: slc(:,:) + real(kind_phys), intent(inout) :: smc(:,:) + real(kind_phys), intent(inout) :: stc(:,:) + real(kind_phys), intent(in) :: tsfcl(:) + real(kind_phys), intent(in) :: snowd(:) + real(kind_phys), intent(in) :: canopy(:) + real(kind_phys), intent(in) :: tg3(:) + real(kind_phys), intent(in) :: stype(:) + real(kind_phys), intent(in) :: con_t0c + + integer, intent(in) :: nthrds character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Local variables - integer :: i, j, ix + integer :: i, j, ix, vegtyp + real(kind_phys) :: rsnow + + !--- Noah MP + integer :: soiltyp, isnow, is, imn + real(kind=kind_phys) :: masslai, masssai, snd + real(kind=kind_phys) :: bexp, ddz, smcmax, smcwlt, dwsat, dksat, psisat + + real(kind=kind_phys), dimension(:), allocatable :: dzsno + real(kind=kind_phys), dimension(:), allocatable :: dzsnso ! Initialize CCPP error handling variables errmsg = '' @@ -153,6 +246,14 @@ subroutine GFS_phys_time_vary_init ( ! hardcoded in module iccn_def.F and GFS_typedefs.F90 endif +!> - Call tau_amf dats for ugwp_v1 + if (do_ugwp_v1) then + call read_tau_amf(me, master, errmsg, errflg) + endif + +!> - Initialize soil vegetation (needed for sncovr calculation further down) + call set_soilveg(me, isot, ivegsrc, nlunit) + !> - Call setindxoz() to initialize ozone data if (ntoz > 0) then call setindxoz (im, xlat_d, jindx1_o3, jindx2_o3, ddy_o3) @@ -178,6 +279,12 @@ subroutine GFS_phys_time_vary_init ( iindx1_ci, iindx2_ci, ddx_ci) endif +!> - Call cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that can resolve GWs + if (do_ugwp_v1) then + call cires_indx_ugwp (im, me, master, xlat_d, jindx1_tau, jindx2_tau, & + ddy_j1tau, ddy_j2tau) + endif + !--- initial calculation of maps local ix -> global i and j ix = 0 do j = 1,ny @@ -188,32 +295,315 @@ subroutine GFS_phys_time_vary_init ( enddo enddo -#if 0 - !Calculate sncovr if it was read in but empty (from FV3/io/FV3GFS_io.F90/sfc_prop_restart_read) - ! if (first_time_step) then - ! if (nint(Sfcprop%sncovr(1)) == -9999) then - ! !--- compute sncovr from existing variables - ! !--- code taken directly from read_fix.f - ! do ix = 1, im - ! Sfcprop%sncovr(ix) = 0.0 - ! if (Sfcprop%slmsk(ix) > 0.001) then - ! vegtyp = Sfcprop%vtype(ix) - ! if (vegtyp == 0) vegtyp = 7 - ! rsnow = 0.001*Sfcprop%weasd(ix)/snupx(vegtyp) - ! if (0.001*Sfcprop%weasd(ix) < snupx(vegtyp)) then - ! Sfcprop%sncovr(ix) = 1.0 - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data)) - ! else - ! Sfcprop%sncovr(ix) = 1.0 - ! endif - ! endif - ! enddo - ! ! DH* 20201104: don't forget snocvr_ice for RUC LSM (see FV3GFS_io.F90) - ! endif - ! endif -#endif + !--- if sncovr does not exist in the restart, need to create it + if (all(sncovr < zero)) then + if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' + !--- compute sncovr from existing variables + !--- code taken directly from read_fix.f + sncovr(:) = zero + do ix=1,im + if (landfrac(ix) >= drythresh .or. fice(ix) >= min_seaice) then + vegtyp = vtype(ix) + if (vegtyp == 0) vegtyp = 7 + rsnow = 0.001_kind_phys*weasd(ix)/snupx(vegtyp) + if (0.001_kind_phys*weasd(ix) < snupx(vegtyp)) then + sncovr(ix) = one - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data)) + else + sncovr(ix) = one + endif + endif + enddo + endif + + !--- For RUC LSM: create sncovr_ice from sncovr + if (lsm == lsm_ruc) then + if (all(sncovr_ice < zero)) then + if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' + sncovr_ice(:) = sncovr(:) + endif + endif + + if (lsm == lsm_noahmp) then + if (all(tvxy <= zero)) then + allocate(dzsno (lsnow_lsm_lbound:lsnow_lsm_ubound)) + allocate(dzsnso(lsnow_lsm_lbound:lsoil) ) + dzsno(:) = missing_value + dzsnso(:) = missing_value + + tvxy(:) = missing_value + tgxy(:) = missing_value + tahxy(:) = missing_value + canicexy(:) = missing_value + canliqxy(:) = missing_value + eahxy(:) = missing_value + cmxy(:) = missing_value + chxy(:) = missing_value + fwetxy(:) = missing_value + sneqvoxy(:) = missing_value + alboldxy(:) = missing_value + qsnowxy(:) = missing_value + wslakexy(:) = missing_value + albdvis(:) = missing_value + albdnir(:) = missing_value + albivis(:) = missing_value + albinir(:) = missing_value + emiss(:) = missing_value + taussxy(:) = missing_value + waxy(:) = missing_value + wtxy(:) = missing_value + zwtxy(:) = missing_value + xlaixy(:) = missing_value + xsaixy(:) = missing_value + + lfmassxy(:) = missing_value + stmassxy(:) = missing_value + rtmassxy(:) = missing_value + woodxy(:) = missing_value + stblcpxy(:) = missing_value + fastcpxy(:) = missing_value + smcwtdxy(:) = missing_value + deeprechxy(:) = missing_value + rechxy(:) = missing_value + + snowxy (:) = missing_value + snicexy(:,:) = missing_value + snliqxy(:,:) = missing_value + tsnoxy (:,:) = missing_value + smoiseq(:,:) = missing_value + zsnsoxy(:,:) = missing_value + + do ix=1,im + if (landfrac(ix) >= drythresh) then + tvxy(ix) = tsfcl(ix) + tgxy(ix) = tsfcl(ix) + tahxy(ix) = tsfcl(ix) + + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tvxy(ix) = con_t0c + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tgxy(ix) = con_t0c + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tahxy(ix) = con_t0c + + canicexy(ix) = 0.0_kind_phys + canliqxy(ix) = canopy(ix) + + eahxy(ix) = 2000.0_kind_phys + + cmxy(ix) = zero + chxy(ix) = zero + fwetxy(ix) = zero + sneqvoxy(ix) = weasd(ix) ! mm + alboldxy(ix) = 0.65_kind_phys + qsnowxy(ix) = zero + +! if (srflag(ix) > 0.001) qsnowxy(ix) = tprcp(ix)/dtp + ! already set to 0.0 + wslakexy(ix) = zero + taussxy(ix) = zero + albdvis(ix) = 0.2_kind_phys + albdnir(ix) = 0.2_kind_phys + albivis(ix) = 0.2_kind_phys + albinir(ix) = 0.2_kind_phys + emiss(ix) = 0.95_kind_phys + + + waxy(ix) = 4900.0_kind_phys + wtxy(ix) = waxy(ix) + zwtxy(ix) = (25.0_kind_phys + 2.0_kind_phys) - waxy(ix) / 1000.0_kind_phys / 0.2_kind_phys + + vegtyp = vtype(ix) + if (vegtyp == 0) vegtyp = 7 + imn = idate(2) + + if ((vegtyp == isbarren_table) .or. (vegtyp == isice_table) .or. (vegtyp == isurban_table) .or. (vegtyp == iswater_table)) then + + xlaixy(ix) = zero + xsaixy(ix) = zero + + lfmassxy(ix) = zero + stmassxy(ix) = zero + rtmassxy(ix) = zero + + woodxy (ix) = zero + stblcpxy (ix) = zero + fastcpxy (ix) = zero + + else + + xlaixy(ix) = max(laim_table(vegtyp, imn),0.05_kind_phys) +! xsaixy(ix) = max(saim_table(vegtyp, imn),0.05) + xsaixy(ix) = max(xlaixy(ix)*0.1_kind_phys,0.05_kind_phys) + + masslai = 1000.0_kind_phys / max(sla_table(vegtyp),one) + lfmassxy(ix) = xlaixy(ix)*masslai + masssai = 1000.0_kind_phys / 3.0_kind_phys + stmassxy(ix) = xsaixy(ix)* masssai + + rtmassxy(ix) = 500.0_kind_phys + + woodxy(ix) = 500.0_kind_phys + stblcpxy(ix) = 1000.0_kind_phys + fastcpxy(ix) = 1000.0_kind_phys + + endif ! non urban ... + + if (vegtyp == isice_table) then + do is = 1,lsoil + stc(ix,is) = min(stc(ix,is),min(tg3(ix),263.15_kind_phys)) + smc(ix,is) = one + slc(ix,is) = zero + enddo + endif + + snd = snowd(ix)/1000.0_kind_phys ! go to m from snwdph + + if (weasd(ix) /= zero .and. snd == zero ) then + snd = weasd(ix)/1000.0 + endif + + if (vegtyp == 15) then ! land ice in MODIS/IGBP + if (weasd(ix) < 0.1_kind_phys) then + weasd(ix) = 0.1_kind_phys + snd = 0.01_kind_phys + endif + endif + + if (snd < 0.025_kind_phys ) then + snowxy(ix) = zero + dzsno(-2:0) = zero + elseif (snd >= 0.025_kind_phys .and. snd <= 0.05_kind_phys ) then + snowxy(ix) = -1.0_kind_phys + dzsno(0) = snd + elseif (snd > 0.05_kind_phys .and. snd <= 0.10_kind_phys ) then + snowxy(ix) = -2.0_kind_phys + dzsno(-1) = 0.5_kind_phys*snd + dzsno(0) = 0.5_kind_phys*snd + elseif (snd > 0.10_kind_phys .and. snd <= 0.25_kind_phys ) then + snowxy(ix) = -2.0_kind_phys + dzsno(-1) = 0.05_kind_phys + dzsno(0) = snd - 0.05_kind_phys + elseif (snd > 0.25_kind_phys .and. snd <= 0.45_kind_phys ) then + snowxy(ix) = -3.0_kind_phys + dzsno(-2) = 0.05_kind_phys + dzsno(-1) = 0.5_kind_phys*(snd-0.05_kind_phys) + dzsno(0) = 0.5_kind_phys*(snd-0.05_kind_phys) + elseif (snd > 0.45_kind_phys) then + snowxy(ix) = -3.0_kind_phys + dzsno(-2) = 0.05_kind_phys + dzsno(-1) = 0.20_kind_phys + dzsno(0) = snd - 0.05_kind_phys - 0.20_kind_phys + else + errmsg = 'Error in GFS_phys_time_vary.fv3.F90: Problem with the logic assigning snow layers in Noah MP initialization' + errflg = 1 + return + endif + +! Now we have the snowxy field +! snice + snliq + tsno allocation and compute them from what we have + + tsnoxy(ix,:) = zero + snicexy(ix,:) = zero + snliqxy(ix,:) = zero + zsnsoxy(ix,:) = zero + + isnow = nint(snowxy(ix))+1 ! snowxy <=0.0, dzsno >= 0.0 + + do is = isnow,0 + tsnoxy(ix,is) = tgxy(ix) + snliqxy(ix,is) = zero + snicexy(ix,is) = one * dzsno(is) * weasd(ix)/snd + enddo +! +!zsnsoxy, all negative ? +! + do is = isnow,0 + dzsnso(is) = -dzsno(is) + enddo + + do is = 1,4 + dzsnso(is) = -dzs(is) + enddo +! +! Assign to zsnsoxy +! + zsnsoxy(ix,isnow) = dzsnso(isnow) + do is = isnow+1,4 + zsnsoxy(ix,is) = zsnsoxy(ix,is-1) + dzsnso(is) + enddo +! +! smoiseq +! Init water table related quantities here +! + soiltyp = stype(ix) + if (soiltyp /= 0) then + bexp = bexp_table(soiltyp) + smcmax = smcmax_table(soiltyp) + smcwlt = smcwlt_table(soiltyp) + dwsat = dwsat_table(soiltyp) + dksat = dksat_table(soiltyp) + psisat = -psisat_table(soiltyp) + endif + + if (vegtyp == isurban_table) then + smcmax = 0.45_kind_phys + smcwlt = 0.40_kind_phys + endif + + if ((bexp > zero) .and. (smcmax > zero) .and. (-psisat > zero)) then + do is = 1, lsoil + if ( is == 1 )then + ddz = -zs(is+1) * 0.5_kind_phys + elseif ( is < lsoil ) then + ddz = ( zs(is-1) - zs(is+1) ) * 0.5_kind_phys + else + ddz = zs(is-1) - zs(is) + endif + smoiseq(ix,is) = min(max(find_eq_smc(bexp, dwsat, dksat, ddz, smcmax),1.e-4_kind_phys),smcmax*0.99_kind_phys) + enddo + else ! bexp <= 0.0 + smoiseq(ix,1:4) = smcmax + endif ! end the bexp condition + + smcwtdxy(ix) = smcmax + deeprechxy(ix) = zero + rechxy(ix) = zero + + endif + + enddo ! ix + + deallocate(dzsno) + deallocate(dzsnso) + + endif + endif !if Noah MP cold start ends is_initialized = .true. + contains + +! +! Use newton-raphson method to find eq soil moisture +! + function find_eq_smc(bexp, dwsat, dksat, ddz, smcmax) result(smc) + implicit none + real(kind=kind_phys), intent(in) :: bexp, dwsat, dksat, ddz, smcmax + real(kind=kind_phys) :: smc + real(kind=kind_phys) :: expon, aa, bb, func, dfunc, dx + integer :: iter + ! + expon = bexp + 1. + aa = dwsat / ddz + bb = dksat / smcmax ** expon + smc = 0.5 * smcmax + ! + do iter = 1,100 + func = (smc - smcmax) * aa + bb * smc ** expon + dfunc = aa + bb * expon * smc ** bexp + dx = func / dfunc + smc = smc - dx + if ( abs (dx) < 1.e-6_kind_phys) return + enddo + end function find_eq_smc + end subroutine GFS_phys_time_vary_init !! @} @@ -228,7 +618,8 @@ subroutine GFS_phys_time_vary_timestep_init ( jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, & - imap, jmap, prsl, seed0, rann, errmsg, errflg) + imap, jmap, prsl, seed0, rann, do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,& + tau_amf, errmsg, errflg) implicit none @@ -252,7 +643,11 @@ subroutine GFS_phys_time_vary_timestep_init ( real(kind_phys), intent(in) :: prsl(:,:) integer, intent(in) :: seed0 real(kind_phys), intent(inout) :: rann(:,:) - ! + + logical, intent(in) :: do_ugwp_v1 + integer, intent(in) :: jindx1_tau(:), jindx2_tau(:) + real(kind_phys), intent(in) :: ddy_j1tau(:), ddy_j2tau(:) + real(kind_phys), intent(inout) :: tau_amf(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -347,6 +742,13 @@ subroutine GFS_phys_time_vary_timestep_init ( levs, prsl, in_nm, ccn_nm) endif +!> - Call cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that resolve GW-activ + if (do_ugwp_v1) then + call tau_amf_interp(me, master, im, idate, fhour, & + jindx1_tau, jindx2_tau, & + ddy_j1tau, ddy_j2tau, tau_amf) + endif + ! Not needed for SCM: !> - Call gcycle() to repopulate specific time-varying surface properties for AMIP/forecast runs !if (nscyc > 0) then @@ -423,6 +825,11 @@ subroutine GFS_phys_time_vary_finalize(errmsg, errflg) if (allocated(ccnin) ) deallocate(ccnin) if (allocated(ci_pres) ) deallocate(ci_pres) + ! Deallocate UGWP-input arrays + if (allocated(ugwp_taulat)) deallocate(ugwp_taulat) + if (allocated(tau_limb )) deallocate(tau_limb) + if (allocated(days_limb )) deallocate(days_limb) + is_initialized = .false. end subroutine GFS_phys_time_vary_finalize diff --git a/physics/GFS_phys_time_vary.scm.meta b/physics/GFS_phys_time_vary.scm.meta index cf0b3afbd..1edaa32c8 100644 --- a/physics/GFS_phys_time_vary.scm.meta +++ b/physics/GFS_phys_time_vary.scm.meta @@ -1,7 +1,8 @@ [ccpp-table-properties] name = GFS_phys_time_vary type = scheme - dependencies = aerclm_def.F,aerinterp.F90,h2o_def.f,h2ointerp.f90,iccn_def.F,iccninterp.F90,machine.F,mersenne_twister.f,namelist_soilveg.f,ozinterp.f90,ozne_def.f + dependencies = aerclm_def.F,aerinterp.F90,h2o_def.f,h2ointerp.f90,iccn_def.F,iccninterp.F90,machine.F,mersenne_twister.f + dependencies = namelist_soilveg.f,set_soilveg.f,ozinterp.f90,ozne_def.f,cires_tauamf_data.F90,noahmp_tables.f90 ######################################################################## [ccpp-arg-table] @@ -307,6 +308,641 @@ type = integer intent = inout optional = F +[do_ugwp_v1] + standard_name = flag_for_ugwp_version_1 + long_name = flag to activate ver 1 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[jindx1_tau] + standard_name = lower_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index1 for weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[jindx2_tau] + standard_name = upper_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index2 for weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[ddy_j1tau] + standard_name = latitude_interpolation_weight_complement_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = inout + kind = kind_phys + optional = F +[ddy_j2tau] + standard_name = latitude_interpolation_weight_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = inout + kind = kind_phys + optional = F +[isot] + standard_name = soil_type_dataset_choice + long_name = soil type dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[ivegsrc] + standard_name = vegetation_type_dataset_choice + long_name = land use dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[nlunit] + standard_name = iounit_namelist + long_name = fortran unit number for file opens + units = none + dimensions = () + type = integer + intent = in + optional = F +[sncovr] + standard_name = surface_snow_area_fraction_over_land + long_name = surface snow area fraction + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sncovr_ice] + standard_name = surface_snow_area_fraction_over_ice + long_name = surface snow area fraction over ice + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lsm] + standard_name = flag_for_land_surface_scheme + long_name = flag for land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[lsm_noahmp] + standard_name = flag_for_noahmp_land_surface_scheme + long_name = flag for NOAH MP land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[lsm_ruc] + standard_name = flag_for_ruc_land_surface_scheme + long_name = flag for RUC land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[min_seaice] + standard_name = sea_ice_minimum + long_name = minimum sea ice value + units = frac + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[fice] + standard_name = sea_ice_concentration + long_name = ice fraction over open water + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[landfrac] + standard_name = land_area_fraction + long_name = fraction of horizontal grid area occupied by land + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[vtype] + standard_name = vegetation_type_classification_real + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[weasd] + standard_name = water_equivalent_accumulated_snow_depth + long_name = water equiv of acc snow depth over land and sea ice + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lsoil] + standard_name = soil_vertical_dimension + long_name = number of soil layers + units = count + dimensions = () + type = integer + intent = in + optional = F +[zs] + standard_name = depth_of_soil_levels_for_land_surface_model + long_name = depth of soil levels for land surface model + units = m + dimensions = (soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = in + optional = F +[dzs] + standard_name = thickness_of_soil_levels_for_land_surface_model + long_name = thickness of soil levels for land surface model + units = m + dimensions = (soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = in + optional = F +[lsnow_lsm_lbound] + standard_name = lower_bound_of_snow_vertical_dimension_for_land_surface_model + long_name = lower bound of of snow-related arrays for land surface model + units = count + dimensions = () + type = integer + intent = in + optional = F +[lsnow_lsm_ubound] + standard_name = upper_bound_of_snow_vertical_dimension_for_land_surface_model + long_name = upper bound of of snow-related arrays for land surface model + units = count + dimensions = () + type = integer + intent = in + optional = F +[tvxy] + standard_name = vegetation_temperature + long_name = vegetation temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tgxy] + standard_name = ground_temperature_for_noahmp + long_name = ground temperature for noahmp + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tahxy] + standard_name = canopy_air_temperature + long_name = canopy air temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[canicexy] + standard_name = canopy_intercepted_ice_mass + long_name = canopy intercepted ice mass + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[canliqxy] + standard_name = canopy_intercepted_liquid_water + long_name = canopy intercepted liquid water + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[eahxy] + standard_name = canopy_air_vapor_pressure + long_name = canopy air vapor pressure + units = Pa + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cmxy] + standard_name = surface_drag_coefficient_for_momentum_for_noahmp + long_name = surface drag coefficient for momentum for noahmp + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[chxy] + standard_name = surface_drag_coefficient_for_heat_and_moisture_for_noahmp + long_name = surface exchange coeff heat & moisture for noahmp + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fwetxy] + standard_name = area_fraction_of_wet_canopy + long_name = area fraction of canopy that is wetted/snowed + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sneqvoxy] + standard_name = snow_mass_at_previous_time_step + long_name = snow mass at previous time step + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[alboldxy] + standard_name = snow_albedo_at_previous_time_step + long_name = snow albedo at previous time step + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qsnowxy] + standard_name = snow_precipitation_rate_at_surface + long_name = snow precipitation rate at surface + units = mm s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[wslakexy] + standard_name = lake_water_storage + long_name = lake water storage + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[taussxy] + standard_name = nondimensional_snow_age + long_name = non-dimensional snow age + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[waxy] + standard_name = water_storage_in_aquifer + long_name = water storage in aquifer + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[wtxy] + standard_name = water_storage_in_aquifer_and_saturated_soil + long_name = water storage in aquifer and saturated soil + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[zwtxy] + standard_name = water_table_depth + long_name = water table depth + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[xlaixy] + standard_name = leaf_area_index + long_name = leaf area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[xsaixy] + standard_name = stem_area_index + long_name = stem area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lfmassxy] + standard_name = leaf_mass + long_name = leaf mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stmassxy] + standard_name = stem_mass + long_name = stem mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rtmassxy] + standard_name = fine_root_mass + long_name = fine root mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[woodxy] + standard_name = wood_mass + long_name = wood mass including woody roots + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stblcpxy] + standard_name = slow_soil_pool_mass_content_of_carbon + long_name = stable carbon in deep soil + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fastcpxy] + standard_name = fast_soil_pool_mass_content_of_carbon + long_name = short-lived carbon in shallow soil + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[smcwtdxy] + standard_name = soil_water_content_between_soil_bottom_and_water_table + long_name = soil water content between the bottom of the soil and the water table + units = m3 m-3 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[deeprechxy] + standard_name = water_table_recharge_when_deep + long_name = recharge to or from the water table when deep + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rechxy] + standard_name = water_table_recharge_when_shallow + long_name = recharge to or from the water table when shallow + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albdvis] + standard_name = surface_albedo_direct_visible + long_name = direct surface albedo visible band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albdnir] + standard_name = surface_albedo_direct_NIR + long_name = direct surface albedo NIR band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albivis] + standard_name = surface_albedo_diffuse_visible + long_name = diffuse surface albedo visible band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albinir] + standard_name = surface_albedo_diffuse_NIR + long_name = diffuse surface albedo NIR band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[emiss] + standard_name = surface_emissivity_lsm + long_name = surface emissivity from lsm + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[snowxy] + standard_name = number_of_snow_layers + long_name = number of snow layers + units = count + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[snicexy] + standard_name = snow_layer_ice + long_name = snow layer ice + units = mm + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[snliqxy] + standard_name = snow_layer_liquid_water + long_name = snow layer liquid water + units = mm + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[tsnoxy] + standard_name = snow_temperature + long_name = snow_temperature + units = K + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[smoiseq] + standard_name = equilibrium_soil_water_content + long_name = equilibrium soil water content + units = m3 m-3 + dimensions = (horizontal_dimension,soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[zsnsoxy] + standard_name = layer_bottom_depth_from_snow_surface + long_name = depth from the top of the snow surface at the bottom of the layer + units = m + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[slc] + standard_name = volume_fraction_of_unfrozen_soil_moisture + long_name = liquid soil moisture + units = frac + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[smc] + standard_name = volume_fraction_of_soil_moisture + long_name = total soil moisture + units = frac + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stc] + standard_name = soil_temperature + long_name = soil temperature + units = K + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tsfcl] + standard_name = surface_skin_temperature_over_land + long_name = surface skin temperature over land + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[snowd] + standard_name = surface_snow_thickness_water_equivalent + long_name = water equivalent snow depth + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[canopy] + standard_name = canopy_water_amount + long_name = canopy water amount + units = kg m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[tg3] + standard_name = deep_soil_temperature + long_name = deep soil temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[stype] + standard_name = soil_type_classification_real + long_name = soil type for lsm + units = index + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[con_t0c] + standard_name = temperature_at_zero_celsius + long_name = temperature at 0 degree Celsius + units = K + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[nthrds] + standard_name = omp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -442,7 +1078,7 @@ [nsswr] standard_name = number_of_timesteps_between_shortwave_radiation_calls long_name = number of timesteps between shortwave radiation calls - units = + units = dimensions = () type = integer intent = in @@ -775,6 +1411,57 @@ kind = kind_phys intent = inout optional = F +[do_ugwp_v1] + standard_name = flag_for_ugwp_version_1 + long_name = flag to activate ver 1 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[jindx1_tau] + standard_name = lower_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index1 for weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[jindx2_tau] + standard_name = upper_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index2 for weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[ddy_j1tau] + standard_name = latitude_interpolation_weight_complement_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = in + kind = kind_phys + optional = F +[ddy_j2tau] + standard_name = latitude_interpolation_weight_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = in + kind = kind_phys + optional = F +[tau_amf] + standard_name = absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = ngw_absolute_momentum_flux + units = various + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -813,4 +1500,4 @@ dimensions = () type = integer intent = out - optional = F \ No newline at end of file + optional = F diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index fee9b8815..0c9eaf3f0 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -24,10 +24,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & imp_physics_zhao_carr_pdf, imp_physics_mg, imp_physics_wsm6, & imp_physics_fer_hires, julian, yearlen, lndp_var_list, lsswr, lslwr, & ltaerosol, lgfdlmprad, uni_cld, effr_in, do_mynnedmf, lmfshal, & - lmfdeep2, fhswr, fhlwr, solhr, sup, eps, epsm1, fvirt, & + lmfdeep2, fhswr, fhlwr, solhr, sup, con_eps, epsm1, fvirt, & rog, rocp, con_rd, xlat_d, xlat, xlon, coslat, sinlat, tsfc, slmsk, & - prsi, prsl, prslk, tgrs, sfc_wts, mg_cld, effrr_in, & - cnvw_in, cnvc_in, qgrs, aer_nm, dx, icloud, & !inputs from here and above + prsi, prsl, prslk, tgrs, sfc_wts, mg_cld, effrr_in, pert_clds,sppt_wts,& + sppt_amp, cnvw_in, cnvc_in, qgrs, aer_nm, dx, icloud, & !inputs from here and above coszen, coszdg, effrl_inout, effri_inout, effrs_inout, & clouds1, clouds2, clouds3, clouds4, clouds5, & !in/out from here and above kd, kt, kb, mtopa, mbota, raddt, tsfg, tsfa, de_lgth, alb1d, delp, dz, & !output from here and below @@ -65,7 +65,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & & profsw_type, NBDSW use module_radlw_parameters, only: topflw_type, sfcflw_type, & & proflw_type, NBDLW - use surface_perturbation, only: cdfnor + use surface_perturbation, only: cdfnor,ppfbet ! For Thompson MP use module_mp_thompson, only: calc_effectRad, Nt_c, & @@ -100,10 +100,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & logical, intent(in) :: lsswr, lslwr, ltaerosol, lgfdlmprad, & uni_cld, effr_in, do_mynnedmf, & - lmfshal, lmfdeep2 + lmfshal, lmfdeep2, pert_clds - real(kind=kind_phys), intent(in) :: fhswr, fhlwr, solhr, sup, julian - real(kind=kind_phys), intent(in) :: eps, epsm1, fvirt, rog, rocp, con_rd + real(kind=kind_phys), intent(in) :: fhswr, fhlwr, solhr, sup, julian, sppt_amp + real(kind=kind_phys), intent(in) :: con_eps, epsm1, fvirt, rog, rocp, con_rd real(kind=kind_phys), dimension(:), intent(in) :: xlat_d, xlat, xlon, & coslat, sinlat, tsfc, & @@ -112,7 +112,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real(kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, prslk, & tgrs, sfc_wts, & mg_cld, effrr_in, & - cnvw_in, cnvc_in + cnvw_in, cnvc_in, & + sppt_wts real(kind=kind_phys), dimension(:,:,:), intent(in) :: qgrs, aer_nm @@ -203,6 +204,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real(kind=kind_phys), dimension(im,lm+LTP,NF_VGAS) :: gasvmr real(kind=kind_phys), dimension(im,lm+LTP,NBDSW,NF_AESW) :: faersw real(kind=kind_phys), dimension(im,lm+LTP,NBDLW,NF_AELW) :: faerlw + + ! for stochastic cloud perturbations + real(kind=kind_phys), dimension(im) :: cldp1d + real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz + integer :: iflag integer :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & @@ -294,7 +300,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & !> - Compute relative humidity. es = min( prsl(i,k2), fpvs( tgrs(i,k2) ) ) ! fpvs and prsl in pa - qs = max( QMIN, eps * es / (prsl(i,k2) + epsm1*es) ) + qs = max( QMIN, con_eps * es / (prsl(i,k2) + epsm1*es) ) rhly(i,k1) = max( 0.0, min( 1.0, max(QMIN, qgrs(i,k2,ntqv))/qs ) ) qstl(i,k1) = qs enddo @@ -637,7 +643,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) - rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + rho (i,k) = con_eps*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+con_eps)) orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) @@ -652,7 +658,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) - rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + rho (i,k) = con_eps*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+con_eps)) orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) @@ -942,7 +948,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ccnd(1:IM,1:LMK,1) = ccnd(1:IM,1:LMK,1) + cnvw(1:IM,1:LMK) endif - if (imp_physics == imp_physics_zhao_carr .or. imp_physics == imp_physics_mg) then ! zhao/moorthi's prognostic cloud scheme ! or unified cloud and/or with MG microphysics @@ -1048,6 +1053,31 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ! endif ! end_if_ntcw +! perturb cld cover + if (pert_clds) then + do i=1,im + tmp_wt= -1*log( ( 2.0 / ( sppt_wts(i,38) ) ) - 1 ) + call cdfnor(tmp_wt,cdfz) + cldp1d(i) = cdfz + enddo + do k = 1, LMK + do i = 1, IM + ! compute beta distribution parameters + m = clouds(i,k,1) + if (m<0.99 .AND. m > 0.01) then + s = sppt_amp*m*(1.-m) + alpha0 = m*m*(1.-m)/(s*s)-m + beta0 = alpha0*(1.-m)/m + ! compute beta distribution value corresponding + ! to the given percentile albPpert to use as new albedo + call ppfbet(cldp1d(i),alpha0,beta0,iflag,cldtmp) + clouds(i,k,1) = cldtmp + else + clouds(i,k,1) = m + endif + enddo ! end_do_i_loop + enddo ! end_do_k_loop + endif do k = 1, LMK do i = 1, IM clouds1(i,k) = clouds(i,k,1) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index eaa878ee7..e26cdeac1 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -395,6 +395,32 @@ type = logical intent = in optional = F +[sppt_wts] + standard_name = weights_for_stochastic_sppt_perturbation + long_name = weights for stochastic sppt perturbation + units = none + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sppt_amp] + standard_name = total_ampltiude_of_sppt_perturbation + long_name = toal ampltidue of stochastic sppt perturbation + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[pert_clds] + standard_name = flag_for_stochastic_cloud_fraction_perturbations + long_name = flag for stochastic cloud fraction physics perturbations + units = flag + dimensions = () + type = logical + intent = in + optional = F [do_mynnedmf] standard_name = do_mynnedmf long_name = flag to activate MYNN-EDMF @@ -455,7 +481,7 @@ kind = kind_phys intent = in optional = F -[eps] +[con_eps] standard_name = ratio_of_dry_air_to_water_vapor_gas_constants long_name = rd/rv units = none diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index c6661d948..b54f27d65 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -41,7 +41,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do i_cldliq, i_cldice, i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, i_cldliq_nc, & i_cldice_nc, i_twa, effr_in, p_lev, p_lay, tv_lay, t_lay, effrin_cldliq, & effrin_cldice, effrin_cldsnow, tracer, qs_lay, q_lay, relhum, cld_frac_mg, con_g, & - con_rd, uni_cld, lmfshal, lmfdeep2, ltaerosol, do_mynnedmf, imfdeepcnv, & + con_rd, con_eps, uni_cld, lmfshal, lmfdeep2, ltaerosol, do_mynnedmf, imfdeepcnv, & imfdeepcnv_gf, doGP_cldoptics_PADE, doGP_cldoptics_LUT, & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, & cld_rerain, precip_frac, errmsg, errflg) @@ -76,7 +76,8 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do doGP_cldoptics_PADE ! (PADE approximation) real(kind_phys), intent(in) :: & con_g, & ! Physical constant: gravitational constant - con_rd ! Physical constant: gas-constant for dry air + con_rd, & ! Physical constant: gas-constant for dry air + con_eps ! Physical constant: gas constant air / gas constant H2O real(kind_phys), dimension(nCol,nLev), intent(in) :: & tv_lay, & ! Virtual temperature (K) @@ -155,7 +156,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do do iLay = 1, nLev do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) - rho = 0.622*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+0.622)) + rho = con_eps*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+con_eps)) orho = 1./rho qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.meta b/physics/GFS_rrtmgp_thompsonmp_pre.meta index 90ec59760..54d266b67 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.meta +++ b/physics/GFS_rrtmgp_thompsonmp_pre.meta @@ -171,7 +171,7 @@ standard_name = effective_radius_of_stratiform_cloud_liquid_water_particle_in_um long_name = eff. radius of cloud liquid water particle in micrometer units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -180,7 +180,7 @@ standard_name = effective_radius_of_stratiform_cloud_ice_particle_in_um long_name = eff. radius of cloud ice water particle in micrometer units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -189,7 +189,7 @@ standard_name = effective_radius_of_stratiform_cloud_snow_particle_in_um long_name = effective radius of cloud snow particle in micrometers units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -247,7 +247,7 @@ standard_name = air_pressure_at_interface_for_RRTMGP_in_hPa long_name = air pressure at vertical interface for radiation calculation units = hPa - dimensions = (horizontal_dimension,vertical_dimension_plus_one) + dimensions = (horizontal_loop_extent,vertical_dimension_plus_one) type = real kind = kind_phys intent = in @@ -256,7 +256,7 @@ standard_name = air_pressure_at_layer_for_RRTMGP_in_hPa long_name = air pressure at vertical layer for radiation calculation units = hPa - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -265,7 +265,7 @@ standard_name = virtual_temperature long_name = layer virtual temperature units = K - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -283,7 +283,7 @@ standard_name = saturation_vapor_pressure long_name = saturation vapor pressure units = Pa - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -292,7 +292,7 @@ standard_name = water_vapor_mixing_ratio long_name = water vaport mixing ratio units = kg/kg - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -301,7 +301,7 @@ standard_name = relative_humidity long_name = layer relative humidity units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -310,7 +310,7 @@ standard_name = chemical_tracers long_name = chemical tracers units = g g-1 - dimensions = (horizontal_dimension,vertical_dimension,number_of_tracers) + dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers) type = real kind = kind_phys intent = in @@ -330,14 +330,23 @@ units = J kg-1 K-1 dimensions = () type = real - kind = kind_phys + kind = kind_phys + intent = in + optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys intent = in - optional = F + optional = F [cld_frac] standard_name = total_cloud_fraction long_name = layer total cloud fraction units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -346,7 +355,7 @@ standard_name = cloud_liquid_water_path long_name = layer cloud liquid water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -355,7 +364,7 @@ standard_name = mean_effective_radius_for_liquid_cloud long_name = mean effective radius for liquid cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -364,7 +373,7 @@ standard_name = cloud_ice_water_path long_name = layer cloud ice water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -373,7 +382,7 @@ standard_name = mean_effective_radius_for_ice_cloud long_name = mean effective radius for ice cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -382,7 +391,7 @@ standard_name = cloud_snow_water_path long_name = layer cloud snow water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -391,7 +400,7 @@ standard_name = mean_effective_radius_for_snow_flake long_name = mean effective radius for snow cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -400,7 +409,7 @@ standard_name = cloud_rain_water_path long_name = layer cloud rain water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -409,7 +418,7 @@ standard_name = mean_effective_radius_for_rain_drop long_name = mean effective radius for rain cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -418,11 +427,11 @@ standard_name = precipitation_fraction_by_layer long_name = precipitation fraction in each layer units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout - optional = F + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/GFS_stochastics.F90 b/physics/GFS_stochastics.F90 index 9b4533cf9..b3dda08da 100644 --- a/physics/GFS_stochastics.F90 +++ b/physics/GFS_stochastics.F90 @@ -26,15 +26,18 @@ end subroutine GFS_stochastics_finalize !! -# defines random seed indices for radiation (in a reproducible way) !! -# interpolates coefficients for prognostic ozone calculation !! -# performs surface data cycling via the GFS gcycle routine - subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, & - do_skeb, do_ca,ca_global,ca1,si,vfact_ca, & + subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblck, & + do_shum ,do_skeb, do_ca,ca_global,ca1,si,vfact_ca, & zmtnblck, sppt_wts, skebu_wts, skebv_wts, shum_wts,& sppt_wts_inv, skebu_wts_inv, skebv_wts_inv, & - shum_wts_inv, diss_est, & - ugrs, vgrs, tgrs, qgrs, gu0, gv0, gt0, gq0, dtdtr, & + shum_wts_inv, diss_est, ugrs, vgrs, tgrs, qgrs_wv, & + qgrs_cw, qgrs_rw, qgrs_sw, qgrs_iw, qgrs_gl, & + gu0, gv0, gt0, gq0_wv, dtdtnp, & + gq0_cw, gq0_rw, gq0_sw, gq0_iw, gq0_gl, & rain, rainc, tprcp, totprcp, cnvprcp, & totprcpb, cnvprcpb, cplflx, & rain_cpl, snow_cpl, drain_cpl, dsnow_cpl, & + ntcw,ntrw,ntsw,ntiw,ntgl, & errmsg, errflg) use machine, only: kind_phys @@ -44,7 +47,9 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, integer, intent(in) :: im integer, intent(in) :: km integer, intent(in) :: kdt + real(kind_phys), intent(in) :: delt logical, intent(in) :: do_sppt + logical, intent(in) :: pert_mp logical, intent(in) :: do_ca logical, intent(in) :: ca_global logical, intent(in) :: use_zmtnblck @@ -67,13 +72,27 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, real(kind_phys), dimension(1:im,1:km), intent(in) :: ugrs real(kind_phys), dimension(1:im,1:km), intent(in) :: vgrs real(kind_phys), dimension(1:im,1:km), intent(in) :: tgrs - real(kind_phys), dimension(1:im,1:km), intent(in) :: qgrs + real(kind_phys), dimension(1:im,1:km), intent(in) :: qgrs_wv + real(kind_phys), dimension(:,:), intent(in) :: qgrs_cw + real(kind_phys), dimension(:,:), intent(in) :: qgrs_rw + real(kind_phys), dimension(:,:), intent(in) :: qgrs_sw + real(kind_phys), dimension(:,:), intent(in) :: qgrs_iw + real(kind_phys), dimension(:,:), intent(in) :: qgrs_gl real(kind_phys), dimension(1:im,1:km), intent(inout) :: gu0 real(kind_phys), dimension(1:im,1:km), intent(inout) :: gv0 real(kind_phys), dimension(1:im,1:km), intent(inout) :: gt0 - real(kind_phys), dimension(1:im,1:km), intent(inout) :: gq0 - ! dtdtr only allocated if do_sppt == .true. - real(kind_phys), dimension(:,:), intent(in) :: dtdtr + real(kind_phys), dimension(1:im,1:km), intent(inout) :: gq0_wv + real(kind_phys), dimension(:,:), intent(inout) :: gq0_cw + real(kind_phys), dimension(:,:), intent(inout) :: gq0_rw + real(kind_phys), dimension(:,:), intent(inout) :: gq0_sw + real(kind_phys), dimension(:,:), intent(inout) :: gq0_iw + real(kind_phys), dimension(:,:), intent(inout) :: gq0_gl + integer, intent(in) :: ntcw + integer, intent(in) :: ntrw + integer, intent(in) :: ntsw + integer, intent(in) :: ntiw + integer, intent(in) :: ntgl + real(kind_phys), dimension(:,:), intent(inout) :: dtdtnp real(kind_phys), dimension(1:im), intent(in) :: rain real(kind_phys), dimension(1:im), intent(in) :: rainc real(kind_phys), dimension(1:im), intent(inout) :: tprcp @@ -130,17 +149,59 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, upert = (gu0(i,k) - ugrs(i,k)) * sppt_wts(i,k) vpert = (gv0(i,k) - vgrs(i,k)) * sppt_wts(i,k) - tpert = (gt0(i,k) - tgrs(i,k) - dtdtr(i,k)) * sppt_wts(i,k) - qpert = (gq0(i,k) - qgrs(i,k)) * sppt_wts(i,k) + tpert = (gt0(i,k) - tgrs(i,k) - (delt*dtdtnp(i,k))) * sppt_wts(i,k) + qpert = (gq0_wv(i,k) - qgrs_wv(i,k)) * sppt_wts(i,k) gu0(i,k) = ugrs(i,k)+upert gv0(i,k) = vgrs(i,k)+vpert !negative humidity check - qnew = qgrs(i,k)+qpert + qnew = qgrs_wv(i,k)+qpert if (qnew >= 1.0e-10) then - gq0(i,k) = qnew - gt0(i,k) = tgrs(i,k) + tpert + dtdtr(i,k) + gq0_wv(i,k) = qnew + gt0(i,k) = tgrs(i,k) + tpert + (delt*dtdtnp(i,k)) + endif + if (pert_mp) then + if (ntcw>0) then + qpert = (gq0_cw(i,k) - qgrs_cw(i,k)) * sppt_wts(i,k) + qnew = qgrs_cw(i,k)+qpert + gq0_cw(i,k) = qnew + if (qnew < 0.0) then + gq0_cw(i,k) = 0.0 + endif + endif + if (ntrw>0) then + qpert = (gq0_rw(i,k) - qgrs_rw(i,k)) * sppt_wts(i,k) + qnew = qgrs_rw(i,k)+qpert + gq0_rw(i,k) = qnew + if (qnew < 0.0) then + gq0_rw(i,k) = 0.0 + endif + endif + if (ntsw>0) then + qpert = (gq0_sw(i,k) - qgrs_sw(i,k)) * sppt_wts(i,k) + qnew = qgrs_sw(i,k)+qpert + gq0_sw(i,k) = qnew + if (qnew < 0.0) then + gq0_sw(i,k) = 0.0 + endif + endif + if (ntiw>0) then + qpert = (gq0_iw(i,k) - qgrs_iw(i,k)) * sppt_wts(i,k) + qnew = qgrs_iw(i,k)+qpert + gq0_iw(i,k) = qnew + if (qnew < 0.0) then + gq0_iw(i,k) = 0.0 + endif + endif + if (ntgl>0) then + qpert = (gq0_gl(i,k) - qgrs_gl(i,k)) * sppt_wts(i,k) + qnew = qgrs_gl(i,k)+qpert + gq0_gl(i,k) = qnew + if (qnew < 0.0) then + gq0_gl(i,k) = 0.0 + endif + endif endif enddo enddo @@ -154,10 +215,12 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, totprcpb(:) = totprcpb(:) + (sppt_wts(:,15) - 1 )*rain(:) cnvprcpb(:) = cnvprcpb(:) + (sppt_wts(:,15) - 1 )*rainc(:) - if (cplflx) then + if (cplflx) then rain_cpl(:) = rain_cpl(:) + (sppt_wts(:,15) - 1.0)*drain_cpl(:) snow_cpl(:) = snow_cpl(:) + (sppt_wts(:,15) - 1.0)*dsnow_cpl(:) - endif + endif + !zero out radiative heating tendency for next physics step + dtdtnp(:,:)=0.0 endif @@ -201,15 +264,57 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, upert = (gu0(i,k) - ugrs(i,k)) * ca(i,k) vpert = (gv0(i,k) - vgrs(i,k)) * ca(i,k) - tpert = (gt0(i,k) - tgrs(i,k) - dtdtr(i,k)) * ca(i,k) - qpert = (gq0(i,k) - qgrs(i,k)) * ca(i,k) + tpert = (gt0(i,k) - tgrs(i,k) - (delt*dtdtnp(i,k))) * ca(i,k) + qpert = (gq0_wv(i,k) - qgrs_wv(i,k)) * ca(i,k) gu0(i,k) = ugrs(i,k)+upert gv0(i,k) = vgrs(i,k)+vpert !negative humidity check - qnew = qgrs(i,k)+qpert + qnew = qgrs_wv(i,k)+qpert if (qnew >= 1.0e-10) then - gq0(i,k) = qnew - gt0(i,k) = tgrs(i,k) + tpert + dtdtr(i,k) + gq0_wv(i,k) = qnew + gt0(i,k) = tgrs(i,k) + tpert + (delt*dtdtnp(i,k)) + endif + if (pert_mp) then + if (ntcw>0) then + qpert = (gq0_cw(i,k) - qgrs_cw(i,k)) * ca(i,k) + qnew = qgrs_cw(i,k)+qpert + gq0_cw(i,k) = qnew + if (qnew < 0.0) then + gq0_cw(i,k) = 0.0 + endif + endif + if (ntrw>0) then + qpert = (gq0_rw(i,k) - qgrs_rw(i,k)) * ca(i,k) + qnew = qgrs_rw(i,k)+qpert + gq0_rw(i,k) = qnew + if (qnew < 0.0) then + gq0_rw(i,k) = 0.0 + endif + endif + if (ntsw>0) then + qpert = (gq0_sw(i,k) - qgrs_sw(i,k)) * ca(i,k) + qnew = qgrs_sw(i,k)+qpert + gq0_sw(i,k) = qnew + if (qnew < 0.0) then + gq0_sw(i,k) = 0.0 + endif + endif + if (ntiw>0) then + qpert = (gq0_iw(i,k) - qgrs_iw(i,k)) * ca(i,k) + qnew = qgrs_iw(i,k)+qpert + gq0_iw(i,k) = qnew + if (qnew < 0.0) then + gq0_iw(i,k) = 0.0 + endif + endif + if (ntgl>0) then + qpert = (gq0_gl(i,k) - qgrs_gl(i,k)) * ca(i,k) + qnew = qgrs_gl(i,k)+qpert + gq0_gl(i,k) = qnew + if (qnew < 0.0) then + gq0_gl(i,k) = 0.0 + endif + endif endif enddo enddo @@ -227,13 +332,15 @@ subroutine GFS_stochastics_run (im, km, kdt, do_sppt, use_zmtnblck, do_shum, rain_cpl(:) = rain_cpl(:) + (ca(:,15) - 1.0)*drain_cpl(:) snow_cpl(:) = snow_cpl(:) + (ca(:,15) - 1.0)*dsnow_cpl(:) endif + !zero out radiative heating tendency for next physics step + dtdtnp(:,:)=0.0 endif if (do_shum) then do k=1,km - gq0(:,k) = gq0(:,k)*(1.0 + shum_wts(:,k)) + gq0_wv(:,k) = gq0_wv(:,k)*(1.0 + shum_wts(:,k)) shum_wts_inv(:,k) = shum_wts(:,k) end do endif diff --git a/physics/GFS_stochastics.meta b/physics/GFS_stochastics.meta index 43c7b2d42..26dcfd47e 100644 --- a/physics/GFS_stochastics.meta +++ b/physics/GFS_stochastics.meta @@ -31,6 +31,55 @@ type = integer intent = in optional = F +[delt] + standard_name = time_step_for_physics + long_name = physics timestep + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[ntcw] + standard_name = index_for_liquid_cloud_condensate + long_name = tracer index for cloud condensate (or liquid water) + units = index + dimensions = () + type = integer + intent = in + optional = F +[ntiw] + standard_name = index_for_ice_cloud_condensate + long_name = tracer index for ice water + units = index + dimensions = () + intent = in + optional = F + type = integer +[ntrw] + standard_name = index_for_rain_water + long_name = tracer index for rain water + units = index + dimensions = () + type = integer + intent = in + optional = F +[ntsw] + standard_name = index_for_snow_water + long_name = tracer index for snow water + units = index + dimensions = () + type = integer + intent = in + optional = F +[ntgl] + standard_name = index_for_graupel + long_name = tracer index for graupel + units = index + dimensions = () + type = integer + intent = in + optional = F [do_sppt] standard_name = flag_for_stochastic_physics_perturbations long_name = flag for stochastic physics perturbations @@ -39,6 +88,14 @@ type = logical intent = in optional = F +[pert_mp] + standard_name = flag_for_stochastic_microphysics_perturbations + long_name = flag for stochastic microphysics physics perturbations + units = flag + dimensions = () + type = logical + intent = in + optional = F [use_zmtnblck] standard_name = flag_for_mountain_blocking long_name = flag for mountain blocking @@ -223,7 +280,7 @@ kind = kind_phys intent = in optional = F -[qgrs] +[qgrs_wv] standard_name = water_vapor_specific_humidity long_name = water vapor specific humidity units = kg kg-1 @@ -232,6 +289,51 @@ kind = kind_phys intent = in optional = F +[qgrs_cw] + standard_name = cloud_condensed_water_mixing_ratio + long_name = ratio of mass of cloud water to mass of dry air plus vapor (without condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qgrs_rw] + standard_name = rain_water_mixing_ratio + long_name = moist mixing ratio of rain + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qgrs_iw] + standard_name = ice_water_mixing_ratio + long_name = moist mixing ratio of cloud ice + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qgrs_sw] + standard_name = snow_water_mixing_ratio + long_name = moist mixing ratio of snow + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qgrs_gl] + standard_name = graupel_mixing_ratio + long_name = moist ratio of mass of graupel to mass of dry air plus vapor (without condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F [gu0] standard_name = x_wind_updated_by_physics long_name = zonal wind updated by physics @@ -259,7 +361,7 @@ kind = kind_phys intent = inout optional = F -[gq0] +[gq0_wv] standard_name = water_vapor_specific_humidity_updated_by_physics long_name = water vapor specific humidity updated by physics units = kg kg-1 @@ -268,10 +370,55 @@ kind = kind_phys intent = inout optional = F -[dtdtr] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_on_physics_time_step - long_name = temp. change due to radiative heating per time step - units = K +[gq0_cw] + standard_name = cloud_condensed_water_mixing_ratio_updated_by_physics + long_name = cloud condensed water mixing ratio updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[gq0_rw] + standard_name = rain_water_mixing_ratio_updated_by_physics + long_name = moist mixing ratio of rain updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[gq0_iw] + standard_name = ice_water_mixing_ratio_updated_by_physics + long_name = moist mixing ratio of cloud ice updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[gq0_sw] + standard_name = snow_water_mixing_ratio_updated_by_physics + long_name = moist mixing ratio of snow updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[gq0_gl] + standard_name = graupel_mixing_ratio_updated_by_physics + long_name = moist ratio of mass of graupel to mass of dry air plus vapor (without condensates) updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[dtdtnp] + standard_name = tendency_of_air_temperature_to_withold_from_sppt + long_name = temp. change from physics that should not be perturbed by sppt + units = K s-1 dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index b51b6f9c8..93106d2de 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -87,7 +87,7 @@ end subroutine GFS_suite_interstitial_1_finalize !! \htmlinclude GFS_suite_interstitial_1_run.html !! subroutine GFS_suite_interstitial_1_run (im, levs, ntrac, dtf, dtp, slmsk, area, dxmin, dxinv, pgr, & - islmsk, work1, work2, psurf, dudt, dvdt, dtdt, dtdtc, dqdt, errmsg, errflg) + islmsk, work1, work2, psurf, dudt, dvdt, dtdt, dqdt, errmsg, errflg) use machine, only: kind_phys @@ -100,7 +100,7 @@ subroutine GFS_suite_interstitial_1_run (im, levs, ntrac, dtf, dtp, slmsk, area, integer, intent(out), dimension(im) :: islmsk real(kind=kind_phys), intent(out), dimension(im) :: work1, work2, psurf - real(kind=kind_phys), intent(out), dimension(im,levs) :: dudt, dvdt, dtdt, dtdtc + real(kind=kind_phys), intent(out), dimension(im,levs) :: dudt, dvdt, dtdt real(kind=kind_phys), intent(out), dimension(im,levs,ntrac) :: dqdt real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys character(len=*), intent(out) :: errmsg @@ -127,7 +127,6 @@ subroutine GFS_suite_interstitial_1_run (im, levs, ntrac, dtf, dtp, slmsk, area, dudt(i,k) = zero dvdt(i,k) = zero dtdt(i,k) = zero - dtdtc(i,k) = zero enddo enddo do n=1,ntrac @@ -156,11 +155,10 @@ end subroutine GFS_suite_interstitial_2_init subroutine GFS_suite_interstitial_2_finalize() end subroutine GFS_suite_interstitial_2_finalize -#if 0 + !> \section arg_table_GFS_suite_interstitial_2_run Argument Table !! \htmlinclude GFS_suite_interstitial_2_run.html !! -#endif subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplflx, flag_cice, shal_cnv, old_monin, mstrat, & do_shoc, frac_grid, imfshalcnv, dtf, xcosz, adjsfcdsw, adjsfcdlw, cice, pgr, ulwsfc_cice, lwhd, htrsw, htrlw, xmu, ctei_rm, & work1, work2, prsi, tgrs, prsl, qgrs_water_vapor, qgrs_cloud_water, cp, hvap, prslk, suntim, adjsfculw, adjsfculw_lnd, & @@ -476,11 +474,9 @@ end subroutine GFS_suite_interstitial_3_init subroutine GFS_suite_interstitial_3_finalize() end subroutine GFS_suite_interstitial_3_finalize -#if 0 !> \section arg_table_GFS_suite_interstitial_3_run Argument Table !! \htmlinclude GFS_suite_interstitial_3_run.html !! -#endif subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & @@ -514,7 +510,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & real(kind=kind_phys), dimension(im, levs), intent(inout) :: rhc, save_qc ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), dimension(:, :), intent(inout) :: save_qi - real(kind=kind_phys), dimension(:, :), intent(inout) :: save_tcp ! ONLY ALLOCATE FOR THOMPSON! TODO + real(kind=kind_phys), dimension(:, :), intent(inout) :: save_tcp real(kind=kind_phys), dimension(im, levs, nn), intent(inout) :: clw character(len=*), intent(out) :: errmsg @@ -653,8 +649,8 @@ end subroutine GFS_suite_interstitial_4_finalize !! subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, dtf, save_qc, save_qi, con_pi, & - gq0, clw, prsl, save_tcp, con_rd, nwfa, spechum, dqdti, errmsg, errflg) + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, & + gq0, clw, prsl, save_tcp, con_rd, con_eps, nwfa, spechum, dqdti, errmsg, errflg) use machine, only: kind_phys use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber @@ -667,7 +663,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf - logical, intent(in) :: ltaerosol, cplchm + logical, intent(in) :: ltaerosol, cplchm, convert_dry_rho real(kind=kind_phys), intent(in) :: con_pi, dtf real(kind=kind_phys), dimension(im,levs), intent(in) :: save_qc @@ -677,7 +673,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to real(kind=kind_phys), dimension(im,levs,ntrac), intent(inout) :: gq0 real(kind=kind_phys), dimension(im,levs,nn), intent(inout) :: clw real(kind=kind_phys), dimension(im,levs), intent(in) :: prsl - real(kind=kind_phys), intent(in) :: con_rd + real(kind=kind_phys), intent(in) :: con_rd, con_eps real(kind=kind_phys), dimension(:,:), intent(in) :: nwfa, save_tcp real(kind=kind_phys), dimension(im,levs), intent(in) :: spechum @@ -740,33 +736,55 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to enddo if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then - do k=1,levs - do i=1,im - !> - Convert specific humidity to dry mixing ratio - qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) - !> - Density of air in kg m-3 and inverse density - rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) - orho = one/rho - if (ntlnc>0) then - !> - Convert moist mixing ratio to dry mixing ratio - qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) - !> - Convert number concentration from moist to dry - nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) - !> - Convert number concentrations from dry to moist - gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) - endif - if (ntinc>0) then - !> - Convert moist mixing ratio to dry mixing ratio - qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) - !> - Convert number concentration from moist to dry - ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) - !> - Convert number concentrations from dry to moist - gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) - endif + if_convert_dry_rho: if (convert_dry_rho) then + do k=1,levs + do i=1,im + !> - Convert specific humidity to dry mixing ratio + qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) + !> - Density of air in kg m-3 and inverse density + rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+con_eps)) + orho = one/rho + if (ntlnc>0) then + !> - Convert moist mixing ratio to dry mixing ratio + qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) + !> - Convert number concentration from moist to dry + nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + !> - Convert number concentrations from dry to moist + gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) + endif + if (ntinc>0) then + !> - Convert moist mixing ratio to dry mixing ratio + qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) + !> - Convert number concentration from moist to dry + ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + !> - Convert number concentrations from dry to moist + gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) + endif + enddo enddo - enddo + else + do k=1,levs + do i=1,im + !> - Density of air in kg m-3 and inverse density + rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(spechum(i,k)+con_eps)) + orho = one/rho + if (ntlnc>0) then + !> - Update cloud water mixing ratio + qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) + !> - Update cloud water number concentration + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif + if (ntinc>0) then + !> - Update cloud ice mixing ratio + qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) + !> - Update cloud ice number concentration + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif + enddo + enddo + end if if_convert_dry_rho endif else diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index fdf1716f1..5b4b0dbf9 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -244,15 +244,6 @@ kind = kind_phys intent = out optional = F -[dtdtc] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_assuming_clear_sky - long_name = clear sky radiative (shortwave + longwave) heating rate at current time - units = K s-1 - dimensions = (horizontal_loop_extent,vertical_dimension) - type = real - kind = kind_phys - intent = out - optional = F [dqdt] standard_name = tendency_of_tracers_due_to_model_physics long_name = updated tendency of the tracers @@ -1778,6 +1769,14 @@ type = integer intent = in optional = F +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in + optional = F [dtf] standard_name = time_step_for_dynamics long_name = dynamics timestep @@ -1859,6 +1858,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [nwfa] standard_name = water_friendly_aerosol_number_concentration long_name = number concentration of water-friendly aerosols diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 72efcea60..483eccdf8 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -26,10 +26,10 @@ end subroutine GFS_surface_generic_pre_finalize !! subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, stype, vtype, slope, & prsik_1, prslk_1, tsfc, phil, con_g, & - sigmaf, soiltyp, vegtype, slopetyp, work3, tsurf, zlvl, do_sppt, ca_global,dtdtr,& + sigmaf, soiltyp, vegtype, slopetyp, work3, tsurf, zlvl, & drain_cpl, dsnow_cpl, rain_cpl, snow_cpl, lndp_type, n_var_lndp, sfc_wts, & lndp_var_list, lndp_prt_list, & - z01d, zt1d, bexp1d, xlai1d, vegf1d, lndp_vgf, & + z01d, zt1d, bexp1d, xlai1d, vegf1d, lndp_vgf, sfc_wts_inv, & cplflx, flag_cice, islmsk_cice, slimskin_cpl, tisfc, tsfco, fice, hice, & wind, u1, v1, cnvwind, smcwlt2, smcref2, errmsg, errflg) @@ -51,8 +51,6 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, real(kind=kind_phys), dimension(im), intent(inout) :: sigmaf, work3, tsurf, zlvl ! Stochastic physics / surface perturbations - logical, intent(in) :: do_sppt, ca_global - real(kind=kind_phys), dimension(im,levs), intent(out) :: dtdtr real(kind=kind_phys), dimension(im), intent(out) :: drain_cpl real(kind=kind_phys), dimension(im), intent(out) :: dsnow_cpl real(kind=kind_phys), dimension(im), intent(in) :: rain_cpl @@ -68,6 +66,7 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, real(kind=kind_phys), dimension(im), intent(out) :: xlai1d real(kind=kind_phys), dimension(im), intent(out) :: vegf1d real(kind=kind_phys), intent(out) :: lndp_vgf + real(kind=kind_phys), dimension(im,n_var_lndp), intent(inout) :: sfc_wts_inv logical, intent(in) :: cplflx real(kind=kind_phys), dimension(im), intent(in) :: slimskin_cpl @@ -99,15 +98,14 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, errmsg = '' errflg = 0 - ! Set initial quantities for stochastic physics deltas - if (do_sppt .or. ca_global) then - dtdtr = 0.0 - endif ! Scale random patterns for surface perturbations with perturbation size ! Turn vegetation fraction pattern into percentile pattern lndp_vgf=-999. + if (lndp_type>0) then + sfc_wts_inv(:,:)=sfc_wts(:,:) + endif if (lndp_type==1) then do k =1,n_var_lndp select case(lndp_var_list(k)) diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 68713ab19..d4c8b1bca 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -188,31 +188,6 @@ kind = kind_phys intent = inout optional = F -[do_sppt] - standard_name = flag_for_stochastic_physics_perturbations - long_name = flag for stochastic physics perturbations - units = flag - dimensions = () - type = logical - intent = in - optional = F -[ca_global] - standard_name = flag_for_global_cellular_automata - long_name = switch for global ca - units = flag - dimensions = () - type = logical - intent = in - optional = F -[dtdtr] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_on_physics_time_step - long_name = temp. change due to radiative heating per time step - units = K - dimensions = (horizontal_loop_extent,vertical_dimension) - type = real - kind = kind_phys - intent = out - optional = F [drain_cpl] standard_name = tendency_of_lwe_thickness_of_precipitation_amount_for_coupling long_name = change in rain_cpl (coupling_type) @@ -274,8 +249,17 @@ kind = kind_phys intent = in optional = F +[sfc_wts_inv] + standard_name = weights_for_stochastic_surface_physics_perturbation_flipped + long_name = weights for stochastic surface physics perturbation, flipped + units = none + dimensions = (horizontal_loop_extent,number_of_land_surface_variables_perturbed) + type = real + kind = kind_phys + intent = inout + optional = F [lndp_prt_list] - standard_name =magnitude_of_perturbations_for_landperts + standard_name = magnitude_of_perturbations_for_landperts long_name = magnitude of perturbations for landperts units = variable dimensions = (number_of_land_surface_variables_perturbed) diff --git a/physics/aerclm_def.F b/physics/aerclm_def.F index 84852a1de..426881fe4 100644 --- a/physics/aerclm_def.F +++ b/physics/aerclm_def.F @@ -2,8 +2,8 @@ module aerclm_def use machine , only : kind_phys implicit none - integer, parameter :: levsaer=50, ntrcaerm=15, timeaer=12 - integer :: latsaer, lonsaer, ntrcaer + integer, parameter :: levsaer=72, ntrcaerm=15, timeaer=12 + integer :: latsaer, lonsaer, ntrcaer, levsw character*10 :: specname(ntrcaerm) real (kind=kind_phys):: aer_time(13) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index e7cd6ca20..bed73c5be 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -11,7 +11,7 @@ module aerinterp private - public :: read_aerdata, setindxaer, aerinterpol + public :: read_aerdata, setindxaer, aerinterpol, read_aerdataf contains @@ -32,11 +32,6 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) logical :: file_exist integer, allocatable :: invardims(:) - real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff - real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx - real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp - real(kind=kind_io8),allocatable,dimension(:) :: aer_lati - real(kind=kind_io8),allocatable,dimension(:) :: aer_loni ! !! =================================================================== if (me == master) then @@ -72,50 +67,62 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) ! specify latsaer, lonsaer, hmx lonsaer = dim1 latsaer = dim2 - hmx = int(dim1/2) ! to swap long from W-E to E-W + levsw = dim3 if(me==master) then print *, 'MERRA2 dim: ',dim1, dim2, dim3 endif ! allocate arrays - if (.not. allocated(aer_loni)) then - allocate (aer_loni(lonsaer)) - allocate (aer_lati(latsaer)) - endif if (.not. allocated(aer_lat)) then allocate(aer_lat(latsaer)) allocate(aer_lon(lonsaer)) - allocate(aerin(lonsaer,latsaer,levsaer,ntrcaerm,timeaer)) - allocate(aer_pres(lonsaer,latsaer,levsaer,timeaer)) endif ! construct lat/lon array call nf_inq_varid(ncid, 'lat', varid) - call nf_get_var(ncid, varid, aer_lati) + call nf_get_var(ncid, varid, aer_lat) call nf_inq_varid(ncid, 'lon', varid) - call nf_get_var(ncid, varid, aer_loni) - - do i = 1, hmx ! flip from (-180,180) to (0,360) - if(aer_loni(i)<0.) aer_loni(i)=aer_loni(i)+360. - aer_lon(i+hmx) = aer_loni(i) - aer_lon(i) = aer_loni(i+hmx) - enddo + call nf_get_var(ncid, varid, aer_lon) + call nf_close(ncid) + END SUBROUTINE read_aerdata +! +!********************************************************************** + SUBROUTINE read_aerdataf (iamin, iamax, jamin, jamax, & + me, master, iflip, idate, errmsg, errflg) + use machine, only: kind_phys, kind_io4, kind_io8 + use aerclm_def + use netcdf - do i = 1, latsaer - aer_lat(i) = aer_lati(i) - enddo +!--- in/out + integer, intent(in) :: me, master, iflip, idate(4) + integer, intent(in) :: iamin, iamax, jamin, jamax + character(len=*), intent(inout) :: errmsg + integer, intent(inout) :: errflg - call nf_close(ncid) +!--- locals + integer :: ncid, varid + integer :: i, j, k, n, ii, imon, klev + character :: fname*50, mn*2, vname*10 + logical :: file_exist + integer, allocatable :: invardims(:) + real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff + real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx + real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp +! + if (.not. allocated(aerin)) then + allocate(aerin(iamin:iamax,jamin:jamax,levsaer,ntrcaerm,timeaer)) + allocate(aer_pres(iamin:iamax,jamin:jamax,levsaer,timeaer)) + endif ! allocate local working arrays if (.not. allocated(buff)) then - allocate (buff(lonsaer, latsaer, dim3)) - allocate (pres_tmp(lonsaer,dim3)) + allocate (buff(lonsaer, latsaer, levsw)) + allocate (pres_tmp(lonsaer,levsw)) endif if (.not. allocated(buffx)) then - allocate (buffx(lonsaer, latsaer, dim3,1)) + allocate (buffx(lonsaer, latsaer, levsw,1)) endif !! =================================================================== @@ -137,11 +144,11 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) call nf_inq_varid(ncid, "DELP", varid) call nf_get_var(ncid, varid, buff) - do j = 1, latsaer - do i = 1, lonsaer + do j = jamin, jamax + do i = iamin, iamax ! constract pres_tmp (top-down), note input is top-down pres_tmp(i,1) = 0. - do k=2, dim3 + do k=2, levsw pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) enddo !k-loop enddo !i-loop (lon) @@ -151,11 +158,10 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) if ( iflip == 0 ) then ! data from toa to sfc klev = k else ! data from sfc to top - klev = ( dim3 - k ) + 1 + klev = ( levsw - k ) + 1 endif - do i = 1, hmx - aer_pres(i+hmx,j,k,imon)= 1.d0*pres_tmp(i,klev) - aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i+hmx,klev) + do i = iamin, iamax + aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i,klev) enddo !i-loop (lon) enddo !k-loop (lev) enddo !j-loop (lat) @@ -168,22 +174,18 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) call nf_inq_varid(ncid, vname, varid) call nf_get_var(ncid, varid, buffx) - do j = 1, latsaer + do j = jamin, jamax do k = 1, levsaer ! input is from toa to sfc if ( iflip == 0 ) then ! data from toa to sfc klev = k else ! data from sfc to top - klev = ( dim3 - k ) + 1 + klev = ( levsw - k ) + 1 endif - do i = 1, hmx - aerin(i+hmx,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1) - if(aerin(i+hmx,j,k,ii,imon)<0.or.aerin(i+hmx,j,k,ii,imon)>1.) then - aerin(i+hmx,j,k,ii,imon) = 0. - end if - aerin(i,j,k,ii,imon) = 1.d0*buffx(i+hmx,j,klev,1) + do i = iamin, iamax + aerin(i,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1) if(aerin(i,j,k,ii,imon)<0.or.aerin(i,j,k,ii,imon)>1.) then - aerin(i,j,k,ii,imon) = 0. + aerin(i,j,k,ii,imon) = 1.e-15 end if enddo !i-loop (lon) enddo !k-loop (lev) @@ -195,13 +197,9 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) call nf_close(ncid) enddo !imon-loop !--- - deallocate (aer_loni, aer_lati) deallocate (buff, pres_tmp) deallocate (buffx) - - END SUBROUTINE read_aerdata -! -!********************************************************************** + END SUBROUTINE read_aerdataf ! SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, & iindx1,iindx2,ddx,me,master) @@ -341,7 +339,6 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & +TEMI*DDY(j)*aer_pres(I1,J2,L,n1)+DDX(j)*TEMJ*aer_pres(I2,J1,L,n1))& +tx2*(TEMI*TEMJ*aer_pres(I1,J1,L,n2)+DDX(j)*DDY(J)*aer_pres(I2,J2,L,n2) & +TEMI*DDY(j)*aer_pres(I1,J2,L,n2)+DDX(j)*TEMJ*aer_pres(I2,J1,L,n2)) - ENDDO ENDDO @@ -369,7 +366,7 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & tx1 = temi/(aerpres(j,i1) - aerpres(j,i2)) tx2 = temj/(aerpres(j,i1) - aerpres(j,i2)) DO ii = 1, ntrcaer - aerout(j,L,ii)= aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2 + aerout(j,L,ii)= aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2 ENDDO endif ENDDO !L-loop diff --git a/physics/cires_ugwpv1_module.F90 b/physics/cires_ugwpv1_module.F90 index 4746d61ff..f5fe7f2ec 100644 --- a/physics/cires_ugwpv1_module.F90 +++ b/physics/cires_ugwpv1_module.F90 @@ -143,8 +143,8 @@ module cires_ugwpv1_module !----------------------------------------------------------------------------------- subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & - con_rerth, fn_nml2, lonr, latr, levs, ak, bk, pref, dtp, & - errmsg, errflg) + con_rerth, fn_nml2, input_nml_file, lonr, latr, levs, ak, bk, & + pref, dtp, errmsg, errflg) ! ! input_nml_file ='input.nml'=fn_nml ..... OLD_namelist and cdmvgwd(4) Corrected Bug Oct 4 ! @@ -177,7 +177,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & real(kind=kind_phys), intent (in) :: con_pi, con_rerth character(len=64), intent (in) :: fn_nml2 - character(len=64), parameter :: fn_nml='input.nml' + character(len=*), intent (in) :: input_nml_file(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -187,7 +187,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & integer :: ios logical :: exists - integer :: ncid, iernc, vid, dimid, status + integer :: ncid, iernc, vid, dimid, status integer :: k integer :: ddd_ugwp, curday_ugwp ! integer :: version @@ -196,21 +196,23 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & errmsg = '' errflg = 0 -! - if (me == master) print *, trim (fn_nml), ' GW-namelist file ' - inquire (file =trim (fn_nml) , exist = exists) +#ifdef INTERNAL_FILE_NML + read (input_nml_file, nml = cires_ugwp_nml) +#else + if (me == master) print *, trim (fn_nml2), ' GW-namelist file ' + inquire (file =trim (fn_nml2) , exist = exists) ! if (.not. exists) then - write(errmsg,'(3a)') 'cires_ugwpv1_init: namelist file: ', trim (fn_nml), ' does not exist' + write(errmsg,'(3a)') 'cires_ugwpv1_init: namelist file: ', trim (fn_nml2), ' does not exist' errflg = 1 return else - open (unit = nlunit, file = trim(fn_nml), action = 'read', status = 'old', iostat = ios) + open (unit = nlunit, file = trim(fn_nml2), action = 'read', status = 'old', iostat = ios) endif rewind (nlunit) read (nlunit, nml = cires_ugwp_nml) close (nlunit) -! +#endif strsolver= knob_ugwp_orosolv diff --git a/physics/dcyc2.f b/physics/dcyc2.f index 389496d07..f29d593a3 100644 --- a/physics/dcyc2.f +++ b/physics/dcyc2.f @@ -48,7 +48,7 @@ end subroutine dcyc2t3_finalize ! im, levs, deltim, fhswr, ! ! dry, icy, wet ! ! input/output: ! -! dtdt,dtdtc, ! +! dtdt,dtdtnp, ! ! outputs: ! ! adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw, ! ! adjsfculw_lnd,adjsfculw_ice,adjsfculw_wat,xmu,xcosz, ! @@ -100,8 +100,7 @@ end subroutine dcyc2t3_finalize ! input/output: ! ! dtdt(im,levs)- real, model time step adjusted total radiation ! ! heating rates ( k/s ) ! -! dtdtc(im,levs)- real, model time step adjusted clear sky radiation! -! heating rates ( k/s ) ! +! dtdtnp(im,levs)- real, heating rate adjustment for SPPT ! ! ! ! outputs: ! ! adjsfcdsw(im)- real, time step adjusted sfc dn sw flux (w/m**2) ! @@ -181,9 +180,10 @@ subroutine dcyc2t3_run & & im, levs, deltim, fhswr, & & dry, icy, wet, & & use_LW_jacobian, sfculw, sfculw_jac, & + & pert_radtend, do_sppt,ca_global, & ! & dry, icy, wet, lprnt, ipr, & ! --- input/output: - & dtdt,dtdtc, & + & dtdt,dtdtnp, & ! --- outputs: & adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw, & & adjsfculw_lnd,adjsfculw_ice,adjsfculw_wat,xmu,xcosz, & @@ -212,7 +212,8 @@ subroutine dcyc2t3_run & ! integer, intent(in) :: ipr ! logical lprnt logical, dimension(im), intent(in) :: dry, icy, wet - logical, intent(in) :: use_LW_jacobian + logical, intent(in) :: use_LW_jacobian, pert_radtend + logical, intent(in) :: do_sppt,ca_global real(kind=kind_phys), intent(in) :: solhr, slag, cdec, sdec, & & deltim, fhswr @@ -232,8 +233,8 @@ subroutine dcyc2t3_run & &, swhc, hlwc ! --- input/output: - real(kind=kind_phys), dimension(im,levs), intent(inout) :: dtdt & - &, dtdtc + real(kind=kind_phys), dimension(im,levs), intent(inout) :: dtdt + real(kind=kind_phys), dimension(:,:), intent(inout) :: dtdtnp ! --- outputs: real(kind=kind_phys), dimension(im), intent(out) :: & @@ -358,9 +359,25 @@ subroutine dcyc2t3_run & do k = 1, levs do i = 1, im dtdt(i,k) = dtdt(i,k) + swh(i,k)*xmu(i) + hlw(i,k) - dtdtc(i,k) = dtdtc(i,k) + swhc(i,k)*xmu(i) + hlwc(i,k) enddo enddo + if (do_sppt .or. ca_global) then + if (pert_radtend) then +! clear sky + do k = 1, levs + do i = 1, im + dtdtnp(i,k) = dtdtnp(i,k) + swhc(i,k)*xmu(i) + hlwc(i,k) + enddo + enddo + else +! all sky + do k = 1, levs + do i = 1, im + dtdtnp(i,k) = dtdtnp(i,k) + swh(i,k)*xmu(i) + hlw(i,k) + enddo + enddo + endif + endif ! return !................................... diff --git a/physics/dcyc2.meta b/physics/dcyc2.meta index efba0a5f5..e4a22477d 100644 --- a/physics/dcyc2.meta +++ b/physics/dcyc2.meta @@ -362,9 +362,33 @@ kind = kind_phys intent = inout optional = F -[dtdtc] - standard_name = tendency_of_air_temperature_due_to_radiative_heating_assuming_clear_sky - long_name = clear sky radiative (shortwave + longwave) heating rate at current time +[pert_radtend] + standard_name = flag_for_stochastic_radiative_heating_perturbations + long_name = flag for stochastic radiative heating perturbations + units = flag + dimensions = () + type = logical + intent = in + optional = F +[do_sppt] + standard_name = flag_for_stochastic_physics_perturbations + long_name = flag for stochastic physics perturbations + units = flag + dimensions = () + type = logical + intent = in + optional = F +[ca_global] + standard_name = flag_for_global_cellular_automata + long_name = switch for global ca + units = flag + dimensions = () + type = logical + intent = in + optional = F +[dtdtnp] + standard_name = tendency_of_air_temperature_to_withold_from_sppt + long_name = temp. change from physics that should not be perturbed by sppt units = K s-1 dimensions = (horizontal_loop_extent,vertical_dimension) type = real diff --git a/physics/gcm_shoc.F90 b/physics/gcm_shoc.F90 index dd7791e18..af7d6db49 100644 --- a/physics/gcm_shoc.F90 +++ b/physics/gcm_shoc.F90 @@ -24,7 +24,7 @@ end subroutine shoc_finalize !! \htmlinclude shoc_run.html !! subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, & - con_pi, con_fvirt, dtp, prsl, delp, phii, phil, u, v, omega, rhc, & + con_pi, con_fvirt, con_eps, dtp, prsl, delp, phii, phil, u, v, omega, rhc, & supice, pcrit, cefac, cesfac, tkef1, dis_opt, hflx, evap, prnum, & gt0, gq0, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc, & cld_sgs, tke, tkh, wthv_sec, errmsg, errflg) @@ -32,7 +32,8 @@ subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_ implicit none integer, intent(in) :: nx, nzm, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc - real(kind=kind_phys), intent(in) :: tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, con_fvirt, & + real(kind=kind_phys), intent(in) :: tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, & + con_rd, con_pi, con_fvirt, con_eps, & dtp, supice, pcrit, cefac, cesfac, tkef1, dis_opt ! real(kind=kind_phys), intent(in), dimension(nx) :: hflx, evap @@ -118,7 +119,8 @@ subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_ rhc, supice, pcrit, cefac, cesfac, tkef1, dis_opt, & cld_sgs, tke, hflx, evap, prnum, tkh, wthv_sec, & ntlnc, ncpl, ncpi, & - con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, con_fvirt) + con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, & + con_fvirt, con_eps) if (ntiw < 0) then ! this is valid only for Zhao-Carr scheme do k=1,nzm @@ -166,13 +168,13 @@ subroutine shoc_work (ix, nx, nzm, nz, dtn, & pcrit, cefac, cesfac, tkef1, dis_opt, & cld_sgs, tke, hflx, evap, prnum, tkh, & wthv_sec, ntlnc, ncpl, ncpi, & - cp, ggr, lcond, lfus, rv, rgas, pi, epsv) + cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps) use funcphys , only : fpvsl, fpvsi, fpvs ! saturation vapor pressure for water & ice implicit none - real, intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv + real, intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps integer, intent(in) :: ix ! max number of points in the physics window in the x integer, intent(in) :: nx ! Number of points in the physics window in the x @@ -219,7 +221,7 @@ subroutine shoc_work (ix, nx, nzm, nz, dtn, & real, intent(in) :: prnum (nx,nzm) ! turbulent Prandtl number real, intent(inout) :: wthv_sec (ix,nzm) ! Buoyancy flux, K*m/s - real, parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, eps=0.622_kp, & + real, parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, & three=3.0_kp, oneb3=one/three, twoby3=two/three, fourb3=twoby3+twoby3 real, parameter :: sqrt2 = sqrt(two), twoby15 = two / 15.0_kp, & nmin = 1.0_kp, RI_cub = 6.4e-14_kp, RL_cub = 1.0e-15_kp, & diff --git a/physics/gcm_shoc.meta b/physics/gcm_shoc.meta index 8cb03727d..047286317 100644 --- a/physics/gcm_shoc.meta +++ b/physics/gcm_shoc.meta @@ -113,6 +113,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [dtp] standard_name = time_step_for_physics long_name = time step for physics diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90 index 1ccedb956..c41323ad5 100644 --- a/physics/gfdl_cloud_microphys.F90 +++ b/physics/gfdl_cloud_microphys.F90 @@ -113,7 +113,7 @@ end subroutine gfdl_cloud_microphys_finalize !! \htmlinclude gfdl_cloud_microphys_run.html !! subroutine gfdl_cloud_microphys_run( & - levs, im, con_g, con_fvirt, con_rd, frland, garea, islmsk, & + levs, im, con_g, con_fvirt, con_rd, con_eps, frland, garea, islmsk, & gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, gq0_ntsw, gq0_ntgl, gq0_ntclamt, & gt0, gu0, gv0, vvl, prsl, phii, del, & rain0, ice0, snow0, graupel0, prcp0, sr, & @@ -134,7 +134,7 @@ subroutine gfdl_cloud_microphys_run( & ! interface variables integer, intent(in ) :: levs, im - real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd + real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd, con_eps real(kind=kind_phys), intent(in ), dimension(1:im) :: frland, garea integer, intent(in ), dimension(1:im) :: islmsk real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, & @@ -295,7 +295,7 @@ subroutine gfdl_cloud_microphys_run( & allocate(den(1:im,1:levs)) do k=1,levs do i=1,im - den(i,k)=0.622*prsl(i,k)/(con_rd*gt0(i,k)*(gq0(i,k)+0.622)) + den(i,k)=con_eps*prsl(i,k)/(con_rd*gt0(i,k)*(gq0(i,k)+con_eps)) enddo enddo call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), & diff --git a/physics/gfdl_cloud_microphys.meta b/physics/gfdl_cloud_microphys.meta index 07847ed17..961a3e33f 100644 --- a/physics/gfdl_cloud_microphys.meta +++ b/physics/gfdl_cloud_microphys.meta @@ -168,6 +168,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [frland] standard_name = land_area_fraction_for_microphysics long_name = land area fraction used in microphysics schemes diff --git a/physics/gfdl_sfc_layer.F90 b/physics/gfdl_sfc_layer.F90 index 93e38c982..008e716e2 100644 --- a/physics/gfdl_sfc_layer.F90 +++ b/physics/gfdl_sfc_layer.F90 @@ -377,7 +377,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -526,7 +526,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -633,7 +633,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), znt_ocn(i), mznt(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -756,7 +756,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, & pert_Cd, ens_random_seed, ens_Cdamp, & upc,vpc,tpc,rpc,dt,wind10,xxfh2,ntsflg,sfenth, & - tzot, errmsg, errflg) + tzot, ep2, errmsg, errflg) !------------------------------------------------------------------------ ! @@ -819,6 +819,8 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: tpc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: rpc + real(kind=kind_phys), intent ( in) :: ep2 + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -1207,7 +1209,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m if(psps1 .EQ. 0.0)then psps1 = .1 endif - rstso(i) = 0.622*estso(i)/psps1 + rstso(i) = ep2*estso(i)/psps1 vrts (i) = 1. + boycon*ecof(i)*rstso(i) enddo @@ -1735,7 +1737,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m if(psps2 .EQ. 0.0)then psps2 = .1 endif - rstsop(i) = 0.622*estsop(i)/psps2 + rstsop(i) = ep2*estsop(i)/psps2 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - & diff --git a/physics/gfdl_sfc_layer.meta b/physics/gfdl_sfc_layer.meta index 77024c813..a9829fec3 100644 --- a/physics/gfdl_sfc_layer.meta +++ b/physics/gfdl_sfc_layer.meta @@ -783,6 +783,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [qss_lnd] standard_name = surface_specific_humidity_over_land long_name = surface air saturation specific humidity over land @@ -790,6 +792,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [qss_ice] standard_name = surface_specific_humidity_over_ice long_name = surface air saturation specific humidity over ice @@ -797,6 +801,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/m_micro.F90 b/physics/m_micro.F90 index 77b51ed62..8e6d6698e 100644 --- a/physics/m_micro.F90 +++ b/physics/m_micro.F90 @@ -542,16 +542,19 @@ subroutine m_micro_run( im, lm, flipv, dt_i & & NCPI(I,K), qc_min) if (rnw(i,k) <= qc_min(1)) then ncpr(i,k) = zero + rnw(i,k) = zero elseif (ncpr(i,k) <= nmin) then ! make sure NL > 0 if Q >0 ncpr(i,k) = max(rnw(i,k) / (fourb3 * PI *RL_cub*997.0_kp), nmin) endif if (snw(i,k) <= qc_min(2)) then ncps(i,k) = zero + snw(i,k) = zero elseif (ncps(i,k) <= nmin) then ncps(i,k) = max(snw(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif if (qgl(i,k) <= qc_min(2)) then ncgl(i,k) = zero + qgl(i,k) = zero elseif (ncgl(i,k) <= nmin) then ncgl(i,k) = max(qgl(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif @@ -1696,16 +1699,19 @@ subroutine m_micro_run( im, lm, flipv, dt_i & QI_TOT(I,K) = QILS(I,K) + QICN(I,K) if (rnw(i,k) <= qc_min(1)) then ncpr(i,k) = zero + rnw(i,k) = zero elseif (ncpr(i,k) <= nmin) then ! make sure NL > 0 if Q >0 ncpr(i,k) = max(rnw(i,k) / (fourb3 * PI *RL_cub*997.0_kp), nmin) endif if (snw(i,k) <= qc_min(2)) then ncps(i,k) = zero + snw(i,k) = zero elseif (ncps(i,k) <= nmin) then ncps(i,k) = max(snw(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif if (qgl(i,k) <= qc_min(2)) then ncgl(i,k) = zero + qgl(i,k) = zero elseif (ncgl(i,k) <= nmin) then ncgl(i,k) = max(qgl(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif @@ -1736,16 +1742,19 @@ subroutine m_micro_run( im, lm, flipv, dt_i & ! if (rnw(i,k) <= qc_min(1)) then ncpr(i,k) = zero + rnw(i,k) = zero elseif (ncpr(i,k) <= nmin) then ! make sure NL > 0 if Q >0 ncpr(i,k) = max(rnw(i,k) / (fourb3 * PI *RL_cub*997.0_kp), nmin) endif if (snw(i,k) <= qc_min(2)) then ncps(i,k) = zero + snw(i,k) = zero elseif (ncps(i,k) <= nmin) then ncps(i,k) = max(snw(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif if (qgl(i,k) <= qc_min(2)) then ncgl(i,k) = zero + qgl(i,k) = zero elseif (ncgl(i,k) <= nmin) then ncgl(i,k) = max(qgl(i,k) / (fourb3 * PI *RL_cub*500.0_kp), nmin) endif diff --git a/physics/micro_mg2_0.F90 b/physics/micro_mg2_0.F90 index 744b46ebc..73c46392d 100644 --- a/physics/micro_mg2_0.F90 +++ b/physics/micro_mg2_0.F90 @@ -1792,7 +1792,7 @@ subroutine micro_mg_tend ( & nnucct(i,k) = ratio * nnucct(i,k) npsacws(i,k) = ratio * npsacws(i,k) nsubc(i,k) = ratio * nsubc(i,k) - end if + endif mnuccri(i,k) = zero nnuccri(i,k) = zero @@ -1800,15 +1800,17 @@ subroutine micro_mg_tend ( & if (do_cldice) then ! freezing of rain to produce ice if mean rain size is smaller than Dcs - if (lamr(i,k) > qsmall .and. one/lamr(i,k) < Dcs) then - mnuccri(i,k) = mnuccr(i,k) - nnuccri(i,k) = nnuccr(i,k) - mnuccr(i,k) = zero - nnuccr(i,k) = zero - end if - end if + if (lamr(i,k) > qsmall) then + if (one/lamr(i,k) < Dcs) then + mnuccri(i,k) = mnuccr(i,k) + nnuccri(i,k) = nnuccr(i,k) + mnuccr(i,k) = zero + nnuccr(i,k) = zero + endif + endif + endif - end do + enddo do i=1,mgncol diff --git a/physics/micro_mg3_0.F90 b/physics/micro_mg3_0.F90 index 636293b86..dde143c4d 100644 --- a/physics/micro_mg3_0.F90 +++ b/physics/micro_mg3_0.F90 @@ -2448,7 +2448,7 @@ subroutine micro_mg_tend ( & nnucct(i,k) = ratio * nnucct(i,k) npsacws(i,k) = ratio * npsacws(i,k) nsubc(i,k) = ratio * nsubc(i,k) - end if + endif mnuccri(i,k) = zero nnuccri(i,k) = zero @@ -2456,15 +2456,17 @@ subroutine micro_mg_tend ( & if (do_cldice) then ! freezing of rain to produce ice if mean rain size is smaller than Dcs - if (lamr(i,k) > qsmall .and. one/lamr(i,k) < Dcs) then - mnuccri(i,k) = mnuccr(i,k) - nnuccri(i,k) = nnuccr(i,k) - mnuccr(i,k) = zero - nnuccr(i,k) = zero - end if - end if + if (lamr(i,k) > qsmall) then + if (one/lamr(i,k) < Dcs) then + mnuccri(i,k) = mnuccr(i,k) + nnuccri(i,k) = nnuccr(i,k) + mnuccr(i,k) = zero + nnuccr(i,k) = zero + endif + endif + endif - end do + enddo do i=1,mgncol diff --git a/physics/module_gfdl_cloud_microphys.F90 b/physics/module_gfdl_cloud_microphys.F90 index 5750d27fd..6cc9275dc 100644 --- a/physics/module_gfdl_cloud_microphys.F90 +++ b/physics/module_gfdl_cloud_microphys.F90 @@ -47,6 +47,9 @@ module gfdl_cloud_microphys_mod public gfdl_cloud_microphys_mod_driver, gfdl_cloud_microphys_mod_init, & gfdl_cloud_microphys_mod_end, cloud_diagnosis +! public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist +! public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d +! public setup_con, wet_bulb real :: missing_value = - 1.e10 @@ -116,7 +119,18 @@ module gfdl_cloud_microphys_mod real, parameter :: sfcrho = 1.2 !< surface air density real, parameter :: rhor = 1.e3 !< density of rain water, lin83 + ! intercept parameters + + real, parameter :: rnzr = 8.0e6 ! lin83 + real, parameter :: rnzs = 3.0e6 ! lin83 + real, parameter :: rnzg = 4.0e6 ! rh84 + real, parameter :: rnzh = 4.0e4 ! lin83 --- lmh 29 sep 17 + + ! density parameters + real, parameter :: rhoh = 0.917e3 ! lin83 --- lmh 29 sep 17 + + public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw !< constants for accretions real :: acco (3, 4) !< constants for accretions real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5) @@ -283,6 +297,7 @@ module gfdl_cloud_microphys_mod logical :: use_ppm = .false. !< use ppm fall scheme logical :: mono_prof = .true. !< perform terminal fall with mono ppm scheme logical :: mp_print = .false. !< cloud microphysics debugging printout + logical :: do_hail = .false. !< use hail parameters instead of graupel ! real :: global_area = - 1. @@ -316,7 +331,7 @@ module gfdl_cloud_microphys_mod rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, & do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, & mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, & - resmin, resmax, regmin, regmax, tintqs + resmin, resmax, regmin, regmax, tintqs, do_hail public & mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, & @@ -330,7 +345,7 @@ module gfdl_cloud_microphys_mod rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, & do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, & mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, & - resmin, resmax, regmin, regmax, tintqs + resmin, resmax, regmin, regmax, tintqs, do_hail contains @@ -1796,9 +1811,11 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & !! threshold from wsm6 scheme, Hong et al. (2004) \cite hong_et_al_2004, !! eq (13) : qi0_crt ~0.8e-4. ! ----------------------------------------------------------------------- - - qim = qi0_crt / den (k) - + if (qi0_crt < 0.) then + qim = - qi0_crt + else + qim = qi0_crt / den (k) + endif ! ----------------------------------------------------------------------- ! assuming linear subgrid vertical distribution of cloud ice ! the mismatch computation following lin et al. 1994, mwr @@ -3280,8 +3297,10 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) real, parameter :: vcons = 6.6280504 real, parameter :: vcong = 87.2382675 + real, parameter :: vconh = vcong * sqrt (rhoh / rhog) real, parameter :: norms = 942477796.076938 real, parameter :: normg = 5026548245.74367 + real, parameter :: normh = pi * rhoh * rnzh real, dimension (ktop:kbot) :: qden, tc, rhof @@ -3346,10 +3365,19 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) ! ----------------------------------------------------------------------- !> - graupel: ! ----------------------------------------------------------------------- - if (const_vg) then vtg (:) = vg_fac ! 2. else + if (do_hail) then + do k = ktop, kbot + if (qg (k) < thg) then + vtg (k) = vf_min + else + vtg (k) = vg_fac * vconh * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normh))) + vtg (k) = min (vg_max, max (vf_min, vtg (k))) + endif + enddo + else do k = ktop, kbot if (qg (k) < thg) then vtg (k) = vf_min @@ -3359,6 +3387,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) endif enddo endif + endif end subroutine fall_speed @@ -3382,9 +3411,9 @@ subroutine setupm ! intercept parameters - real, parameter :: rnzr = 8.0e6 ! lin83 - real, parameter :: rnzs = 3.0e6 ! lin83 - real, parameter :: rnzg = 4.0e6 ! rh84 +! real, parameter :: rnzr = 8.0e6 ! lin83 +! real, parameter :: rnzs = 3.0e6 ! lin83 +! real, parameter :: rnzg = 4.0e6 ! rh84 ! density parameters @@ -3427,8 +3456,13 @@ subroutine setupm cracs = pisq * rnzr * rnzs * rhos csacr = pisq * rnzr * rnzs * rhor - cgacr = pisq * rnzr * rnzg * rhor - cgacs = pisq * rnzg * rnzs * rhos + if (do_hail) then + cgacr = pisq * rnzr * rnzh * rhor + cgacs = pisq * rnzh * rnzs * rhos + else + cgacr = pisq * rnzr * rnzg * rhor + cgacs = pisq * rnzg * rnzs * rhos + endif cgacs = cgacs * c_pgacs ! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ; @@ -3436,7 +3470,11 @@ subroutine setupm act (1) = pie * rnzs * rhos act (2) = pie * rnzr * rhor - act (6) = pie * rnzg * rhog + if (do_hail) then + act (6) = pie * rnzh * rhoh + else + act (6) = pie * rnzg * rhog + endif act (3) = act (2) act (4) = act (1) act (5) = act (2) @@ -3457,7 +3495,11 @@ subroutine setupm craci = pie * rnzr * alin * gam380 / (4. * act (2) ** 0.95) csaci = csacw * c_psaci - cgacw = pie * rnzg * gam350 * gcon / (4. * act (6) ** 0.875) + if (do_hail) then + cgacw = pie * rnzh * gam350 * gcon / (4. * act (6) ** 0.875) + else + cgacw = pie * rnzg * gam350 * gcon / (4. * act (6) ** 0.875) + endif ! cgaci = cgacw * 0.1 ! sjl, may 28, 2012 @@ -3470,7 +3512,11 @@ subroutine setupm ! subl and revp: five constants for three separate processes cssub (1) = 2. * pie * vdifu * tcond * rvgas * rnzs - cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzg + if (do_hail) then + cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzh + else + cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzg + endif crevp (1) = 2. * pie * vdifu * tcond * rvgas * rnzr cssub (2) = 0.78 / sqrt (act (1)) cgsub (2) = 0.78 / sqrt (act (6)) @@ -3498,8 +3544,13 @@ subroutine setupm ! gmlt: five constants - cgmlt (1) = 2. * pie * tcond * rnzg / hltf - cgmlt (2) = 2. * pie * vdifu * rnzg * hltc / hltf + if (do_hail) then + cgmlt (1) = 2. * pie * tcond * rnzh / hltf + cgmlt (2) = 2. * pie * vdifu * rnzh * hltc / hltf + else + cgmlt (1) = 2. * pie * tcond * rnzg / hltf + cgmlt (2) = 2. * pie * vdifu * rnzg * hltc / hltf + endif cgmlt (3) = cgsub (2) cgmlt (4) = cgsub (3) cgmlt (5) = ch2o / hltf diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 5d5e631f5..6031af94a 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -28,9 +28,9 @@ module mp_thompson !! \section arg_table_mp_thompson_init Argument Table !! \htmlinclude mp_thompson_init.html !! - subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & - imp_physics, imp_physics_thompson, & - convert_dry_rho, & + subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & + restart, imp_physics, & + imp_physics_thompson, convert_dry_rho,& spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa2d, nifa2d, & nwfa, nifa, tgrs, prsl, phil, area, & @@ -43,7 +43,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Interface variables integer, intent(in ) :: ncol integer, intent(in ) :: nlev - real(kind_phys), intent(in ) :: con_g, con_rd + real(kind_phys), intent(in ) :: con_g, con_rd, con_eps logical, intent(in ) :: restart integer, intent(in ) :: imp_physics integer, intent(in ) :: imp_physics_thompson @@ -160,7 +160,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & end if ! Density of moist air in kg m-3 and inverse density of air - rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) orho = 1.0/rho ! Ensure we have 1st guess ice number where mass non-zero but no number. @@ -324,7 +324,7 @@ end subroutine mp_thompson_init !>\section gen_thompson_hrrr Thompson MP General Algorithm !>@{ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & - convert_dry_rho, & + con_eps, convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa, nifa, & nwfa2d, nifa2d, & @@ -344,6 +344,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, intent(in ) :: nlev real(kind_phys), intent(in ) :: con_g real(kind_phys), intent(in ) :: con_rd + real(kind_phys), intent(in ) :: con_eps ! Hydrometeors logical, intent(in ) :: convert_dry_rho real(kind_phys), intent(inout) :: spechum(1:ncol,1:nlev) @@ -474,7 +475,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if !> - Density of air in kg m-3 - rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 w = -omega/(rho*con_g) diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index ed54f8d02..7f1e9197e 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -41,6 +41,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [restart] standard_name = flag_for_restart long_name = flag for restart (warmstart) or coldstart @@ -349,6 +358,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [convert_dry_rho] standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air long_name = flag for converting hydrometeors from moist to dry air diff --git a/physics/radiation_aerosols.f b/physics/radiation_aerosols.f index f732c37ef..e1e66b0d9 100644 --- a/physics/radiation_aerosols.f +++ b/physics/radiation_aerosols.f @@ -561,7 +561,7 @@ subroutine aer_init & laswflg= (mod(iaerflg,10) > 0) ! control flag for sw tropospheric aerosol lalwflg= (mod(iaerflg/10,10) > 0) ! control flag for lw tropospheric aerosol - lavoflg= (iaerflg >= 100) ! control flag for stratospheric volcanic aeros + lavoflg= (mod(iaerflg/100,10) >0) ! control flag for stratospheric volcanic aeros !> -# Call wrt_aerlog() to write aerosol parameter configuration to output logs. @@ -4446,6 +4446,8 @@ subroutine aeropt asy1 = f_zero sca1 = f_zero ssa1 = f_zero + asy = f_zero + ssa = f_zero do m = 1, kcm1 cm = max(aerms(k,m),0.0) * dz1(k) ext1 = ext1 + cm*extrhi_grt(m,ib) diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index 8ebbb3ab1..77fd61fcc 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -2946,8 +2946,13 @@ subroutine spcvrtc & else ! for non-conservative scattering za1 = zgam1*zgam4 + zgam2*zgam3 za2 = zgam1*zgam3 + zgam2*zgam4 - zrk = sqrt ( (zgam1 - zgam2) * (zgam1 + zgam2) ) - zrk2= 2.0 * zrk + zrk = (zgam1 - zgam2) * (zgam1 + zgam2) + if (zrk > eps1) then + zrk = sqrt(zrk) + else + zrk = f_zero + endif + zrk2= zrk + zrk zrp = zrk * cosz zrp1 = f_one + zrp @@ -2993,7 +2998,8 @@ subroutine spcvrtc & ze1r45 = zr4*zexp1 + zr5*zexm1 ! ... collimated beam - if (ze1r45>=-eps1 .and. ze1r45<=eps1) then +! if (ze1r45>=-eps1 .and. ze1r45<=eps1) then + if (abs(ze1r45) <= eps1) then zrefb(kp) = eps1 ztrab(kp) = zexm2 else @@ -3005,7 +3011,11 @@ subroutine spcvrtc & endif ! ... diffuse beam - zden1 = zr4 / (ze1r45 * zrkg1) + if (ze1r45 >= f_zero) then + zden1 = zr4 / max(eps1, ze1r45*zrkg1) + else + zden1 = zr4 / min(-eps1, ze1r45*zrkg1) + endif zrefd(kp) = max(f_zero, min(f_one, & & zgam2*(zexp1 - zexm1)*zden1 )) ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 )) @@ -3171,8 +3181,13 @@ subroutine spcvrtc & else ! for non-conservative scattering za1 = zgam1*zgam4 + zgam2*zgam3 za2 = zgam1*zgam3 + zgam2*zgam4 - zrk = sqrt ( (zgam1 - zgam2) * (zgam1 + zgam2) ) - zrk2= 2.0 * zrk + zrk = (zgam1 - zgam2) * (zgam1 + zgam2) + if (zrk > eps1) then + zrk = sqrt(zrk) + else + zrk = f_zero + endif + zrk2= zrk + zrk zrp = zrk * cosz zrp1 = f_one + zrp @@ -3218,7 +3233,8 @@ subroutine spcvrtc & ze1r45 = zr4*zexp1 + zr5*zexm1 ! ... collimated beam - if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then +! if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then + if ( abs(ze1r45) <= eps1 ) then zrefb(kp) = eps1 ztrab(kp) = zexm2 else @@ -3230,7 +3246,11 @@ subroutine spcvrtc & endif ! ... diffuse beam - zden1 = zr4 / (ze1r45 * zrkg1) + if (ze1r45 >= f_zero) then + zden1 = zr4 / max(eps1, ze1r45*zrkg1) + else + zden1 = zr4 / min(-eps1, ze1r45*zrkg1) + endif zrefd(kp) = max(f_zero, min(f_one, & & zgam2*(zexp1 - zexm1)*zden1 )) ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 )) @@ -3723,8 +3743,13 @@ subroutine spcvrtm & else ! for non-conservative scattering za1 = zgam1*zgam4 + zgam2*zgam3 za2 = zgam1*zgam3 + zgam2*zgam4 - zrk = sqrt ( (zgam1 - zgam2) * (zgam1 + zgam2) ) - zrk2= 2.0 * zrk + zrk = (zgam1 - zgam2) * (zgam1 + zgam2) + if (zrk > eps1) then + zrk = sqrt(zrk) + else + zrk = f_zero + endif + zrk2= zrk + zrk zrp = zrk * cosz zrp1 = f_one + zrp @@ -3770,7 +3795,8 @@ subroutine spcvrtm & ze1r45 = zr4*zexp1 + zr5*zexm1 ! ... collimated beam - if (ze1r45>=-eps1 .and. ze1r45<=eps1) then +! if (ze1r45>=-eps1 .and. ze1r45<=eps1) then + if (abs(ze1r45) <= eps1) then zrefb(kp) = eps1 ztrab(kp) = zexm2 else @@ -3782,7 +3808,11 @@ subroutine spcvrtm & endif ! ... diffuse beam - zden1 = zr4 / (ze1r45 * zrkg1) + if (ze1r45 >= f_zero) then + zden1 = zr4 / max(eps1, ze1r45*zrkg1) + else + zden1 = zr4 / min(-eps1, ze1r45*zrkg1) + endif zrefd(kp) = max(f_zero, min(f_one, & & zgam2*(zexp1 - zexm1)*zden1 )) ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 )) @@ -3935,8 +3965,13 @@ subroutine spcvrtm & else ! for non-conservative scattering za1 = zgam1*zgam4 + zgam2*zgam3 za2 = zgam1*zgam3 + zgam2*zgam4 - zrk = sqrt ( (zgam1 - zgam2) * (zgam1 + zgam2) ) - zrk2= 2.0 * zrk + zrk = (zgam1 - zgam2) * (zgam1 + zgam2) + if (zrk > eps1) then + zrk = sqrt(zrk) + else + zrk = f_zero + endif + zrk2= zrk + zrk zrp = zrk * cosz zrp1 = f_one + zrp @@ -3982,7 +4017,8 @@ subroutine spcvrtm & ze1r45 = zr4*zexp1 + zr5*zexm1 ! ... collimated beam - if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then +! if ( ze1r45>=-eps1 .and. ze1r45<=eps1 ) then + if ( abs(ze1r45) <= eps1 ) then zrefb(kp) = eps1 ztrab(kp) = zexm2 else @@ -3994,7 +4030,11 @@ subroutine spcvrtm & endif ! ... diffuse beam - zden1 = zr4 / (ze1r45 * zrkg1) + if (ze1r45 >= f_zero) then + zden1 = zr4 / max(eps1, ze1r45*zrkg1) + else + zden1 = zr4 / min(-eps1, ze1r45*zrkg1) + endif zrefd(kp) = max(f_zero, min(f_one, & & zgam2*(zexp1 - zexm1)*zden1 )) ztrad(kp) = max(f_zero, min(f_one, zrk2*zden1 )) diff --git a/physics/rascnv.F90 b/physics/rascnv.F90 index 1c311e4cf..e78570f34 100644 --- a/physics/rascnv.F90 +++ b/physics/rascnv.F90 @@ -30,21 +30,25 @@ module rascnv &, FOUR_P2=4.0e2_kp, ONE_M10=1.0e-10_kp& &, ONE_M6=1.0e-6_kp, ONE_M5=1.0e-5_kp & &, ONE_M2=1.0e-2_kp, ONE_M1=1.0e-1_kp & - &, oneolog10=one/log(10.0_kp) & - &, facmb = 0.01_kp & ! conversion factor from Pa to hPa (or mb) - &, cmb2pa = 100.0_kp ! Conversion from hPa to Pa -! - real(kind=kind_phys), parameter :: frac=0.5_kp, crtmsf=0.0_kp & - &, rhfacs=0.75_kp, rhfacl=0.75_kp & - &, face=5.0_kp, delx=10000.0_kp& - &, ddfac=face*delx*0.001_kp & - &, max_neg_bouy=0.15_kp & -! &, max_neg_bouy=pt25_kp & - &, testmb=0.1_kp, testmbi=one/testmb & - &, dpd=0.5_kp, rknob=1.0_kp, eknob=1.0_kp + &, oneolog10=one/log(10.0_kp) & + &, rain_min=1.0e-13_kp & + &, facmb=0.01_kp & ! conversion factor from Pa to hPa (or mb) + &, cmb2pa=100.0_kp ! Conversion from hPa to Pa +! +! real (kind=kind_phys), parameter :: frac=0.5_kp, crtmsf=0.0_kp & + real (kind=kind_phys), parameter :: frac=0.1_kp, crtmsf=0.0_kp & + &, tfrac_max=0.15_kp & + &, rhfacs=0.75_kp, rhfacl=0.75_kp & + &, face=5.0_kp, delx=10000.0_kp & + &, ddfac=face*delx*0.001_kp & + &, max_neg_bouy=0.15_kp & +! &, max_neg_bouy=pt25_kp & + &, testmb=0.1_kp, testmbi=one/testmb & + &, dpd=0.5_kp, rknob=1.0_kp, eknob=1.0_kp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - logical, parameter :: do_aw=.true., cumfrc=.true. & +! logical, parameter :: aw_scal=.false., cumfrc=.true. & + logical, parameter :: aw_scal=.true., cumfrc=.true. & &, updret=.false., vsmooth=.false. & &, wrkfun=.false., crtfun=.true. & &, calkbl=.true., botop=.true., revap=.true. & @@ -67,24 +71,24 @@ module rascnv ! ! For Tilting Angle Specification ! - real(kind=kind_phys) REFP(6), REFR(6), TLAC(8), PLAC(8), TLBPL(7) & - &, drdp(5) + real(kind=kind_phys) :: REFP(6), REFR(6), TLAC(8), PLAC(8), & + TLBPL(7), drdp(5) ! DATA PLAC/100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0/ DATA TLAC/ 35.0, 25.0, 20.0, 17.5, 15.0, 12.5, 10.0, 7.5/ DATA REFP/500.0, 300.0, 250.0, 200.0, 150.0, 100.0/ DATA REFR/ 1.0, 2.0, 3.0, 4.0, 6.0, 8.0/ ! - real(kind=kind_phys) AC(16), AD(16) + real(kind=kind_phys) :: AC(16), AD(16) ! integer, parameter :: nqrp=500001 - real(kind=kind_phys) C1XQRP, C2XQRP, TBQRP(NQRP), TBQRA(NQRP) & - &, TBQRB(NQRP) + real(kind=kind_phys) :: C1XQRP, C2XQRP, TBQRP(NQRP), & + TBQRA(NQRP), TBQRB(NQRP) ! integer, parameter :: nvtp=10001 - real(kind=kind_phys) C1XVTP, C2XVTP, TBVTP(NVTP) + real(kind=kind_phys) :: C1XVTP, C2XVTP, TBVTP(NVTP) ! - real(kind=kind_phys) afc, facdt, & + real(kind=kind_phys) :: afc, facdt, & grav, cp, alhl, alhf, rgas, rkap, nu, pi, & t0c, rv, cvap, cliq, csol, ttp, eps, epsm1,& ! @@ -93,7 +97,6 @@ module rascnv deg2rad, PIINV, testmboalhl, & rvi, facw, faci, hsub, tmix, DEN - contains ! ----------------------------------------------------------------------- @@ -118,12 +121,13 @@ subroutine rascnv_init(me, dt, con_g, con_cp, con_rd, & con_g, con_cp, con_rd, con_rv, con_hvap, & con_hfus, con_fvirt, con_t0c, con_cvap, con_cliq, & con_csol, con_ttp, con_eps, con_epsm1 + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! real(kind=kind_phys), parameter :: actp=1.7_kp, facm=1.00_kp ! - real(kind=kind_phys) PH(15), A(15) + real(kind=kind_phys) :: PH(15), A(15) ! DATA PH/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0 & &, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/ @@ -134,8 +138,6 @@ subroutine rascnv_init(me, dt, con_g, con_cp, con_rd, & ! real(kind=kind_phys) tem, actop, tem1, tem2 integer i, l - logical first - data first/.true./ ! ! Initialize CCPP error handling variables errmsg = '' @@ -169,6 +171,12 @@ subroutine rascnv_init(me, dt, con_g, con_cp, con_rd, & ! VTP = 36.34*SQRT(1.2)* (0.001)**0.1364 ! AFC = -(1.01097e-4_kp*DT)*(3600.0_kp/DT)**0.57777778_kp +! + if (fix_ncld_hr) then + facdt = delt_c / dt + else + facdt = one / 3600.0_kp + endif ! grav = con_g ; cp = con_cp ; alhl = con_hvap alhf = con_hfus ; rgas = con_rd @@ -186,9 +194,9 @@ subroutine rascnv_init(me, dt, con_g, con_cp, con_rd, & picon = half*pi*onebg ; zfac = 0.28888889e-4_kp * ONEBG testmboalhl = testmb/alhl ! - rvi = one/rv ; facw=CVAP-CLIQ - faci = CVAP-CSOL ; hsub=alhl+alhf - tmix = TTP-20.0_kp ; DEN=one/(TTP-TMIX) + rvi = one / rv ; facw = CVAP - CLIQ + faci = CVAP - CSOL ; hsub = alhl + alhf + tmix = TTP - 20.0_kp ; DEN = one / (TTP-TMIX) ! if (me == 0) write(0,*) ' NO DOWNDRAFT FOR CLOUD TYPES' & @@ -286,7 +294,7 @@ end subroutine rascnv_finalize !! \section arg_table_rascnv_run Argument Table !! \htmlinclude rascnv_run.html !! - subroutine rascnv_run(IM, k, ntr, dt, dtf & + subroutine rascnv_run(IM, k, ntr, dt, dtf & &, ccwf, area, dxmin, dxinv & &, psauras, prauras, wminras, dlqf, flipv & &, me, rannum, nrcm, mp_phys, mp_phys_mg & @@ -329,10 +337,12 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & &, psauras(2), prauras(2) & &, wminras(2), dlqf(2) ! - real(kind=kind_phys), dimension(im,k) :: tin, qin, uin, vin & - &, prsl, prslk, phil real(kind=kind_phys), dimension(im,k+1) :: prsi, prsik, phii - real(kind=kind_phys), dimension(im,k) :: ud_mf, dd_mf, dt_mf & + + real(kind=kind_phys), dimension(im,k) :: tin, qin, uin, vin & + &, prsl, prslk, phil & + + &, ud_mf, dd_mf, dt_mf & &, rhc, qlcn, qicn, w_upi & &, cnv_mfd & &, cnv_dqldt, clcn & @@ -344,7 +354,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & real(kind=kind_phys) ccin(im,k,ntr+2) real(kind=kind_phys) trcmin(ntr+2) - real(kind=kind_phys) DT, dtf, qw0, qi0 + real(kind=kind_phys) DT, dtf ! ! Added for aerosol scavenging for GOCART ! @@ -377,17 +387,17 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ! integer :: nrcmax ! Maximum # of random clouds per 1200s ! Integer KCR, KFX, NCMX, NC, KTEM, I, ii, Lm1, l & - &, ntrc, ia, ll, km1, kp1, ipt, lv, KBL, n & + &, ntrc, ll, km1, kp1, ipt, lv, KBL, n & &, KRMIN, KRMAX, KFMAX, kblmx, irnd,ib & &, kblmn, ksfc, ncrnd - real(kind=kind_phys) sgcs(k,im) + real(kind=kind_phys) sgcs(k) ! ! Scavenging related parameters ! real fscav_(ntr+2) ! Fraction scavenged per km ! - fscav_ = zero ! By default no scavenging - if (ntr > 0) then + fscav_ = -999.0_kp ! By default no scavenging + if (ntr > 0 .and. fscav(1) > zero) then do i=1,ntr fscav_(i) = fscav(i) enddo @@ -425,7 +435,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & endif ! !!!!! initialization for microphysics ACheng - if(mp_phys == 10) then + if(mp_phys == mp_phys_mg) then do l=1,K do i=1,im QLCN(i,l) = zero @@ -466,7 +476,6 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & c0i = (psauras(1)*tem1 + psauras(2)*tem2) * tem c0 = (prauras(1)*tem1 + prauras(2)*tem2) * tem if (ccwfac == zero) ccwfac = half - ! ! ctei = .false. ! if (ctei_r(ipt) > ctei_rm) ctei = .true. @@ -482,11 +491,12 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & KFMAX = KRMAX kblmx = 1 kblmn = 1 + sgcs(k) = one DO L=1,KM1 ll = l if (flipv) ll = kp1 -l ! Input variables are bottom to top! SGC = prsl(ipt,ll) * tem - sgcs(l,ipt) = sgc + sgcs(l) = sgc IF (SGC <= 0.050_kp) KRMIN = L ! IF (SGC <= 0.700_kp) KRMAX = L ! IF (SGC <= 0.800_kp) KRMAX = L @@ -510,10 +520,8 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/360) + 0.50001 ! & + 0.50001 ! NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * min(1.0,DTF/360) + 0.1 - facdt = delt_c / dt else NCRND = min(nrcmax, (KRMAX-KRMIN+1)) - facdt = one / 3600.0_kp endif NCRND = min(nrcm,max(NCRND, 1)) ! @@ -779,7 +787,6 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ! IB = IC(NC) ! cloud top level index if (ib > kbl-1) cycle - ! !**************************************************************************** ! if (advtvd) then ! TVD flux limiter scheme for updraft @@ -858,12 +865,12 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & rainp = rain CALL CLOUD(K, KP1, IB, ntrc, kblmx, kblmn & - &, FRAC, MAX_NEG_BOUY, vsmooth, do_aw & + &, FRAC, MAX_NEG_BOUY, vsmooth, aw_scal & &, REVAP, WRKFUN, CALKBL, CRTFUN & &, DT, KDT, TLA, DPD & &, ALFINT, rhfacl, rhfacs, area(ipt) & &, ccwfac, CDRAG(ipt), trcfac & - &, alfind, rhc_l, phi_l, phi_h, PRS, PRSM,sgcs(1,ipt) & + &, alfind, rhc_l, phi_l, phi_h, PRS, PRSM,sgcs & &, TOI, QOI, UVI, QLI, QII, KBL, DDVEL(ipt) & &, TCU, QCU, RCU, PCU, FLX, FLXD, RAIN, WFNC, fscav_ & &, trcmin, ntk-2, c0, wminras(1), c0i, wminras(2) & @@ -880,15 +887,15 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ll = kp1 - ib dt_mf(ipt,ll) = dt_mf(ipt,ll) + flx(ib) - if (mp_phys == 10) then ! Anning Cheng for microphysics 11/14/2015 + if (mp_phys == mp_phys_mg) then ! Anning Cheng for microphysics 11/14/2015 CNV_MFD(ipt,ll) = CNV_MFD(ipt,ll) + flx(ib)/dt -! CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll) -! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt +!! CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll) +!! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll) + flx(ib)* & & max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt -! & max(0.,(QLI(ib)+QII(ib)))/dt/3. +!! & max(0.,(QLI(ib)+QII(ib)))/dt/3. if(flx(ib)<0) write(*,*)"AAA666", flx(ib),QLI(ib),QII(ib) & & ,ipt,ll endif @@ -901,16 +908,17 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & enddo dt_mf(ipt,ib) = dt_mf(ipt,ib) + flx(ib) - if (mp_phys == 10) then ! Anning Cheng for microphysics 11/14/2015 + if (mp_phys == mp_phys_mg) then ! Anning Cheng for microphysics 11/14/2015 CNV_MFD(ipt,ib) = CNV_MFD(ipt,ib) + flx(ib)/dt -! CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib) -! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt +!! CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib) +!! & + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib) + flx(ib)* & & max(zero,(QLI(ib)+QII(ib)-qiid-qlid))/dt -! & max(0.,(QLI(ib)+QII(ib)))/dt/3. +!! & max(0.,(QLI(ib)+QII(ib)))/dt/3. if(flx(ib)<0) write(*,*)"AAA666", flx(ib),QLI(ib),QII(ib) & & ,ipt,ib endif + endif ! ! @@ -931,6 +939,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ENDDO ! End of the NC loop! ! RAINC(ipt) = rain * 0.001_kp ! Output rain is in meters + if (rainc(ipt) < rain_min) rainc(ipt) = zero ktop(ipt) = kp1 kbot(ipt) = 0 @@ -944,9 +953,9 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & ! clw(i) = max(clw(i), zero) ! cli(i) = max(cli(i), zero) - if (sgcs(l,ipt) < 0.93_kp .and. abs(tcu(l)) > one_m10) then -! if (sgcs(l,ipt) < 0.90_kp .and. tcu(l) .ne. zero) then -! if (sgcs(l,ipt) < 0.85_kp .and. tcu(l) .ne. zero) then + if (sgcs(l) < 0.93_kp .and. abs(tcu(l)) > one_m10) then +! if (sgcs(l) < 0.90_kp .and. tcu(l) .ne. zero) then +! if (sgcs(l) < 0.85_kp .and. tcu(l) .ne. zero) then kcnv(ipt) = 1 endif ! New test for convective clouds ! added in 08/21/96 @@ -967,7 +976,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & vin(ipt,ll) = uvi(l,ntr+2) ! V momentum !! for 2M microphysics, always output these variables - if (mp_phys == 10) then + if (mp_phys == mp_phys_mg) then if (advcld) then QLCN(ipt,ll) = max(qli(l)-ccin(ipt,ll,2), zero) QICN(ipt,ll) = max(qii(l)-ccin(ipt,ll,1), zero) @@ -1018,7 +1027,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & vin(ipt,l) = uvi(l,ntr+2) ! V momentum !! for 2M microphysics, always output these variables - if (mp_phys == 10) then + if (mp_phys == mp_phys_mg) then if (advcld) then QLCN(ipt,l) = max(qli(l)-ccin(ipt,l,2), zero) QICN(ipt,l) = max(qii(l)-ccin(ipt,l,1), zero) @@ -1071,7 +1080,7 @@ subroutine rascnv_run(IM, k, ntr, dt, dtf & end subroutine rascnv_run SUBROUTINE CLOUD( & & K, KP1, KD, NTRC, KBLMX, kblmn & - &, FRACBL, MAX_NEG_BOUY, vsmooth, do_aw & + &, FRACBL, MAX_NEG_BOUY, vsmooth, aw_scal & &, REVAP, WRKFUN, CALKBL, CRTFUN & &, DT, KDT, TLA, DPD & &, ALFINT, RHFACL, RHFACS, area, ccwf, cd, trcfac & @@ -1145,7 +1154,8 @@ SUBROUTINE CLOUD( & &, RHRAM=0.05_kp & ! PBL RELATIVE HUMIDITY RAMP ! &, RHRAM=0.15_kp !& ! PBL RELATIVE HUMIDITY RAMP &, HCRITD=4000.0_kp & ! Critical Moist Static Energy for Deep clouds - &, HCRITS=2000.0_kp & ! Critical Moist Static Energy for Shallow clouds +! &, HCRITS=2000.0_kp & ! Critical Moist Static Energy for Shallow clouds + &, HCRITS=2500.0_kp & ! Critical Moist Static Energy for Shallow clouds &, pcrit_lcl=250.0_kp & ! Critical pressure difference between boundary layer top ! layer top and lifting condensation level (hPa) ! &, hpert_fac=1.01_kp !& ! Perturbation on hbl when ctei=.true. @@ -1153,7 +1163,9 @@ SUBROUTINE CLOUD( & &, qudfac=quad_lam*half & &, shalfac=3.0_kp & ! &, qudfac=quad_lam*pt25, shalfac=3.0_kp !& ! Yogesh's - &, c0ifac=0.07_kp & ! following Han et al, 2016 MWR +! &, c0ifac=0.07_kp & ! following Han et al, 2016 MWR +! &, c0ifac=0.001_kp & ! following Han et al, 2017 Weather and Forecasting + &, c0ifac=0.0_kp & &, dpnegcr = 150.0_kp ! &, dpnegcr = 100.0_kp ! &, dpnegcr = 200.0_kp @@ -1172,7 +1184,7 @@ SUBROUTINE CLOUD( & ! LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP, ctei LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP - logical vsmooth, do_aw + logical vsmooth, aw_scal INTEGER K, KP1, KD, NTRC, kblmx, kblmn, ntk @@ -1405,7 +1417,8 @@ SUBROUTINE CLOUD( & hmax = hol(kmax) elseif (kmax < k) then do l=kmax+1,k - if (abs(hol(kmax)-hol(l)) > half * hcrit) then +! if (abs(hol(kmax)-hol(l)) > half * hcrit) then + if (abs(hol(kmax)-hol(l)) > hcrit) then kmxb = l - 1 exit endif @@ -1435,7 +1448,6 @@ SUBROUTINE CLOUD( & endif enddo endif - ! klcl = kd1 if (kmax > kd1) then @@ -1446,6 +1458,7 @@ SUBROUTINE CLOUD( & endif enddo endif + ! if (klcl == kd .or. klcl < ktem) return ! This is to handle mid-level convection from quasi-uniform h @@ -1464,7 +1477,6 @@ SUBROUTINE CLOUD( & tem = min(50.0_kp,max(10.0_kp,(prl(kmaxp1)-prl(kd))*0.10_kp)) if (prl(kmaxp1) - prl(ii) > tem .and. ii > kbl) kbl = ii - if (kbl .ne. ii) then if (PRL(kmaxp1)-PRL(KBL) > bldmax) kbl = max(kbl,ii) endif @@ -1503,13 +1515,11 @@ SUBROUTINE CLOUD( & ! KBL = min(kmax,MAX(KBL,KD+2)) KB1 = KBL - 1 -!! - +! if (PRL(Kmaxp1)-PRL(KBL) > bldmax .or. kb1 <= kd ) then ! & .or. PRL(Kmaxp1)-PRL(KBL) < bldmin) then return endif -! ! PRIS = ONE / (PRL(KP1)-PRL(KBL)) PRISM = ONE / (PRL(Kmaxp1)-PRL(KBL)) @@ -1606,7 +1616,7 @@ SUBROUTINE CLOUD( & ENDDO ENDDO ! -! if (ntk > 0 .and. do_aw) then +! if (ntk > 0 .and. aw_scal) then if (ntk > 0) then if (rbl(ntk) > zero) then wcbase = min(two, max(wcbase, sqrt(twoo3*rbl(ntk)))) @@ -1671,7 +1681,8 @@ SUBROUTINE CLOUD( & QLL(KD ) = ALHF * GAF(KD) * QIL(KD) + ONE ! st1 = qil(kd) - st2 = c0i * st1 * exp(c0ifac*min(tol(kd)-t0c,zero)) + st2 = c0i * st1 + if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(kd)-t0c,zero)) tem = c0 * (one-st1) tem2 = st2*qi0 + tem*qw0 ! @@ -1693,7 +1704,8 @@ SUBROUTINE CLOUD( & AKC(L) = one / AKT(L) ! st1 = half * (qil(l)+qil(lp1)) - st2 = c0i * st1 * exp(c0ifac*min(tol(lp1)-t0c,zero)) + st2 = c0i * st1 + if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(lp1)-t0c,zero)) tem = c0 * (one-st1) tem2 = st2*qi0 + tem*qw0 ! @@ -1734,7 +1746,6 @@ SUBROUTINE CLOUD( & ALM = ALHF*QIL(KD) - LTL(KD) * VTF(KD) ! HSU = HST(KD) + LTL(KD) * NU * (QOL(KD)-QST(KD)) - ! !===> VERTICAL INTEGRALS NEEDED TO COMPUTE THE ENTRAINMENT PARAMETER ! @@ -1772,7 +1783,7 @@ SUBROUTINE CLOUD( & ! clp = one st2 = hbl - hsu - +! if (tx2 == zero) then alm = - st2 / tx1 if (alm > almax) alm = -100.0_kp @@ -1790,7 +1801,6 @@ SUBROUTINE CLOUD( & endif endif - ! ! CLIP CASE: ! NON-ENTRAINIG CLOUD DETRAINS IN LOWER HALF OF TOP LAYER. @@ -1887,7 +1897,6 @@ SUBROUTINE CLOUD( & ETAI(L) = one / ETA(L) ENDDO ETAI(KBL) = one - ! !===> CLOUD WORKFUNCTION ! @@ -2046,7 +2055,6 @@ SUBROUTINE CLOUD( & TEM = max(0.05_kp, MIN(CD*200.0_kp, MAX_NEG_BOUY)) IF (.not. cnvflg .and. WFN > ACR .and. & & dpneg < dpnegcr .and. AKM <= TEM) CALCUP = .TRUE. - ! !===> IF NO SOUNDING MEETS THIRD CONDITION, RETURN ! @@ -2145,7 +2153,6 @@ SUBROUTINE CLOUD( & GMS(KD) = (DS + st1 - tem1*det*alhl-tem*alhf) * PRI(KD) GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOS + DH) - ! ! TENDENCY FOR SUSPENDED ENVIRONMENTAL ICE AND/OR LIQUID WATER ! @@ -2185,7 +2192,6 @@ SUBROUTINE CLOUD( & GMH(L) = DH * PRI(L) GMS(L) = DS * PRI(L) - ! GHD(L) = TEM5 * PRI(L) GSD(L) = (TEM5 - ALHL * TEM6) * PRI(L) @@ -2241,7 +2247,6 @@ SUBROUTINE CLOUD( & ! avh = avh + tx1*(prs(l+1)-prs(l)) ENDDO - ! !*********************************************************************** !*********************************************************************** @@ -2303,8 +2308,7 @@ SUBROUTINE CLOUD( & ! hbl = hbl * hpert_fac ! qbl = qbl * hpert_fac ! endif - - + !*********************************************************************** !===> CLOUD WORKFUNCTION FOR MODIFIED SOUNDING, THEN KERNEL (AKM) @@ -2386,7 +2390,6 @@ SUBROUTINE CLOUD( & AMBMAX = (PRL(KMAXP1)-PRL(KBL))*(FRACBL*GRAVCON) AMB = MAX(MIN(AMB, AMBMAX),ZERO) - !*********************************************************************** !*************************RESULTS*************************************** !*********************************************************************** @@ -2400,7 +2403,7 @@ SUBROUTINE CLOUD( & ! tx1 = one - amb * eta(kd) / (rho(kd)*wvl(kd)) ! sigf(kd) = max(zero, min(one, tx1 * tx1)) ! endif - if (do_aw) then + if (aw_scal) then tx1 = (0.2_kp / max(alm, 1.0e-5_kp)) tx2 = one - min(one, pi * tx1 * tx1 / area) @@ -2426,6 +2429,10 @@ SUBROUTINE CLOUD( & else sigf(kd:k) = one endif + + tx1 = max(1.0e-6_kp, abs(gms(kd) * onebcp * sigf(kd))) + amb = min(tx1*amb, tfrac_max*toi(kd)) / tx1 + ! avt = zero avq = zero @@ -2530,7 +2537,7 @@ SUBROUTINE CLOUD( & clfrac = max(ZERO, min(half, rknob*clf(tem)*tem1)) cldfrd = clfrac - +! DO L=KD,KBL ! Testing on 20070926 ! for L=KD,K IF (L >= IDH .AND. DDFT) THEN @@ -2576,7 +2583,6 @@ SUBROUTINE CLOUD( & & TEM4 = POTEVAP * (one - EXP( tx4*TX1**0.57777778_kp )) ACTEVAP = MIN(TX1, TEM4*CLFRAC) - if (tx1 < rainmin*dt) actevap = min(tx1, potevap) ! tem4 = zero @@ -2606,7 +2612,6 @@ SUBROUTINE CLOUD( & ENDDO CUP = CUP + TX1 + DOF * AMB * sigf(kbl) ENDIF - ! ! Convective transport (mixing) of passive tracers ! @@ -3300,6 +3305,7 @@ SUBROUTINE DDRFT( & ! endif ELSE ERRQ = TX2 ! Further iteration ! + ! if (itr == itrmu .and. ERRQ > ERRMIN*10 & ! & .and. ntla == 1) ERRQ = 10.0 ENDIF @@ -3461,9 +3467,7 @@ SUBROUTINE DDRFT( & ! VT(1) = GMS(L-1) * QRP(L-1) ** 0.1364 VT(1) = GMS(L-1) * QRPF(QRP(L-1)) RNT = ROR(L-1) * (WVL(L-1)+VT(1))*QRP(L-1) - ! - ! TEM = MAX(ALM, 2.5E-4) * MAX(ETA(L), 1.0) TEM = MAX(ALM,ONE_M6) * MAX(ETA(L), ONE) ! TEM = MAX(ALM, 1.0E-5) * MAX(ETA(L), 1.0) @@ -3620,6 +3624,7 @@ SUBROUTINE DDRFT( & ENDIF ERRH = HOD(L) - TEM1 ERRQ = ABS(ERRH/HOD(L)) + ABS(ERRE/MAX(ETD(L),ONE_M5)) + DOF = DDZ VT(2) = QQQ ! @@ -3678,6 +3683,7 @@ SUBROUTINE DDRFT( & ! Compute Buoyancy TEM1 = WA(3) + (HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2))) & & * onebcp + TEM1 = TEM1 * (one + NU*QOD(L)) ROR(L) = CMPOR * PRL(L) / TEM1 TEM1 = TEM1 * DOFW @@ -3689,6 +3695,7 @@ SUBROUTINE DDRFT( & TEM1 = WVL(L) WVL(L) = VT(2) * (ETD(L-1)*WVL(L-1) - FACG & & * (BUY(L-1)*QRT(L-1)+BUY(L)*QRB(L-1))) + ! if (wvl(l) < zero) then ! WVL(L) = max(wvl(l), 0.1*tem1) @@ -3707,9 +3714,12 @@ SUBROUTINE DDRFT( & ! ERRQ = ERRQ + ABS(ERRW/MAX(WVL(L),ONE_M5)) + ! IF (ITR >= MIN(ITRMIN,ITRMD/2)) THEN IF (ITR >= MIN(ITRMND,ITRMD/2)) THEN + IF (ETD(L-1) == zero .AND. ERRQ > 0.2_kp) THEN + ROR(L) = BUD(KD) ETD(L) = zero WVL(L) = zero @@ -3878,7 +3888,6 @@ SUBROUTINE DDRFT( & if (tx5 > zero) idnm = idnm + 1 endif ENDIF - ! ! If downdraft properties are not obtainable, (i.e.solution does ! not converge) , no downdraft is assumed @@ -4095,7 +4104,7 @@ end subroutine qrabf SUBROUTINE SETVTP implicit none - real(kind=kind_phys), parameter :: vtpexp=-0.3636_kp, one=1.0_kp + real(kind=kind_phys), parameter :: vtpexp=-0.3636_kp real(kind=kind_phys) xinc,x,xmax,xmin integer jx ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/physics/rascnv.meta b/physics/rascnv.meta index f0ab36f19..4babf620d 100644 --- a/physics/rascnv.meta +++ b/physics/rascnv.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rascnv type = scheme - dependencies = + dependencies = funcphys.f90,machine.F ######################################################################## [ccpp-arg-table] @@ -422,7 +422,7 @@ standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_dimension,tracer_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers_for_convective_transport) type = real kind = kind_phys intent = inout diff --git a/physics/sfc_sice.f b/physics/sfc_sice.f index ab67f849e..081bbf48e 100644 --- a/physics/sfc_sice.f +++ b/physics/sfc_sice.f @@ -287,11 +287,11 @@ subroutine sfc_sice_run & q0 = min(qs1, q0) if (fice(i) < cimin) then - print *,'warning: ice fraction is low:', fice(i) +! print *,'warning: ice fraction is low:', fice(i) fice(i) = cimin tice(i) = tgice tskin(i)= tgice - print *,'fix ice fraction: reset it to:', fice(i) +! print *,'fix ice fraction: reset it to:', fice(i) endif ffw(i) = one - fice(i) @@ -362,9 +362,9 @@ subroutine sfc_sice_run & snowd(i) = min( snowd(i), hsmax ) if (snowd(i) > (2.0_kind_phys*hice(i))) then - print *, 'warning: too much snow :',snowd(i) +! print *, 'warning: too much snow :',snowd(i) snowd(i) = hice(i) + hice(i) - print *,'fix: decrease snow depth to:',snowd(i) +! print *,'fix: decrease snow depth to:',snowd(i) endif endif enddo diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index f473dc9bd..649d994a8 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -230,8 +230,8 @@ subroutine ugwpv1_gsldrag_init ( & if ( do_ugwp_v1 ) then call cires_ugwpv1_init (me, master, nlunit, logunit, jdat, con_pi, & - con_rerth, fn_nml2, lonr, latr, levs, ak, bk, & - con_p0, dtp, errmsg, errflg) + con_rerth, fn_nml2, input_nml_file, lonr, latr, & + levs, ak, bk, con_p0, dtp, errmsg, errflg) if (errflg/=0) return end if