diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a341fd1835..61ab6c93cf 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -483,7 +483,7 @@ subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_ ! Build the new grid and store it in h_new. The old grid is retained as h. ! Both are needed for the subsequent remapping of variables. dzRegrid(:,:,:) = 0.0 - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, & + call regridding_main( CS%remapCS, CS%regridCS, G, GV, US, h, tv, h_new, dzRegrid, & frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell) if (CS%id_dzRegrid>0) then ; if (query_averaging_enabled(CS%diag)) then @@ -497,10 +497,11 @@ end subroutine ALE_regrid !> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have !! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This !! routine builds a grid on the runtime specified vertical coordinate -subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) +subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) type(ALE_CS), pointer :: CS !< Regridding parameters and options type(ocean_grid_type), intent(in ) :: G !< Ocean grid informations type(verticalGrid_type), intent(in ) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure type(tracer_registry_type), pointer :: Reg !< Tracer registry structure @@ -526,7 +527,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) ! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored ! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective ! adjustment right now is not used because it is unclear what to do with vanished layers - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid) + call regridding_main( CS%remapCS, CS%regridCS, G, GV, US, h, tv, h_new, dzRegrid) if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_offline_inputs)") ! Remap all variables from old grid h onto new grid h_new @@ -576,10 +577,11 @@ 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, 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 + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Original thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS) @@ -651,7 +653,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d ! generate new grid if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local) - call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface) + call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface) dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) ! remap from original grid onto new grid diff --git a/src/ALE/MOM_hybgen_regrid.F90 b/src/ALE/MOM_hybgen_regrid.F90 index dc7c90a079..524f9b8ff2 100644 --- a/src/ALE/MOM_hybgen_regrid.F90 +++ b/src/ALE/MOM_hybgen_regrid.F90 @@ -61,7 +61,21 @@ module MOM_hybgen_regrid !> Nominal density of interfaces [R ~> kg m-3] real, allocatable, dimension(:) :: target_density - real :: onem !< Nominally one m in thickness units [H ~> m or kg m-2] + real :: dp_far_from_sfc !< A distance that determines when an interface is suffiently far from + !! the surface that certain adjustments can be made in the Hybgen regridding + !! code [H ~> m or kg m-2]. In Hycom, this is set to tenm (nominally 10 m). + real :: dp_far_from_bot !< A distance that determines when an interface is suffiently far from + !! the bottom that certain adjustments can be made in the Hybgen regridding + !! code [H ~> m or kg m-2]. In Hycom, this is set to onem (nominally 1 m). + real :: h_thin !< A layer thickness below which a layer is considered to be too thin for + !! certain adjustments to be made in the Hybgen regridding code. + !! In Hycom, this is set to onemm (nominally 0.001 m). + + real :: rho_eps !< A small nonzero density that is used to prevent division by zero + !! in several expressions in the Hybgen regridding code [R ~> kg m-3]. + + real :: onem !< Nominally one m in thickness units [H ~> m or kg m-2], used only in + !! certain debugging tests. end type hybgen_regrid_CS @@ -166,6 +180,28 @@ subroutine init_hybgen_regrid(CS, GV, US, param_file) "A bottom boundary layer thickness within which Hybgen is able to move "//& "overlying layers upward to match a target density.", & units="m", default=0.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_FAR_FROM_SURFACE", CS%dp_far_from_sfc, & + "A distance that determines when an interface is suffiently far "//& + "from the surface that certain adjustments can be made in the Hybgen "//& + "regridding code. In Hycom, this is set to tenm (nominally 10 m).", & + units="m", default=10.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_FAR_FROM_BOTTOM", CS%dp_far_from_bot, & + "A distance that determines when an interface is suffiently far "//& + "from the bottom that certain adjustments can be made in the Hybgen "//& + "regridding code. In Hycom, this is set to onem (nominally 1 m).", & + units="m", default=1.0, scale=GV%m_to_H) + call get_param(param_file, mdl, "HYBGEN_H_THIN", CS%h_thin, & + "A layer thickness below which a layer is considered to be too thin for "//& + "certain adjustments to be made in the Hybgen regridding code. "//& + "In Hycom, this is set to onemm (nominally 0.001 m).", & + units="m", default=0.001, scale=GV%m_to_H) + + call get_param(param_file, mdl, "HYBGEN_DENSITY_EPSILON", CS%rho_eps, & + "A small nonzero density that is used to prevent division by zero "//& + "in several expressions in the Hybgen regridding code.", & + units="kg m-3", default=1e-11, scale=US%kg_m3_to_R) + + call get_param(param_file, mdl, "HYBGEN_REMAP_DENSITY_MATCH", CS%hybiso, & "A tolerance between the layer densities and their target, within which "//& "Hybgen determines that remapping uses PCM for a layer.", & @@ -300,12 +336,17 @@ end subroutine get_hybgen_regrid_params !> Modify the input grid to give a new vertical grid based on the HYCOM hybgen code. -subroutine hybgen_regrid(G, GV, US, dp, tv, CS, dzInterface, PCM_cell) +subroutine hybgen_regrid(G, GV, US, dp, nom_depth_H, tv, CS, dzInterface, PCM_cell) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: dp !< Source grid layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure type(hybgen_regrid_CS), intent(in) :: CS !< hybgen control structure real, dimension(SZI_(G),SZJ_(G),CS%nk+1), & @@ -457,7 +498,7 @@ subroutine hybgen_regrid(G, GV, US, dp, tv, CS, dzInterface, PCM_cell) enddo ! The following block of code is used to trigger z* stretching of the targets heights. - nominalDepth = (G%bathyT(i,j) + G%Z_ref)*GV%Z_to_H + nominalDepth = nom_depth_H(i,j) if (h_tot <= CS%min_dilate*nominalDepth) then dilate = CS%min_dilate elseif (h_tot >= CS%max_dilate*nominalDepth) then @@ -482,8 +523,7 @@ subroutine hybgen_regrid(G, GV, US, dp, tv, CS, dzInterface, PCM_cell) enddo !k ! Determine the new layer thicknesses. - call hybgen_column_regrid(CS, nk, CS%thkbot, CS%onem, & - 1.0e-11*US%kg_m3_to_R, Rcv_tgt, fixlay, qhrlx, dp0ij, & + call hybgen_column_regrid(CS, nk, CS%thkbot, Rcv_tgt, fixlay, qhrlx, dp0ij, & dp0cum, Rcv, h_col, dz_int) ! Store the output from hybgenaij_regrid in 3-d arrays. @@ -669,13 +709,11 @@ real function cushn(delp, dp0) end function cushn !> Create a new grid for a column of water using the Hybgen algorithm. -subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & +subroutine hybgen_column_regrid(CS, nk, thkbot, Rcv_tgt, & fixlay, qhrlx, dp0ij, dp0cum, Rcv, h_in, dp_int) type(hybgen_regrid_CS), intent(in) :: CS !< hybgen regridding control structure integer, intent(in) :: nk !< number of layers real, intent(in) :: thkbot !< thickness of bottom boundary layer [H ~> m or kg m-2] - real, intent(in) :: onem !< one m in pressure units [H ~> m or kg m-2] - real, intent(in) :: epsil !< small nonzero density to prevent division by zero [R ~> kg m-3] real, intent(in) :: Rcv_tgt(nk) !< Target potential density [R ~> kg m-3] integer, intent(in) :: fixlay !< deepest fixed coordinate layer real, intent(in) :: qhrlx( nk+1) !< relaxation coefficient per timestep [nondim] @@ -702,20 +740,14 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & real :: h_hat0 ! A first guess at thickness movement upward across the interface ! between layers k and k-1 [H ~> m or kg m-2] real :: dh_cor ! Thickness changes [H ~> m or kg m-2] - real :: tenm ! ten m in pressure units [H ~> m or kg m-2] - real :: onemm ! one mm in pressure units [H ~> m or kg m-2] logical :: trap_errors integer :: k character(len=256) :: mesg ! A string for output messages ! This line needs to be consistent with the parameters set in cushn(). - real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn -! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn -! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn - - !### These hard-coded parameters should be changed to run-time variables. - tenm = 10.0*onem - onemm = 0.001*onem + real, parameter :: qqmn=-4.0, qqmx=2.0 ! shifted range for cushn [nondim] +! real, parameter :: qqmn=-2.0, qqmx=4.0 ! traditional range for cushn [nondim] +! real, parameter :: qqmn=-4.0, qqmx=6.0 ! somewhat wider range for cushn [nondim] trap_errors = .true. @@ -769,26 +801,26 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & ! Remap the non-fixed layers. - ! In the Hycom version, this loop was fused the loop correcting water that is + ! In the Hycom version, this loop was fused with the loop correcting water that is ! too light, and it ran down the water column, but if there are a set of layers ! that are very dense, that structure can lead to all of the water being remapped ! into a single thick layer. Splitting the loops and running the loop upwards - ! (as is done here avoids that catastrophic problem for layers that are far from + ! (as is done here) avoids that catastrophic problem for layers that are far from ! their targets. However, this code is still prone to a thin-thick-thin null mode. do k=nk,fixlay+2,-1 ! This is how the Hycom code would do this loop: do k=fixlay+1,nk ; if (k>fixlay+1) then - if ((Rcv(k) > Rcv_tgt(k) + epsil)) then + if ((Rcv(k) > Rcv_tgt(k) + CS%rho_eps)) then ! Water in layer k is too dense, so try to dilute with water from layer k-1 ! Do not move interface if k = fixlay + 1 if ((Rcv(k-1) >= Rcv_tgt(k-1)) .or. & - (p_int(k) <= dp0cum(k) + onem) .or. & + (p_int(k) <= dp0cum(k) + CS%dp_far_from_bot) .or. & (h_col(k) <= h_col(k-1))) then ! If layer k-1 is too light, there is a conflict in the direction the ! inteface between them should move, so thicken the thinner of the two. - if ((Rcv_tgt(k) - Rcv(k-1)) <= epsil) then + if ((Rcv_tgt(k) - Rcv(k-1)) <= CS%rho_eps) then ! layer k-1 is far too dense, take the entire layer ! If this code is working downward and this branch is repeated in a series ! of successive layers, it can accumulate into a very thick homogenous layers. @@ -814,7 +846,7 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & ! layer (thinner than its minimum thickness) in the interior ocean, ! move interface k-1 (and k-2 if necessary) upward ! Only work on layers that are sufficiently far from the fixed near-surface layers. - if ((h_hat >= 0.0) .and. (k > fixlay+2) .and. (p_int(k-1) > dp0cum(k-1) + tenm)) then + if ((h_hat >= 0.0) .and. (k > fixlay+2) .and. (p_int(k-1) > dp0cum(k-1) + CS%dp_far_from_sfc)) then ! Only act if interface k-1 is near the bottom or layer k-2 could donate water. if ( (p_int(nk+1) - p_int(k-1) < thkbot) .or. & @@ -828,7 +860,7 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2) endif !fixlay+3:else - if (h_hat2 < -onemm) then + if (h_hat2 < -CS%h_thin) then dh_cor = qhrlx(k-1) * max(h_hat2, -h_hat - h_col(k-1)) h_col(k-2) = h_col(k-2) + dh_cor h_col(k-1) = h_col(k-1) - dh_cor @@ -838,9 +870,9 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & h_hat = cushn(h_hat0 + h_col(k-1), dp0ij(k-1)) - h_col(k-1) elseif (k <= fixlay+3) then ! Do nothing. - elseif (p_int(k-2) > dp0cum(k-2) + tenm .and. & - (p_int(nk+1) - p_int(k-2) < thkbot .or. & - h_col(k-3) > qqmx*dp0ij(k-3))) then + elseif ( (p_int(k-2) > dp0cum(k-2) + CS%dp_far_from_sfc) .and. & + ( (p_int(nk+1) - p_int(k-2) < thkbot) .or. & + (h_col(k-3) > qqmx*dp0ij(k-3)) ) ) then ! Determine how much water layer k-3 could supply without becoming too thin. if (k == fixlay+4) then @@ -850,7 +882,7 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & ! Maintain minimum thickess of layer k-3. h_hat3 = cushn(h_col(k-3) + (h_hat0 - h_hat), dp0ij(k-3)) - h_col(k-3) endif !fixlay+4:else - if (h_hat3 < -onemm) then + if (h_hat3 < -CS%h_thin) then ! Water is moved from layer k-3 to k-2, but do not dilute layer k-2 too much. dh_cor = qhrlx(k-2) * max(h_hat3, -h_col(k-2)) h_col(k-3) = h_col(k-3) + dh_cor @@ -860,7 +892,7 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & ! Now layer k-2 might be able donate to layer k-1. h_hat2 = cushn(h_col(k-2) + (h_hat0 - h_hat), dp0ij(k-2)) - h_col(k-2) - if (h_hat2 < -onemm) then + if (h_hat2 < -CS%h_thin) then dh_cor = qhrlx(k-1) * (max(h_hat2, -h_hat - h_col(k-1)) ) h_col(k-2) = h_col(k-2) + dh_cor h_col(k-1) = h_col(k-1) - dh_cor @@ -890,17 +922,17 @@ subroutine hybgen_column_regrid(CS, nk, thkbot, onem, epsil, Rcv_tgt, & enddo do k=fixlay+1,nk - if (Rcv(k) < Rcv_tgt(k) - epsil) then ! layer too light + if (Rcv(k) < Rcv_tgt(k) - CS%rho_eps) then ! layer too light ! Water in layer k is too light, so try to dilute with water from layer k+1. ! Entrainment is not possible if layer k touches the bottom. if (p_int(k+1) < p_int(nk+1)) then ! k 1.0e-13*max(p_int(nk+1), onem)) then + if (abs((h_col(k) - h_in(k)) + (dp_int(K) - dp_int(K+1))) > 1.0e-13*max(p_int(nk+1), CS%onem)) then write(mesg, '("k ",i4," h ",es13.4," h_in ",es13.4, " dp ",2es13.4," err ",es13.4)') & k, h_col(k), h_in(k), dp_int(K), dp_int(K+1), (h_col(k) - h_in(k)) + (dp_int(K) - dp_int(K+1)) call MOM_error(FATAL, "Mismatched thickness changes in hybgen_regrid: "//trim(mesg)) endif - if (h_col(k) < 0.0) then ! Could instead do: -1.0e-15*max(p_int(nk+1), onem)) then + if (h_col(k) < 0.0) then ! Could instead do: -1.0e-15*max(p_int(nk+1), CS%onem)) then write(mesg, '("k ",i4," h ",es13.4," h_in ",es13.4, " dp ",2es13.4, " fixlay ",i4)') & k, h_col(k), h_in(k), dp_int(K), dp_int(K+1), fixlay call MOM_error(FATAL, "Significantly negative final thickness in hybgen_regrid: "//trim(mesg)) endif enddo do K=1,nk+1 - if (abs(dp_int(K) - (p_int(K) - pres_in(K))) > 1.0e-13*max(p_int(nk+1), onem)) then + if (abs(dp_int(K) - (p_int(K) - pres_in(K))) > 1.0e-13*max(p_int(nk+1), CS%onem)) then call MOM_error(FATAL, "Mismatched interface height changes in hybgen_regrid.") endif enddo diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9da4e95b24..c5f5807f66 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -175,7 +175,7 @@ module MOM_regridding character(len=*), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2" !> Default mode for boundary extrapolation logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false. -!> Default minimum thickness for some coordinate generation modes +!> Default minimum thickness for some coordinate generation modes [m] real, parameter, public :: regriddingDefaultMinThickness = 1.e-3 !> Maximum length of parameters @@ -213,7 +213,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m logical :: remap_answers_2018 integer :: remap_answer_date ! The vintage of the remapping expressions to use. integer :: regrid_answer_date ! The vintage of the regridding expressions to use. - real :: tmpReal, P_Ref + real :: tmpReal ! A temporary variable used in setting other variables [various] + real :: P_Ref ! The coordinate variable reference pression [R L2 T-2 ~> Pa] real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z). real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3] @@ -771,7 +772,7 @@ end subroutine end_regridding !------------------------------------------------------------------------------ !> Dispatching regridding routine for orchestrating regridding & remapping -subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, & +subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, & frac_shelf_h, PCM_cell) !------------------------------------------------------------------------------ ! This routine takes care of (1) building a new grid and (2) remapping between @@ -795,45 +796,60 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, & type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after - !! the last time step + !! the last time step [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamical variables (T, S, ...) - real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate - real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage + real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target + !! coordinate [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each + !! interface [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage [nomdim] logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true ! Local variables + real :: nom_depth_H(SZI_(G),SZJ_(G)) !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2] + real :: Z_to_H ! A conversion factor used by some routines to convert coordinate + ! parameters to depth units [H Z-1 ~> nondim or kg m-3] real :: trickGnuCompiler integer :: i, j if (present(PCM_cell)) PCM_cell(:,:,:) = .false. + Z_to_H = US%Z_to_m * GV%m_to_H ! Often this is equivalent to GV%Z_to_H. + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * Z_to_H + ! Consider using the following instead: + ! nom_depth_H(i,j) = max( (G%bathyT(i,j)+G%Z_ref) * Z_to_H , CS%min_nom_depth ) + ! if (G%mask2dT(i,j)==0.) nom_depth_H(i,j) = 0.0 + enddo ; enddo + select case ( CS%regridding_scheme ) case ( REGRIDDING_ZSTAR ) - call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h ) + call build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_h, zScale=Z_to_H ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_SIGMA_SHELF_ZSTAR) - call build_zstar_grid( CS, G, GV, h, dzInterface ) + call build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, zScale=Z_to_H ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_SIGMA ) - call build_sigma_grid( CS, G, GV, h, dzInterface ) + call build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_RHO ) - call build_rho_grid( G, GV, G%US, h, tv, dzInterface, remapCS, CS, frac_shelf_h ) + call build_rho_grid( G, GV, G%US, h, nom_depth_H, tv, dzInterface, remapCS, CS, frac_shelf_h ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ARBITRARY ) - call build_grid_arbitrary( G, GV, h, dzInterface, trickGnuCompiler, CS ) + call build_grid_arbitrary( G, GV, h, nom_depth_H, dzInterface, trickGnuCompiler, CS ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_HYCOM1 ) - call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, remapCS, CS, frac_shelf_h ) + call build_grid_HyCOM1( G, GV, G%US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, & + frac_shelf_h, zScale=Z_to_H ) case ( REGRIDDING_HYBGEN ) - call hybgen_regrid(G, GV, G%US, h, tv, CS%hybgen_CS, dzInterface, PCM_cell) + call hybgen_regrid(G, GV, G%US, h, nom_depth_H, tv, CS%hybgen_CS, dzInterface, PCM_cell) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ADAPTIVE ) - call build_grid_adaptive(G, GV, G%US, h, tv, dzInterface, remapCS, CS) + call build_grid_adaptive(G, GV, G%US, h, nom_depth_H, tv, dzInterface, remapCS, CS) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case default @@ -896,9 +912,12 @@ subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses (arbitrary units) - real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions (same as h) - real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses (same as h) + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses [H ~> m or kg m-2] + !! or other units + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions + !! in the same units as h [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses in the same + !! units as h [H ~> m or kg m-2] ! Local variables integer :: i, j, k, nki @@ -1121,21 +1140,29 @@ end subroutine filtered_grid_motion !> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004). !! z* is defined as !! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H . -subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) +subroutine build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_h, zScale) ! Arguments type(regridding_CS), intent(in) :: CS !< Regridding control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional !! ice shelf coverage [nondim]. + real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate + !! resolution in Z to desired units for zInterface, + !! usually Z_to_H in which case it is in + !! units of [H Z-1 ~> nondim or kg m-3] ! Local variables real :: nominalDepth, minThickness, totalThickness ! Depths and thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ - real :: dh ! [H ~> m or kg m-2] + real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2] #endif real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2] real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2] @@ -1146,13 +1173,13 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) minThickness = CS%min_thickness ice_shelf = present(frac_shelf_h) -!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, & -!$OMP ice_shelf,minThickness) & -!$OMP private(nominalDepth,totalThickness, & + !$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, & + !$OMP ice_shelf,minThickness,zScale,nom_depth_H) & + !$OMP private(nominalDepth,totalThickness, & #ifdef __DO_SAFETY_CHECKS__ -!$OMP dh, & + !$OMP dh, & #endif -!$OMP zNew,zOld) + !$OMP zNew,zOld) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 @@ -1161,8 +1188,8 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) cycle endif - ! Local depth (G%bathyT is positive downward) - nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H + ! Local depth (positive downward) + nominalDepth = nom_depth_H(i,j) ! Determine water column thickness totalThickness = 0.0 @@ -1170,23 +1197,26 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) totalThickness = totalThickness + h(i,j,k) enddo + ! if (GV%Boussinesq) then zOld(nz+1) = - nominalDepth do k = nz,1,-1 zOld(k) = zOld(k+1) + h(i,j,k) enddo + ! else ! Work downward? + ! endif if (ice_shelf) then if (frac_shelf_h(i,j) > 0.) then ! under ice shelf call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, zNew, & - z_rigid_top = totalThickness-nominalDepth, & - eta_orig=zOld(1), zScale=GV%Z_to_H) + z_rigid_top=totalThickness-nominalDepth, & + eta_orig=zOld(1), zScale=zScale) else call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, & - zNew, zScale=GV%Z_to_H) + zNew, zScale=zScale) endif else call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, & - zNew, zScale=GV%Z_to_H) + zNew, zScale=zScale) endif ! Calculate the final change in grid position after blending new and old grids @@ -1225,7 +1255,7 @@ end subroutine build_zstar_grid !------------------------------------------------------------------------------ ! Build sigma grid !> This routine builds a grid based on terrain-following coordinates. -subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) +subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface ) !------------------------------------------------------------------------------ ! This routine builds a grid based on terrain-following coordinates. ! The module parameter coordinateResolution(:) determines the resolution in @@ -1238,18 +1268,22 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] ! Local variables - integer :: i, j, k - integer :: nz - real :: nominalDepth, totalThickness + real :: nominalDepth ! The nominal depth of the sea-floor in thickness units [H ~> m or kg m-2] + real :: totalThickness ! The total thickness of the water column [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ - real :: dh + real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2] #endif real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2] real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2] + integer :: i, j, k, nz nz = GV%ke @@ -1261,28 +1295,35 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) cycle endif - ! The rest of the model defines grids integrating up from the bottom - nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H - ! Determine water column height totalThickness = 0.0 do k = 1,nz totalThickness = totalThickness + h(i,j,k) enddo + ! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that + ! cancels out when determining coordinate motion, so referencing the column postions to + ! the surface is perfectly acceptable, but for preservation of previous answers the + ! referencing is done relative to the bottom when in Boussinesq mode. + ! if (GV%Boussinesq) then + nominalDepth = nom_depth_H(i,j) + ! else + ! nominalDepth = totalThickness + ! endif + call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew) ! Calculate the final change in grid position after blending new and old grids zOld(nz+1) = -nominalDepth do k = nz,1,-1 - zOld(k) = zOld(k+1) + h(i, j, k) + zOld(k) = zOld(k+1) + h(i,j,k) enddo call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ - dh=max(nominalDepth,totalThickness) - if (abs(zNew(1)-zOld(1))>(CS%nk-1)*0.5*epsilon(dh)*dh) then + dh = max(nominalDepth,totalThickness) + if (abs(zNew(1)-zOld(1)) > (CS%nk-1)*0.5*epsilon(dh)*dh) then write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz,CS%nk @@ -1314,11 +1355,11 @@ end subroutine build_sigma_grid ! Build grid based on target interface densities !------------------------------------------------------------------------------ !> This routine builds a new grid based on a given set of target interface densities. -subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shelf_h ) +subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS, frac_shelf_h ) !------------------------------------------------------------------------------ ! This routine builds a new grid based on a given set of target interface ! densities (these target densities are computed by taking the mean value -! of given layer densities). The algorithn operates as follows within each +! of given layer densities). The algorithm operates as follows within each ! column: ! 1. Given T & S within each layer, the layer densities are computed. ! 2. Based on these layer densities, a global density profile is reconstructed @@ -1331,17 +1372,21 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel !------------------------------------------------------------------------------ ! Arguments - type(regridding_CS), intent(in) :: CS !< Regridding control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure - real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth - !! [H ~> m or kg m-2] - type(remapping_CS), intent(in) :: remapCS !< The remapping control structure - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice - !! shelf coverage [nondim] + type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth + !! [H ~> m or kg m-2] + type(remapping_CS), intent(in) :: remapCS !< The remapping control structure + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice + !! shelf coverage [nondim] ! Local variables integer :: nz ! The number of layers in the input grid integer :: i, j, k @@ -1351,7 +1396,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] real :: totalThickness ! Total thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ - real :: dh + real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2] #endif logical :: ice_shelf @@ -1378,15 +1423,22 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel cycle endif - - ! Local depth (G%bathyT is positive downward) - nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H - ! Determine total water column thickness totalThickness = 0.0 do k=1,nz totalThickness = totalThickness + h(i,j,k) enddo + + ! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that + ! cancels out when determining coordinate motion, so referencing the column postions to + ! the surface is perfectly acceptable, but for preservation of previous answers the + ! referencing is done relative to the bottom when in Boussinesq mode. + ! if (GV%Boussinesq) then + nominalDepth = nom_depth_H(i,j) + ! else + ! nominalDepth = totalThickness + ! endif + ! Determine absolute interface positions zOld(nz+1) = - nominalDepth do k = nz,1,-1 @@ -1394,13 +1446,13 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel enddo if (ice_shelf) then - call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & - tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, & - z_rigid_top = totalThickness - nominalDepth, eta_orig = zOld(1), & + call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i,j,:), & + tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, zNew, & + z_rigid_top=totalThickness - nominalDepth, eta_orig = zOld(1), & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) else - call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & - tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, & + call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i,j,:), & + tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, zNew, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) endif @@ -1441,8 +1493,8 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel totalThickness = totalThickness + h(i,j,k) enddo - dh=max(nominalDepth,totalThickness) - if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then + dh = max(nominalDepth, totalThickness) + if (abs(zNew(1)-zOld(1)) > (nz-1)*0.5*epsilon(dh)*dh) then write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness write(0,*) 'zNew(1)-zOld(1) = ',zNew(1)-zOld(1),epsilon(dh),nz @@ -1475,11 +1527,15 @@ 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, tv, h_new, dzInterface, remapCS, CS, frac_shelf_h ) +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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure type(remapping_CS), intent(in) :: remapCS !< The remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure @@ -1487,6 +1543,10 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, remapCS, CS, real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf !! coverage [nondim] + real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate + !! resolution in Z to desired units for zInterface, + !! usually Z_to_H in which case it is in + !! units of [H Z-1 ~> nondim or kg m-3] ! Local variables real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2] @@ -1517,12 +1577,12 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, remapCS, CS, do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then - nominalDepth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H + nominalDepth = nom_depth_H(i,j) if (ice_shelf) then totalThickness = 0.0 do k=1,GV%ke - totalThickness = totalThickness + h(i,j,k) * GV%Z_to_H + totalThickness = totalThickness + h(i,j,k) enddo z_top_col = max(nominalDepth-totalThickness,0.0) else @@ -1538,7 +1598,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, remapCS, CS, call build_hycom1_column(CS%hycom_CS, remapCS, tv%eqn_of_state, GV%ke, nominalDepth, & h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & - z_col, z_col_new, zScale=GV%Z_to_H, & + z_col, z_col_new, zScale=zScale, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) ! Calculate the final change in grid position after blending new and old grids @@ -1562,11 +1622,15 @@ end subroutine build_grid_HyCOM1 !> This subroutine builds an adaptive grid that follows density surfaces where !! possible, subject to constraints on the smoothness of interface heights. -subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) +subroutine build_grid_adaptive(G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables type(regridding_CS), intent(in) :: CS !< Regridding control structure @@ -1576,8 +1640,8 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) ! local variables integer :: i, j, k, nz ! indices and dimension lengths - ! temperature [C ~> degC], salinity [S ~> ppt] and pressure on interfaces - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt ! Temperature on interfaces [C ~> degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: sInt ! Salinity on interfaces [S ~> ppt] ! current interface positions and after tendency term is applied ! positive downward real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2] @@ -1614,7 +1678,7 @@ subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) cycle endif - call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, 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 @@ -1685,7 +1749,7 @@ end subroutine adjust_interface_motion !------------------------------------------------------------------------------ ! Build arbitrary grid !------------------------------------------------------------------------------ -subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) +subroutine build_grid_arbitrary( G, GV, h, nom_depth_H, dzInterface, h_new, CS ) !------------------------------------------------------------------------------ ! This routine builds a grid based on arbitrary rules !------------------------------------------------------------------------------ @@ -1694,6 +1758,10 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Original layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface !! depth [H ~> m or kg m-2] @@ -1718,7 +1786,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) do i = G%isc-1,G%iec+1 ! Local depth - local_depth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H + local_depth = nom_depth_H(i,j) ! Determine water column height total_height = 0.0 @@ -2373,7 +2441,7 @@ function getStaticThickness( CS, SSH, depth ) real, dimension(CS%nk) :: getStaticThickness !< The returned thicknesses in the units of depth ! Local integer :: k - real :: z, dz + 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, & @@ -2407,10 +2475,13 @@ end function getStaticThickness subroutine dz_function1( string, dz ) character(len=*), intent(in) :: string !< String with list of parameters in form !! dz_min, H_total, power, precision - real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses + real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses [m] or other units ! Local variables integer :: nk, k - real :: dz_min, power, prec, H_total + real :: dz_min ! minimum grid spacing [m] or other units + real :: power ! A power to raise the relative position in index space [nondim] + real :: prec ! The precision with which positions are returned [m] or other units + real :: H_total ! The sum of the nominal thicknesses [m] or other units nk = size(dz) ! Number of cells prec = -1024. diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 91df78c021..ee612788c9 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -112,7 +112,7 @@ subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoom if (present(adaptDoMin)) CS%adaptDoMin = adaptDoMin end subroutine set_adapt_params -subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNext) +subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, nom_depth_H, zNext) type(adapt_CS), intent(in) :: CS !< The control structure for this module type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure @@ -125,6 +125,10 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: tInt !< Interface temperatures [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: sInt !< Interface salinities [S ~> ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column + !! relative to mean sea level or another locally + !! valid reference height, converted to thickness + !! units [H ~> m or kg m-2] real, dimension(SZK_(GV)+1), intent(inout) :: zNext !< updated interface positions ! Local variables @@ -144,7 +148,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex zNext(nz+1) = zInt(i,j,nz+1) ! local depth for scaling diffusivity - depth = (G%bathyT(i,j) + G%Z_ref) * GV%Z_to_H + depth = nom_depth_H(i,j) ! initialize del2sigma and the thickness change response to it zero del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0 @@ -244,9 +248,9 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex ! set vertical grid diffusivity kGrid(k) = (CS%adaptTimeRatio * nz**2 * depth) * & - (CS%adaptZoomCoeff / (CS%adaptZoom + 0.5*(zNext(K) + zNext(K+1))) + & - (CS%adaptBuoyCoeff * drdz / CS%adaptDrho0) + & - max(1.0 - CS%adaptZoomCoeff - CS%adaptBuoyCoeff, 0.0) / depth) + ( CS%adaptZoomCoeff / (CS%adaptZoom + 0.5*(zNext(K) + zNext(K+1))) + & + (CS%adaptBuoyCoeff * drdz / CS%adaptDrho0) + & + max(1.0 - CS%adaptZoomCoeff - CS%adaptBuoyCoeff, 0.0) / depth) enddo ! initial denominator (first diagonal element) diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index aa2715eb42..ddc569e45e 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -117,7 +117,8 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ real, dimension(nz+1), intent(in) :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2] real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces [H ~> m or kg m-2] real, optional, intent(in) :: zScale !< Scaling factor from the input coordinate thicknesses in [Z ~> m] - !! to desired units for zInterface, perhaps GV%Z_to_H. + !! to desired units for zInterface, perhaps GV%Z_to_H in which + !! case this has units of [H Z-1 ~> nondim or kg m-3] real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of !! cell reconstruction [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of @@ -220,7 +221,6 @@ subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, depth, h, real, dimension(nz), intent(in) :: S !< Salinity of column [S ~> ppt] real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(nz), intent(in) :: p_col !< Layer pressure [R L2 T-2 ~> Pa] - !! to desired units for zInterface, perhaps GV%Z_to_H. real, dimension(nz), intent(out) :: R !< Layer density [R ~> kg m-3] real, dimension(nz+1), intent(out) :: RiAnom !< The interface density anomaly !! w.r.t. the interface target diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 8454c4be1d..7b6c0e0f8c 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -102,9 +102,9 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(CS%nk+1), & intent(inout) :: z_interface !< Absolute positions of interfaces real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same - !! units as depth) [Z ~> m] or [H ~> m or kg m-2] + !! units as depth) [H ~> m or kg m-2] real, optional, intent(in) :: eta_orig !< The actual original height of the top in the same - !! units as depth) [Z ~> m] or [H ~> m or kg m-2] + !! units as depth) [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose !! of cell reconstructions [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose @@ -119,22 +119,10 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(nz+1) :: xTmp ! Temporary positions [H ~> m or kg m-2] real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2] real, dimension(CS%nk+1) :: x1 ! Interface heights [H ~> m or kg m-2] - real :: z0_top, eta ! Thicknesses or heights [Z ~> m] or [H ~> m or kg m-2] ! Construct source column with vanished layers removed (stored in h_nv) call copy_finite_thicknesses(nz, h, CS%min_thickness, count_nonzero_layers, h_nv, mapping) - z0_top = 0. - eta=0.0 - if (present(z_rigid_top)) then - z0_top = z_rigid_top - eta=z0_top - if (present(eta_orig)) then - eta=eta_orig - endif - endif - - if (count_nonzero_layers > 1) then xTmp(1) = 0.0 do k = 1,count_nonzero_layers diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0321d7511a..802fd33d0f 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -464,7 +464,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (new_sim .and. debug) & call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, scale=GV%H_to_MKS) - call ALE_regrid_accelerated(ALE_CSp, G, GV, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & + call ALE_regrid_accelerated(ALE_CSp, G, GV, US, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & dt=dt, initial=.true.) endif endif @@ -2787,7 +2787,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just call regridding_preadjust_reqs(regridCS, do_conv_adj, ignore) if (do_conv_adj) call convective_adjustment(G, GV_loc, h1, tv_loc) - call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, & + call regridding_main( remapCS, regridCS, G, GV_loc, US, h1, tv_loc, h, dz_interface, & frac_shelf_h=frac_shelf_h ) deallocate( dz_interface ) diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 40dced9b20..0577a12ac5 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -1069,7 +1069,7 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale) call pass_var(h, G%Domain) call pass_var(CS%tv%T, G%Domain) call pass_var(CS%tv%S, G%Domain) - call ALE_offline_inputs(CS%ALE_CSp, G, GV, h, CS%tv, CS%tracer_Reg, CS%uhtr, CS%vhtr, CS%Kd, & + call ALE_offline_inputs(CS%ALE_CSp, G, GV, US, h, CS%tv, CS%tracer_Reg, CS%uhtr, CS%vhtr, CS%Kd, & CS%debug, CS%OBC) if (CS%id_temp_regrid>0) call post_data(CS%id_temp_regrid, CS%tv%T, CS%diag) if (CS%id_salt_regrid>0) call post_data(CS%id_salt_regrid, CS%tv%S, CS%diag)