11 implicit none ;
private
13 #include <MOM_memory.h>
22 real,
allocatable,
dimension(:) :: coordinateresolution
25 real :: adapttimeratio
34 real :: adaptzoomcoeff
37 real :: adaptbuoycoeff
44 logical :: adaptdomin = .false.
54 integer,
intent(in) :: nk
55 real,
dimension(:),
intent(in) :: coordinateresolution
57 real,
optional,
intent(in) :: m_to_h
59 real :: m_to_h_rescale
61 if (
associated(cs))
call mom_error(fatal,
"init_coord_adapt: CS already associated")
63 allocate(cs%coordinateResolution(nk))
65 m_to_h_rescale = 1.0 ;
if (
present(m_to_h)) m_to_h_rescale = m_to_h
68 cs%coordinateResolution(:) = coordinateresolution(:)
71 cs%adaptTimeRatio = 1e-1
73 cs%adaptZoom = 200.0 * m_to_h_rescale
74 cs%adaptZoomCoeff = 0.0
75 cs%adaptBuoyCoeff = 0.0
85 if (.not.
associated(cs))
return
86 deallocate(cs%coordinateResolution)
91 subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, &
92 adaptBuoyCoeff, adaptDrho0, adaptDoMin)
94 real,
optional,
intent(in) :: adapttimeratio
95 real,
optional,
intent(in) :: adaptalpha
97 real,
optional,
intent(in) :: adaptzoom
98 real,
optional,
intent(in) :: adaptzoomcoeff
99 real,
optional,
intent(in) :: adaptbuoycoeff
100 real,
optional,
intent(in) :: adaptdrho0
102 logical,
optional,
intent(in) :: adaptdomin
106 if (.not.
associated(cs))
call mom_error(fatal,
"set_adapt_params: CS not associated")
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
117 subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
119 type(ocean_grid_type),
intent(in) :: g
123 integer,
intent(in) :: i
124 integer,
intent(in) :: j
125 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(in) :: zint
126 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(in) :: tint
127 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(in) :: sint
128 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
129 real,
dimension(SZK_(GV)+1),
intent(inout) :: znext
133 real :: h_up, b1, b_denom_1, d1, depth, drdz, nominal_z, stretching
134 real,
dimension(SZK_(GV)+1) :: alpha, beta, del2sigma
135 real,
dimension(SZK_(GV)) :: kgrid, c1
141 znext(nz+1) = zint(i,j,nz+1)
144 depth = g%bathyT(i,j) * gv%Z_to_H
154 if (g%mask2dT(i,j-1) > 0.)
then
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)
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)))
166 if (g%mask2dT(i,j+1) > 0.)
then
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)
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)))
178 if (g%mask2dT(i-1,j) > 0.)
then
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)
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)))
190 if (g%mask2dT(i+1,j) > 0.)
then
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)
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)))
209 alpha, beta, 1, nz + 1, tv%eqn_of_state)
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)
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))
224 znext(k) = zint(i,j,k) + del2sigma(k)
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))
238 drdz = drdz / (znext(k) - znext(k+1) + gv%H_subroundoff)
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)
256 b_denom_1 = 1. + d1 * kgrid(k-1)
258 b1 = 1.0 / (b_denom_1 + kgrid(k))
261 c1(k) = kgrid(k) * b1
266 znext(k) = b1 * (znext(k) + kgrid(k-1)*znext(k-1))
270 znext(k) = znext(k) + c1(k)*znext(k+1)
273 if (cs%adaptDoMin)
then
275 stretching = zint(i,j,nz+1) / depth
278 nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
280 znext(k) = max(znext(k), nominal_z)
282 znext(k) = min(znext(k), zint(i,j,nz+1))