MOM_tracer_diabatic.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!> This module contains routines that implement physical fluxes of tracers (e.g. due
6!! to surface fluxes or mixing). These are intended to be called from call_tracer_column_fns
7!! in the MOM_tracer_flow_control module.
9
10use mom_grid, only : ocean_grid_type
12use mom_forcing_type, only : forcing
13use mom_error_handler, only : mom_error, fatal, warning
14
15implicit none ; private
16
17#include <MOM_memory.h>
20
21contains
22
23!> This subroutine solves a tridiagonal equation for the final tracer concentrations after the
24!! dual-entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking
25!! is implemented with an fully implicit upwind advection scheme. Alternate time units can be
26!! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent.
27subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
28 sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
29 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
30 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
31 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment
32 !! [H ~> m or kg m-2]
33 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< amount of fluid entrained from the layer
34 !! above [H ~> m or kg m-2]
35 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< amount of fluid entrained from the layer
36 !! below [H ~> m or kg m-2]
37 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU]
38 real, intent(in) :: dt !< amount of time covered by this call [T ~> s]
39 real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of
40 !! [CU R Z T-1 ~> CU kg m-2 s-1] or
41 !! [CU H ~> CU m or CU kg m-2] if
42 !! convert_flux_in is .false.
43 real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the
44 !! tracer in [CU R Z T-1 ~> CU kg m-2 s-1] or
45 !! [CU H ~> CU m or CU kg m-2] if
46 !! convert_flux_in is .false.
47 real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir
48 !! [CU R Z ~> CU kg m-2]
49 real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks
50 !! [Z T-1 ~> m s-1]
51 logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs
52 !! to be integrated in time
53
54 ! local variables
55 real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2].
56 real, dimension(SZI_(G),SZJ_(G)) :: &
57 sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2].
58 btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2].
59 real, dimension(SZI_(G)) :: &
60 b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
61 d1 !! d1=1-c1 is used by the tridiagonal solver [nondim].
62 real :: c1(szi_(g),szk_(gv)) !< c1 is used by the tridiagonal solver [nondim].
63 real :: h_minus_dsink(szi_(g),szk_(gv)) !< The layer thickness minus the
64 !! difference in sinking rates across the layer [H ~> m or kg m-2].
65 !! By construction, 0 <= h_minus_dsink < h_work.
66 real :: sink(szi_(g),szk_(gv)+1) !< The tracer's sinking distances at the
67 !! interfaces, limited to prevent characteristics from
68 !! crossing within a single timestep [H ~> m or kg m-2].
69 real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2].
70 real :: h_tr !< h_tr is h at tracer points with a h_neglect added to
71 !! ensure positive definiteness [H ~> m or kg m-2].
72 real :: h_neglect !< A thickness that is so small it is usually lost
73 !! in roundoff and can be neglected [H ~> m or kg m-2].
74 logical :: convert_flux = .true.
75
76 integer :: i, j, k, is, ie, js, je, nz
77 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
78
79 if (nz == 1) then
80 call mom_error(warning, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//&
81 "with only one vertical level")
82 return
83 endif
84
85 if (present(convert_flux_in)) convert_flux = convert_flux_in
86 h_neglect = gv%H_subroundoff
87 sink_dist = 0.0
88 if (present(sink_rate)) sink_dist = (dt*sink_rate) * gv%Z_to_H
89 !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
90 !$OMP do
91 do j=js,je ; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo
92 if (present(sfc_flux)) then
93 if (convert_flux) then
94 !$OMP do
95 do j=js,je ; do i=is,ie
96 sfc_src(i,j) = (sfc_flux(i,j)*dt) * gv%RZ_to_H
97 enddo ; enddo
98 else
99 !$OMP do
100 do j=js,je ; do i=is,ie
101 sfc_src(i,j) = sfc_flux(i,j)
102 enddo ; enddo
103 endif
104 endif
105 if (present(btm_flux)) then
106 if (convert_flux) then
107 !$OMP do
108 do j=js,je ; do i=is,ie
109 btm_src(i,j) = (btm_flux(i,j)*dt) * gv%RZ_to_H
110 enddo ; enddo
111 else
112 !$OMP do
113 do j=js,je ; do i=is,ie
114 btm_src(i,j) = btm_flux(i,j)
115 enddo ; enddo
116 endif
117 endif
118
119 if (present(sink_rate)) then
120 !$OMP do
121 do j=js,je
122 ! Find the sinking rates at all interfaces, limiting them if necesary
123 ! so that the characteristics do not cross within a timestep.
124 ! If a non-constant sinking rate were used, that would be incorprated
125 ! here.
126 if (present(btm_reservoir)) then
127 do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
128 do k=2,nz ; do i=is,ie
129 sink(i,k) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
130 enddo ; enddo
131 else
132 do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
133 ! Find the limited sinking distance at the interfaces.
134 do k=nz,2,-1 ; do i=is,ie
135 if (sink(i,k+1) >= sink_dist) then
136 sink(i,k) = sink_dist
137 h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,k+1) - sink(i,k))
138 elseif (sink(i,k+1) + h_old(i,j,k) < sink_dist) then
139 sink(i,k) = sink(i,k+1) + h_old(i,j,k)
140 h_minus_dsink(i,k) = 0.0
141 else
142 sink(i,k) = sink_dist
143 h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,k+1)) - sink(i,k)
144 endif
145 enddo ; enddo
146 endif
147 do i=is,ie
148 sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
149 enddo
150
151 ! Now solve the tridiagonal equation for the tracer concentrations.
152 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
153 b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
154 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
155 d1(i) = b_denom_1 * b1(i)
156 h_tr = h_old(i,j,1) + h_neglect
157 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
158 endif ; enddo
159 do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
160 c1(i,k) = eb(i,j,k-1) * b1(i)
161 b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,k)) + &
162 h_neglect
163 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
164 d1(i) = b_denom_1 * b1(i)
165 h_tr = h_old(i,j,k) + h_neglect
166 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
167 (ea(i,j,k) + sink(i,k)) * tr(i,j,k-1))
168 endif ; enddo ; enddo
169 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
170 c1(i,nz) = eb(i,j,nz-1) * b1(i)
171 b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
172 h_neglect
173 b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
174 h_tr = h_old(i,j,nz) + h_neglect
175 tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
176 (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
177 endif ; enddo
178 if (present(btm_reservoir)) then ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
179 btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * gv%H_to_RZ
180 endif ; enddo ; endif
181
182 do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
183 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
184 endif ; enddo ; enddo
185 enddo
186 else
187 !$OMP do
188 do j=js,je
189 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
190 h_tr = h_old(i,j,1) + h_neglect
191 b_denom_1 = h_tr + ea(i,j,1)
192 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
193 d1(i) = h_tr * b1(i)
194 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
195 endif ; enddo
196 do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
197 c1(i,k) = eb(i,j,k-1) * b1(i)
198 h_tr = h_old(i,j,k) + h_neglect
199 b_denom_1 = h_tr + d1(i) * ea(i,j,k)
200 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
201 d1(i) = b_denom_1 * b1(i)
202 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1))
203 endif ; enddo ; enddo
204 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
205 c1(i,nz) = eb(i,j,nz-1) * b1(i)
206 h_tr = h_old(i,j,nz) + h_neglect
207 b_denom_1 = h_tr + d1(i)*ea(i,j,nz)
208 b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz))
209 tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
210 ea(i,j,nz) * tr(i,j,nz-1))
211 endif ; enddo
212 do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
213 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
214 endif ; enddo ; enddo
215 enddo
216 endif
217 !$OMP end parallel
218
219end subroutine tracer_vertdiff
220
221
222!> This subroutine solves a tridiagonal equation for the final tracer concentrations after
223!! Eulerian mixing, and possibly sinking or surface and bottom sources, are applied. The sinking
224!! is implemented with an fully implicit upwind advection scheme. Alternate time units can be
225!! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent.
226subroutine tracer_vertdiff_eulerian(h_old, ent, dt, tr, G, GV, &
227 sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
228 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
229 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
230 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment
231 !! [H ~> m or kg m-2]
232 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< Amount of fluid mixed across interfaces
233 !! [H ~> m or kg m-2]
234 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU]
235 real, intent(in) :: dt !< amount of time covered by this call [T ~> s]
236 real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of
237 !! [CU R Z T-1 ~> CU kg m-2 s-1] or
238 !! [CU H ~> CU m or CU kg m-2] if
239 !! convert_flux_in is .false.
240 real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the
241 !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or
242 !! [CU H ~> CU m or CU kg m-2] if
243 !! convert_flux_in is .false.
244 real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir
245 !! [CU R Z ~> CU kg m-2]
246 real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks
247 !! [Z T-1 ~> m s-1]
248 logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs
249 !! to be integrated in time
250
251 ! local variables
252 real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2].
253 real, dimension(SZI_(G),SZJ_(G)) :: &
254 sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2].
255 btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2].
256 real, dimension(SZI_(G)) :: &
257 b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
258 d1 !! d1=1-c1 is used by the tridiagonal solver [nondim].
259 real :: c1(szi_(g),szk_(gv)) !< c1 is used by the tridiagonal solver [nondim].
260 real :: h_minus_dsink(szi_(g),szk_(gv)) !< The layer thickness minus the
261 !! difference in sinking rates across the layer [H ~> m or kg m-2].
262 !! By construction, 0 <= h_minus_dsink < h_work.
263 real :: sink(szi_(g),szk_(gv)+1) !< The tracer's sinking distances at the
264 !! interfaces, limited to prevent characteristics from
265 !! crossing within a single timestep [H ~> m or kg m-2].
266 real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2].
267 real :: h_tr !< h_tr is h at tracer points with a h_neglect added to
268 !! ensure positive definiteness [H ~> m or kg m-2].
269 real :: h_neglect !< A thickness that is so small it is usually lost
270 !! in roundoff and can be neglected [H ~> m or kg m-2].
271 logical :: convert_flux
272
273 integer :: i, j, k, is, ie, js, je, nz
274 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
275
276 if (nz == 1) then
277 call mom_error(warning, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//&
278 "with only one vertical level")
279 return
280 endif
281
282 convert_flux = .true.
283 if (present(convert_flux_in)) convert_flux = convert_flux_in
284 h_neglect = gv%H_subroundoff
285 sink_dist = 0.0
286 if (present(sink_rate)) sink_dist = (dt*sink_rate) * gv%Z_to_H
287 !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
288 !$OMP do
289 do j=js,je ; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo
290 if (present(sfc_flux)) then
291 if (convert_flux) then
292 !$OMP do
293 do j=js,je ; do i=is,ie
294 sfc_src(i,j) = (sfc_flux(i,j)*dt) * gv%RZ_to_H
295 enddo ; enddo
296 else
297 !$OMP do
298 do j=js,je ; do i=is,ie
299 sfc_src(i,j) = sfc_flux(i,j)
300 enddo ; enddo
301 endif
302 endif
303 if (present(btm_flux)) then
304 if (convert_flux) then
305 !$OMP do
306 do j=js,je ; do i=is,ie
307 btm_src(i,j) = (btm_flux(i,j)*dt) * gv%kg_m2_to_H
308 enddo ; enddo
309 else
310 !$OMP do
311 do j=js,je ; do i=is,ie
312 btm_src(i,j) = btm_flux(i,j)
313 enddo ; enddo
314 endif
315 endif
316
317 if (present(sink_rate)) then
318 !$OMP do
319 do j=js,je
320 ! Find the sinking rates at all interfaces, limiting them if necesary
321 ! so that the characteristics do not cross within a timestep.
322 ! If a non-constant sinking rate were used, that would be incorprated
323 ! here.
324 if (present(btm_reservoir)) then
325 do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
326 do k=2,nz ; do i=is,ie
327 sink(i,k) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
328 enddo ; enddo
329 else
330 do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
331 ! Find the limited sinking distance at the interfaces.
332 do k=nz,2,-1 ; do i=is,ie
333 if (sink(i,k+1) >= sink_dist) then
334 sink(i,k) = sink_dist
335 h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,k+1) - sink(i,k))
336 elseif (sink(i,k+1) + h_old(i,j,k) < sink_dist) then
337 sink(i,k) = sink(i,k+1) + h_old(i,j,k)
338 h_minus_dsink(i,k) = 0.0
339 else
340 sink(i,k) = sink_dist
341 h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,k+1)) - sink(i,k)
342 endif
343 enddo ; enddo
344 endif
345 do i=is,ie
346 sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
347 enddo
348
349 ! Now solve the tridiagonal equation for the tracer concentrations.
350 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
351 b_denom_1 = h_minus_dsink(i,1) + ent(i,j,1) + h_neglect
352 b1(i) = 1.0 / (b_denom_1 + ent(i,j,2))
353 d1(i) = b_denom_1 * b1(i)
354 h_tr = h_old(i,j,1) + h_neglect
355 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
356 endif ; enddo
357 do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
358 c1(i,k) = ent(i,j,k) * b1(i)
359 b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ent(i,j,k) + sink(i,k)) + &
360 h_neglect
361 b1(i) = 1.0 / (b_denom_1 + ent(i,j,k+1))
362 d1(i) = b_denom_1 * b1(i)
363 h_tr = h_old(i,j,k) + h_neglect
364 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
365 (ent(i,j,k) + sink(i,k)) * tr(i,j,k-1))
366 endif ; enddo ; enddo
367 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
368 c1(i,nz) = ent(i,j,nz) * b1(i)
369 b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ent(i,j,nz) + sink(i,nz)) + &
370 h_neglect
371 b1(i) = 1.0 / (b_denom_1 + ent(i,j,nz+1))
372 h_tr = h_old(i,j,nz) + h_neglect
373 tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
374 (ent(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
375 endif ; enddo
376 if (present(btm_reservoir)) then ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
377 btm_reservoir(i,j) = btm_reservoir(i,j) + (sink(i,nz+1)*tr(i,j,nz)) * gv%H_to_RZ
378 endif ; enddo ; endif
379
380 do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
381 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
382 endif ; enddo ; enddo
383 enddo
384 else
385 !$OMP do
386 do j=js,je
387 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
388 h_tr = h_old(i,j,1) + h_neglect
389 b_denom_1 = h_tr + ent(i,j,1)
390 b1(i) = 1.0 / (b_denom_1 + ent(i,j,2))
391 d1(i) = h_tr * b1(i)
392 tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j)
393 endif ; enddo
394 do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
395 c1(i,k) = ent(i,j,k) * b1(i)
396 h_tr = h_old(i,j,k) + h_neglect
397 b_denom_1 = h_tr + d1(i) * ent(i,j,k)
398 b1(i) = 1.0 / (b_denom_1 + ent(i,j,k+1))
399 d1(i) = b_denom_1 * b1(i)
400 tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,k) * tr(i,j,k-1))
401 endif ; enddo ; enddo
402 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
403 c1(i,nz) = ent(i,j,nz) * b1(i)
404 h_tr = h_old(i,j,nz) + h_neglect
405 b_denom_1 = h_tr + d1(i)*ent(i,j,nz)
406 b1(i) = 1.0 / ( b_denom_1 + ent(i,j,nz+1))
407 tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
408 ent(i,j,nz) * tr(i,j,nz-1))
409 endif ; enddo
410 do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
411 tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
412 endif ; enddo ; enddo
413 enddo
414 endif
415 !$OMP end parallel
416
417end subroutine tracer_vertdiff_eulerian
418
419
420!> This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90
421!! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface
422!! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif
423subroutine applytracerboundaryfluxesinout(G, GV, Tr, dt, fluxes, h, evap_CFL_limit, minimum_forcing_depth, &
424 in_flux_optional, out_flux_optional, update_h_opt)
425
426 type(ocean_grid_type), intent(in ) :: g !< Grid structure
427 type(verticalgrid_type), intent(in ) :: gv !< ocean vertical grid structure
428 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< Tracer concentration on T-cell [conc]
429 real, intent(in ) :: dt !< Time-step over which forcing is applied [T ~> s]
430 type(forcing), intent(in ) :: fluxes !< Surface fluxes container
431 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
432 real, intent(in ) :: evap_cfl_limit !< Limit on the fraction of the
433 !! water that can be fluxed out of the top
434 !! layer in a timestep [nondim]
435 real, intent(in ) :: minimum_forcing_depth !< The smallest depth over
436 !! which fluxes can be applied [H ~> m or kg m-2]
437 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: in_flux_optional !< The total time-integrated
438 !! amount of tracer that enters with freshwater
439 !! [conc H ~> conc m or conc kg m-2]
440 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: out_flux_optional !< The total time-integrated
441 !! amount of tracer that leaves with freshwater
442 !! [conc H ~> conc m or conc kg m-2]
443 logical, optional, intent(in) :: update_h_opt !< Optional flag to determine whether
444 !! h should be updated
445
446 integer, parameter :: maxgroundings = 5
447 integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
448 real :: iforcingdepthscale ! The inverse of the scale over which to apply forcing [H-1 ~> m-1 or m2 kg-1]
449 real :: dthickness ! The change in a layer's thickness [H ~> m or kg m-2]
450 real :: dtracer ! The change in the integrated tracer content of a layer [conc H ~> conc m or conc kg m-2]
451 real :: fractionofforcing ! The fraction of the forcing to apply to a layer [nondim]
452 real :: hold ! The layer thickness before surface forcing is applied [H ~> m or kg m-2]
453 real :: ithickness ! The inverse of the new layer thickness [H-1 ~> m-1 or m2 kg-1]
454
455 real :: h2d(szi_(g),szk_(gv)) ! A 2-d work copy of layer thicknesses [H ~> m or kg m-2]
456 real :: tr2d(szi_(g),szk_(gv)) ! A 2-d work copy of tracer concentrations [conc]
457 real :: in_flux(szi_(g),szj_(g)) ! The total time-integrated amount of tracer that
458 ! enters with freshwater [conc H ~> conc m or conc kg m-2]
459 real :: out_flux(szi_(g),szj_(g)) ! The total time-integrated amount of tracer that
460 ! leaves with freshwater [conc H ~> conc m or conc kg m-2]
461 real :: netmassin(szi_(g)) ! The remaining mass entering ocean surface [H ~> m or kg m-2]
462 real :: netmassout(szi_(g)) ! The remaining mass leaving ocean surface [H ~> m or kg m-2]
463 real :: in_flux_1d(szi_(g)) ! The remaining amount of tracer that enters with
464 ! the freshwater [conc H ~> conc m or conc kg m-2]
465 real :: out_flux_1d(szi_(g)) ! The remaining amount of tracer that leaves with
466 ! the freshwater [conc H ~> conc m or conc kg m-2]
467 real :: hgrounding(maxgroundings) ! The remaining fresh water flux that was not able to be
468 ! supplied from a column that grounded out [H ~> m or kg m-2]
469 logical :: update_h
470 integer :: i, j, is, ie, js, je, k, nz
471 character(len=45) :: mesg
472
473 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
474
475 ! If no freshwater fluxes, nothing needs to be done in this routine
476 if ( (.not. associated(fluxes%netMassIn)) .or. (.not. associated(fluxes%netMassOut)) ) return
477
478 in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0
479 if (present(in_flux_optional)) then
480 do j=js,je ; do i=is,ie
481 in_flux(i,j) = in_flux_optional(i,j)
482 enddo ; enddo
483 endif
484 if (present(out_flux_optional)) then
485 do j=js,je ; do i=is,ie
486 out_flux(i,j) = out_flux_optional(i,j)
487 enddo ; enddo
488 endif
489
490 if (present(update_h_opt)) then
491 update_h = update_h_opt
492 else
493 update_h = .true.
494 endif
495
496 numberofgroundings = 0
497
498!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt, &
499!$OMP IforcingDepthScale,minimum_forcing_depth, &
500!$OMP numberOfGroundings,iGround,jGround,update_h, &
501!$OMP in_flux,out_flux,hGrounding,evap_CFL_limit) &
502!$OMP private(h2d,Tr2d,netMassIn,netMassOut, &
503!$OMP in_flux_1d,out_flux_1d,fractionOfForcing, &
504!$OMP dThickness,dTracer,hOld,Ithickness)
505
506 ! Work in vertical slices for efficiency
507 do j=js,je
508
509 ! Copy state into 2D-slice arrays
510 do k=1,nz ; do i=is,ie
511 h2d(i,k) = h(i,j,k)
512 tr2d(i,k) = tr(i,j,k)
513 enddo ; enddo
514
515 do i = is,ie
516 in_flux_1d(i) = in_flux(i,j)
517 out_flux_1d(i) = out_flux(i,j)
518 enddo
519 ! The surface forcing is contained in the fluxes type.
520 ! We aggregate the thermodynamic forcing for a time step into the following:
521 ! These should have been set and stored during a call to applyBoundaryFluxesInOut
522 ! netMassIn = net mass entering at ocean surface over a timestep
523 ! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
524 ! netMassOut < 0 means mass leaves ocean.
525
526 ! Note here that the aggregateFW flag has already been taken care of in the call to
527 ! applyBoundaryFluxesInOut
528 do i=is,ie
529 netmassout(i) = fluxes%netMassOut(i,j)
530 netmassin(i) = fluxes%netMassIn(i,j)
531 enddo
532
533 ! Apply the surface boundary fluxes in three steps:
534 ! A/ update concentration from mass entering the ocean
535 ! B/ update concentration from mass leaving ocean.
536 do i=is,ie
537 if (g%mask2dT(i,j)>0.) then
538
539 ! A/ Update tracer due to incoming mass flux.
540 do k=1,1
541
542 ! Change in state due to forcing
543 dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
544 dtracer = 0.
545
546 ! Update the forcing by the part to be consumed within the present k-layer.
547 ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
548 netmassin(i) = netmassin(i) - dthickness
549 dtracer = dtracer + in_flux_1d(i)
550 in_flux_1d(i) = 0.0
551
552 ! Update state
553 hold = h2d(i,k) ! Keep original thickness in hand
554 h2d(i,k) = h2d(i,k) + dthickness ! New thickness
555 if (h2d(i,k) > 0.0) then
556 ithickness = 1.0/h2d(i,k) ! Inverse new thickness
557 ! The "if"s below avoid changing T/S by roundoff unnecessarily
558 if (dthickness /= 0. .or. dtracer /= 0.) tr2d(i,k) = (hold*tr2d(i,k)+ dtracer)*ithickness
559 endif
560
561 enddo ! k=1,1
562
563 ! B/ Update tracer from mass leaving ocean
564 do k=1,nz
565
566 ! Place forcing into this layer if this layer has nontrivial thickness.
567 ! For layers thin relative to 1/IforcingDepthScale, then distribute
568 ! forcing into deeper layers.
569 iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
570 ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
571 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
572
573 ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
574 ! limit the forcing applied to this cell, leaving the remaining forcing to
575 ! be distributed downwards.
576 if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
577 fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
578 endif
579
580 ! Change in state due to forcing
581 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
582 ! Note this is slightly different to how salt is currently treated
583 dtracer = fractionofforcing*out_flux_1d(i)
584
585 ! Update the forcing by the part to be consumed within the present k-layer.
586 ! If fractionOfForcing = 1, then new netMassOut vanishes.
587 netmassout(i) = netmassout(i) - dthickness
588 out_flux_1d(i) = out_flux_1d(i) - dtracer
589
590 ! Update state by the appropriate increment.
591 hold = h2d(i,k) ! Keep original thickness in hand
592 h2d(i,k) = h2d(i,k) + dthickness ! New thickness
593 if (h2d(i,k) > 0.) then
594 ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
595 tr2d(i,k) = (hold*tr2d(i,k) + dtracer)*ithickness
596 endif
597
598 enddo ! k
599
600 endif
601
602 ! If anything remains after the k-loop, then we have grounded out, which is a problem.
603 if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0) then
604!$OMP critical
605 numberofgroundings = numberofgroundings +1
606 if (numberofgroundings<=maxgroundings) then
607 iground(numberofgroundings) = i ! Record i,j location of event for
608 jground(numberofgroundings) = j ! warning message
609 hgrounding(numberofgroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i))
610 endif
611!$OMP end critical
612 endif
613
614 enddo ! i
615
616 ! Step C/ copy updated tracer concentration from the 2d slice now back into model state.
617 do k=1,nz ; do i=is,ie
618 tr(i,j,k) = tr2d(i,k)
619 enddo ; enddo
620
621 if (update_h) then
622 do k=1,nz ; do i=is,ie
623 h(i,j,k) = h2d(i,k)
624 enddo ; enddo
625 endif
626
627 enddo ! j-loop finish
628
629 if (numberofgroundings>0) then
630 do i = 1, min(numberofgroundings, maxgroundings)
631 write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
632 g%geoLatT( iground(i), jground(i)) , hgrounding(i)
633 call mom_error(warning, "MOM_tracer_diabatic.F90, applyTracerBoundaryFluxesInOut(): "//&
634 "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.)
635 enddo
636
637 if (numberofgroundings - maxgroundings > 0) then
638 write(mesg, '(I0)') numberofgroundings - maxgroundings
639 call mom_error(warning, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
640 trim(mesg) // " groundings remaining", all_print=.true.)
641 endif
642 endif
643
645end module mom_tracer_diabatic