From 73bbc9f1df683e4dab9c22d52d0319c8615ffab2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 26 Nov 2019 08:22:34 -0700 Subject: [PATCH 1/4] Merge gsd/develop into dtc/develop, squashed commit of the following: commit 7f530edd66132aa4d92e042a580c0aebf7e69662 Merge: e0d5f16 b492f2e Author: Dom Heinzeller Date: Thu Nov 21 15:40:20 2019 -0700 Merge pull request #356 from tanyasmirnova/ruc_land_ice_v1 Added the capability to use climatological LAI in RUC LSM commit b492f2efb2a33ce8fbc43518bf1fd6fee44574e2 Merge: bd32702 e0d5f16 Author: tanyasmirnova Date: Wed Nov 20 20:36:42 2019 +0000 Merge branch 'gsd/develop' of https://github.com/NCAR/ccpp-physics into ruc_land_ice_v1 commit bd32702bfd96f2d4bab4b25bfa408c4c7bd098cd Author: tanyasmirnova Date: Wed Nov 20 20:29:42 2019 +0000 Added the capability to use a Leaf Area Index (LAI) climatology in RUC LSM. Variables xlaixy and rdlai are added to the argument list of lsm_ruc_run. If rdlai=.true. in the physics namelist, then the LAI climatology will be passed into the RUC LSM and used instead of look-up table value for a given vegetation type. commit e0d5f16696a64333dc1920b060c43a8dde050c00 Merge: 660ede7 e4d291e Author: Dom Heinzeller Date: Sat Nov 2 05:47:40 2019 +0900 Merge pull request #349 from tanyasmirnova/ruc_land_ice_v1 This commit has a fix for a problem of cloud-radiation coupling with the use of MYNN PBL. commit e4d291e1b08ab68a7820f55921d0b5584d58944b Author: tanyasmirnova Date: Fri Nov 1 16:47:58 2019 +0000 This commit has a fix for a problem of cloud-radiation coupling with the use of MYNN PBL. The problem: the first call to the radiation happens before the first call to MYNN PBL, therefore CLDFRA_BL=0 in the first call to mynnrad_pre, and zero values are sent to array cldcov(:,:). When cloud cover is zero, the RRTMG radiation thinks that there are no clouds at all. The erroneous cloud-free LW and SW downward radiation fluxes affect the first hour of itegration, and cause siginificant cooling in the ploar regions, and too warm land surface temperature from cloud-free SW radiation. The fix: the fist call to mynnrad_pre should be skipped, so that cloud cover - cldcov(:,:) - is not overwritten by zero values of MYNN subgrid-clouds. In this case the initial cloud cover is computed in progcld5 from initial cloud water mixing ratio, relative humidity and specific humidity in the layer. Starting with the second call to the rrtmg radiation, the MYNN subgrid clouds are used. commit 660ede7a9f83a45f6141200cd951446fddf7f15e Merge: 4a17324 db9742d Author: Dom Heinzeller Date: Mon Oct 28 12:38:54 2019 +0900 Merge pull request #344 from tanyasmirnova/ruc_land_ice_v1 Sync RUC LSM code with the version used in RAP/HRRR commit db9742d5609912a7e3db99769665eded32668332 Author: tanyasmirnova Date: Thu Oct 24 22:14:13 2019 +0000 Sync the RUC LSM code with the version in RAPv5/HRRRv4. Some clean-up in sfc_drv_ruc.F90. commit 27eb0898682ca2dce1a8da32826fd7be561a5f68 Merge: fa3c1d3 4a17324 Author: tanyasmirnova Date: Thu Oct 24 22:03:14 2019 +0000 Merge branch 'gsd/develop' of https://github.com/NCAR/ccpp-physics into ruc_land_ice_v1 commit 4a17324ac9c7e9351da6527c541a2d110109f8a5 Merge: 543f640 3a28055 Author: Dom Heinzeller Date: Thu Oct 24 10:53:19 2019 +0900 Merge pull request #338 from haiqinli/gsd/develop-hli "to include GF updates in GSDv0beta4" commit 3a280556fd762f04cd5a2688c861546ce6c097ec Author: Haiqin.Li Date: Wed Oct 23 21:13:25 2019 +0000 "update to pass the ccpp_gsd_noah_repro regression test case" commit 0711b8288eb5825b77f20c4079b28235c03d4c86 Author: Haiqin.Li Date: Sun Oct 20 04:54:18 2019 +0000 "update to pass ccpp_gsd regression test" commit fa3c1d39aa8e5db07f571fff9f3348cd4c0a1423 Author: tanyasmirnova Date: Thu Oct 17 16:28:55 2019 +0000 1. Use fraction of frozen precipitation SR directly from Thompson MP. 2. Bug fix in liquid precipitation and frozen fraction - SRFLAG. This bug was producing 1.e-3 factor maller values of SRFLAG. 3. Modification to comment for precipitation in sfc_drv_ruc.F90 commit a59d416574e7c978da10d2d0f82920e46ec047e0 Author: Haiqin.Li Date: Sun Oct 13 20:40:44 2019 +0000 "clean the code" commit 4ca463ca07c1d381ccb8bd018761bc0adee0e526 Author: Haiqin.Li Date: Sun Oct 13 20:35:36 2019 +0000 "update input of imfdeepcnv following Dom's suggestions" commit 14c1c5bfedc69cb6466d7c8a3a98d0f34454f125 Author: Haiqin.Li Date: Fri Sep 27 18:04:33 2019 +0000 "to include GF updates in GSDv0beta4" --- physics/GFS_suite_interstitial.F90 | 7 +- physics/GFS_suite_interstitial.meta | 8 + physics/cu_gf_deep.F90 | 287 ++++++++++++- physics/cu_gf_driver.F90 | 600 +++++++++++++++------------- physics/cu_gf_driver.meta | 67 ++++ physics/module_MYNNrad_post.F90 | 7 + physics/module_MYNNrad_post.meta | 16 + physics/module_MYNNrad_pre.F90 | 8 +- physics/module_MYNNrad_pre.meta | 16 + physics/module_sf_ruclsm.F90 | 14 +- physics/sfc_drv_ruc.F90 | 40 +- physics/sfc_drv_ruc.meta | 14 + 12 files changed, 774 insertions(+), 310 deletions(-) diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 1df53ff12..20f51f99c 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -659,7 +659,7 @@ end subroutine GFS_suite_interstitial_4_finalize !> \section arg_table_GFS_suite_interstitial_4_run Argument Table !! \htmlinclude GFS_suite_interstitial_4_run.html !! - subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & + subroutine GFS_suite_interstitial_4_run (imfdeepcnv, 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, dqdti, errmsg, errflg) @@ -670,7 +670,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! interface variables - integer, intent(in) :: im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, & + integer, intent(in) :: imfdeepcnv, im, levs, 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 @@ -736,7 +736,8 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to gq0(i,k,ntcw) = clw(i,k,2) ! water enddo enddo - if (imp_physics == imp_physics_thompson) then +! if (imp_physics == imp_physics_thompson) then + if (imp_physics == imp_physics_thompson .and. imfdeepcnv /= 3) then if (ltaerosol) then do k=1,levs do i=1,im diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index 44696dcb0..2fa377c00 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1454,6 +1454,14 @@ [ccpp-arg-table] name = GFS_suite_interstitial_4_run type = scheme +[imfdeepcnv] + standard_name = flag_for_mass_flux_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent diff --git a/physics/cu_gf_deep.F90 b/physics/cu_gf_deep.F90 index 3e865c9ba..4afad80d1 100644 --- a/physics/cu_gf_deep.F90 +++ b/physics/cu_gf_deep.F90 @@ -14,7 +14,7 @@ module cu_gf_deep !> tuning constant for cloudwater/ice detrainment real(kind=kind_phys), parameter:: c1= 0.003 !.002 ! .0005 !> parameter to turn on or off evaporation of rainwater as done in sas - integer, parameter :: irainevap=0 + integer, parameter :: irainevap=1 !> max allowed fractional coverage (frh_thresh) real(kind=kind_phys), parameter::frh_thresh = .9 !> rh threshold. if fractional coverage ~ frh_thres, do not use cupa any further @@ -362,7 +362,7 @@ subroutine cu_gf_deep_run( & c1_max=c1 elocp=xlv/cp el2orc=xlv*xlv/(r_v*cp) - evfact=.2 + evfact=.4 ! .2 evfactl=.2 !evfact=.0 ! for 4F5f !evfactl=.4 @@ -1923,6 +1923,13 @@ subroutine cu_gf_deep_run( & ichoice,imid,ipr,itf,ktf, & its,ite, kts,kte, & dicycle,xf_dicycle ) + +!---------------evap below cloud base + + call rain_evap_below_cloudbase(itf,ktf,its,ite, & + kts,kte,ierr,kbcon,xmb,psur,xland,qo_cup, & + po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy) + k=1 do i=its,itf if(ierr(i).eq.0 .and.pre(i).gt.0.) then @@ -1971,7 +1978,7 @@ subroutine cu_gf_deep_run( & do k = ktop(i), 1, -1 rain = pwo(i,k) + edto(i) * pwdo(i,k) rn(i) = rn(i) + rain * xmb(i) * .001 * dtime - !if(po(i,k).gt.700.)then + if(po(i,k).gt.400.)then if(flg(i))then q1=qo(i,k)+(outq(i,k))*dtime t1=tn(i,k)+(outt(i,k))*dtime @@ -1996,7 +2003,7 @@ subroutine cu_gf_deep_run( & pre(i)=max(pre(i),0.) delqev(i) = delqev(i) + .001*dp*qevap(i)/g endif - !endif ! 700mb + endif ! 400mb endif enddo ! pre(i)=1000.*rn(i)/dtime @@ -2035,6 +2042,271 @@ end subroutine cu_gf_deep_run !> @} !>\ingroup cu_gf_deep_group + + + subroutine fct1d3 (ktop,n,dt,z,tracr,massflx,trflx_in,dellac,g) + +! --- modify a 1-D array of tracer fluxes for the purpose of maintaining +! --- monotonicity (including positive-definiteness) in the tracer field +! --- during tracer transport. + +! --- the underlying transport equation is (d tracr/dt) = - (d trflx/dz) +! --- where dz = |z(k+1)-z(k)| (k=1,...,n) and trflx = massflx * tracr +! --- physical dimensions of tracr,trflx,dz are arbitrary to some extent +! --- but are subject to the constraint dim[trflx] = dim[tracr*(dz/dt)]. + +! --- note: tracr is carried in grid cells while z and fluxes are carried on +! --- interfaces. interface variables at index k are at grid location k-1/2. +! --- sign convention: mass fluxes are considered positive in +k direction. + +! --- massflx and trflx_in must be provided independently to allow the +! --- algorithm to generate an auxiliary low-order (diffusive) tracer flux +! --- as a stepping stone toward the final product trflx_out. + + implicit none + integer,intent(in) :: n,ktop ! number of grid cells + real(kind=kind_phys) ,intent(in) :: dt,g ! transport time step + real(kind=kind_phys) ,intent(in) :: z(n+0) ! location of cell interfaces + real(kind=kind_phys) ,intent(in) :: tracr(n) ! the transported variable + real(kind=kind_phys) ,intent(in) :: massflx(n+0) ! mass flux across interfaces + real(kind=kind_phys) ,intent(in) :: trflx_in(n+0) ! original tracer flux + real(kind=kind_phys) ,intent(out):: dellac(n+0) ! modified tracr flux + real(kind=kind_phys) :: trflx_out(n+0) ! modified tracr flux + integer k,km1,kp1 + logical :: NaN, error=.false., vrbos=.true. + real(kind=kind_phys) dtovdz(n),trmax(n),trmin(n),flx_lo(n+0),antifx(n+0),clipped(n+0), & + soln_hi(n),totlin(n),totlout(n),soln_lo(n),clipin(n),clipout(n),arg + real(kind=kind_phys),parameter :: epsil=1.e-22 ! prevent division by zero + real(kind=kind_phys),parameter :: damp=1. ! damper of antidff flux (1=no damping) + NaN(arg) = .not. (arg.ge.0. .or. arg.le.0.) ! NaN detector + dtovdz(:)=0. + soln_lo(:)=0. + antifx(:)=0. + clipin(:)=0. + totlin(:)=0. + totlout(:)=0. + clipout(:)=0. + flx_lo(:)=0. + trmin(:)=0. + trmax(:)=0. + clipped(:)=0. + trflx_out(:)=0. + do k=1,ktop + dtovdz(k)=.01*dt/abs(z(k+1)-z(k))*g ! time step / grid spacing + if (z(k).eq.z(k+1)) error=.true. + end do +! if (vrbos .or. error) print '(a/(8es10.3))','(fct1d) dtovdz =',dtovdz + + do k=2,ktop + if (massflx(k).ge.0.) then + flx_lo(k)=massflx(k)*tracr(k-1) ! low-order flux, upstream + else + flx_lo(k)=massflx(k)*tracr(k) ! low-order flux, upstream + end if + antifx(k)=trflx_in(k)-flx_lo(k) ! antidiffusive flux + end do + flx_lo( 1)=trflx_in( 1) + flx_lo(ktop+1)=trflx_in(ktop+1) + antifx( 1)=0. + antifx(ktop+1)=0. +! --- clip low-ord fluxes to make sure they don't violate positive-definiteness + do k=1,ktop + totlout(k)=max(0.,flx_lo(k+1))-min(0.,flx_lo(k )) ! total flux out + clipout(k)=min(1.,tracr(k)/max(epsil,totlout(k))/ (1.0001*dtovdz(k))) + end do + + do k=2,ktop + if (massflx(k).ge.0.) then + flx_lo(k)=flx_lo(k)*clipout(k-1) + else + flx_lo(k)=flx_lo(k)*clipout(k) + end if + end do + if (massflx( 1).lt.0.) flx_lo( 1)=flx_lo( 1)*clipout(1) + if (massflx(ktop+1).gt.0.)flx_lo(ktop+1)=flx_lo(ktop+1)*clipout(ktop) + +! --- a positive-definite low-order (diffusive) solution can now be constructed + + do k=1,ktop + soln_lo(k)=tracr(k)-(flx_lo(k+1)-flx_lo(k))*dtovdz(k) ! low-ord solutn + dellac(k)=-(flx_lo(k+1)-flx_lo(k))*dtovdz(k)/dt + !dellac(k)=soln_lo(k) + end do + return + do k=1,ktop + km1=max(1,k-1) + kp1=min(ktop,k+1) + trmax(k)= max(soln_lo(km1),soln_lo(k),soln_lo(kp1), & + tracr (km1),tracr (k),tracr (kp1)) ! upper bound + trmin(k)=max(0.,min(soln_lo(km1),soln_lo(k),soln_lo(kp1), & + tracr (km1),tracr (k),tracr (kp1))) ! lower bound + end do + + do k=1,ktop + totlin (k)=max(0.,antifx(k ))-min(0.,antifx(k+1)) ! total flux in + totlout(k)=max(0.,antifx(k+1))-min(0.,antifx(k )) ! total flux out + + clipin (k)=min(damp,(trmax(k)-soln_lo(k))/max(epsil,totlin (k)) & + / (1.0001*dtovdz(k))) + clipout(k)=min(damp,(soln_lo(k)-trmin(k))/max(epsil,totlout(k)) & + / (1.0001*dtovdz(k))) + + if (NaN(clipin(k))) print *,'(fct1d) error: clipin is NaN, k=',k + if (NaN(clipout(k))) print *,'(fct1d) error: clipout is NaN, k=',k + + if (clipin(k).lt.0.) then +! print 100,'(fct1d) error: clipin < 0 at k =',k, & +! 'clipin',clipin(k),'trmax',trmax(k),'soln_lo',soln_lo(k), & +! 'totlin',totlin(k),'dt/dz',dtovdz(k) + error=.true. + end if + if (clipout(k).lt.0.) then +! print 100,'(fct1d) error: clipout < 0 at k =',k, & +! 'clipout',clipout(k),'trmin',trmin(k),'soln_lo',soln_lo(k), & +! 'totlout',totlout(k),'dt/dz',dtovdz(k) + error=.true. + end if +! 100 format (a,i3/(4(a10,"=",es9.2))) + end do + + do k=2,ktop + if (antifx(k).gt.0.) then + clipped(k)=antifx(k)*min(clipout(k-1),clipin(k)) + else + clipped(k)=antifx(k)*min(clipout(k),clipin(k-1)) + end if + trflx_out(k)=flx_lo(k)+clipped(k) + if (NaN(trflx_out(k))) then + print *,'(fct1d) error: trflx_out is NaN, k=',k + error=.true. + end if + end do + trflx_out( 1)=trflx_in( 1) + trflx_out(ktop+1)=trflx_in(ktop+1) + do k=1,ktop + soln_hi(k)=tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k) + dellac(k)=-g*(trflx_out(k+1)-trflx_out(k))*dtovdz(k)/dt + !dellac(k)=soln_hi(k) + end do + + if (vrbos .or. error) then +! do k=2,ktop +! write(32,99)k, & +! 'tracr(k)', tracr(k), & +! 'flx_in(k)', trflx_in(k), & +! 'flx_in(k+1)', trflx_in(k+1), & +! 'flx_lo(k)', flx_lo(k), & +! 'flx_lo(k+1)', flx_lo(k+1), & +! 'soln_lo(k)', soln_lo(k), & +! 'trmin(k)', trmin(k), & +! 'trmax(k)', trmax(k), & +! 'totlin(k)', totlin(k), & +! 'totlout(k)', totlout(k), & +! 'clipin(k-1)', clipin(k-1), & +! 'clipin(k)', clipin(k), & +! 'clipout(k-1)', clipout(k-1), & +! 'clipout(k)', clipout(k), & +! 'antifx(k)', antifx(k), & +! 'antifx(k+1)', antifx(k+1), & +! 'clipped(k)', clipped(k), & +! 'clipped(k+1)', clipped(k+1), & +! 'flx_out(k)', trflx_out(k), & +! 'flx_out(k+1)', trflx_out(k+1), & +! 'dt/dz(k)', dtovdz(k), & +! 'final', tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k) +! 99 format ('(trc1d) k =',i4/(3(a13,'=',es13.6))) +! end do + if (error) stop '(fct1d error)' + end if + + return + end subroutine fct1d3 + + subroutine rain_evap_below_cloudbase(itf,ktf, its,ite, kts,kte,ierr, & + kbcon,xmb,psur,xland,qo_cup, & + po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy) + + implicit none + real(kind=kind_phys), parameter :: alp1=5.44e-4 & !1/sec + ,alp2=5.09e-3 & !unitless + ,alp3=0.5777 & !unitless + ,c_conv=0.05 !conv fraction area, unitless + + + integer ,intent(in) :: itf,ktf, its,ite, kts,kte + integer, dimension(its:ite) ,intent(in) :: ierr,kbcon + real(kind=kind_phys), dimension(its:ite) ,intent(in) ::psur,xland,pwavo,edto,pwevo,xmb + real(kind=kind_phys), dimension(its:ite,kts:kte),intent(in) :: po_cup,qo_cup,qes_cup + real(kind=kind_phys), dimension(its:ite) ,intent(inout) :: pre + real(kind=kind_phys), dimension(its:ite,kts:kte),intent(inout) :: outt,outq !,outbuoy + + !real, dimension(its:ite) ,intent(out) :: tot_evap_bcb + !real, dimension(its:ite,kts:kte),intent(out) :: evap_bcb,net_prec_bcb + + !-- locals + integer :: i,k + real(kind=kind_phys) :: RH_cr , del_t,del_q,dp,q_deficit + real(kind=kind_phys), dimension(its:ite,kts:kte) :: evap_bcb,net_prec_bcb + real(kind=kind_phys), dimension(its:ite) :: tot_evap_bcb + + do i=its,itf + evap_bcb (i,:)= 0.0 + net_prec_bcb(i,:)= 0.0 + tot_evap_bcb(i) = 0.0 + if(ierr(i) /= 0) cycle + + !-- critical rel humidity + RH_cr=0.9*xland(i)+0.7*(1-xland(i)) + !RH_cr=1. + + !-- net precipitation (after downdraft evap) at cloud base, available to + !evap + k=kbcon(i) + !net_prec_bcb(i,k) = xmb(i)*(pwavo(i)+edto(i)*pwevo(i)) !-- pwevo<0. + net_prec_bcb(i,k) = pre(i) + + do k=kbcon(i)-1, kts, -1 + + q_deficit = max(0.,(RH_cr*qes_cup(i,k) -qo_cup(i,k))) + + if(q_deficit < 1.e-6) then + net_prec_bcb(i,k)= net_prec_bcb(i,k+1) + cycle + endif + + dp = 100.*(po_cup(i,k)-po_cup(i,k+1)) + + !--units here: kg[water]/kg[air}/sec + evap_bcb(i,k) = c_conv * alp1 * q_deficit * & + ( sqrt(po_cup(i,k)/psur(i))/alp2 *net_prec_bcb(i,k+1)/c_conv )**alp3 + + !--units here: kg[water]/kg[air}/sec * kg[air]/m3 * m = kg[water]/m2/sec + evap_bcb(i,k)= evap_bcb(i,k)*dp/g + + if((net_prec_bcb(i,k+1) - evap_bcb(i,k)).lt.0.) cycle + if((pre(i) - evap_bcb(i,k)).lt.0.) cycle + net_prec_bcb(i,k)= net_prec_bcb(i,k+1) - evap_bcb(i,k) + + tot_evap_bcb(i) = tot_evap_bcb(i)+evap_bcb(i,k) + + !-- feedback + del_q = evap_bcb(i,k)*g/dp ! > 0., units: kg[water]/kg[air}/sec + del_t = -evap_bcb(i,k)*g/dp*(xlv/cp) ! < 0., units: K/sec + +! print*,"ebcb2",k,del_q*86400,del_t*86400 + + outq (i,k) = outq (i,k) + del_q + outt (i,k) = outt (i,k) + del_t + !outbuoy(i,k) = outbuoy(i,k) + cp*del_t+xlv*del_q + + pre(i) = pre(i) - evap_bcb(i,k) + enddo + enddo + + end subroutine rain_evap_below_cloudbase + + + subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & pw,ccn,pwev,edtmax,edtmin,edtc,psum2,psumh, & rho,aeroevap,itf,ktf, & @@ -2747,9 +3019,8 @@ subroutine cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2 xff_ens3(12)=0. xff_ens3(13)= 0. xff_ens3(16)= 0. -! closure_n(i)=12. -! hli 05/01/2018 closure_n(i)=12. -! xff_dicycle = 0. +! closure_n(i)=12. +! xff_dicycle = 0. endif !xff0 endif ! ichoice @@ -3682,7 +3953,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & prop_b(kts:kte)=0 iall=0 c0=.002 - clwdet=100. + clwdet=50. bdsp=bdispm ! !--- no precip for small clouds diff --git a/physics/cu_gf_driver.F90 b/physics/cu_gf_driver.F90 index 58a30749a..53e26fb46 100644 --- a/physics/cu_gf_driver.F90 +++ b/physics/cu_gf_driver.F90 @@ -7,8 +7,9 @@ module cu_gf_driver ! DH* TODO: replace constants with arguments to cu_gf_driver_run use physcons , g => con_g, cp => con_cp, xlv => con_hvap, r_v => con_rv use machine , only: kind_phys - use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap + use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3 use cu_gf_sh , only: cu_gf_sh_run + use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber implicit none @@ -68,11 +69,12 @@ end subroutine cu_gf_driver_finalize !! !>\section gen_gf_driver GSD GF Cumulus Scheme General Algorithm !> @{ - subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & - forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d, & - us,vs,t2di,w,qv2di_spechum,p2di,psuri, & - hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & - pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & + subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & + forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d, & + us,vs,t2di,w,qv2di_spechum,p2di,psuri, & + hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & + pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & + nwfa,con_rd,gq0,ntinc,ntlnc,imp_physics,imp_physics_thompson, & errmsg,errflg) !------------------------------------------------------------- implicit none @@ -94,7 +96,7 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & integer :: ishallow_g3 ! depend on imfshalcnv !------------------------------------------------------------- integer :: its,ite, jts,jte, kts,kte - integer, intent(in ) :: im,ix,km + integer, intent(in ) :: im,ix,km,ntracer real(kind=kind_phys), dimension( ix , km ), intent(in ) :: forcet,forceqv_spechum,w,phil real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: t,us,vs @@ -104,16 +106,16 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & real(kind=kind_phys), dimension( ix , km ), intent(out ) :: cnvw_moist,cnvc real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: cliw, clcw -!hj change from ix to im +! change from ix to im integer, dimension (im), intent(inout) :: hbot,htop,kcnv integer, dimension (im), intent(in) :: xland real(kind=kind_phys), dimension (im), intent(in) :: pbl integer, dimension (ix) :: tropics -! ruc variable +! ruc variable real(kind=kind_phys), dimension (im) :: hfx2,qfx2,psuri real(kind=kind_phys), dimension (im,km) :: ud_mf,dd_mf,dt_mf real(kind=kind_phys), dimension (im), intent(inout) :: raincv,cld1d -!hj end change ix to im +! end change ix to im real(kind=kind_phys), dimension (ix,km) :: t2di,p2di ! Specific humidity from FV3 real(kind=kind_phys), dimension (ix,km), intent(in) :: qv2di_spechum @@ -123,80 +125,76 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! real(kind=kind_phys), dimension( im ),intent(in) :: garea real(kind=kind_phys), intent(in ) :: dt + +! additional variables for number concentrations + real(kind=kind_phys), intent(in) :: nwfa(1:im,1:km) + real(kind=kind_phys), intent(in) :: con_rd + real(kind=kind_phys), dimension(im,km,ntracer), intent(inout) :: gq0 + integer, intent(in) :: imp_physics,imp_physics_thompson,ntlnc,ntinc + integer, intent(in ) :: imfshalcnv character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg -!hj define locally for now. - integer, dimension(im),intent(inout) :: cactiv ! hli for gf -!hj change from ix to im +! define locally for now. + integer, dimension(im),intent(inout) :: cactiv integer, dimension(im) :: k22_shallow,kbcon_shallow,ktop_shallow real(kind=kind_phys), dimension(im) :: ht -!hj change -! -!+lxz -!hj real(kind=kind_phys) :: dx real(kind=kind_phys), dimension(im) :: dx -! local vars -!hj change ix to im - real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws - real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm - real(kind=kind_phys), dimension (im,km) :: outts,outqs,outqcs,outu,outv,outus,outvs - real(kind=kind_phys), dimension (im,km) :: outtm,outqm,outqcm,submm,cupclwm - real(kind=kind_phys), dimension (im,km) :: cnvwt,cnvwts,cnvwtm - real(kind=kind_phys), dimension (im,km) :: hco,hcdo,zdo,zdd,hcom,hcdom,zdom - real(kind=kind_phys), dimension (km) :: zh - real(kind=kind_phys), dimension (im) :: tau_ecmwf,edt,edtm,edtd,ter11,aa0,xlandi - real(kind=kind_phys), dimension (im) :: pret,prets,pretm,hexec - real(kind=kind_phys), dimension (im,10) :: forcing,forcing2 -!+lxz - integer, dimension (im) :: kbcon, ktop,ierr,ierrs,ierrm,kpbli - integer, dimension (im) :: k22s,kbcons,ktops,k22,jmin,jminm - integer, dimension (im) :: kbconm,ktopm,k22m -!hj end change ix to im -!.lxz - integer :: iens,ibeg,iend,jbeg,jend,n - integer :: ibegh,iendh,jbegh,jendh - integer :: ibegc,iendc,jbegc,jendc,kstop - real(kind=kind_phys) :: rho_dryar,temp - real(kind=kind_phys) :: pten,pqen,paph,zrho,pahfs,pqhfl,zkhvfl,pgeoh -!hj 10/11/2016: ipn is an input in fim. set it to zero here. - integer, parameter :: ipn = 0 + real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws + real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm + real(kind=kind_phys), dimension (im,km) :: outts,outqs,outqcs,outu,outv,outus,outvs + real(kind=kind_phys), dimension (im,km) :: outtm,outqm,outqcm,submm,cupclwm + real(kind=kind_phys), dimension (im,km) :: cnvwt,cnvwts,cnvwtm + real(kind=kind_phys), dimension (im,km) :: hco,hcdo,zdo,zdd,hcom,hcdom,zdom + real(kind=kind_phys), dimension (km) :: zh + real(kind=kind_phys), dimension (im) :: tau_ecmwf,edt,edtm,edtd,ter11,aa0,xlandi + real(kind=kind_phys), dimension (im) :: pret,prets,pretm,hexec + real(kind=kind_phys), dimension (im,10) :: forcing,forcing2 + + integer, dimension (im) :: kbcon, ktop,ierr,ierrs,ierrm,kpbli + integer, dimension (im) :: k22s,kbcons,ktops,k22,jmin,jminm + integer, dimension (im) :: kbconm,ktopm,k22m + + integer :: iens,ibeg,iend,jbeg,jend,n + integer :: ibegh,iendh,jbegh,jendh + integer :: ibegc,iendc,jbegc,jendc,kstop + real(kind=kind_phys), dimension(im,km) :: rho_dryar + real(kind=kind_phys) :: pten,pqen,paph,zrho,pahfs,pqhfl,zkhvfl,pgeoh + integer, parameter :: ipn = 0 ! ! basic environmental input includes moisture convergence (mconv) ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off ! convection for this call only and at that particular gridpoint ! -!hj 10/11/2016: change ix to im. - real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi - real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg - real(kind=kind_phys), dimension (im) :: ccn,z1,psur,cuten,cutens,cutenm - real(kind=kind_phys), dimension (im) :: umean,vmean,pmean - real(kind=kind_phys), dimension (im) :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv -!hj end change ix to im - - integer :: i,j,k,icldck,ipr,jpr,jpr_deep,ipr_deep - integer :: itf,jtf,ktf,iss,jss,nbegin,nend - integer :: high_resolution - real(kind=kind_phys) :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter - real(kind=kind_phys) :: dsubclw,dsubclws,dsubclwm,ztm,ztq,hfm,qfm,rkbcon,rktop !-lxz -!hj change ix to im - real(kind=kind_phys), dimension (im) :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep - character*50 :: ierrc(im),ierrcm(im) - character*50 :: ierrcs(im) -!hj end change ix to im -! ruc variable -!hj hfx2 -- sensible heat flux (k m/s), positive upward from sfc -!hj qfx2 -- latent heat flux (kg/kg m/s), positive upward from sfc -!hj gf needs them in w/m2. define hfx and qfx after simple unit conversion - real(kind=kind_phys), dimension (im) :: hfx,qfx - real(kind=kind_phys) tem,tem1,tf,tcr,tcrf - - parameter (tf=243.16, tcr=270.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=263.16, tcr=273.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf)) ! as fim - ! initialize ccpp error handling variables + real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi + real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg + real(kind=kind_phys), dimension (im) :: ccn,z1,psur,cuten,cutens,cutenm + real(kind=kind_phys), dimension (im) :: umean,vmean,pmean + real(kind=kind_phys), dimension (im) :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv + + integer :: i,j,k,icldck,ipr,jpr,jpr_deep,ipr_deep + integer :: itf,jtf,ktf,iss,jss,nbegin,nend + integer :: high_resolution + real(kind=kind_phys) :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter + real(kind=kind_phys) :: dsubclw,dsubclws,dsubclwm,dtime_max,ztm,ztq,hfm,qfm,rkbcon,rktop + real(kind=kind_phys), dimension(km) :: massflx,trcflx_in1,clw_in1,clw_ten1,po_cup +! real(kind=kind_phys), dimension(km) :: trcflx_in2,clw_in2,clw_ten2 + real(kind=kind_phys), dimension (im) :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep + character*50 :: ierrc(im),ierrcm(im) + character*50 :: ierrcs(im) +! ruc variable +! hfx2 -- sensible heat flux (k m/s), positive upward from sfc +! qfx2 -- latent heat flux (kg/kg m/s), positive upward from sfc +! gf needs them in w/m2. define hfx and qfx after simple unit conversion + real(kind=kind_phys), dimension (im) :: hfx,qfx + real(kind=kind_phys) tem,tem1,tf,tcr,tcrf + + parameter (tf=243.16, tcr=270.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=263.16, tcr=273.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf)) ! as fim + ! initialize ccpp error handling variables errmsg = '' errflg = 0 ! @@ -212,8 +210,7 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! ! these should be coming in from outside ! -! print*,'hli in gf cactiv',cactiv -! cactiv(:) = 0 +! cactiv(:) = 0 rand_mom(:) = 0. rand_vmas(:) = 0. rand_clos(:,:) = 0. @@ -232,112 +229,113 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! !> - Set tuning constants for radiation coupling ! - tun_rad_shall(:)=.02 - tun_rad_mid(:)=.15 - tun_rad_deep(:)=.13 - edt(:)=0. - edtm(:)=0. - edtd(:)=0. - zdd(:,:)=0. - flux_tun(:)=5. -!hj 10/11/2016 dx and tscl_kf are replaced with input dx(i), is dlength. - ! dx for scale awareness -!hj dx=40075000./float(lonf) -!hj tscl_kf=dx/25000. - ccn(its:ite)=150. - ! - if (imfshalcnv == 3) then - ishallow_g3 = 1 - else - ishallow_g3 = 0 - end if - high_resolution=0 - subcenter=0. - iens=1 + tun_rad_shall(:)=.02 + tun_rad_mid(:)=.15 + tun_rad_deep(:)=.13 + edt(:)=0. + edtm(:)=0. + edtd(:)=0. + zdd(:,:)=0. + flux_tun(:)=5. +! 10/11/2016 dx and tscl_kf are replaced with input dx(i), is dlength. +! dx for scale awareness +! dx=40075000./float(lonf) +! tscl_kf=dx/25000. + ccn(its:ite)=150. + + if (imfshalcnv == 3) then + ishallow_g3 = 1 + else + ishallow_g3 = 0 + end if + high_resolution=0 + subcenter=0. + iens=1 ! ! these can be set for debugging ! - ipr=0 - jpr=0 - ipr_deep=0 - jpr_deep= 0 !53322 ! 528196 !0 ! 1136 !0 !421755 !3536 + ipr=0 + jpr=0 + ipr_deep=0 + jpr_deep= 0 !53322 ! 528196 !0 ! 1136 !0 !421755 !3536 ! ! - ibeg=its - iend=ite - tcrit=258. - - ztm=0. - ztq=0. - hfm=0. - qfm=0. - ud_mf =0. - dd_mf =0. - dt_mf =0. - tau_ecmwf(:)=0. + ibeg=its + iend=ite + tcrit=258. + + ztm=0. + ztq=0. + hfm=0. + qfm=0. + ud_mf =0. + dd_mf =0. + dt_mf =0. + tau_ecmwf(:)=0. ! - j=1 - ht(:)=phil(:,1)/g - do i=its,ite - cld1d(i)=0. - zo(i,:)=phil(i,:)/g - dz8w(i,1)=zo(i,2)-zo(i,1) - zh(1)=0. - kpbli(i)=2 - do k=kts+1,ktf - dz8w(i,k)=zo(i,k+1)-zo(i,k) - enddo - do k=kts+1,ktf - zh(k)=zh(k-1)+dz8w(i,k-1) - if(zh(k).gt.pbl(i))then - kpbli(i)=max(2,k) - exit - endif - enddo - enddo - do i= its,itf - forcing(i,:)=0. - forcing2(i,:)=0. - ccn(i)=100. - hbot(i) =kte - htop(i) =kts - raincv(i)=0. - xlandi(i)=real(xland(i)) -! if(abs(xlandi(i)-1.).le.1.e-3) tun_rad_shall(i)=.15 -! if(abs(xlandi(i)-1.).le.1.e-3) flux_tun(i)=1.5 + j=1 + ht(:)=phil(:,1)/g + do i=its,ite + cld1d(i)=0. + zo(i,:)=phil(i,:)/g + dz8w(i,1)=zo(i,2)-zo(i,1) + zh(1)=0. + kpbli(i)=2 + do k=kts+1,ktf + dz8w(i,k)=zo(i,k+1)-zo(i,k) + enddo + do k=kts+1,ktf + zh(k)=zh(k-1)+dz8w(i,k-1) + if(zh(k).gt.pbl(i))then + kpbli(i)=max(2,k) + exit + endif + enddo enddo + do i= its,itf - mconv(i)=0. + forcing(i,:)=0. + forcing2(i,:)=0. + ccn(i)=100. + hbot(i) =kte + htop(i) =kts + raincv(i)=0. + xlandi(i)=real(xland(i)) +! if(abs(xlandi(i)-1.).le.1.e-3) tun_rad_shall(i)=.15 +! if(abs(xlandi(i)-1.).le.1.e-3) flux_tun(i)=1.5 enddo - do k=kts,kte do i= its,itf - omeg(i,k)=0. - zu(i,k)=0. - zum(i,k)=0. - zus(i,k)=0. - zd(i,k)=0. - zdm(i,k)=0. + mconv(i)=0. enddo + do k=kts,kte + do i= its,itf + omeg(i,k)=0. + zu(i,k)=0. + zum(i,k)=0. + zus(i,k)=0. + zd(i,k)=0. + zdm(i,k)=0. + enddo enddo psur(:)=0.01*psuri(:) do i=its,itf - ter11(i)=max(0.,ht(i)) + ter11(i)=max(0.,ht(i)) enddo do k=kts,kte - do i=its,ite - cnvw(i,k)=0. - cnvc(i,k)=0. - gdc(i,k,1)=0. - gdc(i,k,2)=0. - gdc(i,k,3)=0. - gdc(i,k,4)=0. - gdc(i,k,7)=0. - gdc(i,k,8)=0. - gdc(i,k,9)=0. - gdc(i,k,10)=0. - gdc2(i,k,1)=0. - enddo + do i=its,ite + cnvw(i,k)=0. + cnvc(i,k)=0. + gdc(i,k,1)=0. + gdc(i,k,2)=0. + gdc(i,k,3)=0. + gdc(i,k,4)=0. + gdc(i,k,7)=0. + gdc(i,k,8)=0. + gdc(i,k,9)=0. + gdc(i,k,10)=0. + gdc2(i,k,1)=0. + enddo enddo ierr(:)=0 ierrm(:)=0 @@ -410,88 +408,80 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & subm(:,:)=0. dhdt(:,:)=0. - !print*,'hli t2di',t2di - !print*,'hli forcet',forcet do k=kts,ktf - do i=its,itf - p2d(i,k)=0.01*p2di(i,k) - po(i,k)=p2d(i,k) !*.01 - rhoi(i,k) = 100.*p2d(i,k)/(287.04*(t2di(i,k)*(1.+0.608*qv2di(i,k)))) - qcheck(i,k)=qv(i,k) - tn(i,k)=t(i,k)!+forcet(i,k)*dt - qo(i,k)=max(1.e-16,qv(i,k))!+forceqv(i,k)*dt - t2d(i,k)=t2di(i,k)-forcet(i,k)*dt - !print*,'hli t2di(i,k),forcet(i,k),dt,t2d(i,k)',t2di(i,k),forcet(i,k),dt,t2d(i,k) - q2d(i,k)=max(1.e-16,qv2di(i,k)-forceqv(i,k)*dt) - if(qo(i,k).lt.1.e-16)qo(i,k)=1.e-16 - tshall(i,k)=t2d(i,k) - qshall(i,k)=q2d(i,k) -!hj if(ipn.eq.jpr_deep)then -!hj write(12,123)k,dt,p2d(i,k),t2d(i,k),tn(i,k),q2d(i,k),qo(i,k),forcet(i,k) -!hj endif - enddo + do i=its,itf + p2d(i,k)=0.01*p2di(i,k) + po(i,k)=p2d(i,k) !*.01 + rhoi(i,k) = 100.*p2d(i,k)/(287.04*(t2di(i,k)*(1.+0.608*qv2di(i,k)))) + qcheck(i,k)=qv(i,k) + tn(i,k)=t(i,k)!+forcet(i,k)*dt + qo(i,k)=max(1.e-16,qv(i,k))!+forceqv(i,k)*dt + t2d(i,k)=t2di(i,k)-forcet(i,k)*dt + q2d(i,k)=max(1.e-16,qv2di(i,k)-forceqv(i,k)*dt) + if(qo(i,k).lt.1.e-16)qo(i,k)=1.e-16 + tshall(i,k)=t2d(i,k) + qshall(i,k)=q2d(i,k) + enddo enddo 123 format(1x,i2,1x,2(1x,f8.0),1x,2(1x,f8.3),3(1x,e13.5)) do i=its,itf - do k=kts,kpbli(i) + do k=kts,kpbli(i) tshall(i,k)=t(i,k) qshall(i,k)=max(1.e-16,qv(i,k)) - enddo + enddo enddo ! -!hj converting hfx2 and qfx2 to w/m2 -!hj hfx=cp*rho*hfx2 -!hj qfx=xlv*qfx2 +! converting hfx2 and qfx2 to w/m2 +! hfx=cp*rho*hfx2 +! qfx=xlv*qfx2 do i=its,itf - hfx(i)=hfx2(i)*cp*rhoi(i,1) - qfx(i)=qfx2(i)*xlv*rhoi(i,1) - dx(i) = sqrt(garea(i)) - !print*,'hli dx', dx(i) + hfx(i)=hfx2(i)*cp*rhoi(i,1) + qfx(i)=qfx2(i)*xlv*rhoi(i,1) + dx(i) = sqrt(garea(i)) enddo -!hj write(0,*),'hfx',hfx(3),qfx(3),rhoi(3,1) -!hj + do i=its,itf - do k=kts,kpbli(i) - tn(i,k)=t(i,k) - qo(i,k)=max(1.e-16,qv(i,k)) - enddo + do k=kts,kpbli(i) + tn(i,k)=t(i,k) + qo(i,k)=max(1.e-16,qv(i,k)) + enddo enddo nbegin=0 nend=0 - do i=its,itf - do k=kts,kpbli(i) - dhdt(i,k)=cp*(forcet(i,k)+(t(i,k)-t2di(i,k))/dt) + & - xlv*(forceqv(i,k)+(qv(i,k)-qv2di(i,k))/dt) -! tshall(i,k)=t(i,k) -! qshall(i,k)=qv(i,k) - enddo - enddo - do k= kts+1,ktf-1 - do i = its,itf - if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then - dp=-.5*(p2d(i,k+1)-p2d(i,k-1)) - umean(i)=umean(i)+us(i,k)*dp - vmean(i)=vmean(i)+vs(i,k)*dp - pmean(i)=pmean(i)+dp - endif - enddo + do i=its,itf + do k=kts,kpbli(i) + dhdt(i,k)=cp*(forcet(i,k)+(t(i,k)-t2di(i,k))/dt) + & + xlv*(forceqv(i,k)+(qv(i,k)-qv2di(i,k))/dt) +! tshall(i,k)=t(i,k) +! qshall(i,k)=qv(i,k) enddo - do k=kts,ktf-1 + enddo + do k= kts+1,ktf-1 do i = its,itf - omeg(i,k)= w(i,k) !-g*rhoi(i,k)*w(i,k) -! dq=(q2d(i,k+1)-q2d(i,k)) -! mconv(i)=mconv(i)+omeg(i,k)*dq/g - enddo + if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then + dp=-.5*(p2d(i,k+1)-p2d(i,k-1)) + umean(i)=umean(i)+us(i,k)*dp + vmean(i)=vmean(i)+vs(i,k)*dp + pmean(i)=pmean(i)+dp + endif enddo + enddo + do k=kts,ktf-1 do i = its,itf - if(mconv(i).lt.0.)mconv(i)=0. + omeg(i,k)= w(i,k) !-g*rhoi(i,k)*w(i,k) +! dq=(q2d(i,k+1)-q2d(i,k)) +! mconv(i)=mconv(i)+omeg(i,k)*dq/g enddo + enddo + do i = its,itf + if(mconv(i).lt.0.)mconv(i)=0. + enddo ! !---- call cumulus parameterization ! if(ishallow_g3.eq.1)then -! + do i=its,ite ierrs(i)=0 ierrm(i)=0 @@ -499,14 +489,13 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! !> - Call shallow: cu_gf_sh_run() ! - ! print*,'hli bf shallow t2d',t2d call cu_gf_sh_run (us,vs, & ! input variables, must be supplied zo,t2d,q2d,ter11,tshall,qshall,p2d,psur,dhdt,kpbli, & - rhoi,hfx,qfx,xlandi,ichoice_s,tcrit,dt, & + rhoi,hfx,qfx,xlandi,ichoice_s,tcrit,dt, & ! input variables. ierr should be initialized to zero or larger than zero for ! turning off shallow convection for grid points - zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs, & + zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs, & ! output tendencies outts,outqs,outqcs,outus,outvs,cnvwt,prets,cupclws, & ! dimesnional variables @@ -524,8 +513,8 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ipr=0 jpr_deep=0 !340765 !> - Call cu_gf_deep_run() for middle GF convection - if(imid_gf == 1)then - call cu_gf_deep_run( & + if(imid_gf == 1)then + call cu_gf_deep_run( & itf,ktf,its,ite, kts,kte & ,dicycle_m & ,ichoicem & @@ -594,16 +583,16 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ,jminm,tropics) do i=its,itf - do k=kts,ktf + do k=kts,ktf qcheck(i,k)=qv(i,k) +outqs(i,k)*dt - enddo + enddo enddo !> - Call neg_check() for middle GF convection call neg_check('mid',ipn,dt,qcheck,outqm,outtm,outum,outvm, & outqcm,pretm,its,ite,kts,kte,itf,ktf,ktopm) - endif + endif !> - Call cu_gf_deep_run() for deep GF convection - if(ideep.eq.1)then + if(ideep.eq.1)then call cu_gf_deep_run( & itf,ktf,its,ite, kts,kte & @@ -673,15 +662,15 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & #endif ,k22 & ,jmin,tropics) - jpr=0 - ipr=0 - do i=its,itf - do k=kts,ktf - qcheck(i,k)=qv(i,k) +(outqs(i,k)+outqm(i,k))*dt - enddo - enddo + jpr=0 + ipr=0 + do i=its,itf + do k=kts,ktf + qcheck(i,k)=qv(i,k) +(outqs(i,k)+outqm(i,k))*dt + enddo + enddo !> - Call neg_check() for deep GF convection - call neg_check('deep',ipn,dt,qcheck,outq,outt,outu,outv, & + call neg_check('deep',ipn,dt,qcheck,outq,outt,outu,outv, & outqc,pret,its,ite,kts,kte,itf,ktf,ktop) ! endif @@ -730,6 +719,11 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & enddo ! do i=its,itf + massflx(:)=0. + trcflx_in1(:)=0. + clw_in1(:)=0. + clw_ten1(:)=0. + po_cup(:)=0. kstop=kts if(ktopm(i).gt.kts .or. ktop(i).gt.kts)kstop=max(ktopm(i),ktop(i)) if(ktops(i).gt.kts)kstop=max(kstop,ktops(i)) @@ -738,7 +732,8 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & if(kbcon(i).gt.2 .or. kbconm(i).gt.2)then hbot(i)=max(kbconm(i),kbcon(i)) !jmin(i) endif -!kbcon(i) + + dtime_max=dt do k=kts,kstop cnvc(i,k) = 0.04 * log(1. + 675. * zu(i,k) * xmb(i)) + & 0.04 * log(1. + 675. * zum(i,k) * xmbm(i)) + & @@ -754,66 +749,117 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & us(i,k)=us(i,k)+outu(i,k)*cuten(i)*dt +outum(i,k)*cutenm(i)*dt +outus(i,k)*cutens(i)*dt vs(i,k)=vs(i,k)+outv(i,k)*cuten(i)*dt +outvm(i,k)*cutenm(i)*dt +outvs(i,k)*cutens(i)*dt -!hj 10/11/2016: don't need gdc and gdc2 yet for gsm. -!hli 08/18/2017: couple gdc to radiation - gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i)) ! my mod + gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i)) ! my mod gdc2(i,k,1)=max(0.,tun_rad_deep(i)*(cupclwm(i,k)*cutenm(i)+cupclw(i,k)*cuten(i))) gdc(i,k,2)=(outt(i,k))*86400. gdc(i,k,3)=(outtm(i,k))*86400. gdc(i,k,4)=(outts(i,k))*86400. gdc(i,k,7)=-(gdc(i,k,7)-sqrt(us(i,k)**2 +vs(i,k)**2))/dt - !gdc(i,k,8)=(outq(i,k))*86400.*xlv/cp + !gdc(i,k,8)=(outq(i,k))*86400.*xlv/cp gdc(i,k,8)=(outqm(i,k)+outqs(i,k)+outq(i,k))*86400.*xlv/cp gdc(i,k,9)=gdc(i,k,2)+gdc(i,k,3)+gdc(i,k,4) - if((gdc(i,k,1).ge.0.5).or.(gdc2(i,k,1).ge.0.5))then - print*,'hli gdc(i,k,1),gdc2(i,k,1)',gdc(i,k,1),gdc2(i,k,1) - endif ! !> - Calculate subsidence effect on clw ! - dsubclw=0. - dsubclwm=0. - dsubclws=0. +! dsubclw=0. +! dsubclwm=0. +! dsubclws=0. +! dp=100.*(p2d(i,k)-p2d(i,k+1)) +! if (clcw(i,k) .gt. -999.0 .and. clcw(i,k+1) .gt. -999.0 )then +! clwtot = cliw(i,k) + clcw(i,k) +! clwtot1= cliw(i,k+1) + clcw(i,k+1) +! dsubclw=((-edt(i)*zd(i,k+1)+zu(i,k+1))*clwtot1 & +! -(-edt(i)*zd(i,k) +zu(i,k)) *clwtot )*g/dp +! dsubclwm=((-edtm(i)*zdm(i,k+1)+zum(i,k+1))*clwtot1 & +! -(-edtm(i)*zdm(i,k) +zum(i,k)) *clwtot )*g/dp +! dsubclws=(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp +! dsubclw=dsubclw+(zu(i,k+1)*clwtot1-zu(i,k)*clwtot)*g/dp +! dsubclwm=dsubclwm+(zum(i,k+1)*clwtot1-zum(i,k)*clwtot)*g/dp +! dsubclws=dsubclws+(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp +! endif +! tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & +! +outqcm(i,k)*cutenm(i) & +! +dsubclw*xmb(i)+dsubclws*xmbs(i)+dsubclwm*xmbm(i) & +! ) +! tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) +! if (clcw(i,k) .gt. -999.0) then +! cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice +! clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water +! else +! cliw(i,k) = max(0.,cliw(i,k) + tem) +! endif +! +! enddo + +!> - FCT treats subsidence effect to cloud ice/water (begin) dp=100.*(p2d(i,k)-p2d(i,k+1)) + dtime_max=min(dtime_max,.5*dp) + po_cup(k)=.5*(p2d(i,k)+p2d(i,k+1)) if (clcw(i,k) .gt. -999.0 .and. clcw(i,k+1) .gt. -999.0 )then clwtot = cliw(i,k) + clcw(i,k) + if(clwtot.lt.1.e-32)clwtot=0. clwtot1= cliw(i,k+1) + clcw(i,k+1) - dsubclw=((-edt(i)*zd(i,k+1)+zu(i,k+1))*clwtot1 & - -(-edt(i)*zd(i,k) +zu(i,k)) *clwtot )*g/dp - dsubclwm=((-edtm(i)*zdm(i,k+1)+zum(i,k+1))*clwtot1 & - -(-edtm(i)*zdm(i,k) +zum(i,k)) *clwtot )*g/dp - dsubclws=(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp - dsubclw=dsubclw+(zu(i,k+1)*clwtot1-zu(i,k)*clwtot)*g/dp - dsubclwm=dsubclwm+(zum(i,k+1)*clwtot1-zum(i,k)*clwtot)*g/dp - dsubclws=dsubclws+(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp + if(clwtot1.lt.1.e-32)clwtot1=0. + clw_in1(k)=clwtot + massflx(k)=-(xmb(i) *( zu(i,k)- edt(i)* zd(i,k))) & + -(xmbm(i)*(zdm(i,k)-edtm(i)*zdm(i,k))) & + -(xmbs(i)*zus(i,k)) + trcflx_in1(k)=massflx(k)*.5*(clwtot+clwtot1) endif - tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & + enddo + + massflx (1)=0. + trcflx_in1(1)=0. + call fct1d3 (kstop,kte,dtime_max,po_cup, & + clw_in1,massflx,trcflx_in1,clw_ten1,g) + + do k=1,kstop + tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & +outqcm(i,k)*cutenm(i) & -! +dsubclw*xmb(i)+dsubclws*xmbs(i)+dsubclwm*xmbm(i) & - ) + +clw_ten1(k) & + ) tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) if (clcw(i,k) .gt. -999.0) then cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water - else + else cliw(i,k) = max(0.,cliw(i,k) + tem) - endif + endif - enddo - gdc(i,1,10)=forcing(i,1) - gdc(i,2,10)=forcing(i,2) - gdc(i,3,10)=forcing(i,3) - gdc(i,4,10)=forcing(i,4) - gdc(i,5,10)=forcing(i,5) - gdc(i,6,10)=forcing(i,6) - gdc(i,7,10)=forcing(i,7) - gdc(i,8,10)=forcing(i,8) - gdc(i,10,10)=xmb(i) - gdc(i,11,10)=xmbm(i) - gdc(i,12,10)=xmbs(i) - gdc(i,13,10)=hfx(i) - gdc(i,15,10)=qfx(i) - gdc(i,16,10)=pret(i)*3600. +! +!> calculate cloud water and cloud ice number concentrations +! + rho_dryar(i,k) = p2di(i,k)/(con_rd*t(i,k)) ! Density of dry air in kg m-3 + if (imp_physics == imp_physics_thompson) then + if ((tem*tem1)>1.e-5) then + gq0(i,k,ntinc) = max(0., gq0(i,k,ntinc) + & + make_IceNumber(tem*tem1*rho_dryar(i,k), t(i,k)) * & + (1/rho_dryar(i,k))) + end if + if ((tem*(1-tem1))>1.e-5) then + gq0(i,k,ntlnc) = max(0., gq0(i,k,ntlnc) + & + make_DropletNumber(tem*(1-tem1)*rho_dryar(i,k), nwfa(i,k)) & + * (1/rho_dryar(i,k))) + end if + end if + + enddo + + + gdc(i,1,10)=forcing(i,1) + gdc(i,2,10)=forcing(i,2) + gdc(i,3,10)=forcing(i,3) + gdc(i,4,10)=forcing(i,4) + gdc(i,5,10)=forcing(i,5) + gdc(i,6,10)=forcing(i,6) + gdc(i,7,10)=forcing(i,7) + gdc(i,8,10)=forcing(i,8) + gdc(i,10,10)=xmb(i) + gdc(i,11,10)=xmbm(i) + gdc(i,12,10)=xmbs(i) + gdc(i,13,10)=hfx(i) + gdc(i,15,10)=qfx(i) + gdc(i,16,10)=pret(i)*3600. if(ktop(i).gt.2 .and.pret(i).gt.0.)dt_mf(i,ktop(i)-1)=ud_mf(i,ktop(i)) endif enddo diff --git a/physics/cu_gf_driver.meta b/physics/cu_gf_driver.meta index 1969f9464..d3687a352 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -44,6 +44,14 @@ [ccpp-arg-table] name = cu_gf_driver_run type = scheme +[ntracer] + standard_name = number_of_tracers + long_name = number of tracers + units = count + dimensions = () + type = integer + intent = in + optional = F [garea] standard_name = cell_area long_name = grid cell area @@ -350,6 +358,65 @@ type = integer intent = in optional = F +[nwfa] + standard_name = water_friendly_aerosol_number_concentration + long_name = number concentration of water-friendly aerosols + units = kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[con_rd] + standard_name = gas_constant_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[gq0] + standard_name = tracer_concentration_updated_by_physics + long_name = tracer concentration updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension,number_of_tracers) + type = real + kind = kind_phys + intent = inout + optional = F +[ntinc] + standard_name = index_for_ice_cloud_number_concentration + long_name = tracer index for ice number concentration + units = index + dimensions = () + type = integer + intent = in + optional = F +[ntlnc] + standard_name = index_for_liquid_cloud_number_concentration + long_name = tracer index for liquid number concentration + units = index + dimensions = () + type = integer + intent = in + optional = F +[imp_physics] + standard_name = flag_for_microphysics_scheme + long_name = choice of microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imp_physics_thompson] + standard_name = flag_for_thompson_microphysics_scheme + long_name = choice of Thompson microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/module_MYNNrad_post.F90 b/physics/module_MYNNrad_post.F90 index 7acd2e406..1364db62e 100644 --- a/physics/module_MYNNrad_post.F90 +++ b/physics/module_MYNNrad_post.F90 @@ -22,6 +22,7 @@ end subroutine mynnrad_post_finalize #endif SUBROUTINE mynnrad_post_run( & & ix,im,levs, & + & flag_init,flag_restart, & & qc,qi, & & qc_save, qi_save, & & errmsg, errflg ) @@ -34,6 +35,7 @@ SUBROUTINE mynnrad_post_run( & !------------------------------------------------------------------- integer, intent(in) :: ix, im, levs + logical, intent(in) :: flag_init, flag_restart real(kind=kind_phys), dimension(im,levs), intent(out) :: qc, qi real(kind=kind_phys), dimension(im,levs), intent(in) :: qc_save, qi_save character(len=*), intent(out) :: errmsg @@ -48,6 +50,11 @@ SUBROUTINE mynnrad_post_run( & !write(0,*)"==============================================" !write(0,*)"in mynn rad post" + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_post flag_init = ', flag_init + return + endif + ! Add subgrid cloud information: do k = 1, levs do i = 1, im diff --git a/physics/module_MYNNrad_post.meta b/physics/module_MYNNrad_post.meta index b09abe01e..881a19fff 100644 --- a/physics/module_MYNNrad_post.meta +++ b/physics/module_MYNNrad_post.meta @@ -25,6 +25,22 @@ type = integer intent = in optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F [qc] standard_name = cloud_condensed_water_mixing_ratio long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) diff --git a/physics/module_MYNNrad_pre.F90 b/physics/module_MYNNrad_pre.F90 index 858abebee..95dc95445 100644 --- a/physics/module_MYNNrad_pre.F90 +++ b/physics/module_MYNNrad_pre.F90 @@ -32,6 +32,7 @@ end subroutine mynnrad_pre_finalize !###=================================================================== SUBROUTINE mynnrad_pre_run( & & ix,im,levs, & + & flag_init,flag_restart, & & qc, qi, T3D, & & qc_save, qi_save, & & qc_bl,cldfra_bl, & @@ -50,6 +51,7 @@ SUBROUTINE mynnrad_pre_run( & ! Interface variables real (kind=kind_phys), parameter :: gfac=1.0e5/con_g integer, intent(in) :: ix, im, levs + logical, intent(in) :: flag_init, flag_restart real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc, qi real(kind=kind_phys), dimension(im,levs), intent(in) :: T3D,delp real(kind=kind_phys), dimension(im,levs), intent(inout) :: & @@ -71,13 +73,17 @@ SUBROUTINE mynnrad_pre_run( & !write(0,*)"==============================================" !write(0,*)"in mynn rad pre" + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_pre flag_init = ', flag_init + return + endif ! Add subgrid cloud information: do k = 1, levs do i = 1, im qc_save(i,k) = qc(i,k) qi_save(i,k) = qi(i,k) - clouds1(i,k)=CLDFRA_BL(i,k) + clouds1(i,k) = CLDFRA_BL(i,k) IF (qc(i,k) < 1.E-6 .AND. qi(i,k) < 1.E-8 .AND. CLDFRA_BL(i,k)>0.001) THEN !Partition the BL clouds into water & ice according to a linear diff --git a/physics/module_MYNNrad_pre.meta b/physics/module_MYNNrad_pre.meta index 617ee3f31..3b5943c66 100644 --- a/physics/module_MYNNrad_pre.meta +++ b/physics/module_MYNNrad_pre.meta @@ -25,6 +25,22 @@ type = integer intent = in optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F [qc] standard_name = cloud_condensed_water_mixing_ratio long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index ea5800736..7345f2667 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -67,7 +67,7 @@ SUBROUTINE LSMRUC( & Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & GLW,GSW,EMISS,CHKLOWQ, CHS, & FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, & - Z0,SNOALB,ALBBCK, & !Z0,SNOALB,ALBBCK,LAI, & + Z0,SNOALB,ALBBCK,LAI, & landusef, nlcat, & ! mosaic_lu, mosaic_soil, & soilctop, nscat, & QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, & @@ -218,6 +218,7 @@ SUBROUTINE LSMRUC( & CANWAT, & ! new SNOALB, & ALB, & + LAI, & EMISS, & MAVAIL, & SFCEXC, & @@ -269,7 +270,6 @@ SUBROUTINE LSMRUC( & PC, & SFCRUNOFF, & UDRUNOFF, & - LAI, & EMISSL, & ZNTL, & LMAVAIL, & @@ -431,8 +431,8 @@ SUBROUTINE LSMRUC( & !! or ~100 mm of snow height ! ! snowc(i,j) = min(1.,snow(i,j)/32.) - soilt1(i,j)=soilt(i,j) - if(snow(i,j).le.32.) soilt1(i,j)=tso(i,1,j) +! soilt1(i,j)=soilt(i,j) +! if(snow(i,j).le.32.) soilt1(i,j)=tso(i,1,j) !> - Initializing inside snow temp if it is not defined IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN IF(snow(i,j).gt.32.) THEN @@ -450,7 +450,9 @@ SUBROUTINE LSMRUC( & patmb=P8w(i,kms,j)*1.e-2 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATMB IF((qvg(i,j) .LE. 0.) .or. (qvg(i,j) .GT.0.1)) THEN - qvg (i,j) = QSG(i,j)*mavail(i,j) + !17sept19 - bad approximation with very low mavail. + !qvg(i,j) = QSG(i,j)*mavail(i,j) + qvg (i,j) = qv3d(i,1,j) IF (debug_print ) THEN print *, & 'QVG is initialized in RUCLSM ', qvg(i,j),mavail(i,j),qsg(i,j),i,j @@ -751,7 +753,7 @@ SUBROUTINE LSMRUC( & meltfactor = 0.85 do k=2,nzs - if(zsmain(k).ge.1.0) then + if(zsmain(k).ge.1.1) then NROOT=K goto 111 endif diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90 index fe12b5e17..02061079e 100644 --- a/physics/sfc_drv_ruc.F90 +++ b/physics/sfc_drv_ruc.F90 @@ -139,11 +139,11 @@ end subroutine lsm_ruc_finalize ! DH* TODO - make order of arguments the same as in the metadata table subroutine lsm_ruc_run & ! inputs & ( iter, me, master, kdt, im, nlev, lsoil_ruc, lsoil, zs, & - & t1, q1, qc, soiltyp, vegtype, sigmaf, & + & t1, q1, qc, soiltyp, vegtype, sigmaf, laixy, & & sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, & & prsl1, zf, wind, shdmin, shdmax, alvwf, alnwf, & & snoalb, sfalb, flag_iter, flag_guess, isot, ivegsrc, fice, & - & smc, stc, slc, lsm_ruc, lsm, land, islimsk, & + & smc, stc, slc, lsm_ruc, lsm, land, islimsk, rdlai, & & imp_physics, imp_physics_gfdl, imp_physics_thompson, & & smcwlt2, smcref2, do_mynnsfclay, & & con_cp, con_rv, con_rd, con_g, con_pi, con_hvap, con_fvirt,& ! constants @@ -178,6 +178,8 @@ subroutine lsm_ruc_run & ! inputs & ch, prsl1, wind, shdmin, shdmax, & & snoalb, alvwf, alnwf, zf, qc, q1 + real (kind=kind_phys), dimension(:), intent(in) :: laixy + real (kind=kind_phys), intent(in) :: delt real (kind=kind_phys), intent(in) :: con_cp, con_rv, con_g, & con_pi, con_rd, & @@ -187,6 +189,8 @@ subroutine lsm_ruc_run & ! inputs integer, dimension(im), intent(in) :: islimsk ! sea/land/ice mask (=0/1/2) logical, intent(in) :: do_mynnsfclay + logical, intent(in) :: rdlai + ! --- in/out: integer, dimension(im), intent(inout) :: soiltyp, vegtype real (kind=kind_phys), dimension(lsoil_ruc) :: dzs @@ -317,6 +321,8 @@ subroutine lsm_ruc_run & ! inputs zs, sh2o, smfrkeep, tslb, smois, wetness, & ! out me, master, errmsg, errflg) + xlai = 0. + endif ! flag_init=.true.,iter=1 !-- end of initialization @@ -508,10 +514,10 @@ subroutine lsm_ruc_run & ! inputs ffrozp(i,j) = real(nint(srflag(i)),kind_phys) endif - !tgs - for now set rdlai2d to .false., WRF has LAI maps, and RUC LSM - ! uses rdlai2d = .true. - rdlai2d = .false. - !if( .not. rdlai2d) xlai = lai_data(vtype) + !tgs - rdlai is .false. when the LAI data is not available in the + ! - INPUT/sfc_data.nc + + rdlai2d = rdlai conflx2(i,1,j) = zf(i) * 2. ! factor 2. is needed to get the height of ! atm. forcing inside RUC LSM (inherited @@ -552,13 +558,15 @@ subroutine lsm_ruc_run & ! inputs !prcp(i,j) = rhoh2o * tprcp(i) ! tprcp in [m] - convective plus explicit !raincv(i,j) = rhoh2o * rainc(i) ! total time-step convective precip !rainncv(i,j) = rhoh2o * max(rain(i)-rainc(i),0.0) ! total time-step explicit precip - !graupelncv(i,j) = rhoh2o * graupel(i) - !snowncv(i,j) = rhoh2o * snow(i) - prcp(i,j) = rhoh2o * (rainc(i)+rainnc(i)) ! tprcp in [m] - convective plus explicit - raincv(i,j) = rhoh2o * rainc(i) ! total time-step convective precip - rainncv(i,j) = rhoh2o * rainnc(i) ! total time-step explicit precip + prcp(i,j) = rhoh2o * (rainc(i)+rainnc(i)) ! [mm] - convective plus explicit + raincv(i,j) = rhoh2o * rainc(i) ! [mm] - total time-step convective precip + rainncv(i,j) = rhoh2o * rainnc(i) ! [mm] - total time-step explicit precip graupelncv(i,j) = rhoh2o * graupel(i) snowncv(i,j) = rhoh2o * snow(i) + !if(prcp(i,j) > 0. .and. i==21) then + !print *,'prcp(i,j),rainncv(i,j),graupelncv(i,j),snowncv(i,j),ffrozp(i,j)',i,j, & + ! prcp(i,j),rainncv(i,j),graupelncv(i,j),snowncv(i,j),ffrozp(i,j) + !endif ! ice not used ! precipfr(i,j) = rainncv(i,j) * ffrozp(i,j) @@ -601,6 +609,8 @@ subroutine lsm_ruc_run & ! inputs albbck(i,j) = max(0.01, 0.5 * (alvwf(i) + alnwf(i))) alb(i,j) = sfalb(i) + if(rdlai2d) xlai(i,j) = laixy(i) + tbot(i,j) = tg3(i) !> - 4. history (state) variables (h): @@ -686,7 +696,7 @@ subroutine lsm_ruc_run & ! inputs znt(i,j) = zorl(i)/100. if(debug_print) then - !if(me==0 .and. i==ipr) then + if(me==0 .and. i==ipr) then write (0,*)'before RUC smsoil = ',smsoil(i,:,j), i,j write (0,*)'stsoil = ',stsoil(i,:,j), i,j write (0,*)'soilt = ',soilt(i,j), i,j @@ -780,7 +790,7 @@ subroutine lsm_ruc_run & ! inputs write (0,*)'shdmin1d(i,j) =',i,j,shdmin1d(i,j) write (0,*)'shdmax1d(i,j) =',i,j,shdmax1d(i,j) write (0,*)'rdlai2d =',rdlai2d - !endif + endif endif !> - Call RUC LSM lsmruc(). @@ -796,8 +806,7 @@ subroutine lsm_ruc_run & ! inputs & chs(i,j), flqc(i,j), flhc(i,j), & ! --- input/outputs: & wet(i,j), cmc(i,j), shdfac(i,j), alb(i,j), znt(i,j), & - & z0(i,j), snoalb1d(i,j), albbck(i,j), & -! & z0, snoalb1d, alb, xlai, & + & z0(i,j), snoalb1d(i,j), albbck(i,j), xlai(i,j), & & landusef(i,:,j), nlcat, & ! --- mosaic_lu and mosaic_soil are moved to the namelist ! & mosaic_lu, mosaic_soil, & @@ -820,6 +829,7 @@ subroutine lsm_ruc_run & ! inputs & its,ite, jts,jte, kts,kte ) if(debug_print) then + !if(me==0 .and. i==ipr) then write (0,*)'after sneqv(i,j) =',i,j,sneqv(i,j) write (0,*)'after snowh(i,j) =',i,j,snowh(i,j) write (0,*)'after sncovr(i,j) =',i,j,sncovr(i,j) diff --git a/physics/sfc_drv_ruc.meta b/physics/sfc_drv_ruc.meta index dac459405..3ae9a57a3 100644 --- a/physics/sfc_drv_ruc.meta +++ b/physics/sfc_drv_ruc.meta @@ -198,6 +198,12 @@ type = integer intent = in optional = F +[rdlai] + standard_name = flag_for_reading_leaf_area_index_from_input + long_name = flag for reading leaf area index from initial conditions for RUC LSM + units = flag + dimensions = () + type = logical [zs] standard_name = depth_of_soil_levels_for_land_surface_model long_name = depth of soil levels for land surface model @@ -529,6 +535,14 @@ kind = kind_phys intent = in optional = F +[laixy] + standard_name = leaf_area_index + long_name = leaf area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + optional = F [sfalb] standard_name = surface_diffused_shortwave_albedo long_name = mean surface diffused sw albedo From 988e95a37bd4ea3fa1b420107bcf02c3ed397bd3 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 26 Nov 2019 11:52:21 -0700 Subject: [PATCH 2/4] physics/GFS_suite_interstitial.*: use new imfdeepcnv_gf parameter instead of hard-coded number 3 --- physics/GFS_suite_interstitial.F90 | 10 +++++----- physics/GFS_suite_interstitial.meta | 24 ++++++++++++++++-------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 20f51f99c..1e8545e98 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -659,10 +659,10 @@ end subroutine GFS_suite_interstitial_4_finalize !> \section arg_table_GFS_suite_interstitial_4_run Argument Table !! \htmlinclude GFS_suite_interstitial_4_run.html !! - subroutine GFS_suite_interstitial_4_run (imfdeepcnv, im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & + 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, dqdti, errmsg, errflg) + gq0, clw, dqdti, imfdeepcnv, imfdeepcnv_gf, errmsg, errflg) use machine, only: kind_phys @@ -670,9 +670,9 @@ subroutine GFS_suite_interstitial_4_run (imfdeepcnv, im, levs, ltaerosol, cplchm ! interface variables - integer, intent(in) :: imfdeepcnv, im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, & + integer, intent(in) :: im, levs, 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 + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imfdeepcnv, imfdeepcnv_gf logical, intent(in) :: ltaerosol, cplchm @@ -737,7 +737,7 @@ subroutine GFS_suite_interstitial_4_run (imfdeepcnv, im, levs, ltaerosol, cplchm enddo enddo ! if (imp_physics == imp_physics_thompson) then - if (imp_physics == imp_physics_thompson .and. imfdeepcnv /= 3) then + if (imp_physics == imp_physics_thompson .and. imfdeepcnv /= imfdeepcnv_gf) then if (ltaerosol) then do k=1,levs do i=1,im diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index 2fa377c00..e6e349a2a 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1454,14 +1454,6 @@ [ccpp-arg-table] name = GFS_suite_interstitial_4_run type = scheme -[imfdeepcnv] - standard_name = flag_for_mass_flux_deep_convection_scheme - long_name = flag for mass-flux deep convection scheme - units = flag - dimensions = () - type = integer - intent = in - optional = F [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent @@ -1709,6 +1701,22 @@ kind = kind_phys intent = inout optional = F +[imfdeepcnv] + standard_name = flag_for_mass_flux_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imfdeepcnv_gf] + standard_name = flag_for_gf_deep_convection_scheme + long_name = flag for Grell-Freitas deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 17f585a070cf266859622d18a04cd0106384bf12 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 26 Nov 2019 11:53:00 -0700 Subject: [PATCH 3/4] physics/drag_suite.F90: bugfix, initialize rstoch to zero (since SPP is not used) --- physics/drag_suite.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index c3da28334..269bf0b3a 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -596,6 +596,7 @@ subroutine drag_suite_run( & olss(i) = 0.0 ulow (i) = 0.0 dtfac(i) = 1.0 + rstoch(i) = 0.0 ldrag(i) = .false. icrilv(i) = .false. flag(i) = .true. From 80edd7812cfa70db0c589026a0a65f2cba1b81a9 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 26 Nov 2019 11:53:27 -0700 Subject: [PATCH 4/4] physics/sfc_drv_ruc.F90: remove comment line that was left mistakenly --- physics/sfc_drv_ruc.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90 index 02061079e..3b4b8a118 100644 --- a/physics/sfc_drv_ruc.F90 +++ b/physics/sfc_drv_ruc.F90 @@ -829,7 +829,6 @@ subroutine lsm_ruc_run & ! inputs & its,ite, jts,jte, kts,kte ) if(debug_print) then - !if(me==0 .and. i==ipr) then write (0,*)'after sneqv(i,j) =',i,j,sneqv(i,j) write (0,*)'after snowh(i,j) =',i,j,snowh(i,j) write (0,*)'after sncovr(i,j) =',i,j,sncovr(i,j)