Skip to content

Commit 9cc084a

Browse files
kuanchihwangPeterHjortLauritzennusbaume
authored
Fix initialization of reference pressure for physics when using MPAS dynamical core (ESCOMP#317)
### Tag name (required for release branches): None ### Originator(s): PeterHjortLauritzen, nusbaume, kuanchihwang ### Description (include the issue title, and the keyword ['closes', 'fixes', 'resolves'] followed by the issue number): Presently, MPAS dynamical core is not initializing reference pressure for physics in the same way as others. This PR backports the fix from ESCOMP/CAM#1169. To confirm the fix, observe log entries similar to the following in `atm.log.<job-id>.<date>-<time>`. The reference pressure at surface now starts from 1000 hPa. Previously, it starts from 1013.25 hPa. ``` dyn_debug_print (0): Reference layer information: dyn_debug_print (0): ----- | -------------- | -------------- dyn_debug_print (0): Index | Height (m) | Pressure (hPa) dyn_debug_print (0): ----- | 44999.999819 | 1.431394 dyn_debug_print (0): 1 | 43953.595504 | 1.655836 dyn_debug_print (0): ----- | 42907.191189 | 1.880278 dyn_debug_print (0): 2 | 41877.531046 | 2.177235 dyn_debug_print (0): ----- | 40847.870904 | 2.474192 dyn_debug_print (0): 3 | 39835.232205 | 2.867591 dyn_debug_print (0): ----- | 38822.593507 | 3.260989 ... (SNIPPED) ... dyn_debug_print (0): ----- | 1292.109873 | 856.014954 dyn_debug_print (0): 30 | 997.820649 | 887.686997 dyn_debug_print (0): ----- | 703.531425 | 919.359041 dyn_debug_print (0): 31 | 476.271636 | 945.093654 dyn_debug_print (0): ----- | 249.011847 | 970.828267 dyn_debug_print (0): 32 | 124.505924 | 985.414133 dyn_debug_print (0): ----- | 0.000000 | 1000.000000 ``` Closes ESCOMP#315. ### Describe any changes made to build system: None ### Describe any changes made to the namelist: None ### List any changes to the defaults for the input datasets (e.g. boundary datasets): None ### List all files eliminated and why: None ### List all files added and what they do: None ### List all existing files that have been modified, and describe the changes: * `M src/dynamics/mpas/dyn_comp.F90` * Sort statements * `M src/dynamics/mpas/dyn_grid.F90` * Properly call `std_atm_pres` to initialize reference pressure for physics * Wire up history support * More explicit memory management * Update code comments * `M src/utils/std_atm_profile.F90` * Allow custom surface pressure in `std_atm_pres` * Remove unused variable in `std_atm_pres` --------- Co-authored-by: Peter Hjort Lauritzen <pel@ucar.edu> Co-authored-by: Jesse Nusbaumer <nusbaume@ucar.edu>
1 parent 455003c commit 9cc084a

File tree

3 files changed

+60
-26
lines changed

3 files changed

+60
-26
lines changed

‎src/dynamics/mpas/dyn_comp.F90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -812,7 +812,6 @@ subroutine mark_variable_as_initialized()
812812
call mark_as_initialized('eastward_wind')
813813
call mark_as_initialized('geopotential_height_wrt_surface')
814814
call mark_as_initialized('geopotential_height_wrt_surface_at_interface')
815-
call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure')
816815
call mark_as_initialized('lagrangian_tendency_of_air_pressure')
817816
call mark_as_initialized('ln_air_pressure')
818817
call mark_as_initialized('ln_air_pressure_at_interface')
@@ -821,6 +820,7 @@ subroutine mark_variable_as_initialized()
821820
call mark_as_initialized('northward_wind')
822821
call mark_as_initialized('reciprocal_of_air_pressure_thickness')
823822
call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air')
823+
call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure')
824824
call mark_as_initialized('surface_air_pressure')
825825
call mark_as_initialized('surface_geopotential')
826826
call mark_as_initialized('surface_pressure_of_dry_air')

‎src/dynamics/mpas/dyn_grid.F90‎

Lines changed: 42 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -5,13 +5,14 @@ module dyn_grid
55
use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, &
66
horiz_coord_t, horiz_coord_create, &
77
max_hcoordname_len
8+
use cam_history_support, only: add_vert_coord
89
use cam_initfiles, only: initial_file_get_id
910
use cam_map_utils, only: kind_imap => imap
1011
use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, &
1112
ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, &
1213
ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, &
1314
sphere_radius
14-
use dynconst, only: constant_pi => pi, rad_to_deg, dynconst_init
15+
use dynconst, only: constant_p0 => pref, constant_pi => pi, rad_to_deg, dynconst_init
1516
use physics_column_type, only: kind_pcol, physics_column_t
1617
use physics_grid, only: phys_decomp, phys_grid_init
1718
use ref_pres, only: ref_pres_init
@@ -107,7 +108,7 @@ subroutine model_grid_init()
107108
call endrun('Numbers of vertical layers mismatch', subname, __LINE__)
108109
end if
109110

110-
! Initialize reference pressure.
111+
! Initialize reference pressure for use by physics.
111112
call dyn_debug_print('Calling init_reference_pressure')
112113

113114
call init_reference_pressure()
@@ -123,7 +124,7 @@ subroutine model_grid_init()
123124
call define_cam_grid()
124125
end subroutine model_grid_init
125126

126-
!> Initialize reference pressure by computing necessary variables and calling `ref_pres_init`.
127+
!> Initialize reference pressure for use by physics.
127128
!> (KCW, 2024-03-25)
128129
subroutine init_reference_pressure()
129130
character(*), parameter :: subname = 'dyn_grid::init_reference_pressure'
@@ -139,10 +140,16 @@ subroutine init_reference_pressure()
139140
! `zw` denotes zeta at w-wind levels (i.e., at layer interfaces).
140141
! `dzw` denotes the delta/difference between `zw`.
141142
! `rdzw` denotes the reciprocal of `dzw`.
142-
real(kind_r8), allocatable :: zu(:), zw(:), dzw(:)
143+
real(kind_r8), allocatable :: dzw(:)
143144
real(kind_r8), pointer :: rdzw(:)
145+
real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord`
146+
! just uses pointers to point at it internally.
147+
real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord`
148+
! just uses pointers to point at it internally.
144149

145150
nullify(rdzw)
151+
nullify(zu)
152+
nullify(zw)
146153

147154
! Compute reference height.
148155
call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw')
@@ -169,11 +176,19 @@ subroutine init_reference_pressure()
169176
zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k))
170177
end do
171178

179+
deallocate(dzw)
180+
181+
! Register zeta coordinates with history.
182+
call add_vert_coord('ilev', pverp, 'Height (zeta) level at layer interfaces', 'm', zw, &
183+
positive='up')
184+
call add_vert_coord('lev', pver, 'Height (zeta) level at layer midpoints', 'm', zu, &
185+
positive='up')
186+
172187
! Compute reference pressure from reference height.
173188
allocate(p_ref_int(pverp), stat=ierr)
174189
call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__)
175190

176-
call std_atm_pres(zw, p_ref_int)
191+
call std_atm_pres(zw, p_ref_int, user_specified_ps=constant_p0)
177192

178193
allocate(p_ref_mid(pver), stat=ierr)
179194
call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__)
@@ -201,6 +216,12 @@ subroutine init_reference_pressure()
201216
' | ' // stringify([p_ref_int(k) / 100.0_kind_r8]))
202217

203218
call ref_pres_init(p_ref_int, p_ref_mid, num_pure_p_lev)
219+
220+
deallocate(p_ref_int)
221+
deallocate(p_ref_mid)
222+
223+
nullify(zu)
224+
nullify(zw)
204225
end subroutine init_reference_pressure
205226

206227
!> Initialize physics grid in terms of dynamics decomposition.
@@ -209,7 +230,7 @@ end subroutine init_reference_pressure
209230
subroutine init_physics_grid()
210231
character(*), parameter :: subname = 'dyn_grid::init_physics_grid'
211232
character(max_hcoordname_len), allocatable :: dyn_attribute_name(:)
212-
integer :: hdim1_d, hdim2_d
233+
integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid.
213234
integer :: i
214235
integer :: ierr
215236
integer, pointer :: indextocellid(:) ! Global indexes of cell centers.
@@ -274,6 +295,9 @@ subroutine init_physics_grid()
274295
call check_allocate(ierr, subname, 'dyn_attribute_name(0)', 'dyn_grid', __LINE__)
275296

276297
call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name)
298+
299+
deallocate(dyn_column)
300+
deallocate(dyn_attribute_name)
277301
end subroutine init_physics_grid
278302

279303
!> This subroutine defines and registers four variants of dynamics grids in terms of dynamics decomposition.
@@ -298,25 +322,25 @@ subroutine define_cam_grid()
298322
real(kind_r8), pointer :: lonedge(:) ! Edge node longitudes (radians).
299323
real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians).
300324

301-
! Global grid indexes. CAN be safely deallocated because its values are copied into
302-
! `cam_grid_attribute_*_t` and `horiz_coord_t`.
325+
! Global grid indexes. CAN be safely deallocated because its values are copied internally by
326+
! `cam_grid_attribute_register` and `horiz_coord_create`.
303327
! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
304328
integer(kind_imap), pointer :: global_grid_index(:)
305-
! Global grid maps. CANNOT be safely deallocated because `cam_filemap_t`
306-
! just uses pointers to point at it.
329+
! Global grid maps. CANNOT be safely deallocated because `cam_grid_register`
330+
! just uses pointers to point at it internally.
307331
! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
308332
integer(kind_imap), pointer :: global_grid_map(:, :)
309-
! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_*_t`
310-
! just uses pointers to point at it.
333+
! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_register`
334+
! just uses pointers to point at it internally.
311335
real(kind_r8), pointer :: cell_area(:)
312-
! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_*_t`
313-
! just uses pointers to point at it.
336+
! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_register`
337+
! just uses pointers to point at it internally.
314338
real(kind_r8), pointer :: cell_weight(:)
315-
! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
316-
! just uses pointers to point at it.
339+
! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_register`
340+
! just uses pointers to point at it internally.
317341
type(horiz_coord_t), pointer :: lat_coord
318-
! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
319-
! just uses pointers to point at it.
342+
! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_register`
343+
! just uses pointers to point at it internally.
320344
type(horiz_coord_t), pointer :: lon_coord
321345

322346
nullify(indextocellid, indextoedgeid, indextovertexid)

‎src/utils/std_atm_profile.F90‎

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -52,16 +52,27 @@ module std_atm_profile
5252
CONTAINS
5353
!=========================================================================================
5454

55-
subroutine std_atm_pres(height, pstd)
55+
subroutine std_atm_pres(height, pstd, user_specified_ps)
5656

5757
! arguments
58-
real(r8), intent(in) :: height(:) ! height above sea level in meters
59-
real(r8), intent(out) :: pstd(:) ! std pressure in Pa
58+
real(r8), intent(in) :: height(:) ! height above sea level in meters
59+
real(r8), intent(out) :: pstd(:) ! std pressure in Pa
60+
real(r8), optional, intent(in) :: user_specified_ps
61+
62+
integer :: i, ii, k, nlev
63+
real(r8) :: pb_local(nreg)
6064

61-
integer :: i, ii, k, nlev
6265
character(len=*), parameter :: routine = 'std_atm_pres'
6366
!----------------------------------------------------------------------------
6467

68+
! Initialize local standard pressure values array
69+
pb_local = pb
70+
71+
! Set new surface pressure value if provided by the caller
72+
if (present(user_specified_ps)) then
73+
pb_local(1) = user_specified_ps
74+
end if
75+
6576
nlev = size(height)
6677
do k = 1, nlev
6778
if (height(k) < 0.0_r8) then
@@ -78,13 +89,12 @@ subroutine std_atm_pres(height, pstd)
7889
end if
7990

8091
if (lb(ii) /= 0._r8) then
81-
pstd(k) = pb(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
92+
pstd(k) = pb_local(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
8293
else
83-
pstd(k) = pb(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
94+
pstd(k) = pb_local(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
8495
end if
8596

8697
end do
87-
8898
end subroutine std_atm_pres
8999

90100
!=========================================================================================

0 commit comments

Comments
 (0)