From d51ecac6508072521058bf63d8cd8149568759c5 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Tue, 11 Jul 2023 21:27:51 -0400 Subject: [PATCH 1/9] implement restart for internal tides * add a call to register restart diabatic in MOM * register restart diabatic call register restart internal tides * internal tides restart uses the extra_axes optional argument * support for extra_axes added to MOM_restart/MOM_io --- src/core/MOM.F90 | 11 +- src/framework/MOM_io.F90 | 71 ++++++++---- src/framework/MOM_restart.F90 | 70 ++++++++++-- .../lateral/MOM_internal_tides.F90 | 105 ++++++++++++------ .../vertical/MOM_diabatic_driver.F90 | 36 +++++- 5 files changed, 221 insertions(+), 72 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7b9f2f9d3f..9b1ee12336 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -62,6 +62,7 @@ module MOM use MOM_coord_initialization, only : MOM_initialize_coord, write_vertgrid_file use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end +use MOM_diabatic_driver, only : register_diabatic_restarts use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics @@ -94,6 +95,7 @@ module MOM use MOM_interface_heights, only : find_eta, calc_derived_thermo use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end use MOM_interface_filter, only : interface_filter_CS +use MOM_internal_tides, only : int_tide_CS use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end use MOM_lateral_mixing_coeffs, only : calc_resoln_function, calc_depth_function, VarMix_CS use MOM_MEKE, only : MEKE_alloc_register_restart, step_forward_MEKE @@ -409,6 +411,8 @@ module MOM type(ALE_sponge_CS), pointer :: ALE_sponge_CSp => NULL() !< Pointer to the oda incremental update control structure type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() + !< Pointer to the internal tides control structure + type(int_tide_CS), pointer :: int_tide_CSp => NULL() !< Pointer to the ALE-mode sponge control structure type(ALE_CS), pointer :: ALE_CSp => NULL() !< Pointer to the Arbitrary Lagrangian Eulerian (ALE) vertical coordinate control structure @@ -1911,6 +1915,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(ice_shelf_CS), optional, pointer :: ice_shelf_CSp !< A pointer to an ice shelf control structure type(Wave_parameters_CS), & optional, pointer :: Waves_CSp !< An optional pointer to a wave property CS + ! local variables type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the metric grid use for the run type(ocean_grid_type), pointer :: G_in => NULL() ! Pointer to the input grid @@ -2710,6 +2715,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call stoch_EOS_register_restarts(HI, param_file, CS%stoch_eos_CS, restart_CSp) endif + if ( .not. CS%adiabatic) then + call register_diabatic_restarts(G, US, param_file, CS%int_tide_CSp, restart_CSp) + endif + call callTree_waypoint("restart registration complete (initialize_MOM)") call restart_registry_lock(restart_CSp) @@ -3072,7 +3081,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & else call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, & CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, & - CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp) + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp) endif if (associated(CS%sponge_CSp)) & diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index bebce6f502..b3fbbe648a 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -147,6 +147,20 @@ module MOM_io module procedure read_attribute_int32, read_attribute_int64 end interface read_attribute +!> Type that stores information that can be used to create a non-decomposed axis. +type :: axis_info + character(len=32) :: name = "" !< The name of this axis for use in files + character(len=256) :: longname = "" !< A longer name describing this axis + character(len=48) :: units = "" !< The units of the axis labels + character(len=8) :: cartesian = "N" !< A variable indicating which direction + !! this axis corresponds with. Valid values + !! include 'X', 'Y', 'Z', 'T', and 'N' for none. + integer :: sense = 0 !< This is 1 for axes whose values increase upward, or -1 + !! if they increase downward. The default, 0, is ignored. + integer :: ax_size = 0 !< The number of elements in this axis + real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis [arbitrary] +end type axis_info + !> Type for describing a 3-d variable for output type, public :: vardesc character(len=64) :: name !< Variable name in a NetCDF file @@ -164,22 +178,9 @@ module MOM_io character(len=32) :: dim_names(5) !< The names in the file of the axes for this variable integer :: position = -1 !< An integer encoding the horizontal position, it may !! CENTER, CORNER, EAST_FACE, NORTH_FACE, or 0. + type(axis_info) :: extra_axes(5) end type vardesc -!> Type that stores information that can be used to create a non-decomposed axis. -type :: axis_info ; private - character(len=32) :: name = "" !< The name of this axis for use in files - character(len=256) :: longname = "" !< A longer name describing this axis - character(len=48) :: units = "" !< The units of the axis labels - character(len=8) :: cartesian = "N" !< A variable indicating which direction - !! this axis corresponds with. Valid values - !! include 'X', 'Y', 'Z', 'T', and 'N' for none. - integer :: sense = 0 !< This is 1 for axes whose values increase upward, or -1 - !! if they increase downward. The default, 0, is ignored. - integer :: ax_size = 0 !< The number of elements in this axis - real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis [arbitrary] -end type axis_info - !> Type that stores for a global file attribute type :: attribute_info ; private character(len=:), allocatable :: name !< The name of this attribute @@ -270,7 +271,8 @@ subroutine create_MOM_file(IO_handle, filename, vars, novars, fields, & !! required if the new file uses any !! vertical grid axes. integer(kind=int64), optional, intent(in) :: checksums(:,:) !< checksums of vars - type(axis_info), optional, intent(in) :: extra_axes(:) !< Types with information about + type(axis_info), dimension(:), & + optional, intent(in) :: extra_axes !< Types with information about !! some axes that might be used in this file type(attribute_info), optional, intent(in) :: global_atts(:) !< Global attributes to !! write to this file @@ -1637,7 +1639,8 @@ end subroutine verify_variable_units !! have default values that are empty strings or are appropriate for a 3-d !! tracer field at the tracer cell centers. function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, & - cmor_units, cmor_longname, conversion, caller, position, dim_names, fixed) result(vd) + cmor_units, cmor_longname, conversion, caller, position, dim_names, & + extra_axes, fixed) result(vd) character(len=*), intent(in) :: name !< variable name character(len=*), optional, intent(in) :: units !< variable units character(len=*), optional, intent(in) :: longname !< variable long name @@ -1658,6 +1661,8 @@ function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_na !! NORTH_FACE, and 0 for no horizontal dimensions. character(len=*), dimension(:), & optional, intent(in) :: dim_names !< The names of the dimensions of this variable + type(axis_info), dimension(:), & + optional, intent(in) :: extra_axes logical, optional, intent(in) :: fixed !< If true, this does not evolve with time type(vardesc) :: vd !< vardesc type that is created @@ -1681,7 +1686,8 @@ function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_na call modify_vardesc(vd, units=units, longname=longname, hor_grid=hor_grid, & z_grid=z_grid, t_grid=t_grid, position=position, dim_names=dim_names, & cmor_field_name=cmor_field_name, cmor_units=cmor_units, & - cmor_longname=cmor_longname, conversion=conversion, caller=cllr) + cmor_longname=cmor_longname, conversion=conversion, caller=cllr, & + extra_axes=extra_axes) end function var_desc @@ -1689,7 +1695,8 @@ end function var_desc !> This routine modifies the named elements of a vardesc type. !! All arguments are optional, except the vardesc type to be modified. subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & - cmor_field_name, cmor_units, cmor_longname, conversion, caller, position, dim_names) + cmor_field_name, cmor_units, cmor_longname, conversion, caller, position, dim_names, & + extra_axes) type(vardesc), intent(inout) :: vd !< vardesc type that is modified character(len=*), optional, intent(in) :: name !< name of variable character(len=*), optional, intent(in) :: units !< units of variable @@ -1711,6 +1718,8 @@ subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & !! NORTH_FACE, and 0 for no horizontal dimensions. character(len=*), dimension(:), & optional, intent(in) :: dim_names !< The names of the dimensions of this variable + type(axis_info), dimension(:), & + optional, intent(in) :: extra_axes character(len=120) :: cllr integer :: n @@ -1763,6 +1772,12 @@ subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & endif ; enddo endif + if (present(extra_axes)) then + do n=1,size(extra_axes) ; if (len_trim(extra_axes(n)%name) > 0) then + vd%extra_axes(n) = extra_axes(n) + endif ; enddo + endif + end subroutine modify_vardesc integer function position_from_horgrid(hor_grid) @@ -1906,7 +1921,7 @@ end function cmor_long_std !> This routine queries vardesc subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & cmor_field_name, cmor_units, cmor_longname, conversion, caller, & - position, dim_names) + extra_axes, position, dim_names) type(vardesc), intent(in) :: vd !< vardesc type that is queried character(len=*), optional, intent(out) :: name !< name of variable character(len=*), optional, intent(out) :: units !< units of variable @@ -1921,6 +1936,8 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & !! convert from intensive to extensive !! [various] or [a A-1 ~> 1] character(len=*), optional, intent(in) :: caller !< calling routine? + type(axis_info), dimension(5), & + optional, intent(out) :: extra_axes integer, optional, intent(out) :: position !< A coded integer indicating the horizontal position !! of this variable if it has such dimensions. !! Valid values include CORNER, CENTER, EAST_FACE @@ -1929,7 +1946,8 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & optional, intent(out) :: dim_names !< The names of the dimensions of this variable integer :: n - character(len=120) :: cllr + integer, parameter :: nmax_extraaxes = 5 + character(len=120) :: cllr, varname cllr = "mod_vardesc" if (present(caller)) cllr = trim(caller) @@ -1962,6 +1980,19 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & enddo endif + if (present(extra_axes)) then + ! save_restart expects 5 extra axes (can be empty) + do n=1, nmax_extraaxes + if (vd%extra_axes(n)%ax_size>=1) then + extra_axes(n) = vd%extra_axes(n) + else + ! return an empty axis + write(varname,"('dummy',i1.1)") n + call set_axis_info(extra_axes(n), name=trim(varname), ax_size=1) + endif + enddo + endif + end subroutine query_vardesc diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 75051c32ba..f3284bae24 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -14,6 +14,7 @@ module MOM_restart use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc, get_filename_appendix use MOM_io, only : MULTIPLE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE +use MOM_io, only : axis_info, get_axis_info use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : days_in_month, get_date, set_date @@ -26,6 +27,7 @@ module MOM_restart public restart_registry_lock, restart_init_end, vardesc public restart_files_exist, determine_is_new_run, is_new_run public register_restart_field_as_obsolete, register_restart_pair +public lock_check ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -445,7 +447,7 @@ end subroutine register_restart_pair_ptr4d !> Register a 4-d field for restarts, providing the metadata as individual arguments subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units, conversion, & - hor_grid, z_grid, t_grid) + hor_grid, z_grid, t_grid, extra_axes) real, dimension(:,:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written !! in arbitrary rescaled units [A ~> a] @@ -460,8 +462,22 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent + type(axis_info), dimension(:), & + optional, intent(in) :: extra_axes type(vardesc) :: vd + character(len=32), dimension(:), allocatable :: dim_names + integer :: n, n_extradims + + if (present(extra_axes)) then + n_extradims = size(extra_axes) + allocate( dim_names(n_extradims+2) ) + dim_names(1) = "" + dim_names(2) = "" + do n=3,n_extradims+2 + dim_names(n) = extra_axes(n-2)%name + enddo + endif if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & "register_restart_field_4d: Module must be initialized before "//& @@ -469,8 +485,13 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units call lock_check(CS, name=name) - vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & - z_grid=z_grid, t_grid=t_grid) + if (present(extra_axes)) then + vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid, dim_names=dim_names, extra_axes=extra_axes) + else + vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid) + endif call register_restart_field_ptr4d(f_ptr, vd, mandatory, CS, conversion) @@ -478,7 +499,7 @@ end subroutine register_restart_field_4d !> Register a 3-d field for restarts, providing the metadata as individual arguments subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units, conversion, & - hor_grid, z_grid, t_grid) + hor_grid, z_grid, t_grid, extra_axes) real, dimension(:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written !! in arbitrary rescaled units [A ~> a] @@ -493,8 +514,22 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent + type(axis_info), dimension(:), & + optional, intent(in) :: extra_axes type(vardesc) :: vd + character(len=32), dimension(:), allocatable :: dim_names + integer :: n, n_extradims + + if (present(extra_axes)) then + n_extradims = size(extra_axes) + allocate( dim_names(n_extradims+2) ) + dim_names(1) = "" + dim_names(2) = "" + do n=3,n_extradims+2 + dim_names(n) = extra_axes(n-2)%name + enddo + endif if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & "register_restart_field_3d: Module must be initialized before "//& @@ -502,8 +537,13 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units call lock_check(CS, name=name) - vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & - z_grid=z_grid, t_grid=t_grid) + if (present(extra_axes)) then + vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid, dim_names=dim_names, extra_axes=extra_axes) + else + vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & + z_grid=z_grid, t_grid=t_grid) + endif call register_restart_field_ptr3d(f_ptr, vd, mandatory, CS, conversion) @@ -1309,7 +1349,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ integer :: start_var, next_var ! The starting variables of the ! current and next files. type(MOM_infra_file) :: IO_handle ! The I/O handle of the open fileset - integer :: m, nz + integer :: m, nz, na integer :: num_files ! The number of restart files that will be used. integer :: seconds, days, year, month, hour, minute character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. @@ -1320,9 +1360,13 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ integer(kind=8) :: check_val(CS%max_fields,1) integer :: isL, ieL, jsL, jeL, pos integer :: turns + integer, parameter :: nmax_extradims = 5 + type(axis_info), dimension(:), allocatable :: extra_axes turns = CS%turns + allocate ( extra_axes(nmax_extradims) ) + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "save_restart: Module must be initialized before it is used.") @@ -1361,8 +1405,14 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ do m=start_var,CS%novars call query_vardesc(CS%restart_field(m)%vars, hor_grid=hor_grid, & - z_grid=z_grid, t_grid=t_grid, caller="save_restart") + z_grid=z_grid, t_grid=t_grid, caller="save_restart", & + extra_axes=extra_axes) + var_sz = get_variable_byte_size(hor_grid, z_grid, t_grid, G, nz) + ! factor in size of extra axes, or multiply by 1 + do na=1,nmax_extradims + var_sz = var_sz*extra_axes(na)%ax_size + enddo if ((m==start_var) .OR. (size_in_file < max_file_size-var_sz)) then size_in_file = size_in_file + var_sz @@ -1445,10 +1495,10 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ if (CS%parallel_restartfiles) then call create_MOM_file(IO_handle, trim(restartpath), vars, next_var-start_var, & - fields, MULTIPLE, G=G, GV=GV, checksums=check_val) + fields, MULTIPLE, G=G, GV=GV, checksums=check_val, extra_axes=extra_axes) else call create_MOM_file(IO_handle, trim(restartpath), vars, next_var-start_var, & - fields, SINGLE_FILE, G=G, GV=GV, checksums=check_val) + fields, SINGLE_FILE, G=G, GV=GV, checksums=check_val, extra_axes=extra_axes) endif do m=start_var,next_var-1 diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 8c56107a4f..40a22ba5f6 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -16,8 +16,10 @@ module MOM_internal_tides use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : slasher, MOM_read_data, file_exists +use MOM_io, only : slasher, MOM_read_data, file_exists, axis_info +use MOM_io, only : set_axis_info, get_axis_info use MOM_restart, only : register_restart_field, MOM_restart_CS, restart_init, save_restart +use MOM_restart, only : lock_check, restart_registry_lock use MOM_spatial_means, only : global_area_integral use MOM_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-) use MOM_unit_scaling, only : unit_scale_type @@ -29,12 +31,13 @@ module MOM_internal_tides #include -public propagate_int_tide !, register_int_tide_restarts +public propagate_int_tide, register_int_tide_restarts public internal_tides_init, internal_tides_end public get_lowmode_loss !> This control structure has parameters for the MOM_internal_tides module type, public :: int_tide_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: do_int_tides !< If true, use the internal tide code. integer :: nFreq = 0 !< The number of internal tide frequency bands integer :: nMode = 1 !< The number of internal tide vertical modes @@ -264,6 +267,11 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, & cn(:,:,:) = 0. + ! temporary solution for restart, only works with one freq/1mode + do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,1,1) = CS%En_restart(i,j,a) + enddo ; enddo ; enddo + ! Set properties related to the internal tides, such as the wave speeds, storing some ! of them in the control structure for this module. if (CS%uniform_test_cg > 0.0) then @@ -635,6 +643,11 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, & call post_data(CS%id_En_mode(fr,m), tot_En, CS%diag) endif ; enddo ; enddo + ! save to temporary restart + do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart(i,j,a) = CS%En(i,j,a,1,1) + enddo ; enddo ; enddo + ! Output 3-D (i,j,a) energy density for each frequency and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_En_ang_mode(fr,m) > 0) then call post_data(CS%id_En_ang_mode(fr,m), CS%En(:,:,:,fr,m) , CS%diag) @@ -2258,43 +2271,66 @@ subroutine PPM_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie) enddo ; enddo end subroutine PPM_limit_pos -! subroutine register_int_tide_restarts(G, param_file, CS, restart_CS) -! type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure -! type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters -! type(int_tide_CS), intent(in) :: CS !< Internal tide control structure -! type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure +subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(int_tide_CS), pointer :: CS !< Internal tide control structure + type(MOM_restart_CS), pointer :: restart_CS !< MOM restart control structure -! ! This subroutine is not currently in use!! + ! This subroutine is used to allocate and register any fields in this module + ! that should be written to or read from the restart file. + logical :: use_int_tides + integer :: num_freq, num_angle , num_mode, period_1 + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, i, j, a, fr + character(64) :: var_name, cfr -! ! This subroutine is used to allocate and register any fields in this module -! ! that should be written to or read from the restart file. -! logical :: use_int_tides -! integer :: num_freq, num_angle , num_mode, period_1 -! integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, a -! isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed -! IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + type(axis_info) :: axes_inttides(1) + real, dimension(24) :: angles -! if (associated(CS)) then -! call MOM_error(WARNING, "register_int_tide_restarts called "//& -! "with an associated control structure.") -! return -! endif + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed -! use_int_tides = .false. -! call read_param(param_file, "INTERNAL_TIDES", use_int_tides) -! if (.not.use_int_tides) return + if (associated(CS)) then + call MOM_error(WARNING, "register_int_tide_restarts called "//& + "with an associated control structure.") + return + endif -! allocate(CS) + allocate(CS) -! num_angle = 24 -! call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle) -! allocate(CS%En_restart(isd:ied, jsd:jed, num_angle), source=0.0) + ! write extra axes + call get_param(param_file, "MOM", "INTERNAL_TIDE_ANGLES", num_angle, default=24) + call get_param(param_file, "MOM", "INTERNAL_TIDE_FREQS", num_freq, default=1) + call get_param(param_file, "MOM", "INTERNAL_TIDE_MODES", num_mode, default=1) -! call register_restart_field(CS%En_restart, "En_restart", .false., restart_CS, & -! longname="The internal wave energy density as a function of (i,j,angle,frequency,mode)", & -! units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1') + do a=1,num_angle + angles(a)= a + enddo -! end subroutine register_int_tide_restarts + call set_axis_info(axes_inttides(1), "angle", "", "angle direction", 24, angles, "N", 1) + ! repeat for frequency and vertical modes + + allocate(CS%En_restart(isd:ied, jsd:jed, num_angle), source=0.0) + allocate(CS%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0) + + ! temporary 3d restart + call register_restart_field(CS%En_restart(:,:,:), "IW_energy", .false., restart_CS, & + longname="The internal wave energy density as a function of (i,j,angle,freq,mode)", & + units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + extra_axes=axes_inttides) + ! full 5d restart + !call register_restart_field(CS%En(:,:,:,:,:), "IW_energy", .false., restart_CS, & + ! longname="The internal wave energy density as a function of (i,j,angle,freq,mode)", & + ! units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + ! extra_axes=axes_inttides) + !call restart_registry_lock(restart_CS, unlocked=.false.) + + ! read the restart data back in at init + do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,1,1) = CS%En_restart(i,j,a) + enddo ; enddo ; enddo + +end subroutine register_int_tide_restarts !> This subroutine initializes the internal tides module. subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) @@ -2306,7 +2342,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) !! parameters. type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate !! diagnostic output. - type(int_tide_CS), intent(inout) :: CS !< Internal tide control structure + type(int_tide_CS), pointer :: CS !< Internal tide control structure ! Local variables real :: Angle_size ! size of wedges [rad] @@ -2340,6 +2376,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed nz = GV%ke + CS%initialized = .true. + use_int_tides = .false. call read_param(param_file, "INTERNAL_TIDES", use_int_tides) CS%do_int_tides = use_int_tides @@ -2358,9 +2396,6 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return CS%nFreq = num_freq ; CS%nAngle = num_angle ; CS%nMode = num_mode - ! Allocate energy density array - allocate(CS%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0) - ! Allocate phase speed array allocate(CS%cp(isd:ied, jsd:jed, num_freq, num_mode), source=0.0) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 1bc29ee16f..4274c87b94 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -44,7 +44,7 @@ module MOM_diabatic_driver use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type use MOM_interface_heights, only : find_eta, calc_derived_thermo -use MOM_internal_tides, only : propagate_int_tide +use MOM_internal_tides, only : propagate_int_tide, register_int_tide_restarts use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used use MOM_CVMix_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate @@ -55,6 +55,7 @@ module MOM_diabatic_driver use MOM_opacity, only : absorbRemainingSW, optics_type, optics_nbands use MOM_open_boundary, only : ocean_OBC_type use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS +use MOM_restart, only : MOM_restart_CS use MOM_set_diffusivity, only : set_diffusivity, set_BBL_TKE use MOM_set_diffusivity, only : set_diffusivity_init, set_diffusivity_end use MOM_set_diffusivity, only : set_diffusivity_CS @@ -80,6 +81,7 @@ module MOM_diabatic_driver public extract_diabatic_member public adiabatic public adiabatic_driver_init +public register_diabatic_restarts ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -227,12 +229,13 @@ module MOM_diabatic_driver type(KPP_CS), pointer :: KPP_CSp => NULL() !< Control structure for a child module type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() !< Control structure for a child module type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Control structure for a child module + type(int_tide_CS), pointer :: int_tide_CSp => NULL() !< Control structure for a child module + type(bulkmixedlayer_CS) :: bulkmixedlayer !< Bulk mixed layer control structure type(CVMix_conv_CS) :: CVMix_conv !< CVMix convection control structure type(energetic_PBL_CS) :: ePBL !< Energetic PBL control structure type(entrain_diffusive_CS) :: entrain_diffusive !< Diffusive entrainment control structure type(geothermal_CS) :: geothermal !< Geothermal control structure - type(int_tide_CS) :: int_tide !< Internal tide control structure type(opacity_CS) :: opacity !< Opacity control structure type(regularize_layers_CS) :: regularize_layers !< Regularize layer control structure @@ -387,7 +390,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & CS%int_tide_input_CSp) call propagate_int_tide(h, tv, CS%int_tide_input%TKE_itidal_input, CS%int_tide_input%tideamp, & - CS%int_tide_input%Nb, dt, G, GV, US, CS%int_tide) + CS%int_tide_input%Nb, dt, G, GV, US, CS%int_tide_CSp) if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)") endif ! end CS%use_int_tides @@ -2938,7 +2941,7 @@ end subroutine adiabatic_driver_init !> This routine initializes the diabatic driver module. subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, & ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, & - ALE_sponge_CSp, oda_incupd_CSp) + ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp) type(time_type), target :: Time !< model time type(ocean_grid_type), intent(inout) :: G !< model grid structure type(verticalGrid_type), intent(in) :: GV !< model vertical grid structure @@ -2956,6 +2959,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< pointer to the ALE sponge module control structure type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the ocean data assimilation incremental !! update module control structure + type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tide structure ! Local variables real :: Kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1] @@ -2990,6 +2994,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (associated(sponge_CSp)) CS%sponge_CSp => sponge_CSp if (associated(ALE_sponge_CSp)) CS%ALE_sponge_CSp => ALE_sponge_CSp if (associated(oda_incupd_CSp)) CS%oda_incupd_CSp => oda_incupd_CSp + if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp CS%useALEalgorithm = useALEalgorithm CS%use_bulkmixedlayer = (GV%nkml > 0) @@ -3453,12 +3458,14 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (CS%use_int_tides) then call int_tide_input_init(Time, G, GV, US, param_file, diag, CS%int_tide_input_CSp, & CS%int_tide_input) - call internal_tides_init(Time, G, GV, US, param_file, diag, CS%int_tide) + call internal_tides_init(Time, G, GV, US, param_file, diag, int_tide_CSp) endif + !if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp + physical_OBL_scheme = (CS%use_bulkmixedlayer .or. CS%use_KPP .or. CS%use_energetic_PBL) ! initialize module for setting diffusivities - call set_diffusivity_init(Time, G, GV, US, param_file, diag, CS%set_diff_CSp, CS%int_tide, & + call set_diffusivity_init(Time, G, GV, US, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp, & halo_TS=CS%halo_TS_diff, double_diffuse=CS%double_diffuse, & physical_OBL_scheme=physical_OBL_scheme) @@ -3514,6 +3521,23 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di end subroutine diabatic_driver_init +!> Routine to register restarts, pass-through to children modules +subroutine register_diabatic_restarts(G, US, param_file, int_tide_CSp, restart_CSp) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(int_tide_CS), pointer :: int_tide_CSp !< Internal tide control structure + type(MOM_restart_CS), pointer :: restart_CSp !< MOM restart control structure + + logical :: use_int_tides = .false. + + call read_param(param_file, "INTERNAL_TIDES", use_int_tides) + + if (use_int_tides) then + call register_int_tide_restarts(G, US, param_file, int_tide_CSp, restart_CSp) + endif + +end subroutine register_diabatic_restarts !> Routine to close the diabatic driver module subroutine diabatic_driver_end(CS) From 1e3089749b9538c288e06452d44ce1f0684c6372 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Tue, 22 Aug 2023 08:52:39 -0400 Subject: [PATCH 2/9] implement 4d restart/ parallel_restart not working --- .../lateral/MOM_internal_tides.F90 | 45 ++++++++++--------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 7ce9aa7c9b..acb58fda5b 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -140,7 +140,7 @@ module MOM_internal_tides real, allocatable :: En(:,:,:,:,:) !< The internal wave energy density as a function of (i,j,angle,frequency,mode) !! integrated within an angular and frequency band [R Z3 T-2 ~> J m-2] - real, allocatable :: En_restart(:,:,:) + real, allocatable :: En_restart(:,:,:,:) !< The internal wave energy density as a function of (i,j,angle); temporary for restart real, allocatable, dimension(:) :: frequency !< The frequency of each band [T-1 ~> s-1]. @@ -270,9 +270,9 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, cn(:,:,:) = 0. ! temporary solution for restart, only works with one freq/1mode - do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied - CS%En(i,j,a,1,1) = CS%En_restart(i,j,a) - enddo ; enddo ; enddo + do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) + enddo ; enddo ; enddo ; enddo ! Set properties related to the internal tides, such as the wave speeds, storing some ! of them in the control structure for this module. @@ -655,9 +655,9 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, endif ; enddo ; enddo ! save to temporary restart - do a=1,CS%nAngle ; do j=js,je ; do i=is,ie - CS%En_restart(i,j,a) = CS%En(i,j,a,1,1) - enddo ; enddo ; enddo + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart(i,j,a,fr) = CS%En(i,j,a,fr,1) + enddo ; enddo ; enddo ; enddo ! Output 3-D (i,j,a) energy density for each frequency and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_En_ang_mode(fr,m) > 0) then @@ -2303,8 +2303,8 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, i, j, a, fr character(64) :: var_name, cfr - type(axis_info) :: axes_inttides(1) - real, dimension(24) :: angles + type(axis_info) :: axes_inttides(2) + real, dimension(:), allocatable :: angles, freqs isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2321,19 +2321,22 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) call get_param(param_file, "MOM", "INTERNAL_TIDE_FREQS", num_freq, default=1) call get_param(param_file, "MOM", "INTERNAL_TIDE_MODES", num_mode, default=1) - do a=1,num_angle - angles(a)= a - enddo + allocate ( angles(num_angle) ) + allocate ( freqs(num_freq) ) + + do a=1,num_angle ; angles(a)= a ; enddo + do fr=1,num_freq ; freqs(fr)= fr ; enddo - call set_axis_info(axes_inttides(1), "angle", "", "angle direction", 24, angles, "N", 1) - ! repeat for frequency and vertical modes + call set_axis_info(axes_inttides(1), "angle", "", "angle direction", num_angle, angles, "N", 1) + call set_axis_info(axes_inttides(2), "freq", "", "wave frequency", num_freq, freqs, "N", 1) + ! repeat for vertical modes - allocate(CS%En_restart(isd:ied, jsd:jed, num_angle), source=0.0) + allocate(CS%En_restart(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) allocate(CS%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0) - ! temporary 3d restart - call register_restart_field(CS%En_restart(:,:,:), "IW_energy", .false., restart_CS, & - longname="The internal wave energy density as a function of (i,j,angle,freq,mode)", & + ! temporary 4d restart + call register_restart_field(CS%En_restart(:,:,:,:), "IW_energy", .false., restart_CS, & + longname="The internal wave energy density as a function of (i,j,angle,freq)", & units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & extra_axes=axes_inttides) ! full 5d restart @@ -2344,9 +2347,9 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) !call restart_registry_lock(restart_CS, unlocked=.false.) ! read the restart data back in at init - do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied - CS%En(i,j,a,1,1) = CS%En_restart(i,j,a) - enddo ; enddo ; enddo + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) + enddo ; enddo ; enddo ; enddo end subroutine register_int_tide_restarts From 66c56781c7486e82733c7ac353d18fa9a8034c7c Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Tue, 22 Aug 2023 09:22:01 -0400 Subject: [PATCH 3/9] add proper documentation for extra axes --- src/framework/MOM_io.F90 | 8 ++++---- src/framework/MOM_restart.F90 | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 8b5d525b97..27d244b226 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -179,7 +179,7 @@ module MOM_io character(len=32) :: dim_names(5) !< The names in the file of the axes for this variable integer :: position = -1 !< An integer encoding the horizontal position, it may !! CENTER, CORNER, EAST_FACE, NORTH_FACE, or 0. - type(axis_info) :: extra_axes(5) + type(axis_info) :: extra_axes(5) !< dimensions other than space-time end type vardesc !> Type that stores for a global file attribute @@ -1776,7 +1776,7 @@ function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_na character(len=*), dimension(:), & optional, intent(in) :: dim_names !< The names of the dimensions of this variable type(axis_info), dimension(:), & - optional, intent(in) :: extra_axes + optional, intent(in) :: extra_axes !< dimensions other than space-time logical, optional, intent(in) :: fixed !< If true, this does not evolve with time type(vardesc) :: vd !< vardesc type that is created @@ -1833,7 +1833,7 @@ subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & character(len=*), dimension(:), & optional, intent(in) :: dim_names !< The names of the dimensions of this variable type(axis_info), dimension(:), & - optional, intent(in) :: extra_axes + optional, intent(in) :: extra_axes !< dimensions other than space-time character(len=120) :: cllr integer :: n @@ -2051,7 +2051,7 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & !! [various] or [a A-1 ~> 1] character(len=*), optional, intent(in) :: caller !< calling routine? type(axis_info), dimension(5), & - optional, intent(out) :: extra_axes + optional, intent(out) :: extra_axes !< dimensions other than space-time integer, optional, intent(out) :: position !< A coded integer indicating the horizontal position !! of this variable if it has such dimensions. !! Valid values include CORNER, CENTER, EAST_FACE diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index f3284bae24..24e6da0efd 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -463,7 +463,7 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent type(axis_info), dimension(:), & - optional, intent(in) :: extra_axes + optional, intent(in) :: extra_axes !< dimensions other than space-time type(vardesc) :: vd character(len=32), dimension(:), allocatable :: dim_names @@ -515,7 +515,7 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent type(axis_info), dimension(:), & - optional, intent(in) :: extra_axes + optional, intent(in) :: extra_axes !< dimensions other than space-time type(vardesc) :: vd character(len=32), dimension(:), allocatable :: dim_names From 337182e7ac17f5305c3b33fbadde05558d8baf30 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 23 Aug 2023 16:25:07 -0400 Subject: [PATCH 4/9] fix read restart for 4d array * passes global_file optional arg to allow reading in parallel restart files --- config_src/infra/FMS1/MOM_io_infra.F90 | 57 +++++++++++++++----------- src/framework/MOM_restart.F90 | 2 +- 2 files changed, 35 insertions(+), 24 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index e37e5db3cb..e62216d0b0 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -739,7 +739,7 @@ end subroutine read_field_3d_region !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale, global_file) + timelevel, position, scale, global_file, file_may_be_4d) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The variable name of the data in the file real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data @@ -750,11 +750,14 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. logical, optional, intent(in) :: global_file !< If true, read from a single global file + logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, + !! in which case a more elaborate set of calls + !! is needed to read it due to FMS limitations. ! Local variables character(len=80) :: varname ! The name of a variable in the file type(fieldtype), allocatable :: fields(:) ! An array of types describing all the variables in the file - logical :: file_is_global + logical :: use_fms_read_data, file_is_global integer :: n, unit, ndim, nvar, natt, ntime ! This single call does not work for a 4-d array due to FMS limitations, so multiple calls are @@ -762,32 +765,40 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & ! call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & ! timelevel=timelevel, position=position) + use_fms_read_data = .true. ; if (present(file_may_be_4d)) use_fms_read_data = .not.file_may_be_4d file_is_global = .true. ; if (present(global_file)) file_is_global = global_file - if (file_is_global) then - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) + use_fms_read_data = .false. ! 4d arrays not working with FMS1 + + if (use_fms_read_data) then + call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=position) else - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) - endif - call mpp_get_info(unit, ndim, nvar, natt, ntime) - allocate(fields(nvar)) - call mpp_get_fields(unit, fields(1:nvar)) - do n=1, nvar - call mpp_get_atts(fields(n), name=varname) - if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then - call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) - ! Maybe something should be done depending on the value of ntime. - call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) - exit + if (file_is_global) then + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) + else + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) endif - enddo - if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & - "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) + call mpp_get_info(unit, ndim, nvar, natt, ntime) + allocate(fields(nvar)) + call mpp_get_fields(unit, fields(1:nvar)) + do n=1, nvar + call mpp_get_atts(fields(n), name=varname) + if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then + call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) + ! Maybe something should be done depending on the value of ntime. + call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) + exit + endif + enddo + if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & + "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) - deallocate(fields) - call mpp_close(unit) + deallocate(fields) + call mpp_close(unit) + endif if (present(scale)) then ; if (scale /= 1.0) then call rescale_comp_data(MOM_Domain, data, scale) diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 24e6da0efd..fa1481570a 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -1700,7 +1700,7 @@ subroutine restore_state(filename, directory, day, G, CS) elseif (associated(CS%var_ptr4d(m)%p)) then ! Read a 4d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & - G%Domain, timelevel=1, position=pos, scale=scale) + G%Domain, timelevel=1, position=pos, scale=scale, global_file=unit_is_global(n)) else ! This array is not domain-decomposed. This variant may be under-tested. call MOM_error(FATAL, & "MOM_restart does not support 4-d arrays without domain decomposition.") From 3489833fd312bb0a80a059dca7a4ca8639a8e276 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 23 Aug 2023 16:29:15 -0400 Subject: [PATCH 5/9] write one restart variable per vertical mode vertical modes >= 6 do not propagate hence there is no need for dynamic array size on this dimension --- .../lateral/MOM_internal_tides.F90 | 201 ++++++++++++++---- 1 file changed, 161 insertions(+), 40 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index acb58fda5b..867e854c97 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -140,8 +140,17 @@ module MOM_internal_tides real, allocatable :: En(:,:,:,:,:) !< The internal wave energy density as a function of (i,j,angle,frequency,mode) !! integrated within an angular and frequency band [R Z3 T-2 ~> J m-2] - real, allocatable :: En_restart(:,:,:,:) - !< The internal wave energy density as a function of (i,j,angle); temporary for restart + real, allocatable :: En_restart_mode1(:,:,:,:) + !< The internal wave energy density as a function of (i,j,angle,freq) for mode 1 + real, allocatable :: En_restart_mode2(:,:,:,:) + !< The internal wave energy density as a function of (i,j,angle,freq) for mode 2 + real, allocatable :: En_restart_mode3(:,:,:,:) + !< The internal wave energy density as a function of (i,j,angle,freq) for mode 3 + real, allocatable :: En_restart_mode4(:,:,:,:) + !< The internal wave energy density as a function of (i,j,angle,freq) for mode 4 + real, allocatable :: En_restart_mode5(:,:,:,:) + !< The internal wave energy density as a function of (i,j,angle,freq) for mode 5 + real, allocatable, dimension(:) :: frequency !< The frequency of each band [T-1 ~> s-1]. type(wave_speed_CS) :: wave_speed !< Wave speed control structure @@ -269,11 +278,35 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, cn(:,:,:) = 0. - ! temporary solution for restart, only works with one freq/1mode + ! Rebuild energy density array from multiple restarts do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied - CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) + CS%En(i,j,a,fr,1) = CS%En_restart_mode1(i,j,a,fr) enddo ; enddo ; enddo ; enddo + if (CS%nMode >= 2) then + do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,2) = CS%En_restart_mode2(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 3) then + do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,3) = CS%En_restart_mode3(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 4) then + do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,4) = CS%En_restart_mode4(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 5) then + do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,5) = CS%En_restart_mode5(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + endif + ! Set properties related to the internal tides, such as the wave speeds, storing some ! of them in the control structure for this module. if (CS%uniform_test_cg > 0.0) then @@ -331,7 +364,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -352,7 +385,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, call correct_halo_rotation(CS%En, test, G, CS%nAngle) ! Propagate the waves. - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq CS%TKE_residual_loss(:,:,:,fr,m) = 0. @@ -361,7 +394,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset @@ -377,13 +410,13 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo ; enddo ! Apply the other half of the refraction. - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq call refract(CS%En(:,:,:,fr,m), cn(:,:,m), CS%frequency(fr), 0.5*dt, & G, US, CS%NAngle, CS%use_PPMang) enddo ; enddo ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -402,7 +435,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, .or. (CS%id_tot_En > 0)) then tot_En(:,:) = 0.0 tot_En_mode(:,:,:,:) = 0.0 - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq do j=jsd,jed ; do i=isd,ied ; do a=1,CS%nAngle tot_En(i,j) = tot_En(i,j) + CS%En(i,j,a,fr,m) tot_En_mode(i,j,fr,m) = tot_En_mode(i,j,fr,m) + CS%En(i,j,a,fr,m) @@ -420,7 +453,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo ; enddo ; enddo ; enddo endif ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -461,7 +494,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo ; enddo ; enddo ; enddo endif ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -479,7 +512,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! still need to allow a portion of the extracted energy to go to higher modes. ! First, find velocity profiles if (CS%apply_wave_drag .or. CS%apply_Froude_drag) then - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq ! compute near-bottom and max horizontal baroclinic velocity values at each point do j=jsd,jed ; do i=isd,ied @@ -527,7 +560,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, CS%TKE_itidal_loss, dt, full_halos=.false.) endif ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -543,7 +576,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Extract the energy for mixing due to wave breaking----------------------------- if (CS%apply_Froude_drag) then ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq freq2 = CS%frequency(fr)**2 do j=js,je ; do i=is,ie id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging @@ -594,7 +627,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, enddo ; enddo endif ! Check for En<0 - for debugging, delete later - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle do j=js,je ; do i=is,ie if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset @@ -626,7 +659,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Check for energy conservation on computational domain.************************* - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide') enddo ; enddo @@ -646,7 +679,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, TKE_itidal_input, CS%diag) ! Output 2-D energy density (summed over angles) for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_En_mode(fr,m) > 0) then + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_En_mode(fr,m) > 0) then tot_En(:,:) = 0.0 do a=1,CS%nAngle ; do j=js,je ; do i=is,ie tot_En(i,j) = tot_En(i,j) + CS%En(i,j,a,fr,m) @@ -654,13 +687,37 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, call post_data(CS%id_En_mode(fr,m), tot_En, CS%diag) endif ; enddo ; enddo - ! save to temporary restart + ! split energy array into multiple restarts do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie - CS%En_restart(i,j,a,fr) = CS%En(i,j,a,fr,1) + CS%En_restart_mode1(i,j,a,fr) = CS%En(i,j,a,fr,1) enddo ; enddo ; enddo ; enddo + if (CS%nMode >= 2) then + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart_mode2(i,j,a,fr) = CS%En(i,j,a,fr,2) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 3) then + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart_mode3(i,j,a,fr) = CS%En(i,j,a,fr,3) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 4) then + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart_mode4(i,j,a,fr) = CS%En(i,j,a,fr,4) + enddo ; enddo ; enddo ; enddo + endif + + if (CS%nMode >= 5) then + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + CS%En_restart_mode5(i,j,a,fr) = CS%En(i,j,a,fr,5) + enddo ; enddo ; enddo ; enddo + endif + ! Output 3-D (i,j,a) energy density for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_En_ang_mode(fr,m) > 0) then + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_En_ang_mode(fr,m) > 0) then call post_data(CS%id_En_ang_mode(fr,m), CS%En(:,:,:,fr,m) , CS%diag) endif ; enddo ; enddo @@ -671,7 +728,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, tot_Froude_loss(:,:) = 0.0 tot_residual_loss(:,:) = 0.0 tot_allprocesses_loss(:,:) = 0.0 - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie tot_leak_loss(i,j) = tot_leak_loss(i,j) + CS%TKE_leak_loss(i,j,a,fr,m) tot_quad_loss(i,j) = tot_quad_loss(i,j) + CS%TKE_quad_loss(i,j,a,fr,m) tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + CS%TKE_itidal_loss(i,j,a,fr,m) @@ -709,7 +766,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, endif ! Output 2-D energy loss (summed over angles) for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq + do m=1,CS%nMode ; do fr=1,CS%Nfreq if (CS%id_itidal_loss_mode(fr,m) > 0 .or. CS%id_allprocesses_loss_mode(fr,m) > 0) then itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together @@ -725,37 +782,37 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, endif ; enddo ; enddo ! Output 3-D (i,j,a) energy loss for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_itidal_loss_ang_mode(fr,m) > 0) then + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_itidal_loss_ang_mode(fr,m) > 0) then call post_data(CS%id_itidal_loss_ang_mode(fr,m), CS%TKE_itidal_loss(:,:,:,fr,m) , CS%diag) endif ; enddo ; enddo ! Output 2-D period-averaged horizontal near-bottom mode velocity for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_Ub_mode(fr,m) > 0) then + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_Ub_mode(fr,m) > 0) then call post_data(CS%id_Ub_mode(fr,m), Ub(:,:,fr,m), CS%diag) endif ; enddo ; enddo - do m=1,CS%NMode ; if (CS%id_Ustruct_mode(m) > 0) then + do m=1,CS%nMode ; if (CS%id_Ustruct_mode(m) > 0) then call post_data(CS%id_Ustruct_mode(m), CS%u_struct(:,:,:,m), CS%diag) endif ; enddo - do m=1,CS%NMode ; if (CS%id_Wstruct_mode(m) > 0) then + do m=1,CS%nMode ; if (CS%id_Wstruct_mode(m) > 0) then call post_data(CS%id_Wstruct_mode(m), CS%w_struct(:,:,:,m), CS%diag) endif ; enddo - do m=1,CS%NMode ; if (CS%id_int_w2_mode(m) > 0) then + do m=1,CS%nMode ; if (CS%id_int_w2_mode(m) > 0) then call post_data(CS%id_int_w2_mode(m), CS%int_w2(:,:,m), CS%diag) endif ; enddo - do m=1,CS%NMode ; if (CS%id_int_U2_mode(m) > 0) then + do m=1,CS%nMode ; if (CS%id_int_U2_mode(m) > 0) then call post_data(CS%id_int_U2_mode(m), CS%int_U2(:,:,m), CS%diag) endif ; enddo - do m=1,CS%NMode ; if (CS%id_int_N2w2_mode(m) > 0) then + do m=1,CS%nMode ; if (CS%id_int_N2w2_mode(m) > 0) then call post_data(CS%id_int_N2w2_mode(m), CS%int_N2w2(:,:,m), CS%diag) endif ; enddo ! Output 2-D horizontal phase velocity for each frequency and mode - do m=1,CS%NMode ; do fr=1,CS%Nfreq ; if (CS%id_cp_mode(fr,m) > 0) then + do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_cp_mode(fr,m) > 0) then call post_data(CS%id_cp_mode(fr,m), CS%cp(:,:,fr,m), CS%diag) endif ; enddo ; enddo @@ -2329,16 +2386,80 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) call set_axis_info(axes_inttides(1), "angle", "", "angle direction", num_angle, angles, "N", 1) call set_axis_info(axes_inttides(2), "freq", "", "wave frequency", num_freq, freqs, "N", 1) - ! repeat for vertical modes - allocate(CS%En_restart(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + ! full energy array allocate(CS%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0) - ! temporary 4d restart - call register_restart_field(CS%En_restart(:,:,:,:), "IW_energy", .false., restart_CS, & - longname="The internal wave energy density as a function of (i,j,angle,freq)", & + ! restart strategy: support for 5d restart is not yet available so we split into + ! 4d restarts. Vertical modes >= 6 are dissipated locally and do not propagate + ! so we only allow for 5 vertical modes and each has its own variable + + ! allocate restart arrays + allocate(CS%En_restart_mode1(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + if (num_mode >= 2) allocate(CS%En_restart_mode2(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + if (num_mode >= 3) allocate(CS%En_restart_mode3(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + if (num_mode >= 4) allocate(CS%En_restart_mode4(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + if (num_mode >= 5) allocate(CS%En_restart_mode5(isd:ied, jsd:jed, num_angle, num_freq), source=0.0) + + ! register all 4d restarts and copy into full Energy array when restarting from previous state + call register_restart_field(CS%En_restart_mode1(:,:,:,:), "IW_energy_mode1", .false., restart_CS, & + longname="The internal wave energy density f(i,j,angle,freq) for mode 1", & units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & extra_axes=axes_inttides) + + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,1) = CS%En_restart_mode1(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + + if (num_mode >= 2) then + call register_restart_field(CS%En_restart_mode2(:,:,:,:), "IW_energy_mode2", .false., restart_CS, & + longname="The internal wave energy density f(i,j,angle,freq) for mode 2", & + units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + extra_axes=axes_inttides) + + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,2) = CS%En_restart_mode2(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + + endif + + if (num_mode >= 3) then + call register_restart_field(CS%En_restart_mode3(:,:,:,:), "IW_energy_mode3", .false., restart_CS, & + longname="The internal wave energy density f(i,j,angle,freq) for mode 3", & + units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + extra_axes=axes_inttides) + + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,3) = CS%En_restart_mode3(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + + endif + + if (num_mode >= 4) then + call register_restart_field(CS%En_restart_mode4(:,:,:,:), "IW_energy_mode4", .false., restart_CS, & + longname="The internal wave energy density f(i,j,angle,freq) for mode 4", & + units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + extra_axes=axes_inttides) + + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,4) = CS%En_restart_mode4(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + + endif + + if (num_mode >= 5) then + call register_restart_field(CS%En_restart_mode5(:,:,:,:), "IW_energy_mode5", .false., restart_CS, & + longname="The internal wave energy density f(i,j,angle,freq) for mode 5", & + units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & + extra_axes=axes_inttides) + + do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + CS%En(i,j,a,fr,5) = CS%En_restart_mode5(i,j,a,fr) + enddo ; enddo ; enddo ; enddo + + endif + + ! full 5d restart !call register_restart_field(CS%En(:,:,:,:,:), "IW_energy", .false., restart_CS, & ! longname="The internal wave energy density as a function of (i,j,angle,freq,mode)", & @@ -2347,9 +2468,9 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) !call restart_registry_lock(restart_CS, unlocked=.false.) ! read the restart data back in at init - do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied - CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) - enddo ; enddo ; enddo ; enddo + !do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied + ! CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) + !enddo ; enddo ; enddo ; enddo end subroutine register_int_tide_restarts @@ -2468,7 +2589,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "Inconsistent number of frequencies.") if (CS%NAngle /= num_angle) call MOM_error(FATAL, "Internal_tides_init: "//& "Inconsistent number of angles.") - if (CS%NMode /= num_mode) call MOM_error(FATAL, "Internal_tides_init: "//& + if (CS%nMode /= num_mode) call MOM_error(FATAL, "Internal_tides_init: "//& "Inconsistent number of modes.") if (4*(num_angle/4) /= num_angle) call MOM_error(FATAL, & "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.") From bb281445c1d36e032165251c90d7f12e58b0b087 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 6 Sep 2023 11:35:49 -0400 Subject: [PATCH 6/9] fix minor issues from review --- src/core/MOM.F90 | 2 +- src/framework/MOM_restart.F90 | 14 +++++++--- .../lateral/MOM_internal_tides.F90 | 27 +++++-------------- .../vertical/MOM_diabatic_driver.F90 | 2 +- 4 files changed, 20 insertions(+), 25 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 63a9fc1739..301f563e83 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2763,7 +2763,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call stoch_EOS_register_restarts(HI, param_file, CS%stoch_eos_CS, restart_CSp) endif - if ( .not. CS%adiabatic) then + if (.not. CS%adiabatic) then call register_diabatic_restarts(G, US, param_file, CS%int_tide_CSp, restart_CSp) endif diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index fa1481570a..252f14bfac 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -469,9 +469,13 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units character(len=32), dimension(:), allocatable :: dim_names integer :: n, n_extradims + ! first 2 dimensions in dim_names are reserved for i,j + ! so extra_dimensions are shifted to index 3. + ! this is designed not to break the behavior in SIS2 + ! (see register_restart_field_4d in SIS_restart.F90) if (present(extra_axes)) then n_extradims = size(extra_axes) - allocate( dim_names(n_extradims+2) ) + allocate(dim_names(n_extradims+2)) dim_names(1) = "" dim_names(2) = "" do n=3,n_extradims+2 @@ -521,9 +525,13 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units character(len=32), dimension(:), allocatable :: dim_names integer :: n, n_extradims + ! first 2 dimensions in dim_names are reserved for i,j + ! so extra_dimensions are shifted to index 3. + ! this is designed not to break the behavior in SIS2 + ! (see register_restart_field_4d in SIS_restart.F90) if (present(extra_axes)) then n_extradims = size(extra_axes) - allocate( dim_names(n_extradims+2) ) + allocate(dim_names(n_extradims+2)) dim_names(1) = "" dim_names(2) = "" do n=3,n_extradims+2 @@ -1365,7 +1373,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ turns = CS%turns - allocate ( extra_axes(nmax_extradims) ) + allocate (extra_axes(nmax_extradims)) if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "save_restart: Module must be initialized before it is used.") diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 867e854c97..83910e6690 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -688,30 +688,30 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, endif ; enddo ; enddo ! split energy array into multiple restarts - do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied CS%En_restart_mode1(i,j,a,fr) = CS%En(i,j,a,fr,1) enddo ; enddo ; enddo ; enddo if (CS%nMode >= 2) then - do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied CS%En_restart_mode2(i,j,a,fr) = CS%En(i,j,a,fr,2) enddo ; enddo ; enddo ; enddo endif if (CS%nMode >= 3) then - do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied CS%En_restart_mode3(i,j,a,fr) = CS%En(i,j,a,fr,3) enddo ; enddo ; enddo ; enddo endif if (CS%nMode >= 4) then - do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied CS%En_restart_mode4(i,j,a,fr) = CS%En(i,j,a,fr,4) enddo ; enddo ; enddo ; enddo endif if (CS%nMode >= 5) then - do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie + do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied CS%En_restart_mode5(i,j,a,fr) = CS%En(i,j,a,fr,5) enddo ; enddo ; enddo ; enddo endif @@ -2378,8 +2378,8 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) call get_param(param_file, "MOM", "INTERNAL_TIDE_FREQS", num_freq, default=1) call get_param(param_file, "MOM", "INTERNAL_TIDE_MODES", num_mode, default=1) - allocate ( angles(num_angle) ) - allocate ( freqs(num_freq) ) + allocate (angles(num_angle)) + allocate (freqs(num_freq)) do a=1,num_angle ; angles(a)= a ; enddo do fr=1,num_freq ; freqs(fr)= fr ; enddo @@ -2459,19 +2459,6 @@ subroutine register_int_tide_restarts(G, US, param_file, CS, restart_CS) endif - - ! full 5d restart - !call register_restart_field(CS%En(:,:,:,:,:), "IW_energy", .false., restart_CS, & - ! longname="The internal wave energy density as a function of (i,j,angle,freq,mode)", & - ! units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1', t_grid="s", & - ! extra_axes=axes_inttides) - !call restart_registry_lock(restart_CS, unlocked=.false.) - - ! read the restart data back in at init - !do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied - ! CS%En(i,j,a,fr,1) = CS%En_restart(i,j,a,fr) - !enddo ; enddo ; enddo ; enddo - end subroutine register_int_tide_restarts !> This subroutine initializes the internal tides module. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 4f0baee22f..8053c16265 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3573,7 +3573,7 @@ subroutine register_diabatic_restarts(G, US, param_file, int_tide_CSp, restart_C type(int_tide_CS), pointer :: int_tide_CSp !< Internal tide control structure type(MOM_restart_CS), pointer :: restart_CSp !< MOM restart control structure - logical :: use_int_tides = .false. + logical :: use_int_tides call read_param(param_file, "INTERNAL_TIDES", use_int_tides) From 95cbc78be028b1c4edb4857af25d2830feecd718 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 6 Sep 2023 11:43:28 -0400 Subject: [PATCH 7/9] revert changes to FMS1 infra --- config_src/infra/FMS1/MOM_io_infra.F90 | 98 +++++++------------------- 1 file changed, 24 insertions(+), 74 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index e62216d0b0..c0ccfcbcc8 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -57,7 +57,7 @@ module MOM_io_infra !> Read a data field from a file interface read_field module procedure read_field_4d - module procedure read_field_3d, read_field_3d_region + module procedure read_field_3d module procedure read_field_2d, read_field_2d_region module procedure read_field_1d, read_field_1d_int module procedure read_field_0d, read_field_0d_int @@ -696,50 +696,11 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, & endif ; endif end subroutine read_field_3d -!> This routine uses the fms_io subroutine read_data to read a region from a distributed or -!! global 3-D data field named "fieldname" from file "filename". -subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, & - no_domain, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data - !! should be read - integer, dimension(:), intent(in) :: start !< The starting index to read in each of 4 - !! dimensions. For this 3-d read, the - !! 4th values are always 1. - integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4 - !! dimensions. For this 3-d read, the - !! 4th values are always 1. - type(MOM_domain_type), & - optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - logical, optional, intent(in) :: no_domain !< If present and true, this variable does not - !! use domain decomposion. - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before it is returned. - - if (present(MOM_Domain)) then - call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, & - no_domain=no_domain) - else - call read_data(filename, fieldname, data, start, nread, no_domain=no_domain) - endif - - if (present(scale)) then ; if (scale /= 1.0) then - if (present(MOM_Domain)) then - call rescale_comp_data(MOM_Domain, data, scale) - else - ! Dangerously rescale the whole array - data(:,:,:) = scale*data(:,:,:) - endif - endif ; endif -end subroutine read_field_3d_region - - !> This routine uses the fms_io subroutine read_data to read a distributed !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale, global_file, file_may_be_4d) + timelevel, position, scale, global_file) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The variable name of the data in the file real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data @@ -750,14 +711,11 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. logical, optional, intent(in) :: global_file !< If true, read from a single global file - logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, - !! in which case a more elaborate set of calls - !! is needed to read it due to FMS limitations. ! Local variables character(len=80) :: varname ! The name of a variable in the file type(fieldtype), allocatable :: fields(:) ! An array of types describing all the variables in the file - logical :: use_fms_read_data, file_is_global + logical :: file_is_global integer :: n, unit, ndim, nvar, natt, ntime ! This single call does not work for a 4-d array due to FMS limitations, so multiple calls are @@ -765,40 +723,32 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & ! call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & ! timelevel=timelevel, position=position) - use_fms_read_data = .true. ; if (present(file_may_be_4d)) use_fms_read_data = .not.file_may_be_4d file_is_global = .true. ; if (present(global_file)) file_is_global = global_file - use_fms_read_data = .false. ! 4d arrays not working with FMS1 - - if (use_fms_read_data) then - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=position) + if (file_is_global) then + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) else - if (file_is_global) then - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) - else - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) + endif + call mpp_get_info(unit, ndim, nvar, natt, ntime) + allocate(fields(nvar)) + call mpp_get_fields(unit, fields(1:nvar)) + do n=1, nvar + call mpp_get_atts(fields(n), name=varname) + if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then + call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) + ! Maybe something should be done depending on the value of ntime. + call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) + exit endif - call mpp_get_info(unit, ndim, nvar, natt, ntime) - allocate(fields(nvar)) - call mpp_get_fields(unit, fields(1:nvar)) - do n=1, nvar - call mpp_get_atts(fields(n), name=varname) - if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then - call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) - ! Maybe something should be done depending on the value of ntime. - call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) - exit - endif - enddo - if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & - "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) + enddo + if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & + "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) - deallocate(fields) - call mpp_close(unit) - endif + deallocate(fields) + call mpp_close(unit) if (present(scale)) then ; if (scale /= 1.0) then call rescale_comp_data(MOM_Domain, data, scale) From 0e38327715b3f5db811d66d4ebfd4a376fde3208 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 6 Sep 2023 11:46:25 -0400 Subject: [PATCH 8/9] Revert "revert changes to FMS1 infra" This reverts commit 95cbc78be028b1c4edb4857af25d2830feecd718. --- config_src/infra/FMS1/MOM_io_infra.F90 | 98 +++++++++++++++++++------- 1 file changed, 74 insertions(+), 24 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index c0ccfcbcc8..e62216d0b0 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -57,7 +57,7 @@ module MOM_io_infra !> Read a data field from a file interface read_field module procedure read_field_4d - module procedure read_field_3d + module procedure read_field_3d, read_field_3d_region module procedure read_field_2d, read_field_2d_region module procedure read_field_1d, read_field_1d_int module procedure read_field_0d, read_field_0d_int @@ -696,11 +696,50 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, & endif ; endif end subroutine read_field_3d +!> This routine uses the fms_io subroutine read_data to read a region from a distributed or +!! global 3-D data field named "fieldname" from file "filename". +subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, & + no_domain, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data + !! should be read + integer, dimension(:), intent(in) :: start !< The starting index to read in each of 4 + !! dimensions. For this 3-d read, the + !! 4th values are always 1. + integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4 + !! dimensions. For this 3-d read, the + !! 4th values are always 1. + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + logical, optional, intent(in) :: no_domain !< If present and true, this variable does not + !! use domain decomposion. + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + if (present(MOM_Domain)) then + call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, & + no_domain=no_domain) + else + call read_data(filename, fieldname, data, start, nread, no_domain=no_domain) + endif + + if (present(scale)) then ; if (scale /= 1.0) then + if (present(MOM_Domain)) then + call rescale_comp_data(MOM_Domain, data, scale) + else + ! Dangerously rescale the whole array + data(:,:,:) = scale*data(:,:,:) + endif + endif ; endif +end subroutine read_field_3d_region + + !> This routine uses the fms_io subroutine read_data to read a distributed !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale, global_file) + timelevel, position, scale, global_file, file_may_be_4d) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The variable name of the data in the file real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data @@ -711,11 +750,14 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. logical, optional, intent(in) :: global_file !< If true, read from a single global file + logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, + !! in which case a more elaborate set of calls + !! is needed to read it due to FMS limitations. ! Local variables character(len=80) :: varname ! The name of a variable in the file type(fieldtype), allocatable :: fields(:) ! An array of types describing all the variables in the file - logical :: file_is_global + logical :: use_fms_read_data, file_is_global integer :: n, unit, ndim, nvar, natt, ntime ! This single call does not work for a 4-d array due to FMS limitations, so multiple calls are @@ -723,32 +765,40 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & ! call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & ! timelevel=timelevel, position=position) + use_fms_read_data = .true. ; if (present(file_may_be_4d)) use_fms_read_data = .not.file_may_be_4d file_is_global = .true. ; if (present(global_file)) file_is_global = global_file - if (file_is_global) then - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) + use_fms_read_data = .false. ! 4d arrays not working with FMS1 + + if (use_fms_read_data) then + call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=position) else - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) - endif - call mpp_get_info(unit, ndim, nvar, natt, ntime) - allocate(fields(nvar)) - call mpp_get_fields(unit, fields(1:nvar)) - do n=1, nvar - call mpp_get_atts(fields(n), name=varname) - if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then - call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) - ! Maybe something should be done depending on the value of ntime. - call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) - exit + if (file_is_global) then + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) + else + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) endif - enddo - if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & - "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) + call mpp_get_info(unit, ndim, nvar, natt, ntime) + allocate(fields(nvar)) + call mpp_get_fields(unit, fields(1:nvar)) + do n=1, nvar + call mpp_get_atts(fields(n), name=varname) + if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then + call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) + ! Maybe something should be done depending on the value of ntime. + call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) + exit + endif + enddo + if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & + "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) - deallocate(fields) - call mpp_close(unit) + deallocate(fields) + call mpp_close(unit) + endif if (present(scale)) then ; if (scale /= 1.0) then call rescale_comp_data(MOM_Domain, data, scale) From 5405c3b05ae8aca7b603f2dfba3d51f26bbc0f04 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 6 Sep 2023 11:47:58 -0400 Subject: [PATCH 9/9] remove changes to infra FMS1 --- config_src/infra/FMS1/MOM_io_infra.F90 | 57 ++++++++----------- .../vertical/MOM_diabatic_driver.F90 | 2 + 2 files changed, 25 insertions(+), 34 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index e62216d0b0..e37e5db3cb 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/config_src/infra/FMS1/MOM_io_infra.F90 @@ -739,7 +739,7 @@ end subroutine read_field_3d_region !! 4-D data field named "fieldname" from file "filename". Valid values for !! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale, global_file, file_may_be_4d) + timelevel, position, scale, global_file) character(len=*), intent(in) :: filename !< The name of the file to read character(len=*), intent(in) :: fieldname !< The variable name of the data in the file real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data @@ -750,14 +750,11 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied !! by before it is returned. logical, optional, intent(in) :: global_file !< If true, read from a single global file - logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, - !! in which case a more elaborate set of calls - !! is needed to read it due to FMS limitations. ! Local variables character(len=80) :: varname ! The name of a variable in the file type(fieldtype), allocatable :: fields(:) ! An array of types describing all the variables in the file - logical :: use_fms_read_data, file_is_global + logical :: file_is_global integer :: n, unit, ndim, nvar, natt, ntime ! This single call does not work for a 4-d array due to FMS limitations, so multiple calls are @@ -765,40 +762,32 @@ subroutine read_field_4d(filename, fieldname, data, MOM_Domain, & ! call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & ! timelevel=timelevel, position=position) - use_fms_read_data = .true. ; if (present(file_may_be_4d)) use_fms_read_data = .not.file_may_be_4d file_is_global = .true. ; if (present(global_file)) file_is_global = global_file - use_fms_read_data = .false. ! 4d arrays not working with FMS1 - - if (use_fms_read_data) then - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=position) + if (file_is_global) then + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) else - if (file_is_global) then - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=SINGLE_FILE) !, domain=MOM_Domain%mpp_domain ) - else - call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & - threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) + call mpp_open(unit, trim(filename), form=NETCDF_FILE, action=READONLY_FILE, & + threading=MULTIPLE, fileset=MULTIPLE, domain=MOM_Domain%mpp_domain ) + endif + call mpp_get_info(unit, ndim, nvar, natt, ntime) + allocate(fields(nvar)) + call mpp_get_fields(unit, fields(1:nvar)) + do n=1, nvar + call mpp_get_atts(fields(n), name=varname) + if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then + call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) + ! Maybe something should be done depending on the value of ntime. + call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) + exit endif - call mpp_get_info(unit, ndim, nvar, natt, ntime) - allocate(fields(nvar)) - call mpp_get_fields(unit, fields(1:nvar)) - do n=1, nvar - call mpp_get_atts(fields(n), name=varname) - if (lowercase(trim(varname)) == lowercase(trim(fieldname))) then - call MOM_error(NOTE, "Reading 4-d variable "//trim(fieldname)//" from file "//trim(filename)) - ! Maybe something should be done depending on the value of ntime. - call mpp_read(unit, fields(n), MOM_Domain%mpp_domain, data, timelevel) - exit - endif - enddo - if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & - "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) + enddo + if ((n == nvar+1) .or. (nvar < 1)) call MOM_error(WARNING, & + "read_field apparently did not find 4-d variable "//trim(fieldname)//" in file "//trim(filename)) - deallocate(fields) - call mpp_close(unit) - endif + deallocate(fields) + call mpp_close(unit) if (present(scale)) then ; if (scale /= 1.0) then call rescale_comp_data(MOM_Domain, data, scale) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 8053c16265..1ccd6a7fb2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3575,6 +3575,8 @@ subroutine register_diabatic_restarts(G, US, param_file, int_tide_CSp, restart_C logical :: use_int_tides + use_int_tides=.false. + call read_param(param_file, "INTERNAL_TIDES", use_int_tides) if (use_int_tides) then