dtc/develop: combined PR with bugfixes, new and updated HAFS physics, documentation

physics/cu_ntiedtke.F90

@@ -148,7 +148,7 @@ end subroutine cu_ntiedtke_finalize
 !-----------------------------------------------------------------------
 !          level 1 subroutine 'tiecnvn'
 !-----------------------------------------------------------------
-      subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
+      subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
                                  evap,hfx,zprecc,lmask,lq,ix,km,dt,dx,kbot,ktop,kcnv,&
                                  ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg)
 !-----------------------------------------------------------------
@@ -162,13 +162,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
       integer, dimension( lq ),   intent(in)  :: lmask
       real(kind=kind_phys), dimension( lq ),     intent(in ) :: evap, hfx, dx
       real(kind=kind_phys), dimension( ix , km ),     intent(inout) :: pu, pv, pt, pqv
-      real(kind=kind_phys), dimension( ix , km ),     intent(in )   :: poz, prsl, pomg, pqvf, ptf
+      real(kind=kind_phys), dimension( ix , km ),     intent(in )   :: tdi, qvdi, poz, prsl, pomg, pqvf, ptf
       real(kind=kind_phys), dimension( ix , km+1 ),   intent(in )   :: pzz, prsi
-      ! DH* TODO - check dimensions of clw, ktrac+2 seems to be smaller
-      ! than the actual dimensions (ok as long as only indices 1 and 2
-      ! are accessed here, and as long as these contain what is expected);
-      ! better to expand into the cloud-ice and cloud-water components *DH
-      real(kind=kind_phys), dimension( ix , km, ktrac+2 ),    intent(inout ) ::  clw
+      real(kind=kind_phys), dimension( ix , km, ktrac ),    intent(inout ) ::  clw
 
       integer, dimension( lq ),   intent(out)  :: kbot, ktop, kcnv
       real(kind=kind_phys), dimension( lq ),   intent(out)  :: zprecc
@@ -188,13 +184,13 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
       real(kind=kind_phys) ztp1(lq,km),    zqp1(lq,km),  ztu(lq,km),   zqu(lq,km),&
      &     zlu(lq,km),     zlude(lq,km), zmfu(lq,km),  zmfd(lq,km),  zmfude_rate(lq,km),&
      &     zqsat(lq,km),   zrain(lq)
-      real(kind=kind_phys) pcen(lq,km,ktrac),ptenc(lq,km,ktrac) 
+      real(kind=kind_phys),allocatable ::  pcen(:,:,:),ptenc(:,:,:) 
 
       integer icbot(lq),   ictop(lq),     ktype(lq),   lndj(lq)
       logical locum(lq)
 !
       real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt
-      integer i,j,k,k1,n,km1
+      integer i,j,k,k1,n,km1,ktracer
       real(kind=kind_phys) ztpp1
       real(kind=kind_phys) zew,zqs,zcor
 !
@@ -246,9 +242,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           zqs  = min(0.5,zqs)                                         
           zcor = 1./(1.-vtmpc1*zqs)                                    
           zqsat(j,k1)=zqs*zcor      
-          pqte(j,k1)=pqvf(j,k) 
+          pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst 
           zqq(j,k1) =pqte(j,k1)
-          ptte(j,k1)=ptf(j,k)
+          ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst
           ztt(j,k1) =ptte(j,k1)
           ud_mf(j,k1)=0.
           dd_mf(j,k1)=0.
@@ -258,16 +254,33 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
         end do
       end do
 
-      do n=1,ktrac
-        do k=1,km
-          k1=km-k+1
-          do j=1,lq
-            pcen(j,k1,n) = clw(j,k,n+2)
-            ptenc(j,k1,n)= 0.
+      if(ktrac > 2) then
+        ktracer = ktrac - 2
+        allocate(pcen(lq,km,ktracer))
+        allocate(ptenc(lq,km,ktracer))
+        do n=1,ktracer
+          do k=1,km
+            k1=km-k+1
+            do j=1,lq
+              pcen(j,k1,n) = clw(j,k,n+2)
+              ptenc(j,k1,n)= 0.
+            end do
           end do
         end do
-      end do
-
+      else
+        ktracer = 2
+        allocate(pcen(lq,km,ktracer))
+        allocate(ptenc(lq,km,ktracer))
+        do n=1,ktracer
+          do k=1,km
+            do j=1,lq
+              pcen(j,k,n) = 0.
+              ptenc(j,k,n)= 0.
+            end do
+          end do
+        end do
+      end if
+
 !      print *, "pgeo=",pgeo(1,:)
 !      print *, "pgeoh=",pgeoh(1,:)
 !      print *, "pap=",pap(1,:)
@@ -289,7 +302,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
      &     zqp1,     pum1,     pvm1,     pverv,   zqsat,&
      &     pqhfl,    ztmst,    pap,      paph,    pgeo, &
      &     ptte,     pqte,     pvom,     pvol,    prsfc,&
-     &     pssfc,    locum,    ktrac,    pcen,    ptenc,&
+     &     pssfc,    locum,    ktracer,  pcen,    ptenc,&
      &     ktype,    icbot,    ictop,    ztu,     zqu,  &
      &     zlu,      zlude,    zmfu,     zmfd,    zrain,&
      &     pcte,     phhfl,    lndj,     pgeoh,   zmfude_rate, dx)
@@ -314,7 +327,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst
           pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst
           ud_mf(j,k)= zmfu(j,k1)*ztmst
-          dd_mf(j,k)= zmfd(j,k1)*ztmst
+          dd_mf(j,k)= -zmfd(j,k1)*ztmst
           dt_mf(j,k)= zmfude_rate(j,k1)*ztmst
           cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1))
           cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k))
@@ -343,17 +356,21 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg,   &
           end do
         end do
       endif
+
 !
-      if (ktrac > 0) then
-        do n=1,ktrac
-          do k=1,km
-            k1=km-k+1
-            do j=1,lq
-              clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
-            end do
-          end do
-        end do
-      end if
+! Currently, vertical mixing of tracers are turned off
+!      if(ktrac > 2) then
+!        do n=1,ktrac-2
+!          do k=1,km
+!            k1=km-k+1
+!            do j=1,lq
+!              clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
+!            end do
+!          end do
+!        end do
+!      end if
+      deallocate(pcen)
+      deallocate(ptenc)
 !
       return
       end subroutine cu_ntiedtke_run

physics/cu_ntiedtke.meta

@@ -80,6 +80,24 @@
   kind = kind_phys
   intent = inout
   optional = F
+[tdi]
+  standard_name = air_temperature
+  long_name = mid-layer temperature
+  units = K
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[qvdi]
+  standard_name = water_vapor_specific_humidity
+  long_name = water vapor specific humidity
+  units = kg kg-1
+  dimensions = (horizontal_dimension,vertical_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
 [pqvf]
   standard_name = moisture_tendency_due_to_dynamics
   long_name = moisture tendency due to dynamics only
@@ -254,8 +272,8 @@
   intent = out
   optional = F
 [ktrac]
-  standard_name = number_of_total_tracers
-  long_name = number of total tracers
+  standard_name = number_of_tracers_for_convective_transport
+  long_name = number of tracers for convective transport
   units = count
   dimensions = ()
   type = integer

physics/docs/ccpp_doxyfile

@@ -113,6 +113,7 @@ INPUT                  =  pdftxt/mainpage.txt     \
                           pdftxt/GFS_SFCSICE.txt  \
                           pdftxt/GFS_HEDMF.txt    \
                           pdftxt/GFS_SATMEDMF.txt    \
+                          pdftxt/GFS_SATMEDMFVDIFQ.txt \
                           pdftxt/GFS_GWDPS.txt    \
                           pdftxt/GFS_OZPHYS.txt   \
                           pdftxt/GFS_H2OPHYS.txt  \
@@ -189,6 +190,11 @@ INPUT                  =  pdftxt/mainpage.txt     \
                         ../mfpblt.f         \
                         ../mfscu.f          \
                         ../tridi.f          \
+### satmedmfvdifq
+                        ../satmedmfvdifq.F \
+                        ../mfpbltq.f         \
+                        ../mfscuq.f          \
+                        ../tridi.f          \
 ### Orographic Gravity Wave
                         ../gwdps.f \
 ### Rayleigh Dampling

physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt

@@ -0,0 +1,35 @@
+/**
+\page GFS_SATMEDMFVDIFQ GFS Scale-aware TKE-based Moist Eddy-Diffusion Mass-Flux (EDMF) PBL and Free Atmospheric Turbulence Scheme
+\section des_satmedmfvdifq Description
+
+The current operational \ref GFS_HEDMF uses a hybrid EDMF parameterization for the convective PBL (Han et al. 2016 \cite Han_2016;
+Han et al. 2017 \cite han_et_al_2017), where the EDMF scheme is applied only for the strongly unstable PBL, while the eddy-diffusivity
+counter-gradient(EDCG) scheme is used for the weakly unstable PBL. The new TKE-EDMF is an extended version of \ref GFS_HEDMF  with below enhancement:
+
+-# Eddy diffusivity (K) is now a function of TKE which is prognostically predicted
+
+-# EDMF approach is applied for all the unstable PBL 
+
+-# EDMF approach is also applied to the stratocumulus-top-driven turbulence mixing
+
+-# It includes a moist-adiabatic process when updraft thermal becomes saturated
+
+-# Scale-aware capability
+
+-# It includes interaction between TKE and cumulus convection
+
+The CCPP-compliant subroutine satmedmfvdifq_run() computes subgrid vertical turbulence mixing using scale-aware
+TKE-based moist eddy-diffusion mass-flux paramterization (Han et al. 2019 \cite Han_2019)
+- For the convective boundary layer, the scheme adopts EDMF parameterization (Siebesma et al. (2007)\cite Siebesma_2007)
+to take into account nonlocal transport by large eddies(mfpbltq.f)
+- A new mass-flux paramterization for stratocumulus-top-induced turbulence mixing has been introduced (mfscuq.f; previously,
+it was an eddy diffusion form)
+- For local turbulence mixing, a TKE closure model is used.
+
+\section intra_satmedmfvdifq Intraphysics Communication
+\ref arg_table_satmedmfvdifq_run
+
+\section gen_pbl_satmedmfvdifq  General Algorithm
+\ref gen_satmedmfvdifq
+
+*/

physics/docs/pdftxt/all_shemes_list.txt

@@ -14,6 +14,7 @@ parameterizations in suites.
 - \b PBL \b and \b Turbulence
  - \subpage GFS_HEDMF
  - \subpage GFS_SATMEDMF
+ - \subpage GFS_SATMEDMFVDIFQ
  - \subpage GSD_MYNNEDMF
 
 - \b Land \b Surface \b Model

physics/mfpbltq.f

@@ -3,7 +3,7 @@
 !! updraft parcel properties for thermals driven by surface heating 
 !! for use in the TKE-EDMF PBL scheme (updated version).
 
-!>\ingroup satmedmfq
+!>\ingroup satmedmfvdifq
 !! This subroutine computes mass flux and updraft parcel properties for
 !! thermals driven by surface heating. 
 !!\section mfpbltq_gen GFS mfpblt General Algorithm 

physics/mfscuq.f

@@ -2,7 +2,7 @@
 !! This file contains the mass flux and downdraft parcel preperties
 !! parameterization for stratocumulus-top-driven turbulence (updated version).
 
-!>\ingroup satmedmfq
+!>\ingroup satmedmfvdifq
 !! This subroutine computes mass flux and downdraft parcel properties
 !! for stratocumulus-top-driven turbulence.
 !! \section mfscuq GFS mfscu General Algorithm

physics/module_mp_thompson.F90

@@ -924,30 +924,20 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa,  &
 
       call cpu_time(stime)
 
-!$OMP parallel num_threads(threads)
-
-!$OMP sections
-
-!$OMP section
 !>  - Call qr_acr_qg() to create rain collecting graupel & graupel collecting rain table
       if (mpirank==mpiroot) write(0,*) '  creating rain collecting graupel table'
       call cpu_time(stime)
       call qr_acr_qg
       call cpu_time(etime)
       if (mpirank==mpiroot) print '("Computing rain collecting graupel table took ",f10.3," seconds.")', etime-stime
 
-!$OMP section
 !>  - Call qr_acr_qs() to create rain collecting snow & snow collecting rain table
       if (mpirank==mpiroot) write (*,*) '  creating rain collecting snow table'
       call cpu_time(stime)
       call qr_acr_qs
       call cpu_time(etime)
       if (mpirank==mpiroot) print '("Computing rain collecting snow table took ",f10.3," seconds.")', etime-stime
 
-!$OMP end sections
-
-!$OMP end parallel
-
 !>  - Call freezeh2o() to create cloud water and rain freezing (Bigg, 1953) table
       if (mpirank==mpiroot) write(0,*) '  creating freezing of water drops table'
       call cpu_time(stime)