Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fv_diagnostics update #332

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

fv_diagnostics update #332

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
fc71858
fv_diag added height fluxes, corrected logic for integrated fluxes
Apr 2, 2024
c9bf9d8
fv_diagnostics removed comments
Apr 2, 2024
9691bbc
fv_diagnostics.h removed comments
Apr 2, 2024
d5ac1e4
fv_diagnostics.F90 updated to include sclark edits
Apr 2, 2024
891b3d0
fixed typo in description of vertically integrated temperature diag
Apr 3, 2024
8a3005f
fixed conditional for any_hght to include >0 in all cases
Apr 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 90 additions & 17 deletions tools/fv_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
tmerlis marked this conversation as resolved.
Show resolved Hide resolved
id_any_hght = 1
else
id_any_hght = 0
Expand Down Expand Up @@ -1017,6 +1017,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, &
Expand Down Expand Up @@ -1050,6 +1055,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, &
Expand Down Expand Up @@ -1140,7 +1150,10 @@ 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 )

id_inthght = register_diag_field ( trim(field), 'inthght', axes(1:2), Time, &
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are inthght and inttemp supposed to be proportional to the vertically-integrated potential and internal energies?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that's correct. This could be overt by multiplying by the relevant constants 'grav' and 'cp_air' that are available in the module and changing the unit accordingly.

'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 )
tmerlis marked this conversation as resolved.
Show resolved Hide resolved
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, &
Expand Down Expand Up @@ -1669,6 +1682,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)
Expand Down Expand Up @@ -2263,7 +2282,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, &
Expand All @@ -2274,6 +2293,48 @@ 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
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)
Expand Down Expand Up @@ -2845,6 +2906,18 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
used = send_data(id_intqg, a2*ginv, Time)
endif

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) )
Expand Down Expand Up @@ -3016,6 +3089,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. &
Expand Down Expand Up @@ -3914,95 +3988,95 @@ 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
a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%q(i,j,k,sphum)
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
a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%q(i,j,k,sphum)
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)
endif
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
a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%pt(i,j,k)
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
a4(i,j,k) = Atm(n)%va(i,j,k) * Atm(n)%pt(i,j,k)
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)
endif
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
a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%ua(i,j,k)
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
a4(i,j,k) = Atm(n)%ua(i,j,k) * Atm(n)%va(i,j,k)
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)
Expand Down Expand Up @@ -5845,7 +5919,6 @@ end subroutine eqv_pot

#endif


subroutine compute_brn(ua, va, delp, delz, cape, bd, npz, Time)

type(fv_grid_bounds_type), intent(IN) :: bd
Expand Down
2 changes: 2 additions & 0 deletions tools/fv_diagnostics.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,10 +95,12 @@
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
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

Expand Down