From 356671c202419346414f9c28661b3dab8017ccdc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 23 Mar 2022 15:52:49 -0400 Subject: [PATCH] +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 | 32 ++++------------- src/ALE/MOM_regridding.F90 | 36 +++++++++++++++++-- src/core/MOM.F90 | 20 +++++------ .../MOM_coord_initialization.F90 | 14 +++----- src/ocean_data_assim/MOM_oda_driver.F90 | 3 +- 5 files changed, 56 insertions(+), 49 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index d8eee55c14..fa195d57eb 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/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 diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 008475b16d..4c99ce457d 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/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 ) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5b35c01a42..f3d8869320 100644 --- a/src/core/MOM.F90 +++ b/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) diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index 286dfa7d95..311145bce1 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/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 diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index f183231c88..fc76c23480 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/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)