diff --git a/cime_config/buildlib b/cime_config/buildlib index 54210d74..497b4dd5 100755 --- a/cime_config/buildlib +++ b/cime_config/buildlib @@ -185,8 +185,9 @@ def _setup_mpas(case: Case) -> None: shutil.copytree(mpas_dycore_src_root, mpas_dycore_bld_root, copy_function=_copy2_as_needed, dirs_exist_ok=True) shutil.move(os.path.join(mpas_dycore_bld_root, "Makefile"), os.path.join(mpas_dycore_bld_root, "Makefile.CESM")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile")), os.path.join(mpas_dycore_bld_root, "Makefile")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile.in.CESM")), os.path.join(mpas_dycore_bld_root, "Makefile.in.CESM")) + _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "driver", "dyn_mpas_subdriver.F90")), os.path.join(mpas_dycore_bld_root, "driver")) + _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile")), mpas_dycore_bld_root) + _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile.in.CESM")), mpas_dycore_bld_root) def _copy2_as_needed(src: str, dst: str) -> None: """ @@ -199,6 +200,9 @@ def _copy2_as_needed(src: str, dst: str) -> None: # `src` and `dst` are both files but their modification time or size differ. # Example scenario: User modified some existing source code files. shutil.copy2(src, dst) + elif os.path.isdir(dst) and os.path.isfile(os.path.join(dst, os.path.basename(src))): + dst = os.path.join(dst, os.path.basename(src)) + _copy2_as_needed(src, dst) else: if os.path.isdir(dst) or os.path.isdir(os.path.dirname(dst)): # `src` is a new file that does not exist at `dst`. diff --git a/cime_config/cam_config.py b/cime_config/cam_config.py index 2e2458f1..80bef16e 100644 --- a/cime_config/cam_config.py +++ b/cime_config/cam_config.py @@ -375,6 +375,10 @@ def __init__(self, case, case_log): # MPAS dynamical core relies on its upstream build infrastructure for compilation instead of CIME to take advantage of future upstream changes automatically. self.create_config("dyn_src_dirs", dyn_dirs_desc, ["mpas"], valid_list_type="str") + + # Add XML namelist definition file for MPAS dynamical core. + mpas_dyn_nml_path = os.path.normpath(os.path.join(cime_conf_path, os.pardir, "src", "dynamics", "mpas")) + self._add_xml_nml_file(mpas_dyn_nml_path, "namelist_definition_mpas_dycore.xml") elif dycore == "none": # Source code directories self.create_config("dyn_src_dirs", dyn_dirs_desc, ["none"], diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/Makefile.in.CESM index 7b15a7e0..a0bff7ed 100644 --- a/src/dynamics/mpas/Makefile.in.CESM +++ b/src/dynamics/mpas/Makefile.in.CESM @@ -9,6 +9,7 @@ endif # export CP = cp -afv +export LN = ln -fsv export MKDIR = mkdir -pv export RM = rm -frv @@ -61,22 +62,29 @@ all: @echo ' `make libmpas-clean ESM="CESM" LIBROOT="..."`' .PHONY: libmpas-prepare -libmpas-prepare: libmpas-archiver-script.txt libmpas-no-physics libmpas-preview +libmpas-prepare: libmpas-archiver-script.txt libmpas-no-physics libmpas-prefix-namelist-groups libmpas-preview # Combine multiple static libraries into `libmpas.a` via archiver/MRI script. This requires GNU or GNU-like archiver (`ar`) program. libmpas-archiver-script.txt: - @echo "create libmpas.a" > $(@) - @echo "addlib libdycore.a" >> $(@) - @echo "addlib libframework.a" >> $(@) - @echo "addlib libops.a" >> $(@) - @echo "save" >> $(@) - @echo "end" >> $(@) + @echo "create libmpas.a" > $(@) + @echo "addlib libdycore.a" >> $(@) + @echo "addlib libframework.a" >> $(@) + @echo "addlib libops.a" >> $(@) + @echo "addmod dyn_mpas_subdriver.o" >> $(@) + @echo "save" >> $(@) + @echo "end" >> $(@) # Do not use built-in MPAS/WRF physics. .PHONY: libmpas-no-physics libmpas-no-physics: @sed -E -i -e "s/^ *PHYSICS=.+$$/PHYSICS=/g" core_atmosphere/Makefile +# Prefix `mpas_` to MPAS namelist groups to avoid naming collisions. +.PHONY: libmpas-prefix-namelist-groups +libmpas-prefix-namelist-groups: + @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/Registry.xml + @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/diagnostics/Registry_*.xml + .PHONY: libmpas-preview libmpas-preview: @echo "Previewing build options for $(LIBROOT)/libmpas.a:" @@ -106,13 +114,21 @@ $(LIBROOT)/libmpas.a: libmpas.a $(MKDIR) $(LIBROOT) $(CP) $(<) $(@) -libmpas.a: $(AUTOCLEAN_DEPS) dycore externals frame ops +libmpas.a: $(AUTOCLEAN_DEPS) dycore externals frame ops subdrv $(AR) $(ARFLAGS) < libmpas-archiver-script.txt + @find -P . -name "*.mod" -type f -exec $(LN) "{}" . ";" .PHONY: libmpas-clean libmpas-clean: clean - $(RM) $(LIBROOT)/libmpas.a libmpas.a + $(RM) $(LIBROOT)/libmpas.a libmpas.a *.mod *.o .PHONY: externals externals: $(AUTOCLEAN_DEPS) ( cd external; $(MAKE) FC="$(FC)" SFC="$(SFC)" CC="$(CC)" SCC="$(SCC)" FFLAGS="$(FFLAGS)" CFLAGS="$(CFLAGS)" CPP="$(CPP)" NETCDF="$(NETCDF)" CORE="$(CORE)" ezxml-lib ) + +.PHONY: subdrv +subdrv: driver/dyn_mpas_subdriver.o + +%.o: %.F90 dycore frame ops + ( cd $( The "class" of MPAS dynamical core. + !> Important data structures like states of MPAS dynamical core are encapsulated inside this derived type to prevent misuse. + !> Type-bound procedures provide well-defined APIs for CAM-SIMA to interact with MPAS dynamical core. + type :: mpas_dynamical_core_type + private + + logical, public :: debug_output = .false. + + integer :: log_unit = output_unit + integer :: mpi_comm = mpi_comm_null + integer :: mpi_rank = 0 + logical :: mpi_rank_root = .false. + + ! Actual implementation is supplied at runtime. + procedure(model_error_if), nopass, pointer :: model_error => null() + + type(core_type), pointer :: corelist => null() + type(domain_type), pointer :: domain_ptr => null() + contains + private + + procedure, pass, public :: debug_print => dyn_mpas_debug_print + procedure, pass, public :: init_phase1 => dyn_mpas_init_phase1 + procedure, pass, public :: read_namelist => dyn_mpas_read_namelist + procedure, pass, public :: init_phase2 => dyn_mpas_init_phase2 + end type mpas_dynamical_core_type +contains + !> Print a debug message with optionally the value(s) of a variable. + !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead. + !> (KCW, 2024-02-03) + subroutine dyn_mpas_debug_print(self, message, variable, printer) + class(mpas_dynamical_core_type), intent(in) :: self + character(*), intent(in) :: message + class(*), optional, intent(in) :: variable(:) + integer, optional, intent(in) :: printer + + ! Bail out early if debug output is not requested. + if (.not. self % debug_output) then + return + end if + + if (present(printer)) then + if (self % mpi_rank /= printer) then + return + end if + else + if (.not. self % mpi_rank_root) then + return + end if + end if + + if (present(variable)) then + write(self % log_unit, '(a)') 'dyn_mpas_debug_print (' // stringify([self % mpi_rank]) // '): ' // & + message // stringify(variable) + else + write(self % log_unit, '(a)') 'dyn_mpas_debug_print (' // stringify([self % mpi_rank]) // '): ' // & + message + end if + end subroutine dyn_mpas_debug_print + + !> Convert one or more values of any intrinsic data types to a character string for pretty printing. + !> If `value` contains more than one element, the elements will be stringified, delimited by `separator`, then concatenated. + !> If `value` contains exactly one element, the element will be stringified without using `separator`. + !> If `value` contains zero element or is of unsupported data types, an empty character string is produced. + !> If `separator` is not supplied, it defaults to `, ` (i.e., a comma and a space). + !> (KCW, 2024-02-04) + pure function stringify(value, separator) + use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64 + + class(*), intent(in) :: value(:) + character(*), optional, intent(in) :: separator + character(:), allocatable :: stringify + + integer, parameter :: sizelimit = 1024 + + character(:), allocatable :: buffer, delimiter, format + integer :: i, n, offset + + if (present(separator)) then + delimiter = separator + else + delimiter = ', ' + end if + + n = min(size(value), sizelimit) + + if (n == 0) then + stringify = '' + + return + end if + + select type (value) + type is (character(*)) + allocate(character(len(value) * n + len(delimiter) * (n - 1)) :: buffer) + + buffer(:) = '' + offset = 0 + + do i = 1, n + if (len(delimiter) > 0 .and. i > 1) then + buffer(offset + 1:offset + len(delimiter)) = delimiter + offset = offset + len(delimiter) + end if + + if (len_trim(adjustl(value(i))) > 0) then + buffer(offset + 1:offset + len_trim(adjustl(value(i)))) = trim(adjustl(value(i))) + offset = offset + len_trim(adjustl(value(i))) + end if + end do + type is (integer(int32)) + allocate(character(11 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(17 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(ss, ', n, '(i0, :, "', delimiter, '"))' + write(buffer, format) value + type is (integer(int64)) + allocate(character(20 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(17 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(ss, ', n, '(i0, :, "', delimiter, '"))' + write(buffer, format) value + type is (logical) + allocate(character(1 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(13 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(', n, '(l1, :, "', delimiter, '"))' + write(buffer, format) value + type is (real(real32)) + allocate(character(13 * n + len(delimiter) * (n - 1)) :: buffer) + + if (maxval(abs(value)) < 1.0e5_real32) then + allocate(character(20 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(f13.6, :, "', delimiter, '"))' + else + allocate(character(23 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(es13.6e2, :, "', delimiter, '"))' + end if + + write(buffer, format) value + type is (real(real64)) + allocate(character(13 * n + len(delimiter) * (n - 1)) :: buffer) + + if (maxval(abs(value)) < 1.0e5_real64) then + allocate(character(20 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(f13.6, :, "', delimiter, '"))' + else + allocate(character(23 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(es13.6e2, :, "', delimiter, '"))' + end if + + write(buffer, format) value + class default + stringify = '' + + return + end select + + stringify = trim(buffer) + end function stringify + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_init_phase1 + ! + !> \brief Tracks `mpas_init` up to the point of reading namelist + !> \author Michael Duda + !> \date 19 April 2019 + !> \details + !> This subroutine follows the stand-alone MPAS subdriver up to, but not + !> including, the point where namelist is read. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-02-02) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_unit, mpas_log_unit) + class(mpas_dynamical_core_type), intent(inout) :: self + integer, intent(in) :: mpi_comm + procedure(model_error_if) :: model_error_impl + integer, intent(in) :: log_unit + integer, intent(in) :: mpas_log_unit(2) + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase1' + integer :: ierr + + self % mpi_comm = mpi_comm + self % model_error => model_error_impl + + if (self % mpi_comm == mpi_comm_null) then + call self % model_error('Invalid MPI communicator group', subname, __LINE__) + end if + + call mpi_comm_rank(self % mpi_comm, self % mpi_rank, ierr) + + if (ierr /= mpi_success) then + call self % model_error('Invalid MPI communicator group', subname, __LINE__) + end if + + self % mpi_rank_root = (self % mpi_rank == 0) + self % log_unit = log_unit + + call self % debug_print(subname // ' entered') + + call self % debug_print('Allocating core') + + allocate(self % corelist, stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate corelist', subname, __LINE__) + end if + + nullify(self % corelist % next) + + call self % debug_print('Allocating domain') + + allocate(self % corelist % domainlist, stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate corelist % domainlist', subname, __LINE__) + end if + + nullify(self % corelist % domainlist % next) + + self % domain_ptr => self % corelist % domainlist + self % domain_ptr % core => self % corelist + + call mpas_allocate_domain(self % domain_ptr) + + self % domain_ptr % domainid = 0 + + call self % debug_print('Calling mpas_framework_init_phase1') + + ! Initialize MPAS framework with supplied MPI communicator group. + call mpas_framework_init_phase1(self % domain_ptr % dminfo, mpi_comm=self % mpi_comm) + + call self % debug_print('Setting up core') + + call atm_setup_core(self % corelist) + + call self % debug_print('Setting up domain') + + call atm_setup_domain(self % domain_ptr) + + call self % debug_print('Setting up log') + + ! Set up the log manager as early as possible so we can use it for any errors/messages during subsequent + ! initialization steps. + ! + ! We need: + ! 1) `domain_ptr` to be allocated; + ! 2) `dmpar_init` to be completed for accessing `dminfo`; + ! 3) `*_setup_core` to assign the `setup_log` function pointer. + ierr = self % domain_ptr % core % setup_log(self % domain_ptr % loginfo, self % domain_ptr, unitnumbers=mpas_log_unit) + + if (ierr /= 0) then + call self % model_error('Failed to setup log for MPAS', subname, __LINE__) + end if + + ! At this point, we should be ready to read namelist in `dyn_comp::dyn_readnl`. + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_init_phase1 + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_read_namelist + ! + !> \brief Tracks `mpas_init` where namelist is being read + !> \author Kuan-Chih Wang + !> \date 2024-02-09 + !> \details + !> This subroutine calls upstream MPAS functionality for reading its own + !> namelist. After that, override designated namelist variables according to + !> information provided from CAM-SIMA. + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_read_namelist(self, namelist_path, & + cf_calendar, start_date_time, stop_date_time, run_duration, initial_run) + class(mpas_dynamical_core_type), intent(inout) :: self + character(*), intent(in) :: namelist_path, cf_calendar + integer, intent(in) :: start_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. + stop_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. + run_duration(4) ! DD, hh, mm, ss. + logical, intent(in) :: initial_run + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_read_namelist' + character(strkind) :: mpas_calendar + character(strkind), pointer :: config_value_c => null() + integer :: ierr + logical, pointer :: config_value_l => null() + + call self % debug_print(subname // ' entered') + + call self % debug_print('Reading namelist at ', [namelist_path]) + + ! Override namelist filename so that we can rely on upstream MPAS functionality for reading its own namelist. + ! The case of missing namelist groups (i.e., `iostat == iostat_end` or `iostat == iostat_eor`) will be handled gracefully. + ! All namelist variables will have reasonable default values even if they are missing. + self % domain_ptr % namelist_filename = trim(adjustl(namelist_path)) + + ierr = self % domain_ptr % core % setup_namelist( & + self % domain_ptr % configs, self % domain_ptr % namelist_filename, self % domain_ptr % dminfo) + + if (ierr /= 0) then + call self % model_error('Namelist setup failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ! Override designated namelist variables according to information provided from CAM-SIMA. + ! These include runtime settings that cannot be determined beforehand. + + call self % debug_print('Overriding designated namelist variables') + + ! CAM-SIMA seems to follow "NetCDF Climate and Forecast (CF) Metadata Conventions" for calendar names. See + ! CF-1.11, section "4.4.1. Calendar". + ! However, this is not the case for MPAS. Translate calendar names between CF and MPAS. + select case (trim(adjustl(cf_calendar))) + case ('360_day') + mpas_calendar = '360day' + case ('365_day', 'noleap') + mpas_calendar = 'gregorian_noleap' + case ('gregorian', 'standard') + ! `gregorian` is a deprecated alternative name for `standard`. + mpas_calendar = 'gregorian' + case default + call self % model_error('Unsupported calendar type "' // trim(adjustl(cf_calendar)) // '"', & + subname, __LINE__) + end select + + call mpas_pool_get_config(self % domain_ptr % configs, 'config_calendar_type', config_value_c) + config_value_c = trim(adjustl(mpas_calendar)) + + call self % debug_print('config_calendar_type = ', [config_value_c]) + + nullify(config_value_c) + + ! MPAS represents date and time in ISO 8601 format. However, the separator between date and time is `_` + ! instead of standard `T`. + ! Format in `YYYY-MM-DD_hh:mm:ss` is acceptable. + call mpas_pool_get_config(self % domain_ptr % configs, 'config_start_time', config_value_c) + config_value_c = stringify(start_date_time(1:3), '-') // '_' // stringify(start_date_time(4:6), ':') + + call self % debug_print('config_start_time = ', [config_value_c]) + + nullify(config_value_c) + + call mpas_pool_get_config(self % domain_ptr % configs, 'config_stop_time', config_value_c) + config_value_c = stringify(stop_date_time(1:3), '-') // '_' // stringify(stop_date_time(4:6), ':') + + call self % debug_print('config_stop_time = ', [config_value_c]) + + nullify(config_value_c) + + ! Format in `DD_hh:mm:ss` is acceptable. + call mpas_pool_get_config(self % domain_ptr % configs, 'config_run_duration', config_value_c) + config_value_c = stringify([run_duration(1)]) // '_' // stringify(run_duration(2:4), ':') + + call self % debug_print('config_run_duration = ', [config_value_c]) + + nullify(config_value_c) + + ! Reflect current run type to MPAS. + if (initial_run) then + ! Run type is initial run. + call mpas_pool_get_config(self % domain_ptr % configs, 'config_do_restart', config_value_l) + config_value_l = .false. + else + ! Run type is branch or restart run. + call mpas_pool_get_config(self % domain_ptr % configs, 'config_do_restart', config_value_l) + config_value_l = .true. + end if + + call self % debug_print('config_do_restart = ', [config_value_l]) + + nullify(config_value_l) + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_read_namelist + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_init_phase2 + ! + !> \brief Tracks `mpas_init` after namelist has been read + !> \author Michael Duda + !> \date 19 April 2019 + !> \details + !> This subroutine follows the stand-alone MPAS subdriver from the point + !> where we call the second phase of MPAS framework initialization up + !> to the check on the existence of the streams. file. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-02-07) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_init_phase2(self, pio_iosystem) + class(mpas_dynamical_core_type), intent(inout) :: self + type(iosystem_desc_t), pointer, intent(in) :: pio_iosystem + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase2' + integer :: ierr + logical :: pio_iosystem_active + + call self % debug_print(subname // ' entered') + + call self % debug_print('Checking PIO system descriptor') + + if (.not. associated(pio_iosystem)) then + call self % model_error('Invalid PIO system descriptor', subname, __LINE__) + end if + + call pio_iosystem_is_active(pio_iosystem, pio_iosystem_active) + + if (.not. pio_iosystem_active) then + call self % model_error('Invalid PIO system descriptor', subname, __LINE__) + end if + + call self % debug_print('Calling mpas_framework_init_phase2') + + ! Initialize MPAS framework with supplied PIO system descriptor. + call mpas_framework_init_phase2(self % domain_ptr, io_system=pio_iosystem) + + ierr = self % domain_ptr % core % define_packages(self % domain_ptr % packages) + + if (ierr /= 0) then + call self % model_error('Package definition failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ierr = self % domain_ptr % core % setup_packages( & + self % domain_ptr % configs, self % domain_ptr % packages, self % domain_ptr % iocontext) + + if (ierr /= 0) then + call self % model_error('Package setup failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ierr = self % domain_ptr % core % setup_decompositions(self % domain_ptr % decompositions) + + if (ierr /= 0) then + call self % model_error('Decomposition setup failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ierr = self % domain_ptr % core % setup_clock(self % domain_ptr % clock, self % domain_ptr % configs) + + if (ierr /= 0) then + call self % model_error('Clock setup failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ! At this point, we should be ready to set up decompositions, build halos, allocate blocks, etc. + ! in `dyn_grid::model_grid_init`. + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_init_phase2 +end module dyn_mpas_subdriver diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index b9443d9a..d80a528c 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -1,5 +1,22 @@ module dyn_comp + use dyn_mpas_subdriver, only: mpas_dynamical_core_type + + ! Modules from CAM. + use cam_abortutils, only: endrun + use cam_control_mod, only: initial_run + use cam_instance, only: atm_id + use cam_logfile, only: debug_output, debugout_none, iulog use runtime_obj, only: runtime_options + use spmd_utils, only: mpicom + use time_manager, only: get_start_date, get_stop_date, get_run_duration, timemgr_get_calendar_cf + + ! Modules from CESM Share. + use shr_file_mod, only: shr_file_getunit + use shr_kind_mod, only: shr_kind_cs + use shr_pio_mod, only: shr_pio_getiosys + + ! Modules from external libraries. + use pio, only: iosystem_desc_t implicit none @@ -11,32 +28,93 @@ module dyn_comp public :: dyn_init ! public :: dyn_run ! public :: dyn_final + public :: mpas_dynamical_core - type dyn_import_t + type :: dyn_import_t end type dyn_import_t - type dyn_export_t + type :: dyn_export_t end type dyn_export_t + + !> The "instance/object" of MPAS dynamical core. + type(mpas_dynamical_core_type) :: mpas_dynamical_core contains + !> Read MPAS namelist from supplied path. + !> Additionally, perform early initialization of MPAS dynamical core. + !> (KCW, 2024-02-09) + ! + ! Called by `read_namelist` in `src/control/runtime_opts.F90`. + subroutine dyn_readnl(namelist_path) + character(*), intent(in) :: namelist_path + + character(shr_kind_cs) :: cam_calendar + integer :: log_unit(2) + integer :: start_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. + stop_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. + run_duration(4), & ! DD, hh, mm, ss. + sec_since_midnight ! Second(s) since midnight. + type(iosystem_desc_t), pointer :: pio_iosystem => null() + + ! Enable/disable the debug output of MPAS dynamical core according to the debug verbosity level of CAM-SIMA. + mpas_dynamical_core % debug_output = (debug_output > debugout_none) + + ! Get free units for MPAS so it can write its own log files, e.g., `log.atmosphere.0000.{out,err}`. + log_unit(1) = shr_file_getunit() + log_unit(2) = shr_file_getunit() + + ! Initialize MPAS framework with supplied MPI communicator group and log units. + ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. + call mpas_dynamical_core % init_phase1(mpicom, endrun, iulog, log_unit) + + cam_calendar = timemgr_get_calendar_cf() + + call get_start_date(start_date_time(1), start_date_time(2), start_date_time(3), sec_since_midnight) + start_date_time(4:6) = sec_to_hour_min_sec(sec_since_midnight) + + call get_stop_date(stop_date_time(1), stop_date_time(2), stop_date_time(3), sec_since_midnight) + stop_date_time(4:6) = sec_to_hour_min_sec(sec_since_midnight) + + call get_run_duration(run_duration(1), sec_since_midnight) + run_duration(2:4) = sec_to_hour_min_sec(sec_since_midnight) + + ! Read MPAS-related namelist variables from `namelist_path`, e.g., `atm_in`. + ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. + call mpas_dynamical_core % read_namelist(namelist_path, & + cam_calendar, start_date_time, stop_date_time, run_duration, initial_run) + + pio_iosystem => shr_pio_getiosys(atm_id) + + ! Initialize MPAS framework with supplied PIO system descriptor. + ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. + call mpas_dynamical_core % init_phase2(pio_iosystem) -! Called by `read_namelist` in `src/control/runtime_opts.F90`. -subroutine dyn_readnl(namelist_path) - character(*), intent(in) :: namelist_path -end subroutine dyn_readnl + nullify(pio_iosystem) + contains + !> Convert second(s) to hour(s), minute(s), and second(s). + !> (KCW, 2024-02-07) + pure function sec_to_hour_min_sec(sec) result(hour_min_sec) + integer, intent(in) :: sec + integer :: hour_min_sec(3) -! Called by `cam_init` in `src/control/cam_comp.F90`. -subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) - type(runtime_options), intent(in) :: cam_runtime_opts - type(dyn_import_t), intent(out) :: dyn_in - type(dyn_export_t), intent(out) :: dyn_out -end subroutine dyn_init + ! These are all intended to be integer arithmetics. + hour_min_sec(1) = sec / 3600 + hour_min_sec(2) = sec / 60 - hour_min_sec(1) * 60 + hour_min_sec(3) = sec - hour_min_sec(1) * 3600 - hour_min_sec(2) * 60 + end function sec_to_hour_min_sec + end subroutine dyn_readnl -! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`. -! subroutine dyn_run() -! end subroutine dyn_run + ! Called by `cam_init` in `src/control/cam_comp.F90`. + subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) + type(runtime_options), intent(in) :: cam_runtime_opts + type(dyn_import_t), intent(out) :: dyn_in + type(dyn_export_t), intent(out) :: dyn_out + end subroutine dyn_init -! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`. -! subroutine dyn_final() -! end subroutine dyn_final + ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`. + ! subroutine dyn_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() + ! end subroutine dyn_final end module dyn_comp diff --git a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml new file mode 100644 index 00000000..874b6d86 --- /dev/null +++ b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml @@ -0,0 +1,710 @@ + + + + + + char*256 + mpas + mpas_nhyd_model + + Time integration scheme in MPAS. + + Possible values: `SRK3' + Default: SRK3 + + + SRK3 + + + + integer + mpas + mpas_nhyd_model + + Order for RK time integration in MPAS. + + Possible values: 2 or 3 + Default: 2 + + + 2 + + + + real + mpas + mpas_nhyd_model + + Model time step, seconds in MPAS. + + Possible values: Positive real values + Default: 720.0 + + + 720.0 + + + + logical + mpas + mpas_nhyd_model + + Whether to super-cycle scalar transport in MPAS. + + Possible values: Logical values + Default: true + + + true + + + + integer + mpas + mpas_nhyd_model + + Number of acoustic steps per full RK step in MPAS. + + Possible values: Positive, even integer values, typically 2 or 6 depending + on transport splitting + Default: 2 + + + 2 + + + + integer + mpas + mpas_nhyd_model + + When config_split_dynamics_transport = T, the number of RK steps per + transport step in MPAS. + + Possible values: Positive integer values + Default: 3 + + + 3 + + + + real + mpas + mpas_nhyd_model + + $\nabla^2$ eddy viscosity for horizontal diffusion of momentum in MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of momentum in + MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + $\nabla^2$ eddy viscosity for vertical diffusion of momentum in MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + $\nabla^2$ eddy viscosity for horizontal diffusion of theta in MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of theta in MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + $\nabla^2$ eddy viscosity for vertical diffusion of theta in MPAS. + + Possible values: Positive real values + Default: 0.0 + + + 0.0 + + + + char*256 + mpas + mpas_nhyd_model + + Formulation of horizontal mixing in MPAS. + + Possible values: `2d_fixed' or `2d_smagorinsky' + Default: 2d_smagorinsky + + + 2d_smagorinsky + + + + real + mpas + mpas_nhyd_model + + Horizontal length scale, used by the Smagorinsky formulation of horizontal + diffusion and by 3-d divergence damping in MPAS. + + Possible values: Positive real values. A zero value implies that the + length scale is prescribed by the nominalMinDc value in the input file. + Default: 0.0 + + + 0.0 + + + + real + mpas + mpas_nhyd_model + + Scaling coefficient of $\delta x^3$ to obtain $\nabla^4$ diffusion + coefficient in MPAS. + + Possible values: Non-negative real values + Default: 0.05 + + + 0.05 + + + + real + mpas + mpas_nhyd_model + + Scaling factor for the divergent component of $\nabla^4 u$ calculation in + MPAS. + + Possible values: Positive real values + Default: 10.0 + + + 10.0 + + + + integer + mpas + mpas_nhyd_model + + Horizontal advection order for w in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Horizontal advection order for theta in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Horizontal advection order for scalars in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Vertical advection order for normal velocities (u) in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Vertical advection order for w in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Vertical advection order for theta in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + integer + mpas + mpas_nhyd_model + + Vertical advection order for scalars in MPAS. + + Possible values: 2, 3, or 4 + Default: 3 + + + 3 + + + + logical + mpas + mpas_nhyd_model + + Whether to advect scalar fields in MPAS. + + Possible values: .true. or .false. + Default: true + + + true + + + + logical + mpas + mpas_nhyd_model + + Whether to enable positive-definite advection of scalars in MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + + logical + mpas + mpas_nhyd_model + + Whether to enable monotonic limiter in scalar advection in MPAS. + + Possible values: .true. or .false. + Default: true + + + true + + + + real + mpas + mpas_nhyd_model + + Upwinding coefficient in the 3rd order advection scheme in MPAS. + + Possible values: 0 $\leq$ config_coef_3rd_order $\leq$ 1 + Default: 0.25 + + + 0.25 + + + + real + mpas + mpas_nhyd_model + + Dimensionless empirical parameter relating the strain tensor to the eddy + viscosity in the Smagorinsky turbulence model in MPAS. + + Possible values: Real values typically in the range 0.1 to 0.4 + Default: 0.125 + + + 0.125 + + + + logical + mpas + mpas_nhyd_model + + Mix full $\theta$ and $u$ fields, or mix perturbation from intitial state + in MPAS. + + Possible values: .true. or .false. + Default: true + + + true + + + + real + mpas + mpas_nhyd_model + + Off-centering parameter for the vertically implicit acoustic integration + in MPAS. + + Possible values: Positive real values + Default: 0.1 + + + 0.1 + + + + real + mpas + mpas_nhyd_model + + 3-d divergence damping coefficient in MPAS. + + Possible values: Positive real values + Default: 0.1 + + + 0.1 + + + + real + mpas + mpas_nhyd_model + + Amount of upwinding in APVM in MPAS. + + Possible values: 0 $\leq$ config_apvm_upwinding $\leq$ 1 + Default: 0.5 + + + 0.5 + + + + logical + mpas + mpas_nhyd_model + + Scale eddy viscosities with mesh-density function for horizontal diffusion + in MPAS. + + Possible values: .true. or .false. + Default: true + + + true + + + + real + mpas + mpas_nhyd_model + + Coefficient for the divergent component of the Laplacian filter of + momentum in the relaxation zone in MPAS. + + Possible values: Positive real values + Default: 6.0 + + + 6.0 + + + + real + mpas + mpas_damping + + Height MSL to begin w-damping profile in MPAS. + + Possible values: Positive real values + Default: 22000.0 + + + 22000.0 + + + + real + mpas + mpas_damping + + Maximum w-damping coefficient at model top in MPAS. + + Possible values: 0 $\leq$ config_xnutr $\leq$ 1 + Default: 0.2 + + + 0.2 + + + + real + mpas + mpas_damping + + Coefficient for scaling the 2nd-order horizontal mixing in the mpas_cam + absorbing layer in MPAS. + + Possible values: 0 $\leq$ config_mpas_cam_coef $\leq$ 1, standard value is + 0.2 + Default: 0.0 + + + 0.0 + + + + integer + mpas + mpas_damping + + Number of layers in which to apply cam 2nd-order horizontal filter top of + model; viscosity linearly ramps to zero by layer number from the top in + MPAS. + + Possible values: Positive integer values + Default: 4 + + + 4 + + + + logical + mpas + mpas_damping + + Whether to apply Rayleigh damping on horizontal velocity in the top-most + model levels. The number of levels is specified by the + config_number_rayleigh_damp_u_levels option, and the damping timescale is + specified by the config_rayleigh_damp_u_timescale_days option. in MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + + real + mpas + mpas_damping + + Timescale, in days (86400 s), for the Rayleigh damping on horizontal + velocity in the top-most model levels. in MPAS. + + Possible values: Positive real values + Default: 5.0 + + + 5.0 + + + + integer + mpas + mpas_damping + + Number of layers in which to apply Rayleigh damping on horizontal velocity + at top of model; damping linearly ramps to zero by layer number from the + top in MPAS. + + Possible values: Positive integer values + Default: 6 + + + 6 + + + + logical + mpas + mpas_limited_area + + Whether to apply lateral boundary conditions in MPAS. + + Possible values: true or false; this option must be set to true for + limited-area simulations and false for global simulations + Default: false + + + false + + + + char*256 + mpas + mpas_decomposition + + Prefix of graph decomposition file, to be suffixed with the MPI task count + in MPAS. + + Possible values: Any valid filename + Default: x1.40962.graph.info.part. + + + x1.40962.graph.info.part. + + + + logical + mpas + mpas_restart + + Whether this run of the model is to restart from a previous restart file + or not in MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + + logical + mpas + mpas_printout + + Whether to print the global min/max of horizontal normal velocity and + vertical velocity each timestep in MPAS. + + Possible values: .true. or .false. + Default: true + + + true + + + + logical + mpas + mpas_printout + + Whether to print the global min/max of horizontal normal velocity and + vertical velocity each timestep, along with the location in the domain + where those extrema occurred in MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + + logical + mpas + mpas_printout + + Whether to print the global min/max of scalar fields each timestep in + MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + + logical + mpas + mpas_assimilation + + Whether this run is within the JEDI data assimilation framework; used to + add temperature and specific humidity as diagnostics in MPAS. + + Possible values: .true. or .false. + Default: false + + + false + + + diff --git a/src/utils/string_utils.F90 b/src/utils/string_utils.F90 index 30f730e8..21133bd1 100644 --- a/src/utils/string_utils.F90 +++ b/src/utils/string_utils.F90 @@ -14,6 +14,7 @@ module string_utils public :: increment_string ! increments a string public :: last_sig_char ! Position of last significant character in string public :: to_str ! convert integer to left justified string + public :: stringify ! Convert one or more values of any intrinsic data types to a character string for pretty printing ! Private module variables integer, parameter :: lower_to_upper = iachar("A") - iachar("a") @@ -230,6 +231,109 @@ end function last_sig_char end function to_str +!========================================================================================= + + !> Convert one or more values of any intrinsic data types to a character string for pretty printing. + !> If `value` contains more than one element, the elements will be stringified, delimited by `separator`, then concatenated. + !> If `value` contains exactly one element, the element will be stringified without using `separator`. + !> If `value` contains zero element or is of unsupported data types, an empty character string is produced. + !> If `separator` is not supplied, it defaults to `, ` (i.e., a comma and a space). + !> (KCW, 2024-02-04) + pure function stringify(value, separator) + use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64 + + class(*), intent(in) :: value(:) + character(*), optional, intent(in) :: separator + character(:), allocatable :: stringify + + integer, parameter :: sizelimit = 1024 + + character(:), allocatable :: buffer, delimiter, format + integer :: i, n, offset + + if (present(separator)) then + delimiter = separator + else + delimiter = ', ' + end if + + n = min(size(value), sizelimit) + + if (n == 0) then + stringify = '' + + return + end if + + select type (value) + type is (character(*)) + allocate(character(len(value) * n + len(delimiter) * (n - 1)) :: buffer) + + buffer(:) = '' + offset = 0 + + do i = 1, n + if (len(delimiter) > 0 .and. i > 1) then + buffer(offset + 1:offset + len(delimiter)) = delimiter + offset = offset + len(delimiter) + end if + + if (len_trim(adjustl(value(i))) > 0) then + buffer(offset + 1:offset + len_trim(adjustl(value(i)))) = trim(adjustl(value(i))) + offset = offset + len_trim(adjustl(value(i))) + end if + end do + type is (integer(int32)) + allocate(character(11 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(17 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(ss, ', n, '(i0, :, "', delimiter, '"))' + write(buffer, format) value + type is (integer(int64)) + allocate(character(20 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(17 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(ss, ', n, '(i0, :, "', delimiter, '"))' + write(buffer, format) value + type is (logical) + allocate(character(1 * n + len(delimiter) * (n - 1)) :: buffer) + allocate(character(13 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + + write(format, '(a, i0, 3a)') '(', n, '(l1, :, "', delimiter, '"))' + write(buffer, format) value + type is (real(real32)) + allocate(character(13 * n + len(delimiter) * (n - 1)) :: buffer) + + if (maxval(abs(value)) < 1.0e5_real32) then + allocate(character(20 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(f13.6, :, "', delimiter, '"))' + else + allocate(character(23 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(es13.6e2, :, "', delimiter, '"))' + end if + + write(buffer, format) value + type is (real(real64)) + allocate(character(13 * n + len(delimiter) * (n - 1)) :: buffer) + + if (maxval(abs(value)) < 1.0e5_real64) then + allocate(character(20 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(f13.6, :, "', delimiter, '"))' + else + allocate(character(23 + len(delimiter) + floor(log10(real(n))) + 1) :: format) + write(format, '(a, i0, 3a)') '(ss, ', n, '(es13.6e2, :, "', delimiter, '"))' + end if + + write(buffer, format) value + class default + stringify = '' + + return + end select + + stringify = trim(buffer) + end function stringify + !========================================================================================= end module string_utils diff --git a/src/utils/time_manager.F90 b/src/utils/time_manager.F90 index ff15e614..91114ccd 100644 --- a/src/utils/time_manager.F90 +++ b/src/utils/time_manager.F90 @@ -35,10 +35,12 @@ module time_manager get_curr_date, &! return date components at end of current timestep get_prev_date, &! return date components at beginning of current timestep get_start_date, &! return components of the start date + get_stop_date, &! return components of the stop date get_ref_date, &! return components of the reference date get_perp_date, &! return components of the perpetual date, and current time of day get_curr_time, &! return components of elapsed time since reference date at end of current timestep get_prev_time, &! return components of elapsed time since reference date at beg of current timestep + get_run_duration, &! return run duration in whole days and remaining seconds get_curr_calday, &! return calendar day at end of current timestep get_calday, &! return calendar day from input date is_first_step, &! return true on first step of initial run @@ -712,6 +714,32 @@ subroutine get_start_date(yr, mon, day, tod) end subroutine get_start_date !========================================================================================= +subroutine get_stop_date(yr, mon, day, tod) + +! Return date components valid at end of run. + +! Arguments + integer, intent(out) ::& + yr, &! year + mon, &! month + day, &! day of month + tod ! time of day (seconds past 0Z) + +! Local variables + character(len=*), parameter :: sub = 'get_stop_date' + integer :: rc + type(ESMF_Time) :: date +!----------------------------------------------------------------------------------------- + + call ESMF_ClockGet(tm_clock, stopTime=date, rc=rc) + call chkrc(rc, sub//': error return from ESMF_ClockGet') + + call ESMF_TimeGet(date, yy=yr, mm=mon, dd=day, s=tod, rc=rc) + call chkrc(rc, sub//': error return from ESMF_TimeGet') + +end subroutine get_stop_date +!========================================================================================= + subroutine get_ref_date(yr, mon, day, tod) ! Return date components of the reference date. @@ -792,11 +820,38 @@ subroutine get_prev_time(days, seconds) call chkrc(rc, sub//': error return from ESMF_ClockGet for refTime') diff = date - ref_date call ESMF_TimeIntervalGet( diff, d=days, s=seconds, rc=rc ) - call chkrc(rc, sub//': error return from ESMF_TimeintervalGet') + call chkrc(rc, sub//': error return from ESMF_TimeIntervalGet') end subroutine get_prev_time !========================================================================================= +subroutine get_run_duration(days, seconds) + +! Return run duration in days and seconds. + +! Arguments + integer, intent(out) ::& + days, &! number of whole days in time interval + seconds ! remaining seconds in time interval + +! Local variables + character(len=*), parameter :: sub = 'get_run_duration' + integer :: rc + type(ESMF_Time) :: start_time, stop_time + type(ESMF_TimeInterval) :: diff +!----------------------------------------------------------------------------------------- + + call ESMF_ClockGet(tm_clock, startTime=start_time, stopTime=stop_time, rc=rc) + call chkrc(rc, sub//': error return from ESMF_ClockGet') + + diff = stop_time - start_time + + call ESMF_TimeIntervalGet(diff, d=days, s=seconds, rc=rc) + call chkrc(rc, sub//': error return from ESMF_TimeIntervalGet') + +end subroutine get_run_duration +!========================================================================================= + function get_curr_calday(offset) ! Return calendar day at end of current timestep with optional offset.