From b6dbc3788c753e68fe331fdbed8d11543aa08e23 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Thu, 12 Jun 2025 12:50:45 +0000 Subject: [PATCH 1/5] Add ALE_COORDINATE_CONFIG HYBRID_3D and HYBRID_MAP Add two options supporting 3D HYCOM1 TARGET_DENSITIES and ALE_RESOLUTION. The fully general interface is ALE_COORDINATE_CONFIG=HYBRID_3D which is like HYBRID except that the netCDF variables (e.g. sigma2 and dz) are 3D fields (x-y-z). Like HYBRID, the 2nd variable (e.g. dz) can be replaced by a FNC1 string that sets a dz profile that is used everywhere. The typical use of 3D HYCOM1 fields is to have different TARGET_DENSITIES in semi-enclosed seas for layers that are always below its sill depth. For example, in the OM4 75-layer setup layers 71 amd 72 are only active in the Arctic and layers 74 and 75 are only active in the Mediterranean. So the number of layers can be reduced to 71 by specifying different deep targets in these semi-emclosed seas. In this case ALE_RESOLUTION would be the same everywhere, but the resulting 71 layer configuration can be further improved by using shallow ALE_RESOLUTION in the northern hemisphere and the standard, deep, dz's in the southern hemisphere. This is safe, because deep layers are never in fixed coordinates near the equator. To simplify this typical use case, ALE_COORDINATE_CONFIG=HYBRID_MAP requires three netCDF variables (e.g. map,sigma2,dz) with the last two being 2-D z-index fields containing a small number of profiles and the first a 2-D x-y field of index values indicating which profile to use at each location. Like HYBRID and HYBRID_3D, the 3rd variable (e.g. dz) can be replaced by a FNC1 string that sets a dz profile that is used everywhere. Answers are not changes unless ALE_COORDINATE_CONFIG is set to HYBRID_3D or HYBRID_MAP. --- src/ALE/MOM_ALE.F90 | 10 +- src/ALE/MOM_regridding.F90 | 389 +++++++++++++++--- src/ALE/coord_hycom.F90 | 144 +++++-- src/core/MOM.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 4 +- src/framework/MOM_diag_remap.F90 | 5 +- .../MOM_state_initialization.F90 | 2 +- src/ocean_data_assim/MOM_oda_driver.F90 | 4 +- 8 files changed, 460 insertions(+), 100 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 93c4d35faa..fd07281039 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -165,8 +165,9 @@ module MOM_ALE !! before the main time integration loop to initialize the regridding stuff. !! We read the MOM_input file to register the values of different !! regridding/remapping parameters. -subroutine ALE_init( param_file, GV, US, max_depth, CS) +subroutine ALE_init( param_file, G, GV, US, max_depth, CS) type(param_file_type), intent(in) :: param_file !< Parameter file + type(ocean_grid_type), intent(in) :: G !< 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, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m]. @@ -205,7 +206,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) default=.false.) ! Initialize and configure regridding - call ALE_initRegridding(GV, US, max_depth, param_file, mdl, CS%regridCS) + call ALE_initRegridding(G, 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. @@ -1664,7 +1665,8 @@ end subroutine TS_PPM_edge_values !> Initializes regridding for the main ALE algorithm -subroutine ALE_initRegridding(GV, US, max_depth, param_file, mdl, regridCS) +subroutine ALE_initRegridding(G, GV, US, max_depth, param_file, mdl, regridCS) + type(ocean_grid_type), intent(in) :: G !< 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, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m]. @@ -1680,7 +1682,7 @@ subroutine ALE_initRegridding(GV, US, max_depth, param_file, mdl, regridCS) trim(regriddingCoordinateModeDoc), & default=DEFAULT_COORDINATE_MODE, fail_if_missing=.true.) - call initialize_regridding(regridCS, GV, US, max_depth, param_file, mdl, coord_mode, '', '') + call initialize_regridding(regridCS, G, GV, US, max_depth, param_file, mdl, coord_mode, '', '') end subroutine ALE_initRegridding diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 926f1f741b..92bfda4a16 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -3,9 +3,10 @@ module MOM_regridding ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_error_handler, only : MOM_error, FATAL, WARNING, assert +use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, assert use MOM_file_parser, only : param_file_type, get_param, log_param use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data +use MOM_io, only : read_variable use MOM_io, only : vardesc, var_desc, SINGLE_FILE use MOM_io, only : MOM_netCDF_file, MOM_field use MOM_io, only : create_MOM_file, MOM_write_field @@ -14,6 +15,8 @@ module MOM_regridding use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : EOS_type, calculate_density +use MOM_domains, only : max_across_PEs +!!!use MOM_domains, only : max_across_PEs, pass_var use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real use MOM_remapping, only : remapping_CS @@ -30,6 +33,7 @@ module MOM_regridding use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho use coord_rho, only : old_inflate_layers_1d use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom +use coord_hycom, only : init_3d_coord_hycom 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 @@ -66,6 +70,12 @@ module MOM_regridding !> A flag to indicate that the target_density arrays has been filled with data. logical :: target_density_set = .false. + !> Nominal HYCOM1 3D near-surface resolution [Z ~> m] + real, allocatable, dimension(:,:,:) :: coordinateResolution_3d + + !> Nominal HYCOM1 3D density of interfaces [R ~> kg m-3] + real, allocatable, dimension(:,:,:) :: target_density_3d + !> This array is set by function set_regrid_max_depths() !! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2]. real, dimension(:), allocatable :: max_interface_depths @@ -183,8 +193,10 @@ module MOM_regridding contains !> Initialization and configures a regridding control structure based on customizable run-time parameters -subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix) +subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & + coord_mode, param_prefix, param_suffix) type(regridding_CS), intent(inout) :: CS !< Regridding control structure + 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, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m]. @@ -197,6 +209,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! Local variables integer :: ke ! Number of levels + integer :: np ! Number of profiles, for HYBRID_MAP character(len=80) :: string, string2, varName ! Temporary strings character(len=40) :: coord_units, coord_res_param ! Temporary strings character(len=MAX_PARAM_LENGTH) :: param_name @@ -217,9 +230,16 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m real :: adaptBuoyCoeff, adaptAlpha ! Temporary variables for input parameters [nondim] real :: adaptZoom ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2] real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3] - integer :: k, nzf(4) + integer :: i, j, k, nzf(4) real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m] ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. + real, dimension(:,:), allocatable :: dz_2d ! 2D resolution (thickness) in units of coordinate, which may be [m] + ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. + real, dimension(:,:,:), allocatable :: dz_3d ! 3D resolution (thickness) in units of coordinate, which may be [m] + ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. + real, dimension(:,:), allocatable :: rho_target_2d ! 2D target density used in HYBRID mode [kg m-3] + real, dimension(:,:,:), allocatable :: rho_target_3d ! 3D target density used in HYBRID mode [kg m-3] + real, dimension(:,:), allocatable :: index_map ! Region array of indexes for HYBRID_MAP [nondim] real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2] real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other ! units depending on the coordinate @@ -264,8 +284,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m main_parameters=.false. if (len_trim(param_prefix)==0) main_parameters=.true. - if (main_parameters .and. len_trim(param_suffix)>0) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// & - 'Suffix provided without prefix for parameter names!') + if (main_parameters .and. len_trim(param_suffix)>0) call MOM_error(FATAL,trim(mdl)//& + ' initialize_regridding: Suffix provided without prefix for parameter names!') CS%nk = 0 CS%regridding_scheme = coordinateMode(coord_mode) @@ -356,8 +376,20 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//& " HYBRID:string - read from a file. The string specifies\n"//& " the filename and two variable names, separated\n"//& - " by a comma or space, for sigma-2 and dz. e.g.\n"//& - " HYBRID:vgrid.nc,sigma2,dz",& + " by a comma or space, for sigma-2 and dz.\n"//& + " e.g. HYBRID:vgrid.nc,sigma2,dz\n"//& + " HYBRID_3D:string - read from a file. The string specifies\n"//& + " the filename and two 3D variable names, separated\n"//& + " by a comma or space, for sigma-2 and dz. The\n"//& + " latter can be FNC1:string which is used everywhere.\n"//& + " e.g. HYBRID_3D:vgrid.nc,sigma2,dz\n"//& + " HYBRID_MAP:string - read from a file. The string specifies\n"//& + " the filename and three variable names, separated\n"//& + " by a comma or space, for map, sigma-2 and dz.\n"//& + " Map is a spatial index array with, maxval(map)=N,\n"//& + " and the others are 2D arrays containing N profiles.\n"//& + " Dz can be FNC1:string which is used everywhere.\n"//& + " e.g. HYBRID_MAP:vgrid.nc,map,sigma2,dz",& default=trim(string2)) message = "The distribution of vertical resolution for the target\n"//& "grid used for Eulerian-like coordinates. For example,\n"//& @@ -424,11 +456,13 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m varName=trim(varName(12:)) call verify_variable_units(filename, varName, expected_units, message, ierr, alt_units) if (ierr) call MOM_error(FATAL, trim(mdl)//", initialize_regridding: "//& - "Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message)) + "Unsupported format in grid definition '"//trim(filename)//& + "'. Error message "//trim(message)) call field_size(trim(fileName), trim(varName), nzf) ke = nzf(1)-1 - if (ke < 1) call MOM_error(FATAL, trim(mdl)//" initialize_regridding via Var "//trim(varName)//& - "in FILE "//trim(filename)//" requires at least 2 target interface values.") + if (ke < 1) call MOM_error(FATAL, trim(mdl)//" initialize_regridding via Var "//& + trim(varName)//"in FILE "//trim(filename)//& + " requires at least 2 target interface values.") if (CS%regridding_scheme == REGRIDDING_RHO) then allocate(rho_target(ke+1)) call MOM_read_data(trim(fileName), trim(varName), rho_target) @@ -461,22 +495,30 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! Function used for set target interface densities ke = rho_function1( trim(string(7:)), rho_target ) elseif (index(trim(string),'HYBRID:')==1) then - ke = GV%ke; allocate(dz(ke)) - ! The following assumes the FILE: syntax of above but without "FILE:" in the string + ke = GV%ke + allocate(dz(ke)) allocate(rho_target(ke+1)) + ! The following assumes the FILE: syntax of above but without "FILE:" in the string + varName = trim( extractWord(trim(string(8:)), 3) ) + if (varname == " ") call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID "// & + "Too few arguments in ("//trim(string)//")") fileName = trim( extractWord(trim(string(8:)), 1) ) if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) - if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// & + if (.not. file_exists(fileName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID "// & "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") varName = trim( extractWord(trim(string(8:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// & + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID "// & "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), rho_target) varName = trim( extractWord(trim(string(8:)), 3) ) if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) else ! Read dz from file - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// & + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID "// & "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), dz) endif @@ -486,6 +528,169 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, & 'HYBRID target densities for interfaces', units=coordinateUnits(coord_mode)) endif + elseif (index(trim(string),'HYBRID_3D:')==1) then + ke = GV%ke + allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0) + allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0) + ! The following assumes the FILE: syntax of above but without "FILE:" in the string + varName = trim( extractWord(trim(string(11:)), 3) ) + if (varname == " ") call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Too few arguments in ("//trim(string)//")") + fileName = trim( extractWord(trim(string(11:)), 1) ) + if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) + if (.not. file_exists(fileName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + varName = trim( extractWord(trim(string(11:)), 2) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + call MOM_read_data(trim(fileName), trim(varName), rho_target_3d, G%Domain) +!!!! call pass_var(rho_target_3d, G%Domain, halo=1) + varName = trim( extractWord(trim(string(11:)), 3) ) + if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d + allocate(dz(ke)) + call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) + ! Adjust target grid to be consistent with maximum_depth + tmpReal = sum( dz(:) ) + if (tmpReal < maximum_depth) then + dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) + endif + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + do k=1,ke + dz_3d(i,j,k) = dz(k) + enddo + endif !mask2dT + enddo; enddo + if (main_parameters) then + call log_param(param_file, mdl, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + endif + else ! Read dz from file + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + call MOM_read_data(trim(fileName), trim(varName), dz_3d, G%Domain) +!!!! call pass_var(dz_3d, G%Domain, halo=1) + ! set nominal 1-d dz to UNIFORM + allocate(dz(ke)) + dz(:) = uniformResolution(ke, coord_mode, maximum_depth, & + US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), & + US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) ) + endif !dz + elseif (index(trim(string),'HYBRID_MAP:')==1) then + ke = GV%ke + allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0) + allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0) + allocate(index_map(SZI_(G),SZJ_(G)), source=1.0) + ! The following assumes the FILE: syntax of above but without "FILE:" in the string + varName = trim( extractWord(trim(string(12:)), 4) ) + if (varname == " ") call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Too few arguments in ("//trim(string)//")") + fileName = trim( extractWord(trim(string(12:)), 1) ) + if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) + if (.not. file_exists(fileName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + varName = trim( extractWord(trim(string(12:)), 2) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + call MOM_read_data(trim(fileName), trim(varName), index_map, G%Domain) +!!!! call pass_var(index_map, G%Domain, halo=1) + !find maximum index + tmpReal = 1 + do j=G%jsc, G%jec ; do i=G%isc, G%iec + tmpReal = max(tmpReal,index_map(i,j)) + enddo ; enddo + np = tmpReal + call max_across_PEs(np) + write(string2,"(i3)") np + call MOM_error(NOTE, & + trim(mdl)//", initialize_regridding: HYBRID_MAP NP="//trim(string2)) + if (np<1) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_MAP to small NP from "//trim(varName)) + allocate(dz_2d(ke,np)) + allocate(rho_target_2d(ke+1,np)) + varName = trim( extractWord(trim(string(12:)), 3) ) + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + ! MOM_read_data can't handle this array + call read_variable(trim(fileName), trim(varName), rho_target_2d) + if (main_parameters) then + call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target_2d(:,1), & + 'HYBRID target densities for interfaces', units=coordinateUnits(coord_mode)) + endif + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + if (nint(index_map(i,j))<1 .or. nint(index_map(i,j))>np) then + write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=', index_map(i,j) + call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "index_map out of range") + endif + rho_target_3d(i,j,:) = rho_target_2d(:,nint(index_map(i,j))) + endif !mask2dT + enddo; enddo + varName = trim( extractWord(trim(string(12:)), 4) ) + if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d + allocate(dz(ke)) + call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) + ! Adjust target grid to be consistent with maximum_depth + tmpReal = sum( dz(:) ) + if (tmpReal < maximum_depth) then + dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) + endif + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + do k=1,ke + dz_3d(i,j,k) = dz(k) + enddo + endif !mask2dT + enddo; enddo + if (main_parameters) then + call log_param(param_file, mdl, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + endif + else ! Read dz from file + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + ! MOM_read_data can't handle this array + call read_variable(trim(fileName), trim(varName), dz_2d) + if (main_parameters) then + call log_param(param_file, mdl, "!"//coord_res_param, dz_2d(:,1), & + trim(message), units=coordinateUnits(coord_mode)) + endif + do i=1,np + tmpReal = sum( dz_2d(:,i) ) + if (tmpReal < maximum_depth) then + dz_2d(ke,i) = dz_2d(ke,i) + ( maximum_depth - tmpReal ) + endif + enddo + allocate(dz(ke)) + dz(:) = dz_2d(:,1) + if (main_parameters) then + call log_param(param_file, mdl, "!"//coord_res_param, dz, & + trim(message), units=coordinateUnits(coord_mode)) + endif + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + if (nint(index_map(i,j))<1 .or. nint(index_map(i,j))>np) then + write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=',index_map(i,j) + call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "index_map out of range") + endif + dz_3d(i,j,:) = dz_2d(:,nint(index_map(i,j))) + endif !mask2dT + enddo; enddo + endif !dz + deallocate(index_map) + deallocate(rho_target_2d) + deallocate(dz_2d) elseif (index(trim(string),'WOA09INT')==1) then if (len_trim(string)==8) then ! string=='WOA09INT' tmpReal = 0. ; ke = 0 ; dz_extra = 0. @@ -566,31 +771,45 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m coordinateMode(coord_mode) == REGRIDDING_HYCOM1 .or. & coordinateMode(coord_mode) == REGRIDDING_HYBGEN .or. & coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then - ! Adjust target grid to be consistent with maximum_depth - tmpReal = sum( dz(:) ) - if (tmpReal < maximum_depth) then - dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) - elseif (tmpReal > maximum_depth) then - if ( dz(ke) + ( maximum_depth - tmpReal ) > 0. ) then + if (allocated(dz)) then + ! Adjust target grid to be consistent with maximum_depth + tmpReal = sum( dz(:) ) + if (tmpReal < maximum_depth) then dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) - else - call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) + elseif (tmpReal > maximum_depth) then + if ( dz(ke) + ( maximum_depth - tmpReal ) > 0. ) then + dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) + else + call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & + "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) + endif endif - endif + endif !allocated(dz) endif endif CS%nk=ke ! Target resolution (for fixed coordinates) - allocate( CS%coordinateResolution(CS%nk), source=-1.E30 ) - if (state_dependent(CS%regridding_scheme)) then - ! Target values - allocate( CS%target_density(CS%nk+1), source=-1.E30*US%kg_m3_to_R ) + if (allocated(dz_3d)) then + allocate( CS%coordinateResolution(CS%nk), source=-1.E30 ) + allocate( CS%coordinateResolution_3d(SZI_(G),SZJ_(G),CS%nk), source=-1.E30 ) + allocate( CS%target_density_3d(SZI_(G),SZJ_(G),CS%nk+1), source=-1.E30*US%kg_m3_to_R ) + else + allocate( CS%coordinateResolution(CS%nk), source=-1.E30 ) + if (state_dependent(CS%regridding_scheme)) then + ! Target values + allocate( CS%target_density(CS%nk+1), source=-1.E30*US%kg_m3_to_R ) + endif endif - if (allocated(dz)) then + if (allocated(dz_3d)) then + ! set both 1d and 3d fields + call setCoordinateResolution(dz, CS, scale=US%m_to_Z) + call setCoordinateResolution_3d(dz_3d, CS, scale=US%m_to_Z) + CS%coord_scale = US%Z_to_m + deallocate(dz_3d) + elseif (allocated(dz)) then if (coordinateMode(coord_mode) == REGRIDDING_SIGMA) then call setCoordinateResolution(dz, CS, scale=1.0) elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then @@ -612,11 +831,12 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! ensure CS%ref_pressure is rescaled properly CS%ref_pressure = US%Pa_to_RL2_T2 * CS%ref_pressure - if (allocated(rho_target)) then + if (allocated(rho_target_3d)) then + call set_target_densities_3d(CS, G, US%kg_m3_to_R, rho_target_3d) + deallocate(rho_target_3d) + elseif (allocated(rho_target)) then call set_target_densities(CS, US%kg_m3_to_R*rho_target) deallocate(rho_target) - - ! \todo This line looks like it would overwrite the target densities set just above? elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then call set_target_densities_from_GV(GV, US, CS) call log_param(param_file, mdl, "!TARGET_DENSITIES", US%R_to_kg_m3*CS%target_density(:), & @@ -624,7 +844,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, param_file) + call initCoord(CS, G, GV, US, coord_mode, param_file) if (coord_is_state_dependent) then if (main_parameters) then @@ -634,17 +854,18 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", & units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2) else - call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), P_Ref, & - "The pressure that is used for calculating the diagnostic coordinate "//& - "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//& - "This is only used for the RHO coordinate.", & - units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2) + call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), & + P_Ref, & + "The pressure that is used for calculating the diagnostic coordinate "//& + "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//& + "This is only used for the RHO coordinate.", & + units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2) endif - call get_param(param_file, mdl, create_coord_param(param_prefix, "REGRID_COMPRESSIBILITY_FRACTION", param_suffix), & - tmpReal, & - "When interpolating potential density profiles we can add "//& - "some artificial compressibility solely to make homogeneous "//& - "regions appear stratified.", units="nondim", default=0.) + call get_param(param_file, mdl, create_coord_param(param_prefix, "REGRID_COMPRESSIBILITY_FRACTION", & + param_suffix), tmpReal, & + "When interpolating potential density profiles we can add "//& + "some artificial compressibility solely to make homogeneous "//& + "regions appear stratified.", units="nondim", default=0.) call set_regrid_params(CS, compress_fraction=tmpReal, ref_pressure=P_Ref) endif @@ -660,8 +881,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m if (main_parameters .and. coordinateMode(coord_mode) == REGRIDDING_HYCOM1) then call get_param(param_file, mdl, "HYCOM1_ONLY_IMPROVES", tmpLogical, & - "When regridding, an interface is only moved if this improves the fit to the target density.", & - default=.false.) + "When regridding, an interface is only moved if this improves "//& + "the fit to the target density.", default=.false.) call set_hycom_params(CS%hycom_CS, only_improves=tmpLogical) endif @@ -835,6 +1056,8 @@ subroutine end_regridding(CS) if (associated(CS%hybgen_CS)) call end_hybgen_regrid(CS%hybgen_CS) deallocate( CS%coordinateResolution ) + if (allocated(CS%coordinateResolution_3d)) deallocate( CS%coordinateResolution_3d ) + if (allocated(CS%target_density_3d)) deallocate( CS%target_density_3d ) if (allocated(CS%target_density)) deallocate( CS%target_density ) if (allocated(CS%max_interface_depths) ) deallocate( CS%max_interface_depths ) if (allocated(CS%max_layer_thickness) ) deallocate( CS%max_layer_thickness ) @@ -1448,10 +1671,12 @@ subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface ) write(0,*) k,zOld(nz+1),zNew(k) enddo do k=1,min(nz,CS%nk) - write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k) + write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k), & + CS%coordinateResolution(k) enddo do k=min(nz,CS%nk)+1,CS%nk - write(0,*) k,0.0,zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k),CS%coordinateResolution(k) + write(0,*) k,0.0,zNew(k)-zNew(k+1),totalThickness*CS%coordinateResolution(k), & + CS%coordinateResolution(k) enddo call MOM_error( FATAL, & 'MOM_regridding, build_sigma_grid: top surface has moved!!!' ) @@ -1701,7 +1926,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) enddo - call build_hycom1_column(CS%hycom_CS, remapCS, tv%eqn_of_state, GV%ke, nominalDepth, & + call build_hycom1_column(CS%hycom_CS, remapCS, tv%eqn_of_state, GV%ke, i, j, nominalDepth, & h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & z_col, z_col_new, zScale=zScale, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) @@ -1783,7 +2008,8 @@ subroutine build_grid_adaptive(G, GV, US, h, nom_depth_H, tv, dzInterface, remap cycle endif - call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, nom_depth_H, zNext) + call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, & + nom_depth_H, zNext) call filtered_grid_motion(CS, nz, zInt(i,j,:), zNext, dzInterface(i,j,:)) ! convert from depth to z @@ -1950,7 +2176,8 @@ subroutine convective_adjustment(G, GV, h, tv) tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0 h(i,j,k) = h1 ; h(i,j,k+1) = h0 ! Recompute densities at levels k and k+1 - 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), 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 @@ -2014,13 +2241,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, param_file) +subroutine initCoord(CS, G, GV, US, coord_mode, param_file) type(regridding_CS), intent(inout) :: CS !< Regridding control structure + 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 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)) @@ -2033,8 +2261,14 @@ subroutine initCoord(CS, GV, US, coord_mode, param_file) case (REGRIDDING_RHO) call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS) case (REGRIDDING_HYCOM1) - call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, & - CS%interp_CS) + if (allocated(CS%coordinateResolution_3d)) then + call init_3d_coord_hycom(CS%hycom_CS, G, CS%nk, & + CS%coordinateResolution_3d, CS%target_density_3d, & + CS%interp_CS) + else + call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, & + CS%interp_CS) + endif case (REGRIDDING_HYBGEN) call init_hybgen_regrid(CS%hybgen_CS, GV, US, param_file) case (REGRIDDING_ADAPTIVE) @@ -2065,6 +2299,26 @@ subroutine setCoordinateResolution( dz, CS, scale ) end subroutine setCoordinateResolution +!> Set the 3d fixed resolution data +subroutine setCoordinateResolution_3d( dz_3d, CS, scale ) + real, dimension(:,:,:), intent(in) :: dz_3d !< A vector of vertical grid spacings, in arbitrary coordinate + !! dependent units, such as [m] for a z-coordinate or [kg m-3] + !! for a density coordinate. + type(regridding_CS), intent(inout) :: CS !< Regridding control structure + real, optional, intent(in) :: scale !< A scaling factor converting dz to coordRes [m -> Z] + + if (.not.allocated(CS%coordinateResolution_3d)) & + call MOM_error(FATAL,'setCoordinateResolution_3d: '//& + 'CS%coordinateResolution_3d not allocated.') + + if (present(scale)) then + CS%coordinateResolution_3d(:,:,:) = scale*dz_3d(:,:,:) + else + CS%coordinateResolution_3d(:,:,:) = dz_3d(:,:,:) + endif + +end subroutine setCoordinateResolution_3d + !> Set target densities based on the old Rlay variable subroutine set_target_densities_from_GV( GV, US, CS ) type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -2088,6 +2342,22 @@ subroutine set_target_densities_from_GV( GV, US, CS ) end subroutine set_target_densities_from_GV +!> Set target densities based on vector of interface values +subroutine set_target_densities_3d( CS, G, scale, rho_int_3d ) + type(regridding_CS), intent(inout) :: CS !< Regridding control structure + type(ocean_grid_type),intent(in) :: G !< Ocean grid structure + real, intent(in) :: scale !< A scaling factor converting densities [kg m-3 -> R] + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: rho_int_3d !< Interface densities [kg m-3] + + if (.not.allocated(CS%target_density_3d)) & + call MOM_error(FATAL,'set_target_densities_3d: '//& + 'CS%target_density_3d not allocated.') + + CS%target_density_3d(:,:,:) = scale * rho_int_3d(:,:,:) + CS%target_density_set = .true. + +end subroutine set_target_densities_3d + !> Set target densities based on vector of interface values subroutine set_target_densities( CS, rho_int ) type(regridding_CS), intent(inout) :: CS !< Regridding control structure @@ -2376,7 +2646,8 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & compress_fraction, ref_pressure, & integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, & - adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0) + adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, & + adaptAlpha, adaptDoMin, adaptDrho0) type(regridding_CS), intent(inout) :: CS !< Regridding control structure logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the @@ -2511,8 +2782,8 @@ function getStaticThickness( CS, SSH, depth ) real :: z, dz ! Vertical positions and grid spacing [Z ~> m] select case ( CS%regridding_scheme ) - case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, & - REGRIDDING_ADAPTIVE ) + case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, & + REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE ) if (depth>0.) then z = ssh do k = 1, CS%nk diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index 1e5474770a..e2a2f5aa3a 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -3,14 +3,17 @@ module coord_hycom ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_error_handler, only : MOM_error, FATAL -use MOM_remapping, only : remapping_CS, remapping_core_h +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, NOTE +use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_EOS, only : EOS_type, calculate_density +use MOM_remapping, only : remapping_CS, remapping_core_h use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, regridding_set_ppolys use regrid_interp, only : DEGREE_MAX implicit none ; private +#include + !> Control structure containing required parameters for the HyCOM coordinate type, public :: hycom_CS ; private @@ -32,11 +35,20 @@ module coord_hycom !> If true, an interface only moves if it improves the density fit logical :: only_improves = .false. + !> If true, use 3-D control fields + logical :: use_3d = .false. + + !> Nominal density of interfaces [R ~> kg m-3] + real, allocatable, dimension(:,:,:) :: target_density_3d + + !> Nominal near-surface resolution [Z ~> m] + real, allocatable, dimension(:,:,:) :: coordinateResolution_3d + !> Interpolation control structure type(interp_CS_type) :: interp_CS end type hycom_CS -public init_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom +public init_coord_hycom, init_3d_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom contains @@ -56,18 +68,59 @@ subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp CS%nk = nk CS%coordinateResolution(:) = coordinateResolution(:) CS%target_density(:) = target_density(:) + CS%use_3d = .false. CS%interp_CS = interp_CS + if (is_root_pe()) call MOM_error(NOTE, "init_coord_hycom: use_3d = .false.") + end subroutine init_coord_hycom +!> Initialise a hycom_CS with pointers to parameters +subroutine init_3d_coord_hycom(CS, G, nk, coordinateResolution, target_density, interp_CS) + type(hycom_CS), pointer :: CS !< Unassociated pointer to hold the control structure + type(ocean_grid_type),intent(in) :: G !< Ocean grid structure + integer, intent(in) :: nk !< Number of layers in generated grid + real, dimension(SZI_(G),SZJ_(G),nk), intent(in) :: coordinateResolution !< Nominal near-surface resolution [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),nk+1), intent(in) :: target_density !< Interface target densities [R ~> kg m-3] + type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation + ! Local variables + integer :: i,j,k + + if (associated(CS)) call MOM_error(FATAL, "init_3d_coord_hycom: CS already associated!") + + allocate(CS) + allocate(CS%coordinateResolution_3d(nk,SZI_(G),SZJ_(G)), source=0.0) + allocate(CS%target_density_3d(nk+1,SZI_(G),SZJ_(G)), source=0.0) + + CS%nk = nk + CS%use_3d = .true. + CS%interp_CS = interp_CS + + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + do k= 1,nk + CS%coordinateResolution_3d(k,i,j) = coordinateResolution(i,j,k) + CS%target_density_3d(k,i,j) = target_density(i,j,k) + enddo + CS%target_density_3d(nk+1,i,j) = target_density(i,j,nk+1) + endif !mask2dT + enddo; enddo + + if (is_root_pe()) call MOM_error(NOTE, "init_3d_coord_hycom: use_3d = .true.") + +end subroutine init_3d_coord_hycom + !> This subroutine deallocates memory in the control structure for the coord_hycom module subroutine end_coord_hycom(CS) type(hycom_CS), pointer :: CS !< Coordinate control structure ! nothing to do if (.not. associated(CS)) return - deallocate(CS%coordinateResolution) - deallocate(CS%target_density) + + if (allocated(CS%coordinateResolution)) deallocate(CS%coordinateResolution) + if (allocated(CS%target_density)) deallocate(CS%target_density) + if (allocated(CS%coordinateResolution_3d)) deallocate(CS%coordinateResolution_3d) + if (allocated(CS%target_density_3d)) deallocate(CS%target_density_3d) if (allocated(CS%max_interface_depths)) deallocate(CS%max_interface_depths) if (allocated(CS%max_layer_thickness)) deallocate(CS%max_layer_thickness) deallocate(CS) @@ -103,12 +156,14 @@ subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, only_ end subroutine set_hycom_params !> Build a HyCOM coordinate column -subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_col, & +subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, p_col, & z_col, z_col_new, zScale, h_neglect, h_neglect_edge) type(hycom_CS), intent(in) :: CS !< Coordinate control structure type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure integer, intent(in) :: nz !< Number of levels + integer, intent(in) :: ix !< x direction array index + integer, intent(in) :: jy !< y direction array index real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) real, dimension(nz), intent(in) :: T !< Temperature of column [C ~> degC] real, dimension(nz), intent(in) :: S !< Salinity of column [S ~> ppt] @@ -150,7 +205,7 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ z_scale = 1.0 ; if (present(zScale)) z_scale = zScale if (CS%only_improves .and. nz == CS%nk) then - call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, depth, & + call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, ix, jy, depth, & h, T, S, p_col, rho_col, RiA_ini, h_neglect, h_neglect_edge) else ! Work bottom recording potential density @@ -164,8 +219,13 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ ! Interpolates for the target interface position with the rho_col profile ! Based on global density profile, interpolate to generate a new grid - call build_and_interpolate_grid(CS%interp_CS, rho_col, nz, h(:), z_col, & - CS%target_density, CS%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge) + if (CS%use_3d) then + call build_and_interpolate_grid(CS%interp_CS, rho_col, nz, h(:), z_col, & + CS%target_density_3d(:,ix,jy), CS%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge) + else + call build_and_interpolate_grid(CS%interp_CS, rho_col, nz, h(:), z_col, & + CS%target_density, CS%nk, h_col_new, z_col_new, h_neglect, h_neglect_edge) + endif if (CS%only_improves .and. nz == CS%nk) then ! Only move an interface if it improves the density fit z_1 = 0.5 * ( z_col(1) + z_col(2) ) @@ -177,7 +237,7 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ ! Remap from original h and T,S to get T,S_col_new call remapping_core_h(remapCS, nz, h(:), T, CS%nk, h_col_new, T_col_new) call remapping_core_h(remapCS, nz, h(:), S, CS%nk, h_col_new, S_col_new) - call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, depth, & + call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, ix, jy, depth, & h_col_new, T_col_new, S_col_new, p_col_new, r_col_new, RiA_new, h_neglect, h_neglect_edge) do k= 2,CS%nk if ( abs(RiA_ini(K)) <= abs(RiA_new(K)) .and. z_col(K) > z_col_new(K-1) .and. & @@ -191,11 +251,19 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ ! as deep as a nominal target z* grid nominal_z = 0. stretching = z_col(nz+1) / depth ! Stretches z* to z - do k = 2, CS%nk+1 - nominal_z = nominal_z + (z_scale * CS%coordinateResolution(k-1)) * stretching - z_col_new(k) = max( z_col_new(k), nominal_z ) - z_col_new(k) = min( z_col_new(k), z_col(nz+1) ) - enddo + if (CS%use_3d) then + do k = 2, CS%nk+1 + nominal_z = nominal_z + (z_scale * CS%coordinateResolution_3d(k-1,ix,jy)) * stretching + z_col_new(k) = max( z_col_new(k), nominal_z ) + z_col_new(k) = min( z_col_new(k), z_col(nz+1) ) + enddo + else + do k = 2, CS%nk+1 + nominal_z = nominal_z + (z_scale * CS%coordinateResolution(k-1)) * stretching + z_col_new(k) = max( z_col_new(k), nominal_z ) + z_col_new(k) = min( z_col_new(k), z_col(nz+1) ) + enddo + endif if (maximum_depths_set .and. maximum_h_set) then ; do k=2,CS%nk ! The loop bounds are 2 & nz so the top and bottom interfaces do not move. @@ -210,12 +278,14 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ end subroutine build_hycom1_column !> Calculate interface density anomaly w.r.t. the target. -subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_col, & +subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, p_col, & R, RiAnom, h_neglect, h_neglect_edge) type(hycom_CS), intent(in) :: CS !< Coordinate control structure type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure integer, intent(in) :: nz !< Number of levels + integer, intent(in) :: ix !< x direction array index + integer, intent(in) :: jy !< y direction array index real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) real, dimension(nz), intent(in) :: T !< Temperature of column [C ~> degC] real, dimension(nz), intent(in) :: S !< Salinity of column [S ~> ppt] @@ -247,19 +317,35 @@ subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, depth, h, call regridding_set_ppolys(CS%interp_CS, rho_col, nz, h, ppoly_E, ppoly_S, ppoly_C, & degree, h_neglect, h_neglect_edge) - R(1) = rho_col(1) - RiAnom(1) = ppoly_E(1,1) - CS%target_density(1) - do k= 2,nz - R(k) = rho_col(k) - if (ppoly_E(k-1,2) > CS%target_density(k)) then - RiAnom(k) = ppoly_E(k-1,2) - CS%target_density(k) !interface is heavier than target - elseif (ppoly_E(k,1) < CS%target_density(k)) then - RiAnom(k) = ppoly_E(k,1) - CS%target_density(k) !interface is lighter than target - else - RiAnom(k) = 0.0 !interface spans the target - endif - enddo - RiAnom(nz+1) = ppoly_E(nz,2) - CS%target_density(nz+1) + if (CS%use_3d) then + R(1) = rho_col(1) + RiAnom(1) = ppoly_E(1,1) - CS%target_density_3d(1,ix,jy) + do k= 2,nz + R(k) = rho_col(k) + if (ppoly_E(k-1,2) > CS%target_density_3d(k,ix,jy)) then + RiAnom(k) = ppoly_E(k-1,2) - CS%target_density_3d(k,ix,jy) !interface is heavier than target + elseif (ppoly_E(k,1) < CS%target_density_3d(k,ix,jy)) then + RiAnom(k) = ppoly_E(k,1) - CS%target_density_3d(k,ix,jy) !interface is lighter than target + else + RiAnom(k) = 0.0 !interface spans the target + endif + enddo + RiAnom(nz+1) = ppoly_E(nz,2) - CS%target_density_3d(nz+1,ix,jy) + else + R(1) = rho_col(1) + RiAnom(1) = ppoly_E(1,1) - CS%target_density(1) + do k= 2,nz + R(k) = rho_col(k) + if (ppoly_E(k-1,2) > CS%target_density(k)) then + RiAnom(k) = ppoly_E(k-1,2) - CS%target_density(k) !interface is heavier than target + elseif (ppoly_E(k,1) < CS%target_density(k)) then + RiAnom(k) = ppoly_E(k,1) - CS%target_density(k) !interface is lighter than target + else + RiAnom(k) = 0.0 !interface spans the target + endif + enddo + RiAnom(nz+1) = ppoly_E(nz,2) - CS%target_density(nz+1) + endif !use_3d:else end subroutine build_hycom1_target_anomaly diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fc874c8cbd..0ec56ed74c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3139,7 +3139,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)") if (CS%use_ALE_algorithm) then - call ALE_init(param_file, GV, US, G%max_depth, CS%ALE_CSp) + call ALE_init(param_file, G, GV, US, G%max_depth, CS%ALE_CSp) call callTree_waypoint("returned from ALE_init() (initialize_MOM)") endif diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 3064f633cc..58b7d39a4c 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -506,7 +506,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) do i=1, diag_cs%num_diag_coords ! For each possible diagnostic coordinate - call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, US, param_file) + call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), G, GV, US, param_file) ! Allocate these arrays since the size of the diagnostic array is now known allocate(diag_cs%diag_remap_cs(i)%h(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz)) @@ -704,7 +704,7 @@ subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_n do i=1, diag_cs%num_diag_coords ! For each possible diagnostic coordinate - !call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, param_file) + !call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), G, GV, param_file) ! This vertical coordinate has been configured so can be used. if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 38553a4351..8fa523924c 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -176,8 +176,9 @@ end subroutine diag_remap_set_active !> Configure the vertical axes for a diagnostic remapping control structure. !! Reads a configuration parameters to determine coordinate generation. -subroutine diag_remap_configure_axes(remap_cs, GV, US, param_file) +subroutine diag_remap_configure_axes(remap_cs, G, GV, US, param_file) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remap control structure + type(ocean_grid_type), intent(in) :: G !< The ocean's grid type 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 structure @@ -192,7 +193,7 @@ subroutine diag_remap_configure_axes(remap_cs, GV, US, param_file) layers ! Numerical values for layer vertical coordinates, in unscaled units ! that might be [m], [kg m-3] or [nondim], depending on the coordinate. - call initialize_regridding(remap_cs%regrid_cs, GV, US, GV%max_depth, param_file, mod, & + call initialize_regridding(remap_cs%regrid_cs, G, GV, US, GV%max_depth, param_file, mod, & trim(remap_cs%vertical_coord_name), "DIAG_COORD", trim(remap_cs%diag_coord_name)) call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index bcefb0022e..d0b0a3cd58 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2899,7 +2899,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Build the target grid (and set the model thickness to it) - call ALE_initRegridding( GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS + call ALE_initRegridding( G, GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS if (remap_general) then dz_neglect = set_h_neglect(GV, remap_answer_date, dz_neglect_edge) else diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 882296be99..3b277578a3 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -312,7 +312,7 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) call set_grid_metrics(dG, PF, CS%US) call MOM_initialize_topography(dG%bathyT, dG%max_depth, dG, PF, CS%US) call MOM_initialize_coord(CS%GV, CS%US, PF, tv_dummy, dG%max_depth) - call ALE_init(PF, CS%GV, CS%US, dG%max_depth, CS%ALE_CS) + call ALE_init(PF, CS%G, CS%GV, CS%US, dG%max_depth, CS%ALE_CS) call MOM_grid_init(CS%Grid, PF, global_indexing=.false.) call ALE_updateVerticalGridType(CS%ALE_CS, CS%GV) call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) @@ -337,7 +337,7 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) "If true, use the OM4 remapping-via-subcells algorithm for ODA. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& "We recommend setting this option to false.", default=om4_remap_via_sub_cells) - call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') + call initialize_regridding(CS%regridCS, CS%G, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') h_neglect = set_h_neglect(GV, CS%answer_date, h_neglect_edge) call initialize_remapping(CS%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells, & From 9f77d5555c4dc00906d30a5b2840d621734d40c8 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Sun, 15 Jun 2025 10:40:09 +0000 Subject: [PATCH 2/5] Set the halos of new 3D arrays --- src/ALE/MOM_regridding.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 92bfda4a16..b388dfe495 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -15,8 +15,7 @@ module MOM_regridding use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : EOS_type, calculate_density -use MOM_domains, only : max_across_PEs -!!!use MOM_domains, only : max_across_PEs, pass_var +use MOM_domains, only : max_across_PEs, pass_var use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real use MOM_remapping, only : remapping_CS @@ -547,7 +546,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & trim(mdl)//", initialize_regridding: HYBRID_3D "// & "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), rho_target_3d, G%Domain) -!!!! call pass_var(rho_target_3d, G%Domain, halo=1) + call pass_var(rho_target_3d, G%Domain, halo=1) varName = trim( extractWord(trim(string(11:)), 3) ) if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d allocate(dz(ke)) @@ -573,7 +572,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & trim(mdl)//", initialize_regridding: HYBRID_3D "// & "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), dz_3d, G%Domain) -!!!! call pass_var(dz_3d, G%Domain, halo=1) + call pass_var(dz_3d, G%Domain, halo=1) ! set nominal 1-d dz to UNIFORM allocate(dz(ke)) dz(:) = uniformResolution(ke, coord_mode, maximum_depth, & @@ -600,7 +599,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & trim(mdl)//", initialize_regridding: HYBRID_MAP "// & "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), index_map, G%Domain) -!!!! call pass_var(index_map, G%Domain, halo=1) + call pass_var(index_map, G%Domain, halo=1) !find maximum index tmpReal = 1 do j=G%jsc, G%jec ; do i=G%isc, G%iec From 2b87ecf8e64d7d8384b22c5efc6df62a7bfd74a9 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Tue, 17 Jun 2025 17:14:57 +0000 Subject: [PATCH 3/5] Add SHALLOW_ALE_RESOLUTION SHALLOW_ALE_RESOLUTION implements a HYBGEN-style Z-sigma-Z near surface fixed coordinate for HYCOM1. For example the US Navy's GOFS 3.1 HYCOM setup has 41 layers, with the top 14 layers in a Z-sigma-Z configuration. For MOM6 HYCOM1 this is: SHALLOW_ALE_RESOLUTION = 14*1.0,27*0.0 for 14 1m "shallow" layers. Let N_SIGMA be the number of consecutive non-zero entries, typically < NK. When rest depth is shallower than SUM(SHALLOW_ALE_RESOLUTION(1:N_SIGMA)) use SHALLOW_ALE_RESOLUTION. When rest depth is deeper than SUM(SHALLOW_ALE_RESOLUTION(1:N_SIGMA)) use ALE_RESOLUTION. Otherwise use a linear sum of the two weighted by rest depth. The default of all zeros turns this option off, and when off answers are unchanged. The new parmeter SHALLOW_ALE_RESOLUTION is only present when using HYCOM1. --- src/ALE/MOM_regridding.F90 | 87 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 83 insertions(+), 4 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index b388dfe495..15d1d316f0 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -209,6 +209,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! Local variables integer :: ke ! Number of levels integer :: np ! Number of profiles, for HYBRID_MAP + integer :: n_sigma ! Number of shallow dz's, for HYBRID_MAP or HYBRID_3D character(len=80) :: string, string2, varName ! Temporary strings character(len=40) :: coord_units, coord_res_param ! Temporary strings character(len=MAX_PARAM_LENGTH) :: param_name @@ -225,6 +226,10 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z). real :: dz_extra ! The thickness of an added layer to append to the woa09_dz profile when ! maximum_depth is large [m] (not in Z). + real :: nominalDepth ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2] + real :: depth_q ! A depth scale factor [nondim] + real :: depth_s ! The end of the shallow Z regime (m) + real :: depth_d ! The start of the deep Z regime (m) real :: adaptTimeRatio, adaptZoomCoeff ! Temporary variables for input parameters [nondim] real :: adaptBuoyCoeff, adaptAlpha ! Temporary variables for input parameters [nondim] real :: adaptZoom ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2] @@ -236,6 +241,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. real, dimension(:,:,:), allocatable :: dz_3d ! 3D resolution (thickness) in units of coordinate, which may be [m] ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. + real, dimension(:), allocatable :: dz_shallow ! Shallow resolution (thickness), for HYBRID_MAP or HYBRID_3D [m] real, dimension(:,:), allocatable :: rho_target_2d ! 2D target density used in HYBRID mode [kg m-3] real, dimension(:,:,:), allocatable :: rho_target_3d ! 3D target density used in HYBRID mode [kg m-3] real, dimension(:,:), allocatable :: index_map ! Region array of indexes for HYBRID_MAP [nondim] @@ -525,7 +531,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call log_param(param_file, mdl, "!"//coord_res_param, dz, & trim(message), units=coordinateUnits(coord_mode)) call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, & - 'HYBRID target densities for interfaces', units=coordinateUnits(coord_mode)) + 'HYBRID target densities for interfaces', units="kg m-3") endif elseif (index(trim(string),'HYBRID_3D:')==1) then ke = GV%ke @@ -622,7 +628,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call read_variable(trim(fileName), trim(varName), rho_target_2d) if (main_parameters) then call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target_2d(:,1), & - 'HYBRID target densities for interfaces', units=coordinateUnits(coord_mode)) + 'HYBRID target densities for interfaces', units="kg m-3") endif do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 if (G%mask2dT(i,j)>0.) then @@ -631,7 +637,9 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & "index_map out of range") endif - rho_target_3d(i,j,:) = rho_target_2d(:,nint(index_map(i,j))) + do k=1,ke+1 + rho_target_3d(i,j,k) = rho_target_2d(k,nint(index_map(i,j))) + enddo endif !mask2dT enddo; enddo varName = trim( extractWord(trim(string(12:)), 4) ) @@ -787,6 +795,77 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif endif + if (coordinateMode(coord_mode) == REGRIDDING_HYCOM1) then + allocate(dz_shallow(ke)) + call get_param(param_file, mdl, "SHALLOW_"//trim(coord_res_param), dz_shallow, & + "HYBGEN-style Z-sigma-Z near surface fixed coordinate. "//& + "The default of all zeros turns this option off. "//& + "Let N_SIGMA be the number of consecutive non-zero entries, typically < NK. "//& + "Use SHALLOW_"//trim(coord_res_param)//" when rest depth is shallower than "//& + "SUM(SHALLOW_"//trim(coord_res_param)//"(1:N_SIGMA)). "//& + "Use "//trim(coord_res_param)//" when rest depth is deeper than "//& + "SUM("//trim(coord_res_param)//"(1:N_SIGMA)). "//& + "Otherwise use a linear sum of the two weighted by rest depth.",& + units="m", default=0.0) + n_sigma = ke + depth_s = 0.0 + do k= 1,ke + depth_s = depth_s + dz_shallow(k) + if (dz_shallow(k) == 0.0) then + n_sigma = k-1 + exit + endif + enddo + if (n_sigma > 0) then + if (main_parameters) call log_param(param_file, mdl, "!N_SIGMA", n_sigma, & + "Number of consecutive non-zero entries in SHALLOW_"//& + trim(coord_res_param)//".") + if (.not.allocated(dz_3d)) then + allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0) + allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0) + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + do k=1,ke + dz_3d(i,j,k) = dz(k) + enddo + do k=1,ke+1 + rho_target_3d(i,j,k) = rho_target(k) + enddo + endif !mask2dT + enddo; enddo + endif + do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 + if (G%mask2dT(i,j)>0.) then + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*US%Z_to_m + if (nominalDepth <= depth_s) then + do k= 1,n_sigma + dz_3d(i,j,k) = dz_shallow(k) + enddo + do k= n_sigma+1,ke + dz_3d(i,j,k) = dz_shallow(n_sigma) + enddo + else ! >depth_s + depth_d = 0.0 + do k= 1,n_sigma + depth_d = depth_d + dz_3d(i,j,k) + enddo + ! do nothing if nominalDepth >= depth_d + if (nominalDepth < depth_d) then + depth_q = (nominalDepth - depth_s) / (depth_d - depth_s) + do k= 1,n_sigma + dz_3d(i,j,k) = (1.0-depth_q)*dz_shallow(k) + depth_q*dz_3d(i,j,k) + enddo + do k= n_sigma+1,ke + dz_3d(i,j,k) = (1.0-depth_q)*dz_shallow(n_sigma) + depth_q*dz_3d(i,j,k) + enddo + endif !depth_s + endif !nominalDepth + endif !mask2dT + enddo; enddo + endif !n_sigma + deallocate(dz_shallow) + endif !REGRIDDING_HYCOM1 + CS%nk=ke ! Target resolution (for fixed coordinates) @@ -839,7 +918,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then call set_target_densities_from_GV(GV, US, CS) call log_param(param_file, mdl, "!TARGET_DENSITIES", US%R_to_kg_m3*CS%target_density(:), & - 'RHO target densities for interfaces', units=coordinateUnits(coord_mode)) + 'RHO target densities for interfaces', "kg m-3") endif ! initialise coordinate-specific control structure From 3294a2d1af17e4bfbb1fedd2b8167566681ccc09 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Wed, 27 Aug 2025 17:33:45 +0000 Subject: [PATCH 4/5] Non-integer HYBRID_MAP values The 2-d REAL map array in HYBRID_MAP usually contains integer values each referencing one profile. It can instead contain non-integer values of the form I+frac, which indicate a weighted sum of profiles: (1-frac) p(I) + (frac) p(I+1). The same profile can be used multiple times, e.g. if 1st profile is also 4th can get profiles between 1 and 2 and between 1 and 3. HYBRID_3D is more general, but HYBRID_MAP covers most practical uses. --- src/ALE/MOM_regridding.F90 | 52 +++++++++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 15 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 15d1d316f0..853dc33595 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -208,8 +208,11 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! Local variables integer :: ke ! Number of levels - integer :: np ! Number of profiles, for HYBRID_MAP - integer :: n_sigma ! Number of shallow dz's, for HYBRID_MAP or HYBRID_3D + integer :: n_sigma ! Number of shallow dz's, for HYBRID_MAP or HYBRID_3D + integer :: np ! Number of profiles, for HYBRID_MAP + integer :: nceiling ! ceiling of map index, for HYBRID_MAP + integer :: nfloor ! floor of map index, for HYBRID_MAP + real :: nfrac ! fraction of map index, for HYBRID_MAP character(len=80) :: string, string2, varName ! Temporary strings character(len=40) :: coord_units, coord_res_param ! Temporary strings character(len=MAX_PARAM_LENGTH) :: param_name @@ -393,6 +396,9 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & " by a comma or space, for map, sigma-2 and dz.\n"//& " Map is a spatial index array with, maxval(map)=N,\n"//& " and the others are 2D arrays containing N profiles.\n"//& + " Map typically contains integer values, but it can\n"//& + " contain real values, I+w, which imply using\n"//& + " the weighted sum of profiles I and I+1.\n"//& " Dz can be FNC1:string which is used everywhere.\n"//& " e.g. HYBRID_MAP:vgrid.nc,map,sigma2,dz",& default=trim(string2)) @@ -607,11 +613,10 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call MOM_read_data(trim(fileName), trim(varName), index_map, G%Domain) call pass_var(index_map, G%Domain, halo=1) !find maximum index - tmpReal = 1 + np = 1 do j=G%jsc, G%jec ; do i=G%isc, G%iec - tmpReal = max(tmpReal,index_map(i,j)) + np = max(np,ceiling(index_map(i,j))) enddo ; enddo - np = tmpReal call max_across_PEs(np) write(string2,"(i3)") np call MOM_error(NOTE, & @@ -632,14 +637,24 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 if (G%mask2dT(i,j)>0.) then - if (nint(index_map(i,j))<1 .or. nint(index_map(i,j))>np) then + nfloor = floor(index_map(i,j)) + nceiling = ceiling(index_map(i,j)) + if (nfloor<1 .or. nceiling>np) then write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=', index_map(i,j) call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & "index_map out of range") endif - do k=1,ke+1 - rho_target_3d(i,j,k) = rho_target_2d(k,nint(index_map(i,j))) - enddo + if (nfloor == nceiling) then + do k=1,ke+1 + rho_target_3d(i,j,k) = rho_target_2d(k,nfloor) + enddo + else + nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0 + do k=1,ke+1 + rho_target_3d(i,j,k) = (1.0-nfrac)*rho_target_2d(k,nfloor) + & + nfrac *rho_target_2d(k,nceiling) + enddo + endif !integer:else endif !mask2dT enddo; enddo varName = trim( extractWord(trim(string(12:)), 4) ) @@ -686,12 +701,19 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1 if (G%mask2dT(i,j)>0.) then - if (nint(index_map(i,j))<1 .or. nint(index_map(i,j))>np) then - write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=',index_map(i,j) - call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "index_map out of range") - endif - dz_3d(i,j,:) = dz_2d(:,nint(index_map(i,j))) + nfloor = floor(index_map(i,j)) + nceiling = ceiling(index_map(i,j)) + if (nfloor == nceiling) then + do k=1,ke + dz_3d(i,j,k) = dz_2d(k,nfloor) + enddo + else + nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0 + do k=1,ke + dz_3d(i,j,k) = (1.0-nfrac)*dz_2d(k,nfloor) + & + nfrac *dz_2d(k,nceiling) + enddo + endif !integer:else endif !mask2dT enddo; enddo endif !dz From 0deb6ce171fa31da5c907b1815011de6d16ec956 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Sat, 6 Sep 2025 14:39:24 +0000 Subject: [PATCH 5/5] indent continuations, source code <= 100 chars --- src/ALE/MOM_ALE.F90 | 45 +++++---- src/ALE/MOM_regridding.F90 | 202 ++++++++++++++++++++----------------- src/ALE/coord_hycom.F90 | 12 +-- 3 files changed, 142 insertions(+), 117 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index fd07281039..ab9b7405ee 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -207,7 +207,8 @@ subroutine ALE_init( param_file, G, GV, US, max_depth, CS) ! Initialize and configure regridding call ALE_initRegridding(G, 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) + 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, & @@ -322,12 +323,12 @@ subroutine ALE_init( param_file, G, GV, US, max_depth, CS) default=-0.001, units="m", scale=GV%m_to_H) call get_param(param_file, mdl, "REMAP_VEL_MASK_H_THIN", CS%h_vel_mask, & "A thickness at velocity points below which near-bottom layers are zeroed out "//& - "after remapping, following practice with Hybgen remapping, or a negative value "//& - "to avoid such filtering altogether.", & + "after remapping, following practice with Hybgen remapping, "//& + "or a negative value 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 init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS) call get_param(param_file, mdl, "REMAP_VEL_CONSERVE_KE", CS%conserve_ke, & "If true, a correction is applied to the baroclinic component of velocity "//& @@ -641,7 +642,8 @@ end subroutine ALE_offline_inputs !> For a state-based coordinate, accelerate the process of regridding by !! repeatedly applying the grid calculation algorithm -subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, dt, dzRegrid, initial) +subroutine ALE_regrid_accelerated(CS, G, GV, US, 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 @@ -690,7 +692,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d ! initial total interface displacement due to successive regridding if (CS%remap_uv_using_old_alg) & - dzIntTotal(:,:,:) = 0. + dzIntTotal(:,:,:) = 0. call create_group_pass(pass_T_S_h, T, G%domain) call create_group_pass(pass_T_S_h, S, G%domain) @@ -709,7 +711,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d ! Apply timescale to regridding (for e.g. filtered_grid_motion) if (present(dt)) & - call ALE_update_regrid_weights(dt, CS) + call ALE_update_regrid_weights(dt, CS) do itt = 1, n_itt @@ -723,12 +725,14 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface) if (CS%remap_uv_using_old_alg) & - dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) + dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) ! remap from original grid onto new grid do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:)) - call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:)) + call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), & + tv_local%S(i,j,:)) + call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), & + tv_local%T(i,j,:)) enddo ; enddo ! starting grid for next iteration @@ -1147,7 +1151,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u if (CS%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0. if (((CS%id_remap_delta_integ_u2>0) .or. (CS%id_remap_delta_integ_v2>0)) .and. .not.present(dt))& - call MOM_error(FATAL, "ALE KE diagnostics requires passing dt into ALE_remap_velocities") + call MOM_error(FATAL, "ALE KE diagnostics requires passing dt into ALE_remap_velocities") nz = GV%ke @@ -1213,7 +1217,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u endif if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) & - call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) + call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) ! Copy the column of new velocities back to the 3-d array do k=1,nz @@ -1362,13 +1366,14 @@ subroutine ALE_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val) do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB if ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) > 0.0 ) then - I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1))) + I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + & + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1))) do k=1,nz h_src(k) = ((G%mask2dT(i,j) * h_old(i,j,k) + G%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + & - (G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum + (G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum h_tgt(k) = ((G%mask2dT(i,j) * h_new(i,j,k) + G%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + & - (G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum + (G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum enddo do K=1,nz+1 @@ -1550,7 +1555,8 @@ subroutine ALE_PLM_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b ) do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 slp(1) = 0. do k = 2, GV%ke-1 - slp(k) = PLM_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1)) + slp(k) = PLM_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, & + Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1)) enddo slp(GV%ke) = 0. @@ -1563,7 +1569,8 @@ subroutine ALE_PLM_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b ) mslp = - PLM_extrapolate_slope(h(i,j,2), h(i,j,1), h_neglect, Q(i,j,2), Q(i,j,1)) Q_t(i,j,1) = Q(i,j,1) - 0.5 * mslp Q_b(i,j,1) = Q(i,j,1) + 0.5 * mslp - mslp = PLM_extrapolate_slope(h(i,j,GV%ke-1), h(i,j,GV%ke), h_neglect, Q(i,j,GV%ke-1), Q(i,j,GV%ke)) + mslp = PLM_extrapolate_slope(h(i,j,GV%ke-1), h(i,j,GV%ke), h_neglect, & + Q(i,j,GV%ke-1), Q(i,j,GV%ke)) Q_t(i,j,GV%ke) = Q(i,j,GV%ke) - 0.5 * mslp Q_b(i,j,GV%ke) = Q(i,j,GV%ke) + 0.5 * mslp else @@ -1631,7 +1638,7 @@ subroutine TS_PPM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect, & answer_date=CS%answer_date ) if (bdry_extrap) & - call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) + call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) do k = 1,GV%ke S_t(i,j,k) = ppol_E(k,1) @@ -1652,7 +1659,7 @@ subroutine TS_PPM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect, & answer_date=CS%answer_date ) if (bdry_extrap) & - call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) + call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) do k = 1,GV%ke T_t(i,j,k) = ppol_E(k,1) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 853dc33595..9f36ae9d89 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -25,15 +25,20 @@ module MOM_regridding 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_HYBGEN, REGRIDDING_ADAPTIVE -use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap, set_interp_answer_date +use regrid_interp, only : interp_CS_type +use regrid_interp, only : set_interp_scheme, set_interp_extrap, set_interp_answer_date -use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike -use coord_sigma, only : init_coord_sigma, sigma_CS, set_sigma_params, build_sigma_column, end_coord_sigma +use coord_zlike, only : zlike_CS +use coord_zlike, only : init_coord_zlike, set_zlike_params, build_zstar_column, end_coord_zlike +use coord_sigma, only : sigma_CS +use coord_sigma, only : init_coord_sigma, set_sigma_params, build_sigma_column, end_coord_sigma use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho use coord_rho, only : old_inflate_layers_1d -use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom +use coord_hycom, only : hycom_CS +use coord_hycom, only : init_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom use coord_hycom, only : init_3d_coord_hycom -use coord_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt +use coord_adapt, only : adapt_CS +use coord_adapt, only : init_coord_adapt, 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 @@ -337,11 +342,11 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) call set_regrid_params(CS, remap_answer_date=remap_answer_date) call get_param(param_file, mdl, "REGRIDDING_ANSWER_DATE", regrid_answer_date, & - "The vintage of the expressions and order of arithmetic to use for regridding. "//& - "Values below 20190101 result in the use of older, less accurate expressions "//& - "that were in use at the end of 2018. Higher values result in the use of more "//& - "robust and accurate forms of mathematically equivalent expressions.", & - default=default_answer_date, do_not_log=.not.GV%Boussinesq) + "The vintage of the expressions and order of arithmetic to use for regridding. "//& + "Values below 20190101 result in the use of older, less accurate expressions "//& + "that were in use at the end of 2018. Higher values result in the use of more "//& + "robust and accurate forms of mathematically equivalent expressions.", & + default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) regrid_answer_date = max(regrid_answer_date, 20230701) call set_regrid_params(CS, regrid_answer_date=regrid_answer_date) endif @@ -421,8 +426,8 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif allocate(dz(ke)) dz(:) = uniformResolution(ke, coord_mode, tmpReal, & - US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), & - US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) ) + US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), & + US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) ) if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, & trim(message), units=trim(coord_units)) elseif (trim(string)=='PARAM') then @@ -512,25 +517,25 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! The following assumes the FILE: syntax of above but without "FILE:" in the string varName = trim( extractWord(trim(string(8:)), 3) ) if (varname == " ") call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID "// & - "Too few arguments in ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID "// & + "Too few arguments in ("//trim(string)//")") fileName = trim( extractWord(trim(string(8:)), 1) ) if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) if (.not. file_exists(fileName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") varName = trim( extractWord(trim(string(8:)), 2) ) if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), rho_target) varName = trim( extractWord(trim(string(8:)), 3) ) if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz ) else ! Read dz from file if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), dz) endif if (main_parameters) then @@ -546,17 +551,17 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! The following assumes the FILE: syntax of above but without "FILE:" in the string varName = trim( extractWord(trim(string(11:)), 3) ) if (varname == " ") call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_3D "// & - "Too few arguments in ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Too few arguments in ("//trim(string)//")") fileName = trim( extractWord(trim(string(11:)), 1) ) if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) if (.not. file_exists(fileName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_3D "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") varName = trim( extractWord(trim(string(11:)), 2) ) if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_3D "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), rho_target_3d, G%Domain) call pass_var(rho_target_3d, G%Domain, halo=1) varName = trim( extractWord(trim(string(11:)), 3) ) @@ -581,15 +586,15 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif else ! Read dz from file if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_3D "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), dz_3d, G%Domain) call pass_var(dz_3d, G%Domain, halo=1) ! set nominal 1-d dz to UNIFORM allocate(dz(ke)) dz(:) = uniformResolution(ke, coord_mode, maximum_depth, & - US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), & - US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) ) + US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), & + US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) ) endif !dz elseif (index(trim(string),'HYBRID_MAP:')==1) then ke = GV%ke @@ -599,17 +604,17 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! The following assumes the FILE: syntax of above but without "FILE:" in the string varName = trim( extractWord(trim(string(12:)), 4) ) if (varname == " ") call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_3D "// & - "Too few arguments in ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_3D "// & + "Too few arguments in ("//trim(string)//")") fileName = trim( extractWord(trim(string(12:)), 1) ) if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName ) if (.not. file_exists(fileName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") varName = trim( extractWord(trim(string(12:)), 2) ) if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") call MOM_read_data(trim(fileName), trim(varName), index_map, G%Domain) call pass_var(index_map, G%Domain, halo=1) !find maximum index @@ -620,15 +625,15 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call max_across_PEs(np) write(string2,"(i3)") np call MOM_error(NOTE, & - trim(mdl)//", initialize_regridding: HYBRID_MAP NP="//trim(string2)) + trim(mdl)//", initialize_regridding: HYBRID_MAP NP="//trim(string2)) if (np<1) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_MAP to small NP from "//trim(varName)) + trim(mdl)//", initialize_regridding: HYBRID_MAP to small NP from "//trim(varName)) allocate(dz_2d(ke,np)) allocate(rho_target_2d(ke+1,np)) varName = trim( extractWord(trim(string(12:)), 3) ) if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") ! MOM_read_data can't handle this array call read_variable(trim(fileName), trim(varName), rho_target_2d) if (main_parameters) then @@ -642,7 +647,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & if (nfloor<1 .or. nceiling>np) then write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=', index_map(i,j) call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "index_map out of range") + "index_map out of range") endif if (nfloor == nceiling) then do k=1,ke+1 @@ -679,8 +684,8 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif else ! Read dz from file if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, & - trim(mdl)//", initialize_regridding: HYBRID_MAP "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + trim(mdl)//", initialize_regridding: HYBRID_MAP "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") ! MOM_read_data can't handle this array call read_variable(trim(fileName), trim(varName), dz_2d) if (main_parameters) then @@ -791,7 +796,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & if (ke > size(woa09_dz_approx)) dz(ke) = dz_extra else call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Unrecognized coordinate configuration"//trim(string)) + "Unrecognized coordinate configuration"//trim(string)) endif if (main_parameters) then @@ -810,7 +815,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & dz(ke) = dz(ke) + ( maximum_depth - tmpReal ) else call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) + "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string)) endif endif endif !allocated(dz) @@ -948,7 +953,8 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & if (coord_is_state_dependent) then if (main_parameters) then - call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), P_Ref, & + call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), & + P_Ref, & "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.", & @@ -961,8 +967,8 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & "This is only used for the RHO coordinate.", & units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2) endif - call get_param(param_file, mdl, create_coord_param(param_prefix, "REGRID_COMPRESSIBILITY_FRACTION", & - param_suffix), tmpReal, & + call get_param(param_file, mdl, create_coord_param(param_prefix, & + "REGRID_COMPRESSIBILITY_FRACTION", param_suffix), tmpReal, & "When interpolating potential density profiles we can add "//& "some artificial compressibility solely to make homogeneous "//& "regions appear stratified.", units="nondim", default=0.) @@ -1045,19 +1051,21 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! Otherwise assume we should look for the file in INPUTDIR fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) ) endif - if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)// & + ", initialize_regridding: "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") do_sum = .false. varName = trim( extractWord(trim(string(6:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)// & + ", initialize_regridding: "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") if (len_trim(varName)==0) then if (field_exists(fileName,'z_max')) then; varName = 'z_max' elseif (field_exists(fileName,'dz')) then; varName = 'dz' ; do_sum = .true. elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' ; do_sum = .true. else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") + "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") endif endif if (do_sum) then @@ -1077,7 +1085,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call set_regrid_max_depths(CS, z_max, GV%m_to_H) else call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string)) + "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string)) endif deallocate(z_max) deallocate(dz_max) @@ -1110,17 +1118,19 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & ! Otherwise assume we should look for the file in INPUTDIR fileName = trim(inputdir) // trim( extractWord(trim(longString(6:200)), 1) ) endif - if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(longString)//")") + if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)// & + ", initialize_regridding: "// & + "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(longString)//")") varName = trim( extractWord(trim(longString(6:)), 2) ) - if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(longString)//")") + if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)// & + ", initialize_regridding: "// & + "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(longString)//")") if (len_trim(varName)==0) then if (field_exists(fileName,'h_max')) then; varName = 'h_max' elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") + "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.") endif endif call MOM_read_data(trim(fileName), trim(varName), h_max) @@ -1134,7 +1144,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & call set_regrid_max_thickness(CS, h_max, GV%m_to_H) else call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// & - "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longString)) + "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longString)) endif deallocate(h_max) endif @@ -1394,7 +1404,7 @@ subroutine check_grid_column( nk, h, dzInterface, msg ) write(0,*) 'k,h,hnew=',k,h(k),h_new write(0,*) 'dzI(k+1),dzI(k)=',dzInterface(k+1),dzInterface(k) call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//& - 'Negative layer thickness implied by re-gridding, '//trim(msg)) + 'Negative layer thickness implied by re-gridding, '//trim(msg)) endif total_h_new = total_h_new + h_new @@ -1409,14 +1419,14 @@ subroutine check_grid_column( nk, h, dzInterface, msg ) write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//& - 'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg)) + 'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg)) endif ! Check that the top and bottom are intentionally moving if (dzInterface(1) /= 0.) call MOM_error( FATAL, & - 'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg)) + 'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg)) if (dzInterface(nk+1) /= 0.) call MOM_error( FATAL, & - 'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg)) + 'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg)) end subroutine check_grid_column @@ -1472,11 +1482,11 @@ subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) if (debug) then do k=2,CS%nk+1 if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.") + call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.") enddo do k=2,nk+1 if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.") + call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.") enddo ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0 endif @@ -1550,9 +1560,9 @@ subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) ! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k) ! ! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) & -! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).") +! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).") ! if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) & -! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.") +! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.") ! endif endif @@ -1568,7 +1578,7 @@ subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) enddo do k=2,CS%nk+1 if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_output is tangled.") + call MOM_error(FATAL, "filtered_grid_motion: z_output is tangled.") enddo endif @@ -1960,7 +1970,8 @@ end subroutine build_rho_grid !! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 } -subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, frac_shelf_h, zScale ) +subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, & + frac_shelf_h, zScale ) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -2142,8 +2153,8 @@ subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) ) if (h_new < -3.0*h_err) then write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), & - 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), & - 'h_new=',h_new,'h_err=',h_err + 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), & + 'h_new=',h_new,'h_err=',h_err call MOM_error( FATAL, 'MOM_regridding: adjust_interface_motion() - '//& 'implied h<0 is larger than roundoff!') endif @@ -2154,8 +2165,8 @@ subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) h_new = ( dz_int(k) - dz_int(k+1) ) if (h_new < -3.0*h_err) then write(0,*) 'h<0 at k=',k,'h_old was empty',& - 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), & - 'h_new=',h_new,'h_err=',h_err + 'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), & + 'h_new=',h_new,'h_err=',h_err call MOM_error( FATAL, 'MOM_regridding: adjust_interface_motion() - '//& 'implied h<0 is larger than roundoff!') endif @@ -2164,14 +2175,14 @@ subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) do k = min(CS%nk,nk),2,-1 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) ) if (h_new max_depths(K)) & - call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths!") + do K=1,CS%nk + if (max_depths(K+1) > max_depths(K)) & + call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths!") enddo else - do K=1,CS%nk ; if (max_depths(K+1) < max_depths(K)) & - call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths.") + do K=1,CS%nk + if (max_depths(K+1) < max_depths(K)) & + call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths.") enddo endif @@ -2655,8 +2668,8 @@ function getCoordinateInterfaces( CS, undo_scaling ) ! densities, rather than computing the interfaces based on resolution if (CS%regridding_scheme == REGRIDDING_RHO) then if (.not. CS%target_density_set) & - call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//& - 'target densities not set!') + call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//& + 'target densities not set!') if (unscale) then getCoordinateInterfaces(:) = CS%coord_scale * CS%target_density(:) @@ -2786,10 +2799,13 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri call MOM_error(FATAL,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!') CS%old_grid_weight = old_grid_weight endif - if (present(depth_of_time_filter_shallow)) CS%depth_of_time_filter_shallow = depth_of_time_filter_shallow - if (present(depth_of_time_filter_deep)) CS%depth_of_time_filter_deep = depth_of_time_filter_deep + if (present(depth_of_time_filter_shallow)) CS%depth_of_time_filter_shallow = & + depth_of_time_filter_shallow + if (present(depth_of_time_filter_deep)) CS%depth_of_time_filter_deep = & + depth_of_time_filter_deep if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then - if (CS%depth_of_time_filter_deep Calculate interface density anomaly w.r.t. the target. -subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, p_col, & - R, RiAnom, h_neglect, h_neglect_edge) +subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, ix, jy, depth, h, T, S, & + p_col, R, RiAnom, h_neglect, h_neglect_edge) type(hycom_CS), intent(in) :: CS !< Coordinate control structure type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure