diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index e23e740c9c..b1897aeb2e 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -84,7 +84,7 @@ module MOM_regridding !> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2]. real :: min_thickness - !> Reference pressure for potential density calculations (Pa) + !> Reference pressure for potential density calculations [Pa] real :: ref_pressure = 2.e7 !> Weight given to old coordinate when blending between new and old grids [nondim] @@ -203,6 +203,7 @@ 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_fixed_sfc, Rho_avg_depth, nlay_sfc_int real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha + real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3] integer :: nz_fixed_sfc, 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. @@ -574,26 +575,26 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, & - "Ratio of ALE timestep to grid timescale.", units="nondim", default=1e-1) + "Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1) call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, & - "Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H) + "Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H) call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, & - "Coefficient of near-surface zooming diffusivity.", & - units="nondim", default=0.2) + "Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2) call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptBuoyCoeff, & - "Coefficient of buoyancy diffusivity.", & - units="nondim", default=0.8) + "Coefficient of buoyancy diffusivity.", units="nondim", default=0.8) call get_param(param_file, mdl, "ADAPT_ALPHA", adaptAlpha, & - "Scaling on optimization tendency.", & - units="nondim", default=1.0) + "Scaling on optimization tendency.", units="nondim", default=1.0) call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmpLogical, & - "If true, make a HyCOM-like mixed layer by preventing interfaces "//& - "from being shallower than the depths specified by the regridding coordinate.", & - default=.false.) + "If true, make a HyCOM-like mixed layer by preventing interfaces "//& + "from being shallower than the depths specified by the regridding coordinate.", & + default=.false.) + call get_param(param_file, mdl, "ADAPT_DRHO0", adaptDrho0, & + "Reference density difference for stratification-dependent diffusion.", & + units="kg m-3", default=0.5, scale=US%kg_m3_to_R) call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, & adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, & - adaptDoMin=tmpLogical) + adaptDoMin=tmpLogical, adaptDrho0=US%R_to_kg_m3*adaptDrho0) endif if (main_parameters .and. coord_is_state_dependent) then @@ -1012,9 +1013,9 @@ end subroutine check_grid_column subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) type(regridding_CS), intent(in) :: CS !< Regridding control structure integer, intent(in) :: nk !< Number of cells in source grid - real, dimension(nk+1), intent(in) :: z_old !< Old grid position [m] - real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [m] - real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [m] + real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2] + real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [H ~> m or kg m-2] + real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [H ~> m or kg m-2] ! Local variables real :: sgn ! The sign convention for downward. real :: dz_tgt, zr1, z_old_k @@ -1156,26 +1157,21 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) 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), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2]. - real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage. + real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage [nondim]. ! Local variables - integer :: i, j, k - integer :: nz - real :: nominalDepth, totalThickness, dh - real, dimension(SZK_(GV)+1) :: zOld, zNew - real :: minThickness + real :: nominalDepth, totalThickness, dh ! Depths and thicknesses [H ~> m or kg m-2] + real, dimension(SZK_(GV)+1) :: zOld, zNew ! Coordinate interface heights [H ~> m or kg m-2] + integer :: i, j, k, nz logical :: ice_shelf nz = GV%ke - minThickness = CS%min_thickness ice_shelf = .false. if (present(frac_shelf_h)) then if (associated(frac_shelf_h)) ice_shelf = .true. endif -!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, & -!$OMP ice_shelf,minThickness) & -!$OMP private(nominalDepth,totalThickness, & -!$OMP zNew,dh,zOld) +!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h,ice_shelf) & +!$OMP private(nominalDepth,totalThickness,zNew,dh,zOld) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 @@ -1218,7 +1214,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) #ifdef __DO_SAFETY_CHECKS__ dh=max(nominalDepth,totalThickness) if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then - write(0,*) 'min_thickness=',minThickness + write(0,*) 'min_thickness=',CS%min_thickness write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz do k=1,nz+1 @@ -1350,10 +1346,11 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS ) ! Local variables integer :: nz integer :: i, j, k - real :: nominalDepth, totalThickness + real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: zOld, zNew - real :: h_neglect, h_neglect_edge + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ + real :: totalThickness real :: dh #endif @@ -1539,8 +1536,8 @@ subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt ! current interface positions and after tendency term is applied ! positive downward - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt - real, dimension(SZK_(GV)+1) :: zNext + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2] + real, dimension(SZK_(GV)+1) :: zNext ! New interface depths [H ~> m or kg m-2] nz = GV%ke @@ -2231,7 +2228,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, & nlay_ML_to_interior, fix_haloclines, halocline_filt_len, & halocline_strat_tol, integrate_downward_for_e, remap_answers_2018, & - adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin) + 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 @@ -2266,6 +2263,8 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri logical, optional, intent(in) :: adaptDoMin !< If true, make a HyCOM-like mixed layer by !! preventing interfaces from being shallower than !! the depths specified by the regridding coordinate. + real, optional, intent(in) :: adaptDrho0 !< Reference density difference for stratification-dependent + !! diffusion. [kg m-3] if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme) if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation) @@ -2322,6 +2321,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri if (present(adaptBuoyCoeff)) call set_adapt_params(CS%adapt_CS, adaptBuoyCoeff=adaptBuoyCoeff) if (present(adaptAlpha)) call set_adapt_params(CS%adapt_CS, adaptAlpha=adaptAlpha) if (present(adaptDoMin)) call set_adapt_params(CS%adapt_CS, adaptDoMin=adaptDoMin) + if (present(adaptDrho0)) call set_adapt_params(CS%adapt_CS, adaptDrho0=adaptDrho0) end select end subroutine set_regrid_params diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index 76c346c82e..1686ac51c9 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -15,7 +15,7 @@ module coord_hycom !> Number of layers/levels in generated grid integer :: nk - !> Nominal near-surface resolution + !> Nominal near-surface resolution [Z ~> m] real, allocatable, dimension(:) :: coordinateResolution !> Nominal density of interfaces [R ~> kg m-3] @@ -24,10 +24,10 @@ module coord_hycom !> Density scaling factor [R m3 kg-1 ~> 1] real :: kg_m3_to_R - !> Maximum depths of interfaces + !> Maximum depths of interfaces [H ~> m or kg m-2] real, allocatable, dimension(:) :: max_interface_depths - !> Maximum thicknesses of layers + !> Maximum thicknesses of layers [H ~> m or kg m-2] real, allocatable, dimension(:) :: max_layer_thickness !> Interpolation control structure @@ -42,7 +42,7 @@ module coord_hycom subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS, rho_scale) type(hycom_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in generated grid - real, dimension(nk), intent(in) :: coordinateResolution !< Nominal near-surface resolution [m] + real, dimension(nk), intent(in) :: coordinateResolution !< Nominal near-surface resolution [Z ~> m] real, dimension(nk+1),intent(in) :: target_density !< Interface target densities [R ~> kg m-3] type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density @@ -76,8 +76,8 @@ end subroutine end_coord_hycom !> This subroutine can be used to set the parameters for the coord_hycom module subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, interp_CS) type(hycom_CS), pointer :: CS !< Coordinate control structure - real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces in m - real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers in m + real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2] + real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2] type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation if (.not. associated(CS)) call MOM_error(FATAL, "set_hycom_params: CS not associated") @@ -102,33 +102,31 @@ end subroutine set_hycom_params !> Build a HyCOM coordinate column subroutine build_hycom1_column(CS, eqn_of_state, nz, 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(hycom_CS), intent(in) :: CS !< Coordinate control structure type(EOS_type), pointer :: eqn_of_state !< Equation of state structure - integer, intent(in) :: nz !< Number of levels + integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) - real, dimension(nz), intent(in) :: T !< Temperature of column [degC] - real, dimension(nz), intent(in) :: S !< Salinity of column [ppt] - real, dimension(nz), intent(in) :: h !< Layer thicknesses, in [m] or [H ~> m or kg m-2] + real, dimension(nz), intent(in) :: T !< Temperature of column [degC] + real, dimension(nz), intent(in) :: S !< Salinity of column [ppt] + real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(nz), intent(in) :: p_col !< Layer pressure [Pa] 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 - real, optional, intent(in) :: zScale !< Scaling factor from the input thicknesses in [m] - !! to desired units for zInterface, perhaps m_to_H. - real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h. - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of edge value calculations - !! in the same units as h0. + 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. + 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 + !! edge value calculation [H ~> m or kg m-2] ! Local variables integer :: k real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3] real, dimension(CS%nk) :: h_col_new ! New layer thicknesses - real :: z_scale - real :: stretching ! z* stretching, converts z* to z. - real :: nominal_z ! Nominal depth of interface when using z* [Z ~> m] - real :: hNew + real :: z_scale ! A scaling factor from the input thicknesses to the target thicknesses, + ! perhaps 1 or a factor in [H Z-1 ~> 1 or kg m-3] + real :: stretching ! z* stretching, converts z* to z [nondim]. + real :: nominal_z ! Nominal depth of interface when using z* [H ~> m or kg m-2] logical :: maximum_depths_set ! If true, the maximum depths of interface have been set. logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set. diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 53b83644af..a78b1dd749 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -36,11 +36,6 @@ module coord_rho type(interp_CS_type) :: interp_CS end type rho_CS -!> Maximum number of regridding iterations -integer, parameter :: NB_REGRIDDING_ITERATIONS = 1 -!> Deviation tolerance between succesive grids in regridding iterations -real, parameter :: DEVIATION_TOLERANCE = 1e-10 - public init_coord_rho, set_rho_params, build_rho_column, old_inflate_layers_1d, end_coord_rho contains @@ -50,7 +45,7 @@ subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS, rho_s type(rho_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in the grid real, intent(in) :: ref_pressure !< Coordinate reference pressure [Pa] - real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [kg m-3 or R ~> kg m-3] + real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3] type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density @@ -103,24 +98,24 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) real, intent(in) :: depth !< Depth of ocean bottom (positive in m) real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(nz), intent(in) :: T !< T for source column - real, dimension(nz), intent(in) :: S !< S for source column + real, dimension(nz), intent(in) :: T !< Temperature for source column [degC] + real, dimension(nz), intent(in) :: S !< Salinity for source column [ppt] type(EOS_type), pointer :: eqn_of_state !< Equation of state structure real, dimension(CS%nk+1), & intent(inout) :: z_interface !< Absolute positions of interfaces - real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h0. - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of edge value calculations - !! in the same units as h0. + 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 + !! of edge value calculations [H ~> m or kg m-2] + ! Local variables integer :: k, count_nonzero_layers integer, dimension(nz) :: mapping - real, dimension(nz) :: p, h_nv + real, dimension(nz) :: pres ! Pressures used to calculate density [Pa] + real, dimension(nz) :: h_nv ! Thicknesses of non-vanishing layers [H ~> m or kg m-2] real, dimension(nz) :: densities ! Layer density [R ~> kg m-3] - real, dimension(nz+1) :: xTmp - real, dimension(CS%nk) :: h_new ! New thicknesses + 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 ! Construct source column with vanished layers removed (stored in h_nv) @@ -133,8 +128,8 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & enddo ! Compute densities on source column - p(:) = CS%ref_pressure - call calculate_density(T, S, p, densities, 1, nz, eqn_of_state, scale=CS%kg_m3_to_R) + pres(:) = CS%ref_pressure + call calculate_density(T, S, pres, densities, 1, nz, eqn_of_state, scale=CS%kg_m3_to_R) do k = 1,count_nonzero_layers densities(k) = densities(mapping(k)) enddo @@ -179,6 +174,8 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & end subroutine build_rho_column +!### build_rho_column_iteratively is never used or called. + !> Iteratively build a rho coordinate column !! !! The algorithm operates as follows within each column: @@ -192,7 +189,7 @@ end subroutine build_rho_column !! 5. Return to step 1 until convergence or until the maximum number of !! iterations is reached, whichever comes first. subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, & - zInterface, h_neglect, h_neglect_edge) + zInterface, h_neglect, h_neglect_edge, dev_tol) type(rho_CS), intent(in) :: CS !< Regridding control structure type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options integer, intent(in) :: nz !< Number of levels @@ -208,29 +205,39 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ real, optional, intent(in) :: h_neglect_edge !< A negligibly small width !! for the purpose of edge value calculations !! in the same units as h [Z ~> m] + real, optional, intent(in) :: dev_tol !< The tolerance for the deviation between + !! successive grids for determining when the + !! iterative solver has converged [Z ~> m] + ! Local variables - integer :: k, m - integer :: count_nonzero_layers - real :: deviation ! When iterating to determine the final - ! grid, this is the deviation between two - ! successive grids. - real :: threshold - real, dimension(nz) :: p, densities, T_tmp, S_tmp, Tmp - integer, dimension(nz) :: mapping - real, dimension(nz) :: h0, h1, hTmp - real, dimension(nz+1) :: x0, x1, xTmp + real, dimension(nz+1) :: x0, x1, xTmp ! Temporary interface heights [Z ~> m] + real, dimension(nz) :: pres ! The pressure used in the equation of state [Pa]. + real, dimension(nz) :: densities ! Layer densities [R ~> kg m-3] + real, dimension(nz) :: T_tmp, S_tmp ! A temporary profile of temperature [degC] and salinity [ppt]. + real, dimension(nz) :: Tmp ! A temporary variable holding a remapped variable. + real, dimension(nz) :: h0, h1, hTmp ! Temporary thicknesses [Z ~> m] + real :: deviation ! When iterating to determine the final grid, this is the + ! deviation between two successive grids [Z ~> m]. + real :: deviation_tol ! Deviation tolerance between succesive grids in + ! regridding iterations [Z ~> m] + real :: threshold ! The minimum thickness for a layer to be considered to exist [Z ~> m] + integer, dimension(nz) :: mapping ! The indices of the massive layers in the initial column. + integer :: k, m, count_nonzero_layers + + ! Maximum number of regridding iterations + integer, parameter :: NB_REGRIDDING_ITERATIONS = 1 threshold = CS%min_thickness - p(:) = CS%ref_pressure + pres(:) = CS%ref_pressure T_tmp(:) = T(:) S_tmp(:) = S(:) h0(:) = h(:) ! Start iterations to build grid m = 1 - deviation = 1e10 - do while ( ( m <= NB_REGRIDDING_ITERATIONS ) .and. & - ( deviation > DEVIATION_TOLERANCE ) ) + deviation_tol = 1.0e-15*depth ; if (present(dev_tol)) deviation_tol = dev_tol + + do m=1,NB_REGRIDDING_ITERATIONS ! Construct column with vanished layers removed call copy_finite_thicknesses(nz, h0, threshold, count_nonzero_layers, hTmp, mapping) @@ -245,7 +252,7 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ enddo ! Compute densities within current water column - call calculate_density( T_tmp, S_tmp, p, densities, & + call calculate_density( T_tmp, S_tmp, pres, densities, & 1, nz, eqn_of_state, scale=CS%kg_m3_to_R) do k = 1,count_nonzero_layers @@ -282,11 +289,10 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ enddo deviation = sqrt( deviation / (nz-1) ) - m = m + 1 + if ( deviation <= deviation_tol ) exit ! Copy final grid onto start grid for next iteration h0(:) = h1(:) - enddo ! end regridding iterations if (CS%integrate_downward_for_e) then @@ -309,16 +315,18 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ end subroutine build_rho_column_iteratively !> Copy column thicknesses with vanished layers removed -subroutine copy_finite_thicknesses(nk, h_in, threshold, nout, h_out, mapping) - integer, intent(in) :: nk !< Number of layer for h_in, T_in, S_in - real, dimension(nk), intent(in) :: h_in !< Thickness of input column - real, intent(in) :: threshold !< Thickness threshold defining vanished layers - integer, intent(out) :: nout !< Number of non-vanished layers - real, dimension(nk), intent(out) :: h_out !< Thickness of output column +subroutine copy_finite_thicknesses(nk, h_in, thresh, nout, h_out, mapping) + integer, intent(in) :: nk !< Number of layer for h_in, T_in, S_in + real, dimension(nk), intent(in) :: h_in !< Thickness of input column [H ~> m or kg m-2] or [Z ~> m] + real, intent(in) :: thresh !< Thickness threshold defining vanished + !! layers [H ~> m or kg m-2] or [Z ~> m] + integer, intent(out) :: nout !< Number of non-vanished layers + real, dimension(nk), intent(out) :: h_out !< Thickness of output column [H ~> m or kg m-2] or [Z ~> m] integer, dimension(nk), intent(out) :: mapping !< Index of k-out corresponding to k-in ! Local variables integer :: k, k_thickest - real :: thickness_in_vanished, thickest_h_out + real :: thickness_in_vanished ! Summed thicknesses in discarded layers [H ~> m or kg m-2] or [Z ~> m] + real :: thickest_h_out ! Thickness of the thickest layer [H ~> m or kg m-2] or [Z ~> m] ! Build up new grid nout = 0 @@ -328,7 +336,7 @@ subroutine copy_finite_thicknesses(nk, h_in, threshold, nout, h_out, mapping) do k = 1, nk mapping(k) = nout ! Note k>=nout always h_out(k) = 0. ! Make sure h_out is set everywhere - if (h_in(k) > threshold) then + if (h_in(k) > thresh) then ! For non-vanished layers nout = nout + 1 mapping(nout) = k diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90 index 92de6e1ec3..30f2597090 100644 --- a/src/ALE/coord_slight.F90 +++ b/src/ALE/coord_slight.F90 @@ -76,7 +76,7 @@ subroutine init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_ type(slight_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in the grid real, intent(in) :: ref_pressure !< Coordinate reference pressure [Pa] - real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [kg m-3] + real, dimension(:), intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3] type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation real, optional, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density @@ -202,10 +202,10 @@ subroutine build_slight_column(CS, eqn_of_state, H_to_Pa, H_subroundoff, & real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose !! of edge value calculations [H ~> m or kg m-2]. ! Local variables - real, dimension(nz) :: rho_col ! Layer densities [R ~> kg m-3] - real, dimension(nz) :: T_f, S_f ! Filtered ayer quantities + real, dimension(nz) :: rho_col ! Layer densities [R ~> kg m-3] + real, dimension(nz) :: T_f, S_f ! Filtered layer temperature [degC] and salinity [ppt] logical, dimension(nz+1) :: reliable ! If true, this interface is in a reliable position. - real, dimension(nz+1) :: T_int, S_int ! Temperature and salinity interpolated to interfaces. + real, dimension(nz+1) :: T_int, S_int ! Temperature [degC] and salinity [ppt] interpolated to interfaces. real, dimension(nz+1) :: rho_tmp ! A temporary density [R ~> kg m-3] real, dimension(nz+1) :: drho_dp ! The partial derivative of density with pressure [kg m-3 Pa-1] real, dimension(nz+1) :: p_IS, p_R @@ -224,7 +224,7 @@ subroutine build_slight_column(CS, eqn_of_state, H_to_Pa, H_subroundoff, & real :: z_int_unst real :: dz ! A uniform layer thickness in very shallow water [H ~> m or kg m-2]. real :: dz_ur ! The total thickness of an unstable region [H ~> m or kg m-2]. - real :: wgt, cowgt ! A weight and its complement, nondim. + real :: wgt, cowgt ! A weight and its complement [nondim]. real :: rho_ml_av ! The average potential density in a near-surface region [R ~> kg m-3]. real :: H_ml_av ! A thickness to try to use in taking the near-surface average [H ~> m or kg m-2]. real :: rho_x_z ! A cumulative integral of a density [R H ~> kg m-2 or kg2 m-5]. @@ -492,38 +492,38 @@ end subroutine build_slight_column subroutine rho_interfaces_col(rho_col, h_col, z_col, rho_tgt, nz, z_col_new, & CS, reliable, debug, h_neglect, h_neglect_edge) integer, intent(in) :: nz !< Number of layers - real, dimension(nz), intent(in) :: rho_col !< Initial layer reference densities. - real, dimension(nz), intent(in) :: h_col !< Initial layer thicknesses. - real, dimension(nz+1), intent(in) :: z_col !< Initial interface heights. + real, dimension(nz), intent(in) :: rho_col !< Initial layer reference densities [R ~> kg m-3]. + real, dimension(nz), intent(in) :: h_col !< Initial layer thicknesses [H ~> m or kg m-2]. + real, dimension(nz+1), intent(in) :: z_col !< Initial interface heights [H ~> m or kg m-2]. real, dimension(nz+1), intent(in) :: rho_tgt !< Interface target densities. - real, dimension(nz+1), intent(inout) :: z_col_new !< New interface heights. + real, dimension(nz+1), intent(inout) :: z_col_new !< New interface heights [H ~> m or kg m-2]. type(slight_CS), intent(in) :: CS !< Coordinate control structure logical, dimension(nz+1), intent(inout) :: reliable !< If true, the interface positions !! are well defined from a stable region. - logical, optional, intent(in) :: debug !< If present and true, do debugging checks. - real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h_col. - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of edge value calculations - !! in the same units as h_col. + logical, optional, intent(in) :: debug !< If present and true, do debugging checks. + 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 + !! of edge value calculations [H ~> m or kg m-2] real, dimension(nz+1) :: ru_max_int ! The maximum and minimum densities in - real, dimension(nz+1) :: ru_min_int ! an unstable region around an interface. + real, dimension(nz+1) :: ru_min_int ! an unstable region around an interface [R ~> kg m-3]. real, dimension(nz) :: ru_max_lay ! The maximum and minimum densities in - real, dimension(nz) :: ru_min_lay ! an unstable region containing a layer. - real, dimension(nz,2) :: ppoly_i_E ! Edge value of polynomial - real, dimension(nz,2) :: ppoly_i_S ! Edge slope of polynomial - real, dimension(nz,DEGREE_MAX+1) :: ppoly_i_coefficients ! Coefficients of polynomial + real, dimension(nz) :: ru_min_lay ! an unstable region containing a layer [R ~> kg m-3]. + real, dimension(nz,2) :: ppoly_i_E ! Edge value of polynomial [R ~> kg m-3] + real, dimension(nz,2) :: ppoly_i_S ! Edge slope of polynomial [R H-1 ~> kg m-4 or m-1] + real, dimension(nz,DEGREE_MAX+1) :: ppoly_i_coefficients ! Coefficients of polynomial [R ~> kg m-3] logical, dimension(nz) :: unstable_lay ! If true, this layer is in an unstable region. logical, dimension(nz+1) :: unstable_int ! If true, this interface is in an unstable region. - real :: rt ! The current target density [kg m-3]. - real :: zf ! The fractional z-position within a layer of the target density. - real :: rfn - real :: a(5) ! Coefficients of a local polynomial minus the target density. - real :: zf1, zf2, rfn1, rfn2 - real :: drfn_dzf, sgn, delta_zf, zf_prev - real :: tol + real :: rt ! The current target density [R ~> kg m-3]. + real :: zf ! The fractional z-position within a layer of the target density [nondim]. + real :: rfn ! The target density relative to the interpolated density [R ~> kg m-3] + real :: a(5) ! Coefficients of a local polynomial minus the target density [R ~> kg m-3]. + real :: zf1, zf2 ! Two previous estimates of zf [nondim] + real :: rfn1, rfn2 ! Values of rfn at zf1 and zf2 [R ~> kg m-3] + real :: drfn_dzf ! The partial derivative of rfn with zf [R ~> kg m-3] + real :: sgn, delta_zf, zf_prev ! [nondim] + real :: tol ! The tolerance for convergence of zf [nondim] logical :: k_found ! If true, the position has been found. integer :: k_layer ! The index of the stable layer containing an interface. integer :: ppoly_degree diff --git a/src/ALE/coord_zlike.F90 b/src/ALE/coord_zlike.F90 index 1f4949431d..f2ed7f0035 100644 --- a/src/ALE/coord_zlike.F90 +++ b/src/ALE/coord_zlike.F90 @@ -63,17 +63,21 @@ end subroutine set_zlike_params subroutine build_zstar_column(CS, depth, total_thickness, zInterface, & z_rigid_top, eta_orig, zScale) type(zlike_CS), intent(in) :: CS !< Coordinate control structure - real, intent(in) :: depth !< Depth of ocean bottom (positive in the output units) - real, intent(in) :: total_thickness !< Column thickness (positive in the same units as depth) + real, intent(in) :: depth !< Depth of ocean bottom (positive downward in the + !! output units), units may be [Z ~> m] or [H ~> m or kg m-2] + real, intent(in) :: total_thickness !< Column thickness (positive definite in the same + !! units as depth) [Z ~> m] or [H ~> m or kg m-2] real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces - real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (negative in the - !! same units as depth) - real, optional, intent(in) :: eta_orig !< The actual original height of the top in the - !! same units as depth + 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] + 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] real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate resolution !! in Z to desired units for zInterface, perhaps Z_to_H ! Local variables - real :: eta, stretching, dh, min_thickness, z0_top, z_star, z_scale + real :: eta ! Free surface height [Z ~> m] or [H ~> m or kg m-2] + real :: stretching ! A stretching factor for the coordinate [nondim] + real :: dh, min_thickness, z0_top, z_star, z_scale ! Thicknesses or heights [Z ~> m] or [H ~> m or kg m-2] integer :: k logical :: new_zstar_def diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 19082292be..5a1d151487 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -58,13 +58,13 @@ module regrid_interp !> When the N-R algorithm produces an estimate that lies outside [0,1], the !! estimate is set to be equal to the boundary location, 0 or 1, plus or minus -!! an offset, respectively, when the derivative is zero at the boundary. +!! an offset, respectively, when the derivative is zero at the boundary [nondim]. real, public, parameter :: NR_OFFSET = 1e-6 !> Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are !! used to build the new grid by finding the coordinates associated with !! target densities and interpolations of degree larger than 1. integer, public, parameter :: NR_ITERATIONS = 8 -!> Tolerance for Newton-Raphson iterations (stop when increment falls below this) +!> Tolerance for Newton-Raphson iterations (stop when increment falls below this) [nondim] real, public, parameter :: NR_TOLERANCE = 1e-12 contains @@ -79,17 +79,17 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & ppoly0_coefs, degree, h_neglect, h_neglect_edge) type(interp_CS_type), intent(in) :: CS !< Interpolation control structure integer, intent(in) :: n0 !< Number of cells on source grid - real, dimension(n0), intent(in) :: densities !< Actual cell densities - real, dimension(n0), intent(in) :: h0 !< cell widths on source grid - real, dimension(n0,2), intent(inout) :: ppoly0_E !< Edge value of polynomial - real, dimension(n0,2), intent(inout) :: ppoly0_S !< Edge slope of polynomial - real, dimension(n0,DEGREE_MAX+1), intent(inout) :: ppoly0_coefs !< Coefficients of polynomial + real, dimension(n0), intent(in) :: densities !< Actual cell densities [A] + real, dimension(n0), intent(in) :: h0 !< cell widths on source grid [H] + real, dimension(n0,2), intent(inout) :: ppoly0_E !< Edge value of polynomial [A] + real, dimension(n0,2), intent(inout) :: ppoly0_S !< Edge slope of polynomial [A H-1] + real, dimension(n0,DEGREE_MAX+1), intent(inout) :: ppoly0_coefs !< Coefficients of polynomial [A] integer, intent(inout) :: degree !< The degree of the polynomials real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions + !! purpose of cell reconstructions [H] !! in the same units as h0. real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of edge value calculations + !! for the purpose of edge value calculations [H] !! in the same units as h0. ! Local variables logical :: extrapolate @@ -271,15 +271,15 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, & target_values, degree, n1, h1, x1, answers_2018 ) integer, intent(in) :: n0 !< Number of points on source grid integer, intent(in) :: n1 !< Number of points on target grid - real, dimension(n0), intent(in) :: h0 !< Thicknesses of source grid cells - real, dimension(n0+1), intent(in) :: x0 !< Source interface positions - real, dimension(n0,2), intent(in) :: ppoly0_E !< Edge values of interpolating polynomials + real, dimension(n0), intent(in) :: h0 !< Thicknesses of source grid cells [H] + real, dimension(n0+1), intent(in) :: x0 !< Source interface positions [H] + real, dimension(n0,2), intent(in) :: ppoly0_E !< Edge values of interpolating polynomials [A] real, dimension(n0,DEGREE_MAX+1), & - intent(in) :: ppoly0_coefs !< Coefficients of interpolating polynomials - real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces + intent(in) :: ppoly0_coefs !< Coefficients of interpolating polynomials [A] + real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [A] integer, intent(in) :: degree !< Degree of interpolating polynomials - real, dimension(n1), intent(inout) :: h1 !< Thicknesses of target grid cells - real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions + real, dimension(n1), intent(inout) :: h1 !< Thicknesses of target grid cells [H] + real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions [H] logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. ! Local variables @@ -309,21 +309,22 @@ subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, type(interp_CS_type), intent(in) :: CS !< A control structure for regrid_interp integer, intent(in) :: n0 !< The number of points on the input grid integer, intent(in) :: n1 !< The number of points on the output grid - real, dimension(n0), intent(in) :: densities !< Input cell densities [kg m-3] - real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces - real, dimension(n0), intent(in) :: h0 !< Initial cell widths - real, dimension(n0+1), intent(in) :: x0 !< Source interface positions - real, dimension(n1), intent(inout) :: h1 !< Output cell widths - real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions + real, dimension(n0), intent(in) :: densities !< Input cell densities [R ~> kg m-3] + real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [R ~> kg m-3] + real, dimension(n0), intent(in) :: h0 !< Initial cell widths [H] + real, dimension(n0+1), intent(in) :: x0 !< Source interface positions [H] + real, dimension(n1), intent(inout) :: h1 !< Output cell widths [H] + real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions [H] real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions + !! purpose of cell reconstructions [H] !! in the same units as h0. real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of edge value calculations + !! for the purpose of edge value calculations [H] !! in the same units as h0. - real, dimension(n0,2) :: ppoly0_E, ppoly0_S - real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C + real, dimension(n0,2) :: ppoly0_E ! Polynomial edge values [R ~> kg m-3] + real, dimension(n0,2) :: ppoly0_S ! Polynomial edge slopes [R H-1] + real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3] integer :: degree call regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, ppoly0_C, & @@ -352,28 +353,28 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & target_value, degree, answers_2018 ) result ( x_tgt ) ! Arguments integer, intent(in) :: N !< Number of grid cells - real, dimension(N), intent(in) :: h !< Grid cell thicknesses - real, dimension(N+1), intent(in) :: x_g !< Grid interface locations - real, dimension(N,2), intent(in) :: ppoly_E !< Edge values of interpolating polynomials - real, dimension(N,DEGREE_MAX+1), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials - real, intent(in) :: target_value !< Target value to find position for + real, dimension(N), intent(in) :: h !< Grid cell thicknesses [H] + real, dimension(N+1), intent(in) :: x_g !< Grid interface locations [H] + real, dimension(N,2), intent(in) :: ppoly_E !< Edge values of interpolating polynomials [A] + real, dimension(N,DEGREE_MAX+1), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials [A] + real, intent(in) :: target_value !< Target value to find position for [A] integer, intent(in) :: degree !< Degree of the interpolating polynomials logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. - real :: x_tgt !< The position of x_g at which target_value is found. + real :: x_tgt !< The position of x_g at which target_value is found [H] + ! Local variables - integer :: i, k ! loop indices - integer :: k_found ! index of target cell - integer :: iter - real :: xi0 ! normalized target coordinate - real, dimension(DEGREE_MAX) :: a ! polynomial coefficients + real :: xi0 ! normalized target coordinate [nondim] + real, dimension(DEGREE_MAX) :: a ! polynomial coefficients [A] real :: numerator real :: denominator - real :: delta ! Newton-Raphson increment - real :: x ! global target coordinate - real :: eps ! offset used to get away from - ! boundaries - real :: grad ! gradient during N-R iterations - logical :: use_2018_answers ! If true use older, less acccurate expressions. + real :: delta ! Newton-Raphson increment [nondim] +! real :: x ! global target coordinate + real :: eps ! offset used to get away from boundaries [nondim] + real :: grad ! gradient during N-R iterations [A] + integer :: i, k, iter ! loop indices + integer :: k_found ! index of target cell + character(len=200) :: mesg + logical :: use_2018_answers ! If true use older, less acccurate expressions. eps = NR_OFFSET k_found = -1 @@ -390,11 +391,9 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & ! Since discontinuous edge values are allowed, we check whether the target ! value lies between two discontinuous edge values at interior interfaces do k = 2,N - if ( ( target_value >= ppoly_E(k-1,2) ) .AND. & - ( target_value <= ppoly_E(k,1) ) ) then + if ( ( target_value >= ppoly_E(k-1,2) ) .AND. ( target_value <= ppoly_E(k,1) ) ) then x_tgt = x_g(k) return ! return because there is no need to look further - exit endif enddo @@ -412,8 +411,7 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & ! contains the target value. The variable k_found holds the index value ! of the cell where the taregt value lies. do k = 1,N - if ( ( target_value > ppoly_E(k,1) ) .AND. & - ( target_value < ppoly_E(k,2) ) ) then + if ( ( target_value > ppoly_E(k,1) ) .AND. ( target_value < ppoly_E(k,2) ) ) then k_found = k exit endif @@ -425,12 +423,10 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & ! means there is a major problem with the interpolant. This needs to be ! reported. if ( k_found == -1 ) then - write(*,*) target_value, ppoly_E(1,1), ppoly_E(N,2) - write(*,*) 'Could not find target coordinate in ' //& - '"get_polynomial_coordinate". This is caused by an '//& - 'inconsistent interpolant (perhaps not monotonically '//& - 'increasing)' - call MOM_error( FATAL, 'Aborting execution' ) + write(mesg,*) 'Could not find target coordinate', target_value, 'in get_polynomial_coordinate. This is '//& + 'caused by an inconsistent interpolant (perhaps not monotonically increasing):', & + target_value, ppoly_E(1,1), ppoly_E(N,2) + call MOM_error( FATAL, mesg ) endif ! Reset all polynomial coefficients to 0 and copy those pertaining to @@ -440,18 +436,11 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & a(i) = ppoly_coefs(k_found,i) enddo - ! Guess value to start Newton-Raphson iterations (middle of cell) + ! Guess the middle of the cell to start Newton-Raphson iterations xi0 = 0.5 - iter = 1 - delta = 1e10 ! Newton-Raphson iterations - do - ! break if converged or too many iterations taken - if ( ( iter > NR_ITERATIONS ) .OR. & - ( abs(delta) < NR_TOLERANCE ) ) then - exit - endif + do iter = 1,NR_ITERATIONS if (use_2018_answers) then numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & @@ -487,7 +476,8 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & if ( grad == 0.0 ) xi0 = xi0 - eps endif - iter = iter + 1 + ! break if converged or too many iterations taken + if ( abs(delta) < NR_TOLERANCE ) exit enddo ! end Newton-Raphson iterations x_tgt = x_g(k_found) + xi0 * h(k_found) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 421c23cf68..ceb782ce4b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -145,7 +145,8 @@ module MOM_diag_mediator !> Stores all the remapping grids and the model's native space thicknesses type, public :: diag_grid_storage integer :: num_diag_coords !< Number of target coordinates - real, dimension(:,:,:), allocatable :: h_state !< Layer thicknesses in native space + real, dimension(:,:,:), allocatable :: h_state !< Layer thicknesses in native + !! space [H ~> m or kg m-2] type(diag_grids_type), dimension(:), allocatable :: diag_grids !< Primarily empty, except h field end type diag_grid_storage @@ -312,9 +313,9 @@ module MOM_diag_mediator !>@} ! Pointer to H, G and T&S needed for remapping - real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping - real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping - real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping + real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping [H ~> m or kg m-2] + real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping [degC] + real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [ppt] type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type type(ocean_grid_type), pointer :: G => null() !< The ocean grid type type(verticalGrid_type), pointer :: GV => null() !< The model's vertical ocean grid @@ -324,7 +325,7 @@ module MOM_diag_mediator integer :: volume_cell_measure_dm_id = -1 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - ! Keep a copy of h so that we know whether it has changed. If it has then + ! Keep a copy of h so that we know whether it has changed [H ~> m or kg m-2]. If it has then ! need the target grid for vertical remapping needs to have been updated. real, dimension(:,:,:), allocatable :: h_old #endif @@ -1525,8 +1526,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) - call diag_remap_do_remap(diag_cs%diag_remap_cs( & - diag%axes%vertical_coordinate_number), & + call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & diag_cs%G, diag_cs%GV, h_diag, staggered_in_x, staggered_in_y, & diag%axes%mask3d, diag_cs%missing_value, field, remapped_field) if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) @@ -3202,7 +3202,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) !! the current salinity ! Local variables integer :: i - real, dimension(:,:,:), pointer :: h_diag => NULL() + real, dimension(:,:,:), pointer :: h_diag => NULL() ! The layer thickneses for diagnostics [H ~> m or kg m-2] real, dimension(:,:,:), pointer :: T_diag => NULL(), S_diag => NULL() if (present(alt_h)) then @@ -3231,9 +3231,8 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) endif do i=1, diag_cs%num_diag_coords - call diag_remap_update(diag_cs%diag_remap_cs(i), & - diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state) + call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, & + h_diag, T_diag, S_diag, diag_cs%eqn_of_state) enddo #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) @@ -3518,7 +3517,7 @@ end subroutine diag_grid_storage_init !> Copy from the main diagnostic arrays to the grid storage as well as the native thicknesses subroutine diag_copy_diag_to_storage(grid_storage, h_state, diag) type(diag_grid_storage), intent(inout) :: grid_storage !< Structure containing a snapshot of the target grids - real, dimension(:,:,:), intent(in) :: h_state !< Current model thicknesses + real, dimension(:,:,:), intent(in) :: h_state !< Current model thicknesses [H ~> m or kg m-2] type(diag_ctrl), intent(in) :: diag !< Diagnostic control structure used as the contructor integer :: m diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index b61c10eb7e..cadd74950a 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -276,14 +276,14 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) type(ocean_grid_type), pointer :: 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 - real, dimension(:, :, :), intent(in) :: h !< New thickness - real, dimension(:, :, :), intent(in) :: T !< New T - real, dimension(:, :, :), intent(in) :: S !< New S + real, dimension(:, :, :), intent(in) :: h !< New thickness [H ~> m or kg m-2] + real, dimension(:, :, :), intent(in) :: T !< New temperatures [degC] + real, dimension(:, :, :), intent(in) :: S !< New salinities [ppt] type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state ! Local variables - real, dimension(remap_cs%nz + 1) :: zInterfaces - real :: h_neglect, h_neglect_edge + real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] integer :: i, j, k, nz ! Note that coordinateMode('LAYER') is never 'configured' so will @@ -326,16 +326,17 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then +!### I think that the conversion factor in the 2nd line should be GV%Z_to_H call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, & US%Z_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then ! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & -! US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!") elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then ! call build_hycom1_column(remap_cs%regrid_cs, nz, & -! US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") endif remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1)