Skip to content

Commit

Permalink
Merge branch 'NOAA-GFDL:main' into diss_est_fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@JosephMouallem
JosephMouallem authored Oct 10, 2023
2 parents 9f507e2 + c59f36d commit 03128ba
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 93 deletions.
28 changes: 14 additions & 14 deletions model/fv_mapz.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)) )
Expand Down Expand Up @@ -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)) )
Expand Down Expand Up @@ -366,21 +366,21 @@ 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

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))
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)

Expand All @@ -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)) )
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
80 changes: 1 addition & 79 deletions tools/fv_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand Down

0 comments on commit 03128ba

Please sign in to comment.