From 52349450f645bef796ad7582914e25ef27edd352 Mon Sep 17 00:00:00 2001 From: Peter Hjort Lauritzen Date: Wed, 9 Oct 2024 13:29:42 -0600 Subject: [PATCH 1/7] Allow custom surface pressure in `std_atm_pres` Co-authored-by: Jesse Nusbaumer --- src/utils/std_atm_profile.F90 | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/utils/std_atm_profile.F90 b/src/utils/std_atm_profile.F90 index f0821d1f..d7ddc558 100644 --- a/src/utils/std_atm_profile.F90 +++ b/src/utils/std_atm_profile.F90 @@ -52,16 +52,28 @@ 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 + integer :: ierr + real(r8) :: pb_local(nreg) - integer :: i, ii, k, nlev character(len=*), parameter :: routine = 'std_atm_pres' !---------------------------------------------------------------------------- + ! Initialize local standard pressure values array + pb_local = pb + + ! Set new surface pressure value if provided by the caller + 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 @@ -78,13 +90,12 @@ 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 - end subroutine std_atm_pres !========================================================================================= From 89a183f4af664c068bff6c60e8091cdc52652a7b Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 23 Oct 2024 16:44:04 -0600 Subject: [PATCH 2/7] Remove unused variable in `std_atm_pres` --- src/utils/std_atm_profile.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/utils/std_atm_profile.F90 b/src/utils/std_atm_profile.F90 index d7ddc558..aff962e4 100644 --- a/src/utils/std_atm_profile.F90 +++ b/src/utils/std_atm_profile.F90 @@ -60,7 +60,6 @@ subroutine std_atm_pres(height, pstd, user_specified_ps) real(r8), optional, intent(in) :: user_specified_ps integer :: i, ii, k, nlev - integer :: ierr real(r8) :: pb_local(nreg) character(len=*), parameter :: routine = 'std_atm_pres' From 446f987eb1ef9c7dc39019fa65c5d214ac647612 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 23 Oct 2024 16:47:32 -0600 Subject: [PATCH 3/7] Properly call `std_atm_pres` to initialize reference pressure for physics --- src/dynamics/mpas/dyn_grid.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 26b0e345..565db35f 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -11,7 +11,7 @@ module dyn_grid ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, & ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & sphere_radius - use dynconst, only: constant_pi => pi, rad_to_deg, dynconst_init + use dynconst, only: constant_p0 => pref, constant_pi => pi, rad_to_deg, dynconst_init use physics_column_type, only: kind_pcol, physics_column_t use physics_grid, only: phys_decomp, phys_grid_init use ref_pres, only: ref_pres_init @@ -107,7 +107,7 @@ subroutine model_grid_init() call endrun('Numbers of vertical layers mismatch', subname, __LINE__) end if - ! Initialize reference pressure. + ! Initialize reference pressure for use by physics. call dyn_debug_print('Calling init_reference_pressure') call init_reference_pressure() @@ -123,7 +123,7 @@ subroutine model_grid_init() call define_cam_grid() end subroutine model_grid_init - !> Initialize reference pressure by computing necessary variables and calling `ref_pres_init`. + !> Initialize reference pressure for use by physics. !> (KCW, 2024-03-25) subroutine init_reference_pressure() character(*), parameter :: subname = 'dyn_grid::init_reference_pressure' @@ -173,7 +173,7 @@ subroutine init_reference_pressure() allocate(p_ref_int(pverp), stat=ierr) call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__) - call std_atm_pres(zw, p_ref_int) + call std_atm_pres(zw, p_ref_int, user_specified_ps=constant_p0) allocate(p_ref_mid(pver), stat=ierr) call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__) From 9c0998f2cec83e5822ebec8b07097a3cfef1a5c6 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sun, 27 Oct 2024 06:25:31 -0600 Subject: [PATCH 4/7] Wire up history support --- src/dynamics/mpas/dyn_grid.F90 | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 565db35f..39e1e5f4 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -5,6 +5,7 @@ module dyn_grid use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, & horiz_coord_t, horiz_coord_create, & max_hcoordname_len + use cam_history_support, only: add_vert_coord use cam_initfiles, only: initial_file_get_id use cam_map_utils, only: kind_imap => imap use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & @@ -139,10 +140,16 @@ subroutine init_reference_pressure() ! `zw` denotes zeta at w-wind levels (i.e., at layer interfaces). ! `dzw` denotes the delta/difference between `zw`. ! `rdzw` denotes the reciprocal of `dzw`. - real(kind_r8), allocatable :: zu(:), zw(:), dzw(:) + real(kind_r8), allocatable :: dzw(:) real(kind_r8), pointer :: rdzw(:) + real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord` + ! just uses pointers to point at it internally. + real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord` + ! just uses pointers to point at it internally. nullify(rdzw) + nullify(zu) + nullify(zw) ! Compute reference height. call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw') @@ -169,6 +176,12 @@ subroutine init_reference_pressure() zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k)) end do + ! Register zeta coordinates with history. + call add_vert_coord('ilev', pverp, 'Height (zeta) level at layer interfaces', 'm', zw, & + positive='up') + call add_vert_coord('lev', pver, 'Height (zeta) level at layer midpoints', 'm', zu, & + positive='up') + ! Compute reference pressure from reference height. allocate(p_ref_int(pverp), stat=ierr) call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__) From 8cda5343a8dfc705973399ba683eefc3de641d18 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sun, 27 Oct 2024 06:28:21 -0600 Subject: [PATCH 5/7] More explicit memory management --- src/dynamics/mpas/dyn_grid.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 39e1e5f4..5889ff67 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -176,6 +176,8 @@ subroutine init_reference_pressure() zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k)) end do + deallocate(dzw) + ! Register zeta coordinates with history. call add_vert_coord('ilev', pverp, 'Height (zeta) level at layer interfaces', 'm', zw, & positive='up') @@ -214,6 +216,12 @@ subroutine init_reference_pressure() ' | ' // stringify([p_ref_int(k) / 100.0_kind_r8])) call ref_pres_init(p_ref_int, p_ref_mid, num_pure_p_lev) + + deallocate(p_ref_int) + deallocate(p_ref_mid) + + nullify(zu) + nullify(zw) end subroutine init_reference_pressure !> Initialize physics grid in terms of dynamics decomposition. @@ -287,6 +295,9 @@ subroutine init_physics_grid() call check_allocate(ierr, subname, 'dyn_attribute_name(0)', 'dyn_grid', __LINE__) call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name) + + deallocate(dyn_column) + deallocate(dyn_attribute_name) end subroutine init_physics_grid !> This subroutine defines and registers four variants of dynamics grids in terms of dynamics decomposition. From 35d99c6bb6042fd61eddab10cc2deea4eef3f293 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sun, 27 Oct 2024 06:29:14 -0600 Subject: [PATCH 6/7] Update code comments --- src/dynamics/mpas/dyn_grid.F90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 5889ff67..07277b1c 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -230,7 +230,7 @@ end subroutine init_reference_pressure subroutine init_physics_grid() character(*), parameter :: subname = 'dyn_grid::init_physics_grid' character(max_hcoordname_len), allocatable :: dyn_attribute_name(:) - integer :: hdim1_d, hdim2_d + integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid. integer :: i integer :: ierr integer, pointer :: indextocellid(:) ! Global indexes of cell centers. @@ -322,25 +322,25 @@ subroutine define_cam_grid() real(kind_r8), pointer :: lonedge(:) ! Edge node longitudes (radians). real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians). - ! Global grid indexes. CAN be safely deallocated because its values are copied into - ! `cam_grid_attribute_*_t` and `horiz_coord_t`. + ! Global grid indexes. CAN be safely deallocated because its values are copied internally by + ! `cam_grid_attribute_register` and `horiz_coord_create`. ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`. integer(kind_imap), pointer :: global_grid_index(:) - ! Global grid maps. CANNOT be safely deallocated because `cam_filemap_t` - ! just uses pointers to point at it. + ! Global grid maps. CANNOT be safely deallocated because `cam_grid_register` + ! just uses pointers to point at it internally. ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`. integer(kind_imap), pointer :: global_grid_map(:, :) - ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_*_t` - ! just uses pointers to point at it. + ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_register` + ! just uses pointers to point at it internally. real(kind_r8), pointer :: cell_area(:) - ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_*_t` - ! just uses pointers to point at it. + ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_register` + ! just uses pointers to point at it internally. real(kind_r8), pointer :: cell_weight(:) - ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_t` - ! just uses pointers to point at it. + ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_register` + ! just uses pointers to point at it internally. type(horiz_coord_t), pointer :: lat_coord - ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_t` - ! just uses pointers to point at it. + ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_register` + ! just uses pointers to point at it internally. type(horiz_coord_t), pointer :: lon_coord nullify(indextocellid, indextoedgeid, indextovertexid) From 3428563edc0f4e8f2bd639149763d135cf80bdee Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 28 Oct 2024 13:00:46 -0600 Subject: [PATCH 7/7] Sort statements The standard name was updated. It should have been sorted again at that time. --- src/dynamics/mpas/dyn_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 30e4b2a5..b1892f47 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -812,7 +812,6 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('eastward_wind') call mark_as_initialized('geopotential_height_wrt_surface') call mark_as_initialized('geopotential_height_wrt_surface_at_interface') - call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure') call mark_as_initialized('lagrangian_tendency_of_air_pressure') call mark_as_initialized('ln_air_pressure') call mark_as_initialized('ln_air_pressure_at_interface') @@ -821,6 +820,7 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('northward_wind') call mark_as_initialized('reciprocal_of_air_pressure_thickness') call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air') + call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure') call mark_as_initialized('surface_air_pressure') call mark_as_initialized('surface_geopotential') call mark_as_initialized('surface_pressure_of_dry_air')