Skip to content

Commit

Permalink
fix issue ESCOMP#1168 (MPAS reference pressures)
Browse files Browse the repository at this point in the history
  • Loading branch information
PeterHjortLauritzen committed Oct 9, 2024
1 parent 1abe4a8 commit 2d2f08a
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 8 deletions.
2 changes: 1 addition & 1 deletion src/dynamics/mpas/dyn_grid.F90
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ subroutine dyn_grid_init()
ierr = pio_get_att(fh_ini, pio_global, 'sphere_radius', sphere_radius)

! Compute reference pressures from reference heights.
call std_atm_pres(zw, pref_edge)
call std_atm_pres(zw, pref_edge, user_specified_ps=1.0E5_r8)
pref_mid = (pref_edge(1:plev) + pref_edge(2:plevp)) * 0.5_r8

num_pr_lev = 0
Expand Down
21 changes: 14 additions & 7 deletions src/utils/std_atm_profile.F90
Original file line number Diff line number Diff line change
Expand Up @@ -52,16 +52,23 @@ module std_atm_profile
CONTAINS
!=========================================================================================

subroutine std_atm_pres(height, pstd)
subroutine std_atm_pres(height, pstd, user_specified_ps)

! arguments
real(r8), intent(in) :: height(:) ! height above sea level in meters
real(r8), intent(out) :: pstd(:) ! std pressure in Pa
real(r8), intent(in) :: height(:) ! height above sea level in meters
real(r8), intent(out) :: pstd(:) ! std pressure in Pa
real(r8), optional, intent(in) :: user_specified_ps

integer :: i, ii, k, nlev
character(len=*), parameter :: routine = 'std_atm_pres'
real(r8), allocatable :: pb_local(:)
!----------------------------------------------------------------------------

allocate(pb_local(size(height)))
pb_local = pb
if (present(user_specified_ps)) then
pb_local(1) = user_specified_ps
end if

nlev = size(height)
do k = 1, nlev
if (height(k) < 0.0_r8) then
Expand All @@ -78,13 +85,13 @@ subroutine std_atm_pres(height, pstd)
end if

if (lb(ii) /= 0._r8) then
pstd(k) = pb(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
pstd(k) = pb_local(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
else
pstd(k) = pb(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
pstd(k) = pb_local(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
end if

end do

deallocate(pb_local)
end subroutine std_atm_pres

!=========================================================================================
Expand Down

0 comments on commit 2d2f08a

Please sign in to comment.