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)