diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 93c4d35faa..ab9b7405ee 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,8 +206,9 @@ 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 regridding_preadjust_reqs(CS%regridCS, CS%do_conv_adj, CS%use_hybgen_unmix, hybgen_CS=hybgen_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. call get_param(param_file, mdl, "REMAPPING_SCHEME", string, & @@ -321,12 +323,12 @@ subroutine ALE_init( param_file, 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 "//& @@ -640,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 @@ -689,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) @@ -708,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 @@ -722,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 @@ -1146,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 @@ -1212,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 @@ -1361,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 @@ -1549,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. @@ -1562,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 @@ -1630,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) @@ -1651,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) @@ -1664,7 +1672,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 +1689,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..9f36ae9d89 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,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, pass_var use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real use MOM_remapping, only : remapping_CS @@ -23,14 +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_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt +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 : 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 @@ -66,6 +74,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 +197,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 +213,11 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! Local variables integer :: ke ! Number of levels + 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 @@ -213,13 +234,25 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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] 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 :: 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] 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 +297,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) @@ -309,11 +342,11 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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 @@ -356,8 +389,23 @@ 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"//& + " 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)) message = "The distribution of vertical resolution for the target\n"//& "grid used for Eulerian-like coordinates. For example,\n"//& @@ -378,8 +426,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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 @@ -424,11 +472,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,31 +511,220 @@ 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 "// & - "Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")") + 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 "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + 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 "// & - "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")") + 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 if (main_parameters) then 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 + 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 + np = 1 + do j=G%jsc, G%jec ; do i=G%isc, G%iec + np = max(np,ceiling(index_map(i,j))) + enddo ; enddo + 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="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 + 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 + 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) ) + 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 + 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 + 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. @@ -557,7 +796,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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 @@ -566,31 +805,116 @@ 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 + 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) - 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,39 +936,42 @@ 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(:), & - 'RHO target densities for interfaces', units=coordinateUnits(coord_mode)) + 'RHO target densities for interfaces', "kg m-3") 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 - 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.", & 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 +987,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 @@ -724,19 +1051,21 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! 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 @@ -756,7 +1085,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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) @@ -789,17 +1118,19 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m ! 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) @@ -813,7 +1144,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m 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 @@ -835,6 +1166,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 ) @@ -1071,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 @@ -1086,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 @@ -1149,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 @@ -1227,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 @@ -1245,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 @@ -1448,10 +1781,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!!!' ) @@ -1635,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 @@ -1701,7 +2037,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 +2119,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 @@ -1816,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 @@ -1828,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 @@ -1838,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 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 +2372,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 +2410,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 +2453,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 @@ -2120,12 +2501,14 @@ subroutine set_regrid_max_depths( CS, max_depths, units_to_H ) ! Check for sign reversals in the depths. if (max_depths(CS%nk+1) < max_depths(1)) then - 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 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 @@ -2285,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(:) @@ -2376,7 +2759,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 @@ -2415,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 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..f5062d6f68 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) @@ -85,14 +138,14 @@ subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, only_ if (present(max_interface_depths)) then if (size(max_interface_depths) /= CS%nk+1) & - call MOM_error(FATAL, "set_hycom_params: max_interface_depths inconsistent size") + call MOM_error(FATAL, "set_hycom_params: max_interface_depths inconsistent size") allocate(CS%max_interface_depths(CS%nk+1)) CS%max_interface_depths(:) = max_interface_depths(:) endif if (present(max_layer_thickness)) then if (size(max_layer_thickness) /= CS%nk) & - call MOM_error(FATAL, "set_hycom_params: max_layer_thickness inconsistent size") + call MOM_error(FATAL, "set_hycom_params: max_layer_thickness inconsistent size") allocate(CS%max_layer_thickness(CS%nk)) CS%max_layer_thickness(:) = max_layer_thickness(:) endif @@ -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,20 +219,25 @@ 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) ) z_nz = 0.5 * ( z_col(nz) + z_col(nz+1) ) do k = 1,CS%nk - p_col_new(k) = p_col(1) + ( 0.5 * ( z_col_new(K) + z_col_new(K+1) ) - z_1 ) / ( z_nz - z_1 ) * & - ( p_col(nz) - p_col(1) ) + p_col_new(k) = p_col(1) + ( 0.5 * ( z_col_new(K) + z_col_new(K+1) ) - z_1 ) & + / ( z_nz - z_1 ) * ( p_col(nz) - p_col(1) ) enddo ! 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, & - 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 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, &