Merge pull request #13 from dustinswales/ncar-main-PR967

NCAR/main PR NCAR#967 (Bugfix and optimization of prognostic closure for the P8 physics suit)

physics/GFS_MP_generic_pre.F90

@@ -9,14 +9,14 @@ module GFS_MP_generic_pre
 !> \section arg_table_GFS_MP_generic_pre_run Argument Table
 !! \htmlinclude GFS_MP_generic_pre_run.html
 !!
-      subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl, &
+      subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, progsigma, ntcw, nncl, &
                                         ntrac, gt0, gq0, save_t, save_q, num_dfi_radar, errmsg, errflg)
 !
       use machine,               only: kind_phys
 
       implicit none
       integer,                                intent(in) :: im, levs, ntcw, nncl, ntrac, num_dfi_radar
-      logical,                                intent(in) :: ldiag3d, qdiag3d, do_aw
+      logical,                                intent(in) :: ldiag3d, qdiag3d, do_aw, progsigma
       real(kind=kind_phys), dimension(:,:),   intent(in) :: gt0
       real(kind=kind_phys), dimension(:,:,:), intent(in) :: gq0
 
@@ -39,7 +39,7 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl,
           enddo
         enddo
       endif
-      if (ldiag3d .or. do_aw) then
+      if (ldiag3d .or. do_aw .or. progsigma) then
         if(qdiag3d) then
            do n=1,ntrac
               do k=1,levs
@@ -48,7 +48,7 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl,
                  enddo
               enddo
            enddo
-        else if(do_aw) then
+        else if(do_aw .or. progsigma) then
            ! if qdiag3d, all q are saved already
            save_q(1:im,:,1) = gq0(1:im,:,1)
            do n=ntcw,ntcw+nncl-1
@@ -59,4 +59,4 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl,
 
       end subroutine GFS_MP_generic_pre_run
 
-      end module GFS_MP_generic_pre
+      end module GFS_MP_generic_pre

physics/GFS_MP_generic_pre.meta

@@ -42,6 +42,13 @@
   dimensions = ()
   type = logical
   intent = in
+[progsigma]
+  standard_name = do_prognostic_updraft_area_fraction
+  long_name = flag for prognostic area fraction in cumulus convection
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [ntcw]
   standard_name = index_of_cloud_liquid_water_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for cloud condensate (or liquid water)

physics/progsigma_calc.f90

@@ -14,8 +14,8 @@
 !> @{ 
 
       subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
-           del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,                 &
-           delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout,       &
+           flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,    &
+           delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout,           &
            sigmab,errmsg,errflg)
 !                                                           
 !                                                                                                                                             
@@ -30,7 +30,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
       real(kind=kind_phys), intent(in)  :: prevsq(im,km), q(im,km),del(im,km),    &
            qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km),           &
            omega_u(im,km),zeta(im,km)
-      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im)
+      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im),flag_shallow
       real(kind=kind_phys), intent(in) :: sigmain(im,km)
 
 !     intent out
@@ -48,14 +48,15 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 
       real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2,   &
                           fdqb,dtdyn,dxlim,rmulacvg,tem,     &
-                          DEN,betascu,dp1,invdelt
+                          DEN,betascu,betadcu,dp1,invdelt
 
      !Parameters
       gcvalmx = 0.1
       rmulacvg=10.
       epsilon=1.E-11
       km1=km-1
-      betascu = 3.0
+      betadcu = 2.0
+      betascu = 3.6
       invdelt = 1./delt
 
      !Initialization 2D
@@ -212,12 +213,19 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 
       !Reduce area fraction before coupling back to mass-flux computation. 
       !This tuning could be addressed in updraft velocity equation instead.
-      do i= 1, im
-         if(cnvflg(i)) then
-            sigmab(i)=sigmab(i)/betascu
-         endif
-      enddo
-
+      if(flag_shallow)then
+         do i= 1, im
+            if(cnvflg(i)) then
+               sigmab(i)=sigmab(i)/betascu
+            endif
+         enddo
+      else
+         do i= 1, im
+            if(cnvflg(i)) then
+               sigmab(i)=sigmab(i)/betadcu
+            endif
+         enddo
+      endif
 
 
      end subroutine progsigma_calc

physics/samfdeepcnv.f

@@ -212,6 +212,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km),
      &     omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im)
       real(kind=kind_phys) gravinv
+      logical flag_shallow
 c  physical parameters
 !     parameter(grav=grav,asolfac=0.958)
 !     parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
@@ -2882,7 +2883,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 
 !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget
       if(progsigma)then
-         call progsigma_calc(im,km,first_time_step,restart,
+         flag_shallow = .false.
+         call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
      &        del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
      &        prevsq,q,kbcon1,ktcon,cnvflg,
      &        sigmain,sigmaout,sigmab,errmsg,errflg)

physics/samfshalcnv.f

@@ -160,7 +160,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
      &                     omegac(im),zeta(im,km),dbyo1(im,km),
      &                     sigmab(im)
       real(kind=kind_phys) gravinv,dxcrtas
-
+      logical flag_shallow
 c  physical parameters
 !     parameter(g=grav,asolfac=0.89)
 !     parameter(g=grav)
@@ -1929,7 +1929,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
 c
 c Prognostic closure
       if(progsigma)then
-         call progsigma_calc(im,km,first_time_step,restart,
+         flag_shallow = .true.
+         call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
      &        del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
      &        prevsq,q,kbcon1,ktcon,cnvflg,
      &        sigmain,sigmaout,sigmab,errmsg,errflg)