Skip to content

Fix mpas_atm compute_full_pressure #429

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

Closed
wants to merge 7 commits into from
Closed

Fix mpas_atm compute_full_pressure #429

Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ef4ac15
Bug fix for full pressure with the additional term for qv^2.
syha Nov 17, 2022
9705191
rvordm1 was missed, which is now back in.
syha Nov 17, 2022
b405613
Parameters are updated for the variable transformation, to be consist…
syha Nov 17, 2022
9a069ff
Comment was incorrect for rv. Now fixed.
syha Nov 17, 2022
669fc2c
Merge branch 'main' into fix-mpas-atm-pressure
hkershaw-brown May 10, 2023
92f54a3
Merge branch 'main' into fix-mpas-atm-pressure
hkershaw-brown Jul 27, 2023
41b423f
Merge branch 'main' into fix-mpas-atm-pressure
hkershaw-brown Jul 27, 2023
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
20 changes: 11 additions & 9 deletions models/mpas_atm/model_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,14 +209,15 @@ module model_mod
integer, parameter :: TIMELEN = 19

! Real (physical) constants as defined exactly in MPAS.
! redefined here for consistency with the model.
real(r8), parameter :: rgas = 287.0_r8
real(r8), parameter :: rv = 461.6_r8
real(r8), parameter :: cp = 1003.0_r8
real(r8), parameter :: cv = 716.0_r8
! redefined here for consistency with the model (MPAS/src/framework/mpas_constants.F).
real(r8), parameter :: rgas = 287.0_r8 ! = R_d (Gas constant for dry air [J kg-1 K-1])
real(r8), parameter :: rv = 461.6_r8 ! = R_v (Gas constant for water varpor [J kg-1 K-1])
real(r8), parameter :: cp = 7.*rgas/2. ! = 1004.5
Copy link
Contributor

Choose a reason for hiding this comment

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

Another comment might be helpful here: "For an ideal diatomic gas, cp:cv:R is 7:5:2"

real(r8), parameter :: cv = cp-rgas ! = 717.5
real(r8), parameter :: p0 = 100000.0_r8
real(r8), parameter :: rcv = rgas/(cp-rgas)
real(r8), parameter :: rvord = rv/rgas
real(r8), parameter :: rvord = rv/rgas ! = 1.6083623693379792
real(r8), parameter :: rvordm1 = rv/rgas-1.0_r8 ! = 0.6083623693379792

! earth radius; needed to convert lat/lon to x,y,z cartesian coords.
! for the highest accuracy this should match what the model uses.
Expand Down Expand Up @@ -7429,7 +7430,7 @@ function theta_to_tk (ens_size, theta, rho, qv, istatus)

integer, intent(in) :: ens_size
real(r8), dimension(ens_size), intent(in) :: theta ! potential temperature [K]
real(r8), dimension(ens_size), intent(in) :: rho ! dry density
real(r8), dimension(ens_size), intent(in) :: rho ! dry air density [kg/m3]
real(r8), dimension(ens_size), intent(in) :: qv ! water vapor mixing ratio [kg/kg]
integer, dimension(ens_size), intent(inout) :: istatus
real(r8), dimension(ens_size) :: theta_to_tk ! sensible temperature [K]
Expand Down Expand Up @@ -7457,7 +7458,7 @@ function theta_to_tk (ens_size, theta, rho, qv, istatus)
endif
where (istatus == 0)

theta_m = (1.0_r8 + 1.61_r8 * qv_nonzero)*theta
theta_m = (1.0_r8 + rvord * qv_nonzero)*theta

where (theta_m > 0.0_r8 .and. rho > 0.0_r8) ! Check if all the input are positive

Expand Down Expand Up @@ -7506,7 +7507,8 @@ subroutine compute_full_pressure(ens_size, theta, rho, qv, pressure, tk, istatus
tk = theta_to_tk(ens_size, theta, rho, qv_nonzero, istatus)

where (istatus == 0) ! We only take non-missing tk here
pressure = rho * rgas * tk * (1.0_r8 + 1.61_r8 * qv_nonzero)
! Soyoung (Nov-2022): The last term is added for full pressure P = (rho_d*R_d + rho_v*R_v)*tk
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't see yet where the qv**2 term comes from.
From the ideal gas law, as written in the comment here,
P = rho_d*R_d*tk * [1 + (rho_v/rho_d)(R_v/R_d)]
Translating into the variables and parameters defined in this module:
P = rho *rgas*tk * [1 + qv_nonzero * rvord]
which is the unmodified code. It would be helpful to confirm that
MPAS uses mixing ratio (mass_v/mass_d),
not specific humidity (mass_v/mass_(v+d))?

pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero + rvordm1 * qv_nonzero**2)
end where

if ( debug > 1 ) then
Expand Down