Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 33 additions & 33 deletions src/ALE/MOM_regridding.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
46 changes: 22 additions & 24 deletions src/ALE/coord_hycom.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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")
Expand All @@ -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.

Expand Down
Loading