diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ce001483ff..d112de07b7 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, thickness_to_dz 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 @@ -1964,6 +1968,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & 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 @@ -2774,6 +2779,10 @@ 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 + 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) @@ -3139,7 +3148,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & 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 220a7d6bcf..27d244b226 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -148,6 +148,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 @@ -165,22 +179,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) !< dimensions other than space-time 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 @@ -271,7 +272,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 @@ -1751,7 +1753,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 @@ -1772,6 +1775,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 !< 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 @@ -1795,7 +1800,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 @@ -1803,7 +1809,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 @@ -1825,6 +1832,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 !< dimensions other than space-time character(len=120) :: cllr integer :: n @@ -1877,6 +1886,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) @@ -2020,7 +2035,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 @@ -2035,6 +2050,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 !< 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 @@ -2043,7 +2060,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) @@ -2076,6 +2094,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..252f14bfac 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,26 @@ 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 !< dimensions other than space-time type(vardesc) :: vd + 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)) + 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 +489,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 +503,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 +518,26 @@ 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 !< dimensions other than space-time type(vardesc) :: vd + 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)) + 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 +545,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 +1357,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 +1368,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 +1413,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 +1503,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 @@ -1650,7 +1708,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.") diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 788e922ff2..83910e6690 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 @@ -137,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 @@ -266,6 +278,35 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, cn(:,:,:) = 0. + ! 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_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 @@ -323,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 @@ -344,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. @@ -353,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 @@ -369,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 @@ -394,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) @@ -412,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 @@ -453,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 @@ -471,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 @@ -519,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 @@ -535,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 @@ -586,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 @@ -618,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 @@ -638,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) @@ -646,8 +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 + ! split energy array into multiple restarts + 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=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=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=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=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 + ! 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 @@ -658,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) @@ -696,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 @@ -712,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 @@ -2276,43 +2346,120 @@ 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 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 + + type(axis_info) :: axes_inttides(2) + real, dimension(:), allocatable :: angles, freqs + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + if (associated(CS)) then + call MOM_error(WARNING, "register_int_tide_restarts called "//& + "with an associated control structure.") + return + endif + + allocate(CS) + + ! 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) + + 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", num_angle, angles, "N", 1) + call set_axis_info(axes_inttides(2), "freq", "", "wave frequency", num_freq, freqs, "N", 1) + + ! full energy array + allocate(CS%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0) + + ! 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) -! ! This subroutine is not currently in use!! + 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 -! ! 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 + 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) -! if (associated(CS)) then -! call MOM_error(WARNING, "register_int_tide_restarts called "//& -! "with an associated control structure.") -! return -! endif + 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 -! use_int_tides = .false. -! call read_param(param_file, "INTERNAL_TIDES", use_int_tides) -! if (.not.use_int_tides) return + endif -! allocate(CS) + 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) -! 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) + 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 -! 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') + endif -! end subroutine register_int_tide_restarts +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) @@ -2324,7 +2471,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] @@ -2358,6 +2505,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 @@ -2376,9 +2525,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) @@ -2430,7 +2576,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.") diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b2b8527819..1ccd6a7fb2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -45,7 +45,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, thickness_to_dz -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 @@ -56,6 +56,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 @@ -81,6 +82,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 @@ -226,12 +228,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 @@ -386,7 +389,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, CS%int_tide_input%Rho_bot, dt, G, GV, US, CS%int_tide) + CS%int_tide_input%Nb, CS%int_tide_input%Rho_bot, 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 @@ -2980,7 +2983,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 @@ -2998,6 +3001,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] @@ -3032,6 +3036,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) @@ -3497,12 +3502,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) @@ -3558,6 +3565,25 @@ 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 + + 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)