+Refactor writing the vertical_coordinate file

  Rearrange how the vertical_coordinate file is written, to avoid an unnecessary
repetition of the writing of this file, and to prepare for the Hybgen grid
generator code to write its own version of this file.  The specific changes
include:

 - Added the new subroutine write_regrid_file() to write the vertical_coordinate
   file when in ALE mode.

 - Relocated the code to write the vertical_coordinate file in ALE mode
   from ALE_writeCoordinateFile to write_regrid_file()

 - Moved the call to write_vertgrid_file() into initialize_MOM(), and do not
   call it when USE_REGRIDDING=True, since in this case this file will be
   rewritten again when write_regrid_file() is called.

 - Eliminated the write_geom and output_dir arguments to MOM_initialize_coord()

All answers and output are bitwise identical.

src/ALE/MOM_ALE.F90

@@ -21,8 +21,6 @@ module MOM_ALE
 use MOM_error_handler,    only : callTree_showQuery
 use MOM_error_handler,    only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser,      only : get_param, param_file_type, log_param
-use MOM_io,               only : vardesc, var_desc, fieldtype, SINGLE_FILE
-use MOM_io,               only : create_file, write_field, close_file, file_type
 use MOM_interface_heights,only : find_eta
 use MOM_open_boundary,    only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
 use MOM_open_boundary,    only : OBC_DIRECTION_N, OBC_DIRECTION_S
@@ -34,8 +32,8 @@ module MOM_ALE
 use MOM_regridding,       only : regriddingInterpSchemeDoc, regriddingDefaultInterpScheme
 use MOM_regridding,       only : regriddingDefaultBoundaryExtrapolation
 use MOM_regridding,       only : regriddingDefaultMinThickness
-use MOM_regridding,       only : regridding_CS, set_regrid_params
-use MOM_regridding,       only : getCoordinateInterfaces, getCoordinateResolution
+use MOM_regridding,       only : regridding_CS, set_regrid_params, write_regrid_file
+use MOM_regridding,       only : getCoordinateInterfaces
 use MOM_regridding,       only : getCoordinateUnits, getCoordinateShortName
 use MOM_regridding,       only : getStaticThickness
 use MOM_remapping,        only : initialize_remapping, end_remapping
@@ -333,7 +331,7 @@ end subroutine ALE_register_diags
 !! We read the MOM_input file to register the values of different
 !! regridding/remapping parameters.
 subroutine adjustGridForIntegrity( CS, G, GV, h )
-  type(ALE_CS),                              pointer       :: CS  !< Regridding parameters and options
+  type(ALE_CS),                              intent(in)    :: CS  !< Regridding parameters and options
   type(ocean_grid_type),                     intent(in)    :: G   !< Ocean grid informations
   type(verticalGrid_type),                   intent(in)    :: GV  !< Ocean vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h   !< Current 3D grid thickness that
@@ -1413,26 +1411,10 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory )
   character(len=*),        intent(in)  :: directory  !< directory for writing grid info
 
   character(len=240) :: filepath
-  type(vardesc)      :: vars(2)
-  type(fieldtype)    :: fields(2)
-  type(file_type)    :: IO_handle ! The I/O handle of the fileset
-  real               :: ds(GV%ke), dsi(GV%ke+1)
-
-  filepath    = trim(directory) // trim("Vertical_coordinate")
-  ds(:)       = getCoordinateResolution( CS%regridCS, undo_scaling=.true. )
-  dsi(1)      = 0.5*ds(1)
-  dsi(2:GV%ke) = 0.5*( ds(1:GV%ke-1) + ds(2:GV%ke) )
-  dsi(GV%ke+1) = 0.5*ds(GV%ke)
-
-  vars(1) = var_desc('ds', getCoordinateUnits( CS%regridCS ), &
-                    'Layer Coordinate Thickness','1','L','1')
-  vars(2) = var_desc('ds_interface', getCoordinateUnits( CS%regridCS ), &
-                    'Layer Center Coordinate Separation','1','i','1')
-
-  call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV)
-  call write_field(IO_handle, fields(1), ds)
-  call write_field(IO_handle, fields(2), dsi)
-  call close_file(IO_handle)
+
+  filepath = trim(directory) // trim("Vertical_coordinate")
+
+  call write_regrid_file(CS%regridCS, GV, filepath)
 
 end subroutine ALE_writeCoordinateFile
 

src/ALE/MOM_regridding.F90

@@ -6,6 +6,8 @@ module MOM_regridding
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, assert
 use MOM_file_parser,   only : param_file_type, get_param, log_param
 use MOM_io,            only : file_exists, field_exists, field_size, MOM_read_data
+use MOM_io,            only : vardesc, var_desc, fieldtype, SINGLE_FILE
+use MOM_io,            only : create_file, MOM_write_field, close_file, file_type
 use MOM_io,            only : verify_variable_units, slasher
 use MOM_unit_scaling,  only : unit_scale_type
 use MOM_variables,     only : ocean_grid_type, thermo_var_ptrs
@@ -129,9 +131,8 @@ module MOM_regridding
 ! The following routines are visible to the outside world
 public initialize_regridding, end_regridding, regridding_main
 public inflate_vanished_layers_old, check_remapping_grid, check_grid_column
-public set_regrid_params, get_regrid_size
+public set_regrid_params, get_regrid_size, write_regrid_file
 public uniformResolution, setCoordinateResolution
-public build_rho_column
 public set_target_densities_from_GV, set_target_densities
 public set_regrid_max_depths, set_regrid_max_thickness
 public getCoordinateResolution, getCoordinateInterfaces
@@ -2107,6 +2108,37 @@ subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
 end subroutine set_regrid_max_thickness
 
 
+!> Write the vertical coordinate information into a file.
+!! This subroutine writes out a file containing any available data related
+!! to the vertical grid used by the MOM ocean model when in ALE mode.
+subroutine write_regrid_file( CS, GV, filepath )
+  type(regridding_CS),     intent(in) :: CS        !< Regridding control structure
+  type(verticalGrid_type), intent(in) :: GV        !< ocean vertical grid structure
+  character(len=*),        intent(in) :: filepath  !< The full path to the file to write
+
+  type(vardesc)      :: vars(2)
+  type(fieldtype)    :: fields(2)
+  type(file_type)    :: IO_handle ! The I/O handle of the fileset
+  real               :: ds(GV%ke), dsi(GV%ke+1)
+
+  ds(:)       = CS%coord_scale * CS%coordinateResolution(:)
+  dsi(1)      = 0.5*ds(1)
+  dsi(2:GV%ke) = 0.5*( ds(1:GV%ke-1) + ds(2:GV%ke) )
+  dsi(GV%ke+1) = 0.5*ds(GV%ke)
+
+  vars(1) = var_desc('ds', getCoordinateUnits( CS ), &
+                     'Layer Coordinate Thickness', '1', 'L', '1')
+  vars(2) = var_desc('ds_interface', getCoordinateUnits( CS ), &
+                     'Layer Center Coordinate Separation', '1', 'i', '1')
+
+  call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV)
+  call MOM_write_field(IO_handle, fields(1), ds)
+  call MOM_write_field(IO_handle, fields(2), dsi)
+  call close_file(IO_handle)
+
+end subroutine write_regrid_file
+
+
 !------------------------------------------------------------------------------
 !> Query the fixed resolution data
 function getCoordinateResolution( CS, undo_scaling )

src/core/MOM.F90

@@ -57,7 +57,7 @@ module MOM
 use MOM_barotropic,            only : Barotropic_CS
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
 use MOM_check_scaling,         only : check_MOM6_scaling_factors
-use MOM_coord_initialization,  only : MOM_initialize_coord
+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_stochastics,           only : stochastics_init, update_stochastics, stochastic_CS
@@ -2484,8 +2484,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   ! Initialize dynamically evolving fields, perhaps from restart files.
   call cpu_clock_begin(id_clock_MOM_init)
-  call MOM_initialize_coord(GV, US, param_file, write_geom_files, &
-                            dirs%output_directory, CS%tv, G%max_depth)
+  call MOM_initialize_coord(GV, US, param_file, CS%tv, G%max_depth)
   call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
 
   if (CS%use_ALE_algorithm) then
@@ -2648,8 +2647,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     ! \todo This block exists for legacy reasons and we should phase it out of
     ! all examples. !###
     if (CS%debug) then
-      call uvchksum("Pre ALE adjust init cond [uv]", &
-                    CS%u, CS%v, G%HI, haloshift=1)
+      call uvchksum("Pre ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1)
       call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m)
     endif
     call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
@@ -2712,19 +2710,21 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   ! must be defined.
   call set_masks_for_axes(G, diag)
 
-  ! Diagnose static fields AND associate areas/volumes with axes
-  call write_static_fields(G, GV, US, CS%tv, CS%diag)
-  call callTree_waypoint("static fields written (initialize_MOM)")
-
   ! Register the volume cell measure (must be one of first diagnostics)
   call register_cell_measure(G, CS%diag, Time)
 
   call cpu_clock_begin(id_clock_MOM_init)
+  ! Diagnose static fields AND associate areas/volumes with axes
+  call write_static_fields(G, GV, US, CS%tv, CS%diag)
+  call callTree_waypoint("static fields written (initialize_MOM)")
+
   if (CS%use_ALE_algorithm) then
     call ALE_writeCoordinateFile( CS%ALE_CSp, GV, dirs%output_directory )
+    call callTree_waypoint("ALE initialized (initialize_MOM)")
+  elseif (write_geom_files) then
+    call write_vertgrid_file(GV, US, param_file, dirs%output_directory)
   endif
   call cpu_clock_end(id_clock_MOM_init)
-  call callTree_waypoint("ALE initialized (initialize_MOM)")
 
   CS%useMEKE = MEKE_init(Time, G, US, param_file, diag, CS%MEKE_CSp, CS%MEKE, restart_CSp)
 

src/initialization/MOM_coord_initialization.F90

@@ -9,8 +9,7 @@ module MOM_coord_initialization
 use MOM_error_handler,    only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser,      only : get_param, read_param, log_param, param_file_type, log_version
 use MOM_io,               only : close_file, create_file, file_type, fieldtype, file_exists
-use MOM_io,               only : MOM_read_data, MOM_write_field, vardesc, var_desc
-use MOM_io,               only : SINGLE_FILE, MULTIPLE
+use MOM_io,               only : MOM_read_data, MOM_write_field, vardesc, var_desc, SINGLE_FILE
 use MOM_string_functions, only : slasher, uppercase
 use MOM_unit_scaling,     only : unit_scale_type
 use MOM_variables,        only : thermo_var_ptrs
@@ -20,7 +19,7 @@ module MOM_coord_initialization
 
 implicit none ; private
 
-public MOM_initialize_coord
+public MOM_initialize_coord, write_vertgrid_file
 
 ! 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
@@ -33,13 +32,11 @@ module MOM_coord_initialization
 
 !> MOM_initialize_coord sets up time-invariant quantities related to MOM6's
 !!   vertical coordinate.
-subroutine MOM_initialize_coord(GV, US, PF, write_geom, output_dir, tv, max_depth)
+subroutine MOM_initialize_coord(GV, US, PF, tv, max_depth)
   type(verticalGrid_type), intent(inout) :: GV         !< Ocean vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US         !< A dimensional unit scaling type
   type(param_file_type),   intent(in)    :: PF         !< A structure indicating the open file
                                                        !! to parse for model parameter values.
-  logical,                 intent(in)    :: write_geom !< If true, write grid geometry files.
-  character(len=*),        intent(in)    :: output_dir !< The directory into which to write files.
   type(thermo_var_ptrs),   intent(inout) :: tv         !< The thermodynamic variable structure.
   real,                    intent(in)    :: max_depth  !< The ocean's maximum depth [Z ~> m].
   ! Local
@@ -107,12 +104,9 @@ subroutine MOM_initialize_coord(GV, US, PF, write_geom, output_dir, tv, max_dept
   if (debug) call chksum(US%m_to_Z*US%L_to_m**2*US%s_to_T**2*GV%g_prime(:), "MOM_initialize_coord: g_prime ", 1, nz)
   call setVerticalGridAxes( GV%Rlay, GV, scale=US%R_to_kg_m3 )
 
-! Copy the maximum depth across from the input argument
+  ! Copy the maximum depth across from the input argument
   GV%max_depth = max_depth
 
-! Write out all of the grid data used by this run.
-  if (write_geom) call write_vertgrid_file(GV, US, PF, output_dir)
-
   call callTree_leave('MOM_initialize_coord()')
 
 end subroutine MOM_initialize_coord

src/ocean_data_assim/MOM_oda_driver.F90

@@ -280,8 +280,7 @@ subroutine init_oda(Time, G, GV, diag_CS, CS)
   call clone_MOM_domain(CS%Grid%Domain, dG%Domain,symmetric=.false.)
   call set_grid_metrics(dG, PF, CS%US)
   call MOM_initialize_topography(dG%bathyT, dG%max_depth, dG, PF, CS%US)
-  call MOM_initialize_coord(CS%GV, CS%US, PF, .false., &
-           dirs%output_directory, tv_dummy, dG%max_depth)
+  call MOM_initialize_coord(CS%GV, CS%US, PF, tv_dummy, dG%max_depth)
   call ALE_init(PF, CS%GV, CS%US, dG%max_depth, CS%ALE_CS)
   call MOM_grid_init(CS%Grid, PF, global_indexing=.false.)
   call ALE_updateVerticalGridType(CS%ALE_CS, CS%GV)