From 7fac3a9edffee7b18f885b49aee86fb244a3d369 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 11 Sep 2024 11:17:25 -0600 Subject: [PATCH 01/22] Sort `use` statements and adjust indentation --- src/dynamics/mpas/stepon.F90 | 97 ++++++++++++++++++------------------ 1 file changed, 49 insertions(+), 48 deletions(-) diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index 9d53ca43..7dfdbc23 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,9 +1,12 @@ module stepon - use camsrfexch, only: cam_out_t - use dyn_comp, only: dyn_import_t, dyn_export_t + ! Modules from CAM-SIMA. + use camsrfexch, only: cam_out_t + use dyn_comp, only: dyn_import_t, dyn_export_t use physics_types, only: physics_state, physics_tend - use runtime_obj, only: runtime_options - use shr_kind_mod, only: r8 => shr_kind_r8 + use runtime_obj, only: runtime_options + + ! Modules from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 implicit none @@ -15,48 +18,46 @@ module stepon public :: stepon_run3 public :: stepon_final contains - -! Called by `cam_init` in `src/control/cam_comp.F90`. -subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out) - type(runtime_options), intent(in) :: cam_runtime_opts - type(dyn_import_t), intent(in) :: dyn_in - type(dyn_export_t), intent(in) :: dyn_out -end subroutine stepon_init - -! Called by `cam_timestep_init` in `src/control/cam_comp.F90`. -subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) - real(r8), intent(out) :: dtime_phys - type(runtime_options), intent(in) :: cam_runtime_opts - type(physics_state), intent(inout) :: phys_state - type(physics_tend), intent(inout) :: phys_tend - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out -end subroutine stepon_timestep_init - -! Called by `cam_run2` in `src/control/cam_comp.F90`. -subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) - type(runtime_options), intent(in) :: cam_runtime_opts - type(physics_state), intent(inout) :: phys_state - type(physics_tend), intent(inout) :: phys_tend - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out -end subroutine stepon_run2 - -! Called by `cam_run3` in `src/control/cam_comp.F90`. -subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out) - real(r8), intent(in) :: dtime_phys - type(runtime_options), intent(in) :: cam_runtime_opts - type(cam_out_t), intent(inout) :: cam_out - type(physics_state), intent(inout) :: phys_state - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out -end subroutine stepon_run3 - -! Called by `cam_final` in `src/control/cam_comp.F90`. -subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) - type(runtime_options), intent(in) :: cam_runtime_opts - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out -end subroutine stepon_final - + ! Called by `cam_init` in `src/control/cam_comp.F90`. + subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out) + type(runtime_options), intent(in) :: cam_runtime_opts + type(dyn_import_t), intent(in) :: dyn_in + type(dyn_export_t), intent(in) :: dyn_out + end subroutine stepon_init + + ! Called by `cam_timestep_init` in `src/control/cam_comp.F90`. + subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + real(kind_r8), intent(out) :: dtime_phys + type(runtime_options), intent(in) :: cam_runtime_opts + type(physics_state), intent(inout) :: phys_state + type(physics_tend), intent(inout) :: phys_tend + type(dyn_import_t), intent(inout) :: dyn_in + type(dyn_export_t), intent(inout) :: dyn_out + end subroutine stepon_timestep_init + + ! Called by `cam_run2` in `src/control/cam_comp.F90`. + subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + type(runtime_options), intent(in) :: cam_runtime_opts + type(physics_state), intent(inout) :: phys_state + type(physics_tend), intent(inout) :: phys_tend + type(dyn_import_t), intent(inout) :: dyn_in + type(dyn_export_t), intent(inout) :: dyn_out + end subroutine stepon_run2 + + ! Called by `cam_run3` in `src/control/cam_comp.F90`. + subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out) + real(kind_r8), intent(in) :: dtime_phys + type(runtime_options), intent(in) :: cam_runtime_opts + type(cam_out_t), intent(inout) :: cam_out + type(physics_state), intent(inout) :: phys_state + type(dyn_import_t), intent(inout) :: dyn_in + type(dyn_export_t), intent(inout) :: dyn_out + end subroutine stepon_run3 + + ! Called by `cam_final` in `src/control/cam_comp.F90`. + subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) + type(runtime_options), intent(in) :: cam_runtime_opts + type(dyn_import_t), intent(inout) :: dyn_in + type(dyn_export_t), intent(inout) :: dyn_out + end subroutine stepon_final end module stepon From be3cef986250802e2bc433b5d495e51d8a4499ea Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 11:55:48 -0600 Subject: [PATCH 02/22] Implement MPAS subdriver MPAS dynamical core can now integrate the dynamical states with time. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 121 +++++++++++++++++- 1 file changed, 116 insertions(+), 5 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index f6bb502e..72643196 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -21,7 +21,7 @@ module dyn_mpas_subdriver pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr ! Modules from MPAS. - use atm_core, only: atm_mpas_init_block + use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep, atm_mpas_init_block use atm_core_interface, only: atm_setup_core, atm_setup_domain use atm_time_integration, only: mpas_atm_dynamics_init use mpas_atm_dimensions, only: mpas_atm_set_dims @@ -34,7 +34,7 @@ module dyn_mpas_subdriver mpas_pool_type, mpas_pool_field_info_type, & mpas_pool_character, mpas_pool_real, mpas_pool_integer, & mpas_stream_type, mpas_stream_noerr, & - mpas_time_type, mpas_start_time, & + mpas_time_type, mpas_timeinterval_type, mpas_now, mpas_start_time, & mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & field0dchar, field1dchar, & field0dinteger, field1dinteger, field2dinteger, field3dinteger, & @@ -51,11 +51,13 @@ module dyn_mpas_subdriver mpas_pool_get_array, & mpas_pool_add_dimension, mpas_pool_get_dimension, & mpas_pool_get_field, mpas_pool_get_field_info, & - mpas_pool_initialize_time_levels + mpas_pool_initialize_time_levels, mpas_pool_shift_time_levels use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex use mpas_string_utils, only: mpas_string_replace - use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time + use mpas_timekeeping, only: mpas_advance_clock, mpas_get_clock_time, mpas_get_time, & + mpas_set_timeinterval, & + operator(+), operator(<) use mpas_vector_operations, only: mpas_initialize_vectors implicit none @@ -106,7 +108,8 @@ end subroutine model_error_if logical, allocatable :: is_water_species(:) ! Initialized by `dyn_mpas_init_phase4`. - integer :: coupling_time_interval + integer :: coupling_time_interval = 0 + integer :: number_of_time_steps = 0 contains private @@ -123,6 +126,7 @@ end subroutine model_error_if procedure, pass, public :: compute_unit_vector => dyn_mpas_compute_unit_vector procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4 + procedure, pass, public :: run => dyn_mpas_run ! Accessor subroutines for users to access internal states of MPAS dynamical core. @@ -2673,6 +2677,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) end if self % coupling_time_interval = coupling_time_interval + self % number_of_time_steps = 0 call self % debug_print('Coupling time interval is ' // stringify([real(self % coupling_time_interval, rkind)]) // & ' seconds') @@ -2867,6 +2872,112 @@ pure function almost_equal(a, b, absolute_tolerance, relative_tolerance) end function almost_equal end subroutine dyn_mpas_init_phase4 + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_run + ! + !> \brief Integrates the dynamical states with time + !> \author Michael Duda + !> \date 29 February 2020 + !> \details + !> This subroutine calls MPAS dynamical solver in a loop, with each iteration + !> of the loop advancing the dynamical states forward by one time step, until + !> the coupling time interval is reached. + !> Essentially, it closely follows what is done in `atm_core_run`, but without + !> any calls to MPAS diagnostics manager or MPAS stream manager. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-06-21) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_run(self) + class(mpas_dynamical_core_type), intent(inout) :: self + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_run' + character(strkind) :: date_time + integer :: ierr + real(rkind), pointer :: config_dt + type(mpas_pool_type), pointer :: mpas_pool_diag, mpas_pool_mesh, mpas_pool_state + type(mpas_time_type) :: mpas_time_end, mpas_time_now ! This derived type is analogous to `ESMF_Time`. + type(mpas_timeinterval_type) :: mpas_time_interval ! This derived type is analogous to `ESMF_TimeInterval`. + + call self % debug_print(subname // ' entered') + + nullify(config_dt) + nullify(mpas_pool_diag, mpas_pool_mesh, mpas_pool_state) + + call self % get_variable_pointer(config_dt, 'cfg', 'config_dt') + call self % get_pool_pointer(mpas_pool_diag, 'diag') + call self % get_pool_pointer(mpas_pool_mesh, 'mesh') + call self % get_pool_pointer(mpas_pool_state, 'state') + + mpas_time_now = mpas_get_clock_time(self % domain_ptr % clock, mpas_now, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__) + end if + + call mpas_get_time(mpas_time_now, datetimestring=date_time, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__) + end if + + call self % debug_print('Time integration of MPAS dynamical core begins at ' // trim(adjustl(date_time))) + + call mpas_set_timeinterval(mpas_time_interval, s=self % coupling_time_interval, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to set coupling time interval', subname, __LINE__) + end if + + ! The `+` operator is overloaded here. + mpas_time_end = mpas_time_now + mpas_time_interval + + ! Integrate until the coupling time interval is reached. + ! The `<` operator is overloaded here. + do while (mpas_time_now < mpas_time_end) + ! Number of time steps that has been completed in this MPAS dynamical core instance. + self % number_of_time_steps = self % number_of_time_steps + 1 + + ! Advance the dynamical states forward in time by `config_dt` seconds. + ! Current states are in time level 1. Upon exit, time level 2 will contain updated states. + call atm_do_timestep(self % domain_ptr, config_dt, self % number_of_time_steps) + + ! MPAS `state` pool has two time levels. + ! Swap them after advancing a time step. + call mpas_pool_shift_time_levels(mpas_pool_state) + + call mpas_advance_clock(self % domain_ptr % clock, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to advance clock', subname, __LINE__) + end if + + mpas_time_now = mpas_get_clock_time(self % domain_ptr % clock, mpas_now, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__) + end if + + call self % debug_print('Time step ' // stringify([self % number_of_time_steps]) // ' completed') + end do + + call mpas_get_time(mpas_time_now, datetimestring=date_time, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_now"', subname, __LINE__) + end if + + call self % debug_print('Time integration of MPAS dynamical core ends at ' // trim(adjustl(date_time))) + + ! Compute diagnostic variables like `pressure`, `rho` and `theta` by calling upstream MPAS functionality. + call atm_compute_output_diagnostics(mpas_pool_state, 1, mpas_pool_diag, mpas_pool_mesh) + + nullify(config_dt) + nullify(mpas_pool_diag, mpas_pool_mesh, mpas_pool_state) + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_run + !------------------------------------------------------------------------------- ! function dyn_mpas_get_constituent_name ! From 6cfeab03aa8764031a02849c34776cfc450dffff Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 10 Sep 2024 14:17:25 -0600 Subject: [PATCH 03/22] Support for computing edge wind tendency in MPAS subdriver This functionality is intended for use by dynamics-physics coupling. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 53 +++++++++++++++---- src/dynamics/mpas/dyn_comp.F90 | 2 +- 2 files changed, 44 insertions(+), 11 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 72643196..eb1ad3c4 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -42,6 +42,7 @@ module dyn_mpas_subdriver use mpas_dmpar, only: mpas_dmpar_exch_halo_field, & mpas_dmpar_max_int, mpas_dmpar_sum_int use mpas_domain_routines, only: mpas_allocate_domain + use mpas_field_routines, only: mpas_allocate_scratch_field use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2 use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, & mpas_readstream, mpas_writestream, mpas_writestreamatt @@ -2529,19 +2530,21 @@ end subroutine dyn_mpas_compute_unit_vector !> \date 16 January 2020 !> \details !> This subroutine computes the edge-normal wind vectors at edge points - !> (i.e., `u` in MPAS `state` pool) from wind components at cell points + !> (i.e., `u` in MPAS `state` pool) from the wind components at cell points !> (i.e., `uReconstruct{Zonal,Meridional}` in MPAS `diag` pool). In MPAS, the !> former are PROGNOSTIC variables, while the latter are DIAGNOSTIC variables !> that are "reconstructed" from the former. This subroutine is essentially the !> inverse function of that reconstruction. The purpose is to provide an !> alternative way for MPAS to initialize from zonal and meridional wind - !> components at cell points. + !> components at cell points. If `wind_tendency` is `.true.`, this subroutine + !> operates on the wind tendency due to physics instead. !> \addenda !> Ported and refactored for CAM-SIMA. (KCW, 2024-05-08) ! !------------------------------------------------------------------------------- - subroutine dyn_mpas_compute_edge_wind(self) + subroutine dyn_mpas_compute_edge_wind(self, wind_tendency) class(mpas_dynamical_core_type), intent(in) :: self + logical, intent(in) :: wind_tendency character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_edge_wind' integer :: cell1, cell2, i @@ -2565,14 +2568,24 @@ subroutine dyn_mpas_compute_edge_wind(self) nullify(uedge) ! Make sure halo layers are up-to-date before computation. - call self % exchange_halo('uReconstructZonal') - call self % exchange_halo('uReconstructMeridional') + if (wind_tendency) then + call self % exchange_halo('tend_uzonal') + call self % exchange_halo('tend_umerid') + else + call self % exchange_halo('uReconstructZonal') + call self % exchange_halo('uReconstructMeridional') + end if ! Input. call self % get_variable_pointer(nedges, 'dim', 'nEdges') - call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') - call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') + if (wind_tendency) then + call self % get_variable_pointer(ucellzonal, 'tend_physics', 'tend_uzonal') + call self % get_variable_pointer(ucellmeridional, 'tend_physics', 'tend_umerid') + else + call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') + call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') + end if call self % get_variable_pointer(cellsonedge, 'mesh', 'cellsOnEdge') call self % get_variable_pointer(east, 'mesh', 'east') @@ -2580,7 +2593,11 @@ subroutine dyn_mpas_compute_edge_wind(self) call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors') ! Output. - call self % get_variable_pointer(uedge, 'state', 'u', time_level=1) + if (wind_tendency) then + call self % get_variable_pointer(uedge, 'tend_physics', 'tend_ru_physics') + else + call self % get_variable_pointer(uedge, 'state', 'u', time_level=1) + end if do i = 1, nedges cell1 = cellsonedge(1, i) @@ -2611,7 +2628,11 @@ subroutine dyn_mpas_compute_edge_wind(self) nullify(uedge) ! Make sure halo layers are up-to-date after computation. - call self % exchange_halo('u') + if (wind_tendency) then + call self % exchange_halo('tend_ru_physics') + else + call self % exchange_halo('u') + end if call self % debug_print(subname // ' completed') end subroutine dyn_mpas_compute_edge_wind @@ -2644,6 +2665,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) logical, pointer :: config_do_restart real(rkind), pointer :: config_dt type(field0dreal), pointer :: field_0d_real + type(field2dreal), pointer :: field_2d_real type(mpas_pool_type), pointer :: mpas_pool type(mpas_time_type) :: mpas_time @@ -2655,6 +2677,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) nullify(config_do_restart) nullify(config_dt) nullify(field_0d_real) + nullify(field_2d_real) nullify(mpas_pool) if (coupling_time_interval <= 0) then @@ -2819,6 +2842,16 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) ! Prepare dynamics for time integration. call mpas_atm_dynamics_init(self % domain_ptr) + ! Some additional "scratch" fields are needed for interoperability with CAM-SIMA, but they are not initialized by + ! `mpas_atm_dynamics_init`. Initialize them below. + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_uzonal', field_2d_real, timelevel=1) + call mpas_allocate_scratch_field(field_2d_real) + nullify(field_2d_real) + + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_umerid', field_2d_real, timelevel=1) + call mpas_allocate_scratch_field(field_2d_real) + nullify(field_2d_real) + call self % debug_print(subname // ' completed') call self % debug_print('Successful initialization of MPAS dynamical core') @@ -3259,7 +3292,7 @@ subroutine dyn_mpas_get_pool_pointer(self, pool_pointer, pool_name) pool_pointer => self % domain_ptr % configs case ('dim') pool_pointer => self % domain_ptr % blocklist % dimensions - case ('diag', 'mesh', 'state', 'tend') + case ('diag', 'mesh', 'state', 'tend', 'tend_physics') call mpas_pool_get_subpool(self % domain_ptr % blocklist % allstructs, trim(adjustl(pool_name)), pool_pointer) case default call self % model_error('Unsupported pool name "' // trim(adjustl(pool_name)) // '"', subname, __LINE__) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 5c2d8f0b..dcf954a4 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -474,7 +474,7 @@ subroutine set_mpas_state_u() nullify(ucellzonal, ucellmeridional) - call mpas_dynamical_core % compute_edge_wind() + call mpas_dynamical_core % compute_edge_wind(.false.) end subroutine set_mpas_state_u !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). From c3d7968d5d1e67476909d293e7d68704da0ef350 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 30 Jul 2024 15:40:24 -0600 Subject: [PATCH 04/22] Implement `dyn_run` --- src/dynamics/mpas/dyn_comp.F90 | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index dcf954a4..6f9fefc9 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -49,7 +49,7 @@ module dyn_comp public :: dyn_export_t public :: dyn_readnl public :: dyn_init - ! public :: dyn_run + public :: dyn_run ! public :: dyn_final public :: dyn_debug_print @@ -834,9 +834,14 @@ subroutine mark_variable_as_initialized() end do end subroutine mark_variable_as_initialized - ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`. - ! subroutine dyn_run() - ! end subroutine dyn_run + !> Run MPAS dynamical core to integrate the dynamical states with time. + !> (KCW, 2024-07-11) + subroutine dyn_run() + character(*), parameter :: subname = 'dyn_comp::dyn_run' + + ! MPAS dynamical core will run until the coupling time interval is reached. + call mpas_dynamical_core % run() + end subroutine dyn_run ! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`. ! subroutine dyn_final() From 6d805f65947c052c7245e2eedf3dfe843588d5f4 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 30 Jul 2024 15:41:56 -0600 Subject: [PATCH 05/22] Wire up `dyn_run` --- src/dynamics/mpas/stepon.F90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index 7dfdbc23..ab517b47 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,7 +1,7 @@ module stepon ! Modules from CAM-SIMA. use camsrfexch, only: cam_out_t - use dyn_comp, only: dyn_import_t, dyn_export_t + use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run use physics_types, only: physics_state, physics_tend use runtime_obj, only: runtime_options @@ -31,8 +31,8 @@ subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_t type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(inout) :: phys_state type(physics_tend), intent(inout) :: phys_tend - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out + type(dyn_import_t), intent(in) :: dyn_in + type(dyn_export_t), intent(in) :: dyn_out end subroutine stepon_timestep_init ! Called by `cam_run2` in `src/control/cam_comp.F90`. @@ -40,8 +40,8 @@ subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(inout) :: phys_state type(physics_tend), intent(inout) :: phys_tend - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out + type(dyn_import_t), intent(in) :: dyn_in + type(dyn_export_t), intent(in) :: dyn_out end subroutine stepon_run2 ! Called by `cam_run3` in `src/control/cam_comp.F90`. @@ -50,14 +50,16 @@ subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in type(runtime_options), intent(in) :: cam_runtime_opts type(cam_out_t), intent(inout) :: cam_out type(physics_state), intent(inout) :: phys_state - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out + type(dyn_import_t), intent(in) :: dyn_in + type(dyn_export_t), intent(in) :: dyn_out + + call dyn_run() end subroutine stepon_run3 ! Called by `cam_final` in `src/control/cam_comp.F90`. subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) type(runtime_options), intent(in) :: cam_runtime_opts - type(dyn_import_t), intent(inout) :: dyn_in - type(dyn_export_t), intent(inout) :: dyn_out + type(dyn_import_t), intent(in) :: dyn_in + type(dyn_export_t), intent(in) :: dyn_out end subroutine stepon_final end module stepon From dc786e4fbe6d311414b7bbb39dedffefac3ad56b Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 27 Sep 2024 12:45:51 -0600 Subject: [PATCH 06/22] Sort `use` statements and add comments --- src/dynamics/mpas/dyn_comp.F90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 6f9fefc9..09a74cbe 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -28,16 +28,16 @@ module dyn_comp use time_manager, only: get_start_date, get_stop_date, get_step_size, get_run_duration, timemgr_get_calendar_cf use vert_coord, only: pver, pverp - ! Modules from CESM Share. - use shr_file_mod, only: shr_file_getunit - use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8 - use shr_pio_mod, only: shr_pio_getiosys - ! Modules from CCPP. use cam_ccpp_cap, only: cam_constituents_array use ccpp_kinds, only: kind_phys use phys_vars_init_check, only: mark_as_initialized, std_name_len + ! Modules from CESM Share. + use shr_file_mod, only: shr_file_getunit + use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8 + use shr_pio_mod, only: shr_pio_getiosys + ! Modules from external libraries. use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open @@ -357,10 +357,10 @@ subroutine set_analytic_initial_condition() integer, allocatable :: global_grid_index(:) real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :) real(kind_r8), allocatable :: lat_rad(:), lon_rad(:) - real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow CAM-SIMA convention. - real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow MPAS convention. + real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow CAM-SIMA convention. + real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow MPAS convention. call init_shared_variable() @@ -577,6 +577,8 @@ subroutine set_mpas_state_rho_theta() real(kind_r8), allocatable :: qv_mid_col(:) ! Water vapor mixing ratio (kg/kg) at layer midpoints of each column. real(kind_r8), allocatable :: t_mid(:, :) ! Temperature (K) at layer midpoints. real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column. + ! Be advised that it is not virtual temperature. + ! See doi:10.5065/1DFH-6P97 and doi:10.1175/MWR-D-11-00215.1 for details. real(kind_r8), pointer :: rho(:, :) real(kind_r8), pointer :: theta(:, :) real(kind_r8), pointer :: scalars(:, :, :) From bb8750894be9926efbebe039ac29c7fab3195672 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 9 Aug 2024 13:37:38 -0600 Subject: [PATCH 07/22] Factor out variable finalization --- src/dynamics/mpas/dyn_comp.F90 | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 09a74cbe..6582c8ed 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -370,11 +370,7 @@ subroutine set_analytic_initial_condition() call set_mpas_state_rho_theta() call set_mpas_state_rho_base_theta_base() - deallocate(global_grid_index) - deallocate(lat_rad, lon_rad) - deallocate(z_int) - - nullify(zgrid) + call final_shared_variable() contains !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. !> (KCW, 2024-05-13) @@ -434,6 +430,18 @@ subroutine init_shared_variable() end do end subroutine init_shared_variable + !> Finalize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. + !> (KCW, 2024-05-13) + subroutine final_shared_variable() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::final_shared_variable' + + deallocate(global_grid_index) + deallocate(lat_rad, lon_rad) + deallocate(z_int) + + nullify(zgrid) + end subroutine final_shared_variable + !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_u() From 10613a1189b6ba58b2b12876cf0b97854c928aab Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 5 Aug 2024 15:41:47 -0600 Subject: [PATCH 08/22] Implement `reverse` This helper function reverses the order of elements in a 1-d array. --- src/dynamics/mpas/dyn_comp.F90 | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 6582c8ed..d26fe34d 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -53,6 +53,7 @@ module dyn_comp ! public :: dyn_final public :: dyn_debug_print + public :: reverse public :: mpas_dynamical_core public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels public :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max @@ -856,4 +857,22 @@ end subroutine dyn_run ! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`. ! subroutine dyn_final() ! end subroutine dyn_final + + !> Helper function for reversing the order of elements in `array`. + !> (KCW, 2024-07-17) + pure function reverse(array) + real(kind_r8), intent(in) :: array(:) + real(kind_r8) :: reverse(size(array)) + + integer :: n + + n = size(array) + + ! There is nothing to reverse. + if (n == 0) then + return + end if + + reverse(:) = array(n:1:-1) + end function reverse end module dyn_comp From e30582c1a813f5371ac8a7833802e75eb9714cdb Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 5 Aug 2024 15:46:26 -0600 Subject: [PATCH 09/22] Refactor assignments between CAM-SIMA and MPAS in terms of `reverse` This makes the assignments less error-prone and much more intuitive. --- src/dynamics/mpas/dyn_comp.F90 | 48 +++++++++++++++++----------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index d26fe34d..5d594c8d 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -378,7 +378,7 @@ subroutine set_analytic_initial_condition() subroutine init_shared_variable() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variable' integer :: ierr - integer :: k + integer :: i integer, pointer :: indextocellid(:) real(kind_r8), pointer :: lat_deg(:), lon_deg(:) @@ -426,8 +426,8 @@ subroutine init_shared_variable() call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pverp - z_int(:, k) = zgrid(pverp - k + 1, 1:ncells_solve) + do i = 1, ncells_solve + z_int(i, :) = reverse(zgrid(:, i)) end do end subroutine init_shared_variable @@ -448,7 +448,7 @@ end subroutine final_shared_variable subroutine set_mpas_state_u() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u' integer :: ierr - integer :: k + integer :: i real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) call dyn_debug_print('Setting MPAS state "u"') @@ -466,8 +466,8 @@ subroutine set_mpas_state_u() call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, u=buffer_2d_real) ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pver - ucellzonal(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1) + do i = 1, ncells_solve + ucellzonal(:, i) = reverse(buffer_2d_real(i, :)) end do buffer_2d_real(:, :) = 0.0_kind_r8 @@ -475,8 +475,8 @@ subroutine set_mpas_state_u() call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, v=buffer_2d_real) ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pver - ucellmeridional(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1) + do i = 1, ncells_solve + ucellmeridional(:, i) = reverse(buffer_2d_real(i, :)) end do deallocate(buffer_2d_real) @@ -514,7 +514,7 @@ subroutine set_mpas_state_scalars() 'water_vapor_mixing_ratio_wrt_dry_air' character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_scalars' - integer :: i, k + integer :: i, j integer :: ierr integer, allocatable :: constituent_index(:) integer, pointer :: index_qv @@ -540,12 +540,12 @@ subroutine set_mpas_state_scalars() call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, q=buffer_3d_real, & m_cnst=constituent_index) - ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index. - do i = 1, num_advected - ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pver - scalars(i, k, 1:ncells_solve) = & - buffer_3d_real(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i)) + do i = 1, ncells_solve + ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index. + do j = 1, num_advected + ! Vertical index order is reversed between CAM-SIMA and MPAS. + scalars(j, :, i) = & + reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j))) end do end do @@ -617,8 +617,8 @@ subroutine set_mpas_state_rho_theta() call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, t=buffer_2d_real) ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pver - t_mid(k, :) = buffer_2d_real(:, pver - k + 1) + do i = 1, ncells_solve + t_mid(:, i) = reverse(buffer_2d_real(i, :)) end do deallocate(buffer_2d_real) @@ -770,7 +770,7 @@ end subroutine set_analytic_initial_condition !> (KCW, 2024-07-09) subroutine set_default_constituent() character(*), parameter :: subname = 'dyn_comp::set_default_constituent' - integer :: i, k + integer :: i, j real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. @@ -787,12 +787,12 @@ subroutine set_default_constituent() call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) - ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index. - do i = 1, num_advected - ! Vertical index order is reversed between CAM-SIMA and MPAS. - do k = 1, pver - scalars(i, k, 1:ncells_solve) = & - constituents(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i)) + do i = 1, ncells_solve + ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index. + do j = 1, num_advected + ! Vertical index order is reversed between CAM-SIMA and MPAS. + scalars(j, :, i) = & + reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))) end do end do From 759768fa4767d35a06e0970479aa15c3b7ef9365 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 9 Aug 2024 14:01:07 -0600 Subject: [PATCH 10/22] Implement `dyn_exchange_constituent_state` This subroutine provides a streamlined mechanism for exchanging constituent states between CAM-SIMA and MPAS. --- src/dynamics/mpas/dyn_comp.F90 | 154 ++++++++++++++++++++++++++++++++- 1 file changed, 153 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 5d594c8d..4c5006a6 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -9,7 +9,7 @@ module dyn_comp thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, & thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: const_name, const_is_water_species, num_advected, readtrace + use cam_constituents, only: const_name, const_is_dry, const_is_water_species, num_advected, readtrace use cam_control_mod, only: initial_run use cam_field_read, only: cam_read_field use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id @@ -32,6 +32,7 @@ module dyn_comp use cam_ccpp_cap, only: cam_constituents_array use ccpp_kinds, only: kind_phys use phys_vars_init_check, only: mark_as_initialized, std_name_len + use physics_types, only: phys_state ! Modules from CESM Share. use shr_file_mod, only: shr_file_getunit @@ -53,6 +54,7 @@ module dyn_comp ! public :: dyn_final public :: dyn_debug_print + public :: dyn_exchange_constituent_state public :: reverse public :: mpas_dynamical_core public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels @@ -766,6 +768,156 @@ pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0) end function theta_by_poisson_equation end subroutine set_analytic_initial_condition + !> Exchange and/or convert constituent states between CAM-SIMA and MPAS. + !> If `exchange` is `.true.` and `direction` is "e" or "export", set MPAS state `scalars` from physics state `constituents`. + !> If `exchange` is `.true.` and `direction` is "i" or "import", set physics state `constituents` from MPAS state `scalars`. + !> Think of it as "exporting/importing constituent states in CAM-SIMA to/from MPAS". + !> Otherwise, if `exchange` is `.false.`, no exchange is performed at all. + !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratio that has different + !> definitions between CAM-SIMA and MPAS (i.e., dry/moist). + !> Otherwise, if `conversion` is `.false.`, no conversion is performed at all. + !> This subroutine is intentionally designed to have these elaborate controls due to complications in CAM-SIMA. + !> Some procedures in CAM-SIMA expect constituent states to be dry, while the others expect them to be moist. + !> (KCW, 2024-09-26) + subroutine dyn_exchange_constituent_state(direction, exchange, conversion) + character(*), intent(in) :: direction + logical, intent(in) :: exchange + logical, intent(in) :: conversion + + character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_state' + integer :: i, j + integer :: ierr + integer, allocatable :: is_water_species_index(:) + logical, allocatable :: is_conversion_needed(:) + logical, allocatable :: is_water_species(:) + real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. + real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water mixing ratio. + real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. + + select case (trim(adjustl(direction))) + case ('e', 'export') + if (exchange) then + call dyn_debug_print('Setting MPAS state "scalars" from physics state "constituents"') + end if + + if (conversion) then + call dyn_debug_print('Converting MPAS state "scalars"') + end if + case ('i', 'import') + if (exchange) then + call dyn_debug_print('Setting physics state "constituents" from MPAS state "scalars"') + end if + + if (conversion) then + call dyn_debug_print('Converting physics state "constituents"') + end if + case default + call endrun('Unsupported exchange direction "' // trim(adjustl(direction)) // '"', subname, __LINE__) + end select + + nullify(constituents) + nullify(scalars) + + allocate(is_conversion_needed(num_advected), stat=ierr) + call check_allocate(ierr, subname, & + 'is_conversion_needed(num_advected)', & + 'dyn_comp', __LINE__) + + allocate(is_water_species(num_advected), stat=ierr) + call check_allocate(ierr, subname, & + 'is_water_species(num_advected)', & + 'dyn_comp', __LINE__) + + do j = 1, num_advected + ! All constituent mixing ratio in MPAS is dry. + ! Therefore, conversion in between is needed for any constituent mixing ratio that is not dry in CAM-SIMA. + is_conversion_needed(j) = .not. const_is_dry(j) + is_water_species(j) = const_is_water_species(j) + end do + + allocate(is_water_species_index(count(is_water_species)), stat=ierr) + call check_allocate(ierr, subname, & + 'is_water_species_index(count(is_water_species))', & + 'dyn_comp', __LINE__) + + allocate(sigma_all_q(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'sigma_all_q(pver)', & + 'dyn_comp', __LINE__) + + constituents => cam_constituents_array() + + if (.not. associated(constituents)) then + call endrun('Failed to find variable "constituents"', subname, __LINE__) + end if + + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + + if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then + do i = 1, ncells_solve + if (conversion .and. any(is_conversion_needed)) then + ! The summation term of equation 8 in doi:10.1029/2017MS001257. + ! Using equation 7 here is not possible because it requires all constituent mixing ratio to be moist + ! on the RHS of it. There is no such guarantee in CAM-SIMA. + sigma_all_q(:) = phys_state % pdel(i, :) / phys_state % pdeldry(i, :) + end if + + ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index. + do j = 1, num_advected + if (exchange) then + ! Vertical index order is reversed between CAM-SIMA and MPAS. + scalars(j, :, i) = & + reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))) + end if + + if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then + ! Equation 8 in doi:10.1029/2017MS001257. + scalars(j, :, i) = & + scalars(j, :, i) * reverse(sigma_all_q) + end if + end do + end do + else + is_water_species_index(:) = & + pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species) + + do i = 1, ncells_solve + if (conversion .and. any(is_conversion_needed)) then + ! The summation term of equation 8 in doi:10.1029/2017MS001257. + sigma_all_q(:) = 1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1) + end if + + ! `j` is indexing into `constituents`, so it is regarded as constituent index. + do j = 1, num_advected + if (exchange) then + ! Vertical index order is reversed between CAM-SIMA and MPAS. + constituents(i, :, j) = & + reverse(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i)) + end if + + if (conversion .and. is_conversion_needed(j)) then + ! Equation 8 in doi:10.1029/2017MS001257. + constituents(i, :, j) = & + constituents(i, :, j) / reverse(sigma_all_q) + end if + end do + end do + end if + + deallocate(is_conversion_needed) + deallocate(is_water_species) + deallocate(is_water_species_index) + deallocate(sigma_all_q) + + nullify(constituents) + nullify(scalars) + + if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('scalars') + end if + end subroutine dyn_exchange_constituent_state + !> Set default MPAS state `scalars` (i.e., constituents) in accordance with CCPP, which is a component of CAM-SIMA. !> (KCW, 2024-07-09) subroutine set_default_constituent() From 201fc4c628c87d4bec0f56d2afa4b60c4a5cea4e Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 9 Aug 2024 14:06:37 -0600 Subject: [PATCH 11/22] Switch to use `dyn_exchange_constituent_state` Use the new `dyn_exchange_constituent_state` subroutine to perform default initialization for all constituents instead. --- src/dynamics/mpas/dyn_comp.F90 | 41 ++-------------------------------- 1 file changed, 2 insertions(+), 39 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 4c5006a6..eec89488 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -254,9 +254,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Perform default initialization for all constituents. ! Subsequently, they can be overridden depending on the namelist option (below) and ! the actual availability (checked and handled by MPAS). - call dyn_debug_print('Calling set_default_constituent') + call dyn_debug_print('Calling dyn_exchange_constituent_state') - call set_default_constituent() + call dyn_exchange_constituent_state('e', .true., .false.) ! Namelist option that controls if constituents are to be read from the file. if (readtrace) then @@ -918,43 +918,6 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) end if end subroutine dyn_exchange_constituent_state - !> Set default MPAS state `scalars` (i.e., constituents) in accordance with CCPP, which is a component of CAM-SIMA. - !> (KCW, 2024-07-09) - subroutine set_default_constituent() - character(*), parameter :: subname = 'dyn_comp::set_default_constituent' - integer :: i, j - real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. - real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. - - call dyn_debug_print('Setting default MPAS state "scalars"') - - nullify(constituents) - nullify(scalars) - - constituents => cam_constituents_array() - - if (.not. associated(constituents)) then - call endrun('Failed to find variable "constituents"', subname, __LINE__) - end if - - call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) - - do i = 1, ncells_solve - ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index. - do j = 1, num_advected - ! Vertical index order is reversed between CAM-SIMA and MPAS. - scalars(j, :, i) = & - reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))) - end do - end do - - nullify(constituents) - nullify(scalars) - - ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. - call mpas_dynamical_core % exchange_halo('scalars') - end subroutine set_default_constituent - !> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized !> to prevent physics from attempting to read them from a file. These variables are to be exchanged later !> during dynamics-physics coupling. From 936306e295153f0cda127284506a5fd684ead167 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 22 Jul 2024 13:22:29 -0600 Subject: [PATCH 12/22] Implement dynamics-physics coupling and vice versa --- src/dynamics/mpas/dyn_coupling.F90 | 627 +++++++++++++++++++++++++++++ 1 file changed, 627 insertions(+) create mode 100644 src/dynamics/mpas/dyn_coupling.F90 diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 new file mode 100644 index 00000000..2858e9ff --- /dev/null +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -0,0 +1,627 @@ +module dyn_coupling + ! Modules from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_constituents, only: const_is_water_species, const_qmin, num_advected + use cam_thermo, only: cam_thermo_update + use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, reverse, mpas_dynamical_core, & + ncells_solve + use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, & + constant_rd => rair, constant_rv => rh2o + use runtime_obj, only: cam_runtime_opts + use vert_coord, only: pver, pverp + + ! Modules from CCPP. + use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_kinds, only: kind_phys + use geopotential_temp, only: geopotential_temp_run + use physics_types, only: cappav, cpairv, rairv, zvirv, & + dtime_phys, lagrangian_vertical, & + phys_state, phys_tend + use qneg, only: qneg_run + use static_energy, only: update_dry_static_energy_run + + ! Modules from CESM Share. + use shr_kind_mod, only: kind_cx => shr_kind_cx, kind_r8 => shr_kind_r8 + + implicit none + + private + ! Provide APIs required by CAM-SIMA. + public :: dynamics_to_physics_coupling + public :: physics_to_dynamics_coupling +contains + !> Perform one-way coupling from the dynamics output states to the physics input states. + !> The other coupling direction is implemented by its counterpart, `physics_to_dynamics_coupling`. + !> (KCW, 2024-07-31) + subroutine dynamics_to_physics_coupling() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling' + integer :: column_index + integer, allocatable :: is_water_species_index(:) + integer, pointer :: index_qv + ! Variable name suffixes have the following meanings: + ! `*_col`: Variable is of each column. + ! `*_int`: Variable is at layer interfaces. + ! `*_mid`: Variable is at layer midpoints. + real(kind_r8), allocatable :: pd_int_col(:), & ! Dry hydrostatic air pressure (Pa). + pd_mid_col(:), & ! Dry hydrostatic air pressure (Pa). + p_int_col(:), & ! Full hydrostatic air pressure (Pa). + p_mid_col(:), & ! Full hydrostatic air pressure (Pa). + z_int_col(:) ! Geometric height (m). + real(kind_r8), allocatable :: dpd_col(:), & ! Dry air pressure difference (Pa) between layer interfaces. + dp_col(:), & ! Full air pressure difference (Pa) between layer interfaces. + dz_col(:) ! Geometric height difference (m) between layer interfaces. + real(kind_r8), allocatable :: qv_mid_col(:), & ! Water vapor mixing ratio (kg kg-1). + sigma_all_q_mid_col(:) ! Summation of all water mixing ratio (kg kg-1). + real(kind_r8), allocatable :: rhod_mid_col(:), & ! Dry air density (kg m-3). + rho_mid_col(:) ! Full air density (kg m-3). + real(kind_r8), allocatable :: t_mid_col(:), & ! Temperature (K). + tm_mid_col(:), & ! Modified "moist" temperature (K). + tv_mid_col(:) ! Virtual temperature (K). + real(kind_r8), allocatable :: u_mid_col(:), & ! Eastward wind velocity (m s-1). + v_mid_col(:), & ! Northward wind velocity (m s-1). + omega_mid_col(:) ! Vertical wind velocity (Pa s-1). + real(kind_r8), pointer :: exner(:, :) + real(kind_r8), pointer :: rho_zz(:, :) + real(kind_r8), pointer :: scalars(:, :, :) + real(kind_r8), pointer :: theta_m(:, :) + real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :) + real(kind_r8), pointer :: zgrid(:, :) + real(kind_r8), pointer :: zz(:, :) + + call init_shared_variable() + + call dyn_exchange_constituent_state('i', .true., .false.) + + call dyn_debug_print('Setting physics state variables column by column') + + ! Set variables in the `physics_state` derived type column by column. + ! This way, peak memory usage can be reduced. + do column_index = 1, ncells_solve + call update_shared_variable(column_index) + call set_physics_state_column(column_index) + end do + + call set_physics_state_external() + + call final_shared_variable() + contains + !> Initialize variables that are shared and repeatedly used by the `update_shared_variable` and + !> `set_physics_state_column` internal subroutines. + !> (KCW, 2024-07-20) + subroutine init_shared_variable() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variable' + integer :: i + integer :: ierr + logical, allocatable :: is_water_species(:) + + call dyn_debug_print('Preparing for dynamics-physics coupling') + + nullify(index_qv) + nullify(exner) + nullify(rho_zz) + nullify(scalars) + nullify(theta_m) + nullify(ucellzonal, ucellmeridional, w) + nullify(zgrid) + nullify(zz) + + allocate(is_water_species(num_advected), stat=ierr) + call check_allocate(ierr, subname, & + 'is_water_species(num_advected)', & + 'dyn_coupling', __LINE__) + + do i = 1, num_advected + is_water_species(i) = const_is_water_species(i) + end do + + allocate(is_water_species_index(count(is_water_species)), stat=ierr) + call check_allocate(ierr, subname, & + 'is_water_species_index(count(is_water_species))', & + 'dyn_coupling', __LINE__) + + is_water_species_index(:) = & + pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species) + + deallocate(is_water_species) + + allocate(pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp), stat=ierr) + call check_allocate(ierr, subname, & + 'pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp)', & + 'dyn_coupling', __LINE__) + + allocate(dpd_col(pver), dp_col(pver), dz_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'dpd_col(pver), dp_col(pver), dz_col(pver)', & + 'dyn_coupling', __LINE__) + + allocate(qv_mid_col(pver), sigma_all_q_mid_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'qv_mid_col(pver), sigma_all_q_mid_col(pver)', & + 'dyn_coupling', __LINE__) + + allocate(rhod_mid_col(pver), rho_mid_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'rhod_mid_col(pver), rho_mid_col(pver)', & + 'dyn_coupling', __LINE__) + + allocate(t_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 't_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver)', & + 'dyn_coupling', __LINE__) + + allocate(u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver)', & + 'dyn_coupling', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') + call mpas_dynamical_core % get_variable_pointer(exner, 'diag', 'exner') + call mpas_dynamical_core % get_variable_pointer(rho_zz, 'state', 'rho_zz', time_level=1) + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1) + call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') + call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') + call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) + call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') + call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz') + end subroutine init_shared_variable + + !> Finalize variables that are shared and repeatedly used by the `update_shared_variable` and + !> `set_physics_state_column` internal subroutines. + !> (KCW, 2024-07-20) + subroutine final_shared_variable() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::final_shared_variable' + + deallocate(is_water_species_index) + deallocate(pd_int_col, pd_mid_col, p_int_col, p_mid_col, z_int_col) + deallocate(dpd_col, dp_col, dz_col) + deallocate(qv_mid_col, sigma_all_q_mid_col) + deallocate(rhod_mid_col, rho_mid_col) + deallocate(t_mid_col, tm_mid_col, tv_mid_col) + deallocate(u_mid_col, v_mid_col, omega_mid_col) + + nullify(index_qv) + nullify(exner) + nullify(rho_zz) + nullify(scalars) + nullify(theta_m) + nullify(ucellzonal, ucellmeridional, w) + nullify(zgrid) + nullify(zz) + end subroutine final_shared_variable + + !> Update variables for the specific column, indicated by `i`. This subroutine and `set_physics_state_column` + !> should be called in pairs. + !> (KCW, 2024-07-30) + subroutine update_shared_variable(i) + integer, intent(in) :: i + + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variable' + integer :: k + + ! The summation term of equation 5 in doi:10.1029/2017MS001257. + sigma_all_q_mid_col(:) = 1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1) + + ! Compute thermodynamic variables. + + ! By definition. + z_int_col(:) = zgrid(:, i) + dz_col(:) = z_int_col(2:pverp) - z_int_col(1:pver) + qv_mid_col(:) = scalars(index_qv, :, i) + rhod_mid_col(:) = rho_zz(:, i) * zz(:, i) + + ! Equation 5 in doi:10.1029/2017MS001257. + rho_mid_col(:) = rhod_mid_col(:) * sigma_all_q_mid_col(:) + + ! Hydrostatic equation. + dpd_col(:) = -rhod_mid_col(:) * constant_g * dz_col(:) + dp_col(:) = -rho_mid_col(:) * constant_g * dz_col(:) + + ! By definition of Exner function. Also see below. + tm_mid_col(:) = theta_m(:, i) * exner(:, i) + + ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. + ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. + t_mid_col(:) = tm_mid_col(:) / & + (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:)) + + ! Equation 16 in doi:10.1029/2017MS001257. + ! The numerator terms are just `tm_mid_col` here (i.e., modified "moist" temperature). + tv_mid_col(:) = tm_mid_col(:) / sigma_all_q_mid_col(:) + + ! Hydrostatic equation with equation of state plugged in and arranging for pressure. + pd_mid_col(:) = -constant_rd * t_mid_col(:) * dpd_col(:) / (constant_g * dz_col(:)) + p_mid_col(:) = -constant_rd * tv_mid_col(:) * dp_col(:) / (constant_g * dz_col(:)) + + ! By definition. + p_int_col(pverp) = p_mid_col(pver) + 0.5_kind_r8 * dp_col(pver) + + ! Assume no water at top of model. + pd_int_col(pverp) = p_int_col(pverp) + + ! Integrate downward. + do k = pver, 1, -1 + pd_int_col(k) = pd_int_col(k + 1) - dpd_col(k) + p_int_col(k) = p_int_col(k + 1) - dp_col(k) + end do + + ! Compute momentum variables. + + ! By definition. + u_mid_col(:) = ucellzonal(:, i) + v_mid_col(:) = ucellmeridional(:, i) + omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * (w(1:pver, i) + w(2:pverp, i)) + end subroutine update_shared_variable + + !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type. + !> This subroutine and `update_shared_variable` should be called in pairs. + !> (KCW, 2024-07-30) + subroutine set_physics_state_column(i) + integer, intent(in) :: i + + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_column' + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + ! Always call `reverse` when assigning anything to/from the `physics_state` derived type. + + phys_state % u(i, :) = reverse(u_mid_col) + phys_state % v(i, :) = reverse(v_mid_col) + phys_state % omega(i, :) = reverse(omega_mid_col) + + phys_state % psdry(i) = pd_int_col(1) + phys_state % pintdry(i, :) = reverse(pd_int_col) + phys_state % pmiddry(i, :) = reverse(pd_mid_col) + phys_state % pdeldry(i, :) = reverse(-dpd_col) + phys_state % lnpintdry(i, :) = log(phys_state % pintdry(i, :)) + phys_state % lnpmiddry(i, :) = log(phys_state % pmiddry(i, :)) + phys_state % rpdeldry(i, :) = 1.0_kind_r8 / phys_state % pdeldry(i, :) + + phys_state % ps(i) = p_int_col(1) + phys_state % pint(i, :) = reverse(p_int_col) + phys_state % pmid(i, :) = reverse(p_mid_col) + phys_state % pdel(i, :) = reverse(-dp_col) + phys_state % lnpint(i, :) = log(phys_state % pint(i, :)) + phys_state % lnpmid(i, :) = log(phys_state % pmid(i, :)) + phys_state % rpdel(i, :) = 1.0_kind_r8 / phys_state % pdel(i, :) + + phys_state % t(i, :) = reverse(t_mid_col) + + phys_state % phis(i) = constant_g * z_int_col(1) + end subroutine set_physics_state_column + + !> Set variables in the `physics_state` derived type by calling external procedures. + !> (KCW, 2024-07-30) + subroutine set_physics_state_external() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external' + character(kind_cx) :: cerr + integer :: i + integer :: ierr + real(kind_phys), allocatable :: minimum_constituents(:) + real(kind_phys), pointer :: constituents(:, :, :) + type(ccpp_constituent_prop_ptr_t), pointer :: constituent_properties(:) + + call dyn_debug_print('Setting physics state variables externally') + + nullify(constituents) + nullify(constituent_properties) + + allocate(minimum_constituents(num_advected), stat=ierr) + call check_allocate(ierr, subname, & + 'minimum_constituents(num_advected)', & + 'dyn_coupling', __LINE__) + + do i = 1, num_advected + minimum_constituents(i) = const_qmin(i) + end do + + constituents => cam_constituents_array() + + if (.not. associated(constituents)) then + call endrun('Failed to find variable "constituents"', subname, __LINE__) + end if + + constituent_properties => cam_model_const_properties() + + if (.not. associated(constituent_properties)) then + call endrun('Failed to find variable "constituent_properties"', subname, __LINE__) + end if + + ! Update `cappav`, `cpairv`, `rairv`, `zvirv`, etc. as needed by calling `cam_thermo_update`. + ! Note that this subroutine expects constituents to be dry. + call cam_thermo_update( & + constituents, phys_state % t, ncells_solve, cam_runtime_opts % update_thermodynamic_variables()) + + ! This variable name is really misleading. It actually represents the reciprocal of Exner function + ! with respect to surface pressure. This definition is sometimes used for boundary layer work. See + ! the paragraph below equation 1.5.1c in doi:10.1007/978-94-009-3027-8. + ! Also note that `cappav` is updated externally by `cam_thermo_update`. + do i = 1, ncells_solve + phys_state % exner(i, :) = (phys_state % ps(i) / phys_state % pmid(i, :)) ** cappav(i, :) + end do + + ! Note that constituents become moist after this. + call dyn_exchange_constituent_state('i', .false., .true.) + + ! Impose minimum limits on constituents. + call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr) + + if (ierr /= 0) then + call endrun('Failed to impose minimum limits on constituents externally', subname, __LINE__) + end if + + ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints). + ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_update`. + call geopotential_temp_run( & + pver, lagrangian_vertical, pver, 1, pverp, 1, num_advected, & + phys_state % lnpint, phys_state % pint, phys_state % pmid, phys_state % pdel, phys_state % rpdel, phys_state % t, & + constituents(:, :, mpas_dynamical_core % map_constituent_index(index_qv)), constituents, & + constituent_properties, rairv, constant_g, zvirv, phys_state % zi, phys_state % zm, ncells_solve, ierr, cerr) + + if (ierr /= 0) then + call endrun('Failed to set variable "zi" and "zm" externally', subname, __LINE__) + end if + + ! Set `dse` (i.e., dry static energy). + ! Note that `cpairv` is updated externally by `cam_thermo_update`. + call update_dry_static_energy_run( & + pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr) + + if (ierr /= 0) then + call endrun('Failed to set variable "dse" externally', subname, __LINE__) + end if + + deallocate(minimum_constituents) + + nullify(constituents) + nullify(constituent_properties) + end subroutine set_physics_state_external + end subroutine dynamics_to_physics_coupling + + !> Perform one-way coupling from the physics output states to the dynamics input states. + !> The other coupling direction is implemented by its counterpart, `dynamics_to_physics_coupling`. + !> (KCW, 2024-09-20) + subroutine physics_to_dynamics_coupling() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling' + integer, pointer :: index_qv + real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1) + ! before being updated by physics. + real(kind_r8), pointer :: rho_zz(:, :) + real(kind_r8), pointer :: scalars(:, :, :) + real(kind_r8), pointer :: zz(:, :) + + call init_shared_variable() + + call dyn_exchange_constituent_state('e', .true., .true.) + + call set_mpas_physics_tendency_ru() + call set_mpas_physics_tendency_rho() + call set_mpas_physics_tendency_rtheta() + + call final_shared_variable() + contains + !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. + !> (KCW, 2024-09-13) + subroutine init_shared_variable() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variable' + integer :: ierr + + call dyn_debug_print('Preparing for physics-dynamics coupling') + + nullify(index_qv) + nullify(rho_zz) + nullify(scalars) + nullify(zz) + + call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') + call mpas_dynamical_core % get_variable_pointer(rho_zz, 'state', 'rho_zz', time_level=1) + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz') + + allocate(qv_prev(pver, ncells_solve), stat=ierr) + call check_allocate(ierr, subname, & + 'qv_prev(pver, ncells_solve)', & + 'dyn_coupling', __LINE__) + + ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta` + ! needs it. This must be done before calling `dyn_exchange_constituent_state`. + qv_prev(:, :) = scalars(index_qv, :, 1:ncells_solve) + end subroutine init_shared_variable + + !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. + !> (KCW, 2024-09-13) + subroutine final_shared_variable() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::final_shared_variable' + + deallocate(qv_prev) + + nullify(index_qv) + nullify(rho_zz) + nullify(scalars) + nullify(zz) + end subroutine final_shared_variable + + !> Set MPAS physics tendency `tend_ru_physics` (i.e., "coupled" tendency of horizontal velocity at edge interfaces + !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. + !> (KCW, 2024-09-11) + subroutine set_mpas_physics_tendency_ru() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru' + integer :: i + real(kind_r8), pointer :: u_tendency(:, :), v_tendency(:, :) + + call dyn_debug_print('Setting MPAS physics tendency "tend_ru_physics"') + + nullify(u_tendency, v_tendency) + + call mpas_dynamical_core % get_variable_pointer(u_tendency, 'tend_physics', 'tend_uzonal') + call mpas_dynamical_core % get_variable_pointer(v_tendency, 'tend_physics', 'tend_umerid') + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type. + do i = 1, ncells_solve + u_tendency(:, i) = reverse(phys_tend % dudt_total(i, :)) * rho_zz(:, i) + v_tendency(:, i) = reverse(phys_tend % dvdt_total(i, :)) * rho_zz(:, i) + end do + + nullify(u_tendency, v_tendency) + + call mpas_dynamical_core % compute_edge_wind(.true.) + end subroutine set_mpas_physics_tendency_ru + + !> Set MPAS physics tendency `tend_rho_physics` (i.e., "coupled" tendency of dry air density due to physics). + !> In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. + !> (KCW, 2024-09-11) + subroutine set_mpas_physics_tendency_rho() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho' + real(kind_r8), pointer :: rho_tendency(:, :) + + call dyn_debug_print('Setting MPAS physics tendency "tend_rho_physics"') + + nullify(rho_tendency) + + call mpas_dynamical_core % get_variable_pointer(rho_tendency, 'tend_physics', 'tend_rho_physics') + + ! The material derivative of `rho` (i.e., dry air density) is zero for incompressible fluid. + rho_tendency(:, 1:ncells_solve) = 0.0_kind_r8 + + nullify(rho_tendency) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('tend_rho_physics') + end subroutine set_mpas_physics_tendency_rho + + !> Set MPAS physics tendency `tend_rtheta_physics` (i.e., "coupled" tendency of modified "moist" potential temperature + !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. + !> (KCW, 2024-09-19) + subroutine set_mpas_physics_tendency_rtheta() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rtheta' + integer :: i + integer :: ierr + ! Variable name suffixes have the following meanings: + ! `*_col`: Variable is of each column. + ! `*_prev`: Variable is before being updated by physics. + ! `*_curr`: Variable is after being updated by physics. + real(kind_r8), allocatable :: qv_col_prev(:), qv_col_curr(:) ! Water vapor mixing ratio (kg kg-1). + real(kind_r8), allocatable :: rhod_col(:) ! Dry air density (kg m-3). + real(kind_r8), allocatable :: t_col_prev(:), t_col_curr(:) ! Temperature (K). + real(kind_r8), allocatable :: theta_col_prev(:), theta_col_curr(:) ! Potential temperature (K). + real(kind_r8), allocatable :: thetam_col_prev(:), thetam_col_curr(:) ! Modified "moist" potential temperature (K). + real(kind_r8), pointer :: theta_m(:, :) + real(kind_r8), pointer :: theta_m_tendency(:, :) + + call dyn_debug_print('Setting MPAS physics tendency "tend_rtheta_physics"') + + nullify(theta_m) + nullify(theta_m_tendency) + + allocate(qv_col_prev(pver), qv_col_curr(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'qv_col_prev(pver), qv_col_curr(pver)', & + 'dyn_coupling', __LINE__) + + allocate(rhod_col(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'rhod_col(pver)', & + 'dyn_coupling', __LINE__) + + allocate(t_col_prev(pver), t_col_curr(pver), stat=ierr) + call check_allocate(ierr, subname, & + 't_col_prev(pver), t_col_curr(pver)', & + 'dyn_coupling', __LINE__) + + allocate(theta_col_prev(pver), theta_col_curr(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'theta_col_prev(pver), theta_col_curr(pver)', & + 'dyn_coupling', __LINE__) + + allocate(thetam_col_prev(pver), thetam_col_curr(pver), stat=ierr) + call check_allocate(ierr, subname, & + 'thetam_col_prev(pver), thetam_col_curr(pver)', & + 'dyn_coupling', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1) + call mpas_dynamical_core % get_variable_pointer(theta_m_tendency, 'tend_physics', 'tend_rtheta_physics') + + ! Set `theta_m_tendency` column by column. This way, peak memory usage can be reduced. + do i = 1, ncells_solve + qv_col_curr(:) = scalars(index_qv, :, i) + qv_col_prev(:) = qv_prev(:, i) + rhod_col(:) = rho_zz(:, i) * zz(:, i) + + thetam_col_prev(:) = theta_m(:, i) + theta_col_prev(:) = thetam_col_prev(:) / (1.0_kind_r8 + constant_rv / constant_rd * qv_col_prev(:)) + t_col_prev(:) = t_of_theta_rhod_qv(theta_col_prev, rhod_col, qv_col_prev) + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type. + t_col_curr(:) = t_col_prev(:) + reverse(phys_tend % dtdt_total(i, :)) * dtime_phys + theta_col_curr(:) = theta_of_t_rhod_qv(t_col_curr, rhod_col, qv_col_curr) + thetam_col_curr(:) = theta_col_curr(:) * (1.0_kind_r8 + constant_rv / constant_rd * qv_col_curr(:)) + + theta_m_tendency(:, i) = (thetam_col_curr(:) - thetam_col_prev(:)) * rho_zz(:, i) / dtime_phys + end do + + deallocate(qv_col_prev, qv_col_curr) + deallocate(rhod_col) + deallocate(t_col_prev, t_col_curr) + deallocate(theta_col_prev, theta_col_curr) + deallocate(thetam_col_prev, thetam_col_curr) + + nullify(theta_m) + nullify(theta_m_tendency) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('tend_rtheta_physics') + end subroutine set_mpas_physics_tendency_rtheta + + !> Compute temperature `t` as a function of potential temperature `theta`, dry air density `rhod` and water vapor + !> mixing ratio `qv`. The formulation comes from Poisson equation with equation of state plugged in and arranging + !> for temperature. This function is the exact inverse of `theta_of_t_rhod_qv`, which means that: + !> `t == t_of_theta_rhod_qv(theta_of_t_rhod_qv(t, rhod, qv), rhod, qv)`. + !> (KCW, 2024-09-13) + pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t) + real(kind_r8), intent(in) :: theta, rhod, qv + real(kind_r8) :: t + + real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume. + + ! Mayer's relation. + constant_cvd = constant_cpd - constant_rd + + ! Poisson equation with equation of state plugged in and arranging for temperature. For equation of state, + ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that + ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is + ! described herein: + ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. + ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. + t = (theta ** (constant_cpd / constant_cvd)) * & + (((rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv)) / constant_p0) ** & + (constant_rd / constant_cvd)) + end function t_of_theta_rhod_qv + + !> Compute potential temperature `theta` as a function of temperature `t`, dry air density `rhod` and water vapor + !> mixing ratio `qv`. The formulation comes from Poisson equation with equation of state plugged in and arranging + !> for potential temperature. This function is the exact inverse of `t_of_theta_rhod_qv`, which means that: + !> `theta == theta_of_t_rhod_qv(t_of_theta_rhod_qv(theta, rhod, qv), rhod, qv)`. + !> (KCW, 2024-09-13) + pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta) + real(kind_r8), intent(in) :: t, rhod, qv + real(kind_r8) :: theta + + real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume. + + ! Mayer's relation. + constant_cvd = constant_cpd - constant_rd + + ! Poisson equation with equation of state plugged in and arranging for potential temperature. For equation of state, + ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that + ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is + ! described herein: + ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. + ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. + theta = (t ** (constant_cvd / constant_cpd)) * & + ((constant_p0 / (rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv))) ** & + (constant_rd / constant_cpd)) + end function theta_of_t_rhod_qv + end subroutine physics_to_dynamics_coupling +end module dyn_coupling From e048618ce099d4ba585fcb89a839e8e0b0cee11d Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 30 Jul 2024 12:55:46 -0600 Subject: [PATCH 13/22] Wire up dynamics-physics coupling --- src/dynamics/mpas/stepon.F90 | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index ab517b47..d089636e 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -2,11 +2,13 @@ module stepon ! Modules from CAM-SIMA. use camsrfexch, only: cam_out_t use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run + use dyn_coupling, only: dynamics_to_physics_coupling, physics_to_dynamics_coupling use physics_types, only: physics_state, physics_tend use runtime_obj, only: runtime_options + use time_manager, only: get_step_size - ! Modules from CESM Share. - use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Modules from CCPP. + use ccpp_kinds, only: kind_phys implicit none @@ -27,29 +29,36 @@ end subroutine stepon_init ! Called by `cam_timestep_init` in `src/control/cam_comp.F90`. subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) - real(kind_r8), intent(out) :: dtime_phys + real(kind_phys), intent(out) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts - type(physics_state), intent(inout) :: phys_state - type(physics_tend), intent(inout) :: phys_tend + type(physics_state), intent(in) :: phys_state + type(physics_tend), intent(in) :: phys_tend type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out + + ! Set timestep for physics. + dtime_phys = real(get_step_size(), kind_phys) + + call dynamics_to_physics_coupling() end subroutine stepon_timestep_init ! Called by `cam_run2` in `src/control/cam_comp.F90`. subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) type(runtime_options), intent(in) :: cam_runtime_opts - type(physics_state), intent(inout) :: phys_state - type(physics_tend), intent(inout) :: phys_tend + type(physics_state), intent(in) :: phys_state + type(physics_tend), intent(in) :: phys_tend type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out + + call physics_to_dynamics_coupling() end subroutine stepon_run2 ! Called by `cam_run3` in `src/control/cam_comp.F90`. subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out) - real(kind_r8), intent(in) :: dtime_phys + real(kind_phys), intent(in) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts - type(cam_out_t), intent(inout) :: cam_out - type(physics_state), intent(inout) :: phys_state + type(cam_out_t), intent(in) :: cam_out + type(physics_state), intent(in) :: phys_state type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out From f557372c84a190ce5e759abda13be24479ced37c Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:08:46 -0700 Subject: [PATCH 14/22] Add more detailed comments about an MPAS subroutine --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index eb1ad3c4..ca8e4159 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -3002,7 +3002,8 @@ subroutine dyn_mpas_run(self) call self % debug_print('Time integration of MPAS dynamical core ends at ' // trim(adjustl(date_time))) - ! Compute diagnostic variables like `pressure`, `rho` and `theta` by calling upstream MPAS functionality. + ! Compute diagnostic variables like `pressure`, `rho` and `theta` from time level 1 of MPAS `state` pool + ! by calling upstream MPAS functionality. call atm_compute_output_diagnostics(mpas_pool_state, 1, mpas_pool_diag, mpas_pool_mesh) nullify(config_dt) From b66368d5dfe003b59c0b54f28635ca1029026fdd Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:16:05 -0700 Subject: [PATCH 15/22] Fix grammar in comments --- src/dynamics/mpas/dyn_comp.F90 | 6 +++--- src/dynamics/mpas/dyn_coupling.F90 | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index eec89488..3e512ad7 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -791,7 +791,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) logical, allocatable :: is_conversion_needed(:) logical, allocatable :: is_water_species(:) real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. - real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water mixing ratio. + real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water species mixing ratios. real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. select case (trim(adjustl(direction))) @@ -829,8 +829,8 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) 'dyn_comp', __LINE__) do j = 1, num_advected - ! All constituent mixing ratio in MPAS is dry. - ! Therefore, conversion in between is needed for any constituent mixing ratio that is not dry in CAM-SIMA. + ! All constituent mixing ratios in MPAS are dry. + ! Therefore, conversion in between is needed for any constituent mixing ratios that are not dry in CAM-SIMA. is_conversion_needed(j) = .not. const_is_dry(j) is_water_species(j) = const_is_water_species(j) end do diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 2858e9ff..8d54e481 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -52,7 +52,7 @@ subroutine dynamics_to_physics_coupling() dp_col(:), & ! Full air pressure difference (Pa) between layer interfaces. dz_col(:) ! Geometric height difference (m) between layer interfaces. real(kind_r8), allocatable :: qv_mid_col(:), & ! Water vapor mixing ratio (kg kg-1). - sigma_all_q_mid_col(:) ! Summation of all water mixing ratio (kg kg-1). + sigma_all_q_mid_col(:) ! Summation of all water species mixing ratios (kg kg-1). real(kind_r8), allocatable :: rhod_mid_col(:), & ! Dry air density (kg m-3). rho_mid_col(:) ! Full air density (kg m-3). real(kind_r8), allocatable :: t_mid_col(:), & ! Temperature (K). From ee5cedb18d3f9107fb894b2976640047875b23b1 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:20:33 -0700 Subject: [PATCH 16/22] Only call `reverse` once per cell when performing conversion --- src/dynamics/mpas/dyn_comp.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 3e512ad7..6b228ece 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -859,7 +859,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) ! The summation term of equation 8 in doi:10.1029/2017MS001257. ! Using equation 7 here is not possible because it requires all constituent mixing ratio to be moist ! on the RHS of it. There is no such guarantee in CAM-SIMA. - sigma_all_q(:) = phys_state % pdel(i, :) / phys_state % pdeldry(i, :) + sigma_all_q(:) = reverse(phys_state % pdel(i, :) / phys_state % pdeldry(i, :)) end if ! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index. @@ -873,7 +873,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then ! Equation 8 in doi:10.1029/2017MS001257. scalars(j, :, i) = & - scalars(j, :, i) * reverse(sigma_all_q) + scalars(j, :, i) * sigma_all_q(:) end if end do end do @@ -884,7 +884,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) do i = 1, ncells_solve if (conversion .and. any(is_conversion_needed)) then ! The summation term of equation 8 in doi:10.1029/2017MS001257. - sigma_all_q(:) = 1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1) + sigma_all_q(:) = reverse(1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1)) end if ! `j` is indexing into `constituents`, so it is regarded as constituent index. @@ -898,7 +898,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (conversion .and. is_conversion_needed(j)) then ! Equation 8 in doi:10.1029/2017MS001257. constituents(i, :, j) = & - constituents(i, :, j) / reverse(sigma_all_q) + constituents(i, :, j) / sigma_all_q(:) end if end do end do From 122c64f642f3ddb211990a4f5fb6bc93f8617e98 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:28:33 -0700 Subject: [PATCH 17/22] Pass arguments by keywords for better clarity --- src/dynamics/mpas/dyn_comp.F90 | 2 +- src/dynamics/mpas/dyn_coupling.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 6b228ece..c421b0e0 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -256,7 +256,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! the actual availability (checked and handled by MPAS). call dyn_debug_print('Calling dyn_exchange_constituent_state') - call dyn_exchange_constituent_state('e', .true., .false.) + call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.false.) ! Namelist option that controls if constituents are to be read from the file. if (readtrace) then diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 8d54e481..afa24f27 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -71,7 +71,7 @@ subroutine dynamics_to_physics_coupling() call init_shared_variable() - call dyn_exchange_constituent_state('i', .true., .false.) + call dyn_exchange_constituent_state(direction='i', exchange=.true., conversion=.false.) call dyn_debug_print('Setting physics state variables column by column') @@ -341,7 +341,7 @@ subroutine set_physics_state_external() end do ! Note that constituents become moist after this. - call dyn_exchange_constituent_state('i', .false., .true.) + call dyn_exchange_constituent_state(direction='i', exchange=.false., conversion=.true.) ! Impose minimum limits on constituents. call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr) @@ -392,7 +392,7 @@ subroutine physics_to_dynamics_coupling() call init_shared_variable() - call dyn_exchange_constituent_state('e', .true., .true.) + call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.true.) call set_mpas_physics_tendency_ru() call set_mpas_physics_tendency_rho() From 679e9147e9c8f5472f6a4d8635ecc265c90292a8 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:34:09 -0700 Subject: [PATCH 18/22] Rename internal subroutines --- src/dynamics/mpas/dyn_comp.F90 | 16 +++++------ src/dynamics/mpas/dyn_coupling.F90 | 46 +++++++++++++++--------------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index c421b0e0..2800c74b 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -365,7 +365,7 @@ subroutine set_analytic_initial_condition() real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. ! Dimension and vertical index orders follow MPAS convention. - call init_shared_variable() + call init_shared_variables() call set_mpas_state_u() call set_mpas_state_w() @@ -373,12 +373,12 @@ subroutine set_analytic_initial_condition() call set_mpas_state_rho_theta() call set_mpas_state_rho_base_theta_base() - call final_shared_variable() + call final_shared_variables() contains !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. !> (KCW, 2024-05-13) - subroutine init_shared_variable() - character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variable' + subroutine init_shared_variables() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variables' integer :: ierr integer :: i integer, pointer :: indextocellid(:) @@ -431,19 +431,19 @@ subroutine init_shared_variable() do i = 1, ncells_solve z_int(i, :) = reverse(zgrid(:, i)) end do - end subroutine init_shared_variable + end subroutine init_shared_variables !> Finalize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. !> (KCW, 2024-05-13) - subroutine final_shared_variable() - character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::final_shared_variable' + subroutine final_shared_variables() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::final_shared_variables' deallocate(global_grid_index) deallocate(lat_rad, lon_rad) deallocate(z_int) nullify(zgrid) - end subroutine final_shared_variable + end subroutine final_shared_variables !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). !> (KCW, 2024-05-13) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index afa24f27..3ebc0818 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -69,7 +69,7 @@ subroutine dynamics_to_physics_coupling() real(kind_r8), pointer :: zgrid(:, :) real(kind_r8), pointer :: zz(:, :) - call init_shared_variable() + call init_shared_variables() call dyn_exchange_constituent_state(direction='i', exchange=.true., conversion=.false.) @@ -78,19 +78,19 @@ subroutine dynamics_to_physics_coupling() ! Set variables in the `physics_state` derived type column by column. ! This way, peak memory usage can be reduced. do column_index = 1, ncells_solve - call update_shared_variable(column_index) + call update_shared_variables(column_index) call set_physics_state_column(column_index) end do call set_physics_state_external() - call final_shared_variable() + call final_shared_variables() contains - !> Initialize variables that are shared and repeatedly used by the `update_shared_variable` and + !> Initialize variables that are shared and repeatedly used by the `update_shared_variables` and !> `set_physics_state_column` internal subroutines. !> (KCW, 2024-07-20) - subroutine init_shared_variable() - character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variable' + subroutine init_shared_variables() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variables' integer :: i integer :: ierr logical, allocatable :: is_water_species(:) @@ -165,13 +165,13 @@ subroutine init_shared_variable() call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz') - end subroutine init_shared_variable + end subroutine init_shared_variables - !> Finalize variables that are shared and repeatedly used by the `update_shared_variable` and + !> Finalize variables that are shared and repeatedly used by the `update_shared_variables` and !> `set_physics_state_column` internal subroutines. !> (KCW, 2024-07-20) - subroutine final_shared_variable() - character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::final_shared_variable' + subroutine final_shared_variables() + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::final_shared_variables' deallocate(is_water_species_index) deallocate(pd_int_col, pd_mid_col, p_int_col, p_mid_col, z_int_col) @@ -189,15 +189,15 @@ subroutine final_shared_variable() nullify(ucellzonal, ucellmeridional, w) nullify(zgrid) nullify(zz) - end subroutine final_shared_variable + end subroutine final_shared_variables !> Update variables for the specific column, indicated by `i`. This subroutine and `set_physics_state_column` !> should be called in pairs. !> (KCW, 2024-07-30) - subroutine update_shared_variable(i) + subroutine update_shared_variables(i) integer, intent(in) :: i - character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variable' + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variables' integer :: k ! The summation term of equation 5 in doi:10.1029/2017MS001257. @@ -252,10 +252,10 @@ subroutine update_shared_variable(i) u_mid_col(:) = ucellzonal(:, i) v_mid_col(:) = ucellmeridional(:, i) omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * (w(1:pver, i) + w(2:pverp, i)) - end subroutine update_shared_variable + end subroutine update_shared_variables !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type. - !> This subroutine and `update_shared_variable` should be called in pairs. + !> This subroutine and `update_shared_variables` should be called in pairs. !> (KCW, 2024-07-30) subroutine set_physics_state_column(i) integer, intent(in) :: i @@ -390,7 +390,7 @@ subroutine physics_to_dynamics_coupling() real(kind_r8), pointer :: scalars(:, :, :) real(kind_r8), pointer :: zz(:, :) - call init_shared_variable() + call init_shared_variables() call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.true.) @@ -398,12 +398,12 @@ subroutine physics_to_dynamics_coupling() call set_mpas_physics_tendency_rho() call set_mpas_physics_tendency_rtheta() - call final_shared_variable() + call final_shared_variables() contains !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. !> (KCW, 2024-09-13) - subroutine init_shared_variable() - character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variable' + subroutine init_shared_variables() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variables' integer :: ierr call dyn_debug_print('Preparing for physics-dynamics coupling') @@ -426,12 +426,12 @@ subroutine init_shared_variable() ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta` ! needs it. This must be done before calling `dyn_exchange_constituent_state`. qv_prev(:, :) = scalars(index_qv, :, 1:ncells_solve) - end subroutine init_shared_variable + end subroutine init_shared_variables !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. !> (KCW, 2024-09-13) - subroutine final_shared_variable() - character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::final_shared_variable' + subroutine final_shared_variables() + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::final_shared_variables' deallocate(qv_prev) @@ -439,7 +439,7 @@ subroutine final_shared_variable() nullify(rho_zz) nullify(scalars) nullify(zz) - end subroutine final_shared_variable + end subroutine final_shared_variables !> Set MPAS physics tendency `tend_ru_physics` (i.e., "coupled" tendency of horizontal velocity at edge interfaces !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. From 8c3725cb062cdba0b0160cdf9d539c8e6c85d47a Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:41:22 -0700 Subject: [PATCH 19/22] Include error codes and messages returned by external procedures --- src/dynamics/mpas/dyn_coupling.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 3ebc0818..d777789c 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -20,6 +20,7 @@ module dyn_coupling phys_state, phys_tend use qneg, only: qneg_run use static_energy, only: update_dry_static_energy_run + use string_utils, only: stringify ! Modules from CESM Share. use shr_kind_mod, only: kind_cx => shr_kind_cx, kind_r8 => shr_kind_r8 @@ -347,7 +348,9 @@ subroutine set_physics_state_external() call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr) if (ierr /= 0) then - call endrun('Failed to impose minimum limits on constituents externally', subname, __LINE__) + call endrun('Failed to impose minimum limits on constituents externally' // new_line('') // & + 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), & + subname, __LINE__) end if ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints). @@ -359,7 +362,9 @@ subroutine set_physics_state_external() constituent_properties, rairv, constant_g, zvirv, phys_state % zi, phys_state % zm, ncells_solve, ierr, cerr) if (ierr /= 0) then - call endrun('Failed to set variable "zi" and "zm" externally', subname, __LINE__) + call endrun('Failed to set variable "zi" and "zm" externally' // new_line('') // & + 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), & + subname, __LINE__) end if ! Set `dse` (i.e., dry static energy). @@ -368,7 +373,9 @@ subroutine set_physics_state_external() pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr) if (ierr /= 0) then - call endrun('Failed to set variable "dse" externally', subname, __LINE__) + call endrun('Failed to set variable "dse" externally' // new_line('') // & + 'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), & + subname, __LINE__) end if deallocate(minimum_constituents) From 54378f0a29d7a3833661fdb229a064a3a19f340e Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:43:22 -0700 Subject: [PATCH 20/22] Just use equation of state --- src/dynamics/mpas/dyn_coupling.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index d777789c..158e0913 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -231,9 +231,9 @@ subroutine update_shared_variables(i) ! The numerator terms are just `tm_mid_col` here (i.e., modified "moist" temperature). tv_mid_col(:) = tm_mid_col(:) / sigma_all_q_mid_col(:) - ! Hydrostatic equation with equation of state plugged in and arranging for pressure. - pd_mid_col(:) = -constant_rd * t_mid_col(:) * dpd_col(:) / (constant_g * dz_col(:)) - p_mid_col(:) = -constant_rd * tv_mid_col(:) * dp_col(:) / (constant_g * dz_col(:)) + ! Equation of state. + pd_mid_col(:) = rhod_mid_col(:) * constant_rd * t_mid_col(:) + p_mid_col(:) = rho_mid_col(:) * constant_rd * tv_mid_col(:) ! By definition. p_int_col(pverp) = p_mid_col(pver) + 0.5_kind_r8 * dp_col(pver) From cf5ce82133fe883f17655415c0b6daaa1d01fe11 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:47:09 -0700 Subject: [PATCH 21/22] Do not assume no water at top of model This assumption does not always hold. --- src/dynamics/mpas/dyn_coupling.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 158e0913..f469dd63 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -236,11 +236,9 @@ subroutine update_shared_variables(i) p_mid_col(:) = rho_mid_col(:) * constant_rd * tv_mid_col(:) ! By definition. + pd_int_col(pverp) = pd_mid_col(pver) + 0.5_kind_r8 * dpd_col(pver) p_int_col(pverp) = p_mid_col(pver) + 0.5_kind_r8 * dp_col(pver) - ! Assume no water at top of model. - pd_int_col(pverp) = p_int_col(pverp) - ! Integrate downward. do k = pver, 1, -1 pd_int_col(k) = pd_int_col(k + 1) - dpd_col(k) From 81331442873166f48a2d569934440afdac43d6fd Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Nov 2024 14:49:56 -0700 Subject: [PATCH 22/22] Add more detailed comments about equation derivation --- src/dynamics/mpas/dyn_coupling.F90 | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index f469dd63..6757158f 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -599,6 +599,15 @@ pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t) ! described herein: ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. + ! + ! In short, solve the below equation set for $T$ in terms of $\theta$, $\rho_d$ and $q_v$: + ! \begin{equation*} + ! \begin{cases} + ! \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\ + ! P &= \rho_d R_d T_m \\ + ! T_m &= T (1 + \frac{R_v}{R_d} q_v) + ! \end{cases} + ! \end{equation*} t = (theta ** (constant_cpd / constant_cvd)) * & (((rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv)) / constant_p0) ** & (constant_rd / constant_cvd)) @@ -624,6 +633,15 @@ pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta) ! described herein: ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. + ! + ! In short, solve the below equation set for $\theta$ in terms of $T$, $\rho_d$ and $q_v$: + ! \begin{equation*} + ! \begin{cases} + ! \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\ + ! P &= \rho_d R_d T_m \\ + ! T_m &= T (1 + \frac{R_v}{R_d} q_v) + ! \end{cases} + ! \end{equation*} theta = (t ** (constant_cvd / constant_cpd)) * & ((constant_p0 / (rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv))) ** & (constant_rd / constant_cpd))