From ef4ac155fdadbc66d9ebbdd2333d8c548144fbfc Mon Sep 17 00:00:00 2001 From: syha Date: Thu, 17 Nov 2022 15:59:01 -0700 Subject: [PATCH 1/5] Bug fix for full pressure with the additional term for qv^2. --- models/mpas_atm/model_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index eb42c3ec5c..d1aa1e5ff7 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -7506,7 +7506,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 + pressure = rho * rgas * tk * (1.0_r8 + 1.61_r8 * qv_nonzero + rvordm1 * qv_nonzero**2) end where if ( debug > 1 ) then From 970519166ca141b29fb5519e7cb5a95437bf2e69 Mon Sep 17 00:00:00 2001 From: syha Date: Thu, 17 Nov 2022 16:03:48 -0700 Subject: [PATCH 2/5] rvordm1 was missed, which is now back in. --- models/mpas_atm/model_mod.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index d1aa1e5ff7..3eee0a36d0 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -216,7 +216,8 @@ module model_mod real(r8), parameter :: cv = 716.0_r8 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. From b405613a0c5ac0aa31596a52547b8e6b0d29f6d6 Mon Sep 17 00:00:00 2001 From: syha Date: Thu, 17 Nov 2022 16:08:01 -0700 Subject: [PATCH 3/5] Parameters are updated for the variable transformation, to be consistent with the model. --- models/mpas_atm/model_mod.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index 3eee0a36d0..0aa54cea42 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -209,11 +209,11 @@ 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_d (Gas constant for dry air [J kg-1 K-1]) +real(r8), parameter :: cp = 7.*rgas/2. ! = 1004.5 +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 ! = 1.6083623693379792 @@ -7430,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] @@ -7458,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 @@ -7508,7 +7508,7 @@ subroutine compute_full_pressure(ens_size, theta, rho, qv, pressure, tk, istatus where (istatus == 0) ! We only take non-missing tk here ! Soyoung (Nov-2022): The last term is added for full pressure P = (rho_d*R_d + rho_v*R_v)*tk - pressure = rho * rgas * tk * (1.0_r8 + 1.61_r8 * qv_nonzero + rvordm1 * qv_nonzero**2) + pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero + rvordm1 * qv_nonzero**2) end where if ( debug > 1 ) then From 9a069ff38738c23d19aec5bc47f6f73ea768777f Mon Sep 17 00:00:00 2001 From: syha Date: Thu, 17 Nov 2022 16:20:44 -0700 Subject: [PATCH 4/5] Comment was incorrect for rv. Now fixed. --- models/mpas_atm/model_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index 0aa54cea42..9950b05b88 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -211,7 +211,7 @@ module model_mod ! Real (physical) constants as defined exactly in MPAS. ! 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_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 real(r8), parameter :: cv = cp-rgas ! = 717.5 real(r8), parameter :: p0 = 100000.0_r8 From b138219e29b0c5de545f0fe0796841b25ea00bd1 Mon Sep 17 00:00:00 2001 From: syha Date: Thu, 27 Jul 2023 15:47:56 -0600 Subject: [PATCH 5/5] Parameters update to be consistent with the corresponding constants in the mpas model. --- models/mpas_atm/model_mod.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/models/mpas_atm/model_mod.f90 b/models/mpas_atm/model_mod.f90 index c93168f339..e8037e417a 100644 --- a/models/mpas_atm/model_mod.f90 +++ b/models/mpas_atm/model_mod.f90 @@ -7458,8 +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 - ! Soyoung (Nov-2022): The last term is added for full pressure P = (rho_d*R_d + rho_v*R_v)*tk - pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero + rvordm1 * qv_nonzero**2) + pressure = rho * rgas * tk * (1.0_r8 + rvord * qv_nonzero) end where if ( debug > 1 ) then