From fc7185833989f1a817114833cd9e89e8d34f63aa Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Tue, 2 Apr 2024 14:05:21 -0400 Subject: [PATCH 1/6] fv_diag added height fluxes, corrected logic for integrated fluxes --- tools/fv_diagnostics.F90 | 190 ++++++++++++++++++++++++++++++++++----- tools/fv_diagnostics.h | 6 +- 2 files changed, 170 insertions(+), 26 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 3697ff1ff..a104f2a5e 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -872,7 +872,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) if (id_t_dt_phys_plev_ave > 0 .and. .not. allocated(Atm(n)%phys_diag%phys_t_dt) ) allocate(Atm(n)%phys_diag%phys_t_dt(isc:iec,jsc:jec,npz)) ! flag for calculation of geopotential - if ( any(id_h > 0) .or. id_h_plev>0 .or. id_hght3d>0) then + if ( any(id_h > 0) .or. id_h_plev>0 .or. id_hght3d>0 .or. id_inthght .or. id_uz .or. id_iuz .or. id_vz .or. id_ivz ) then id_any_hght = 1 else id_any_hght = 0 @@ -960,10 +960,6 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) ! 'Relative Humidity', '%', missing_value=missing_value, range=rhrange ) id_delp = register_diag_field ( trim(field), 'delp', axes(1:3), Time, & 'pressure thickness'//massdef_str, 'pa', missing_value=missing_value ) -#ifdef GFS_PHYS - id_delp_total = register_diag_field ( trim(field), 'delp_total', axes(1:3), Time, & - 'FV3 pressure thickness (dry air + all water species)', 'pa', missing_value=missing_value ) -#endif if ( .not. Atm(n)%flagstruct%hydrostatic ) & id_delz = register_diag_field ( trim(field), 'delz', axes(1:3), Time, & 'height thickness', 'm', missing_value=missing_value ) @@ -1017,6 +1013,11 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) id_vt = register_diag_field ( trim(field), 'vt', axes(1:3), Time, & 'meridional heat flux', 'K*m/sec', missing_value=missing_value ) + id_uz = register_diag_field ( trim(field), 'uz', axes(1:3), Time, & + 'zonal height flux', 'K*m/sec', missing_value=missing_value ) + id_vz = register_diag_field ( trim(field), 'vz', axes(1:3), Time, & + 'meridional height flux', 'K*m/sec', missing_value=missing_value ) + id_uu = register_diag_field ( trim(field), 'uu', axes(1:3), Time, & 'zonal flux of zonal wind', '(m/sec)^2', missing_value=missing_value ) id_uv = register_diag_field ( trim(field), 'uv', axes(1:3), Time, & @@ -1050,6 +1051,11 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) id_ivt = register_diag_field ( trim(field), 'vt_vi', axes(1:2), Time, & 'vertical integral of vt', 'K*m/sec*Pa', missing_value=missing_value ) + id_iuz = register_diag_field ( trim(field), 'uz_vi', axes(1:2), Time, & + 'vertical integral of uz', 'm*m/sec*Pa', missing_value=missing_value ) + id_ivz = register_diag_field ( trim(field), 'vz_vi', axes(1:2), Time, & + 'vertical integral of vz', 'm*m/sec*Pa', missing_value=missing_value ) + id_iuu = register_diag_field ( trim(field), 'uu_vi', axes(1:2), Time, & 'vertical integral of uu', '(m/sec)^2*Pa', missing_value=missing_value ) id_iuv = register_diag_field ( trim(field), 'uv_vi', axes(1:2), Time, & @@ -1074,7 +1080,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/m/s^2', missing_value=missing_value ) + 'Total Energy', 'J/kg', 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 ) @@ -1140,7 +1146,11 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) 'Vertically Integrated Snow', 'kg/m**2', missing_value=missing_value ) id_intqg = register_diag_field ( trim(field), 'intqg', axes(1:2), Time, & 'Vertically Integrated Graupel', 'kg/m**2', missing_value=missing_value ) - + ! Tim M. mod + id_inthght = register_diag_field ( trim(field), 'inthght', axes(1:2), Time, & + 'Vertically Integrated Height', 'kg/m', missing_value=missing_value ) + id_inttemp = register_diag_field ( trim(field), 'inttemp', axes(1:2), Time, & + 'Vertically Integrated Temperautre', 'K kg/m**2', missing_value=missing_value ) id_acl = register_diag_field ( trim(field), 'acl', axes(1:2), Time, & 'Column-averaged Cl mixing ratio', 'kg/kg', missing_value=missing_value ) id_acl2 = register_diag_field ( trim(field), 'acl2', axes(1:2), Time, & @@ -1669,6 +1679,12 @@ 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) @@ -2263,7 +2279,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) - if( id_slp>0 .or. id_tm>0 .or. id_any_hght>0 .or. id_hght3d>0 .or. id_c15>0 .or. id_ctz>0 ) then + if( id_slp>0 .or. id_tm>0 .or. id_any_hght>0 .or. id_c15>0 .or. id_ctz>0) then allocate ( wz(isc:iec,jsc:jec,npz+1) ) call get_height_field(isc, iec, jsc, jec, ngc, npz, Atm(n)%flagstruct%hydrostatic, Atm(n)%delz, & @@ -2274,6 +2290,49 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) if (id_hght3d > 0) then used = send_data(id_hght3d, 0.5*(wz(isc:iec,jsc:jec,1:npz)+wz(isc:iec,jsc:jec,2:npz+1)), Time) endif + ! Tim M. addition of integrated hght for MSE + if ( id_inthght>0 ) then + a2 = 0. + do k=1,npz + do j=jsc,jec + do i=isc,iec + a2(i,j) = a2(i,j) + 0.5*(wz(i,j,k)+wz(i,j,k+1))*Atm(n)%delp(i,j,k) + enddo + enddo + enddo + used = send_data(id_inthght, a2*ginv, Time) + endif + allocate ( a4(isc:iec,jsc:jec,npz) ) + ! zonal height flux + if(id_uz > 0 .or. id_iuz > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%ua(i,j,k) * 0.5*(wz(i,j,k)+wz(i,j,k+1)) + enddo + enddo + enddo + if(id_uz > 0) used=send_data(id_uz, a4, Time) + if(id_iuz > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(id_iuz, a2, Time) + endif + endif + ! meridional height flux + if(id_vz > 0 .or. id_ivz > 0) then + do k=1,npz + do j=jsc,jec + do i=isc,iec + a4(i,j,k) = Atm(n)%va(i,j,k) * 0.5*(wz(i,j,k)+wz(i,j,k+1)) + enddo + enddo + enddo + if(id_vz > 0) used=send_data(id_vz, a4, Time) + if(id_ivz > 0) then + call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) + used=send_data(id_ivz, a2, Time) + endif + endif if(id_slp > 0) then ! Cumpute SLP (pressure at height=0) @@ -2845,6 +2904,19 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) used = send_data(id_intqg, a2*ginv, Time) endif + ! Tim M. addition of integrated temperature for MSE + if ( id_inttemp>0 ) then + a2 = 0. + do k=1,npz + do j=jsc,jec + do i=isc,iec + a2(i,j) = a2(i,j) + Atm(n)%pt(i,j,k)*Atm(n)%delp(i,j,k) + enddo + enddo + enddo + used = send_data(id_inttemp, a2*ginv, Time) + endif + ! Cloud top temperature & cloud top press: if ( (id_ctt>0 .or. id_ctp>0 .or. id_ctz>0).and. Atm(n)%flagstruct%nwat==6) then allocate ( var1(isc:iec,jsc:jec) ) @@ -3016,6 +3088,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif endif + #ifdef GFS_PHYS if(id_delp > 0 .or. id_cape > 0 .or. id_cin > 0 .or. & ((.not. Atm(n)%flagstruct%hydrostatic) .and. (id_pfnh > 0 .or. id_ppnh > 0)) .or. & @@ -3029,7 +3102,6 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo if (id_delp > 0) used=send_data(id_delp, wk, Time) endif - if(id_delp_total > 0) used=send_data(id_delp_total, Atm(n)%delp(isc:iec,jsc:jec,:), Time) #else if(id_delp > 0) used=send_data(id_delp, Atm(n)%delp(isc:iec,jsc:jec,:), Time) #endif @@ -3846,7 +3918,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) do itrac=1, Atm(n)%ncnst call get_tracer_names (MODEL_ATMOS, itrac, tname) - if (itrac.gt.nq) then + if (id_tracer(itrac) > 0 .and. itrac.gt.nq) then used = send_data (id_tracer(itrac), Atm(n)%qdiag(isc:iec,jsc:jec,:,itrac), Time ) else used = send_data (id_tracer(itrac), Atm(n)%q(isc:iec,jsc:jec,:,itrac), Time ) @@ -3914,10 +3986,10 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) !---------------------------------- ! compute 3D flux terms !---------------------------------- - allocate ( a4(isc:iec,jsc:jec,npz) ) + if (.not. allocated(a4)) allocate ( a4(isc:iec,jsc:jec,npz) ) ! zonal moisture flux - if(id_uq > 0) then + if(id_uq > 0 .or. id_iuq > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -3925,14 +3997,14 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_uq, a4, Time) + if(id_uq > 0) used=send_data(id_uq, a4, Time) if(id_iuq > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_iuq, a2, Time) endif endif ! meridional moisture flux - if(id_vq > 0) then + if(id_vq > 0 .or. id_ivq > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -3940,7 +4012,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_vq, a4, Time) + if(id_vq > 0) used=send_data(id_vq, a4, Time) if(id_ivq > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_ivq, a2, Time) @@ -3948,7 +4020,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif ! zonal heat flux - if(id_ut > 0) then + if(id_ut > 0 .or. id_iut > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -3956,14 +4028,14 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_ut, a4, Time) + if(id_ut > 0) used=send_data(id_ut, a4, Time) if(id_iut > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_iut, a2, Time) endif endif ! meridional heat flux - if(id_vt > 0) then + if(id_vt > 0 .or. id_ivt > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -3971,7 +4043,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_vt, a4, Time) + if(id_vt > 0) used=send_data(id_vt, a4, Time) if(id_ivt > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_ivt, a2, Time) @@ -3979,7 +4051,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) endif ! zonal flux of u - if(id_uu > 0) then + if(id_uu > 0 .or. id_iuu > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -3987,14 +4059,14 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_uu, a4, Time) + if(id_uu > 0) used=send_data(id_uu, a4, Time) if(id_iuu > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_iuu, a2, Time) endif endif ! zonal flux of v - if(id_uv > 0) then + if(id_uv > 0 .or. id_iuv > 0) then do k=1,npz do j=jsc,jec do i=isc,iec @@ -4002,7 +4074,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo enddo enddo - used=send_data(id_uv, a4, Time) + if(id_uv > 0) used=send_data(id_uv, a4, Time) if(id_iuv > 0) then call z_sum(isc, iec, jsc, jec, npz, 0, Atm(n)%delp(isc:iec,jsc:jec,1:npz), a4, a2) used=send_data(id_iuv, a2, Time) @@ -5845,6 +5917,78 @@ 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) diff --git a/tools/fv_diagnostics.h b/tools/fv_diagnostics.h index d745cd2e6..108b8ee92 100644 --- a/tools/fv_diagnostics.h +++ b/tools/fv_diagnostics.h @@ -95,17 +95,17 @@ integer :: id_qr_dt_phys, id_qg_dt_phys, id_qs_dt_phys integer :: id_liq_wat_dt_phys, id_ice_wat_dt_phys integer :: id_intqv, id_intql, id_intqi, id_intqr, id_intqs, id_intqg + ! Tim M. mod + integer :: id_inthght, id_inttemp ! ESM/CM 3-D diagostics integer :: id_uq, id_vq, id_wq, id_iuq, id_ivq, id_iwq, & ! moisture flux & vertical integral id_ut, id_vt, id_wt, id_iut, id_ivt, id_iwt, & ! heat flux + id_uz, id_vz, id_iuz, id_ivz, & ! height flux id_uu, id_uv, id_vv, id_ww, & ! momentum flux id_iuu, id_iuv, id_iuw, id_ivv, id_ivw, id_iww ! vertically integral of momentum flux integer :: id_uw, id_vw integer :: id_lagrangian_tendency_of_hydrostatic_pressure integer :: id_t_dt_nudge, id_ps_dt_nudge, id_delp_dt_nudge, id_u_dt_nudge, id_v_dt_nudge -#ifdef GFS_PHYS - integer :: id_delp_total -#endif #endif _FV_DIAG__ From c9bf9d8623e7412827b312a5f4ce0167d9f89c3b Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Tue, 2 Apr 2024 14:06:33 -0400 Subject: [PATCH 2/6] fv_diagnostics removed comments --- tools/fv_diagnostics.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index a104f2a5e..fe66b1464 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -1146,7 +1146,6 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) 'Vertically Integrated Snow', 'kg/m**2', missing_value=missing_value ) id_intqg = register_diag_field ( trim(field), 'intqg', axes(1:2), Time, & 'Vertically Integrated Graupel', 'kg/m**2', missing_value=missing_value ) - ! Tim M. mod id_inthght = register_diag_field ( trim(field), 'inthght', axes(1:2), Time, & 'Vertically Integrated Height', 'kg/m', missing_value=missing_value ) id_inttemp = register_diag_field ( trim(field), 'inttemp', axes(1:2), Time, & @@ -2290,7 +2289,6 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) if (id_hght3d > 0) then used = send_data(id_hght3d, 0.5*(wz(isc:iec,jsc:jec,1:npz)+wz(isc:iec,jsc:jec,2:npz+1)), Time) endif - ! Tim M. addition of integrated hght for MSE if ( id_inthght>0 ) then a2 = 0. do k=1,npz @@ -2904,7 +2902,6 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) used = send_data(id_intqg, a2*ginv, Time) endif - ! Tim M. addition of integrated temperature for MSE if ( id_inttemp>0 ) then a2 = 0. do k=1,npz From 9691bbc533488bcd6bab9741358c84ce71bd0bfd Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Tue, 2 Apr 2024 14:08:58 -0400 Subject: [PATCH 3/6] fv_diagnostics.h removed comments --- tools/fv_diagnostics.h | 1 - 1 file changed, 1 deletion(-) diff --git a/tools/fv_diagnostics.h b/tools/fv_diagnostics.h index 108b8ee92..1462075f0 100644 --- a/tools/fv_diagnostics.h +++ b/tools/fv_diagnostics.h @@ -95,7 +95,6 @@ integer :: id_qr_dt_phys, id_qg_dt_phys, id_qs_dt_phys integer :: id_liq_wat_dt_phys, id_ice_wat_dt_phys integer :: id_intqv, id_intql, id_intqi, id_intqr, id_intqs, id_intqg - ! Tim M. mod integer :: id_inthght, id_inttemp ! ESM/CM 3-D diagostics From d5ac1e40d64ac7bdac40c45edb9f3a2bff0334fb Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Tue, 2 Apr 2024 14:21:33 -0400 Subject: [PATCH 4/6] fv_diagnostics.F90 updated to include sclark edits --- tools/fv_diagnostics.F90 | 82 ++++------------------------------------ tools/fv_diagnostics.h | 3 ++ 2 files changed, 10 insertions(+), 75 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index fe66b1464..64997e3d3 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -960,6 +960,10 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) ! 'Relative Humidity', '%', missing_value=missing_value, range=rhrange ) id_delp = register_diag_field ( trim(field), 'delp', axes(1:3), Time, & 'pressure thickness'//massdef_str, 'pa', missing_value=missing_value ) +#ifdef GFS_PHYS + id_delp_total = register_diag_field ( trim(field), 'delp_total', axes(1:3), Time, & + 'FV3 pressure thickness (dry air + all water species)', 'pa', missing_value=missing_value ) +#endif if ( .not. Atm(n)%flagstruct%hydrostatic ) & id_delz = register_diag_field ( trim(field), 'delz', axes(1:3), Time, & 'height thickness', 'm', missing_value=missing_value ) @@ -1080,7 +1084,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 ) @@ -3099,6 +3103,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) enddo if (id_delp > 0) used=send_data(id_delp, wk, Time) endif + if(id_delp_total > 0) used=send_data(id_delp_total, Atm(n)%delp(isc:iec,jsc:jec,:), Time) #else if(id_delp > 0) used=send_data(id_delp, Atm(n)%delp(isc:iec,jsc:jec,:), Time) #endif @@ -3915,7 +3920,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) do itrac=1, Atm(n)%ncnst call get_tracer_names (MODEL_ATMOS, itrac, tname) - if (id_tracer(itrac) > 0 .and. itrac.gt.nq) then + if (itrac.gt.nq) then used = send_data (id_tracer(itrac), Atm(n)%qdiag(isc:iec,jsc:jec,:,itrac), Time ) else used = send_data (id_tracer(itrac), Atm(n)%q(isc:iec,jsc:jec,:,itrac), Time ) @@ -5914,79 +5919,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) type(fv_grid_bounds_type), intent(IN) :: bd diff --git a/tools/fv_diagnostics.h b/tools/fv_diagnostics.h index 1462075f0..a8deda83e 100644 --- a/tools/fv_diagnostics.h +++ b/tools/fv_diagnostics.h @@ -107,4 +107,7 @@ integer :: id_uw, id_vw integer :: id_lagrangian_tendency_of_hydrostatic_pressure integer :: id_t_dt_nudge, id_ps_dt_nudge, id_delp_dt_nudge, id_u_dt_nudge, id_v_dt_nudge +#ifdef GFS_PHYS + integer :: id_delp_total +#endif #endif _FV_DIAG__ From 891b3d04ef38f716fc98080602bd1e4851861a0f Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Wed, 3 Apr 2024 15:44:30 -0400 Subject: [PATCH 5/6] fixed typo in description of vertically integrated temperature diag --- tools/fv_diagnostics.F90 | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 64997e3d3..cd5b90661 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -1153,7 +1153,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) id_inthght = register_diag_field ( trim(field), 'inthght', axes(1:2), Time, & 'Vertically Integrated Height', 'kg/m', missing_value=missing_value ) id_inttemp = register_diag_field ( trim(field), 'inttemp', axes(1:2), Time, & - 'Vertically Integrated Temperautre', 'K kg/m**2', missing_value=missing_value ) + 'Vertically Integrated Temperature', 'K kg/m**2', missing_value=missing_value ) id_acl = register_diag_field ( trim(field), 'acl', axes(1:2), Time, & 'Column-averaged Cl mixing ratio', 'kg/kg', missing_value=missing_value ) id_acl2 = register_diag_field ( trim(field), 'acl2', axes(1:2), Time, & @@ -1682,12 +1682,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) From 8a3005f529ca88e924265671c6823641c208669b Mon Sep 17 00:00:00 2001 From: Timothy Merlis Date: Tue, 9 Apr 2024 10:22:32 -0400 Subject: [PATCH 6/6] fixed conditional for any_hght to include >0 in all cases --- tools/fv_diagnostics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index cd5b90661..d646a3e73 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -872,7 +872,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref) if (id_t_dt_phys_plev_ave > 0 .and. .not. allocated(Atm(n)%phys_diag%phys_t_dt) ) allocate(Atm(n)%phys_diag%phys_t_dt(isc:iec,jsc:jec,npz)) ! flag for calculation of geopotential - if ( any(id_h > 0) .or. id_h_plev>0 .or. id_hght3d>0 .or. id_inthght .or. id_uz .or. id_iuz .or. id_vz .or. id_ivz ) then + if ( any(id_h > 0) .or. id_h_plev>0 .or. id_hght3d>0 .or. id_inthght>0 .or. id_uz>0 .or. id_iuz>0 .or. id_vz>0 .or. id_ivz>0 ) then id_any_hght = 1 else id_any_hght = 0 @@ -2287,7 +2287,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) if (id_hght3d > 0) then used = send_data(id_hght3d, 0.5*(wz(isc:iec,jsc:jec,1:npz)+wz(isc:iec,jsc:jec,2:npz+1)), Time) endif - if ( id_inthght>0 ) then + if (id_inthght > 0) then a2 = 0. do k=1,npz do j=jsc,jec