MOM6
coord_adapt Module Reference

Detailed Description

Regrid columns for the adaptive coordinate.

Data Types

type  adapt_cs
 Control structure for adaptive coordinates (coord_adapt). More...
 

Functions/Subroutines

subroutine, public init_coord_adapt (CS, nk, coordinateResolution, m_to_H)
 Initialise an adapt_CS with parameters. More...
 
subroutine, public end_coord_adapt (CS)
 Clean up the coordinate control structure. More...
 
subroutine, public set_adapt_params (CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptDrho0, adaptDoMin)
 This subtroutine can be used to set the parameters for coord_adapt module. More...
 
subroutine, public build_adapt_column (CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
 

Function/Subroutine Documentation

◆ build_adapt_column()

subroutine, public coord_adapt::build_adapt_column ( type(adapt_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(thermo_var_ptrs), intent(in)  tv,
integer, intent(in)  i,
integer, intent(in)  j,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)  zInt,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)  tInt,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke+1), intent(in)  sInt,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  h,
real, dimension( gv %ke+1), intent(inout)  zNext 
)
Parameters
[in]csThe control structure for this module
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]tvA structure pointing to various thermodynamic variables
[in]iThe i-index of the column to work on
[in]jThe j-index of the column to work on
[in]zintInterface heights [H ~> m or kg m-2].
[in]tintInterface temperatures [degC]
[in]sintInterface salinities [ppt]
[in]hLayer thicknesses [H ~> m or kg m-2]
[in,out]znextupdated interface positions

Definition at line 118 of file coord_adapt.F90.

118  type(adapt_CS), intent(in) :: CS !< The control structure for this module
119  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
120  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
121  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
122  !! thermodynamic variables
123  integer, intent(in) :: i !< The i-index of the column to work on
124  integer, intent(in) :: j !< The j-index of the column to work on
125  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: zInt !< Interface heights [H ~> m or kg m-2].
126  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: tInt !< Interface temperatures [degC]
127  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: sInt !< Interface salinities [ppt]
128  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
129  real, dimension(SZK_(GV)+1), intent(inout) :: zNext !< updated interface positions
130 
131  ! Local variables
132  integer :: k, nz
133  real :: h_up, b1, b_denom_1, d1, depth, drdz, nominal_z, stretching
134  real, dimension(SZK_(GV)+1) :: alpha, beta, del2sigma ! drho/dT and drho/dS
135  real, dimension(SZK_(GV)) :: kGrid, c1 ! grid diffusivity on layers, and tridiagonal work array
136 
137  nz = cs%nk
138 
139  ! set bottom and surface of zNext
140  znext(1) = 0.
141  znext(nz+1) = zint(i,j,nz+1)
142 
143  ! local depth for scaling diffusivity
144  depth = g%bathyT(i,j) * gv%Z_to_H
145 
146  ! initialize del2sigma to zero
147  del2sigma(:) = 0.
148 
149  ! calculate del-squared of neutral density by a
150  ! stencilled finite difference
151  ! TODO: this needs to be adjusted to account for vanished layers near topography
152 
153  ! up (j-1)
154  if (g%mask2dT(i,j-1) > 0.) then
155  call calculate_density_derivs( &
156  0.5 * (tint(i,j,2:nz) + tint(i,j-1,2:nz)), &
157  0.5 * (sint(i,j,2:nz) + sint(i,j-1,2:nz)), &
158  0.5 * (zint(i,j,2:nz) + zint(i,j-1,2:nz)) * gv%H_to_Pa, &
159  alpha, beta, 2, nz - 1, tv%eqn_of_state)
160 
161  del2sigma(2:nz) = del2sigma(2:nz) + &
162  (alpha(2:nz) * (tint(i,j-1,2:nz) - tint(i,j,2:nz)) + &
163  beta(2:nz) * (sint(i,j-1,2:nz) - sint(i,j,2:nz)))
164  endif
165  ! down (j+1)
166  if (g%mask2dT(i,j+1) > 0.) then
167  call calculate_density_derivs( &
168  0.5 * (tint(i,j,2:nz) + tint(i,j+1,2:nz)), &
169  0.5 * (sint(i,j,2:nz) + sint(i,j+1,2:nz)), &
170  0.5 * (zint(i,j,2:nz) + zint(i,j+1,2:nz)) * gv%H_to_Pa, &
171  alpha, beta, 2, nz - 1, tv%eqn_of_state)
172 
173  del2sigma(2:nz) = del2sigma(2:nz) + &
174  (alpha(2:nz) * (tint(i,j+1,2:nz) - tint(i,j,2:nz)) + &
175  beta(2:nz) * (sint(i,j+1,2:nz) - sint(i,j,2:nz)))
176  endif
177  ! left (i-1)
178  if (g%mask2dT(i-1,j) > 0.) then
179  call calculate_density_derivs( &
180  0.5 * (tint(i,j,2:nz) + tint(i-1,j,2:nz)), &
181  0.5 * (sint(i,j,2:nz) + sint(i-1,j,2:nz)), &
182  0.5 * (zint(i,j,2:nz) + zint(i-1,j,2:nz)) * gv%H_to_Pa, &
183  alpha, beta, 2, nz - 1, tv%eqn_of_state)
184 
185  del2sigma(2:nz) = del2sigma(2:nz) + &
186  (alpha(2:nz) * (tint(i-1,j,2:nz) - tint(i,j,2:nz)) + &
187  beta(2:nz) * (sint(i-1,j,2:nz) - sint(i,j,2:nz)))
188  endif
189  ! right (i+1)
190  if (g%mask2dT(i+1,j) > 0.) then
191  call calculate_density_derivs( &
192  0.5 * (tint(i,j,2:nz) + tint(i+1,j,2:nz)), &
193  0.5 * (sint(i,j,2:nz) + sint(i+1,j,2:nz)), &
194  0.5 * (zint(i,j,2:nz) + zint(i+1,j,2:nz)) * gv%H_to_Pa, &
195  alpha, beta, 2, nz - 1, tv%eqn_of_state)
196 
197  del2sigma(2:nz) = del2sigma(2:nz) + &
198  (alpha(2:nz) * (tint(i+1,j,2:nz) - tint(i,j,2:nz)) + &
199  beta(2:nz) * (sint(i+1,j,2:nz) - sint(i,j,2:nz)))
200  endif
201 
202  ! at this point, del2sigma contains the local neutral density curvature at
203  ! h-points, on interfaces
204  ! we need to divide by drho/dz to give an interfacial displacement
205  !
206  ! a positive curvature means we're too light relative to adjacent columns,
207  ! so del2sigma needs to be positive too (push the interface deeper)
208  call calculate_density_derivs(tint(i,j,:), sint(i,j,:), zint(i,j,:) * gv%H_to_Pa, &
209  alpha, beta, 1, nz + 1, tv%eqn_of_state)
210  do k = 2, nz
211  ! TODO make lower bound here configurable
212  del2sigma(k) = del2sigma(k) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &
213  max(alpha(k) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + &
214  beta(k) * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20)
215 
216  ! don't move the interface so far that it would tangle with another
217  ! interface in the direction we're moving (or exceed a Nyquist limit
218  ! that could cause oscillations of the interface)
219  h_up = merge(h(i,j,k), h(i,j,k-1), del2sigma(k) > 0.)
220  del2sigma(k) = 0.5 * cs%adaptAlpha * &
221  sign(min(abs(del2sigma(k)), 0.5 * h_up), del2sigma(k))
222 
223  ! update interface positions so we can diffuse them
224  znext(k) = zint(i,j,k) + del2sigma(k)
225  enddo
226 
227  ! solve diffusivity equation to smooth grid
228  ! upper diagonal coefficients: -kGrid(2:nz)
229  ! lower diagonal coefficients: -kGrid(1:nz-1)
230  ! diagonal coefficients: 1 + (kGrid(1:nz-1) + kGrid(2:nz))
231  !
232  ! first, calculate the diffusivities within layers
233  do k = 1, nz
234  ! calculate the dr bit of drdz
235  drdz = 0.5 * (alpha(k) + alpha(k+1)) * (tint(i,j,k+1) - tint(i,j,k)) + &
236  0.5 * (beta(k) + beta(k+1)) * (sint(i,j,k+1) - sint(i,j,k))
237  ! divide by dz from the new interface positions
238  drdz = drdz / (znext(k) - znext(k+1) + gv%H_subroundoff)
239  ! don't do weird stuff in unstably-stratified regions
240  drdz = max(drdz, 0.)
241 
242  ! set vertical grid diffusivity
243  kgrid(k) = (cs%adaptTimeRatio * nz**2 * depth) * &
244  (cs%adaptZoomCoeff / (cs%adaptZoom + 0.5*(znext(k) + znext(k+1))) + &
245  (cs%adaptBuoyCoeff * drdz / cs%adaptDrho0) + &
246  max(1.0 - cs%adaptZoomCoeff - cs%adaptBuoyCoeff, 0.0) / depth)
247  enddo
248 
249  ! initial denominator (first diagonal element)
250  b1 = 1.0
251  ! initial Q_1 = 1 - q_1 = 1 - 0/1
252  d1 = 1.0
253  ! work on all interior interfaces
254  do k = 2, nz
255  ! calculate numerator of Q_k
256  b_denom_1 = 1. + d1 * kgrid(k-1)
257  ! update denominator for k
258  b1 = 1.0 / (b_denom_1 + kgrid(k))
259 
260  ! calculate q_k
261  c1(k) = kgrid(k) * b1
262  ! update Q_k = 1 - q_k
263  d1 = b_denom_1 * b1
264 
265  ! update RHS
266  znext(k) = b1 * (znext(k) + kgrid(k-1)*znext(k-1))
267  enddo
268  ! final substitution
269  do k = nz, 2, -1
270  znext(k) = znext(k) + c1(k)*znext(k+1)
271  enddo
272 
273  if (cs%adaptDoMin) then
274  nominal_z = 0.
275  stretching = zint(i,j,nz+1) / depth
276 
277  do k = 2, nz+1
278  nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
279  ! take the deeper of the calculated and nominal positions
280  znext(k) = max(znext(k), nominal_z)
281  ! interface can't go below topography
282  znext(k) = min(znext(k), zint(i,j,nz+1))
283  enddo
284  endif

Referenced by mom_regridding::build_grid_adaptive().

Here is the caller graph for this function:

◆ end_coord_adapt()

subroutine, public coord_adapt::end_coord_adapt ( type(adapt_cs), pointer  CS)

Clean up the coordinate control structure.

Parameters
csThe control structure for this module

Definition at line 82 of file coord_adapt.F90.

82  type(adapt_CS), pointer :: CS !< The control structure for this module
83 
84  ! nothing to do
85  if (.not. associated(cs)) return
86  deallocate(cs%coordinateResolution)
87  deallocate(cs)

Referenced by mom_regridding::end_regridding().

Here is the caller graph for this function:

◆ init_coord_adapt()

subroutine, public coord_adapt::init_coord_adapt ( type(adapt_cs), pointer  CS,
integer, intent(in)  nk,
real, dimension(:), intent(in)  coordinateResolution,
real, intent(in), optional  m_to_H 
)

Initialise an adapt_CS with parameters.

Parameters
csUnassociated pointer to hold the control structure
[in]nkNumber of layers in the grid
[in]coordinateresolutionNominal near-surface resolution [m] or other units specified with m_to_H
[in]m_to_hA conversion factor from m to the units of thicknesses

Definition at line 53 of file coord_adapt.F90.

53  type(adapt_CS), pointer :: CS !< Unassociated pointer to hold the control structure
54  integer, intent(in) :: nk !< Number of layers in the grid
55  real, dimension(:), intent(in) :: coordinateResolution !< Nominal near-surface resolution [m] or
56  !! other units specified with m_to_H
57  real, optional, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses
58 
59  real :: m_to_H_rescale ! A unit conversion factor.
60 
61  if (associated(cs)) call mom_error(fatal, "init_coord_adapt: CS already associated")
62  allocate(cs)
63  allocate(cs%coordinateResolution(nk))
64 
65  m_to_h_rescale = 1.0 ; if (present(m_to_h)) m_to_h_rescale = m_to_h
66 
67  cs%nk = nk
68  cs%coordinateResolution(:) = coordinateresolution(:)
69 
70  ! Set real parameter default values
71  cs%adaptTimeRatio = 1e-1 ! Nondim.
72  cs%adaptAlpha = 1.0 ! Nondim.
73  cs%adaptZoom = 200.0 * m_to_h_rescale
74  cs%adaptZoomCoeff = 0.0 ! Nondim.
75  cs%adaptBuoyCoeff = 0.0 ! Nondim.
76  cs%adaptDrho0 = 0.5 ! [kg m-3]
77 

References mom_error_handler::mom_error().

Referenced by mom_regridding::initcoord().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_adapt_params()

subroutine, public coord_adapt::set_adapt_params ( type(adapt_cs), pointer  CS,
real, intent(in), optional  adaptTimeRatio,
real, intent(in), optional  adaptAlpha,
real, intent(in), optional  adaptZoom,
real, intent(in), optional  adaptZoomCoeff,
real, intent(in), optional  adaptBuoyCoeff,
real, intent(in), optional  adaptDrho0,
logical, intent(in), optional  adaptDoMin 
)

This subtroutine can be used to set the parameters for coord_adapt module.

Parameters
csThe control structure for this module
[in]adapttimeratioRatio of optimisation and diffusion timescales
[in]adaptalphaNondimensional coefficient determining how much optimisation to apply
[in]adaptzoomNear-surface zooming depth [H ~> m or kg m-2]
[in]adaptzoomcoeffNear-surface zooming coefficient
[in]adaptbuoycoeffStratification-dependent diffusion coefficient
[in]adaptdrho0Reference density difference for stratification-dependent diffusion
[in]adaptdominIf true, form a HYCOM1-like mixed layer by preventing interfaces from becoming shallower than the depths set by coordinateResolution

Definition at line 93 of file coord_adapt.F90.

93  type(adapt_CS), pointer :: CS !< The control structure for this module
94  real, optional, intent(in) :: adaptTimeRatio !< Ratio of optimisation and diffusion timescales
95  real, optional, intent(in) :: adaptAlpha !< Nondimensional coefficient determining
96  !! how much optimisation to apply
97  real, optional, intent(in) :: adaptZoom !< Near-surface zooming depth [H ~> m or kg m-2]
98  real, optional, intent(in) :: adaptZoomCoeff !< Near-surface zooming coefficient
99  real, optional, intent(in) :: adaptBuoyCoeff !< Stratification-dependent diffusion coefficient
100  real, optional, intent(in) :: adaptDrho0 !< Reference density difference for
101  !! stratification-dependent diffusion
102  logical, optional, intent(in) :: adaptDoMin !< If true, form a HYCOM1-like mixed layer by
103  !! preventing interfaces from becoming shallower than
104  !! the depths set by coordinateResolution
105 
106  if (.not. associated(cs)) call mom_error(fatal, "set_adapt_params: CS not associated")
107 
108  if (present(adapttimeratio)) cs%adaptTimeRatio = adapttimeratio
109  if (present(adaptalpha)) cs%adaptAlpha = adaptalpha
110  if (present(adaptzoom)) cs%adaptZoom = adaptzoom
111  if (present(adaptzoomcoeff)) cs%adaptZoomCoeff = adaptzoomcoeff
112  if (present(adaptbuoycoeff)) cs%adaptBuoyCoeff = adaptbuoycoeff
113  if (present(adaptdrho0)) cs%adaptDrho0 = adaptdrho0
114  if (present(adaptdomin)) cs%adaptDoMin = adaptdomin

References mom_error_handler::mom_error().

Referenced by mom_regridding::set_regrid_params().

Here is the call graph for this function:
Here is the caller graph for this function: