Skip to content

Update params #521

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 8 commits into from
+10 −9
Closed

Update params #521

Changes from all commits
Commits
Show all changes
8 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
b138219
Parameters update to be consistent with the corresponding constants i…
syha 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
19 changes: 10 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
Member

Choose a reason for hiding this comment

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

Suggested change
real(r8), parameter :: cp = 7.*rgas/2. ! = 1004.5
real(r8), parameter :: cp = 7.0_r8*rgas/2.0_r8 ! = 1004.5

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Looks better to me.

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 @@ -7380,7 +7381,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]
Copy link
Member

Choose a reason for hiding this comment

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

There are a bunch of subroutine arguments labelled as dry density. Should these all be dry air density?

Copy link
Member

Choose a reason for hiding this comment

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

edit: it is one in compute_full_pressure

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You are right. Here, rho is all 'dry air density'. Nice catch!

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 @@ -7408,7 +7409,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 @@ -7457,7 +7458,7 @@ 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)
pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero)
Copy link
Collaborator

Choose a reason for hiding this comment

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

i wanted to leave a comment on line 7445 - the 'dry density' comment should be 'dry air density' to be consistent. i think that's the last place 'dry density' is in a comment (that i could find) that needs changing. if all these suggested changes are ok, should we go ahead and accept them?

Copy link
Member

Choose a reason for hiding this comment

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

@nancycollins I'm going to close this pull request and put the constant fix in cleanly.
As Kevin says, it is probably best not to have the Jedi equations in the commits, and we need to add documentation adding about the constant changes.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@hkershaw-brown, just for clarification, there are no jedi equations or comments here at all. In this PR, I only updated the constants based on the model. I'm curious how you can make the fix cleaner?

Copy link
Member

Choose a reason for hiding this comment

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

So the Jedi stuff is in the commit history:

Screen Shot 2023-08-02 at 4 31 25 PM

Copy link
Member

Choose a reason for hiding this comment

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

@syha I'm aware you're busy with other things, so I'll put the fix in (I'll add you as a joint commit so you show up in the contributions). Thanks, Helen

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ah, the commit history, not the code itself. Got it. Thanks!

end where

if ( debug > 1 ) then
Expand Down