coord_adapt.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Regrid columns for the adaptive coordinate
6module coord_adapt
7
8use mom_eos, only : calculate_density_derivs
9use mom_error_handler, only : mom_error, fatal
11use mom_variables, only : ocean_grid_type, thermo_var_ptrs
13
14implicit none ; private
15
16#include <MOM_memory.h>
17
18!> Control structure for adaptive coordinates (coord_adapt).
19type, public :: adapt_cs ; private
20
21 !> Number of layers/levels
22 integer :: nk
23
24 !> Nominal near-surface resolution [H ~> m or kg m-2]
25 real, allocatable, dimension(:) :: coordinateresolution
26
27 !> Ratio of optimisation and diffusion timescales [nondim]
28 real :: adapttimeratio
29
30 !> Nondimensional coefficient determining how much optimisation to apply [nondim]
31 real :: adaptalpha
32
33 !> Near-surface zooming depth [H ~> m or kg m-2]
34 real :: adaptzoom
35
36 !> Near-surface zooming coefficient [nondim]
37 real :: adaptzoomcoeff
38
39 !> Stratification-dependent diffusion coefficient [nondim]
40 real :: adaptbuoycoeff
41
42 !> Reference density difference for stratification-dependent diffusion [R ~> kg m-3]
43 real :: adaptdrho0
44
45 !> If true, form a HYCOM1-like mixed layet by preventing interfaces
46 !! from becoming shallower than the depths set by coordinateResolution
47 logical :: adaptdomin = .false.
48end type adapt_cs
49
51
52contains
53
54!> Initialise an adapt_CS with parameters
55subroutine init_coord_adapt(CS, nk, coordinateResolution, m_to_H, kg_m3_to_R)
56 type(adapt_cs), pointer :: cs !< Unassociated pointer to hold the control structure
57 integer, intent(in) :: nk !< Number of layers in the grid
58 real, dimension(:), intent(in) :: coordinateresolution !< Nominal near-surface resolution [m] or
59 !! other units specified with m_to_H
60 real, intent(in) :: m_to_h !< A conversion factor from m to the units of thicknesses,
61 !! perhaps in units of [H m-1 ~> 1 or kg m-3]
62 real, intent(in) :: kg_m3_to_r !< A conversion factor from kg m-3 to the units of density,
63 !! perhaps in units of [R m3 kg-1 ~> 1]
64
65 if (associated(cs)) call mom_error(fatal, "init_coord_adapt: CS already associated")
66 allocate(cs)
67 allocate(cs%coordinateResolution(nk))
68
69 cs%nk = nk
70 cs%coordinateResolution(:) = coordinateresolution(:)
71
72 ! Set real parameter default values
73 cs%adaptTimeRatio = 1e-1 ! Nondim.
74 cs%adaptAlpha = 1.0 ! Nondim.
75 cs%adaptZoom = 200.0 * m_to_h ! [H ~> m or kg m-2]
76 cs%adaptZoomCoeff = 0.0 ! Nondim.
77 cs%adaptBuoyCoeff = 0.0 ! Nondim.
78 cs%adaptDrho0 = 0.5 * kg_m3_to_r ! [R ~> kg m-3]
79
80end subroutine init_coord_adapt
81
82!> Clean up the coordinate control structure
83subroutine end_coord_adapt(CS)
84 type(adapt_cs), pointer :: cs !< The control structure for this module
85
86 ! nothing to do
87 if (.not. associated(cs)) return
88 deallocate(cs%coordinateResolution)
89 deallocate(cs)
90end subroutine end_coord_adapt
91
92!> This subtroutine can be used to set the parameters for coord_adapt module
93subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, &
94 adaptBuoyCoeff, adaptDrho0, adaptDoMin)
95 type(adapt_cs), pointer :: cs !< The control structure for this module
96 real, optional, intent(in) :: adapttimeratio !< Ratio of optimisation and diffusion timescales [nondim]
97 real, optional, intent(in) :: adaptalpha !< Nondimensional coefficient determining
98 !! how much optimisation to apply [nondim]
99 real, optional, intent(in) :: adaptzoom !< Near-surface zooming depth [H ~> m or kg m-2]
100 real, optional, intent(in) :: adaptzoomcoeff !< Near-surface zooming coefficient [nondim]
101 real, optional, intent(in) :: adaptbuoycoeff !< Stratification-dependent diffusion coefficient [nondim]
102 real, optional, intent(in) :: adaptdrho0 !< Reference density difference for
103 !! stratification-dependent diffusion [R ~> kg m-3]
104 logical, optional, intent(in) :: adaptdomin !< If true, form a HYCOM1-like mixed layer by
105 !! preventing interfaces from becoming shallower than
106 !! the depths set by coordinateResolution
107
108 if (.not. associated(cs)) call mom_error(fatal, "set_adapt_params: CS not associated")
109
110 if (present(adapttimeratio)) cs%adaptTimeRatio = adapttimeratio
111 if (present(adaptalpha)) cs%adaptAlpha = adaptalpha
112 if (present(adaptzoom)) cs%adaptZoom = adaptzoom
113 if (present(adaptzoomcoeff)) cs%adaptZoomCoeff = adaptzoomcoeff
114 if (present(adaptbuoycoeff)) cs%adaptBuoyCoeff = adaptbuoycoeff
115 if (present(adaptdrho0)) cs%adaptDrho0 = adaptdrho0
116 if (present(adaptdomin)) cs%adaptDoMin = adaptdomin
117end subroutine set_adapt_params
118
119subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, nom_depth_H, zNext)
120 type(adapt_cs), intent(in) :: cs !< The control structure for this module
121 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
122 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
123 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
124 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
125 !! thermodynamic variables
126 integer, intent(in) :: i !< The i-index of the column to work on
127 integer, intent(in) :: j !< The j-index of the column to work on
128 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: zint !< Interface heights [H ~> m or kg m-2].
129 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: tint !< Interface temperatures [C ~> degC]
130 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: sint !< Interface salinities [S ~> ppt]
131 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
132 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_h !< The bathymetric depth of this column
133 !! relative to mean sea level or another locally
134 !! valid reference height, converted to thickness
135 !! units [H ~> m or kg m-2]
136 real, dimension(SZK_(GV)+1), intent(inout) :: znext !< updated interface positions [H ~> m or kg m-2]
137
138 ! Local variables
139 integer :: k, nz
140 real :: h_up ! The upwind source grid thickness based on the direction of the
141 ! adjustive fluxes [H ~> m or kg m-2]
142 real :: b1 ! The inverse of the tridiagonal denominator [nondim]
143 real :: b_denom_1 ! The leading term in the tridiagonal denominator [nondim]
144 real :: d1 ! A term in the tridiagonal expressions [nondim]
145 real :: depth ! Depth in thickness units [H ~> m or kg m-2]
146 real :: nominal_z ! A nominal interface position in thickness units [H ~> m or kg m-2]
147 real :: stretching ! A stretching factor for the water column [nondim]
148 real :: drdz ! The vertical density gradient [R H-1 ~> kg m-4 or m-1]
149 real, dimension(SZK_(GV)+1) :: alpha ! drho/dT [R C-1 ~> kg m-3 degC-1]
150 real, dimension(SZK_(GV)+1) :: beta ! drho/dS [R S-1 ~> kg m-3 ppt-1]
151 real, dimension(SZK_(GV)+1) :: del2sigma ! Laplacian of in situ density times grid spacing [R ~> kg m-3]
152 real, dimension(SZK_(GV)+1) :: dh_d2s ! Thickness change in response to del2sigma [H ~> m or kg m-2]
153 real, dimension(SZK_(GV)) :: kgrid ! grid diffusivity on layers [nondim]
154 real, dimension(SZK_(GV)) :: c1 ! A tridiagonal work array [nondim]
155
156 nz = cs%nk
157
158 ! set bottom and surface of zNext
159 znext(1) = 0.
160 znext(nz+1) = zint(i,j,nz+1)
161
162 ! local depth for scaling diffusivity
163 depth = nom_depth_h(i,j)
164
165 ! initialize del2sigma and the thickness change response to it zero
166 del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0
167
168 ! calculate del-squared of neutral density by a
169 ! stencilled finite difference
170 ! TODO: this needs to be adjusted to account for vanished layers near topography
171
172 ! up (j-1)
173 if (g%mask2dT(i,j-1) > 0.0) then
174 call calculate_density_derivs( &
175 0.5 * (tint(i,j,2:nz) + tint(i,j-1,2:nz)), &
176 0.5 * (sint(i,j,2:nz) + sint(i,j-1,2:nz)), &
177 0.5 * (zint(i,j,2:nz) + zint(i,j-1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
178 alpha, beta, tv%eqn_of_state, (/2,nz/) )
179
180 del2sigma(2:nz) = del2sigma(2:nz) + &
181 (alpha(2:nz) * (tint(i,j-1,2:nz) - tint(i,j,2:nz)) + &
182 beta(2:nz) * (sint(i,j-1,2:nz) - sint(i,j,2:nz)))
183 endif
184 ! down (j+1)
185 if (g%mask2dT(i,j+1) > 0.0) then
186 call calculate_density_derivs( &
187 0.5 * (tint(i,j,2:nz) + tint(i,j+1,2:nz)), &
188 0.5 * (sint(i,j,2:nz) + sint(i,j+1,2:nz)), &
189 0.5 * (zint(i,j,2:nz) + zint(i,j+1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
190 alpha, beta, tv%eqn_of_state, (/2,nz/) )
191
192 del2sigma(2:nz) = del2sigma(2:nz) + &
193 (alpha(2:nz) * (tint(i,j+1,2:nz) - tint(i,j,2:nz)) + &
194 beta(2:nz) * (sint(i,j+1,2:nz) - sint(i,j,2:nz)))
195 endif
196 ! left (i-1)
197 if (g%mask2dT(i-1,j) > 0.0) then
198 call calculate_density_derivs( &
199 0.5 * (tint(i,j,2:nz) + tint(i-1,j,2:nz)), &
200 0.5 * (sint(i,j,2:nz) + sint(i-1,j,2:nz)), &
201 0.5 * (zint(i,j,2:nz) + zint(i-1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
202 alpha, beta, tv%eqn_of_state, (/2,nz/) )
203
204 del2sigma(2:nz) = del2sigma(2:nz) + &
205 (alpha(2:nz) * (tint(i-1,j,2:nz) - tint(i,j,2:nz)) + &
206 beta(2:nz) * (sint(i-1,j,2:nz) - sint(i,j,2:nz)))
207 endif
208 ! right (i+1)
209 if (g%mask2dT(i+1,j) > 0.0) then
210 call calculate_density_derivs( &
211 0.5 * (tint(i,j,2:nz) + tint(i+1,j,2:nz)), &
212 0.5 * (sint(i,j,2:nz) + sint(i+1,j,2:nz)), &
213 0.5 * (zint(i,j,2:nz) + zint(i+1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
214 alpha, beta, tv%eqn_of_state, (/2,nz/) )
215
216 del2sigma(2:nz) = del2sigma(2:nz) + &
217 (alpha(2:nz) * (tint(i+1,j,2:nz) - tint(i,j,2:nz)) + &
218 beta(2:nz) * (sint(i+1,j,2:nz) - sint(i,j,2:nz)))
219 endif
220
221 ! at this point, del2sigma contains the local neutral density curvature at
222 ! h-points, on interfaces
223 ! we need to divide by drho/dz to give an interfacial displacement
224 !
225 ! a positive curvature means we're too light relative to adjacent columns,
226 ! so del2sigma needs to be positive too (push the interface deeper)
227 call calculate_density_derivs(tint(i,j,:), sint(i,j,:), zint(i,j,:) * (gv%H_to_RZ * gv%g_Earth), &
228 alpha, beta, tv%eqn_of_state, (/1,nz+1/) )
229 do k = 2, nz
230 ! TODO make lower bound here configurable
231 dh_d2s(k) = del2sigma(k) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &
232 max(alpha(k) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + &
233 beta(k) * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20*us%kg_m3_to_R)
234
235 ! don't move the interface so far that it would tangle with another
236 ! interface in the direction we're moving (or exceed a Nyquist limit
237 ! that could cause oscillations of the interface)
238 h_up = merge(h(i,j,k), h(i,j,k-1), dh_d2s(k) > 0.)
239 dh_d2s(k) = 0.5 * cs%adaptAlpha * &
240 sign(min(abs(del2sigma(k)), 0.5 * h_up), dh_d2s(k))
241
242 ! update interface positions so we can diffuse them
243 znext(k) = zint(i,j,k) + dh_d2s(k)
244 enddo
245
246 ! solve diffusivity equation to smooth grid
247 ! upper diagonal coefficients: -kGrid(2:nz)
248 ! lower diagonal coefficients: -kGrid(1:nz-1)
249 ! diagonal coefficients: 1 + (kGrid(1:nz-1) + kGrid(2:nz))
250 !
251 ! first, calculate the diffusivities within layers
252 do k = 1, nz
253 ! calculate the dr bit of drdz
254 drdz = 0.5 * (alpha(k) + alpha(k+1)) * (tint(i,j,k+1) - tint(i,j,k)) + &
255 0.5 * (beta(k) + beta(k+1)) * (sint(i,j,k+1) - sint(i,j,k))
256 ! divide by dz from the new interface positions
257 drdz = drdz / (znext(k) - znext(k+1) + gv%H_subroundoff)
258 ! don't do weird stuff in unstably-stratified regions
259 drdz = max(drdz, 0.)
260
261 ! set vertical grid diffusivity
262 kgrid(k) = (cs%adaptTimeRatio * nz**2 * depth) * &
263 ( cs%adaptZoomCoeff / (cs%adaptZoom + 0.5*(znext(k) + znext(k+1))) + &
264 (cs%adaptBuoyCoeff * drdz / cs%adaptDrho0) + &
265 max(1.0 - cs%adaptZoomCoeff - cs%adaptBuoyCoeff, 0.0) / depth)
266 enddo
267
268 ! initial denominator (first diagonal element)
269 b1 = 1.0
270 ! initial Q_1 = 1 - q_1 = 1 - 0/1
271 d1 = 1.0
272 ! work on all interior interfaces
273 do k = 2, nz
274 ! calculate numerator of Q_k
275 b_denom_1 = 1. + d1 * kgrid(k-1)
276 ! update denominator for k
277 b1 = 1.0 / (b_denom_1 + kgrid(k))
278
279 ! calculate q_k
280 c1(k) = kgrid(k) * b1
281 ! update Q_k = 1 - q_k
282 d1 = b_denom_1 * b1
283
284 ! update RHS
285 znext(k) = b1 * (znext(k) + kgrid(k-1)*znext(k-1))
286 enddo
287 ! final substitution
288 do k = nz, 2, -1
289 znext(k) = znext(k) + c1(k)*znext(k+1)
290 enddo
291
292 if (cs%adaptDoMin) then
293 nominal_z = 0.
294 stretching = zint(i,j,nz+1) / depth
295
296 do k = 2, nz+1
297 nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
298 ! take the deeper of the calculated and nominal positions
299 znext(k) = max(znext(k), nominal_z)
300 ! interface can't go below topography
301 znext(k) = min(znext(k), zint(i,j,nz+1))
302 enddo
303 endif
304end subroutine build_adapt_column
305
306end module coord_adapt