implement restart for internal tides

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,44 +750,55 @@ 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
   ! needed.
   ! 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)

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
@@ -1943,6 +1947,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
@@ -2758,6 +2763,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)
 
@@ -3123,7 +3132,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)) &

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,15 +1800,17 @@ 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
 
 
 !> 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