diff --git a/.testing/Makefile b/.testing/Makefile index f330c92e3b..4096436f30 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -265,9 +265,9 @@ build/openmp/Makefile: MOM_ACFLAGS=--enable-openmp build/target/Makefile: MOM_ACFLAGS= build/opt/Makefile: MOM_ACFLAGS= build/opt_target/Makefile: MOM_ACFLAGS= -build/coupled/Makefile: MOM_ACFLAGS=--with-driver=coupled_driver -build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_driver -build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_driver +build/coupled/Makefile: MOM_ACFLAGS=--with-driver=FMS_cap +build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_cap +build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_cap # Fetch regression target source code build/target/Makefile: | $(TARGET_CODEBASE) diff --git a/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 b/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 index e50f2ccf0b..b7ee7de684 100644 --- a/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 +++ b/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 @@ -11,9 +11,9 @@ module FMS_coupler_util !> Get element and index of a boundary condition subroutine extract_coupler_values(BC_struc, BC_index, BC_element, array_out, ilb, jlb, & is, ie, js, je, conversion) - real, dimension(ilb:,jlb:),intent(out) :: array_out !< The array being filled with the input values integer, intent(in) :: ilb !< Lower bounds integer, intent(in) :: jlb !< Lower bounds + real, dimension(ilb:,jlb:),intent(out) :: array_out !< The array being filled with the input values type(coupler_2d_bc_type), intent(in) :: BC_struc !< The type from which the data is being extracted integer, intent(in) :: BC_index !< The boundary condition number being extracted integer, intent(in) :: BC_element !< The element of the boundary condition being extracted @@ -27,9 +27,9 @@ end subroutine extract_coupler_values !> Set element and index of a boundary condition subroutine set_coupler_values(array_in, BC_struc, BC_index, BC_element, ilb, jlb,& is, ie, js, je, conversion) - real, dimension(ilb:,jlb:), intent(in) :: array_in !< The array containing the values to load into the BC integer, intent(in) :: ilb !< Lower bounds integer, intent(in) :: jlb !< Lower bounds + real, dimension(ilb:,jlb:), intent(in) :: array_in !< The array containing the values to load into the BC type(coupler_2d_bc_type), intent(inout) :: BC_struc !< The type into which the data is being loaded integer, intent(in) :: BC_index !< The boundary condition number being set integer, intent(in) :: BC_element !< The element of the boundary condition being set diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 index 6bd445ae8b..42c386497a 100644 --- a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 +++ b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 @@ -69,13 +69,13 @@ end subroutine generic_tracer_coupler_accumulate subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dtts,& grid_dat,model_time,nbands,max_wavelength_band,sw_pen_band,opacity_band,internal_heat,& frunoff,grid_ht, current_wave_stress, sosga) + integer, intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain + integer, intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain real, dimension(ilb:,jlb:,:), intent(in) :: Temp !< Potential temperature [deg C] real, dimension(ilb:,jlb:,:), intent(in) :: Salt !< Salinity [psu] real, dimension(ilb:,jlb:,:), intent(in) :: rho_dzt !< Mass per unit area of each layer [kg m-2] real, dimension(ilb:,jlb:,:), intent(in) :: dzt !< Ocean layer thickness [m] real, dimension(ilb:,jlb:), intent(in) :: hblt_depth !< Boundary layer depth [m] - integer, intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain - integer, intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain integer, intent(in) :: tau !< Time step index of %field real, intent(in) :: dtts !< The time step for this call [s] real, dimension(ilb:,jlb:), intent(in) :: grid_dat !< Grid cell areas [m2] diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index fa195d57eb..876c2c684b 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -20,6 +20,8 @@ module MOM_ALE use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_error_handler, only : callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint +use MOM_hybgen_unmix, only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS +use MOM_hybgen_regrid, only : hybgen_regrid_CS use MOM_file_parser, only : get_param, param_file_type, log_param use MOM_interface_heights,only : find_eta use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W @@ -27,6 +29,7 @@ module MOM_ALE use MOM_regridding, only : initialize_regridding, regridding_main, end_regridding use MOM_regridding, only : uniformResolution use MOM_regridding, only : inflate_vanished_layers_old +use MOM_regridding, only : regridding_preadjust_reqs, convective_adjustment use MOM_regridding, only : set_target_densities_from_GV, set_target_densities use MOM_regridding, only : regriddingCoordinateModeDoc, DEFAULT_COORDINATE_MODE use MOM_regridding, only : regriddingInterpSchemeDoc, regriddingDefaultInterpScheme @@ -73,6 +76,11 @@ module MOM_ALE type(remapping_CS) :: remapCS !< Remapping parameters and work arrays type(remapping_CS) :: vel_remapCS !< Remapping parameters for velocities and work arrays + type(hybgen_unmix_CS), pointer :: hybgen_unmixCS => NULL() !< Parameters for hybgen remapping + + logical :: use_hybgen_unmix !< If true, use the hybgen unmixing code before regridding + logical :: do_conv_adj !< If true, do convective adjustment before regridding + integer :: nk !< Used only for queries, not directly by this module real :: BBL_h_vel_mask !< The thickness of a bottom boundary layer within which velocities in !! thin layers are zeroed out after remapping, following practice with @@ -119,7 +127,6 @@ module MOM_ALE public ALE_main_offline public ALE_offline_inputs public ALE_offline_tracer_final -public ALE_build_grid public ALE_regrid_accelerated public ALE_remap_scalar public ALE_PLM_edge_values @@ -165,6 +172,8 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) logical :: force_bounds_in_subcell logical :: local_logical logical :: remap_boundary_extrap + type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding + ! for sharing parameters. if (associated(CS)) then call MOM_error(WARNING, "ALE_init called with an associated "// & @@ -184,6 +193,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) ! Initialize and configure regridding call ALE_initRegridding(GV, US, max_depth, param_file, mdl, CS%regridCS) + call regridding_preadjust_reqs(CS%regridCS, CS%do_conv_adj, CS%use_hybgen_unmix, hybgen_CS=hybgen_regridCS) ! Initialize and configure remapping that is orchestrated by ALE. call get_param(param_file, mdl, "REMAPPING_SCHEME", string, & @@ -277,6 +287,9 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) "to avoid such filtering altogether.", & default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=(CS%BBL_h_vel_mask<=0.0)) + if (CS%use_hybgen_unmix) & + call init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS) + call get_param(param_file, "MOM", "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) @@ -350,6 +363,7 @@ subroutine ALE_end(CS) ! Deallocate memory used for the regridding call end_remapping( CS%remapCS ) + if (CS%use_hybgen_unmix) call end_hybgen_unmix( CS%hybgen_unmixCS ) call end_regridding( CS%regridCS ) deallocate(CS) @@ -379,9 +393,9 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2] logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) !< If true, PCM remapping should be used in a cell. - integer :: nk, i, j, k, isc, iec, jsc, jec, ntr + integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90") @@ -401,10 +415,19 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) endif dzRegrid(:,:,:) = 0.0 + ! If necessary, do some preparatory work to clean up the model state before regridding. + + ! This adjusts the input thicknesses prior to remapping, based on the verical coordinate. + if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv) + if (CS%use_hybgen_unmix) then + ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr + call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) + endif + ! Build new grid. The new grid is stored in h_new. The old grid is h. ! Both are needed for the subsequent remapping of variables. - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, & - frac_shelf_h, PCM_cell=PCM_cell) + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false., & + frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell) call check_grid( G, GV, h, 0. ) @@ -459,9 +482,9 @@ subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, OBC, dt) ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2] - integer :: nk, i, j, k, isc, iec, jsc, jec + integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90") @@ -470,9 +493,16 @@ subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, OBC, dt) endif dzRegrid(:,:,:) = 0.0 + ! This adjusts the input state prior to remapping, depending on the verical coordinate. + if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv) + if (CS%use_hybgen_unmix) then + ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr + call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) + endif + ! Build new grid. The new grid is stored in h_new. The old grid is h. ! Both are needed for the subsequent remapping of variables. - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid ) + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false. ) call check_grid( G, GV, h, 0. ) @@ -519,7 +549,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: temp_vec ! Transports on the destination grid [H L2 ~> m3 or kg] - nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke dzRegrid(:,:,:) = 0.0 h_new(:,:,:) = 0.0 @@ -595,13 +625,18 @@ subroutine ALE_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS, OBC) real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid !< The change in grid interface positions real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new !< Regridded target thicknesses - integer :: nk, i, j, k, isc, iec, jsc, jec + integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90") ! Need to make sure that h_target is consistent with the current offline ALE confiuration - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid ) + if (CS%do_conv_adj) call convective_adjustment(G, GV, h_target, tv) + if (CS%use_hybgen_unmix) then + ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr + call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) + endif + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid, conv_adjust=.false. ) call check_grid( G, GV, h_target, 0. ) @@ -651,52 +686,17 @@ subroutine check_grid( G, GV, h, threshold ) end subroutine check_grid -!### This routine does not appear to be used. -!> Generates new grid -subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h ) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(regridding_CS), intent(in) :: regridCS !< Regridding parameters and options - type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variable structure - real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the - !! last time step [H ~> m or kg m-2] - logical, optional, intent(in) :: debug !< If true, show the call tree - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in):: frac_shelf_h !< Fractional ice shelf coverage [nondim] - ! Local variables - integer :: nk, i, j, k - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! The new grid thicknesses - logical :: show_call_tree - - show_call_tree = .false. - if (present(debug)) show_call_tree = debug - if (show_call_tree) call callTree_enter("ALE_build_grid(), MOM_ALE.F90") - - ! Build new grid. The new grid is stored in h_new. The old grid is h. - ! Both are needed for the subsequent remapping of variables. - call regridding_main( remapCS, regridCS, G, GV, h, tv, h_new, dzRegrid, frac_shelf_h ) - - ! Override old grid with new one. The new grid 'h_new' is built in - ! one of the 'build_...' routines above. -!$OMP parallel do default(none) shared(G,h,h_new) - do j = G%jsc,G%jec ; do i = G%isc,G%iec - if (G%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:) - enddo ; enddo - - if (show_call_tree) call callTree_leave("ALE_build_grid()") -end subroutine ALE_build_grid !> For a state-based coordinate, accelerate the process of regridding by !! repeatedly applying the grid calculation algorithm -subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial) +subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, dzRegrid, initial) type(ALE_CS), pointer :: CS !< ALE control structure type(ocean_grid_type), intent(inout) :: G !< Ocean grid type(verticalGrid_type), intent(in) :: GV !< Vertical grid real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Original thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS) - integer, intent(in) :: n !< Number of times to regrid + integer, intent(in) :: n_itt !< Number of times to regrid real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -711,7 +711,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzReg !! routine (and expect diagnostics to work) ! Local variables - integer :: i, j, k, nz + integer :: i, j, itt, nz type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt type(group_pass_type) :: pass_T_S_h ! group pass if the coordinate has a stencil real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig ! A working copy of layer thicknesses @@ -744,11 +744,12 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzReg if (present(dt)) & call ALE_update_regrid_weights(dt, CS) - do k = 1, n + do itt = 1, n_itt call do_group_pass(pass_T_S_h, G%domain) ! generate new grid - call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface) + if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local) + call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface, conv_adjust=.false.) dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) ! remap from original grid onto new grid diff --git a/src/ALE/MOM_hybgen_regrid.F90 b/src/ALE/MOM_hybgen_regrid.F90 new file mode 100644 index 0000000000..783820829f --- /dev/null +++ b/src/ALE/MOM_hybgen_regrid.F90 @@ -0,0 +1,983 @@ +!> This module contains the hybgen regridding routines from HYCOM, with minor +!! modifications to follow the MOM6 coding conventions +module MOM_hybgen_regrid + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_EOS, only : EOS_type, calculate_density +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, assert +use MOM_file_parser, only : get_param, param_file_type, log_param +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, SINGLE_FILE +use MOM_string_functions, only : slasher +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : ocean_grid_type, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +!> Control structure containing required parameters for the hybgen coordinate generator +type, public :: hybgen_regrid_CS ; private + + real :: min_thickness !< Minimum thickness allowed for layers [H ~> m or kg m-2] + + integer :: nk !< Number of layers on the target grid + + !> Reference pressure for density calculations [R L2 T-2 ~> Pa] + real :: ref_pressure + + !> Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3] + real :: hybiso + !> Number of sigma levels used by HYBGEN + integer :: nsigma + + real :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2] + real :: qhybrlx !< Fractional relaxation within a regridding step [nondim] + + real, allocatable, dimension(:) :: & + dp0k, & !< minimum deep z-layer separation [H ~> m or kg m-2] + ds0k !< minimum shallow z-layer separation [H ~> m or kg m-2] + + real :: coord_scale = 1.0 !< A scaling factor to restores the depth coordinates to values in m + real :: Rho_coord_scale = 1.0 !< A scaling factor to restores the denesity coordinates to values in kg m-3 + + real :: dpns !< depth to start terrain following [H ~> m or kg m-2] + real :: dsns !< depth to stop terrain following [H ~> m or kg m-2] + real :: min_dilate !< The minimum amount of dilation that is permitted when converting target + !! coordinates from z to z* [nondim]. This limit applies when wetting occurs. + real :: max_dilate !< The maximum amount of dilation that is permitted when converting target + !! coordinates from z to z* [nondim]. This limit applies when drying occurs. + + real :: thkbot !< Thickness of a bottom boundary layer, within which hybgen does + !! something different. [H ~> m or kg m-2] + + !> Shallowest depth for isopycnal layers [H ~> m or kg m-2] + real :: topiso_const + ! real, dimension(:,:), allocatable :: topiso + + !> Nominal density of interfaces [R ~> kg m-3] + real, allocatable, dimension(:) :: target_density + + real :: onem !< Nominally one m in thickness units [H ~> m or kg m-2] + +end type hybgen_regrid_CS + + +public hybgen_regrid, init_hybgen_regrid, end_hybgen_regrid +public hybgen_column_init, get_hybgen_regrid_params, write_Hybgen_coord_file + +contains + +!> Initialise a hybgen_regrid_CS control structure and store its parameters +subroutine init_hybgen_regrid(CS, GV, US, param_file) + type(hybgen_regrid_CS), pointer :: CS !< Unassociated pointer to hold the control structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file + + character(len=40) :: mdl = "MOM_hybgen_regrid" ! This module's name. + real :: hybrlx ! The number of remappings over which to move toward the target coordinate [timesteps] + character(len=40) :: dp0_coord_var, ds0_coord_var, rho_coord_var + character(len=200) :: filename, coord_file, inputdir ! Strings for file/path + logical :: use_coord_file + integer :: k + + if (associated(CS)) call MOM_error(FATAL, "init_hybgen_regrid: CS already associated!") + allocate(CS) + + CS%nk = GV%ke + + allocate(CS%target_density(CS%nk)) + allocate(CS%dp0k(CS%nk), source=0.0) ! minimum deep z-layer separation + allocate(CS%ds0k(CS%nk), source=0.0) ! minimum shallow z-layer separation + + do k=1,CS%nk ; CS%target_density(k) = GV%Rlay(k) ; enddo + + call get_param(param_file, mdl, "P_REF", CS%ref_pressure, & + "The pressure that is used for calculating the coordinate "//& + "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//& + "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", & + units="Pa", default=2.0e7, scale=US%kg_m3_to_R*US%m_s_to_L_T**2) + + call get_param(param_file, mdl, "HYBGEN_MIN_THICKNESS", CS%min_thickness, & + "The minimum layer thickness allowed when regridding with Hybgen.", & + units="m", default=0.0, scale=GV%m_to_H ) + + call get_param(param_file, mdl, "HYBGEN_N_SIGMA", CS%nsigma, & + "The number of sigma-coordinate (terrain-following) layers with Hybgen regridding.", & + default=0) + call get_param(param_file, mdl, "HYBGEN_COORD_FILE", coord_file, & + "The file from which the Hybgen profile is read, or blank to use a list of "//& + "real input parameters from the MOM_input file.", default="") + + use_coord_file = (len_trim(coord_file) > 0) + call get_param(param_file, mdl, "HYBGEN_DEEP_DZ_PR0FILE", CS%dp0k, & + "The layerwise list of deep z-level minimum thicknesses for Hybgen (dp0k in Hycom).", & + units="m", default=0.0, scale=GV%m_to_H, do_not_log=use_coord_file) + call get_param(param_file, mdl, "HYBGEN_SHALLOW_DZ_PR0FILE", CS%ds0k, & + "The layerwise list of shallow z-level minimum thicknesses for Hybgen (ds0k in Hycom).", & + units="m", default=0.0, scale=GV%m_to_H, do_not_log=use_coord_file) + + if (use_coord_file) then + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + filename = trim(inputdir)//trim(coord_file) + call log_param(param_file, mdl, "INPUTDIR/HYBGEN_COORD_FILE", filename) + if (.not.file_exists(filename)) call MOM_error(FATAL, & + " set_coord_from_file: Unable to open "//trim(filename)) + + call get_param(param_file, mdl, "HYBGEN_DEEP_DZ_VAR", dp0_coord_var, & + "The variable in HYBGEN_COORD_FILE that is to be used for the "//& + "deep z-level minimum thicknesses for Hybgen (dp0k in Hycom).", & + default="dp0") + call get_param(param_file, mdl, "HYBGEN_SHALLOW_DZ_VAR", ds0_coord_var, & + "The variable in HYBGEN_COORD_FILE that is to be used for the "//& + "shallow z-level minimum thicknesses for Hybgen (ds0k in Hycom).", & + default="ds0") + call get_param(param_file, mdl, "HYBGEN_TGT_DENSITY_VAR", rho_coord_var, & + "The variable in HYBGEN_COORD_FILE that is to be used for the Hybgen "//& + "target layer densities, or blank to reuse the values in GV%Rlay.", & + default="") + + call MOM_read_data(filename, dp0_coord_var, CS%dp0k, scale=GV%m_to_H) + + call MOM_read_data(filename, ds0_coord_var, CS%ds0k, scale=GV%m_to_H) + + if (len_trim(rho_coord_var) > 0) & + call MOM_read_data(filename, rho_coord_var, CS%target_density, scale=US%kg_m3_to_R) + endif + + call get_param(param_file, mdl, "HYBGEN_ISOPYCNAL_DZ_MIN", CS%dp00i, & + "The Hybgen deep isopycnal spacing minimum thickness (dp00i in Hycom)", & + units="m", default=0.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_MIN_ISO_DEPTH", CS%topiso_const, & + "The Hybgen shallowest depth for isopycnal layers (isotop in Hycom)", & + units="m", default=0.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_RELAX_PERIOD", hybrlx, & + "The Hybgen coordinate relaxation period in timesteps, or 1 to move to "//& + "the new target coordinates in a single step. This must be >= 1.", & + units="timesteps", default=1.0) + if (hybrlx < 1.0) call MOM_error(FATAL, "init_hybgen_regrid: HYBGEN_RELAX_PERIOD must be at least 1.") + CS%qhybrlx = 1.0 / hybrlx + call get_param(param_file, mdl, "HYBGEN_BBL_THICKNESS", CS%thkbot, & + "A bottom boundary layer thickness within which Hybgen is able to move "//& + "overlying layers upward to match a target density.", & + units="m", default=0.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_REMAP_DENSITY_MATCH", CS%hybiso, & + "A tolerance between the layer densities and their target, within which "//& + "Hybgen determines that remapping uses PCM for a layer.", & + units="kg m-3", default=0.0, scale=US%kg_m3_to_R) + call get_param(param_file, mdl, "HYBGEN_REMAP_MIN_ZSTAR_DILATE", CS%min_dilate, & + "The maximum amount of dilation that is permitted when converting target "//& + "coordinates from z to z* [nondim]. This limit applies when drying occurs.", & + default=0.5) + call get_param(param_file, mdl, "HYBGEN_REMAP_MAX_ZSTAR_DILATE", CS%max_dilate, & + "The maximum amount of dilation that is permitted when converting target "//& + "coordinates from z to z* [nondim]. This limit applies when drying occurs.", & + default=2.0) + + CS%onem = 1.0 * GV%m_to_H + + do k=1,CS%nk ; CS%dp0k(k) = max(CS%dp0k(k), CS%min_thickness) ; enddo + CS%dp00i = max(CS%dp00i, CS%min_thickness) + + ! Determine the depth range over which to use a sigma (terrain-following) coordinate. + ! --- terrain following starts at depth dpns and ends at depth dsns + if (CS%nsigma == 0) then + CS%dpns = CS%dp0k(1) + CS%dsns = 0.0 + else + CS%dpns = 0.0 + CS%dsns = 0.0 + do k=1,CS%nsigma + CS%dpns = CS%dpns + CS%dp0k(k) + CS%dsns = CS%dsns + CS%ds0k(k) + enddo !k + endif !nsigma + + CS%coord_scale = GV%H_to_m + CS%Rho_coord_scale = US%R_to_kg_m3 + +end subroutine init_hybgen_regrid + +!> Writes out a file containing any available data related +!! to the vertical grid used by the MOM ocean model. +subroutine write_Hybgen_coord_file(GV, CS, filepath) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(hybgen_regrid_CS), intent(in) :: CS !< Control structure for this module + character(len=*), intent(in) :: filepath !< The full path to the file to write + ! Local variables + type(vardesc) :: vars(3) + type(fieldtype) :: fields(3) + type(file_type) :: IO_handle ! The I/O handle of the fileset + + vars(1) = var_desc("dp0", "meter", "Deep z-level minimum thicknesses for Hybgen", '1', 'L', '1') + vars(2) = var_desc("ds0", "meter", "Shallow z-level minimum thicknesses for Hybgen", '1', 'L', '1') + vars(3) = var_desc("Rho_tgt", "kg m-3", "Target coordinate potential densities for Hybgen", '1', 'L', '1') + call create_file(IO_handle, trim(filepath), vars, 3, fields, SINGLE_FILE, GV=GV) + + call MOM_write_field(IO_handle, fields(1), CS%dp0k, scale=CS%coord_scale) + call MOM_write_field(IO_handle, fields(2), CS%ds0k, scale=CS%coord_scale) + call MOM_write_field(IO_handle, fields(3), CS%target_density, scale=CS%Rho_coord_scale) + + call close_file(IO_handle) + +end subroutine write_Hybgen_coord_file + +!> This subroutine deallocates memory in the control structure for the hybgen module +subroutine end_hybgen_regrid(CS) + type(hybgen_regrid_CS), pointer :: CS !< Coordinate control structure + + ! nothing to do + if (.not. associated(CS)) return + + deallocate(CS%target_density) + deallocate(CS%dp0k, CS%ds0k) + deallocate(CS) +end subroutine end_hybgen_regrid + +!> This subroutine can be used to retrieve the parameters for the hybgen regrid module +subroutine get_hybgen_regrid_params(CS, nk, ref_pressure, hybiso, nsigma, dp00i, qhybrlx, & + dp0k, ds0k, dpns, dsns, min_dilate, max_dilate, & + thkbot, topiso_const, target_density) + type(hybgen_regrid_CS), pointer :: CS !< Coordinate regridding control structure + integer, optional, intent(out) :: nk !< Number of layers on the target grid + real, optional, intent(out) :: ref_pressure !< Reference pressure for density calculations [R L2 T-2 ~> Pa] + real, optional, intent(out) :: hybiso !< Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3] + integer, optional, intent(out) :: nsigma !< Number of sigma levels used by HYBGEN + real, optional, intent(out) :: dp00i !< Deep isopycnal spacing minimum thickness (m) + real, optional, intent(out) :: qhybrlx !< Fractional relaxation amount per timestep, 0 < qyhbrlx <= 1 [nondim] + real, optional, intent(out) :: dp0k(:) !< minimum deep z-layer separation [H ~> m or kg m-2] + real, optional, intent(out) :: ds0k(:) !< minimum shallow z-layer separation [H ~> m or kg m-2] + real, optional, intent(out) :: dpns !< depth to start terrain following [H ~> m or kg m-2] + real, optional, intent(out) :: dsns !< depth to stop terrain following [H ~> m or kg m-2] + real, optional, intent(out) :: min_dilate !< The minimum amount of dilation that is permitted when + !! converting target coordinates from z to z* [nondim]. + !! This limit applies when wetting occurs. + real, optional, intent(out) :: max_dilate !< The maximum amount of dilation that is permitted when + !! converting target coordinates from z to z* [nondim]. + !! This limit applies when drying occurs. + real, optional, intent(out) :: thkbot !< Thickness of a bottom boundary layer, within which + !! hybgen does something different. [H ~> m or kg m-2] + real, optional, intent(out) :: topiso_const !< Shallowest depth for isopycnal layers [H ~> m or kg m-2] + ! real, dimension(:,:), allocatable :: topiso + real, optional, intent(out) :: target_density(:) !< Nominal density of interfaces [R ~> kg m-3] + + if (.not. associated(CS)) call MOM_error(FATAL, "get_hybgen_params: CS not associated") + + if (present(nk)) nk = CS%nk + if (present(ref_pressure)) ref_pressure = CS%ref_pressure + if (present(hybiso)) hybiso = CS%hybiso + if (present(nsigma)) nsigma = CS%nsigma + if (present(dp00i)) dp00i = CS%dp00i + if (present(qhybrlx)) qhybrlx = CS%qhybrlx + if (present(dp0k)) then + if (size(dp0k) < CS%nk) call MOM_error(FATAL, "get_hybgen_regrid_params: "//& + "The dp0k argument is not allocated with enough space.") + dp0k(1:CS%nk) = CS%dp0k(1:CS%nk) + endif + if (present(ds0k)) then + if (size(ds0k) < CS%nk) call MOM_error(FATAL, "get_hybgen_regrid_params: "//& + "The ds0k argument is not allocated with enough space.") + ds0k(1:CS%nk) = CS%ds0k(1:CS%nk) + endif + if (present(dpns)) dpns = CS%dpns + if (present(dsns)) dsns = CS%dsns + if (present(min_dilate)) min_dilate = CS%min_dilate + if (present(max_dilate)) max_dilate = CS%max_dilate + if (present(thkbot)) thkbot = CS%thkbot + if (present(topiso_const)) topiso_const = CS%topiso_const + if (present(target_density)) then + if (size(target_density) < CS%nk) call MOM_error(FATAL, "get_hybgen_regrid_params: "//& + "The target_density argument is not allocated with enough space.") + target_density(1:CS%nk) = CS%target_density(1:CS%nk) + endif + +end subroutine get_hybgen_regrid_params + + +!> Modify the input grid to give a new vertical grid based on the HYCOM hybgen code. +subroutine hybgen_regrid(G, GV, US, dp, tv, CS, dzInterface, PCM_cell) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: dp !< Source grid layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + type(hybgen_regrid_CS), intent(in) :: CS !< hybgen control structure + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), & + intent(inout) :: dzInterface !< The change in height of each interface, + !! using a sign convention opposite to the change + !! in pressure [H ~> m or kg m-2] + logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(inout) :: PCM_cell !< If true, PCM remapping should be used in a cell. + !! This is effectively intent out, but values in wide + !! halo regions and land points are reused. + + ! --- ------------------------------------- + ! --- hybrid grid generator from HYCOM + ! --- ------------------------------------- + + ! These notes on the parameters for the hybrid grid generator are inhereted from the + ! Hycom source code for these algorithms. + ! + ! From blkdat.input (units may have changed from m to pressure): + ! + ! --- 'nsigma' = number of sigma levels + ! --- 'dp0k ' = layer k deep z-level spacing minimum thickness (m) + ! --- k=1,nk + ! --- 'ds0k ' = layer k shallow z-level spacing minimum thickness (m) + ! --- k=1,nsigma + ! --- 'dp00i' = deep isopycnal spacing minimum thickness (m) + ! --- 'isotop' = shallowest depth for isopycnal layers (m) + ! now in topiso(:,:) + ! --- 'sigma ' = isopycnal layer target densities (sigma units) + ! --- now in theta(:,:,1:nk) + ! + ! --- the above specifies a vertical coord. that is isopycnal or: + ! --- near surface z in deep water, based on dp0k + ! --- near surface z in shallow water, based on ds0k and nsigma + ! --- terrain-following between them, based on ds0k and nsigma + ! + ! --- terrain following starts at depth dpns=sum(dp0k(k),k=1,nsigma) and + ! --- ends at depth dsns=sum(ds0k(k),k=1,nsigma), and the depth of the + ! --- k-th layer interface varies linearly with total depth between + ! --- these two reference depths, i.e. a z-sigma-z fixed coordinate. + ! + ! --- near the surface (i.e. shallower than isotop), layers are always + ! --- fixed depth (z or sigma). + ! -- layer 1 is always fixed, so isotop=0.0 is not realizable. + ! --- near surface layers can also be forced to be fixed depth + ! --- by setting target densities (sigma(k)) very small. + ! + ! --- away from the surface, the minimum layer thickness is dp00i. + ! + ! --- for fixed depth targets to be: + ! --- z-only set nsigma=0, + ! --- sigma-z (shallow-deep) use a very small ds0k(:), + ! --- sigma-only set nsigma=nk, dp0k large, and ds0k small. + + ! These arrays work with the input column + real :: p_col(GV%ke) ! A column of reference pressures [R L2 T-2 ~> Pa] + real :: temp_in(GV%ke) ! A column of input potential temperatures [degC] + real :: saln_in(GV%ke) ! A column of input layer salinities [ppt] + real :: Rcv_in(GV%ke) ! An input column of coordinate potential density [R ~> kg m-3] + real :: dp_in(GV%ke) ! The input column of layer thicknesses [H ~> m or kg m-2] + logical :: PCM_lay(GV%ke) ! If true for a layer, use PCM remapping for that layer + + ! These arrays are on the target grid. + real :: Rcv_tgt(CS%nk) ! Target potential density [R ~> kg m-3] + real :: Rcv(CS%nk) ! Initial values of coordinate potential density on the target grid [R ~> kg m-3] + real :: h_col(CS%nk) ! A column of layer thicknesses [H ~> m or kg m-2] + real :: dz_int(CS%nk+1) ! The change in interface height due to remapping [H ~> m or kg m-2] + real :: Rcv_integral ! Integrated coordinate potential density in a layer [R H ~> kg m-2 or kg2 m-5] + + real :: qhrlx(CS%nk+1) ! Fractional relaxation within a timestep (between 0 and 1) [nondim] + real :: dp0ij(CS%nk) ! minimum layer thickness [H ~> m or kg m-2] + real :: dp0cum(CS%nk+1) ! minimum interface depth [H ~> m or kg m-2] + + real :: h_tot ! Total thickness of the water column [H ~> m or kg m-2] + real :: nominalDepth ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2] + real :: dilate ! A factor by which to dilate the target positions from z to z* [nondim] + integer :: fixlay ! Deepest fixed coordinate layer + integer, dimension(0:CS%nk) :: k_end ! The index of the deepest source layer that contributes to + ! each target layer, in the unusual case where the the input grid is + ! larger than the new grid. This situation only occurs during certain + ! types of initialization or when generating output diagnostics. + integer :: i, j, k, nk, m, k2, nk_in + + nk = CS%nk + + p_col(:) = CS%ref_pressure + + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then + + ! Store one-dimensional arrays of thicknesses for the 'old' vertical grid before regridding + h_tot = 0.0 + do K=1,GV%ke + temp_in(k) = tv%T(i,j,k) + saln_in(k) = tv%S(i,j,k) + dp_in(k) = dp(i,j,k) + h_tot = h_tot + dp_in(k) + enddo + + ! This sets the input column's coordinate potential density from T and S. + call calculate_density(temp_in, saln_in, p_col, Rcv_in, tv%eqn_of_state) + + ! Set the initial properties on the new grid from the old grid. + nk_in = GV%ke + if (GV%ke > CS%nk) then ; do k=GV%ke,CS%nk+1,-1 + ! Remove any excess massless layers from the bottom of the input column. + if (dp_in(k) > 0.0) exit + nk_in = k-1 + enddo ; endif + + if (CS%nk >= nk_in) then + ! Simply copy over the common layers. This is the usual case. + do k=1,min(CS%nk,GV%ke) + h_col(k) = dp_in(k) + Rcv(k) = Rcv_in(k) + enddo + if (CS%nk > GV%ke) then + ! Pad out the input column with additional massless layers with the bottom properties. + ! This case only occurs during initialization or perhaps when writing diagnostics. + do k=GV%ke+1,CS%nk + Rcv(k) = Rcv_in(GV%ke) + h_col(k) = 0.0 + enddo + endif + else ! (CS%nk < nk_in) + ! The input column has more data than the output. For now, combine layers to + ! make them the same size, but there may be better approaches that should be taken. + ! This case only occurs during initialization or perhaps when writing diagnostics. + ! This case was not handled by the original Hycom code in hybgen.F90. + do k=0,CS%nk ; k_end(k) = (k * nk_in) / CS%nk ; enddo + do k=1,CS%nk + h_col(k) = 0.0 ; Rcv_integral = 0.0 + do k2=k_end(k-1) + 1,k_end(k) + h_col(k) = h_col(k) + dp_in(k2) + Rcv_integral = Rcv_integral + dp_in(k2)*Rcv_in(k2) + enddo + if (h_col(k) > GV%H_subroundoff) then + ! Take the volume-weighted average properties. + Rcv(k) = Rcv_integral / h_col(k) + else ! Take the properties of the topmost source layer that contributes. + Rcv(k) = Rcv_in(k_end(k-1) + 1) + endif + enddo + endif + + ! Set the target densities for the new layers. + do k=1,CS%nk + ! Rcv_tgt(k) = theta(i,j,k) ! If a 3-d target density were set up in theta, use that here. + Rcv_tgt(k) = CS%target_density(k) ! MOM6 does not yet support 3-d target densities. + enddo + + ! The following block of code is used to trigger z* stretching of the targets heights. + nominalDepth = (G%bathyT(i,j) + G%Z_ref)*GV%Z_to_H + if (h_tot <= CS%min_dilate*nominalDepth) then + dilate = CS%min_dilate + elseif (h_tot >= CS%max_dilate*nominalDepth) then + dilate = CS%max_dilate + else + dilate = h_tot / nominalDepth + endif + + ! Convert the regridding parameters into specific constraints for this column. + call hybgen_column_init(nk, CS%nsigma, CS%dp0k, CS%ds0k, CS%dp00i, & + CS%topiso_const, CS%qhybrlx, CS%dpns, CS%dsns, h_tot, dilate, & + h_col, fixlay, qhrlx, dp0ij, dp0cum) + + ! Determine whether to require the use of PCM remapping from each source layer. + do k=1,GV%ke + if (CS%hybiso > 0.0) then + ! --- thin or isopycnal source layers are remapped with PCM. + PCM_lay(k) = (k > fixlay) .and. (abs(Rcv(k) - Rcv_tgt(k)) < CS%hybiso) + else ! hybiso==0.0, so purely isopycnal layers use PCM + PCM_lay(k) = .false. + endif ! hybiso + enddo !k + + ! Determine the new layer thicknesses. + call hybgen_column_regrid(CS, nk, CS%thkbot, CS%onem, & + 1.0e-11*US%kg_m3_to_R, Rcv_tgt, fixlay, qhrlx, dp0ij, & + dp0cum, Rcv, h_col, dz_int) + + ! Store the output from hybgenaij_regrid in 3-d arrays. + if (present(PCM_cell)) then ; do k=1,GV%ke + PCM_cell(i,j,k) = PCM_lay(k) + enddo ; endif + + do K=1,nk+1 + ! Note that dzInterface uses the opposite sign convention from the change in p. + dzInterface(i,j,K) = -dz_int(K) + enddo + + else + if (present(PCM_cell)) then ; do k=1,GV%ke + PCM_cell(i,j,k) = .false. + enddo ; endif + do k=1,CS%nk+1 ; dzInterface(i,j,k) = 0.0 ; enddo + endif ; enddo ; enddo !i & j. + +end subroutine hybgen_regrid + +!> Initialize some of the variables that are used for regridding or unmixing, including the +!! stretched contraits on where the new interfaces can be. +subroutine hybgen_column_init(nk, nsigma, dp0k, ds0k, dp00i, topiso_i_j, & + qhybrlx, dpns, dsns, h_tot, dilate, h_col, & + fixlay, qhrlx, dp0ij, dp0cum) + integer, intent(in) :: nk !< The number of layers in the new grid + integer, intent(in) :: nsigma !< The number of sigma levels + real, intent(in) :: dp0k(nk) !< Layer deep z-level spacing minimum thicknesses [H ~> m or kg m-2] + real, intent(in) :: ds0k(nsigma) !< Layer shallow z-level spacing minimum thicknesses [H ~> m or kg m-2] + real, intent(in) :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2] + real, intent(in) :: topiso_i_j !< Shallowest depth for isopycnal layers [H ~> m or kg m-2] + real, intent(in) :: qhybrlx !< Fractional relaxation amount per timestep, 0 < qyhbrlx <= 1 [nondim] + real, intent(in) :: h_tot !< The sum of the initial layer thicknesses [H ~> m or kg m-2] + real, intent(in) :: dilate !< A factor by which to dilate the target positions + !! from z to z* [nondim] + real, intent(in) :: h_col(nk) !< Initial layer thicknesses [H ~> m or kg m-2] + real, intent(in) :: dpns !< Vertical sum of dp0k [H ~> m or kg m-2] + real, intent(in) :: dsns !< Vertical sum of ds0k [H ~> m or kg m-2] + integer, intent(out) :: fixlay !< Deepest fixed coordinate layer + real, intent(out) :: qhrlx(nk+1) !< Fractional relaxation within a timestep (between 0 and 1) [nondim] + real, intent(out) :: dp0ij(nk) !< minimum layer thickness [H ~> m or kg m-2] + real, intent(out) :: dp0cum(nk+1) !< minimum interface depth [H ~> m or kg m-2] + + ! --- -------------------------------------------------------------- + ! --- hybrid grid generator, single column - initialization. + ! --- -------------------------------------------------------------- + + ! Local variables + character(len=256) :: mesg ! A string for output messages + real :: hybrlx ! The relaxation rate in the hybrid region [timestep-1]? + real :: qdep ! Total water column thickness as a fraction of dp0k (vs ds0k) [nondim] + real :: q ! A portion of the thickness that contributes to the new cell [H ~> m or kg m-2] + real :: p_int(nk+1) ! Interface depths [H ~> m or kg m-2] + integer :: k, fixall + + ! --- dpns = sum(dp0k(k),k=1,nsigma) + ! --- dsns = sum(ds0k(k),k=1,nsigma) + ! --- terrain following starts (on the deep side) at depth dpns and ends (on the + ! --- shallow side) at depth dsns and the depth of the k-th layer interface varies + ! --- linearly with total depth between these two reference depths. + if ((h_tot >= dilate * dpns) .or. (dpns <= dsns)) then + qdep = 1.0 ! Not terrain following - this column is too thick or terrain following is disabled. + elseif (h_tot <= dilate * dsns) then + qdep = 0.0 ! Not terrain following - this column is too thin + else + qdep = (h_tot - dilate * dsns) / (dilate * (dpns - dsns)) + endif + + if (qdep < 1.0) then + ! Terrain following or shallow fixed coordinates, qhrlx=1 and ignore dp00 + p_int( 1) = 0.0 + dp0cum(1) = 0.0 + qhrlx( 1) = 1.0 + dp0ij( 1) = dilate * (qdep*dp0k(1) + (1.0-qdep)*ds0k(1)) + + dp0cum(2) = dp0cum(1) + dp0ij(1) + qhrlx( 2) = 1.0 + p_int( 2) = p_int(1) + h_col(1) + do k=2,nk + qhrlx( k+1) = 1.0 + dp0ij( k) = dilate * (qdep*dp0k(k) + (1.0-qdep)*ds0k(k)) + dp0cum(k+1) = dp0cum(k) + dp0ij(k) + p_int( k+1) = p_int(k) + h_col(k) + enddo !k + else + ! Not terrain following + p_int( 1) = 0.0 + dp0cum(1) = 0.0 + qhrlx( 1) = 1.0 !no relaxation in top layer + dp0ij( 1) = dilate * dp0k(1) + + dp0cum(2) = dp0cum(1) + dp0ij(1) + qhrlx( 2) = 1.0 !no relaxation in top layer + p_int( 2) = p_int(1) + h_col(1) + do k=2,nk + if ((dp0k(k) <= dp00i) .or. (dilate * dp0k(k) >= p_int(k) - dp0cum(k))) then + ! This layer is in fixed surface coordinates. + dp0ij(k) = dp0k(k) + qhrlx(k+1) = 1.0 + else + q = dp0k(k) * (dilate * dp0k(k) / ( p_int(k) - dp0cum(k)) ) ! A fraction between 0 and 1 of dp0 to use here. + if (dp00i >= q) then + ! This layer is much deeper than the fixed surface coordinates. + dp0ij(k) = dp00i + qhrlx(k+1) = qhybrlx + else + ! This layer spans the margins of the fixed surface coordinates. + ! In this case dp00i < q < dp0k. + dp0ij(k) = dilate * q + qhrlx(k+1) = qhybrlx * (dp0k(k) - dp00i) / & + ((dp0k(k) - q) + (q - dp00i)*qhybrlx) ! 1 at dp0k, qhybrlx at dp00i + endif + + ! The old equivalent code is: + ! hybrlx = 1.0 / qhybrlx + ! q = max( dp00i, dp0k(k) * (dp0k(k) / max(dp0k( k), p_int(k) - dp0cum(k)) ) ) + ! qts = 1.0 - (q-dp00i) / (dp0k(k) - dp00i) !0 at q = dp0k, 1 at q=dp00i + ! qhrlx( k+1) = 1.0 / (1.0 + qts*(hybrlx-1.0)) !1 at dp0k, qhybrlx at dp00i + endif + dp0cum(k+1) = dp0cum(k) + dp0ij(k) + p_int(k+1) = p_int(k) + h_col(k) + enddo !k + endif !qdep<1:else + + ! Identify the current fixed coordinate layers + fixlay = 1 !layer 1 always fixed + do k=2,nk + if (dp0cum(k) >= dilate * topiso_i_j) then + exit !layers k to nk might be isopycnal + endif + ! Top of layer is above topiso, i.e. always fixed coordinate layer + qhrlx(k+1) = 1.0 !no relaxation in fixed layers + fixlay = fixlay+1 + enddo !k + + fixall = fixlay + do k=fixall+1,nk + if (p_int(k+1) > dp0cum(k+1) + 0.1*dp0ij(k)) then + if ( (fixlay > fixall) .and. (p_int(k) > dp0cum(k)) ) then + ! --- The previous layer should remain fixed. + fixlay = fixlay-1 + endif + exit !layers k to nk might be isopycnal + endif + ! Sometimes fixed coordinate layer + qhrlx(k) = 1.0 !no relaxation in fixed layers + fixlay = fixlay+1 + enddo !k + +end subroutine hybgen_column_init + +!> The cushion function from Bleck & Benjamin, 1992, which returns a smoothly varying +!! but limited value that goes between dp0 and delp +real function cushn(delp, dp0) + real, intent(in) :: delp ! A thickness change [H ~> m or kg m-2] + real, intent(in) :: dp0 ! A non-negative reference thickness [H ~> m or kg m-2] + + real :: qq ! A limited ratio of delp/dp0 [nondim] + + ! These are the nondimensional parameters that define the cushion function. + real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn [nondim] +! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn [nondim] +! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn [nondim] + ! These are derivative nondimensional parameters. + ! real, parameter :: cusha = qqmn**2 * (qqmx-1.0) / (qqmx-qqmn)**2 + ! real, parameter :: I_qqmn = 1.0 / qqmn + real, parameter :: qq_scale = (qqmx-1.0) / (qqmx-qqmn)**2 + real, parameter :: I_qqmx = 1.0 / qqmx + + ! --- if delp >= qqmx*dp0 >> dp0, cushn returns delp. + ! --- if delp <= qqmn*dp0 << -dp0, cushn returns dp0. + + ! This is the original version from Hycom. + ! qq = max(qqmn, min(qqmx, delp/dp0)) + ! cushn = dp0 * (1.0 + cusha * (1.0-I_qqmn*qq)**2) * max(1.0, delp/(dp0*qqmx)) + + ! This is mathematically equivalent, has one fewer divide, and works as intended even if dp0 = 0. + if (delp >= qqmx*dp0) then + cushn = delp + elseif (delp < qqmn*dp0) then + cushn = max(dp0, delp * I_qqmx) + else + cushn = max(dp0, delp * I_qqmx) * (1.0 + qq_scale * ((delp / dp0) - qqmn)**2) + endif + +end function cushn + +!> Create a new grid for a column of water using the Hybgen algorithm. +subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & + fixlay, qhrlx, dp0ij, dp0cum, Rcv, h_in, dp_int) + type(hybgen_regrid_CS), intent(in) :: CS !< hybgen regridding control structure + integer, intent(in) :: nk !< number of layers + real, intent(in) :: thkbot !< thickness of bottom boundary layer [H ~> m or kg m-2] + real, intent(in) :: onem !< one m in pressure units [H ~> m or kg m-2] + real, intent(in) :: epsil !< small nonzero density to prevent division by zero [R ~> kg m-3] + real, intent(in) :: Rcv_tgt(nk) !< Target potential density [R ~> kg m-3] + integer, intent(in) :: fixlay !< deepest fixed coordinate layer + real, intent(in) :: qhrlx( nk+1) !< relaxation coefficient per timestep [nondim] + real, intent(in) :: dp0ij( nk) !< minimum layer thickness [H ~> m or kg m-2] + real, intent(in) :: dp0cum(nk+1) !< minimum interface depth [H ~> m or kg m-2] + real, intent(in) :: Rcv(nk) !< Coordinate potential density [R ~> kg m-3] + real, intent(in) :: h_in(nk) !< Layer thicknesses [H ~> m or kg m-2] + real, intent(out) :: dp_int(nk+1) !< The change in interface positions [H ~> m or kg m-2] + + ! --- ------------------------------------------------------ + ! --- hybrid grid generator, single column - regrid. + ! --- ------------------------------------------------------ + + ! Local variables + real :: p_new ! A new interface position [H ~> m or kg m-2] + real :: pres_in(nk+1) ! layer interface positions [H ~> m or kg m-2] + real :: p_int(nk+1) ! layer interface positions [H ~> m or kg m-2] + real :: h_col(nk) ! Updated layer thicknesses [H ~> m or kg m-2] + real :: q_frac ! A fraction of a layer to entrain [nondim] + real :: h_min ! The minimum layer thickness [H ~> m or kg m-2] + real :: h_hat3 ! Thickness movement upward across the interface between layers k-2 and k-3 [H ~> m or kg m-2] + real :: h_hat2 ! Thickness movement upward across the interface between layers k-1 and k-2 [H ~> m or kg m-2] + real :: h_hat ! Thickness movement upward across the interface between layers k and k-1 [H ~> m or kg m-2] + real :: h_hat0 ! A first guess at thickness movement upward across the interface + ! between layers k and k-1 [H ~> m or kg m-2] + real :: dh_cor ! Thickness changes [H ~> m or kg m-2] + real :: h1_tgt ! A target thickness for the top layer [H ~> m or kg m-2] + real :: tenm ! ten m in pressure units [H ~> m or kg m-2] + real :: onemm ! one mm in pressure units [H ~> m or kg m-2] + logical :: trap_errors + integer :: k + character(len=256) :: mesg ! A string for output messages + + ! This line needs to be consistent with the parameters set in cushn(). + real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn +! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn +! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn + + !### These hard-coded parameters should be changed to run-time variables. + tenm = 10.0*onem + onemm = 0.001*onem + + trap_errors = .true. + + do K=1,nk+1 ; dp_int(K) = 0.0 ; enddo + + p_int(1) = 0.0 + do k=1,nk + h_col(k) = max(h_in(k), 0.0) + p_int(K+1) = p_int(K) + h_col(k) + enddo + h_min = min( CS%min_thickness, p_int(nk+1)/real(CS%nk) ) + + if (trap_errors) then + do K=1,nk+1 ; pres_in(K) = p_int(K) ; enddo + endif + + ! Try to restore isopycnic conditions by moving layer interfaces + ! qhrlx(k) are relaxation amounts per timestep. + + ! Maintain prescribed thickness in layer k <= fixlay + ! There may be massless layers at the bottom, so work upwards. + do k=min(nk-1,fixlay),1,-1 + p_new = min(dp0cum(k+1), p_int(nk+1) - (nk-k)*h_min) ! This could be positive or negative. + dh_cor = p_new - p_int(K+1) + if (k= h_min) exit ! usually get here quickly + dh_cor = h_min - h_col(k) ! This is positive. + h_col(k) = h_min ! = h_col(k) + dh_cor + h_col(k+1) = h_col(k+1) - dh_cor + dp_int(k+1) = dp_int(k+1) + dh_cor + p_int(k+1) = p_int(fixlay+1) + enddo + if (h_col(nk) < h_min) then ! This should be uncommon, and should only arise at the level of roundoff. + do k=nk,2,-1 + if (h_col(k) >= h_min) exit + dh_cor = h_col(k) - h_min ! dh_cor is negative. + h_col(k-1) = h_col(k-1) + dh_cor + h_col(k) = h_min ! = h_col(k) - dh_cor + dp_int(k) = dp_int(k) + dh_cor + p_int(k) = p_int(k) + dh_cor + enddo + endif + + ! Remap the non-fixed layers. + + ! In the Hycom version, this loop was fused the loop correcting water that is + ! too light, and it ran down the water column, but if there are a set of layers + ! that are very dense, that structure can lead to all of the water being remapped + ! into a single thick layer. Splitting the loops and running the loop upwards + ! (as is done here avoids that catastrophic problem for layers that are far from + ! their targets. However, this code is still prone to a thin-thick-thin null mode. + do k=nk,fixlay+2,-1 + ! This is how the Hycom code would do this loop: do k=fixlay+1,nk ; if (k>fixlay+1) then + + if ((Rcv(k) > Rcv_tgt(k) + epsil)) then + ! Water in layer k is too dense, so try to dilute with water from layer k-1 + ! Do not move interface if k = fixlay + 1 + + if ((Rcv(k-1) >= Rcv_tgt(k-1)) .or. & + (p_int(k) <= dp0cum(k) + onem) .or. & + (h_col(k) <= h_col(k-1))) then + ! If layer k-1 is too light, there is a conflict in the direction the + ! inteface between them should move, so thicken the thinner of the two. + + if ((Rcv_tgt(k) - Rcv(k-1)) <= epsil) then + ! layer k-1 is far too dense, take the entire layer + ! If this code is working downward and this branch is repeated in a series + ! of successive layers, it can accumulate into a very thick homogenous layers. + h_hat0 = 0.0 ! This line was not in the Hycom version of hybgen.F90. + h_hat = dp0ij(k-1) - h_col(k-1) + else + ! Entrain enough from the layer above to bring layer k to its target density. + q_frac = (Rcv_tgt(k) - Rcv(k)) / (Rcv_tgt(k) - Rcv(k-1)) ! -1 <= q_frac < 0 + h_hat0 = q_frac*h_col(k) ! -h_col(k-1) <= h_hat0 < 0 + if (k == fixlay+2) then + ! Treat layer k-1 as fixed. + h_hat = max(h_hat0, dp0ij(k-1) - h_col(k-1)) + else + ! Maintain the minimum thickess of layer k-1. + h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1) + endif !fixlay+2:else + endif + ! h_hat is usually negative, so this check may be unnecessary if the values of + ! dp0ij are limited to not be below the seafloor? + h_hat = min(h_hat, p_int(nk+1) - p_int(k)) + + ! If isopycnic conditions cannot be achieved because of a blocking + ! layer (thinner than its minimum thickness) in the interior ocean, + ! move interface k-1 (and k-2 if necessary) upward + ! Only work on layers that are sufficiently far from the fixed near-surface layers. + if ((h_hat >= 0.0) .and. (k > fixlay+2) .and. (p_int(k-1) > dp0cum(k-1) + tenm)) then + + ! Only act if interface k-1 is near the bottom or layer k-2 could donate water. + if ( (p_int(nk+1) - p_int(k-1) < thkbot) .or. & + (h_col(k-2) > qqmx*dp0ij(k-2)) ) then + ! Determine how much water layer k-2 could supply without becoming too thin. + if (k == fixlay+3) then + ! Treat layer k-2 as fixed. + h_hat2 = max(h_hat0 - h_hat, dp0ij(k-2) - h_col(k-2)) + else + ! Maintain minimum thickess of layer k-2. + h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2) + endif !fixlay+3:else + + if (h_hat2 < -onemm) then + dh_cor = qhrlx(k-1) * max(h_hat2, -h_hat - h_col(k-1)) + h_col(k-2) = h_col(k-2) + dh_cor + h_col(k-1) = h_col(k-1) - dh_cor + dp_int(k-1) = dp_int(k-1) + dh_cor + p_int(k-1) = p_int(k-1) + dh_cor + ! Recalculate how much layer k-1 could donate to layer k. + h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1) + elseif (k <= fixlay+3) then + ! Do nothing. + elseif (p_int(k-2) > dp0cum(k-2) + tenm .and. & + (p_int(nk+1) - p_int(k-2) < thkbot .or. & + h_col(k-3) > qqmx*dp0ij(k-3))) then + + ! Determine how much water layer k-3 could supply without becoming too thin. + if (k == fixlay+4) then + ! Treat layer k-3 as fixed. + h_hat3 = max(h_hat0 - h_hat, dp0ij(k-3) - h_col(k-3)) + else + ! Maintain minimum thickess of layer k-3. + h_hat3 = cushn(h_col(k-3) + (h_hat0 - h_hat), dp0ij(k-3)) - h_col(k-3) + endif !fixlay+4:else + if (h_hat3 < -onemm) then + ! Water is moved from layer k-3 to k-2, but do not dilute layer k-2 too much. + dh_cor = qhrlx(k-2) * max(h_hat3, -h_col(k-2)) + h_col(k-3) = h_col(k-3) + dh_cor + h_col(k-2) = h_col(k-2) - dh_cor + dp_int(k-2) = dp_int(k-2) + dh_cor + p_int(k-2) = p_int(k-2) + dh_cor + + ! Now layer k-2 might be able donate to layer k-1. + h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2) + if (h_hat2 < -onemm) then + dh_cor = qhrlx(k-1) * (max(h_hat2, -h_hat - h_col(k-1)) ) + h_col(k-2) = h_col(k-2) + dh_cor + h_col(k-1) = h_col(k-1) - dh_cor + dp_int(k-1) = dp_int(k-1) + dh_cor + p_int(k-1) = p_int(k-1) + dh_cor + ! Recalculate how much layer k-1 could donate to layer k. + h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1) + endif !h_hat2 + endif !h_hat3 + endif !h_hat2:blocking + endif ! Layer k-2 could move. + endif ! blocking, i.e., h_hat >= 0, and far enough from the fixed layers to permit a change. + + if (h_hat < 0.0) then + ! entrain layer k-1 water into layer k, move interface up. + dh_cor = qhrlx(k) * h_hat + h_col(k-1) = h_col(k-1) + dh_cor + h_col(k) = h_col(k) - dh_cor + dp_int(k) = dp_int(k) + dh_cor + p_int(k) = p_int(k) + dh_cor + endif !entrain + + endif !too-dense adjustment + endif + + ! In the original Hycom version, there is not a break between these two loops. + enddo + + do k=fixlay+1,nk + if (Rcv(k) < Rcv_tgt(k) - epsil) then ! layer too light + ! Water in layer k is too light, so try to dilute with water from layer k+1. + ! Entrainment is not possible if layer k touches the bottom. + if (p_int(k+1) < p_int(nk+1)) then ! k dp0ij(k) + dp0ij(k+1)) then + h_hat = h_col(k+1) - cushn(h_col(k+1) - h_hat, dp0ij(k+1)) + endif + ! Try to bring layer layer k up to its minimum thickness. + h_hat = max(h_hat, dp0ij(k) - h_col(k)) + ! Do not drive layer k+1 below its minimum thickness or take more than half of it. + h_hat = min(h_hat, max(0.5*h_col(k+1), h_col(k+1) - dp0ij(k+1)) ) + else + ! Layers that touch the bottom can lose their entire contents. + h_hat = min(h_col(k+1), h_hat) + endif !p.k+2 0.0) then + ! Entrain layer k+1 water into layer k. + dh_cor = qhrlx(k+1) * h_hat + h_col(k) = h_col(k) + dh_cor + h_col(k+1) = h_col(k+1) - dh_cor + dp_int(k+1) = dp_int(k+1) + dh_cor + p_int(k+1) = p_int(k+1) + dh_cor + endif !entrain + + endif !too-light adjustment + endif !above bottom + endif !too light + + ! If layer above is still too thin, move interface down. + dh_cor = min(qhrlx(k-1) * min(dp0ij(k-1) - h_col(k-1), p_int(nk+1) - p_int(k)), h_col(k)) + if (dh_cor > 0.0) then + h_col(k-1) = h_col(k-1) + dh_cor + h_col(k) = h_col(k) - dh_cor + dp_int(k) = dp_int(k) + dh_cor + p_int(k) = p_int(k) + dh_cor + endif + + enddo !k Hybrid vertical coordinate relocation moving interface downward + + if (trap_errors) then + ! Verify that everything is consistent. + do k=1,nk + if (abs((h_col(k) - h_in(k)) + (dp_int(K) - dp_int(K+1))) > 1.0e-13*max(p_int(nk+1), onem)) then + write(mesg, '("k ",i4," h ",es13.4," h_in ",es13.4, " dp ",2es13.4," err ",es13.4)') & + k, h_col(k), h_in(k), dp_int(K), dp_int(K+1), (h_col(k) - h_in(k)) + (dp_int(K) - dp_int(K+1)) + call MOM_error(FATAL, "Mismatched thickness changes in hybgen_regrid: "//trim(mesg)) + endif + if (h_col(k) < 0.0) then ! Could instead do: -1.0e-15*max(p_int(nk+1), onem)) then + write(mesg, '("k ",i4," h ",es13.4," h_in ",es13.4, " dp ",2es13.4, " fixlay ",i4)') & + k, h_col(k), h_in(k), dp_int(K), dp_int(K+1), fixlay + call MOM_error(FATAL, "Significantly negative final thickness in hybgen_regrid: "//trim(mesg)) + endif + enddo + do K=1,nk+1 + if (abs(dp_int(K) - (p_int(K) - pres_in(K))) > 1.0e-13*max(p_int(nk+1), onem)) then + call MOM_error(FATAL, "Mismatched interface height changes in hybgen_regrid.") + endif + enddo + endif + +end subroutine hybgen_column_regrid + +end module MOM_hybgen_regrid + +! This code was translated in 2022 from the HYCOM hybgen code, which was primarily developed +! between 2000 and 2015, with some minor subsequent changes and bug fixes. diff --git a/src/ALE/MOM_hybgen_unmix.F90 b/src/ALE/MOM_hybgen_unmix.F90 new file mode 100644 index 0000000000..6e204f856b --- /dev/null +++ b/src/ALE/MOM_hybgen_unmix.F90 @@ -0,0 +1,502 @@ +!> This module contains the hybgen unmixing routines from HYCOM, with +!! modifications to follow the MOM6 coding conventions and several bugs fixed +module MOM_hybgen_unmix + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_EOS, only : EOS_type, calculate_density, calculate_density_derivs +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, param_file_type, log_param +use MOM_hybgen_regrid, only : hybgen_column_init +use MOM_hybgen_regrid, only : hybgen_regrid_CS, get_hybgen_regrid_params +use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : ocean_grid_type, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +!> Control structure containing required parameters for the hybgen coordinate generator +type, public :: hybgen_unmix_CS ; private + + integer :: nsigma !< Number of sigma levels used by HYBGEN + real :: hybiso !< Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3] + + real :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2] + real :: qhybrlx !< Hybgen relaxation amount per thermodynamic time steps [nondim] + + real, allocatable, dimension(:) :: & + dp0k, & !< minimum deep z-layer separation [H ~> m or kg m-2] + ds0k !< minimum shallow z-layer separation [H ~> m or kg m-2] + + real :: dpns !< depth to start terrain following [H ~> m or kg m-2] + real :: dsns !< depth to stop terrain following [H ~> m or kg m-2] + real :: min_dilate !< The minimum amount of dilation that is permitted when converting target + !! coordinates from z to z* [nondim]. This limit applies when wetting occurs. + real :: max_dilate !< The maximum amount of dilation that is permitted when converting target + !! coordinates from z to z* [nondim]. This limit applies when drying occurs. + + real :: topiso_const !< Shallowest depth for isopycnal layers [H ~> m or kg m-2] + ! real, dimension(:,:), allocatable :: topiso + + real :: ref_pressure !< Reference pressure for density calculations [R L2 T-2 ~> Pa] + real, allocatable, dimension(:) :: target_density !< Nominal density of interfaces [R ~> kg m-3] + +end type hybgen_unmix_CS + +public hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix +public set_hybgen_unmix_params + +contains + +!> Initialise a hybgen_unmix_CS control structure and store its parameters +subroutine init_hybgen_unmix(CS, GV, US, param_file, hybgen_regridCS) + type(hybgen_unmix_CS), pointer :: CS !< Unassociated pointer to hold the control structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file + type(hybgen_regrid_CS), pointer :: hybgen_regridCS !< Control structure for hybgen + !! regridding for sharing parameters. + + character(len=40) :: mdl = "MOM_hybgen" ! This module's name. + integer :: k + + if (associated(CS)) call MOM_error(FATAL, "init_hybgen_unmix: CS already associated!") + allocate(CS) + allocate(CS%target_density(GV%ke)) + + allocate(CS%dp0k(GV%ke), source=0.0) ! minimum deep z-layer separation + allocate(CS%ds0k(GV%ke), source=0.0) ! minimum shallow z-layer separation + + ! Set the parameters for the hybgen unmixing from a hybgen regridding control structure. + call get_hybgen_regrid_params(hybgen_regridCS, ref_pressure=CS%ref_pressure, & + nsigma=CS%nsigma, dp0k=CS%dp0k, ds0k=CS%ds0k, & + dp00i=CS%dp00i, topiso_const=CS%topiso_const, qhybrlx=CS%qhybrlx, & + hybiso=CS%hybiso, min_dilate=CS%min_dilate, max_dilate=CS%max_dilate, & + target_density=CS%target_density) + + ! Determine the depth range over which to use a sigma (terrain-following) coordinate. + ! --- terrain following starts at depth dpns and ends at depth dsns + if (CS%nsigma == 0) then + CS%dpns = CS%dp0k(1) + CS%dsns = 0.0 + else + CS%dpns = 0.0 + CS%dsns = 0.0 + do k=1,CS%nsigma + CS%dpns = CS%dpns + CS%dp0k(k) + CS%dsns = CS%dsns + CS%ds0k(k) + enddo !k + endif !nsigma + +end subroutine init_hybgen_unmix + +!> This subroutine deallocates memory in the control structure for the hybgen unmixing module +subroutine end_hybgen_unmix(CS) + type(hybgen_unmix_CS), pointer :: CS !< Coordinate control structure + + ! nothing to do + if (.not. associated(CS)) return + + deallocate(CS%target_density) + deallocate(CS%dp0k, CS%ds0k) + deallocate(CS) +end subroutine end_hybgen_unmix + +!> This subroutine can be used to set the parameters for the hybgen module +subroutine set_hybgen_unmix_params(CS, min_thickness) + type(hybgen_unmix_CS), pointer :: CS !< Coordinate unmixing control structure + real, optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2] + + if (.not. associated(CS)) call MOM_error(FATAL, "set_hybgen_params: CS not associated") + +! if (present(min_thickness)) CS%min_thickness = min_thickness +end subroutine set_hybgen_unmix_params + + +!> Unmix the properties in the lowest layer with mass if it is too light, and make +!! any other changes to the water column to prepare for regridding. +subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(hybgen_unmix_CS), intent(in) :: CS !< hybgen control structure + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + type(tracer_registry_type), pointer :: Reg !< Tracer registry structure + integer, intent(in) :: ntr !< The number of tracers in the registry, or + !! 0 if the registry is not in use. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + +! --- -------------------------------------------- +! --- hybrid grid generator, single j-row (part A). +! --- -------------------------------------------- + + character(len=256) :: mesg ! A string for output messages + integer :: fixlay ! deepest fixed coordinate layer + real :: qhrlx( GV%ke+1) ! relaxation coefficient per timestep [nondim] + real :: dp0ij( GV%ke) ! minimum layer thickness [H ~> m or kg m-2] + real :: dp0cum(GV%ke+1) ! minimum interface depth [H ~> m or kg m-2] + + real :: Rcv_tgt(GV%ke) ! Target potential density [R ~> kg m-3] + real :: temp(GV%ke) ! A column of potential temperature [degC] + real :: saln(GV%ke) ! A column of salinity [ppt] + real :: Rcv(GV%ke) ! A column of coordinate potential density [R ~> kg m-3] + real :: h_col(GV%ke) ! A column of layer thicknesses [H ~> m or kg m-2] + real :: p_col(GV%ke) ! A column of reference pressures [R L2 T-2 ~> Pa] + real :: tracer(GV%ke,max(ntr,1)) ! Columns of each tracer [Conc] + real :: h_tot ! Total thickness of the water column [H ~> m or kg m-2] + real :: nominalDepth ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2] + real :: h_thin ! A negligibly small thickness to identify essentially + ! vanished layers [H ~> m or kg m-2] + real :: dilate ! A factor by which to dilate the target positions from z to z* [nondim] + + real :: Th_tot_in, Th_tot_out ! Column integrated temperature [degC H ~> degC m or degC kg m-2] + real :: Sh_tot_in, Sh_tot_out ! Column integrated salinity [ppt H ~> ppt m or ppt kg m-2] + real :: Trh_tot_in(max(ntr,1)) ! Initial column integrated tracer amounts [conc H ~> conc m or conc kg m-2] + real :: Trh_tot_out(max(ntr,1)) ! Final column integrated tracer amounts [conc H ~> conc m or conc kg m-2] + + logical :: debug_conservation ! If true, test for non-conservation. + logical :: terrain_following ! True if this column is terrain following. + integer :: trcflg(max(ntr,1)) ! Hycom tracer type flag for each tracer + integer :: i, j, k, nk, m + + nk = GV%ke + + ! Set all tracers to be passive. Setting this to 2 treats a tracer like temperature. + trcflg(:) = 3 + + h_thin = 1e-6*GV%m_to_H + debug_conservation = .false. ! Set this to true for debugging + + p_col(:) = CS%ref_pressure + + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then + + h_tot = 0.0 + do k=1,nk + ! Rcv_tgt(k) = theta(i,j,k) ! If a 3-d target density were set up in theta, use that here. + Rcv_tgt(k) = CS%target_density(k) ! MOM6 does not yet support 3-d target densities. + h_col(k) = h(i,j,k) + h_tot = h_tot + h_col(k) + temp(k) = tv%T(i,j,k) + saln(k) = tv%S(i,j,k) + enddo + + ! This sets the potential density from T and S. + call calculate_density(temp, saln, p_col, Rcv, tv%eqn_of_state) + + do m=1,ntr ; do k=1,nk + tracer(k,m) = Reg%Tr(m)%t(i,j,k) + enddo ; enddo + + ! Store original amounts to test for conservation of temperature, salinity, and tracers. + if (debug_conservation) then + Th_tot_in = 0.0 ; Sh_tot_in = 0.0 ; Trh_tot_in(:) = 0.0 + do k=1,nk + Sh_tot_in = Sh_tot_in + h_col(k)*saln(k) + Th_tot_in = Th_tot_in + h_col(k)*temp(k) + enddo + do m=1,ntr ; do k=1,nk + Trh_tot_in(m) = Trh_tot_in(m) + h_col(k)*tracer(k,m) + enddo ; enddo + endif + + ! The following block of code is used to trigger z* stretching of the targets heights. + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H + if (h_tot <= CS%min_dilate*nominalDepth) then + dilate = CS%min_dilate + elseif (h_tot >= CS%max_dilate*nominalDepth) then + dilate = CS%max_dilate + else + dilate = h_tot / nominalDepth + endif + + terrain_following = (h_tot < dilate*CS%dpns) .and. (CS%dpns >= CS%dsns) + + ! Convert the regridding parameters into specific constraints for this column. + call hybgen_column_init(nk, CS%nsigma, CS%dp0k, CS%ds0k, CS%dp00i, & + CS%topiso_const, CS%qhybrlx, CS%dpns, CS%dsns, h_tot, dilate, & + h_col, fixlay, qhrlx, dp0ij, dp0cum) + + ! Do any unmixing of the column that is needed to move the layer properties toward their targets. + call hybgen_column_unmix(CS, nk, Rcv_tgt, temp, saln, Rcv, tv%eqn_of_state, & + ntr, tracer, trcflg, fixlay, qhrlx, h_col, & + terrain_following, h_thin) + + ! Store the output from hybgen_unmix in the 3-d arrays. + do k=1,nk + h(i,j,k) = h_col(k) + enddo + ! Note that temperature and salinity are among the tracers unmixed here. + do m=1,ntr ; do k=1,nk + Reg%Tr(m)%t(i,j,k) = tracer(k,m) + enddo ; enddo + ! However, temperature and salinity may have been treated differently from other tracers. + do k=1,nk + tv%T(i,j,k) = temp(k) + tv%S(i,j,k) = saln(k) + enddo + + ! Test for conservation of temperature, salinity, and tracers. + if (debug_conservation) then + Th_tot_out = 0.0 ; Sh_tot_out = 0.0 ; Trh_tot_out(:) = 0.0 + do k=1,nk + Sh_tot_out = Sh_tot_out + h_col(k)*saln(k) + Th_tot_out = Th_tot_out + h_col(k)*temp(k) + enddo + do m=1,ntr ; do k=1,nk + Trh_tot_out(m) = Trh_tot_out(m) + h_col(k)*tracer(k,m) + enddo ; enddo + if (abs(Sh_tot_in - Sh_tot_out) > 1.e-15*(abs(Sh_tot_in) + abs(Sh_tot_out))) then + write(mesg, '("i,j=",2i8,"Sh_tot = ",2es17.8," err = ",es13.4)') & + i, j, Sh_tot_in, Sh_tot_out, (Sh_tot_in - Sh_tot_out) + call MOM_error(FATAL, "Mismatched column salinity in hybgen_unmix: "//trim(mesg)) + endif + if (abs(Th_tot_in - Th_tot_out) > 1.e-10*(abs(Th_tot_in) + abs(Th_tot_out))) then + write(mesg, '("i,j=",2i8,"Th_tot = ",2es17.8," err = ",es13.4)') & + i, j, Th_tot_in, Th_tot_out, (Th_tot_in - Th_tot_out) + call MOM_error(FATAL, "Mismatched column temperature in hybgen_unmix: "//trim(mesg)) + endif + do m=1,ntr + if (abs(Trh_tot_in(m) - Trh_tot_out(m)) > 1.e-10*(abs(Trh_tot_in(m)) + abs(Trh_tot_out(m)))) then + write(mesg, '("i,j=",2i8,"Trh_tot(",i2,") = ",2es17.8," err = ",es13.4)') & + i, j, m, Trh_tot_in(m), Trh_tot_out(m), (Trh_tot_in(m) - Trh_tot_out(m)) + call MOM_error(FATAL, "Mismatched column tracer in hybgen_unmix: "//trim(mesg)) + endif + enddo + endif + endif ; enddo ; enddo !i & j. + +end subroutine hybgen_unmix + + +!> Unmix the properties in the lowest layer if it is too light. +subroutine hybgen_column_unmix(CS, nk, Rcv_tgt, temp, saln, Rcv, eqn_of_state, & + ntr, tracer, trcflg, fixlay, qhrlx, h_col, & + terrain_following, h_thin) + type(hybgen_unmix_CS), intent(in) :: CS !< hybgen unmixing control structure + integer, intent(in) :: nk !< The number of layers + integer, intent(in) :: fixlay !< deepest fixed coordinate layer + real, intent(in) :: qhrlx(nk+1) !< Relaxation fraction per timestep [nondim], < 1. + real, intent(in) :: Rcv_tgt(nk) !< Target potential density [R ~> kg m-3] + real, intent(inout) :: temp(nk) !< A column of potential temperature [degC] + real, intent(inout) :: saln(nk) !< A column of salinity [ppt] + real, intent(inout) :: Rcv(nk) !< Coordinate potential density [R ~> kg m-3] + type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure + integer, intent(in) :: ntr !< The number of registered passive tracers + real, intent(inout) :: tracer(nk, max(ntr,1)) !< Columns of the passive tracers [Conc] + integer, intent(in) :: trcflg(max(ntr,1)) !< Hycom tracer type flag for each tracer + real, intent(inout) :: h_col(nk+1) !< Layer thicknesses [H ~> m or kg m-2] + logical, intent(in) :: terrain_following !< True if this column is terrain following + real, intent(in) :: h_thin !< A negligibly small thickness to identify + !! essentially vanished layers [H ~> m or kg m-2] + +! +! --- ------------------------------------------------------------------ +! --- hybrid grid generator, single column - ummix lowest massive layer. +! --- ------------------------------------------------------------------ +! + ! Local variables + real :: h_hat ! A portion of a layer to move across an interface [H ~> m or kg m-2] + real :: delt, deltm ! Temperature differences between successive layers [degC] + real :: dels, delsm ! Salinity differences between successive layers [ppt] + real :: abs_dRdT ! The absolute value of the derivative of the coordinate density + ! with temperature [R degC-1 ~> kg m-3 degC-1] + real :: abs_dRdS ! The absolute value of the derivative of the coordinate density + ! with salinity [R ppt-1 ~> kg m-3 ppt-1] + real :: q, qts ! Nondimensional fractions in the range of 0 to 1 [nondim] + real :: frac_dts ! The fraction of the temperature or salinity difference between successive + ! layers by which the source layer's property changes by the loss of water + ! that matches the destination layers properties via unmixing [nondim]. + real :: qtr ! The fraction of the water that will come from the layer below, + ! used for updating the concentration of passive tracers [nondim] + real :: swap_T ! A swap variable for temperature [degC] + real :: swap_S ! A swap variable for salinity [ppt] + real :: swap_R ! A swap variable for the coordinate potential density [R ~> kg m-3] + real :: swap_tr ! A temporary swap variable for the tracers [conc] + logical, parameter :: lunmix=.true. ! unmix a too light deepest layer + integer :: k, ka, kp, kt, m + + ! --- identify the deepest layer kp with significant thickness (> h_thin) + kp = 2 !minimum allowed value + do k=nk,3,-1 + if (h_col(k) >= h_thin) then + kp = k + exit + endif + enddo !k + + k = kp !at least 2 + ka = max(k-2,1) !k might be 2 +! + if ( ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth + (h_col(k-1) >= h_thin) .and. & ! layer above not too thin + (Rcv_tgt(k) > Rcv(k)) .and. & ! layer is lighter than its target + ((Rcv(k-1) > Rcv(k)) .and. (Rcv(ka) > Rcv(k))) ) then +! +! --- water in the deepest inflated layer with significant thickness +! --- (kp) is too light, and it is lighter than the two layers above. +! --- +! --- this should only occur when relaxing or nudging layer thickness +! --- and is a bug (bad interaction with tsadvc) even in those cases +! --- +! --- entrain the entire layer into the one above +!--- note the double negative in T=T-q*(T-T'), equiv. to T=T+q*(T'-T) + q = h_col(k) / (h_col(k) + h_col(k-1)) + temp(k-1) = temp(k-1) - q*(temp(k-1) - temp(k)) + saln(k-1) = saln(k-1) - q*(saln(k-1) - saln(k)) + call calculate_density(temp(k-1), saln(k-1), CS%ref_pressure, Rcv(k-1), eqn_of_state) + + do m=1,ntr + tracer(k-1,m) = tracer(k-1,m) - q*(tracer(k-1,m) - tracer(k,m) ) + enddo !m +! --- entrained the entire layer into the one above, so now kp=kp-1 + h_col(k-1) = h_col(k-1) + h_col(k) + h_col(k) = 0.0 + kp = k-1 + elseif ( ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth + (h_col(k-1) >= h_thin) .and. & ! layer above not too thin + (Rcv_tgt(k) > Rcv(k)) .and. & ! layer is lighter than its target + (Rcv(k-1) > Rcv(k)) ) then +! --- water in the deepest inflated layer with significant thickness +! --- (kp) is too light, and it is lighter than the layer above, but not the layer two above. +! --- +! --- swap the entire layer with the one above. + if (h_col(k) <= h_col(k-1)) then + ! The bottom layer is thinner; swap the entire bottom layer with a portion of the layer above. + q = h_col(k) / h_col(k-1) !<=1.0 + + swap_T = temp(k-1) + temp(k-1) = temp(k-1) + q*(temp(k) - temp(k-1)) + temp(k) = swap_T + + swap_S = saln(k-1) + saln(k-1) = saln(k-1) + q*(saln(k) - saln(k-1)) + saln(k) = swap_S + + Rcv(k) = Rcv(k-1) + call calculate_density(temp(k-1), saln(k-1), CS%ref_pressure, Rcv(k-1), eqn_of_state) + + do m=1,ntr + swap_tr = tracer(k-1,m) + tracer(k-1,m) = tracer(k-1,m) - q * (tracer(k-1,m) - tracer(k,m)) + tracer(k,m) = swap_tr + enddo !m + else + ! The bottom layer is thicker; swap the entire layer above with a portion of the bottom layer. + q = h_col(k-1) / h_col(k) !<1.0 + + swap_T = temp(k) + temp(k) = temp(k) + q*(temp(k-1) - temp(k)) + temp(k-1) = swap_T + + swap_S = saln(k) + saln(k) = saln(k) + q*(saln(k-1) - saln(k)) + saln(k-1) = swap_S + + Rcv(k-1) = Rcv(k) + call calculate_density(temp(k), saln(k), CS%ref_pressure, Rcv(k), eqn_of_state) + + do m=1,ntr + swap_tr = tracer(k,m) + tracer(k,m) = tracer(k,m) + q * (tracer(k-1,m) - tracer(k,m)) + tracer(k-1,m) = swap_tr + enddo !m + endif !bottom too light + endif + + k = kp !at least 2 + ka = max(k-2,1) !k might be 2 + + if ( lunmix .and. & ! usually .true. + ((k > fixlay+1) .and. (.not.terrain_following)) .and. & ! layer not fixed depth + (h_col(k-1) >= h_thin) .and. & ! layer above not too thin + (Rcv(k) < Rcv_tgt(k)) .and. & ! layer is lighter than its target + (Rcv(k) > Rcv_tgt(k-1)) .and. & ! layer is denser than the target above + (abs(Rcv_tgt(k-1) - Rcv(k-1)) < CS%hybiso) .and. & ! layer above is near its target + (Rcv(k) - Rcv(k-1) > 0.001*(Rcv_tgt(k) - Rcv_tgt(k-1))) ) then +! +! --- water in the deepest inflated layer with significant thickness (kp) is too +! --- light but denser than the layer above, with the layer above near-isopycnal +! --- +! --- split layer into 2 sublayers, one near the desired density +! --- and one exactly matching the T&S properties of layer k-1. +! --- To prevent "runaway" T or S, the result satisfies either +! --- abs(T.k - T.k-1) <= abs(T.k-N - T.k-1) or +! --- abs(S.k - S.k-1) <= abs(S.k-N - S.k-1) where +! --- Rcv.k-1 - Rcv.k-N is at least Rcv_tgt(k-1) - Rcv_tgt(k-2) +! --- It is also limited to a 50% change in layer thickness. + + ka = 1 + do kt=k-2,2,-1 + if ( Rcv(k-1) - Rcv(kt) >= Rcv_tgt(k-1) - Rcv_tgt(k-2) ) then + ka = kt !usually k-2 + exit + endif + enddo + + delsm = abs(saln(ka) - saln(k-1)) + dels = abs(saln(k-1) - saln(k)) + deltm = abs(temp(ka) - temp(k-1)) + delt = abs(temp(k-1) - temp(k)) + + call calculate_density_derivs(temp(k-1), saln(k-1), CS%ref_pressure, abs_dRdT, abs_dRdS, eqn_of_state) + ! Bound deltm and delsm based on the equation of state and density differences between layers. + abs_dRdT = abs(abs_dRdT) ; abs_dRdS = abs(abs_dRdS) + if (abs_dRdT * deltm > Rcv_tgt(k)-Rcv_tgt(k-1)) deltm = (Rcv_tgt(k)-Rcv_tgt(k-1)) / abs_dRdT + if (abs_dRdS * delsm > Rcv_tgt(k)-Rcv_tgt(k-1)) delsm = (Rcv_tgt(k)-Rcv_tgt(k-1)) / abs_dRdS + + qts = 0.0 + if (qts*dels < min(delsm-dels, dels)) qts = min(delsm-dels, dels) / dels + if (qts*delt < min(deltm-delt, delt)) qts = min(deltm-delt, delt) / delt + + ! Note that Rcv_tgt(k) > Rcv(k) > Rcv(k-1), and 0 <= qts <= 1. + ! qhrlx is relaxation coefficient (inverse baroclinic time steps), 0 <= qhrlx <= 1. + ! This takes the minimum of the two estimates. + if ((1.0+qts) * (Rcv_tgt(k)-Rcv(k)) < qts * (Rcv_tgt(k)-Rcv(k-1))) then + q = qhrlx(k) * ((Rcv_tgt(k)-Rcv(k)) / (Rcv_tgt(k)-Rcv(k-1))) + else + q = qhrlx(k) * (qts / (1.0+qts)) ! upper sublayer <= 50% of total + endif + frac_dts = q / (1.0-q) ! 0 <= q <= 0.5, so 0 <= frac_dts <= 1 + + h_hat = q * h_col(k) + h_col(k-1) = h_col(k-1) + h_hat + h_col(k) = h_col(k) - h_hat + + temp(k) = temp(k) + frac_dts * (temp(k) - temp(k-1)) + saln(k) = saln(k) + frac_dts * (saln(k) - saln(k-1)) + call calculate_density(temp(k), saln(k), CS%ref_pressure, Rcv(k), eqn_of_state) + + if ((ntr > 0) .and. (h_hat /= 0.0)) then + ! qtr is the fraction of the new upper layer from the old lower layer. + ! The nonconservative original from Hycom: qtr = h_hat / max(h_hat, h_col(k)) !between 0 and 1 + qtr = h_hat / h_col(k-1) ! Between 0 and 1, noting the h_col(k-1) = h_col(k-1) + h_hat above. + do m=1,ntr + if (trcflg(m) == 2) then !temperature tracer + tracer(k,m) = tracer(k,m) + frac_dts * (tracer(k,m) - tracer(k-1,m)) + else !standard tracer - not split into two sub-layers + tracer(k-1,m) = tracer(k-1,m) + qtr * (tracer(k,m) - tracer(k-1,m)) + endif !trcflg + enddo !m + endif !tracers + endif !too light + +! ! Fill properties of massless or near-massless (thickness < h_thin) layers +! ! This was in the Hycom verion, but it appears to be unnecessary in MOM6. +! do k=kp+1,nk +! ! --- fill thin and massless layers on sea floor with fluid from above +! Rcv(k) = Rcv(k-1) +! do m=1,ntr +! tracer(k,m) = tracer(k-1,m) +! enddo !m +! saln(k) = saln(k-1) +! temp(k) = temp(k-1) +! enddo !k + +end subroutine hybgen_column_unmix + +end module MOM_hybgen_unmix diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 4c99ce457d..78159146a2 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -21,7 +21,7 @@ module MOM_regridding use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR -use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE +use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike @@ -31,6 +31,8 @@ module MOM_regridding use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom use coord_slight, only : init_coord_slight, slight_CS, set_slight_params, build_slight_column, end_coord_slight use coord_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt +use MOM_hybgen_regrid, only : hybgen_regrid, hybgen_regrid_CS, init_hybgen_regrid, end_hybgen_regrid +use MOM_hybgen_regrid, only : write_Hybgen_coord_file implicit none ; private @@ -119,17 +121,21 @@ module MOM_regridding !! If false, use more robust forms of the same remapping expressions. logical :: remap_answers_2018 = .true. + logical :: use_hybgen_unmix = .false. !< If true, use the hybgen unmixing code before remapping + type(zlike_CS), pointer :: zlike_CS => null() !< Control structure for z-like coordinate generator type(sigma_CS), pointer :: sigma_CS => null() !< Control structure for sigma coordinate generator type(rho_CS), pointer :: rho_CS => null() !< Control structure for rho coordinate generator type(hycom_CS), pointer :: hycom_CS => null() !< Control structure for hybrid coordinate generator type(slight_CS), pointer :: slight_CS => null() !< Control structure for Slight-coordinate generator type(adapt_CS), pointer :: adapt_CS => null() !< Control structure for adaptive coordinate generator + type(hybgen_regrid_CS), pointer :: hybgen_CS => NULL() !< Control structure for hybgen regridding end type ! The following routines are visible to the outside world public initialize_regridding, end_regridding, regridding_main +public regridding_preadjust_reqs, convective_adjustment public inflate_vanished_layers_old, check_remapping_grid, check_grid_column public set_regrid_params, get_regrid_size, write_regrid_file public uniformResolution, setCoordinateResolution @@ -148,6 +154,7 @@ module MOM_regridding " SIGMA - terrain following coordinates\n"//& " RHO - continuous isopycnal\n"//& " HYCOM1 - HyCOM-like hybrid coordinate\n"//& + " HYBGEN - Hybrid coordinate from the Hycom hybgen code\n"//& " SLIGHT - stretched coordinates above continuous isopycnal\n"//& " ADAPTIVE - optimize for smooth neutral density surfaces" @@ -458,6 +465,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! This is a work around to apparently needed to work with the from_Z initialization... ??? if (coordinateMode(coord_mode) == REGRIDDING_ZSTAR .or. & coordinateMode(coord_mode) == REGRIDDING_HYCOM1 .or. & + coordinateMode(coord_mode) == REGRIDDING_HYBGEN .or. & coordinateMode(coord_mode) == REGRIDDING_SLIGHT .or. & coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then ! Adjust target grid to be consistent with maximum_depth @@ -511,7 +519,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m endif ! initialise coordinate-specific control structure - call initCoord(CS, GV, US, coord_mode) + call initCoord(CS, GV, US, coord_mode, param_file) if (main_parameters .and. coord_is_state_dependent) then call get_param(param_file, mdl, "P_REF", P_Ref, & @@ -536,6 +544,13 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call set_regrid_params(CS, min_thickness=0.) endif + CS%use_hybgen_unmix = .false. + if (coordinateMode(coord_mode) == REGRIDDING_HYBGEN) then + call get_param(param_file, mdl, "USE_HYBGEN_UNMIX", CS%use_hybgen_unmix, & + "If true, use hybgen unmixing code before regridding.", & + default=.false.) + endif + if (coordinateMode(coord_mode) == REGRIDDING_SLIGHT) then ! Set SLight-specific regridding parameters. call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, & @@ -744,6 +759,7 @@ subroutine end_regridding(CS) if (associated(CS%hycom_CS)) call end_coord_hycom(CS%hycom_CS) if (associated(CS%slight_CS)) call end_coord_slight(CS%slight_CS) if (associated(CS%adapt_CS)) call end_coord_adapt(CS%adapt_CS) + if (associated(CS%hybgen_CS)) call end_hybgen_regrid(CS%hybgen_CS) deallocate( CS%coordinateResolution ) if (allocated(CS%target_density)) deallocate( CS%target_density ) @@ -754,8 +770,8 @@ end subroutine end_regridding !------------------------------------------------------------------------------ !> Dispatching regridding routine for orchestrating regridding & remapping -subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, & - conv_adjust, PCM_cell) +subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, conv_adjust, & + frac_shelf_h, PCM_cell) !------------------------------------------------------------------------------ ! This routine takes care of (1) building a new grid and (2) remapping between ! the old grid and the new grid. The creation of the new grid can be based @@ -778,22 +794,29 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after !! the last time step - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...) + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamical variables (T, S, ...) real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface + logical, intent(in ) :: conv_adjust !< If true, regridding_main should do + !! convective adjustment, but because it no + !! longer does convective adjustment this must + !! be false. This argument has been retained to + !! trap inconsistent code, but will eventually + !! be eliminated. real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage - logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true ! Local variables real :: trickGnuCompiler - logical :: do_convective_adjustment - do_convective_adjustment = .true. - if (present(conv_adjust)) do_convective_adjustment = conv_adjust + if (conv_adjust) call MOM_error(FATAL, & + "regridding_main: convective adjustment no longer is done inside of regridding_main. "//& + "The code needs to be modified to call regridding_main() with conv_adjust=.false, "//& + "and a call to convective_adjustment added before calling regridding_main() "//& + "if regridding_preadjust_reqs() indicates that this is necessary.") if (present(PCM_cell)) PCM_cell(:,:,:) = .false. select case ( CS%regridding_scheme ) @@ -808,7 +831,6 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ call build_sigma_grid( CS, G, GV, h, dzInterface ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_RHO ) - if (do_convective_adjustment) call convective_adjustment(G, GV, h, tv) call build_rho_grid( G, GV, G%US, h, tv, dzInterface, remapCS, CS, frac_shelf_h ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ARBITRARY ) @@ -816,6 +838,9 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_HYCOM1 ) call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) + case ( REGRIDDING_HYBGEN ) + call hybgen_regrid(G, GV, G%US, h, tv, CS%hybgen_CS, dzInterface, PCM_cell) + call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_SLIGHT ) call build_grid_SLight( G, GV, G%US, h, tv, dzInterface, CS ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) @@ -835,6 +860,39 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ end subroutine regridding_main +!------------------------------------------------------------------------------ +!> This routine returns flags indicating which pre-remapping state adjustments +!! are needed depending on the coordinate mode in use. +subroutine regridding_preadjust_reqs(CS, do_conv_adj, do_hybgen_unmix, hybgen_CS) + + ! Arguments + type(regridding_CS), intent(in) :: CS !< Regridding control structure + logical, intent(out) :: do_conv_adj !< Convective adjustment should be done + logical, intent(out) :: do_hybgen_unmix !< Hybgen unmixing should be done + type(hybgen_regrid_CS), pointer, & + optional, intent(out) :: hybgen_CS !< Control structure for hybgen regridding for sharing parameters. + + + do_conv_adj = .false. ; do_hybgen_unmix = .false. + select case ( CS%regridding_scheme ) + + case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_SIGMA, REGRIDDING_ARBITRARY, & + REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) + do_conv_adj = .false. ; do_hybgen_unmix = .false. + case ( REGRIDDING_RHO ) + do_conv_adj = .true. ; do_hybgen_unmix = .false. + case ( REGRIDDING_HYBGEN ) + do_conv_adj = .false. ; do_hybgen_unmix = CS%use_hybgen_unmix + case default + call MOM_error(FATAL,'MOM_regridding, regridding_preadjust_reqs: '//& + 'Unknown regridding scheme selected!') + end select ! type of grid + + if (present(hybgen_CS) .and. do_hybgen_unmix) hybgen_CS => CS%hybgen_CS + +end subroutine regridding_preadjust_reqs + + !> Calculates h_new from h + delta_k dzInterface subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) type(regridding_CS), intent(in) :: CS !< Regridding control structure @@ -1869,14 +1927,17 @@ subroutine convective_adjustment(G, GV, h, tv) !------------------------------------------------------------------------------ ! Local variables - integer :: i, j, k - real :: T0, T1 ! temperatures - real :: S0, S1 ! salinities - real :: r0, r1 ! densities - real :: h0, h1 + real :: T0, T1 ! temperatures of two layers [degC] + real :: S0, S1 ! salinities of two layers [ppt] + real :: r0, r1 ! densities of two layers [R ~> kg m-3] + real :: h0, h1 ! Layer thicknesses [H ~> m or kg m-2] + real, dimension(GV%ke) :: p_col ! A column of zero pressures [R L2 T-2 ~> Pa] + real, dimension(GV%ke) :: densities ! Densities in the column [R ~> kg m-3] logical :: stratified - real, dimension(GV%ke) :: p_col, densities + integer :: i, j, k + !### Doing convective adjustment based on potential densities with zero pressure seems + ! questionable, although it does avoid ambiguous sorting. -RWH p_col(:) = 0. ! Loop on columns @@ -1905,6 +1966,8 @@ subroutine convective_adjustment(G, GV, h, tv) call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state) call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), & densities(k+1), tv%eqn_of_state ) + ! Because p_col is has uniform values, these calculate_density calls are equivalent to + ! densities(k) = r1 ; densities(k+1) = r0 stratified = .false. endif enddo ! k @@ -1940,8 +2003,8 @@ function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) scheme = coordinateMode(coordMode) select case ( scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_SIGMA_SHELF_ZSTAR, & - REGRIDDING_ADAPTIVE ) + case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, & + REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_ADAPTIVE ) uniformResolution(:) = maxDepth / real(nk) case ( REGRIDDING_RHO ) @@ -1960,13 +2023,14 @@ end function uniformResolution !> Initialize the coordinate resolutions by calling the appropriate initialization !! routine for the specified coordinate mode. -subroutine initCoord(CS, GV, US, coord_mode) +subroutine initCoord(CS, GV, US, coord_mode, param_file) type(regridding_CS), intent(inout) :: CS !< Regridding control structure character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode. !! See the documentation for regrid_consts !! for the recognized values. type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file select case (coordinateMode(coord_mode)) case (REGRIDDING_ZSTAR) @@ -1980,6 +2044,8 @@ subroutine initCoord(CS, GV, US, coord_mode) case (REGRIDDING_HYCOM1) call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, & CS%interp_CS) + case (REGRIDDING_HYBGEN) + call init_hybgen_regrid(CS%hybgen_CS, GV, US, param_file) case (REGRIDDING_SLIGHT) call init_coord_slight(CS%slight_CS, CS%nk, CS%ref_pressure, CS%target_density, & CS%interp_CS, GV%m_to_H) @@ -2121,6 +2187,11 @@ subroutine write_regrid_file( CS, GV, filepath ) type(file_type) :: IO_handle ! The I/O handle of the fileset real :: ds(GV%ke), dsi(GV%ke+1) + if (CS%regridding_scheme == REGRIDDING_HYBGEN) then + call write_Hybgen_coord_file(GV, CS%hybgen_CS, filepath) + return + endif + 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) ) @@ -2209,7 +2280,8 @@ function getCoordinateUnits( CS ) character(len=20) :: getCoordinateUnits select case ( CS%regridding_scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) + case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_SLIGHT, & + REGRIDDING_ADAPTIVE ) getCoordinateUnits = 'meter' case ( REGRIDDING_SIGMA_SHELF_ZSTAR ) getCoordinateUnits = 'meter/fraction' @@ -2247,6 +2319,8 @@ function getCoordinateShortName( CS ) getCoordinateShortName = 'coordinate' case ( REGRIDDING_HYCOM1 ) getCoordinateShortName = 'z-rho' + case ( REGRIDDING_HYBGEN ) + getCoordinateShortName = 'hybrid' case ( REGRIDDING_SLIGHT ) getCoordinateShortName = 's-rho' case ( REGRIDDING_ADAPTIVE ) @@ -2341,6 +2415,8 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri case (REGRIDDING_HYCOM1) if (associated(CS%hycom_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) & call set_hycom_params(CS%hycom_CS, interp_CS=CS%interp_CS) + case (REGRIDDING_HYBGEN) + ! Do nothing for now. case (REGRIDDING_SLIGHT) if (present(min_thickness)) call set_slight_params(CS%slight_CS, min_thickness=min_thickness) if (present(dz_min_surface)) call set_slight_params(CS%slight_CS, dz_ml_min=dz_min_surface) @@ -2409,7 +2485,8 @@ function getStaticThickness( CS, SSH, depth ) real :: z, dz select case ( CS%regridding_scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) + case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, & + REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE ) if (depth>0.) then z = ssh do k = 1, CS%nk diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 82087eea24..ebe0f81743 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -516,13 +516,13 @@ subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, & u_min = min(u_l, u_c) u_max = max(u_l, u_c) if (ppoly_r_E(i0,1) < u_min) then - write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge undershoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & - 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_min + write(0,'(a,i4,5(1x,a,1pe24.16))') 'Left edge undershoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & + 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_min problem_detected = .true. endif if (ppoly_r_E(i0,1) > u_max) then - write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge overshoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & - 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_max + write(0,'(a,i4,5(1x,a,1pe24.16))') 'Left edge overshoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & + 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_max problem_detected = .true. endif endif @@ -530,27 +530,27 @@ subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, & u_min = min(u_c, u_r) u_max = max(u_c, u_r) if (ppoly_r_E(i0,2) < u_min) then - write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge undershoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & - 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_min + write(0,'(a,i4,5(1x,a,1pe24.16))') 'Right edge undershoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & + 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_min problem_detected = .true. endif if (ppoly_r_E(i0,2) > u_max) then - write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge overshoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & - 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_max + write(0,'(a,i4,5(1x,a,1pe24.16))') 'Right edge overshoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & + 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_max problem_detected = .true. endif endif if (i0 > 1) then if ( (u_c-u_l)*(ppoly_r_E(i0,1)-ppoly_r_E(i0-1,2)) < 0.) then - write(0,'(a,i4,5(x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & - 'right edge=',ppoly_r_E(i0-1,2),'left edge=',ppoly_r_E(i0,1) - write(0,'(5(a,1pe24.16,x))') 'u(i0)-u(i0-1)',u_c-u_l,'edge diff=',ppoly_r_E(i0,1)-ppoly_r_E(i0-1,2) + write(0,'(a,i4,5(1x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & + 'right edge=',ppoly_r_E(i0-1,2),'left edge=',ppoly_r_E(i0,1) + write(0,'(5(a,1pe24.16,1x))') 'u(i0)-u(i0-1)',u_c-u_l,'edge diff=',ppoly_r_E(i0,1)-ppoly_r_E(i0-1,2) problem_detected = .true. endif endif if (problem_detected) then write(0,'(a,1p9e24.16)') 'Polynomial coeffs:',ppoly_r_coefs(i0,:) - write(0,'(3(a,1pe24.16,x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r + write(0,'(3(a,1pe24.16,1x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r write(0,'(a4,10a24)') 'i0','h0(i0)','u0(i0)','left edge','right edge','Polynomial coefficients' do n = 1, n0 write(0,'(i4,1p10e24.16)') n,h0(n),u0(n),ppoly_r_E(n,1),ppoly_r_E(n,2),ppoly_r_coefs(n,:) @@ -1960,7 +1960,7 @@ logical function test_answer(verbose, n, u, u_true, label, tol) if (abs(u(k) - u_true(k)) > tolerance) test_answer = .true. enddo if (test_answer .or. verbose) then - write(stdout,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label + write(stdout,'(a4,2a24,1x,a)') 'k','Calculated value','Correct value',label do k = 1, n if (abs(u(k) - u_true(k)) > tolerance) then write(stdout,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong' diff --git a/src/ALE/regrid_consts.F90 b/src/ALE/regrid_consts.F90 index 7e8edea344..9fe638dd5b 100644 --- a/src/ALE/regrid_consts.F90 +++ b/src/ALE/regrid_consts.F90 @@ -21,6 +21,7 @@ module regrid_consts integer, parameter :: REGRIDDING_SIGMA_SHELF_ZSTAR = 8 !< Identifiered for z* coordinates at the bottom, !! sigma-near the top integer, parameter :: REGRIDDING_ADAPTIVE = 9 !< Adaptive coordinate mode identifier +integer, parameter :: REGRIDDING_HYBGEN = 10 !< Hybgen coordinates identifier character(len=*), parameter :: REGRIDDING_LAYER_STRING = "LAYER" !< Layer string character(len=*), parameter :: REGRIDDING_ZSTAR_STRING_OLD = "Z*" !< z* string (legacy name) @@ -29,6 +30,7 @@ module regrid_consts character(len=*), parameter :: REGRIDDING_SIGMA_STRING = "SIGMA" !< Sigma string character(len=*), parameter :: REGRIDDING_ARBITRARY_STRING = "ARB" !< Arbitrary coordinates character(len=*), parameter :: REGRIDDING_HYCOM1_STRING = "HYCOM1" !< Hycom string +character(len=*), parameter :: REGRIDDING_HYBGEN_STRING = "HYBGEN" !< Hybgen string character(len=*), parameter :: REGRIDDING_SLIGHT_STRING = "SLIGHT" !< Hybrid S-rho string character(len=*), parameter :: REGRIDDING_SIGMA_SHELF_ZSTAR_STRING = "SIGMA_SHELF_ZSTAR" !< Hybrid z*/sigma character(len=*), parameter :: REGRIDDING_ADAPTIVE_STRING = "ADAPTIVE" !< Adaptive coordinate string @@ -60,6 +62,7 @@ function coordinateMode(string) case (trim(REGRIDDING_RHO_STRING)); coordinateMode = REGRIDDING_RHO case (trim(REGRIDDING_SIGMA_STRING)); coordinateMode = REGRIDDING_SIGMA case (trim(REGRIDDING_HYCOM1_STRING)); coordinateMode = REGRIDDING_HYCOM1 + case (trim(REGRIDDING_HYBGEN_STRING)); coordinateMode = REGRIDDING_HYBGEN case (trim(REGRIDDING_SLIGHT_STRING)); coordinateMode = REGRIDDING_SLIGHT case (trim(REGRIDDING_ARBITRARY_STRING)); coordinateMode = REGRIDDING_ARBITRARY case (trim(REGRIDDING_SIGMA_SHELF_ZSTAR_STRING)); coordinateMode = REGRIDDING_SIGMA_SHELF_ZSTAR @@ -81,6 +84,7 @@ function coordinateUnitsI(coordMode) case (REGRIDDING_RHO); coordinateUnitsI = "kg m^-3" case (REGRIDDING_SIGMA); coordinateUnitsI = "Non-dimensional" case (REGRIDDING_HYCOM1); coordinateUnitsI = "m" + case (REGRIDDING_HYBGEN); coordinateUnitsI = "m" case (REGRIDDING_SLIGHT); coordinateUnitsI = "m" case (REGRIDDING_ADAPTIVE); coordinateUnitsI = "m" case default ; call MOM_error(FATAL, "coordinateUnts: "//& @@ -116,6 +120,7 @@ logical function state_dependent_int(mode) case (REGRIDDING_RHO); state_dependent_int = .true. case (REGRIDDING_SIGMA); state_dependent_int = .false. case (REGRIDDING_HYCOM1); state_dependent_int = .true. + case (REGRIDDING_HYBGEN); state_dependent_int = .true. case (REGRIDDING_SLIGHT); state_dependent_int = .true. case (REGRIDDING_ADAPTIVE); state_dependent_int = .true. case default ; call MOM_error(FATAL, "state_dependent: "//& diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f3d8869320..0703e45696 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3551,7 +3551,7 @@ subroutine extract_surface_state(CS, sfc_state_in) ig = i + G%HI%idg_offset ! Global i-index jg = j + G%HI%jdg_offset ! Global j-index if (use_temperature) then - write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),8(a,es11.4,x))') & + write(msg(1:240),'(2(a,i4,1x),4(a,f8.3,1x),8(a,es11.4,1x))') & 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), & @@ -3560,7 +3560,7 @@ subroutine extract_surface_state(CS, sfc_state_in) 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) else - write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') & + write(msg(1:240),'(2(a,i4,1x),4(a,f8.3,1x),6(a,es11.4))') & 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), & @@ -3577,7 +3577,7 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo ; enddo call sum_across_PEs(numberOfErrors) if (numberOfErrors>0) then - write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberOfErrors, & + write(msg(1:240),'(3(a,i9,1x))') 'There were a total of ',numberOfErrors, & 'locations detected with extreme surface values!' call MOM_error(FATAL, trim(msg)) endif diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 9a58dddd0f..085b4e02c3 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -154,7 +154,7 @@ module MOM_dynamics_unsplit !> A pointer to the set_visc control structure type(set_visc_CS), pointer :: set_visc_CSp => NULL() !> A pointer to the tidal forcing control structure - type(tidal_forcing_CS), pointer :: tides_CSp => NULL() + type(tidal_forcing_CS) :: tides_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index ec4a1aa843..91aaca508e 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -156,7 +156,7 @@ module MOM_dynamics_unsplit_RK2 !> A pointer to the set_visc control structure type(set_visc_CS), pointer :: set_visc_CSp => NULL() !> A pointer to the tidal forcing control structure - type(tidal_forcing_CS), pointer :: tides_CSp => NULL() + type(tidal_forcing_CS) :: tides_CSp !> A pointer to the ALE control structure. type(ALE_CS), pointer :: ALE_CSp => NULL() diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 22b7ee5c2b..6de46f7738 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -704,7 +704,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & fluxes%num_msg = fluxes%num_msg + 1 write(mesg,'("Penetrating shortwave of ",1pe17.10, & &" exceeds total shortwave of ",1pe17.10,& - &" at ",1pg11.4,"E, "1pg11.4,"N.")') & + &" at ",1pg11.4,",E,",1pg11.4,"N.")') & Pen_SW_tot(i), I_Cp_Hconvert*scale*dt * fluxes%sw(i,j), & G%geoLonT(i,j), G%geoLatT(i,j) call MOM_error(WARNING,mesg) @@ -3125,7 +3125,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & if (present(waves)) then; if (waves) then; if (.not. present(num_stk_bands)) then call MOM_error(FATAL,"Requested to & - initialize with waves, but no waves are present.") + &initialize with waves, but no waves are present.") endif if (num_stk_bands > 0) then if (.not.associated(forces%ustkb)) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 6ce0940ddb..bda0b2df38 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4777,7 +4777,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) if (color(i,j) /= color2(i,j)) then fatal_error = .True. write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", & - "the masking of the outside grid points.")') i, j + &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM register_tracer: "//mesg, all_print=.true.) endif if (color(i,j) == cout) G%bathyT(i,j) = min_depth diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index a4badaf8e7..e52feec697 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -147,7 +147,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday) write (file,'(/,"--------------------------")') write (file,'(/,"Time ",i5,i4,F6.2," U-velocity violation at ",I4,": ",2(I3), & - & " (",F7.2," E "F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') & + & " (",F7.2," E ",F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') & yr, yearday, (REAL(sec)/3600.0), pe_here(), I, j, & G%geoLonCu(I,j), G%geoLatCu(I,j), ks, ke, US%T_to_s*dt @@ -156,174 +156,174 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st if ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*Angstrom) do_k(k) = .true. enddo - write(file,'(/,"Layers:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k) ; enddo - write(file,'(/,"u(m): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*um(I,j,k)) ; enddo + write(file,'(/,"Layers:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(I10," ")', advance='no') (k) ; enddo + write(file,'(/,"u(m): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*um(I,j,k)) ; enddo if (prev_avail) then - write(file,'(/,"u(mp): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%u_prev(I,j,k)) ; enddo + write(file,'(/,"u(mp): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%u_prev(I,j,k)) ; enddo endif - write(file,'(/,"u(3): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%u_av(I,j,k)) ; enddo + write(file,'(/,"u(3): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%u_av(I,j,k)) ; enddo - write(file,'(/,"CFL u: ",$)') + write(file,'(/,"CFL u: ")', advance='no') do k=ks,ke ; if (do_k(k)) then CFL = abs(um(I,j,k)) * dt * G%dy_Cu(I,j) if (um(I,j,k) < 0.0) then ; CFL = CFL * G%IareaT(i+1,j) else ; CFL = CFL * G%IareaT(i,j) ; endif - write(file,'(ES10.3," ",$)') CFL + write(file,'(ES10.3," ")', advance='no') CFL endif ; enddo - write(file,'(/,"CFL0 u:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"CFL0 u:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & abs(um(I,j,k)) * dt * G%IdxCu(I,j) ; enddo if (prev_avail) then - write(file,'(/,"du: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"du: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*(um(I,j,k)-CS%u_prev(I,j,k))) ; enddo endif - write(file,'(/,"CAu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%CAu(I,j,k)) ; enddo - write(file,'(/,"PFu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%PFu(I,j,k)) ; enddo - write(file,'(/,"diffu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%diffu(I,j,k)) ; enddo + write(file,'(/,"CAu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%CAu(I,j,k)) ; enddo + write(file,'(/,"PFu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%PFu(I,j,k)) ; enddo + write(file,'(/,"diffu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%diffu(I,j,k)) ; enddo if (associated(ADp%gradKEu)) then - write(file,'(/,"KEu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"KEu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*ADp%gradKEu(I,j,k)) ; enddo endif if (associated(ADp%rv_x_v)) then - write(file,'(/,"Coru: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"Coru: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & vel_scale*dt*(ADp%CAu(I,j,k)-ADp%rv_x_v(I,j,k)) ; enddo endif if (associated(ADp%du_dt_visc)) then - write(file,'(/,"ubv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"ubv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & vel_scale*(um(I,j,k) - dt*ADp%du_dt_visc(I,j,k)) ; enddo - write(file,'(/,"duv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"duv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*ADp%du_dt_visc(I,j,k)) ; enddo endif if (associated(ADp%du_other)) then - write(file,'(/,"du_other: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"du_other: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*ADp%du_other(I,j,k)) ; enddo endif if (present(a)) then - write(file,'(/,"a: ",$)') - do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') (US%Z_to_m*a(I,j,k)*dt) ; enddo + write(file,'(/,"a: ")', advance='no') + do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,k)*dt) ; enddo endif if (present(hv)) then - write(file,'(/,"hvel: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hv(I,j,k) ; enddo + write(file,'(/,"hvel: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(I,j,k) ; enddo endif write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) if (associated(CS%u_accel_bt)) then - write(file,'("dubt: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'("dubt: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*CS%u_accel_bt(I,j,k)) ; enddo write(file,'(/)') endif - write(file,'(/,"h--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j-1,k)) ; enddo - write(file,'(/,"h+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j-1,k)) ; enddo - write(file,'(/,"h-0: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j,k)) ; enddo - write(file,'(/,"h+0: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j,k)) ; enddo - write(file,'(/,"h-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j+1,k)) ; enddo - write(file,'(/,"h++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j+1,k)) ; enddo + write(file,'(/,"h--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i,j-1,k)) ; enddo + write(file,'(/,"h+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i+1,j-1,k)) ; enddo + write(file,'(/,"h-0: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i,j,k)) ; enddo + write(file,'(/,"h+0: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i+1,j,k)) ; enddo + write(file,'(/,"h-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i,j+1,k)) ; enddo + write(file,'(/,"h++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i+1,j+1,k)) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo - write(file,'(/,"e-: ",$)') - write(file,'(ES10.3," ",$)') e(ks) - do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K) ; enddo + write(file,'(/,"e-: ")', advance='no') + write(file,'(ES10.3," ")', advance='no') e(ks) + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i+1,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i+1,j,k) ; enddo - write(file,'(/,"e+: ",$)') - write(file,'(ES10.3," ",$)') e(ks) - do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K) ; enddo + write(file,'(/,"e+: ")', advance='no') + write(file,'(ES10.3," ")', advance='no') e(ks) + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo if (associated(CS%T)) then - write(file,'(/,"T-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j,k) ; enddo - write(file,'(/,"T+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i+1,j,k) ; enddo + write(file,'(/,"T-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j,k) ; enddo + write(file,'(/,"T+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i+1,j,k) ; enddo endif if (associated(CS%S)) then - write(file,'(/,"S-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j,k) ; enddo - write(file,'(/,"S+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i+1,j,k) ; enddo + write(file,'(/,"S-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j,k) ; enddo + write(file,'(/,"S+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i+1,j,k) ; enddo endif if (prev_avail) then - write(file,'(/,"v--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_prev(i,J-1,k)) ; enddo - write(file,'(/,"v-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_prev(i,J,k)) ; enddo - write(file,'(/,"v+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_prev(i+1,J-1,k)) ; enddo - write(file,'(/,"v++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_prev(i+1,J,k)) ; enddo - endif - - write(file,'(/,"vh--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"v--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_prev(i,J-1,k)) ; enddo + write(file,'(/,"v-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_prev(i,J,k)) ; enddo + write(file,'(/,"v+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_prev(i+1,J-1,k)) ; enddo + write(file,'(/,"v++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_prev(i+1,J,k)) ; enddo + endif + + write(file,'(/,"vh--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%vh(i,J-1,k)*G%IdxCv(i,J-1)) ; enddo - write(file,'(/," vhC--:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhC--:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_av(i,j-1,k)*uh_scale*(hin(i,j-1,k) + hin(i,j,k))) ; enddo if (prev_avail) then - write(file,'(/," vhCp--:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhCp--:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_prev(i,j-1,k)*uh_scale*(hin(i,j-1,k) + hin(i,j,k))) ; enddo endif - write(file,'(/,"vh-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"vh-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%vh(i,J,k)*G%IdxCv(i,J)) ; enddo - write(file,'(/," vhC-+:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhC-+:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_av(i,J,k)*uh_scale*(hin(i,j,k) + hin(i,j+1,k))) ; enddo if (prev_avail) then - write(file,'(/," vhCp-+:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhCp-+:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_prev(i,J,k)*uh_scale*(hin(i,j,k) + hin(i,j+1,k))) ; enddo endif - write(file,'(/,"vh+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"vh+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%vh(i+1,J-1,k)*G%IdxCv(i+1,J-1)) ; enddo - write(file,'(/," vhC+-:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhC+-:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_av(i+1,J-1,k)*uh_scale*(hin(i+1,j-1,k) + hin(i+1,j,k))) ; enddo if (prev_avail) then - write(file,'(/," vhCp+-:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhCp+-:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_prev(i+1,J-1,k)*uh_scale*(hin(i+1,j-1,k) + hin(i+1,j,k))) ; enddo endif - write(file,'(/,"vh++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"vh++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%vh(i+1,J,k)*G%IdxCv(i+1,J)) ; enddo - write(file,'(/," vhC++:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhC++:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_av(i+1,J,k)*uh_scale*(hin(i+1,j,k) + hin(i+1,j+1,k))) ; enddo if (prev_avail) then - write(file,'(/," vhCp++:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," vhCp++:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (0.5*CS%v_av(i+1,J,k)*uh_scale*(hin(i+1,j,k) + hin(i+1,j+1,k))) ; enddo endif @@ -337,48 +337,48 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st Inorm(k) = 1.0 / du enddo - write(file,'(2/,"Norm: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (vel_scale / Inorm(k)) ; enddo + write(file,'(2/,"Norm: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') (vel_scale / Inorm(k)) ; enddo - write(file,'(/,"du: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"du: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & ((um(I,j,k)-CS%u_prev(I,j,k)) * Inorm(k)) ; enddo - write(file,'(/,"CAu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"CAu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%CAu(I,j,k) * Inorm(k)) ; enddo - write(file,'(/,"PFu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"PFu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%PFu(I,j,k) * Inorm(k)) ; enddo - write(file,'(/,"diffu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"diffu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%diffu(I,j,k) * Inorm(k)) ; enddo if (associated(ADp%gradKEu)) then - write(file,'(/,"KEu: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"KEu: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%gradKEu(I,j,k) * Inorm(k)) ; enddo endif if (associated(ADp%rv_x_v)) then - write(file,'(/,"Coru: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"Coru: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*(ADp%CAu(I,j,k)-ADp%rv_x_v(I,j,k)) * Inorm(k)) ; enddo endif if (associated(ADp%du_dt_visc)) then - write(file,'(/,"duv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"duv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%du_dt_visc(I,j,k) * Inorm(k)) ; enddo endif if (associated(ADp%du_other)) then - write(file,'(/,"du_other: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"du_other: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (ADp%du_other(I,j,k) * Inorm(k)) ; enddo endif if (associated(CS%u_accel_bt)) then - write(file,'(/,"dubt: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"dubt: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*CS%u_accel_bt(I,j,k) * Inorm(k)) ; enddo endif endif @@ -487,178 +487,178 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st if ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*Angstrom) do_k(k) = .true. enddo - write(file,'(/,"Layers:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k) ; enddo - write(file,'(/,"v(m): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*vm(i,J,k)) ; enddo + write(file,'(/,"Layers:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(I10," ")', advance='no') (k) ; enddo + write(file,'(/,"v(m): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*vm(i,J,k)) ; enddo if (prev_avail) then - write(file,'(/,"v(mp): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_prev(i,J,k)) ; enddo + write(file,'(/,"v(mp): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_prev(i,J,k)) ; enddo endif - write(file,'(/,"v(3): ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*CS%v_av(i,J,k)) ; enddo - write(file,'(/,"CFL v: ",$)') + write(file,'(/,"v(3): ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*CS%v_av(i,J,k)) ; enddo + write(file,'(/,"CFL v: ")', advance='no') do k=ks,ke ; if (do_k(k)) then CFL = abs(vm(i,J,k)) * dt * G%dx_Cv(i,J) if (vm(i,J,k) < 0.0) then ; CFL = CFL * G%IareaT(i,j+1) else ; CFL = CFL * G%IareaT(i,j) ; endif - write(file,'(ES10.3," ",$)') CFL + write(file,'(ES10.3," ")', advance='no') CFL endif ; enddo - write(file,'(/,"CFL0 v:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"CFL0 v:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & abs(vm(i,J,k)) * dt * G%IdyCv(i,J) ; enddo if (prev_avail) then - write(file,'(/,"dv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"dv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*(vm(i,J,k)-CS%v_prev(i,J,k))) ; enddo endif - write(file,'(/,"CAv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%CAv(i,J,k)) ; enddo + write(file,'(/,"CAv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%CAv(i,J,k)) ; enddo - write(file,'(/,"PFv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%PFv(i,J,k)) ; enddo + write(file,'(/,"PFv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%PFv(i,J,k)) ; enddo - write(file,'(/,"diffv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vel_scale*dt*ADp%diffv(i,J,k)) ; enddo + write(file,'(/,"diffv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (vel_scale*dt*ADp%diffv(i,J,k)) ; enddo if (associated(ADp%gradKEv)) then - write(file,'(/,"KEv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"KEv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*ADp%gradKEv(i,J,k)) ; enddo endif if (associated(ADp%rv_x_u)) then - write(file,'(/,"Corv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"Corv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & vel_scale*dt*(ADp%CAv(i,J,k)-ADp%rv_x_u(i,J,k)) ; enddo endif if (associated(ADp%dv_dt_visc)) then - write(file,'(/,"vbv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"vbv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & vel_scale*(vm(i,J,k) - dt*ADp%dv_dt_visc(i,J,k)) ; enddo - write(file,'(/,"dvv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"dvv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*ADp%dv_dt_visc(i,J,k)) ; enddo endif if (associated(ADp%dv_other)) then - write(file,'(/,"dv_other: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"dv_other: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*ADp%dv_other(i,J,k)) ; enddo endif if (present(a)) then - write(file,'(/,"a: ",$)') - do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') (US%Z_to_m*a(i,j,k)*dt) ; enddo + write(file,'(/,"a: ")', advance='no') + do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,j,k)*dt) ; enddo endif if (present(hv)) then - write(file,'(/,"hvel: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hv(i,J,k) ; enddo + write(file,'(/,"hvel: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(i,J,k) ; enddo endif write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0) if (associated(CS%v_accel_bt)) then - write(file,'("dvbt: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'("dvbt: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (vel_scale*dt*CS%v_accel_bt(i,J,k)) ; enddo write(file,'(/)') endif - write(file,'("h--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i-1,j,k) ; enddo - write(file,'(/,"h0-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i,j,k) ; enddo - write(file,'(/,"h+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i+1,j,k) ; enddo - write(file,'(/,"h-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i-1,j+1,k) ; enddo - write(file,'(/,"h0+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i,j+1,k) ; enddo - write(file,'(/,"h++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i+1,j+1,k) ; enddo + write(file,'("h--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i-1,j,k) ; enddo + write(file,'(/,"h0-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i,j,k) ; enddo + write(file,'(/,"h+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i+1,j,k) ; enddo + write(file,'(/,"h-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i-1,j+1,k) ; enddo + write(file,'(/,"h0+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i,j+1,k) ; enddo + write(file,'(/,"h++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i+1,j+1,k) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo - write(file,'(/,"e-: ",$)') - write(file,'(ES10.3," ",$)') e(ks) - do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K) ; enddo + write(file,'(/,"e-: ")', advance='no') + write(file,'(ES10.3," ")', advance='no') e(ks) + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo e(nz+1) = -US%Z_to_m*(G%bathyT(i,j+1) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j+1,k) ; enddo - write(file,'(/,"e+: ",$)') - write(file,'(ES10.3," ",$)') e(ks) - do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K) ; enddo + write(file,'(/,"e+: ")', advance='no') + write(file,'(ES10.3," ")', advance='no') e(ks) + do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo if (associated(CS%T)) then - write(file,'(/,"T-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j,k) ; enddo - write(file,'(/,"T+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j+1,k) ; enddo + write(file,'(/,"T-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j,k) ; enddo + write(file,'(/,"T+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j+1,k) ; enddo endif if (associated(CS%S)) then - write(file,'(/,"S-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j,k) ; enddo - write(file,'(/,"S+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j+1,k) ; enddo + write(file,'(/,"S-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j,k) ; enddo + write(file,'(/,"S+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j+1,k) ; enddo endif if (prev_avail) then - write(file,'(/,"u--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') vel_scale*CS%u_prev(I-1,j,k) ; enddo - write(file,'(/,"u-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') vel_scale*CS%u_prev(I-1,j+1,k) ; enddo - write(file,'(/,"u+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') vel_scale*CS%u_prev(I,j,k) ; enddo - write(file,'(/,"u++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') vel_scale*CS%u_prev(I,j+1,k) ; enddo - endif - - write(file,'(/,"uh--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"u--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') vel_scale*CS%u_prev(I-1,j,k) ; enddo + write(file,'(/,"u-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') vel_scale*CS%u_prev(I-1,j+1,k) ; enddo + write(file,'(/,"u+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') vel_scale*CS%u_prev(I,j,k) ; enddo + write(file,'(/,"u++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') vel_scale*CS%u_prev(I,j+1,k) ; enddo + endif + + write(file,'(/,"uh--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%uh(I-1,j,k)*G%IdyCu(I-1,j)) ; enddo - write(file,'(/," uhC--: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhC--: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_av(I-1,j,k) * uh_scale*0.5*(hin(i-1,j,k) + hin(i,j,k))) ; enddo if (prev_avail) then - write(file,'(/," uhCp--:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhCp--:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_prev(I-1,j,k) * uh_scale*0.5*(hin(i-1,j,k) + hin(i,j,k))) ; enddo endif - write(file,'(/,"uh-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"uh-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%uh(I-1,j+1,k)*G%IdyCu(I-1,j+1)) ; enddo - write(file,'(/," uhC-+: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhC-+: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_av(I-1,j+1,k) * uh_scale*0.5*(hin(i-1,j+1,k) + hin(i,j+1,k))) ; enddo if (prev_avail) then - write(file,'(/," uhCp-+:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhCp-+:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_prev(I-1,j+1,k) * uh_scale*0.5*(hin(i-1,j+1,k) + hin(i,j+1,k))) ; enddo endif - write(file,'(/,"uh+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"uh+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%uh(I,j,k)*G%IdyCu(I,j)) ; enddo - write(file,'(/," uhC+-: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhC+-: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_av(I,j,k) * uh_scale*0.5*(hin(i,j,k) + hin(i+1,j,k))) ; enddo if (prev_avail) then - write(file,'(/," uhCp+-:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhCp+-:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_prev(I,j,k) * uh_scale*0.5*(hin(i,j,k) + hin(i+1,j,k))) ; enddo endif - write(file,'(/,"uh++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/,"uh++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (uh_scale*CDp%uh(I,j+1,k)*G%IdyCu(I,j+1)) ; enddo - write(file,'(/," uhC++: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhC++: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_av(I,j+1,k) * uh_scale*0.5*(hin(i,j+1,k) + hin(i+1,j+1,k))) ; enddo if (prev_avail) then - write(file,'(/," uhCp++:",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') & + write(file,'(/," uhCp++:")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') & (CS%u_prev(I,j+1,k) * uh_scale*0.5*(hin(i,j+1,k) + hin(i+1,j+1,k))) ; enddo endif @@ -672,44 +672,44 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st Inorm(k) = 1.0 / dv enddo - write(file,'(2/,"Norm: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (vel_scale / Inorm(k)) ; enddo - write(file,'(/,"dv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(2/,"Norm: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') (vel_scale / Inorm(k)) ; enddo + write(file,'(/,"dv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & ((vm(i,J,k)-CS%v_prev(i,J,k)) * Inorm(k)) ; enddo - write(file,'(/,"CAv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"CAv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%CAv(i,J,k) * Inorm(k)) ; enddo - write(file,'(/,"PFv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"PFv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%PFv(i,J,k) * Inorm(k)) ; enddo - write(file,'(/,"diffv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"diffv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%diffv(i,J,k) * Inorm(k)) ; enddo if (associated(ADp%gradKEu)) then - write(file,'(/,"KEv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"KEv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%gradKEv(i,J,k) * Inorm(k)) ; enddo endif if (associated(ADp%rv_x_u)) then - write(file,'(/,"Corv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"Corv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*(ADp%CAv(i,J,k)-ADp%rv_x_u(i,J,k)) * Inorm(k)) ; enddo endif if (associated(ADp%dv_dt_visc)) then - write(file,'(/,"dvv: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"dvv: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*ADp%dv_dt_visc(i,J,k) * Inorm(k)) ; enddo endif if (associated(ADp%dv_other)) then - write(file,'(/,"dv_other: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"dv_other: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (ADp%dv_other(i,J,k) * Inorm(k)) ; enddo endif if (associated(CS%v_accel_bt)) then - write(file,'(/,"dvbt: ",$)') - do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') & + write(file,'(/,"dvbt: ")', advance='no') + do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ")', advance='no') & (dt*CS%v_accel_bt(i,J,k) * Inorm(k)) ; enddo endif endif diff --git a/src/diagnostics/MOM_debugging.F90 b/src/diagnostics/MOM_debugging.F90 index fda5a97d69..3b24a3871b 100644 --- a/src/diagnostics/MOM_debugging.F90 +++ b/src/diagnostics/MOM_debugging.F90 @@ -196,7 +196,7 @@ subroutine check_redundant_vC2d(mesg, u_comp, v_comp, G, is, ie, js, je, & if (v_resym(i,j) /= v_comp(i,j) .and. & redundant_prints(3) < max_redundant_prints) then write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", & - & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') & + & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)," on pe ",i4)') & v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, & G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() write(0,'(A155)') trim(mesg)//trim(mesg2) @@ -386,7 +386,7 @@ subroutine check_redundant_vB2d(mesg, u_comp, v_comp, G, is, ie, js, je, & if (v_resym(i,j) /= v_comp(i,j) .and. & redundant_prints(2) < max_redundant_prints) then write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", & - & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') & + & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)," on pe ",i4)') & v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, & G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() write(0,'(A155)') trim(mesg)//trim(mesg2) @@ -549,7 +549,7 @@ subroutine check_redundant_vT2d(mesg, u_comp, v_comp, G, is, ie, js, je, & if (v_nonsym(i,j) /= v_comp(i,j) .and. & redundant_prints(1) < max_redundant_prints) then write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", & - & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') & + & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)," on pe ",i4)') & v_comp(i,j), v_nonsym(i,j),v_comp(i,j)-v_nonsym(i,j),i,j, & G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() write(0,'(A155)') trim(mesg)//trim(mesg2) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index a7cae98620..48097b40c4 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -848,13 +848,13 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci endif do m=1,nTr_stocks - write(stdout,'(" Total ",a,": ",ES24.16,X,a)') & + write(stdout,'(" Total ",a,": ",ES24.16,1X,a)') & trim(Tr_names(m)), Tr_stocks(m), trim(Tr_units(m)) if (Tr_minmax_avail(m)) then - write(stdout,'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & + write(stdout,'(64X,"Global Min:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & Tr_min(m),Tr_min_x(m),Tr_min_y(m),Tr_min_z(m) - write(stdout,'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & + write(stdout,'(64X,"Global Max:",ES24.16,1X,"at: (",f7.2,",",f7.2,",",f8.2,")" )') & Tr_max(m),Tr_max_x(m),Tr_max_y(m),Tr_max_z(m) endif diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index d1a8102fc1..15db6abce9 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -2196,7 +2196,7 @@ subroutine chk_sum_msg1(fmsg, bc0, mesg, iounit) integer, intent(in) :: iounit !< Checksum logger IO unit if (is_root_pe()) & - write(iounit, '(A,1(A,I10,X),A)') fmsg, " c=", bc0, trim(mesg) + write(iounit, '(a,1(a,i10,1x),a)') fmsg, " c=", bc0, trim(mesg) end subroutine chk_sum_msg1 !> Write a message including checksums of non-shifted and diagonally shifted arrays diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 92d53330a7..1ed89ec47e 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3732,7 +3732,7 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string mesg = '"'//trim(field_name)//'" [Unused]' endif if (len(trim((comment)))>0) then - write(diag_CS%available_diag_doc_unit, '(a,x,"(",a,")")') trim(mesg),trim(comment) + write(diag_CS%available_diag_doc_unit, '(a,1x,"(",a,")")') trim(mesg),trim(comment) else write(diag_CS%available_diag_doc_unit, '(a)') trim(mesg) endif @@ -3754,7 +3754,7 @@ subroutine log_chksum_diag(docunit, description, chksum) character(len=*), intent(in) :: description !< Name of the diagnostic module integer, intent(in) :: chksum !< chksum of the diagnostic - write(docunit, '(a,x,i9.8)') description, chksum + write(docunit, '(a,1x,i9.8)') description, chksum flush(docunit) end subroutine log_chksum_diag diff --git a/src/framework/MOM_diag_vkernels.F90 b/src/framework/MOM_diag_vkernels.F90 index 3d6e3e3f65..24ff07e7d0 100644 --- a/src/framework/MOM_diag_vkernels.F90 +++ b/src/framework/MOM_diag_vkernels.F90 @@ -311,8 +311,8 @@ logical function test_interp(verbose, missing_value, msg, nsrc, h_src, u_src, nd if (error==0.) then write(stdout,'(i3,3(1pe24.16))') k,u_dest(k),u_true(k),u_dest(k)-u_true(k) else - write(stdout,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' - write(stderr,'(i3,3(1pe24.16),x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' + write(stdout,'(i3,3(1pe24.16),1x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' + write(stderr,'(i3,3(1pe24.16),1x,a)') k,u_dest(k),u_true(k),u_dest(k)-u_true(k),'<--- WRONG!' endif enddo endif @@ -350,8 +350,8 @@ logical function test_reintegrate(verbose, missing_value, msg, nsrc, h_src, uh_s if (error==0.) then write(stdout,'(i3,3(1pe24.16))') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k) else - write(stdout,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' - write(stderr,'(i3,3(1pe24.16),x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' + write(stdout,'(i3,3(1pe24.16),1x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' + write(stderr,'(i3,3(1pe24.16),1x,a)') k,uh_dest(k),uh_true(k),uh_dest(k)-uh_true(k),'<--- WRONG!' endif enddo endif diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index dc6c0a8996..00c42727e1 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -343,7 +343,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & ! Idiot check that fewer PEs than columns have been requested if (layout(1)*layout(2) > n_global(1)*n_global(2)) then - write(mesg,'(a,2(i5,x,a))') 'You requested to use',layout(1)*layout(2), & + write(mesg,'(a,2(i5,1x,a))') 'You requested to use', layout(1)*layout(2), & 'PEs but there are only', n_global(1)*n_global(2), 'columns in the model' call MOM_error(FATAL, mesg) endif diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index de511688a9..5a125d5abd 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -66,7 +66,7 @@ subroutine myStats(array, missing, is, ie, js, je, k, mesg) call min_across_PEs(minA) call max_across_PEs(maxA) if (is_root_pe()) then - write(lMesg(1:120),'(2(a,es12.4),a,i3,x,a)') & + write(lMesg(1:120),'(2(a,es12.4),a,i3,1x,a)') & 'init_from_Z: min=',minA,' max=',maxA,' Level=',k,trim(mesg) call MOM_mesg(lMesg,2) endif diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index bef78a433a..cf649c32f7 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -7,6 +7,7 @@ module MOM_random use MOM_time_manager, only : time_type, set_date, get_date use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use iso_fortran_env, only : int32 implicit none ; private @@ -20,11 +21,13 @@ module MOM_random public :: random_unit_tests ! Private period parameters for the Mersenne Twister -integer, parameter :: blockSize = 624, & !< Size of the state vector - M = 397, & !< Pivot element in state vector - MATRIX_A = -1727483681, & !< constant vector a (0x9908b0dfUL) - UMASK = -2147483648_8, & !< most significant w-r bits (0x80000000UL) - LMASK = 2147483647 !< least significant r bits (0x7fffffffUL) +integer, parameter :: & + blockSize = 624, & !< Size of the state vector + M = 397, & !< Pivot element in state vector + MATRIX_A = -1727483681, & !< constant vector a (0x9908b0dfUL) + UMASK = ibset(0, 31), & !< most significant w-r bits (0x80000000UL) + LMASK = 2147483647 !< least significant r bits (0x7fffffffUL) + ! Private tempering parameters for the Mersenne Twister integer, parameter :: TMASKB= -1658038656, & !< (0x9d2c5680UL) TMASKC= -272236544 !< (0xefc60000UL) diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index 9df994448b..cc43ec19ea 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -200,7 +200,7 @@ subroutine write_cputime(day, n, CS, nmax, call_end) (CS%startup_cputime / CLOCKS_PER_SEC), num_pes() write(CS%fileCPU_ascii,*)" Day, Step number, CPU time, CPU time change" endif - write(CS%fileCPU_ascii,'(F12.3,", "I11,", ", F12.3,", ", F12.3)') & + write(CS%fileCPU_ascii,'(F12.3,", ",I11,", ",F12.3,", ",F12.3)') & reday, n, (CS%cputime2 / real(CLOCKS_PER_SEC)), & d_cputime / real(CLOCKS_PER_SEC) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 608a672a46..e1ecd45347 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -2901,10 +2901,10 @@ end subroutine bilinear_shape_fn_grid subroutine bilinear_shape_functions_subgrid(Phisub, nsub) + integer, intent(in) :: nsub !< The number of subgridscale quadrature locations in each direction real, dimension(nsub,nsub,2,2,2,2), & intent(inout) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] - integer, intent(in) :: nsub !< The number of subgridscale quadrature locations in each direction ! this subroutine is a helper for interpolation of floatation condition ! for the purposes of evaluating the terms \int (u,v) \phi_i dx dy in a cell that is diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 4e8dd4f186..ee3052f166 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -88,12 +88,10 @@ module MOM_state_initialization use dumbbell_initialization, only : dumbbell_initialize_sponges use MOM_tracer_Z_init, only : tracer_Z_init_array, determine_temperature use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord -use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid, ALE_regrid_accelerated -use MOM_ALE, only : TS_PLM_edge_values +use MOM_ALE, only : ALE_remap_scalar, ALE_regrid_accelerated, TS_PLM_edge_values use MOM_regridding, only : regridding_CS, set_regrid_params, getCoordinateResolution -use MOM_regridding, only : regridding_main -use MOM_remapping, only : remapping_CS, initialize_remapping -use MOM_remapping, only : remapping_core_h +use MOM_regridding, only : regridding_main, regridding_preadjust_reqs, convective_adjustment +use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field @@ -212,7 +210,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & use_EOS = associated(tv%eqn_of_state) use_OBC = associated(OBC) if (use_EOS) eos => tv%eqn_of_state - use_ice_shelf=PRESENT(frac_shelf_h) + use_ice_shelf = PRESENT(frac_shelf_h) !==================================================================== ! Initialize temporally evolving fields, either as initial @@ -2439,18 +2437,21 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! data arrays real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m] real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3] - real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z + real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [degC] + real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [ppt] + real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim] real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor ! relative to the surface [Z ~> m]. - integer, dimension(SZI_(G),SZJ_(G)) :: nlevs + integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa]. ! Local variables for ALE remapping real, dimension(:), allocatable :: hTarget ! Target thicknesses [Z ~> m]. - real, dimension(:,:,:), allocatable, target :: tmpT1dIn, tmpS1dIn - real, dimension(:,:,:), allocatable :: tmp_mask_in + real, dimension(:,:,:), allocatable, target :: tmpT1dIn ! Input temperatures on a model-sized grid [degC] + real, dimension(:,:,:), allocatable, target :: tmpS1dIn ! Input salinities on a model-sized grid [ppt] + real, dimension(:,:,:), allocatable :: tmp_mask_in ! The valid data mask on a model-sized grid [nondim] real, dimension(:,:,:), allocatable :: h1 ! Thicknesses [H ~> m or kg m-2]. real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to regridding real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m]. @@ -2459,7 +2460,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg logical :: answers_2018, default_2018_answers, hor_regrid_answers_2018 - logical :: use_ice_shelf logical :: pre_gridded logical :: separate_mixed_layer ! If true, handle the mixed layers differently. logical :: density_extrap_bug ! If true use an expression with a vertical indexing bug for @@ -2467,7 +2467,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! from data when finding the initial interface locations in ! layered mode from a dataset of T and S. character(len=64) :: remappingScheme - real :: tempAvg, saltAvg + real :: tempAvg ! Spatially averaged temperatures on a layer [degC] + real :: saltAvg ! Spatially averaged salinities on a layer [ppt] + logical :: do_conv_adj, ignore integer :: nPoints, ans integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE @@ -2493,8 +2495,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just reentrant_x = .false. ; call get_param(PF, mdl, "REENTRANT_X", reentrant_x, default=.true.) tripolar_n = .false. ; call get_param(PF, mdl, "TRIPOLAR_N", tripolar_n, default=.false.) - use_ice_shelf = present(frac_shelf_h) - call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE", filename, & "The name of the z-space input file used to initialize "//& "temperatures (T) and salinities (S). If T and S are not "//& @@ -2722,11 +2722,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just GV_loc = GV GV_loc%ke = nkd allocate( dz_interface(isd:ied,jsd:jed,nkd+1) ) ! Need for argument to regridding_main() but is not used - if (use_ice_shelf) then - call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, frac_shelf_h ) - else - call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface ) - endif + + call regridding_preadjust_reqs(regridCS, do_conv_adj, ignore) + if (do_conv_adj) call convective_adjustment(G, GV_loc, h1, tv_loc) + call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, conv_adjust=.false., & + frac_shelf_h=frac_shelf_h ) + deallocate( dz_interface ) endif call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, & diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index d8da752a1f..e1dff4f5bd 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -178,7 +178,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2]. real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer ! densities [R L2 T-2 ~> Pa]. - real :: aFac, bFac ! Nondimensional ratios [nondim] + real :: aFac, bFac ! Nondimensional ratios [nondim] real :: ddRho ! A density difference [R ~> kg m-3] real :: hAtVel ! Thickness at the velocity points [H ~> m or kg m-2] real :: zpa ! Fractional position within the mixed layer of the interface above a layer [nondim] @@ -190,15 +190,6 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - real :: PSI, PSI1, z, BOTTOP, XP, DD ! For the following statement functions [nondim] - ! Stream function as a function of non-dimensional position within mixed-layer (F77 statement function) - !PSI1(z) = max(0., (1. - (2.*z+1.)**2 ) ) - PSI1(z) = max(0., (1. - (2.*z+1.)**2 ) * (1. + (5./21.)*(2.*z+1.)**2) ) - BOTTOP(z) = 0.5*(1.-SIGN(1.,z+0.5)) ! =0 for z>-0.5, =1 for z<-0.5 - XP(z) = max(0., min(1., (-z-0.5)*2./(1.+2.*CS%MLE_tail_dh) ) ) - DD(z) = (1.-3.*(XP(z)**2)+2.*(XP(z)**3))**(1.+2.*CS%MLE_tail_dh) - PSI(z) = max( PSI1(z), DD(z)*BOTTOP(z) ) ! Combines original PSI1 with tail - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -557,6 +548,22 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var ! This needs to happen after the H update and before the next post_data. call diag_update_remap_grids(CS%diag) +contains + !> Stream function as a function of non-dimensional position within mixed-layer + real function psi(z) + real, intent(in) :: z !< Fractional mixed layer depth [nondim] + real :: psi1, bottop, xp, dd + + !psi1 = max(0., (1. - (2.*z + 1.)**2)) + psi1 = max(0., (1. - (2.*z + 1.)**2) * (1. + (5./21.)*(2.*z + 1.)**2)) + + xp = max(0., min(1., (-z - 0.5)*2. / (1. + 2.*CS%MLE_tail_dh))) + dd = (1. - 3.*(xp**2) + 2.*(xp**3))**(1. + 2.*CS%MLE_tail_dh) + bottop = 0.5*(1. - sign(1., z + 0.5)) ! =0 for z>-0.5, =1 for z<-0.5 + + psi = max(psi1, dd*bottop) ! Combines original psi1 with tail + end function psi + end subroutine mixedlayer_restrat_general diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index a128390ae6..446cdb8b45 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -691,7 +691,7 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, CS%fldno = CS%fldno + 1 if (CS%fldno > MAX_FIELDS_) then write(mesg, '("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease "//& - "the number of fields to be damped in the call to initialize_ALE_sponge." )') CS%fldno + &"the number of fields to be damped in the call to initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif ! get a unique time interp id for this field. If sponge data is on-grid, then setup diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 227e3ffb06..7296640560 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -343,8 +343,8 @@ end subroutine swap !> Receives a 1D array x and sorts it into ascending order. subroutine sort(x, n) - real, dimension(n), intent(inout) :: x !< 1D array to be sorted integer, intent(in ) :: n !< # of pts in the array + real, dimension(n), intent(inout) :: x !< 1D array to be sorted ! local variables integer :: i, location @@ -1012,8 +1012,8 @@ logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_a write(stdout,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans endif - 20 format(A,"=",i3,X,A,"=",i3) - 30 format(A,"=",f20.16,X,A,"=",f20.16) + 20 format(A,"=",i3,1X,A,"=",i3) + 30 format(A,"=",f20.16,1X,A,"=",f20.16) end function test_boundary_k_range diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 6bafb064b8..08361ba9af 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -1364,7 +1364,7 @@ subroutine find_neutral_surface_positions_discontinuous(CS, nk, & if (CS%debug) then write(stdout,'(A,I2)') "Searching left layer ", kl_left - write(stdout,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right + write(stdout,'(A,I2,1X,I2)') "Searching from right: ", kl_right, ki_right write(stdout,*) "Temp/Salt Reference: ", Tr(kl_right,ki_right), Sr(kl_right,ki_right) write(stdout,*) "Temp/Salt Top L: ", Tl(kl_left,1), Sl(kl_left,1) write(stdout,*) "Temp/Salt Bot L: ", Tl(kl_left,2), Sl(kl_left,2) @@ -1387,7 +1387,7 @@ subroutine find_neutral_surface_positions_discontinuous(CS, nk, & if (CS%debug) then write(stdout,'(A,I2)') "Searching right layer ", kl_right - write(stdout,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left + write(stdout,'(A,I2,1X,I2)') "Searching from left: ", kl_left, ki_left write(stdout,*) "Temp/Salt Reference: ", Tl(kl_left,ki_left), Sl(kl_left,ki_left) write(stdout,*) "Temp/Salt Top L: ", Tr(kl_right,1), Sr(kl_right,1) write(stdout,*) "Temp/Salt Bot L: ", Tr(kl_right,2), Sr(kl_right,2) @@ -2651,9 +2651,9 @@ logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, ti if (test_fv_diff) stdunit = stderr ! In case of wrong results, write to error stream write(stdunit,'(a)') title if (test_fv_diff) then - write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' + write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' else - write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue + write(stdunit,'(2(1x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue endif endif @@ -2683,9 +2683,9 @@ logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue if (test_fvlsq_slope) stdunit = stderr ! In case of wrong results, write to error stream write(stdunit,'(a)') title if (test_fvlsq_slope) then - write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' + write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' else - write(stdunit,'(2(x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue + write(stdunit,'(2(1x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue endif endif @@ -2713,10 +2713,10 @@ logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title) if (test_ifndp) stdunit = stderr ! In case of wrong results, write to error stream write(stdunit,'(a)') title if (test_ifndp) then - write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') & + write(stdunit,'(4(1x,a,f20.16),2(1x,a,1pe22.15),1x,a)') & 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue,'WRONG!' else - write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') & + write(stdunit,'(4(1x,a,f20.16),2(1x,a,1pe22.15))') & 'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue endif endif @@ -2746,11 +2746,11 @@ logical function test_data1d(verbose, nk, Po, Ptrue, title) do k = 1,nk if (Po(k) /= Ptrue(k)) then test_data1d = .true. - write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') & + write(stdunit,'(a,i2,2(1x,a,f20.16),1x,a,1pe22.15,1x,a)') & 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k),'WRONG!' else if (verbose) & - write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') & + write(stdunit,'(a,i2,2(1x,a,f20.16),1x,a,1pe22.15)') & 'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k) endif enddo @@ -2781,10 +2781,10 @@ logical function test_data1di(verbose, nk, Po, Ptrue, title) do k = 1,nk if (Po(k) /= Ptrue(k)) then test_data1di = .true. - write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k),'WRONG!' + write(stdunit,'(a,i2,2(1x,a,i5),1x,a)') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k),'WRONG!' else if (verbose) & - write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k) + write(stdunit,'(a,i2,2(1x,a,i5))') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k) endif enddo endif diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 79fd4f3cbf..42e09b574e 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -326,6 +326,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & is, ie, js, je, k, G, GV, US, usePPM, useHuynh) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + integer, intent(in) :: ntr !< The number of tracers type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous !! tracer change [H L2 ~> m3 or kg] @@ -337,7 +338,6 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & logical, dimension(SZJ_(G),SZK_(GV)), intent(inout) :: domore_u !< If true, there is more advection to be !! done in this u-row real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1] - integer, intent(in) :: ntr !< The number of tracers integer, intent(in) :: is !< The starting tracer i-index to work on integer, intent(in) :: ie !< The ending tracer i-index to work on integer, intent(in) :: js !< The starting tracer j-index to work on @@ -698,6 +698,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & is, ie, js, je, k, G, GV, US, usePPM, useHuynh) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + integer, intent(in) :: ntr !< The number of tracers type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous !! tracer change [H L2 ~> m3 or kg] @@ -709,7 +710,6 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & logical, dimension(SZJB_(G),SZK_(GV)), intent(inout) :: domore_v !< If true, there is more advection to be !! done in this v-row real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1] - integer, intent(in) :: ntr !< The number of tracers integer, intent(in) :: is !< The starting tracer i-index to work on integer, intent(in) :: ie !< The ending tracer i-index to work on integer, intent(in) :: js !< The starting tracer j-index to work on