MOM_offline_aux.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!> Contains routines related to offline transport of tracers. These routines are likely to be called from
6!> the MOM_offline_main module
8
10use mom_domains, only : pass_var, pass_vector, to_all
11use mom_diag_mediator, only : post_data
12use mom_error_handler, only : calltree_enter, calltree_leave, mom_error, fatal, warning, is_root_pe
13use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_forcing_type, only : forcing
15use mom_grid, only : ocean_grid_type
16use mom_io, only : mom_read_data, mom_read_vector, center
17use mom_opacity, only : optics_type
18use mom_time_manager, only : time_type, operator(-)
22use astronomy_mod, only : orbital_time, diurnal_solar, daily_mean_solar
23
24implicit none ; private
25
37
38#include "MOM_memory.h"
39
40contains
41
42!> This updates thickness based on the convergence of horizontal mass fluxes
43!! NOTE: Only used in non-ALE mode
44subroutine update_h_horizontal_flux(G, GV, uhtr, vhtr, h_pre, h_new)
45 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
46 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
47 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
48 intent(in) :: uhtr !< Accumulated mass flux through zonal face [H L2 ~> m3 or kg]
49 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
50 intent(in) :: vhtr !< Accumulated mass flux through meridional face [H L2 ~> m3 or kg]
51 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
52 intent(in) :: h_pre !< Previous layer thicknesses [H ~> m or kg m-2]
53 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
54 intent(inout) :: h_new !< Updated layer thicknesses [H ~> m or kg m-2]
55
56 ! Local variables
57 integer :: i, j, k, is, ie, js, je, nz
58 ! Set index-related variables for fields on T-grid
59 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
60
61 do k=1,nz
62 do i=is-1,ie+1 ; do j=js-1,je+1
63
64 h_new(i,j,k) = max(0.0, g%areaT(i,j)*h_pre(i,j,k) + &
65 ((uhtr(i-1,j,k) - uhtr(i,j,k)) + (vhtr(i,j-1,k) - vhtr(i,j,k))))
66
67 ! Convert back to thickness
68 h_new(i,j,k) = max(gv%Angstrom_H, h_new(i,j,k) * g%IareaT(i,j))
69
70 enddo ; enddo
71 enddo
72end subroutine update_h_horizontal_flux
73
74!> Updates layer thicknesses due to vertical mass transports
75!! NOTE: Only used in non-ALE configuration
76subroutine update_h_vertical_flux(G, GV, ea, eb, h_pre, h_new)
77 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
78 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
79 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
80 intent(in) :: ea !< Mass of fluid entrained from the layer
81 !! above within this timestep [H ~> m or kg m-2]
82 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
83 intent(in) :: eb !< Mass of fluid entrained from the layer
84 !! below within this timestep [H ~> m or kg m-2]
85 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
86 intent(in) :: h_pre !< Layer thicknesses at the end of the previous
87 !! step [H ~> m or kg m-2]
88 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
89 intent(inout) :: h_new !< Updated layer thicknesses [H ~> m or kg m-2]
90
91 ! Local variables
92 integer :: i, j, k, is, ie, js, je, nz
93 ! Set index-related variables for fields on T-grid
94 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
95
96 ! Update h_new with convergence of vertical mass transports
97 do j=js-1,je+1
98 do i=is-1,ie+1
99 ! Top layer
100 h_new(i,j,1) = max(0.0, h_pre(i,j,1) + ((eb(i,j,1) - ea(i,j,2)) + ea(i,j,1)))
101
102 ! Bottom layer
103 h_new(i,j,nz) = max(0.0, h_pre(i,j,nz) + ((ea(i,j,nz) - eb(i,j,nz-1)) + eb(i,j,nz)))
104 enddo
105
106 ! Interior layers
107 do k=2,nz-1 ; do i=is-1,ie+1
108 h_new(i,j,k) = max(0.0, h_pre(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
109 (eb(i,j,k) - ea(i,j,k+1))))
110 enddo ; enddo
111 enddo
112
113end subroutine update_h_vertical_flux
114
115!> This routine limits the mass fluxes so that the a layer cannot be completely depleted.
116!! NOTE: Only used in non-ALE mode
117subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre)
118 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
119 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
120 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
121 intent(inout) :: uh !< Mass flux through zonal face [H L2 ~> m3 or kg]
122 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
123 intent(inout) :: vh !< Mass flux through meridional face [H L2 ~> m3 or kg]
124 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
125 intent(inout) :: ea !< Mass of fluid entrained from the layer
126 !! above within this timestep [H ~> m or kg m-2]
127 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
128 intent(inout) :: eb !< Mass of fluid entrained from the layer
129 !! below within this timestep [H ~> m or kg m-2]
130 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
131 intent(in) :: h_pre !< Layer thicknesses at the end of the previous
132 !! step [H ~> m or kg m-2]
133
134 ! Local variables
135 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: top_flux ! Net upward fluxes through the layer
136 ! top [H ~> m or kg m-2]
137 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: bottom_flux ! Net downward fluxes through the layer
138 ! bottom [H ~> m or kg m-2]
139 real :: pos_flux ! Net flux out of cell [H L2 ~> m3 or kg]
140 real :: hvol ! Cell volume [H L2 ~> m3 or kg]
141 real :: scale_factor ! A nondimensional rescaling factor between 0 and 1 [nondim]
142 real :: max_off_cfl ! The maximum permitted fraction that can leave in a timestep [nondim]
143 integer :: i, j, k, is, ie, js, je, nz
144
145 max_off_cfl = 0.5
146
147 ! In this subroutine, fluxes out of the box are scaled away if they deplete
148 ! the layer, note that we define the positive direction as flux out of the box.
149 ! Hence, uh(I-1) is multipled by negative one, but uh(I) is not
150
151 ! Set index-related variables for fields on T-grid
152 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
153
154 ! Calculate top and bottom fluxes from ea and eb. Note the explicit negative signs
155 ! to enforce the positive out convention
156 k = 1
157 do j=js-1,je+1 ; do i=is-1,ie+1
158 top_flux(i,j,k) = -ea(i,j,k)
159 bottom_flux(i,j,k) = -(eb(i,j,k)-ea(i,j,k+1))
160 enddo ; enddo
161
162 do k=2,nz-1 ; do j=js-1,je+1 ; do i=is-1,ie+1
163 top_flux(i,j,k) = -(ea(i,j,k)-eb(i,j,k-1))
164 bottom_flux(i,j,k) = -(eb(i,j,k)-ea(i,j,k+1))
165 enddo ; enddo ; enddo
166
167 k=nz
168 do j=js-1,je+1 ; do i=is-1,ie+1
169 top_flux(i,j,k) = -(ea(i,j,k)-eb(i,j,k-1))
170 bottom_flux(i,j,k) = -eb(i,j,k)
171 enddo ; enddo
172
173
174 ! Calculate sum of positive fluxes (negatives applied to enforce convention)
175 ! in a given cell and scale it back if it would deplete a layer
176 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
177
178 hvol = h_pre(i,j,k) * g%areaT(i,j)
179 pos_flux = ((max(0.0, -uh(i-1,j,k)) + max(0.0, uh(i,j,k))) + &
180 (max(0.0, -vh(i,j-1,k)) + max(0.0, vh(i,j,k)))) + &
181 (max(0.0, top_flux(i,j,k)) + max(0.0, bottom_flux(i,j,k))) * g%areaT(i,j)
182
183 if ((pos_flux > hvol) .and. (pos_flux > 0.0)) then
184 scale_factor = (hvol / pos_flux) * max_off_cfl
185 else ! Don't scale
186 scale_factor = 1.0
187 endif
188
189 ! Scale horizontal fluxes
190 if (-uh(i-1,j,k) > 0.0) uh(i-1,j,k) = uh(i-1,j,k) * scale_factor
191 if (uh(i,j,k) > 0.0) uh(i,j,k) = uh(i,j,k) * scale_factor
192 if (-vh(i,j-1,k) > 0.0) vh(i,j-1,k) = vh(i,j-1,k) * scale_factor
193 if (vh(i,j,k) > 0.0) vh(i,j,k) = vh(i,j,k) * scale_factor
194
195 ! Scale the flux across the interface atop a layer if it is upward
196 if (top_flux(i,j,k) > 0.0) then
197 ea(i,j,k) = ea(i,j,k) * scale_factor
198 if (k > 1) &
199 eb(i,j,k-1) = eb(i,j,k-1) * scale_factor
200 endif
201 ! Scale the flux across the interface atop a layer if it is downward
202 if (bottom_flux(i,j,k) > 0.0) then
203 eb(i,j,k) = eb(i,j,k) * scale_factor
204 if (k < nz) &
205 ea(i,j,k+1) = ea(i,j,k+1) * scale_factor
206 endif
207 enddo ; enddo ; enddo
208
209end subroutine limit_mass_flux_3d
210
211!> In the case where offline advection has failed to converge, redistribute the u-flux
212!! into remainder of the water column as a barotropic equivalent
213subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh)
214 type(ocean_grid_type), intent(in ) :: g !< ocean grid structure
215 type(verticalgrid_type), intent(in ) :: gv !< ocean vertical grid structure
216 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
217 intent(in ) :: hvol !< Mass of water in the cells at the end
218 !! of the previous timestep [H L2 ~> m3 or kg]
219 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
220 intent(inout) :: uh !< Zonal mass transport within a timestep [H L2 ~> m3 or kg]
221
222 ! Local variables
223 real, dimension(SZIB_(G),SZK_(GV)) :: uh2d ! A 2-d slice of transports [H L2 ~> m3 or kg]
224 real, dimension(SZIB_(G)) :: uh2d_sum ! Vertically summed transports [H L2 ~> m3 or kg]
225 real, dimension(SZI_(G),SZK_(GV)) :: h2d ! A 2-d slice of cell volumes [H L2 ~> m3 or kg]
226 real, dimension(SZI_(G)) :: h2d_sum ! Vertically summed cell volumes [H L2 ~> m3 or kg]
227
228 real :: abs_uh_sum ! The vertical sum of the absolute value of the transports [H L2 ~> m3 or kg]
229 real :: new_uh_sum ! The vertically summed transports after redistribution [H L2 ~> m3 or kg]
230 real :: uh_neglect ! A negligible transport [H L2 ~> m3 or kg]
231 integer :: i, j, k, is, ie, js, je, nz
232
233 ! Set index-related variables for fields on T-grid
234 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
235
236 do j=js,je
237 uh2d_sum(:) = 0.0
238 ! Copy over uh to a working array and sum up the remaining fluxes in a column
239 do k=1,nz ; do i=is-1,ie
240 uh2d(i,k) = uh(i,j,k)
241 uh2d_sum(i) = uh2d_sum(i) + uh2d(i,k)
242 enddo ; enddo
243
244 ! Copy over h to a working array and calculate total column volume
245 h2d_sum(:) = 0.0
246 do k=1,nz ; do i=is-1,ie+1
247 h2d(i,k) = hvol(i,j,k)
248 if (hvol(i,j,k)>0.) then
249 h2d_sum(i) = h2d_sum(i) + h2d(i,k)
250 else
251 h2d(i,k) = gv%H_subroundoff * g%areaT(i,j)
252 endif
253 enddo ; enddo
254
255 ! Distribute flux. Note min/max is intended to make sure that the mass transport
256 ! does not deplete a cell
257 do i=is-1,ie
258 if ( uh2d_sum(i)>0.0 ) then
259 do k=1,nz
260 uh2d(i,k) = uh2d_sum(i)*(h2d(i,k)/h2d_sum(i))
261 enddo
262 elseif (uh2d_sum(i)<0.0) then
263 do k=1,nz
264 uh2d(i,k) = uh2d_sum(i)*(h2d(i+1,k)/h2d_sum(i+1))
265 enddo
266 else
267 do k=1,nz
268 uh2d(i,k) = 0.0
269 enddo
270 endif
271
272 ! Check that column integrated transports match the original to within roundoff.
273 uh_neglect = gv%Angstrom_H * min(g%areaT(i,j), g%areaT(i+1,j))
274 abs_uh_sum = 0.0 ; new_uh_sum = 0.0
275 do k=1,nz
276 abs_uh_sum = abs_uh_sum + abs(uh2d(j,k))
277 new_uh_sum = new_uh_sum + uh2d(j,k)
278 enddo
279 if ( abs(new_uh_sum - uh2d_sum(j)) > max(uh_neglect, (5.0e-16*nz)*abs_uh_sum) ) &
280 call mom_error(warning, "Column integral of uh does not match after "//&
281 "barotropic redistribution")
282 enddo
283
284 do k=1,nz ; do i=is-1,ie
285 uh(i,j,k) = uh2d(i,k)
286 enddo ; enddo
287 enddo
288
290
291!> Redistribute the v-flux as a barotropic equivalent
292subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh)
293 type(ocean_grid_type), intent(in ) :: g !< ocean grid structure
294 type(verticalgrid_type), intent(in ) :: gv !< ocean vertical grid structure
295 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
296 intent(in ) :: hvol !< Mass of water in the cells at the end
297 !! of the previous timestep [H L2 ~> m3 or kg]
298 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
299 intent(inout) :: vh !< Meridional mass transport within a timestep [H L2 ~> m3 or kg]
300
301 ! Local variables
302 real, dimension(SZJB_(G),SZK_(GV)) :: vh2d ! A 2-d slice of transports [H L2 ~> m3 or kg]
303 real, dimension(SZJB_(G)) :: vh2d_sum ! Vertically summed transports [H L2 ~> m3 or kg]
304 real, dimension(SZJ_(G),SZK_(GV)) :: h2d ! A 2-d slice of cell volumes [H L2 ~> m3 or kg]
305 real, dimension(SZJ_(G)) :: h2d_sum ! Vertically summed cell volumes [H L2 ~> m3 or kg]
306
307 real :: abs_vh_sum ! The vertical sum of the absolute value of the transports [H L2 ~> m3 or kg]
308 real :: new_vh_sum ! The vertically summed transports after redistribution [H L2 ~> m3 or kg]
309 real :: vh_neglect ! A negligible transport [H L2 ~> m3 or kg]
310 integer :: i, j, k, is, ie, js, je, nz
311
312 ! Set index-related variables for fields on T-grid
313 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
314
315 do i=is,ie
316 vh2d_sum(:) = 0.0
317 ! Copy over uh to a working array and sum up the remaining fluxes in a column
318 do k=1,nz ; do j=js-1,je
319 vh2d(j,k) = vh(i,j,k)
320 vh2d_sum(j) = vh2d_sum(j) + vh2d(j,k)
321 enddo ; enddo
322
323 ! Copy over h to a working array and calculate column volume
324 h2d_sum(:) = 0.0
325 do k=1,nz ; do j=js-1,je+1
326 h2d(j,k) = hvol(i,j,k)
327 if (hvol(i,j,k)>0.) then
328 h2d_sum(j) = h2d_sum(j) + h2d(j,k)
329 else
330 h2d(i,k) = gv%H_subroundoff * g%areaT(i,j)
331 endif
332 enddo ; enddo
333
334 ! Distribute flux evenly throughout a column
335 do j=js-1,je
336 if ( vh2d_sum(j)>0.0 ) then
337 do k=1,nz
338 vh2d(j,k) = vh2d_sum(j)*(h2d(j,k)/h2d_sum(j))
339 enddo
340 elseif (vh2d_sum(j)<0.0) then
341 do k=1,nz
342 vh2d(j,k) = vh2d_sum(j)*(h2d(j+1,k)/h2d_sum(j+1))
343 enddo
344 else
345 do k=1,nz
346 vh2d(j,k) = 0.0
347 enddo
348 endif
349
350 ! Check that column integrated transports match the original to within roundoff.
351 vh_neglect = gv%Angstrom_H * min(g%areaT(i,j), g%areaT(i,j+1))
352 abs_vh_sum = 0.0 ; new_vh_sum = 0.0
353 do k=1,nz
354 abs_vh_sum = abs_vh_sum + abs(vh2d(j,k))
355 new_vh_sum = new_vh_sum + vh2d(j,k)
356 enddo
357 if ( abs(new_vh_sum - vh2d_sum(j)) > max(vh_neglect, (5.0e-16*nz)*abs_vh_sum) ) &
358 call mom_error(warning, "Column integral of vh does not match after "//&
359 "barotropic redistribution")
360 enddo
361
362 do k=1,nz ; do j=js-1,je
363 vh(i,j,k) = vh2d(j,k)
364 enddo ; enddo
365 enddo
366
368
369!> In the case where offline advection has failed to converge, redistribute the u-flux
370!! into layers above
371subroutine distribute_residual_uh_upwards(G, GV, hvol, uh)
372 type(ocean_grid_type), intent(in ) :: g !< ocean grid structure
373 type(verticalgrid_type), intent(in ) :: gv !< ocean vertical grid structure
374 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
375 intent(in ) :: hvol !< Mass of water in the cells at the end
376 !! of the previous timestep [H L2 ~> m3 or kg]
377 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
378 intent(inout) :: uh !< Zonal mass transport within a timestep [H L2 ~> m3 or kg]
379
380 ! Local variables
381 real, dimension(SZIB_(G),SZK_(GV)) :: uh2d ! A slice of transports [H L2 ~> m3 or kg]
382 real, dimension(SZI_(G),SZK_(GV)) :: h2d ! A slice of updated cell volumes [H L2 ~> m3 or kg]
383
384 real :: uh_neglect, uh_remain, uh_sum, uh_col ! Transports [H L2 ~> m3 or kg]
385 real :: hup, hlos ! Various cell volumes [H L2 ~> m3 or kg]
386 real :: min_h ! A minimal layer thickness [H ~> m or kg m-2]
387 integer :: i, j, k, is, ie, js, je, nz, k_rev
388
389 ! Set index-related variables for fields on T-grid
390 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
391
392 min_h = gv%Angstrom_H*0.1
393
394 do j=js,je
395 ! Copy over uh and cell volume to working arrays
396 do k=1,nz ; do i=is-2,ie+1
397 uh2d(i,k) = uh(i,j,k)
398 enddo ; enddo
399 do k=1,nz ; do i=is-1,ie+1
400 ! Subtract just a little bit of thickness to avoid roundoff errors
401 h2d(i,k) = hvol(i,j,k) - min_h * g%areaT(i,j)
402 enddo ; enddo
403
404 do i=is-1,ie
405 uh_col = sum(uh2d(i,:)) ! Store original column-integrated transport
406 do k=1,nz
407 uh_remain = uh2d(i,k)
408 uh2d(i,k) = 0.0
409 if (abs(uh_remain) > 0.0) then
410 do k_rev = k,1,-1
411 uh_sum = uh_remain + uh2d(i,k_rev)
412 if (uh_sum<0.0) then ! Transport to the left
413 hup = h2d(i+1,k_rev)
414 hlos = max(0.0,uh2d(i+1,k_rev))
415 if ((((hup - hlos) + uh_sum) < 0.0) .and. &
416 ((0.5*hup + uh_sum) < 0.0)) then
417 uh2d(i,k_rev) = min(-0.5*hup,-hup+hlos,0.0)
418 uh_remain = uh_sum - uh2d(i,k_rev)
419 else
420 uh2d(i,k_rev) = uh_sum
421 uh_remain = 0.0
422 exit
423 endif
424 else ! Transport to the right
425 hup = h2d(i,k_rev)
426 hlos = max(0.0,-uh2d(i-1,k_rev))
427 if ((((hup - hlos) - uh_sum) < 0.0) .and. &
428 ((0.5*hup - uh_sum) < 0.0)) then
429 uh2d(i,k_rev) = max(0.5*hup,hup-hlos,0.0)
430 uh_remain = uh_sum - uh2d(i,k_rev)
431 else
432 uh2d(i,k_rev) = uh_sum
433 uh_remain = 0.0
434 exit
435 endif
436 endif
437 enddo ! k_rev
438 endif
439
440 if (abs(uh_remain) > 0.0) then
441 if (k<nz) then
442 uh2d(i,k+1) = uh2d(i,k+1) + uh_remain
443 else
444 uh2d(i,k) = uh2d(i,k) + uh_remain
445 call mom_error(warning,"Water column cannot accommodate UH redistribution. Tracer may not be conserved")
446 endif
447 endif
448 enddo ! k-loop
449
450 ! Calculate and check that column integrated transports match the original to
451 ! within the tolerance limit
452 uh_neglect = gv%Angstrom_H * min(g%areaT(i,j), g%areaT(i+1,j))
453 if (abs(uh_col - sum(uh2d(i,:))) > uh_neglect) then
454 call mom_error(warning,"Column integral of uh does not match after upwards redistribution")
455 endif
456
457 enddo ! i-loop
458
459 do k=1,nz ; do i=is-1,ie
460 uh(i,j,k) = uh2d(i,k)
461 enddo ; enddo
462 enddo
463
465
466!> In the case where offline advection has failed to converge, redistribute the u-flux
467!! into layers above
468subroutine distribute_residual_vh_upwards(G, GV, hvol, vh)
469 type(ocean_grid_type), intent(in ) :: g !< ocean grid structure
470 type(verticalgrid_type), intent(in ) :: gv !< ocean vertical grid structure
471 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
472 intent(in ) :: hvol !< Mass of water in the cells at the end
473 !! of the previous timestep [H L2 ~> m3 or kg]
474 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
475 intent(inout) :: vh !< Meridional mass transport within a timestep [H L2 ~> m3 or kg]
476
477 ! Local variables
478 real, dimension(SZJB_(G),SZK_(GV)) :: vh2d ! A slice of transports [H L2 ~> m3 or kg]
479 real, dimension(SZJ_(G),SZK_(GV)) :: h2d ! A slice of updated cell volumes [H L2 ~> m3 or kg]
480
481 real :: vh_neglect, vh_remain, vh_col, vh_sum ! Transports [H L2 ~> m3 or kg]
482 real :: hup, hlos ! Various cell volumes [H L2 ~> m3 or kg]
483 real :: min_h ! A minimal layer thickness [H ~> m or kg m-2]
484 integer :: i, j, k, is, ie, js, je, nz, k_rev
485
486 ! Set index-related variables for fields on T-grid
487 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
488
489 min_h = 0.1*gv%Angstrom_H
490
491 do i=is,ie
492 ! Copy over uh and cell volume to working arrays
493 do k=1,nz ; do j=js-2,je+1
494 vh2d(j,k) = vh(i,j,k)
495 enddo ; enddo
496 do k=1,nz ; do j=js-1,je+1
497 h2d(j,k) = hvol(i,j,k) - min_h * g%areaT(i,j)
498 enddo ; enddo
499
500 do j=js-1,je
501 vh_col = sum(vh2d(j,:))
502 do k=1,nz
503 vh_remain = vh2d(j,k)
504 vh2d(j,k) = 0.0
505 if (abs(vh_remain) > 0.0) then
506 do k_rev = k,1,-1
507 vh_sum = vh_remain + vh2d(j,k_rev)
508 if (vh_sum<0.0) then ! Transport to the left
509 hup = h2d(j+1,k_rev)
510 hlos = max(0.0,vh2d(j+1,k_rev))
511 if ((((hup - hlos) + vh_sum) < 0.0) .and. &
512 ((0.5*hup + vh_sum) < 0.0)) then
513 vh2d(j,k_rev) = min(-0.5*hup,-hup+hlos,0.0)
514 vh_remain = vh_sum - vh2d(j,k_rev)
515 else
516 vh2d(j,k_rev) = vh_sum
517 vh_remain = 0.0
518 exit
519 endif
520 else ! Transport to the right
521 hup = h2d(j,k_rev)
522 hlos = max(0.0,-vh2d(j-1,k_rev))
523 if ((((hup - hlos) - vh_sum) < 0.0) .and. &
524 ((0.5*hup - vh_sum) < 0.0)) then
525 vh2d(j,k_rev) = max(0.5*hup,hup-hlos,0.0)
526 vh_remain = vh_sum - vh2d(j,k_rev)
527 else
528 vh2d(j,k_rev) = vh_sum
529 vh_remain = 0.0
530 exit
531 endif
532 endif
533
534 enddo ! k_rev
535 endif
536
537 if (abs(vh_remain) > 0.0) then
538 if (k<nz) then
539 vh2d(j,k+1) = vh2d(j,k+1) + vh_remain
540 else
541 vh2d(j,k) = vh2d(j,k) + vh_remain
542 call mom_error(warning,"Water column cannot accommodate VH redistribution. Tracer will not be conserved")
543 endif
544 endif ! k-loop
545 enddo
546
547 ! Calculate and check that column integrated transports match the original to
548 ! within the tolerance limit
549 vh_neglect = gv%Angstrom_H * min(g%areaT(i,j), g%areaT(i,j+1))
550 if ( abs(vh_col-sum(vh2d(j,:))) > vh_neglect) then
551 call mom_error(warning,"Column integral of vh does not match after "//&
552 "upwards redistribution")
553 endif
554 enddo
555
556 do k=1,nz ; do j=js-1,je
557 vh(i,j,k) = vh2d(j,k)
558 enddo ; enddo
559 enddo
560
562
563!> add_diurnal_SW adjusts the shortwave fluxes in an forcying_type variable
564!! to add a synthetic diurnal cycle. Adapted from SIS2
565subroutine offline_add_diurnal_sw(fluxes, G, Time_start, Time_end)
566 type(forcing), intent(inout) :: fluxes !< The type with atmospheric fluxes to be adjusted.
567 type(ocean_grid_type), intent(in) :: g !< The ocean lateral grid type.
568 type(time_type), intent(in) :: time_start !< The start time for this step.
569 type(time_type), intent(in) :: time_end !< The ending time for this step.
570
571 real :: diurnal_factor ! A scaling factor to insert a synthetic diurnal cycle [nondim]
572 real :: time_since_ae ! Time since the autumnal equinox expressed as a fraction of a year times 2 pi [nondim]
573 real :: rad ! A conversion factor from degrees to radians = pi/180 degrees [nondim]
574 real :: fracday_dt ! Daylight fraction averaged over a timestep [nondim]
575 real :: fracday_day ! Daylight fraction averaged over a day [nondim]
576 real :: cosz_day ! Cosine of the solar zenith angle averaged over a day [nondim]
577 real :: cosz_dt ! Cosine of the solar zenith angle averaged over a timestep [nondim]
578 real :: rrsun_day ! Earth-Sun distance (r) relative to the semi-major axis of
579 ! the orbital ellipse averaged over a day [nondim]
580 real :: rrsun_dt ! Earth-Sun distance (r) relative to the semi-major axis of
581 ! the orbital ellipse averaged over a timestep [nondim]
582 type(time_type) :: dt_here ! The time increment covered by this call
583
584 integer :: i, j, i2, j2, isc, iec, jsc, jec, i_off, j_off
585
586 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
587 i_off = lbound(fluxes%sens,1) - g%isc ; j_off = lbound(fluxes%sens,2) - g%jsc
588
589 ! Orbital_time extracts the time of year relative to the northern
590 ! hemisphere autumnal equinox from a time_type variable.
591 time_since_ae = orbital_time(time_start)
592 dt_here = time_end - time_start
593 rad = acos(-1.)/180.
594
595 !$OMP parallel do default(shared) private(i,j,i2,j2,cosz_dt,fracday_dt,rrsun_dt, &
596 !$OMP fracday_day,cosz_day,rrsun_day,diurnal_factor)
597 do j=jsc,jec ; do i=isc,iec
598! Per Rick Hemler:
599! Call diurnal_solar with dtime=dt_here to get cosz averaged over dt_here.
600! Call daily_mean_solar to get cosz averaged over a day. Then
601! diurnal_factor = cosz_dt_ice*fracday_dt_ice*rrsun_dt_ice /
602! cosz_day*fracday_day*rrsun_day
603
604 call diurnal_solar(g%geoLatT(i,j)*rad, g%geoLonT(i,j)*rad, time_start, cosz=cosz_dt, &
605 fracday=fracday_dt, rrsun=rrsun_dt, dt_time=dt_here)
606 call daily_mean_solar(g%geoLatT(i,j)*rad, time_since_ae, cosz_day, fracday_day, rrsun_day)
607 diurnal_factor = cosz_dt*fracday_dt*rrsun_dt / &
608 max(1e-30, cosz_day*fracday_day*rrsun_day)
609
610 i2 = i+i_off ; j2 = j+j_off
611 fluxes%sw(i2,j2) = fluxes%sw(i2,j2) * diurnal_factor
612 fluxes%sw_vis_dir(i2,j2) = fluxes%sw_vis_dir(i2,j2) * diurnal_factor
613 fluxes%sw_vis_dif(i2,j2) = fluxes%sw_vis_dif(i2,j2) * diurnal_factor
614 fluxes%sw_nir_dir(i2,j2) = fluxes%sw_nir_dir(i2,j2) * diurnal_factor
615 fluxes%sw_nir_dif(i2,j2) = fluxes%sw_nir_dif(i2,j2) * diurnal_factor
616 enddo ; enddo
617
618end subroutine offline_add_diurnal_sw
619
620!> Controls the reading in 3d mass fluxes, diffusive fluxes, and other fields stored
621!! in a previous integration of the online model
622subroutine update_offline_from_files(G, GV, US, nk_input, mean_file, sum_file, snap_file, &
623 surf_file, h_end, uhtr, vhtr, temp_mean, salt_mean, mld, Kd, fluxes, &
624 ridx_sum, ridx_snap, read_mld, read_sw, read_ts_uvh, do_ale_in)
625
626 type(ocean_grid_type), intent(inout) :: g !< Horizontal grid type
627 type(verticalgrid_type), intent(in ) :: gv !< Vertical grid type
628 type(unit_scale_type), intent(in ) :: us !< A dimensional unit scaling type
629 integer, intent(in ) :: nk_input !< Number of levels in input file
630 character(len=*), intent(in ) :: mean_file !< Name of file with averages fields
631 character(len=*), intent(in ) :: sum_file !< Name of file with summed fields
632 character(len=*), intent(in ) :: snap_file !< Name of file with snapshot fields
633 character(len=*), intent(in ) :: surf_file !< Name of file with surface fields
634 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
635 intent(inout) :: h_end !< End of timestep layer thickness [H ~> m or kg m-2]
636 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
637 intent(inout) :: uhtr !< Zonal mass fluxes [H L2 ~> m3 or kg]
638 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
639 intent(inout) :: vhtr !< Meridional mass fluxes [H L2 ~> m3 or kg]
640 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
641 intent(inout) :: temp_mean !< Averaged temperature [C ~> degC]
642 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
643 intent(inout) :: salt_mean !< Averaged salinity [S ~> ppt]
644 real, dimension(SZI_(G),SZJ_(G)), &
645 intent(inout) :: mld !< Averaged mixed layer depth [Z ~> m]
646 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
647 intent(inout) :: kd !< Diapycnal diffusivities at interfaces
648 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
649 type(forcing), intent(inout) :: fluxes !< Fields with surface fluxes
650 integer, intent(in ) :: ridx_sum !< Read index for sum, mean, and surf files
651 integer, intent(in ) :: ridx_snap !< Read index for snapshot file
652 logical, intent(in ) :: read_mld !< True if reading in MLD
653 logical, intent(in ) :: read_sw !< True if reading in radiative fluxes
654 logical, intent(in ) :: read_ts_uvh !< True if reading in uh, vh, and h
655 logical, optional, intent(in ) :: do_ale_in !< True if using ALE algorithms
656
657 logical :: do_ale
658 real :: convert_to_h ! A scale conversion factor from the thickness units in the
659 ! file to H [H m-1 ~> 1] or [H m2 kg-1 ~> 1]
660 integer :: i, j, k, is, ie, js, je, nz
661
662 do_ale = .false.
663 if (present(do_ale_in)) do_ale = do_ale_in
664
665 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
666
667 if (gv%Boussinesq) then
668 convert_to_h = gv%m_to_H
669 else
670 convert_to_h = gv%kg_m2_to_H
671 endif
672
673 ! Check if reading in temperature, salinity, transports and ending thickness
674 if (read_ts_uvh) then
675 h_end(:,:,:) = 0.0
676 temp_mean(:,:,:) = 0.0
677 salt_mean(:,:,:) = 0.0
678 uhtr(:,:,:) = 0.0
679 vhtr(:,:,:) = 0.0
680 ! Time-summed fields
681 call mom_read_vector(sum_file, 'uhtr_sum', 'vhtr_sum', uhtr(:,:,1:nk_input), &
682 vhtr(:,:,1:nk_input), g%Domain, timelevel=ridx_sum, &
683 scale=us%m_to_L**2*gv%kg_m2_to_H)
684 call mom_read_data(snap_file, 'h_end', h_end(:,:,1:nk_input), g%Domain, &
685 timelevel=ridx_snap, position=center, scale=convert_to_h)
686 call mom_read_data(mean_file, 'temp', temp_mean(:,:,1:nk_input), g%Domain, &
687 timelevel=ridx_sum, position=center, scale=us%degC_to_C)
688 call mom_read_data(mean_file, 'salt', salt_mean(:,:,1:nk_input), g%Domain, &
689 timelevel=ridx_sum, position=center, scale=us%ppt_to_S)
690
691 ! Fill temperature and salinity downward from the deepest input data.
692 do k=nk_input+1,nz ; do j=js,je ; do i=is,ie
693 if (g%mask2dT(i,j)>0.) then
694 temp_mean(i,j,k) = temp_mean(i,j,nk_input)
695 salt_mean(i,j,k) = salt_mean(i,j,nk_input)
696 endif
697 enddo ; enddo ; enddo
698 endif
699
700 ! Check if reading vertical diffusivities or entrainment fluxes
701 call mom_read_data( mean_file, 'Kd_interface', kd(:,:,1:nk_input+1), g%Domain, &
702 timelevel=ridx_sum, position=center, scale=gv%m2_s_to_HZ_T)
703
704 ! This block makes sure that the fluxes control structure, which may not be used in the solo_driver,
705 ! contains netMassIn and netMassOut which is necessary for the applyTracerBoundaryFluxesInOut routine
706 if (do_ale) then
707 if (.not. associated(fluxes%netMassOut)) &
708 allocate(fluxes%netMassOut(g%isd:g%ied,g%jsd:g%jed), source=0.0)
709 if (.not. associated(fluxes%netMassIn)) &
710 allocate(fluxes%netMassIn(g%isd:g%ied,g%jsd:g%jed), source=0.0)
711
712 fluxes%netMassOut(:,:) = 0.0
713 fluxes%netMassIn(:,:) = 0.0
714 call mom_read_data(surf_file,'massout_flux_sum',fluxes%netMassOut, g%Domain, &
715 timelevel=ridx_sum, scale=gv%kg_m2_to_H)
716 call mom_read_data(surf_file,'massin_flux_sum', fluxes%netMassIn, g%Domain, &
717 timelevel=ridx_sum, scale=gv%kg_m2_to_H)
718
719 do j=js,je ; do i=is,ie
720 if (g%mask2dT(i,j)<1.0) then
721 fluxes%netMassOut(i,j) = 0.0
722 fluxes%netMassIn(i,j) = 0.0
723 endif
724 enddo ; enddo
725
726 endif
727
728 if (read_mld) then
729 call mom_read_data(surf_file, 'ePBL_h_ML', mld, g%Domain, timelevel=ridx_sum, scale=us%m_to_Z)
730 endif
731
732 if (read_sw) then
733 ! Shortwave radiation is only needed for offline mode with biogeochemistry but without the coupler.
734 ! Need to double check, but set_opacity seems to only need the sum of the diffuse and
735 ! direct fluxes in the visible and near-infrared bands. For convenience, we store the
736 ! sum of the direct and diffuse fluxes in the 'dir' field and set the 'dif' fields to zero
737 call mom_read_data(mean_file,'sw_vis', fluxes%sw_vis_dir, g%Domain, &
738 timelevel=ridx_sum, scale=us%W_m2_to_QRZ_T)
739 call mom_read_data(mean_file,'sw_nir', fluxes%sw_nir_dir, g%Domain, &
740 timelevel=ridx_sum, scale=us%W_m2_to_QRZ_T)
741 fluxes%sw_vis_dir(:,:) = fluxes%sw_vis_dir(:,:)*0.5
742 fluxes%sw_vis_dif(:,:) = fluxes%sw_vis_dir(:,:)
743 fluxes%sw_nir_dir(:,:) = fluxes%sw_nir_dir(:,:)*0.5
744 fluxes%sw_nir_dif(:,:) = fluxes%sw_nir_dir(:,:)
745 fluxes%sw = (fluxes%sw_vis_dir + fluxes%sw_vis_dif) + (fluxes%sw_nir_dir + fluxes%sw_nir_dif)
746 do j=js,je ; do i=is,ie
747 if (g%mask2dT(i,j)<1.0) then
748 fluxes%sw(i,j) = 0.0
749 fluxes%sw_vis_dir(i,j) = 0.0
750 fluxes%sw_nir_dir(i,j) = 0.0
751 fluxes%sw_vis_dif(i,j) = 0.0
752 fluxes%sw_nir_dif(i,j) = 0.0
753 endif
754 enddo ; enddo
755 call pass_var(fluxes%sw,g%Domain)
756 call pass_var(fluxes%sw_vis_dir,g%Domain)
757 call pass_var(fluxes%sw_vis_dif,g%Domain)
758 call pass_var(fluxes%sw_nir_dir,g%Domain)
759 call pass_var(fluxes%sw_nir_dif,g%Domain)
760 endif
761
762end subroutine update_offline_from_files
763
764!> Fields for offline transport are copied from the stored arrays read during initialization
765subroutine update_offline_from_arrays(G, GV, nk_input, ridx_sum, mean_file, sum_file, snap_file, uhtr, vhtr, &
766 hend, uhtr_all, vhtr_all, hend_all, temp, salt, temp_all, salt_all )
767 type(ocean_grid_type), intent(inout) :: g !< Horizontal grid type
768 type(verticalgrid_type), intent(in ) :: gv !< Vertical grid type
769 integer, intent(in ) :: nk_input !< Number of levels in input file
770 integer, intent(in ) :: ridx_sum !< Index to read from
771 character(len=200), intent(in ) :: mean_file !< Name of file with averages fields
772 character(len=200), intent(in ) :: sum_file !< Name of file with summed fields
773 character(len=200), intent(in ) :: snap_file !< Name of file with snapshot fields
774 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Zonal mass fluxes [H L2 ~> m3 or kg]
775 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Meridional mass fluxes [H L2 ~> m3 or kg]
776 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hend !< End of timestep layer thickness
777 !! [H ~> m or kg m-2]
778 real, dimension(:,:,:,:), allocatable, intent(inout) :: uhtr_all !< Zonal mass fluxes [H L2 ~> m3 or kg]
779 real, dimension(:,:,:,:), allocatable, intent(inout) :: vhtr_all !< Meridional mass fluxes [H L2 ~> m3 or kg]
780 real, dimension(:,:,:,:), allocatable, intent(inout) :: hend_all !< End of timestep layer thickness
781 !! [H ~> m or kg m-2]
782 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: temp !< Temperature array [C ~> degC]
783 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: salt !< Salinity array [S ~> ppt]
784 real, dimension(:,:,:,:), allocatable, intent(inout) :: temp_all !< Temperature array [C ~> degC]
785 real, dimension(:,:,:,:), allocatable, intent(inout) :: salt_all !< Salinity array [S ~> ppt]
786
787 integer :: i, j, k, is, ie, js, je, nz
788 real, parameter :: fill_value = 0. ! The fill value for input arrays [various]
789 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
790
791 ! Check that all fields are allocated (this is a redundant check)
792 if (.not. allocated(uhtr_all)) &
793 call mom_error(fatal, "uhtr_all not allocated before call to update_transport_from_arrays")
794 if (.not. allocated(vhtr_all)) &
795 call mom_error(fatal, "vhtr_all not allocated before call to update_transport_from_arrays")
796 if (.not. allocated(hend_all)) &
797 call mom_error(fatal, "hend_all not allocated before call to update_transport_from_arrays")
798 if (.not. allocated(temp_all)) &
799 call mom_error(fatal, "temp_all not allocated before call to update_transport_from_arrays")
800 if (.not. allocated(salt_all)) &
801 call mom_error(fatal, "salt_all not allocated before call to update_transport_from_arrays")
802
803 ! Copy uh, vh, h_end, temp, and salt
804 do k=1,nk_input ; do j=js,je ; do i=is,ie
805 uhtr(i,j,k) = uhtr_all(i,j,k,ridx_sum)
806 vhtr(i,j,k) = vhtr_all(i,j,k,ridx_sum)
807 hend(i,j,k) = hend_all(i,j,k,ridx_sum)
808 temp(i,j,k) = temp_all(i,j,k,ridx_sum)
809 salt(i,j,k) = salt_all(i,j,k,ridx_sum)
810 enddo ; enddo ; enddo
811
812 ! Fill the rest of the arrays with 0s (fill_value could probably be changed to a runtime parameter)
813 do k=nk_input+1,nz ; do j=js,je ; do i=is,ie
814 uhtr(i,j,k) = fill_value
815 vhtr(i,j,k) = fill_value
816 hend(i,j,k) = fill_value
817 temp(i,j,k) = fill_value
818 salt(i,j,k) = fill_value
819 enddo ; enddo ; enddo
820
821end subroutine update_offline_from_arrays
822
823!> Calculates the next timelevel to read from the input fields. This allows the 'looping'
824!! of the fields
825function next_modulo_time(inidx, numtime)
826 ! Returns the next time interval to be read
827 integer :: numtime ! Number of time levels in input fields
828 integer :: inidx ! The current time index
829
830 integer :: read_index ! The index in the input files that corresponds
831 ! to the current timestep
832
833 integer :: next_modulo_time
834
835 read_index = mod(inidx+1,numtime)
836 if (read_index < 0) read_index = inidx-read_index
837 if (read_index == 0) read_index = numtime
838
839 next_modulo_time = read_index
840
841end function next_modulo_time
842
843end module mom_offline_aux
844