Add an option to use liquid potential temperature in temperature

equation of GFS PBL scheme. (from @WeiguoWang-NOAA)

Author: BinLiu-NOAA

physics/PBL/SATMEDMF/satmedmfvdifq.F

@@ -82,7 +82,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
      &     prsi,del,prsl,prslk,phii,phil,delt,                          &
      &     dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh,           &
      &     kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr,          &
-     &     rlmx,elmx,sfc_rlm,tc_pbl,                                    &
+     &     rlmx,elmx,sfc_rlm,tc_pbl,use_lpt,                            &
      &     ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind,       &
      &     index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d,       &
      &     errmsg,errflg)
@@ -97,6 +97,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
      &                        ntke, ntqv
       integer, intent(in)  :: sfc_rlm
       integer, intent(in)  :: tc_pbl
+      integer, intent(in)  :: use_lpt
       integer, intent(in)  :: kinver(:)
       integer, intent(out) :: kpbl(:)
       logical, intent(in)  :: gen_tend,ldiag3d
@@ -259,6 +260,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
       real(kind=kind_phys) h1
 
       real(kind=kind_phys) bfac, mffac
+
+      real(kind=kind_phys) qice(im,km),qliq(im,km)
 !!
       parameter(bfac=100.)
       parameter(wfac=7.0,cfac=4.5)
@@ -467,15 +470,20 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
         do i=1,im
           pix(i,k)   = psk(i) / prslk(i,k)
           theta(i,k) = t1(i,k) * pix(i,k)
+          qice(i,k) = 0.0
+          qliq(i,k) = 0.0
           if(ntiw > 0) then
             tem = max(q1(i,k,ntcw),qlmin)
             tem1 = max(q1(i,k,ntiw),qlmin)
+            qice(i,k) = tem1
+            qliq(i,k) = tem
             qlx(i,k) = tem + tem1
             ptem = hvap*tem + (hvap+hfus)*tem1
             slx(i,k)   = cp * t1(i,k) + phil(i,k) - ptem
           else
             qlx(i,k) = max(q1(i,k,ntcw),qlmin)
             slx(i,k)   = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k)
+            qliq(i,k) = qlx(i,k)
           endif
           tem2       = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k)
           thvx(i,k)  = theta(i,k) * tem2
@@ -1818,6 +1826,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,               &
           rdz     = rdzt(i,k)
           tem1    = dsig * dkt(i,k) * rdz
           dsdzt   = tem1 * gocp
+          if (use_lpt > 0) then
+            dsdzt = dsdzt-tem1*elocp*(qliq(i,k+1)-qliq(i,k))*rdz
+     &        -(1+0.33/2.5)*tem1*elocp*(qice(i,k+1)-qice(i,k))*rdz
+          endif
           dsdz2   = tem1 * rdz
           au(i,k) = -dtodsd*dsdz2
           al(i,k) = -dtodsu*dsdz2

physics/PBL/SATMEDMF/satmedmfvdifq.meta

@@ -619,6 +619,13 @@
   dimensions = ()
   type = integer
   intent = in
+[use_lpt]
+  standard_name = control_for_using_LPT_for_TC_applications_in_the_PBL_scheme
+  long_name = control for using LPT in TC applications in the PBL scheme
+  units = none
+  dimensions = ()
+  type = integer
+  intent = in
 [ntqv]
   standard_name = index_of_specific_humidity_in_tracer_concentration_array
   long_name = tracer index for water vapor (specific humidity)