From 1452c3fcb2ef5ed0d78ed24a2a34f7febb20fff4 Mon Sep 17 00:00:00 2001 From: kaiyuan-cheng <74800123+kaiyuan-cheng@users.noreply.github.com> Date: Thu, 5 Oct 2023 10:46:33 -0400 Subject: [PATCH 1/2] Update fv_diagnostics.F90 Fix the unit of total energy diagnostic and clean up unused nh_total_energy(). --- tools/fv_diagnostics.F90 | 80 +--------------------------------------- 1 file changed, 1 insertion(+), 79 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index e4b346e74..40a3d7af9 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -1070,7 +1070,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) ! Total energy (only when moist_phys = .T.) idiag%id_te = register_diag_field ( trim(field), 'te', axes(1:2), Time, & - 'Total Energy', 'J/kg', missing_value=missing_value ) + 'Total Energy', 'J/m/s^2', missing_value=missing_value ) ! Total Kinetic energy id_ke = register_diag_field ( trim(field), 'ke', axes(1:2), Time, & 'Total KE', 'm^2/s^2', missing_value=missing_value ) @@ -1665,12 +1665,6 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif - if ( .not. Atm(n)%flagstruct%hydrostatic ) & - call nh_total_energy(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, & - Atm(n)%w, Atm(n)%delz, Atm(n)%pt, Atm(n)%delp, & - Atm(n)%q, Atm(n)%phis, Atm(n)%gridstruct%area, Atm(n)%domain, & - sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, Atm(n)%flagstruct%nwat, & - Atm(n)%ua, Atm(n)%va, Atm(n)%flagstruct%moist_phys, a2) #endif call prt_mxm('UA_Top (m/s): ', Atm(n)%ua(isc:iec,jsc:jec,1), & isc, iec, jsc, jec, 0, 1, 1., Atm(n)%gridstruct%area_64, Atm(n)%domain) @@ -5847,78 +5841,6 @@ end subroutine eqv_pot #endif - subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, & - w, delz, pt, delp, q, hs, area, domain, & - sphum, liq_wat, rainwat, ice_wat, & - snowwat, graupel, nwat, ua, va, moist_phys, te) -!------------------------------------------------------ -! Compute vertically integrated total energy per column -!------------------------------------------------------ -! !INPUT PARAMETERS: - integer, intent(in):: km, is, ie, js, je, isd, ied, jsd, jed - integer, intent(in):: nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel - real, intent(in), dimension(isd:ied,jsd:jed,km):: ua, va, pt, delp, w - real, intent(in), dimension(is:ie,js:je,km) :: delz - real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q - real, intent(in):: hs(isd:ied,jsd:jed) ! surface geopotential - real, intent(in):: area(isd:ied, jsd:jed) - logical, intent(in):: moist_phys - type(domain2d), intent(INOUT) :: domain - real, intent(out):: te(is:ie,js:je) ! vertically integrated TE -! Local - real(kind=R_Grid) :: area_l(isd:ied, jsd:jed) - real, parameter:: cv_vap = cp_vapor - rvgas ! 1384.5 - real phiz(is:ie,km+1) - real, dimension(is:ie):: cvm, qc - real cv_air, psm - integer i, j, k - - area_l = area - cv_air = cp_air - rdgas - -!$OMP parallel do default(none) shared(te,nwat,is,ie,js,je,isd,ied,jsd,jed,km,ua,va, & -!$OMP w,q,pt,delp,delz,hs,cv_air,moist_phys,sphum,liq_wat,rainwat,ice_wat,snowwat,graupel) & -!$OMP private(phiz,cvm, qc) - do j=js,je - - do i=is,ie - te(i,j) = 0. - phiz(i,km+1) = hs(i,j) - enddo - - do i=is,ie - do k=km,1,-1 - phiz(i,k) = phiz(i,k+1) - grav*delz(i,j,k) - enddo - enddo - - if ( moist_phys ) then - do k=1,km - call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, qc, cvm) - do i=is,ie - te(i,j) = te(i,j) + delp(i,j,k)*( cvm(i)*pt(i,j,k) + hlv*q(i,j,k,sphum) + & - 0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) ) - enddo - enddo - else - do k=1,km - do i=is,ie - te(i,j) = te(i,j) + delp(i,j,k)*( cv_air*pt(i,j,k) + & - 0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) ) - enddo - enddo - endif -! Unit: kg*(m/s)^2/m^2 = Joule/m^2 - do i=is,ie - te(i,j) = te(i,j)/grav - enddo - enddo - - psm = g_sum(domain, te, is, ie, js, je, 3, area_l, 1) - if( master ) write(*,*) 'Total_Energy (J/m**2 * E9) = ', psm * 1.E-9 - - end subroutine nh_total_energy subroutine compute_brn(ua, va, delp, delz, cape, bd, npz, Time) From 17961a9b2806598e483b95f9e91e1ed9bc16cef5 Mon Sep 17 00:00:00 2001 From: JosephMouallem <75327266+JosephMouallem@users.noreply.github.com> Date: Tue, 10 Oct 2023 08:02:26 -0400 Subject: [PATCH 2/2] fix gz dimension to remap v (#292) --- model/fv_mapz.F90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index bdaac733c..8952d23f8 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -155,7 +155,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis real, dimension(isd:ied,jsd:jed,km):: pe4 real, dimension(is:ie+1,km+1):: pe0, pe3 - real, dimension(is:ie):: gz, cvm + real, dimension(is:ie+1):: gz, cvm real, dimension(isd:ied,jsd:jed,km):: qnl, qni real rcp, rg, rrg, bkh, dtmp, k1k, dlnp, tpe @@ -220,7 +220,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & do k=1,km #ifdef MOIST_CAPPA call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + ice_wat, snowwat, graupel, q, gz(is:ie), cvm(is:ie)) do i=is,ie q_con(i,j,k) = gz(i) cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) @@ -272,7 +272,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & enddo #ifdef MOIST_CAPPA call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + ice_wat, snowwat, graupel, q, gz(is:ie), cvm(is:ie)) do i=is,ie q_con(i,j,k) = gz(i) cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) @@ -366,7 +366,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & km, pe2, te, & is, ie, j, isd, ied, jsd, jed, akap, T_VAR=1, conserv=.true.) else - call map_scalar(km, peln(is,1,j), te, gz, & + call map_scalar(km, peln(is,1,j), te, gz(is:ie), & km, pn2, te, & is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), cp_air*t_min) endif @@ -374,13 +374,13 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & else if ( kord_tm<0 ) then ! Map t using logp - call map_scalar(km, peln(is,1,j), pt, gz, & + call map_scalar(km, peln(is,1,j), pt, gz(is:ie), & km, pn2, pt, & is, ie, j, isd, ied, jsd, jed, & 1, abs(kord_tm), t_min) else ! Map pt using pe - call map1_ppm (km, pe1, pt, gz, & + call map1_ppm (km, pe1, pt, gz(is:ie), & km, pe2, pt, & is, ie, j, isd, ied, jsd, jed, & 1, abs(kord_tm)) @@ -423,7 +423,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & -2, abs(kord_wz)) endif ! Remap delz for hybrid sigma-p coordinate - call map1_ppm (km, pe1, delz, gz, & ! works + call map1_ppm (km, pe1, delz, gz(is:ie), & ! works km, pe2, delz, & is, ie, j, is, ie, js, je, & 1, abs(kord_tm)) @@ -527,8 +527,8 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & ! Note: pt at this stage is T_v or T_m , unless kord_tm > 0 do k=1,km #ifdef MOIST_CAPPA - call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + call moist_cv(is,ie+1,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & + ice_wat, snowwat, graupel, q, gz(is:ie+1), cvm(is:ie+1)) if ( kord_tm < 0 ) then do i=is,ie q_con(i,j,k) = gz(i) @@ -613,7 +613,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & enddo enddo - call map1_ppm( km, pe0(is:ie,:), u, gz, & + call map1_ppm( km, pe0(is:ie,:), u, gz(is:ie), & km, pe3(is:ie,:), u, & is, ie, j, isd, ied, jsd, jed+1, & -1, kord_mt) @@ -634,7 +634,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & enddo enddo - call map1_ppm (km, pe0, v, gz, & + call map1_ppm (km, pe0, v, gz(is:ie+1), & km, pe3, v, is, ie+1, & j, isd, ied+1, jsd, jed, -1, kord_mt) @@ -660,7 +660,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & do k=km,1,-1 #ifdef MOIST_CAPPA call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + ice_wat, snowwat, graupel, q, gz(is:ie), cvm(is:ie)) do i=is,ie q_con(i,j,k) = gz(i) cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) @@ -759,7 +759,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & do k=1,km #ifdef USE_COND call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + ice_wat, snowwat, graupel, q, gz(is:ie), cvm(is:ie)) do i=is,ie ! KE using 3D winds: q_con(i,j,k) = gz(i) @@ -870,7 +870,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & else #ifdef USE_COND call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, & - ice_wat, snowwat, graupel, q, gz, cvm) + ice_wat, snowwat, graupel, q, gz(is:ie), cvm(is:ie)) do i=is,ie pt(i,j,k) = (pt(i,j,k)+dtmp/cvm(i)*pkz(i,j,k))/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i))) enddo