Fix virtual potential temperature in the total energy calculation for…

… hydrostatic=.T. and -DUSE_COND compiling flag on (#326)
NOAA-GFDL · Feb 29, 2024 · 5f24469 · 5f24469
1 parent 5cabcf1
commit 5f24469
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 4 deletions.
diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90
@@ -361,7 +361,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot,  ng, bdt, consv_te, fill,
 !---------------------
       if ( (consv_te > 0. .or. idiag%id_te>0)  .and. (.not.do_adiabatic_init) ) then
            call compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, npz,        &
-                                     u, v, w, delz, pt, delp, q, dp1, pe, peln, phis, &
+                                     u, v, w, delz, pt, delp, q, dp1, q_con, pe, peln, phis, &
                                      gridstruct%rsin2, gridstruct%cosa_s, &
                                      zvir, cp_air, rdgas, hlv, te_2d, ua, va, teq,        &
                                      flagstruct%moist_phys, nwat, sphum, liq_wat, rainwat,   &

diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90
@@ -900,7 +900,7 @@ end subroutine Lagrangian_to_Eulerian
 
 
  subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km,       &
-                                 u, v, w, delz, pt, delp, q, qc, pe, peln, hs, &
+                                 u, v, w, delz, pt, delp, q, qc, q_con, pe, peln, hs, &
                                  rsin2_l, cosa_s_l, &
                                  r_vir,  cp, rg, hlv, te_2d, ua, va, teq, &
                                  moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
@@ -913,7 +913,7 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km,       &
    real, intent(inout), dimension(isd:ied,jsd:jed,km):: ua, va
    real, intent(in), dimension(isd:ied,jsd:jed,km):: pt, delp
    real, intent(in), dimension(isd:ied,jsd:jed,km,*):: q
-   real, intent(in), dimension(isd:ied,jsd:jed,km):: qc
+   real, intent(in), dimension(isd:ied,jsd:jed,km):: qc, q_con
    real, intent(inout)::  u(isd:ied,  jsd:jed+1,km)
    real, intent(inout)::  v(isd:ied+1,jsd:jed,  km)
    real, intent(in)::  w(isd:,jsd:,1:)   ! vertical velocity (m/s)
@@ -939,7 +939,7 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km,       &
 !----------------------
 !  call cubed_to_latlon(u, v, ua, va, dx, dy, rdxa, rdya, km, flagstruct%c2l_ord)
 
-!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,hydrostatic,hs,pt,qc,rg,peln,te_2d, &
+!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,hydrostatic,hs,pt,qc,q_con,rg,peln,te_2d, &
 !$OMP                                  pe,delp,cp,rsin2_l,u,v,cosa_s_l,delz,moist_phys,w, &
 !$OMP                                  q,nwat,liq_wat,rainwat,ice_wat,snowwat,graupel,sphum)   &
 !$OMP                          private(phiz, tv, cvm, qd)
@@ -952,7 +952,11 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km,       &
         enddo
         do k=km,1,-1
            do i=is,ie
+#ifdef USE_COND
+                tv(i,k) = pt(i,j,k)*(1.+qc(i,j,k))*(1-q_con(i,j,k))
+#else
                 tv(i,k) = pt(i,j,k)*(1.+qc(i,j,k))
+#endif
               phiz(i,k) = phiz(i,k+1) + rg*tv(i,k)*(peln(i,k+1,j)-peln(i,k,j))
            enddo
         enddo