MOM_offline_main.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!> The routines here implement the offline tracer algorithm used in MOM6. These are called from step_offline
6!! Some routines called here can be found in the MOM_offline_aux module.
8
11use mom_ale, only : ale_remap_tracers
12use mom_checksums, only : hchksum, uvchksum
13use mom_coms, only : reproducing_sum
14use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
15use mom_cpu_clock, only : clock_component, clock_subcomponent
16use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
21use mom_domains, only : pass_var, pass_vector
22use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
24use mom_file_parser, only : read_param, get_param, log_version, param_file_type
25use mom_forcing_type, only : forcing
26use mom_grid, only : ocean_grid_type
27use mom_interface_heights, only : calc_derived_thermo, thickness_to_dz
28use mom_io, only : mom_read_data, mom_read_vector, center
36use mom_time_manager, only : time_type, real_to_time
40use mom_tracer_registry, only : tracer_registry_type, mom_tracer_chksum, mom_tracer_chkinv
44
45implicit none ; private
46
47#include "MOM_memory.h"
48
49!> The control structure for the offline transport module
50type, public :: offline_transport_cs ; private
51
52 ! Pointers to relevant fields from the main MOM control structure
53 type(ale_cs), pointer :: ale_csp => null()
54 !< A pointer to the ALE control structure
55 type(diabatic_cs), pointer :: diabatic_csp => null()
56 !< A pointer to the diabatic control structure
57 type(diag_ctrl), pointer :: diag => null()
58 !< Structure that regulates diagnostic output
59 type(ocean_obc_type), pointer :: obc => null()
60 !< A pointer to the open boundary condition control structure
61 type(tracer_advect_cs), pointer :: tracer_adv_csp => null()
62 !< A pointer to the tracer advection control structure
63 type(opacity_cs), pointer :: opacity_csp => null()
64 !< A pointer to the opacity control structure
65 type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null()
66 !< A pointer to control structure that orchestrates the calling of tracer packages
67 type(tracer_registry_type), pointer :: tracer_reg => null()
68 !< A pointer to the tracer registry
69 type(thermo_var_ptrs), pointer :: tv => null()
70 !< A structure pointing to various thermodynamic variables
71 type(optics_type), pointer :: optics => null()
72 !< Pointer to the optical properties type
73 type(diabatic_aux_cs), pointer :: diabatic_aux_csp => null()
74 !< Pointer to the diabatic_aux control structure
75
76 !> Variables related to reading in fields from online run
77 integer :: start_index !< Timelevel to start
78 integer :: iter_no !< Timelevel to start
79 integer :: numtime !< How many timelevels in the input fields
80 type(time_type) :: accumulated_time !< Length of time accumulated in the current offline interval
81 type(time_type) :: vertical_time !< The next value of accumulate_time at which to apply vertical processes
82 ! Index of each of the variables to be read in with separate indices for each variable if they
83 ! are set off from each other in time
84 integer :: ridx_sum = -1 !< Read index offset of the summed variables
85 integer :: ridx_snap = -1 !< Read index offset of the snapshot variables
86 integer :: nk_input !< Number of input levels in the input fields
87 character(len=200) :: offlinedir !< Directory where offline fields are stored
88 character(len=200) :: & ! Names of input files
89 surf_file, & !< Contains surface fields (2d arrays)
90 snap_file, & !< Snapshotted fields (layer thicknesses)
91 sum_file, & !< Fields which are accumulated over time
92 mean_file !< Fields averaged over time
93 character(len=20) :: redistribute_method !< 'barotropic' if evenly distributing extra flow
94 !! throughout entire watercolumn, 'upwards',
95 !! if trying to do it just in the layers above
96 !! 'both' if both methods are used
97 character(len=20) :: mld_var_name !< Name of the mixed layer depth variable to use
98 logical :: fields_are_offset !< True if the time-averaged fields and snapshot fields are
99 !! offset by one time level
100 logical :: x_before_y !< Which horizontal direction is advected first
101 logical :: print_adv_offline !< Prints out some updates each advection sub interation
102 logical :: skip_diffusion !< Skips horizontal diffusion of tracers
103 logical :: read_sw !< Read in averaged values for shortwave radiation
104 logical :: read_mld !< Check to see whether mixed layer depths should be read in
105 logical :: diurnal_sw !< Adds a synthetic diurnal cycle on shortwave radiation
106 logical :: debug !< If true, write verbose debugging messages
107 logical :: redistribute_barotropic !< Redistributes column-summed residual transports throughout
108 !! a column weighted by thickness
109 logical :: redistribute_upwards !< Redistributes remaining fluxes only in layers above
110 !! the current one based as the max allowable transport
111 !! in that cell
112 logical :: read_all_ts_uvh !< If true, then all timelevels of temperature, salinity, mass transports, and
113 !! Layer thicknesses are read during initialization
114 !! Variables controlling some of the numerical considerations of offline transport
115 integer :: num_off_iter !< Number of advection iterations per offline step
116 integer :: num_vert_iter !< Number of vertical iterations per offline step
117 integer :: off_ale_mod !< Sets how frequently the ALE step is done during the advection
118 real :: dt_offline !< Timestep used for offline tracers [T ~> s]
119 real :: dt_offline_vertical !< Timestep used for calls to tracer vertical physics [T ~> s]
120 real :: evap_cfl_limit !< Limit on the fraction of the water that can be fluxed out of the top
121 !! layer in a timestep [nondim]. This is Copied from diabatic_CS controlling
122 !! how tracers follow freshwater fluxes
123 real :: minimum_forcing_depth !< The smallest depth over which fluxes can be applied [H ~> m or kg m-2].
124 !! This is copied from diabatic_CS controlling how tracers follow freshwater fluxes
125
126 real :: kd_max !< Runtime parameter specifying the maximum value of vertical diffusivity
127 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
128 real :: min_residual !< The minimum amount of total mass flux before exiting the main advection
129 !! routine [H L2 ~> m3 or kg]
130 !>@{ Diagnostic manager IDs for some fields that may be of interest when doing offline transport
131 integer :: &
132 id_uhr = -1, &
133 id_vhr = -1, &
134 id_ear = -1, &
135 id_ebr = -1, &
136 id_hr = -1, &
137 id_hdiff = -1, &
138 id_uhr_redist = -1, &
139 id_vhr_redist = -1, &
140 id_uhr_end = -1, &
141 id_vhr_end = -1, &
142 id_eta_pre_distribute = -1, &
143 id_eta_post_distribute = -1, &
144 id_h_redist = -1, &
145 id_eta_diff_end = -1
146
147 ! Diagnostic IDs for the regridded/remapped input fields
148 integer :: &
149 id_uhtr_regrid = -1, &
150 id_vhtr_regrid = -1, &
151 id_temp_regrid = -1, &
152 id_salt_regrid = -1, &
153 id_h_regrid = -1
154 !>@}
155
156 ! IDs for timings of various offline components
157 integer :: id_clock_read_fields = -1 !< A CPU time clock
158 integer :: id_clock_offline_diabatic = -1 !< A CPU time clock
159 integer :: id_clock_offline_adv = -1 !< A CPU time clock
160 integer :: id_clock_redistribute = -1 !< A CPU time clock
161
162 !> Zonal transport that may need to be stored between calls to step_MOM [H L2 ~> m3 or kg]
163 real, allocatable, dimension(:,:,:) :: uhtr
164 !> Meridional transport that may need to be stored between calls to step_MOM [H L2 ~> m3 or kg]
165 real, allocatable, dimension(:,:,:) :: vhtr
166
167 ! Fields at T-point
168 real, allocatable, dimension(:,:,:) :: eatr
169 !< Amount of fluid entrained from the layer above within
170 !! one time step [H ~> m or kg m-2]
171 real, allocatable, dimension(:,:,:) :: ebtr
172 !< Amount of fluid entrained from the layer below within
173 !! one time step [H ~> m or kg m-2]
174 ! Fields at T-points on interfaces
175 real, allocatable, dimension(:,:,:) :: kd !< Vertical diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
176 real, allocatable, dimension(:,:,:) :: h_end !< Thicknesses at the end of offline timestep [H ~> m or kg m-2]
177
178 real, allocatable, dimension(:,:) :: mld !< Mixed layer depths at thickness points [Z ~> m]
179
180 ! Allocatable arrays to read in entire fields during initialization
181 real, allocatable, dimension(:,:,:,:) :: uhtr_all !< Entire field of zonal transport [H L2 ~> m3 or kg]
182 real, allocatable, dimension(:,:,:,:) :: vhtr_all !< Entire field of meridional transport [H L2 ~> m3 or kg]
183 real, allocatable, dimension(:,:,:,:) :: hend_all !< Entire field of layer thicknesses [H ~> m or kg m-2]
184 real, allocatable, dimension(:,:,:,:) :: temp_all !< Entire field of temperatures [C ~> degC]
185 real, allocatable, dimension(:,:,:,:) :: salt_all !< Entire field of salinities [S ~> ppt]
186
188
202
203contains
204
205!> 3D advection is done by doing flux-limited nonlinear horizontal advection interspersed with an ALE
206!! regridding/remapping step. The loop in this routine is exited if remaining residual transports are below
207!! a runtime-specified value or a maximum number of iterations is reached.
208subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, CS, id_clock_ale, &
209 h_pre, uhtr, vhtr, converged)
210 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
211 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
212 real, intent(in) :: time_interval !< time interval covered by this call [T ~> s]
213 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
214 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
215 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
216 type(offline_transport_cs), pointer :: cs !< control structure for offline module
217 integer, intent(in) :: id_clock_ale !< Clock for ALE routines
218 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
219 intent(inout) :: h_pre !< layer thicknesses before advection
220 !! [H ~> m or kg m-2]
221 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
222 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
223 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
224 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
225 logical, intent( out) :: converged !< True if the iterations have converged
226
227 ! Local variables
228 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhtr_sub ! Substep zonal mass transports [H L2 ~> m3 or kg]
229 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub ! Substep meridional mass transports [H L2 ~> m3 or kg]
230
231 real :: prev_tot_residual, tot_residual ! Used to keep track of how close to convergence we are [H L2 ~> m3 or kg]
232
233 ! Variables used to keep track of layer thicknesses at various points in the code
234 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
235 h_new, & ! Updated layer thicknesses [H ~> m or kg m-2]
236 h_post_remap, & ! Layer thicknesses after remapping [H ~> m or kg m-2]
237 h_vol ! Layer volumes [H L2 ~> m3 or kg]
238 real :: dzregrid(szi_(g),szj_(g),szk_(gv)+1) ! The change in grid interface positions due to regridding,
239 ! in the same units as thicknesses [H ~> m or kg m-2]
240 integer :: niter, iter
241 real :: inum_iter ! The inverse of the number of iterations [nondim]
242 character(len=256) :: mesg ! The text of an error message
243 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
244 integer :: isdb, iedb, jsdb, jedb
245 logical :: x_before_y
246 real :: evap_cfl_limit ! Limit on the fraction of the water that can be fluxed out of the
247 ! top layer in a timestep [nondim]
248 real :: minimum_forcing_depth ! The smallest depth over which fluxes can be applied [H ~> m or kg m-2]
249 real :: dt_iter ! The timestep to use for each iteration [T ~> s]
250 real :: hl2_to_kg_scale ! Unit conversion factors to cell mass [kg H-1 L-2 ~> kg m-3 or 1]
251 character(len=20) :: debug_msg
252 call cpu_clock_begin(cs%id_clock_offline_adv)
253
254 ! Grid-related pointer assignments
255
256 x_before_y = cs%x_before_y
257
258 ! Initialize some shorthand variables from other structures
259 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
260 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
261 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
262
263 evap_cfl_limit = cs%evap_CFL_limit
264 minimum_forcing_depth = cs%minimum_forcing_depth
265
266 hl2_to_kg_scale = us%L_to_m**2*gv%H_to_kg_m2
267 niter = cs%num_off_iter
268 inum_iter = 1./real(niter)
269 dt_iter = cs%dt_offline*inum_iter
270
271 ! Initialize working arrays
272 h_new(:,:,:) = 0.0
273 h_vol(:,:,:) = 0.0
274 uhtr_sub(:,:,:) = 0.0
275 vhtr_sub(:,:,:) = 0.0
276
277 ! converged should only be true if there are no remaining mass fluxes
278 converged = .false.
279
280 ! Tracers are transported using the stored mass fluxes. Where possible, operators are Strang-split around
281 ! the call to
282 ! 1) Using the layer thicknesses and tracer concentrations from the previous timestep,
283 ! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to tracer_column_fns.
284 ! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline
285 ! 2) Half of the accumulated surface freshwater fluxes are applied
286 !! START ITERATION
287 ! 3) Accumulated mass fluxes are used to do horizontal transport. The number of iterations used in
288 ! advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are stored for later use
289 ! and resulting layer thicknesses fed into the next step
290 ! 4) Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for layers which might
291 ! 'vanish' because of horizontal mass transport to be 'reinflated'
292 ! 5) Check that transport is done if the remaining mass fluxes equals 0 or if the max number of iterations
293 ! has been reached
294 !! END ITERATION
295 ! 6) Repeat steps 1 and 2
296 ! 7) Force a remapping to the stored layer thicknesses that correspond to the snapshot of the online model
297 ! 8) Reset T/S and h to their stored snapshotted values to prevent model drift
298
299 ! Copy over the horizontal mass fluxes from the total mass fluxes
300 do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
301 uhtr_sub(i,j,k) = uhtr(i,j,k)
302 enddo ; enddo ; enddo
303 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
304 vhtr_sub(i,j,k) = vhtr(i,j,k)
305 enddo ; enddo ; enddo
306 do k=1,nz ; do j=js,je ; do i=is,ie
307 h_new(i,j,k) = h_pre(i,j,k)
308 enddo ; enddo ; enddo
309
310 if (cs%debug) then
311 call hchksum(h_pre, "h_pre before transport", g%HI, unscale=gv%H_to_MKS)
312 call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, g%HI, unscale=hl2_to_kg_scale)
313 endif
314 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
315 if (cs%print_adv_offline) then
316 write(mesg,'(A,ES24.16)') "Main advection starting transport: ", tot_residual*hl2_to_kg_scale
317 call mom_mesg(mesg)
318 endif
319
320 ! This loop does essentially a flux-limited, nonlinear advection scheme until all mass fluxes
321 ! are used. ALE is done after the horizontal advection.
322 do iter=1,cs%num_off_iter
323
324 do k=1,nz ; do j=js,je ; do i=is,ie
325 h_vol(i,j,k) = h_new(i,j,k) * g%areaT(i,j)
326 h_pre(i,j,k) = h_new(i,j,k)
327 enddo ; enddo ; enddo
328
329 if (cs%debug) then
330 call hchksum(h_vol, "h_vol before advect", g%HI, unscale=hl2_to_kg_scale)
331 call uvchksum("[uv]htr_sub before advect", uhtr_sub, vhtr_sub, g%HI, unscale=hl2_to_kg_scale)
332 write(debug_msg, '(A,I4.4)') 'Before advect ', iter
333 call mom_tracer_chkinv(debug_msg, g, gv, h_pre, cs%tracer_reg)
334 endif
335
336 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, cs%dt_offline, g, gv, us, &
337 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
338 max_iter_in=1, update_vol_prev=.true., uhr_out=uhtr, vhr_out=vhtr)
339
340 ! Switch the direction every iteration
341 x_before_y = .not. x_before_y
342
343 ! Update the new layer thicknesses after one round of advection has happened
344 do k=1,nz ; do j=js,je ; do i=is,ie
345 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
346 enddo ; enddo ; enddo
347
348 if (modulo(iter,cs%off_ale_mod)==0) then
349 ! Do ALE remapping/regridding to allow for more advection to occur in the next iteration
350 call pass_var(h_new,g%Domain)
351 if (cs%debug) then
352 call hchksum(h_new,"h_new before ALE", g%HI, unscale=gv%H_to_MKS)
353 write(debug_msg, '(A,I4.4)') 'Before ALE ', iter
354 call mom_tracer_chkinv(debug_msg, g, gv, h_new, cs%tracer_reg)
355 endif
356 call cpu_clock_begin(id_clock_ale)
357
358 call ale_update_regrid_weights(cs%dt_offline, cs%ALE_CSp)
359 call pre_ale_adjustments(g, gv, us, h_new, cs%tv, cs%tracer_Reg, cs%ALE_CSp)
360 ! Uncomment this to adjust the target grids for diagnostics, if there have been thickness
361 ! adjustments, but the offline tracer code does not yet have the other corresponding calls
362 ! that would be needed to support remapping its output.
363 ! call diag_update_remap_grids(CS%diag, alt_h=h_new)
364
365 call ale_regrid(g, gv, us, h_new, h_post_remap, dzregrid, cs%tv, cs%ALE_CSp)
366
367 ! Remap all variables from the old grid h_new onto the new grid h_post_remap
368 call ale_remap_tracers(cs%ALE_CSp, g, gv, h_new, h_post_remap, cs%tracer_Reg, &
369 cs%debug, dt=cs%dt_offline)
370 if (allocated(cs%tv%SpV_avg)) cs%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
371
372 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
373 h_new(i,j,k) = h_post_remap(i,j,k)
374 enddo ; enddo ; enddo
375 call cpu_clock_end(id_clock_ale)
376
377 if (cs%debug) then
378 call hchksum(h_new, "h_new after ALE", g%HI, unscale=gv%H_to_MKS)
379 write(debug_msg, '(A,I4.4)') 'After ALE ', iter
380 call mom_tracer_chkinv(debug_msg, g, gv, h_new, cs%tracer_reg)
381 endif
382 endif
383
384 do k=1,nz ; do j=js,je ; do i=is,ie
385 uhtr_sub(i,j,k) = uhtr(i,j,k)
386 vhtr_sub(i,j,k) = vhtr(i,j,k)
387 enddo ; enddo ; enddo
388 call pass_var(h_new, g%Domain)
389 call pass_vector(uhtr_sub, vhtr_sub, g%Domain)
390
391 ! Check for whether we've used up all the advection, or if we need to move on because
392 ! advection has stalled
393 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
394 if (cs%print_adv_offline) then
395 write(mesg,'(A,ES24.16)') "Main advection remaining transport: ", tot_residual*hl2_to_kg_scale
396 call mom_mesg(mesg)
397 endif
398 ! If all the mass transports have been used u, then quit
399 if (tot_residual == 0.0) then
400 write(mesg,*) "Converged after iteration ", iter
401 call mom_mesg(mesg)
402 converged = .true.
403 exit
404 endif
405 ! If advection has stalled or the remaining residual is less than a specified amount, quit
406 if ( (tot_residual == prev_tot_residual) .or. (tot_residual<cs%min_residual) ) then
407 converged = .false.
408 exit
409 endif
410
411 prev_tot_residual = tot_residual
412
413 enddo
414
415 ! Make sure that uhtr and vhtr halos are updated
416 h_pre(:,:,:) = h_new(:,:,:)
417 call pass_vector(uhtr, vhtr, g%Domain)
418
419 if (cs%debug) then
420 call hchksum(h_pre, "h after offline_advection_ale", g%HI, unscale=gv%H_to_MKS)
421 call uvchksum("[uv]htr after offline_advection_ale", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
422 call mom_tracer_chkinv("After offline_advection_ale", g, gv, h_pre, cs%tracer_reg)
423 endif
424
425 call cpu_clock_end(cs%id_clock_offline_adv)
426
427end subroutine offline_advection_ale
428
429!> In the case where the main advection routine did not converge, something needs to be done with the remaining
430!! transport. Two different ways are offered, 'barotropic' means that the residual is distributed equally
431!! throughout the water column. 'upwards' attempts to redistribute the transport in the layers above and will
432!! eventually work down the entire water column
433subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, converged)
434 type(offline_transport_cs), pointer :: cs !< control structure from initialize_MOM
435 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
436 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
437 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
438 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
439 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
440 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
441 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
442 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
443 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
444 logical, intent(in ) :: converged !< True if the iterations have converged
445
446 logical :: x_before_y
447 ! Variables used to keep track of layer thicknesses at various points in the code
448 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
449 h_new, & ! New layer thicknesses [H ~> m or kg m-2]
450 h_vol ! Cell volume [H L2 ~> m3 or kg]
451
452 ! Used to calculate the eta diagnostics
453 real, dimension(SZI_(G),SZJ_(G)) :: eta_work ! The total column thickness [H ~> m or kg m-2]
454 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhr !< Remaining zonal mass transport [H L2 ~> m3 or kg]
455 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhr !< Remaining meridional mass transport [H L2 ~> m3 or kg]
456
457 character(len=256) :: mesg ! The text of an error message
458 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, iter
459 real :: hl2_to_kg_scale ! Unit conversion factors to cell mass [kg H-1 L-2 ~> kg m-3 or 1]
460 real :: prev_tot_residual, tot_residual ! The absolute value of the remaining transports [H L2 ~> m3 or kg]
461
462 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
463 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
464
465 x_before_y = cs%x_before_y
466 hl2_to_kg_scale = us%L_to_m**2*gv%H_to_kg_m2
467
468 if (cs%id_eta_pre_distribute>0) then
469 eta_work(:,:) = 0.0
470 do k=1,nz ; do j=js,je ; do i=is,ie
471 if (h_pre(i,j,k) > gv%Angstrom_H) then
472 eta_work(i,j) = eta_work(i,j) + h_pre(i,j,k)
473 endif
474 enddo ; enddo ; enddo
475 call post_data(cs%id_eta_pre_distribute, eta_work, cs%diag)
476 endif
477
478 ! These are used to find out how much will be redistributed in this routine
479 if (cs%id_h_redist>0) call post_data(cs%id_h_redist, h_pre, cs%diag)
480 if (cs%id_uhr_redist>0) call post_data(cs%id_uhr_redist, uhtr, cs%diag)
481 if (cs%id_vhr_redist>0) call post_data(cs%id_vhr_redist, vhtr, cs%diag)
482
483 if (converged) return
484
485 if (cs%debug) then
486 call mom_tracer_chkinv("Before redistribute ", g, gv, h_pre, cs%tracer_reg)
487 endif
488
489 call cpu_clock_begin(cs%id_clock_redistribute)
490
491 if (cs%redistribute_upwards .or. cs%redistribute_barotropic) then
492 do iter = 1, cs%num_off_iter
493
494 ! Perform upwards redistribution
495 if (cs%redistribute_upwards) then
496
497 ! Calculate the layer volumes at beginning of redistribute
498 do k=1,nz ; do j=js,je ; do i=is,ie
499 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
500 enddo ; enddo ; enddo
501 call pass_var(h_vol,g%Domain)
502 call pass_vector(uhtr, vhtr, g%Domain)
503
504 if (cs%debug) then
505 call mom_tracer_chksum("Before upwards redistribute ", cs%tracer_Reg, g)
506 call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
507 endif
508
509 if (x_before_y) then
510 call distribute_residual_uh_upwards(g, gv, h_vol, uhtr)
511 call distribute_residual_vh_upwards(g, gv, h_vol, vhtr)
512 else
513 call distribute_residual_vh_upwards(g, gv, h_vol, vhtr)
514 call distribute_residual_uh_upwards(g, gv, h_vol, uhtr)
515 endif
516
517 call advect_tracer(h_pre, uhtr, vhtr, cs%OBC, cs%dt_offline, g, gv, us, &
518 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
519 max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr)
520
521 if (cs%debug) then
522 call mom_tracer_chksum("After upwards redistribute ", cs%tracer_Reg, g)
523 endif
524
525 ! Convert h_new back to layer thickness for ALE remapping
526 do k=1,nz ; do j=js,je ; do i=is,ie
527 uhtr(i,j,k) = uhr(i,j,k)
528 vhtr(i,j,k) = vhr(i,j,k)
529 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
530 h_pre(i,j,k) = h_new(i,j,k)
531 enddo ; enddo ; enddo
532
533 endif ! redistribute upwards
534
535 ! Perform barotropic redistribution
536 if (cs%redistribute_barotropic) then
537
538 ! Calculate the layer volumes at beginning of redistribute
539 do k=1,nz ; do j=js,je ; do i=is,ie
540 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
541 enddo ; enddo ; enddo
542 call pass_var(h_vol, g%Domain)
543 call pass_vector(uhtr, vhtr, g%Domain)
544
545 if (cs%debug) then
546 call mom_tracer_chksum("Before barotropic redistribute ", cs%tracer_Reg, g)
547 call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
548 endif
549
550 if (x_before_y) then
551 call distribute_residual_uh_barotropic(g, gv, h_vol, uhtr)
552 call distribute_residual_vh_barotropic(g, gv, h_vol, vhtr)
553 else
554 call distribute_residual_vh_barotropic(g, gv, h_vol, vhtr)
555 call distribute_residual_uh_barotropic(g, gv, h_vol, uhtr)
556 endif
557
558 call advect_tracer(h_pre, uhtr, vhtr, cs%OBC, cs%dt_offline, g, gv, us, &
559 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
560 max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr)
561
562 if (cs%debug) then
563 call mom_tracer_chksum("After barotropic redistribute ", cs%tracer_Reg, g)
564 endif
565
566 ! Convert h_new back to layer thickness for ALE remapping
567 do k=1,nz ; do j=js,je ; do i=is,ie
568 uhtr(i,j,k) = uhr(i,j,k)
569 vhtr(i,j,k) = vhr(i,j,k)
570 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
571 h_pre(i,j,k) = h_new(i,j,k)
572 enddo ; enddo ; enddo
573
574 endif ! redistribute barotropic
575
576 ! Check to see if all transport has been exhausted
577 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
578 if (cs%print_adv_offline) then
579 write(mesg,'(A,ES24.16)') "Residual advection remaining transport: ", tot_residual*hl2_to_kg_scale
580 call mom_mesg(mesg)
581 endif
582 ! If the remaining residual is 0, then this return is done
583 if (tot_residual==0.0 ) then
584 exit
585 endif
586
587 if ( (tot_residual == prev_tot_residual) .or. (tot_residual<cs%min_residual) ) exit
588 prev_tot_residual = tot_residual
589
590 enddo ! Redistribution iteration
591 endif ! If one of the redistribution routines is requested
592
593 if (cs%id_eta_post_distribute>0) then
594 eta_work(:,:) = 0.0
595 do k=1,nz ; do j=js,je ; do i=is,ie
596 if (h_pre(i,j,k)>gv%Angstrom_H) then
597 eta_work(i,j) = eta_work(i,j) + h_pre(i,j,k)
598 endif
599 enddo ; enddo ; enddo
600 call post_data(cs%id_eta_post_distribute, eta_work, cs%diag)
601 endif
602
603 if (cs%id_uhr>0) call post_data(cs%id_uhr, uhtr, cs%diag)
604 if (cs%id_vhr>0) call post_data(cs%id_vhr, vhtr, cs%diag)
605
606 if (cs%debug) then
607 call hchksum(h_pre, "h_pre after redistribute", g%HI, unscale=gv%H_to_MKS)
608 call uvchksum("uhtr after redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
609 call mom_tracer_chkinv("after redistribute ", g, gv, h_new, cs%tracer_Reg)
610 endif
611
612 call cpu_clock_end(cs%id_clock_redistribute)
613
615
616!> Returns the sums of any non-negligible remaining transport [H L2 ~> m3 or kg] to check for advection convergence
617real function remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new)
618 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
619 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
620 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
621 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
622 intent(in ) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
623 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
624 intent(in ) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
625 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
626 intent(in ) :: h_new !< Layer thicknesses [H ~> m or kg m-2]
627
628 ! Local variables
629 real, dimension(SZI_(G),SZJ_(G)) :: trans_rem_col !< The vertical sum of the absolute value of
630 !! transports through the faces of a column [R Z L2 ~> kg].
631 real :: trans_cell !< The sum of the absolute value of the remaining transports through the faces
632 !! of a tracer cell [H L2 ~> m3 or kg]
633 integer :: i, j, k, is, ie, js, je, nz
634
635 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
636
637 trans_rem_col(:,:) = 0.0
638 do k=1,nz ; do j=js,je ; do i=is,ie
639 trans_cell = (abs(uhtr(i-1,j,k)) + abs(uhtr(i,j,k))) + &
640 (abs(vhtr(i,j-1,k)) + abs(vhtr(i,j,k)))
641 if (trans_cell > max(1.0e-16*h_new(i,j,k), gv%H_subroundoff) * g%areaT(i,j)) &
642 trans_rem_col(i,j) = trans_rem_col(i,j) + gv%H_to_RZ * trans_cell
643 enddo ; enddo ; enddo
644
645 ! The factor of 0.5 here is to avoid double-counting because two cells share a face.
646 remaining_transport_sum = 0.5 * gv%RZ_to_H * reproducing_sum(trans_rem_col, &
647 is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd), unscale=us%RZL2_to_kg)
648
649end function remaining_transport_sum
650
651!> The vertical/diabatic driver for offline tracers. First the eatr/ebtr associated with the interpolated
652!! vertical diffusivities are calculated and then any tracer column functions are done which can include
653!! vertical diffuvities and source/sink terms.
654subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, G, GV, US, CS, h_pre, tv, eatr, ebtr)
655
656 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
657 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
658 type(time_type), intent(in) :: time_end !< ending time of a segment, as a time type
659 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
660 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
661 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
662 type(offline_transport_cs), pointer :: cs !< control structure from initialize_MOM
663 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
664 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
665 type(thermo_var_ptrs), intent(in ) :: tv !< A structure pointing to various thermodynamic variables
666 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
667 intent(inout) :: eatr !< Entrainment from layer above [H ~> m or kg m-2]
668 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
669 intent(inout) :: ebtr !< Entrainment from layer below [H ~> m or kg m-2]
670
671 ! Local variables
672 real, dimension(SZI_(G),SZJ_(G)) :: &
673 sw, sw_vis, sw_nir !< Save old values of shortwave radiation [Q R Z T-1 ~> W m-2]
674 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Vertical distance across layers [Z ~> m]
675 real :: i_dzval ! An inverse distance between layer centers [Z-1 ~> m-1]
676 integer :: i, j, k, is, ie, js, je, nz
677 integer :: k_nonzero
678 real :: kd_bot ! Near-bottom diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
679 nz = gv%ke
680 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
681
682 call cpu_clock_begin(cs%id_clock_offline_diabatic)
683
684 call mom_mesg("Applying tracer source, sinks, and vertical mixing")
685
686 if (cs%debug) then
687 call hchksum(h_pre, "h_pre before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
688 call hchksum(eatr, "eatr before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
689 call hchksum(ebtr, "ebtr before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
690 call mom_tracer_chkinv("Before offline_diabatic_ale", g, gv, h_pre, cs%tracer_reg)
691 endif
692
693 call thickness_to_dz(h_pre, tv, dz, g, gv, us)
694
695 eatr(:,:,:) = 0.
696 ebtr(:,:,:) = 0.
697 ! Calculate eatr and ebtr if vertical diffusivity is read
698 ! Because the saved remapped diagnostics from the online run assume a zero minimum thickness
699 ! but ALE may have a minimum thickness. Flood the diffusivities for all layers with the value
700 ! of Kd closest to the bottom which is non-zero
701 do j=js,je ; do i=is,ie
702 k_nonzero = nz+1
703 ! Find the nonzero bottom Kd
704 do k=nz+1,1,-1
705 if (cs%Kd(i,j,k)>0.) then
706 kd_bot = cs%Kd(i,j,k)
707 k_nonzero = k
708 exit
709 endif
710 enddo
711 ! Flood the bottom interfaces
712 do k=k_nonzero,nz+1
713 cs%Kd(i,j,k) = kd_bot
714 enddo
715 enddo ; enddo
716
717 do j=js,je ; do i=is,ie
718 eatr(i,j,1) = 0.
719 enddo ; enddo
720 do k=2,nz ; do j=js,je ; do i=is,ie
721 i_dzval = 1.0 / (gv%dZ_subroundoff + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
722 eatr(i,j,k) = cs%dt_offline_vertical * i_dzval * cs%Kd(i,j,k)
723 ebtr(i,j,k-1) = eatr(i,j,k)
724 enddo ; enddo ; enddo
725 do j=js,je ; do i=is,ie
726 ebtr(i,j,nz) = 0.
727 enddo ; enddo
728
729 ! Add diurnal cycle for shortwave radiation (only used if run in ocean-only mode)
730 if (cs%diurnal_SW .and. cs%read_sw) then
731 sw(:,:) = fluxes%sw(:,:)
732 sw_vis(:,:) = fluxes%sw_vis_dir(:,:)
733 sw_nir(:,:) = fluxes%sw_nir_dir(:,:)
734 call offline_add_diurnal_sw(fluxes, g, time_start, time_end)
735 endif
736
737 if (associated(cs%optics)) &
738 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, &
739 cs%opacity_CSp, cs%tracer_flow_CSp)
740
741 ! Note that tracerBoundaryFluxesInOut within this subroutine should NOT be called
742 ! as the freshwater fluxes have already been accounted for
743 call call_tracer_column_fns(h_pre, h_pre, eatr, ebtr, fluxes, cs%MLD, cs%dt_offline_vertical, &
744 g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
745
746 if (cs%diurnal_SW .and. cs%read_sw) then
747 fluxes%sw(:,:) = sw(:,:)
748 fluxes%sw_vis_dir(:,:) = sw_vis(:,:)
749 fluxes%sw_nir_dir(:,:) = sw_nir(:,:)
750 endif
751
752 if (cs%debug) then
753 call hchksum(h_pre, "h_pre after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
754 call hchksum(eatr, "eatr after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
755 call hchksum(ebtr, "ebtr after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
756 call mom_tracer_chkinv("After offline_diabatic_ale", g, gv, h_pre, cs%tracer_reg)
757 endif
758
759 call cpu_clock_end(cs%id_clock_offline_diabatic)
760
761end subroutine offline_diabatic_ale
762
763!> Apply positive freshwater fluxes (into the ocean) and update netMassOut with only the negative
764!! (out of the ocean) fluxes
765subroutine offline_fw_fluxes_into_ocean(G, GV, CS, fluxes, h, in_flux_optional)
766 type(offline_transport_cs), intent(inout) :: cs !< Offline control structure
767 type(ocean_grid_type), intent(in) :: g !< Grid structure
768 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
769 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
770 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
771 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
772 real, dimension(SZI_(G),SZJ_(G)), &
773 optional, intent(in) :: in_flux_optional !< The total time-integrated amount
774 !! of tracer that leaves with freshwater
775 !! [CU H ~> Conc m or Conc kg m-2]
776
777 integer :: i, j, m
778 real, dimension(SZI_(G),SZJ_(G)) :: negative_fw !< store all negative fluxes [H ~> m or kg m-2]
779 logical :: update_h !< Flag for whether h should be updated
780
781 if ( present(in_flux_optional) ) &
782 call mom_error(warning, "Positive freshwater fluxes with non-zero tracer concentration not supported yet")
783
784 ! Set all fluxes to 0
785 negative_fw(:,:) = 0.
786
787 ! Sort fluxes into positive and negative
788 do j=g%jsc,g%jec ; do i=g%isc,g%iec
789 if (fluxes%netMassOut(i,j)<0.0) then
790 negative_fw(i,j) = fluxes%netMassOut(i,j)
791 fluxes%netMassOut(i,j) = 0.
792 endif
793 enddo ; enddo
794
795 if (cs%debug) then
796 call hchksum(h, "h before fluxes into ocean", g%HI, unscale=gv%H_to_MKS)
797 call mom_tracer_chkinv("Before fluxes into ocean", g, gv, h, cs%tracer_reg)
798 endif
799 do m = 1,cs%tracer_reg%ntr
800 ! Layer thicknesses should only be updated after the last tracer is finished
801 update_h = ( m == cs%tracer_reg%ntr )
802 call applytracerboundaryfluxesinout(g, gv, cs%tracer_reg%tr(m)%t, cs%dt_offline, fluxes, h, &
803 cs%evap_CFL_limit, cs%minimum_forcing_depth, update_h_opt=update_h)
804 enddo
805 if (cs%debug) then
806 call hchksum(h, "h after fluxes into ocean", g%HI, unscale=gv%H_to_MKS)
807 call mom_tracer_chkinv("After fluxes into ocean", g, gv, h, cs%tracer_reg)
808 endif
809
810 ! Now that fluxes into the ocean are done, save the negative fluxes for later
811 fluxes%netMassOut(:,:) = negative_fw(:,:)
812
813end subroutine offline_fw_fluxes_into_ocean
814
815!> Apply negative freshwater fluxes (out of the ocean)
816subroutine offline_fw_fluxes_out_ocean(G, GV, CS, fluxes, h, out_flux_optional)
817 type(offline_transport_cs), intent(inout) :: cs !< Offline control structure
818 type(ocean_grid_type), intent(in) :: g !< Grid structure
819 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
820 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
821 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
822 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
823 real, dimension(SZI_(G),SZJ_(G)), &
824 optional, intent(in) :: out_flux_optional !< The total time-integrated amount
825 !! of tracer that leaves with freshwater
826 !! [CU H ~> Conc m or Conc kg m-2]
827
828 integer :: m
829 logical :: update_h !< Flag for whether h should be updated
830
831 if ( present(out_flux_optional) ) &
832 call mom_error(warning, "Negative freshwater fluxes with non-zero tracer concentration not supported yet")
833
834 if (cs%debug) then
835 call hchksum(h, "h before fluxes out of ocean", g%HI, unscale=gv%H_to_MKS)
836 call mom_tracer_chkinv("Before fluxes out of ocean", g, gv, h, cs%tracer_reg)
837 endif
838 do m = 1, cs%tracer_reg%ntr
839 ! Layer thicknesses should only be updated after the last tracer is finished
840 update_h = ( m == cs%tracer_reg%ntr )
841 call applytracerboundaryfluxesinout(g, gv, cs%tracer_reg%tr(m)%t, cs%dt_offline, fluxes, h, &
842 cs%evap_CFL_limit, cs%minimum_forcing_depth, update_h_opt = update_h)
843 enddo
844 if (cs%debug) then
845 call hchksum(h, "h after fluxes out of ocean", g%HI, unscale=gv%H_to_MKS)
846 call mom_tracer_chkinv("Before fluxes out of ocean", g, gv, h, cs%tracer_reg)
847 endif
848
849end subroutine offline_fw_fluxes_out_ocean
850
851!> When in layer mode, 3D horizontal advection using stored mass fluxes must be used. Horizontal advection is
852!! done via tracer_advect, whereas the vertical component is actually handled by vertdiff in tracer_column_fns
853subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, CS, h_pre, eatr, ebtr, uhtr, vhtr)
854 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
855 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
856 real, intent(in) :: time_interval !< Offline transport time interval [T ~> s]
857 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
858 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
859 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
860 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
861 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
862 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
863 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
864 intent(inout) :: eatr !< Entrainment from layer above [H ~> m or kg m-2]
865 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
866 intent(inout) :: ebtr !< Entrainment from layer below [H ~> m or kg m-2]
867 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
868 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
869 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
870 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
871
872 ! Local variables
873
874 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhtr_sub ! Remaining zonal mass transports [H L2 ~> m3 or kg]
875 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub ! Remaining meridional mass transports [H L2 ~> m3 or kg]
876
877 real, dimension(SZI_(G),SZJB_(G)) :: rem_col_flux ! The summed absolute value of the remaining
878 ! mass fluxes through the faces of a column or within a column [R Z L2 ~> kg]
879 real :: sum_flux ! Globally summed absolute value of fluxes [R Z L2 ~> kg], which is
880 ! used to keep track of how close to convergence we are.
881
882 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
883 eatr_sub, & ! Layer entrainment rate from above for this sub-cycle [H ~> m or kg m-2]
884 ebtr_sub ! Layer entrainment rate from below for this sub-cycle [H ~> m or kg m-2]
885 ! Variables used to keep track of layer thicknesses at various points in the code
886 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
887 h_new, & ! Updated thicknesses [H ~> m or kg m-2]
888 h_vol ! Cell volumes [H L2 ~> m3 or kg]
889 ! Work arrays for temperature and salinity
890 integer :: iter
891 real :: dt_iter ! The timestep of each iteration [T ~> s]
892 character(len=160) :: mesg ! The text of an error message
893 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
894 integer :: isdb, iedb, jsdb, jedb
895 logical :: z_first, x_before_y
896
897 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
898 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
899 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
900
901 dt_iter = time_interval / real(max(1, cs%num_off_iter))
902 x_before_y = cs%x_before_y
903
904 do iter=1,cs%num_off_iter
905
906 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
907 eatr_sub(i,j,k) = eatr(i,j,k)
908 ebtr_sub(i,j,k) = ebtr(i,j,k)
909 enddo ; enddo ; enddo
910
911 do k=1,nz ; do j=js-1,je+1 ; do i=is-2,ie+1
912 uhtr_sub(i,j,k) = uhtr(i,j,k)
913 enddo ; enddo ; enddo
914
915 do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
916 vhtr_sub(i,j,k) = vhtr(i,j,k)
917 enddo ; enddo ; enddo
918
919 ! Calculate 3d mass transports to be used in this iteration
920 call limit_mass_flux_3d(g, gv, uhtr_sub, vhtr_sub, eatr_sub, ebtr_sub, h_pre)
921
922 if (z_first) then
923 ! First do vertical advection
924 call update_h_vertical_flux(g, gv, eatr_sub, ebtr_sub, h_pre, h_new)
925 call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, &
926 fluxes, cs%mld, dt_iter, g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
927 ! We are now done with the vertical mass transports, so now h_new is h_sub
928 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
929 h_pre(i,j,k) = h_new(i,j,k)
930 enddo ; enddo ; enddo
931 call pass_var(h_pre,g%Domain)
932
933 ! Second zonal and meridional advection
934 call update_h_horizontal_flux(g, gv, uhtr_sub, vhtr_sub, h_pre, h_new)
935 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
936 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
937 enddo ; enddo ; enddo
938 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, dt_iter, g, gv, us, &
939 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30)
940
941 ! Done with horizontal so now h_pre should be h_new
942 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
943 h_pre(i,j,k) = h_new(i,j,k)
944 enddo ; enddo ; enddo
945
946 endif
947
948 if (.not. z_first) then
949
950 ! First zonal and meridional advection
951 call update_h_horizontal_flux(g, gv, uhtr_sub, vhtr_sub, h_pre, h_new)
952 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
953 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
954 enddo ; enddo ; enddo
955 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, dt_iter, g, gv, us, cs%tracer_adv_CSp, &
956 cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30)
957
958 ! Done with horizontal so now h_pre should be h_new
959 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
960 h_pre(i,j,k) = h_new(i,j,k)
961 enddo ; enddo ; enddo
962
963 ! Second vertical advection
964 call update_h_vertical_flux(g, gv, eatr_sub, ebtr_sub, h_pre, h_new)
965 call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, &
966 fluxes, cs%mld, dt_iter, g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
967 ! We are now done with the vertical mass transports, so now h_new is h_sub
968 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
969 h_pre(i,j,k) = h_new(i,j,k)
970 enddo ; enddo ; enddo
971
972 endif
973
974 ! Update remaining transports
975 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
976 eatr(i,j,k) = eatr(i,j,k) - eatr_sub(i,j,k)
977 ebtr(i,j,k) = ebtr(i,j,k) - ebtr_sub(i,j,k)
978 enddo ; enddo ; enddo
979
980 do k=1,nz ; do j=js-1,je+1 ; do i=is-2,ie+1
981 uhtr(i,j,k) = uhtr(i,j,k) - uhtr_sub(i,j,k)
982 enddo ; enddo ; enddo
983
984 do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
985 vhtr(i,j,k) = vhtr(i,j,k) - vhtr_sub(i,j,k)
986 enddo ; enddo ; enddo
987
988 call pass_var(eatr,g%Domain)
989 call pass_var(ebtr,g%Domain)
990 call pass_var(h_pre,g%Domain)
991 call pass_vector(uhtr,vhtr,g%Domain)
992
993 ! Calculate how close we are to converging by summing the remaining fluxes at each point
994 rem_col_flux(:,:) = 0.0
995 do k=1,nz ; do j=js,je ; do i=is,ie
996 rem_col_flux(i,j) = rem_col_flux(i,j) + gv%H_to_RZ * &
997 ( (abs(eatr(i,j,k)) + abs(ebtr(i,j,k))) + &
998 ((abs(uhtr(i-1,j,k)) + abs(uhtr(i,j,k))) + &
999 (abs(vhtr(i,j-1,k)) + abs(vhtr(i,j,k))) ) )
1000 enddo ; enddo ; enddo
1001 sum_flux = reproducing_sum(rem_col_flux, is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd), &
1002 unscale=us%RZL2_to_kg)
1003
1004 if (sum_flux==0) then
1005 write(mesg,*) 'offline_advection_layer: Converged after iteration', iter
1006 call mom_mesg(mesg)
1007 exit
1008 else
1009 write(mesg,*) "offline_advection_layer: Iteration ", iter, " remaining total fluxes: ", sum_flux*us%RZL2_to_kg
1010 call mom_mesg(mesg)
1011 endif
1012
1013 ! Switch order of Strang split every iteration
1014 z_first = .not. z_first
1015 x_before_y = .not. x_before_y
1016 enddo
1017
1018end subroutine offline_advection_layer
1019
1020!> Update fields used in this round of offline transport. First fields are updated from files or from arrays
1021!! read during initialization. Then if in an ALE-dependent coordinate, regrid/remap fields.
1022subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
1023 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
1024 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
1025 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1026 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1027 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1028 intent(inout) :: h !< The regridded layer thicknesses [H ~> m or kg m-2]
1029 type(forcing), intent(inout) :: fluxes !< Pointers to forcing fields
1030 logical, intent(in ) :: do_ale !< True if using ALE
1031 ! Local variables
1032 integer :: stencil
1033 integer :: i, j, k, is, ie, js, je, nz
1034 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_start ! Initial thicknesses [H ~> m or kg m-2]
1035 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1036
1037 call cpu_clock_begin(cs%id_clock_read_fields)
1038 call calltree_enter("update_offline_fields, MOM_offline_main.F90")
1039
1040 if (cs%debug) then
1041 call uvchksum("[uv]htr before update_offline_fields", cs%uhtr, cs%vhtr, g%HI, &
1042 unscale=us%L_to_m**2*gv%H_to_kg_m2)
1043 call hchksum(cs%h_end, "h_end before update_offline_fields", g%HI, unscale=gv%H_to_MKS)
1044 call hchksum(cs%tv%T, "Temp before update_offline_fields", g%HI, unscale=us%C_to_degC)
1045 call hchksum(cs%tv%S, "Salt before update_offline_fields", g%HI, unscale=us%S_to_ppt)
1046 endif
1047
1048 ! Store a copy of the layer thicknesses before ALE regrid/remap
1049 h_start(:,:,:) = h(:,:,:)
1050
1051 ! Most fields will be read in from files
1052 call update_offline_from_files( g, gv, us, cs%nk_input, cs%mean_file, cs%sum_file, cs%snap_file, &
1053 cs%surf_file, cs%h_end, cs%uhtr, cs%vhtr, cs%tv%T, cs%tv%S, &
1054 cs%mld, cs%Kd, fluxes, cs%ridx_sum, cs%ridx_snap, cs%read_mld, &
1055 cs%read_sw, .not.cs%read_all_ts_uvh, do_ale)
1056 ! If uh, vh, h_end, temp, salt were read in at the beginning, fields are copied from those arrays
1057 if (cs%read_all_ts_uvh) then
1058 call update_offline_from_arrays(g, gv, cs%nk_input, cs%ridx_sum, cs%mean_file, cs%sum_file, &
1059 cs%snap_file, cs%uhtr, cs%vhtr, cs%h_end, cs%uhtr_all, cs%vhtr_all, &
1060 cs%hend_all, cs%tv%T, cs%tv%S, cs%temp_all, cs%salt_all)
1061 endif
1062 if (cs%debug) then
1063 call uvchksum("[uv]h after update offline from files and arrays", cs%uhtr, cs%vhtr, g%HI, &
1064 unscale=us%L_to_m**2*gv%H_to_kg_m2)
1065 call hchksum(cs%tv%T, "Temp after update offline from files and arrays", g%HI, unscale=us%C_to_degC)
1066 call hchksum(cs%tv%S, "Salt after update offline from files and arrays", g%HI, unscale=us%S_to_ppt)
1067 endif
1068
1069 ! If using an ALE-dependent vertical coordinate, fields will need to be remapped
1070 if (do_ale) then
1071 ! These halo passes are necessary because u, v fields will need information 1 step into the halo
1072 call pass_var(h, g%Domain)
1073 call pass_var(cs%tv%T, g%Domain)
1074 call pass_var(cs%tv%S, g%Domain)
1075 call ale_offline_inputs(cs%ALE_CSp, g, gv, us, h, cs%tv, cs%tracer_Reg, cs%uhtr, cs%vhtr, cs%Kd, &
1076 cs%debug, cs%OBC)
1077 if (cs%id_temp_regrid>0) call post_data(cs%id_temp_regrid, cs%tv%T, cs%diag)
1078 if (cs%id_salt_regrid>0) call post_data(cs%id_salt_regrid, cs%tv%S, cs%diag)
1079 if (cs%id_uhtr_regrid>0) call post_data(cs%id_uhtr_regrid, cs%uhtr, cs%diag)
1080 if (cs%id_vhtr_regrid>0) call post_data(cs%id_vhtr_regrid, cs%vhtr, cs%diag)
1081 if (cs%id_h_regrid>0) call post_data(cs%id_h_regrid, h, cs%diag)
1082 if (cs%debug) then
1083 call uvchksum("[uv]htr after ALE regridding/remapping of inputs", cs%uhtr, cs%vhtr, g%HI, &
1084 unscale=us%L_to_m**2*gv%H_to_kg_m2)
1085 call hchksum(h_start,"h_start after ALE regridding/remapping of inputs", g%HI, unscale=gv%H_to_MKS)
1086 endif
1087 endif
1088
1089 ! Update halos for some
1090 call pass_var(cs%h_end, g%Domain)
1091 call pass_var(cs%tv%T, g%Domain)
1092 call pass_var(cs%tv%S, g%Domain)
1093
1094 ! Update derived thermodynamic quantities.
1095 if (allocated(cs%tv%SpV_avg)) then
1096 stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1097 call calc_derived_thermo(cs%tv, cs%h_end, g, gv, us, halo=stencil)
1098 endif
1099
1100 ! Update the read indices
1101 cs%ridx_snap = next_modulo_time(cs%ridx_snap,cs%numtime)
1102 cs%ridx_sum = next_modulo_time(cs%ridx_sum,cs%numtime)
1103
1104 ! Apply masks/factors at T, U, and V points
1105 do k=1,nz ; do j=js,je ; do i=is,ie
1106 if (g%mask2dT(i,j)<1.0) then
1107 cs%h_end(i,j,k) = gv%Angstrom_H
1108 endif
1109 enddo ; enddo ; enddo
1110
1111 do k=1,nz+1 ; do j=js,je ; do i=is,ie
1112 cs%Kd(i,j,k) = max(0.0, cs%Kd(i,j,k))
1113 if (cs%Kd_max>0.) then
1114 cs%Kd(i,j,k) = min(cs%Kd_max, cs%Kd(i,j,k))
1115 endif
1116 enddo ; enddo ; enddo
1117
1118 do k=1,nz ; do j=js-1,je ; do i=is,ie
1119 if (g%mask2dCv(i,j)<1.0) then
1120 cs%vhtr(i,j,k) = 0.0
1121 endif
1122 enddo ; enddo ; enddo
1123
1124 do k=1,nz ; do j=js,je ; do i=is-1,ie
1125 if (g%mask2dCu(i,j)<1.0) then
1126 cs%uhtr(i,j,k) = 0.0
1127 endif
1128 enddo ; enddo ; enddo
1129
1130 if (cs%debug) then
1131 call uvchksum("[uv]htr after update_offline_fields", cs%uhtr, cs%vhtr, g%HI, &
1132 unscale=us%L_to_m**2*gv%H_to_kg_m2)
1133 call hchksum(cs%h_end, "h_end after update_offline_fields", g%HI, unscale=gv%H_to_MKS)
1134 call hchksum(cs%tv%T, "Temp after update_offline_fields", g%HI, unscale=us%C_to_degC)
1135 call hchksum(cs%tv%S, "Salt after update_offline_fields", g%HI, unscale=us%S_to_ppt)
1136 endif
1137
1138 call calltree_leave("update_offline_fields")
1139 call cpu_clock_end(cs%id_clock_read_fields)
1140
1141end subroutine update_offline_fields
1142
1143!> Initialize additional diagnostics required for offline tracer transport
1144subroutine register_diags_offline_transport(Time, diag, CS, GV, US)
1145
1146 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
1147 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1148 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1149 type(time_type), intent(in) :: time !< current model time
1150 type(diag_ctrl), intent(in) :: diag !< Structure that regulates diagnostic output
1151
1152 ! U-cell fields
1153 cs%id_uhr = register_diag_field('ocean_model', 'uhr', diag%axesCuL, time, &
1154 'Zonal thickness fluxes remaining at end of advection', &
1155 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1156 cs%id_uhr_redist = register_diag_field('ocean_model', 'uhr_redist', diag%axesCuL, time, &
1157 'Zonal thickness fluxes to be redistributed vertically', &
1158 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1159 cs%id_uhr_end = register_diag_field('ocean_model', 'uhr_end', diag%axesCuL, time, &
1160 'Zonal thickness fluxes at end of offline step', &
1161 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1162
1163 ! V-cell fields
1164 cs%id_vhr = register_diag_field('ocean_model', 'vhr', diag%axesCvL, time, &
1165 'Meridional thickness fluxes remaining at end of advection', &
1166 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1167 cs%id_vhr_redist = register_diag_field('ocean_model', 'vhr_redist', diag%axesCvL, time, &
1168 'Meridional thickness to be redistributed vertically', &
1169 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1170 cs%id_vhr_end = register_diag_field('ocean_model', 'vhr_end', diag%axesCvL, time, &
1171 'Meridional thickness at end of offline step', &
1172 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1173
1174 ! T-cell fields
1175 cs%id_hdiff = register_diag_field('ocean_model', 'hdiff', diag%axesTL, time, &
1176 'Difference between the stored and calculated layer thickness', &
1177 'm', conversion=gv%H_to_m)
1178 cs%id_hr = register_diag_field('ocean_model', 'hr', diag%axesTL, time, &
1179 'Layer thickness at end of offline step', 'm', conversion=gv%H_to_m)
1180 cs%id_ear = register_diag_field('ocean_model', 'ear', diag%axesTL, time, &
1181 'Remaining thickness entrained from above', 'm')
1182 cs%id_ebr = register_diag_field('ocean_model', 'ebr', diag%axesTL, time, &
1183 'Remaining thickness entrained from below', 'm')
1184 cs%id_eta_pre_distribute = register_diag_field('ocean_model','eta_pre_distribute', &
1185 diag%axesT1, time, 'Total water column height before residual transport redistribution', &
1186 'm', conversion=gv%H_to_m)
1187 cs%id_eta_post_distribute = register_diag_field('ocean_model','eta_post_distribute', &
1188 diag%axesT1, time, 'Total water column height after residual transport redistribution', &
1189 'm', conversion=gv%H_to_m)
1190 cs%id_eta_diff_end = register_diag_field('ocean_model','eta_diff_end', diag%axesT1, time, &
1191 'Difference in total water column height from online and offline ' // &
1192 'at the end of the offline timestep', 'm', conversion=gv%H_to_m)
1193 cs%id_h_redist = register_diag_field('ocean_model','h_redist', diag%axesTL, time, &
1194 'Layer thicknesses before redistribution of mass fluxes', &
1195 get_thickness_units(gv), conversion=gv%H_to_MKS)
1196
1197 ! Regridded/remapped input fields
1198 cs%id_uhtr_regrid = register_diag_field('ocean_model', 'uhtr_regrid', diag%axesCuL, time, &
1199 'Zonal mass transport regridded/remapped onto offline grid', &
1200 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1201 cs%id_vhtr_regrid = register_diag_field('ocean_model', 'vhtr_regrid', diag%axesCvL, time, &
1202 'Meridional mass transport regridded/remapped onto offline grid', &
1203 'kg', conversion=us%L_to_m**2*gv%H_to_kg_m2)
1204 cs%id_temp_regrid = register_diag_field('ocean_model', 'temp_regrid', diag%axesTL, time, &
1205 'Temperature regridded/remapped onto offline grid',&
1206 'C', conversion=us%C_to_degC)
1207 cs%id_salt_regrid = register_diag_field('ocean_model', 'salt_regrid', diag%axesTL, time, &
1208 'Salinity regridded/remapped onto offline grid', &
1209 'g kg-1', conversion=us%S_to_ppt)
1210 cs%id_h_regrid = register_diag_field('ocean_model', 'h_regrid', diag%axesTL, time, &
1211 'Layer thicknesses regridded/remapped onto offline grid', &
1212 'm', conversion=gv%H_to_m)
1213
1215
1216!> Posts diagnostics related to offline convergence diagnostics
1217subroutine post_offline_convergence_diags(G, GV, CS, h_off, h_end, uhtr, vhtr)
1218 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1219 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1220 type(offline_transport_cs), intent(in ) :: cs !< Offline control structure
1221 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1222 intent(inout) :: h_off !< Thicknesses at end of offline step [H ~> m or kg m-2]
1223 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1224 intent(inout) :: h_end !< Stored thicknesses [H ~> m or kg m-2]
1225 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1226 intent(inout) :: uhtr !< Remaining zonal mass transport [H L2 ~> m3 or kg]
1227 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1228 intent(inout) :: vhtr !< Remaining meridional mass transport [H L2 ~> m3 or kg]
1229
1230 real, dimension(SZI_(G),SZJ_(G)) :: eta_diff ! Differences in column thickness [H ~> m or kg m-2]
1231 integer :: i, j, k
1232
1233 if (cs%id_eta_diff_end>0) then
1234 ! Calculate difference in column thickness
1235 eta_diff = 0.
1236 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1237 eta_diff(i,j) = eta_diff(i,j) + h_off(i,j,k)
1238 enddo ; enddo ; enddo
1239 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1240 eta_diff(i,j) = eta_diff(i,j) - h_end(i,j,k)
1241 enddo ; enddo ; enddo
1242
1243 call post_data(cs%id_eta_diff_end, eta_diff, cs%diag)
1244 endif
1245
1246 if (cs%id_hdiff>0) call post_data(cs%id_hdiff, h_off-h_end, cs%diag)
1247 if (cs%id_hr>0) call post_data(cs%id_hr, h_off, cs%diag)
1248 if (cs%id_uhr_end>0) call post_data(cs%id_uhr_end, uhtr, cs%diag)
1249 if (cs%id_vhr_end>0) call post_data(cs%id_vhr_end, vhtr, cs%diag)
1250
1251end subroutine post_offline_convergence_diags
1252
1253!> Extracts members of the offline main control structure. All arguments are optional except
1254!! the control structure itself
1255subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, vertical_time, &
1256 dt_offline, dt_offline_vertical, skip_diffusion)
1257 type(offline_transport_cs), target, intent(in ) :: cs !< Offline control structure
1258 ! Returned optional arguments
1259 real, dimension(:,:,:), optional, pointer :: uhtr !< Remaining zonal mass transport [H L2 ~> m3 or kg]
1260 real, dimension(:,:,:), optional, pointer :: vhtr !< Remaining meridional mass transport [H L2 ~> m3 or kg]
1261 real, dimension(:,:,:), optional, pointer :: eatr !< Amount of fluid entrained from the layer above within
1262 !! one time step [H ~> m or kg m-2]
1263 real, dimension(:,:,:), optional, pointer :: ebtr !< Amount of fluid entrained from the layer below within
1264 !! one time step [H ~> m or kg m-2]
1265 real, dimension(:,:,:), optional, pointer :: h_end !< Thicknesses at the end of offline timestep
1266 !! [H ~> m or kg m-2]
1267 type(time_type), optional, pointer :: accumulated_time !< Length of time accumulated in the
1268 !! current offline interval
1269 type(time_type), optional, pointer :: vertical_time !< The next value of accumulate_time at which to
1270 !! vertical processes
1271 real, optional, intent( out) :: dt_offline !< Timestep used for offline tracers [T ~> s]
1272 real, optional, intent( out) :: dt_offline_vertical !< Timestep used for calls to tracer
1273 !! vertical physics [T ~> s]
1274 logical, optional, intent( out) :: skip_diffusion !< Skips horizontal diffusion of tracers
1275
1276 ! Pointers to 3d members
1277 if (present(uhtr)) uhtr => cs%uhtr
1278 if (present(vhtr)) vhtr => cs%vhtr
1279 if (present(eatr)) eatr => cs%eatr
1280 if (present(ebtr)) ebtr => cs%ebtr
1281 if (present(h_end)) h_end => cs%h_end
1282
1283 ! Pointers to integer members which need to be modified
1284 if (present(accumulated_time)) accumulated_time => cs%accumulated_time
1285 if (present(vertical_time)) vertical_time => cs%vertical_time
1286
1287 ! Return value of non-modified integers
1288 if (present(dt_offline)) dt_offline = cs%dt_offline
1289 if (present(dt_offline_vertical)) dt_offline_vertical = cs%dt_offline_vertical
1290 if (present(skip_diffusion)) skip_diffusion = cs%skip_diffusion
1291
1292end subroutine extract_offline_main
1293
1294!> Inserts (assigns values to) members of the offline main control structure. All arguments
1295!! are optional except for the CS itself
1296subroutine insert_offline_main(CS, ALE_CSp, diabatic_CSp, diag, OBC, tracer_adv_CSp, &
1297 tracer_flow_CSp, tracer_Reg, tv, x_before_y, debug)
1298 type(offline_transport_cs), intent(inout) :: cs !< Offline control structure
1299 ! Inserted optional arguments
1300 type(ale_cs), &
1301 target, optional, intent(in ) :: ale_csp !< A pointer to the ALE control structure
1302 type(diabatic_cs), &
1303 target, optional, intent(in ) :: diabatic_csp !< A pointer to the diabatic control structure
1304 type(diag_ctrl), &
1305 target, optional, intent(in ) :: diag !< A pointer to the structure that regulates diagnostic output
1306 type(ocean_obc_type), &
1307 target, optional, intent(in ) :: obc !< A pointer to the open boundary condition control structure
1308 type(tracer_advect_cs), &
1309 target, optional, intent(in ) :: tracer_adv_csp !< A pointer to the tracer advection control structure
1310 type(tracer_flow_control_cs), &
1311 target, optional, intent(in ) :: tracer_flow_csp !< A pointer to the tracer flow control control structure
1312 type(tracer_registry_type), &
1313 target, optional, intent(in ) :: tracer_reg !< A pointer to the tracer registry
1314 type(thermo_var_ptrs), &
1315 target, optional, intent(in ) :: tv !< A structure pointing to various thermodynamic variables
1316 logical, optional, intent(in ) :: x_before_y !< Indicates which horizontal direction is advected first
1317 logical, optional, intent(in ) :: debug !< If true, write verbose debugging messages
1318
1319
1320 if (present(ale_csp)) cs%ALE_CSp => ale_csp
1321 if (present(diabatic_csp)) cs%diabatic_CSp => diabatic_csp
1322 if (present(diag)) cs%diag => diag
1323 if (present(obc)) cs%OBC => obc
1324 if (present(tracer_adv_csp)) cs%tracer_adv_CSp => tracer_adv_csp
1325 if (present(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
1326 if (present(tracer_reg)) cs%tracer_Reg => tracer_reg
1327 if (present(tv)) cs%tv => tv
1328 if (present(x_before_y)) cs%x_before_y = x_before_y
1329 if (present(debug)) cs%debug = debug
1330
1331end subroutine insert_offline_main
1332
1333!> Initializes the control structure for offline transport and reads in some of the
1334! run time parameters from MOM_input
1335subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US)
1336
1337 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1338 type(offline_transport_cs), pointer :: cs !< Offline control structure
1339 type(diabatic_cs), intent(in) :: diabatic_csp !< The diabatic control structure
1340 type(ocean_grid_type), target, intent(in) :: g !< ocean grid structure
1341 type(verticalgrid_type), target, intent(in) :: gv !< ocean vertical grid structure
1342 type(unit_scale_type), target, intent(in) :: us !< A dimensional unit scaling type
1343
1344 character(len=40) :: mdl = "offline_transport"
1345 character(len=20) :: redistribute_method
1346 ! This include declares and sets the variable "version".
1347# include "version_variable.h"
1348 integer :: is, ie, js, je, isd, ied, jsd, jed, nz
1349 integer :: isdb, iedb, jsdb, jedb
1350
1351 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1352 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1353 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1354
1355 call calltree_enter("offline_transport_init, MOM_offline_control.F90")
1356
1357 if (associated(cs)) then
1358 call mom_error(warning, "offline_transport_init called with an associated control structure.")
1359 return
1360 endif
1361 allocate(cs)
1362 call log_version(param_file, mdl, version, "This module allows for tracers to be run offline")
1363
1364 ! Parse MOM_input for offline control
1365 call get_param(param_file, mdl, "OFFLINEDIR", cs%offlinedir, &
1366 "Input directory where the offline fields can be found", fail_if_missing=.true.)
1367 call get_param(param_file, mdl, "OFF_SUM_FILE", cs%sum_file, &
1368 "Filename where the accumulated fields can be found", fail_if_missing=.true.)
1369 call get_param(param_file, mdl, "OFF_SNAP_FILE", cs%snap_file, &
1370 "Filename where snapshot fields can be found", fail_if_missing=.true.)
1371 call get_param(param_file, mdl, "OFF_MEAN_FILE", cs%mean_file, &
1372 "Filename where averaged fields can be found", fail_if_missing=.true.)
1373 call get_param(param_file, mdl, "OFF_SURF_FILE", cs%surf_file, &
1374 "Filename where averaged fields can be found", fail_if_missing=.true.)
1375 call get_param(param_file, mdl, "NUMTIME", cs%numtime, &
1376 "Number of timelevels in offline input files", fail_if_missing=.true.)
1377 call get_param(param_file, mdl, "NK_INPUT", cs%nk_input, &
1378 "Number of vertical levels in offline input files", default=nz)
1379 call get_param(param_file, mdl, "DT_OFFLINE", cs%dt_offline, &
1380 "Length of time between reading in of input fields", units='s', scale=us%s_to_T, fail_if_missing=.true.)
1381 call get_param(param_file, mdl, "DT_OFFLINE_VERTICAL", cs%dt_offline_vertical, &
1382 "Length of the offline timestep for tracer column sources/sinks " //&
1383 "This should be set to the length of the coupling timestep for " //&
1384 "tracers which need shortwave fluxes", units="s", scale=us%s_to_T, fail_if_missing=.true.)
1385 call get_param(param_file, mdl, "START_INDEX", cs%start_index, &
1386 "Which time index to start from", default=1)
1387 call get_param(param_file, mdl, "FIELDS_ARE_OFFSET", cs%fields_are_offset, &
1388 "True if the time-averaged fields and snapshot fields "//&
1389 "are offset by one time level", default=.false.)
1390 call get_param(param_file, mdl, "REDISTRIBUTE_METHOD", redistribute_method, &
1391 "Redistributes any remaining horizontal fluxes throughout " //&
1392 "the rest of water column. Options are 'barotropic' which " //&
1393 "evenly distributes flux throughout the entire water column, " //&
1394 "'upwards' which adds the maximum of the remaining flux in " //&
1395 "each layer above, both which first applies upwards and then " //&
1396 "barotropic, and 'none' which does no redistribution", &
1397 default='barotropic')
1398 call get_param(param_file, mdl, "NUM_OFF_ITER", cs%num_off_iter, &
1399 "Number of iterations to subdivide the offline tracer advection and diffusion", &
1400 default=60)
1401 call get_param(param_file, mdl, "OFF_ALE_MOD", cs%off_ale_mod, &
1402 "Sets how many horizontal advection steps are taken before an ALE "//&
1403 "remapping step is done. 1 would be x->y->ALE, 2 would be x->y->x->y->ALE", default=1)
1404 call get_param(param_file, mdl, "PRINT_ADV_OFFLINE", cs%print_adv_offline, &
1405 "Print diagnostic output every advection subiteration", default=.false.)
1406 call get_param(param_file, mdl, "SKIP_DIFFUSION_OFFLINE", cs%skip_diffusion, &
1407 "Do not do horizontal diffusion", default=.false.)
1408 call get_param(param_file, mdl, "READ_SW", cs%read_sw, &
1409 "Read in shortwave radiation field instead of using values from the coupler "//&
1410 "when in offline tracer mode", default=.false.)
1411 call get_param(param_file, mdl, "READ_MLD", cs%read_mld, &
1412 "Read in mixed layer depths for tracers which exchange with the atmosphere "//&
1413 "when in offline tracer mode", default=.false.)
1414 call get_param(param_file, mdl, "MLD_VAR_NAME", cs%mld_var_name, &
1415 "Name of the variable containing the depth of active mixing", default='ePBL_h_ML')
1416 call get_param(param_file, mdl, "OFFLINE_ADD_DIURNAL_SW", cs%diurnal_sw, &
1417 "Adds a synthetic diurnal cycle in the same way that the ice "//&
1418 "model would have when time-averaged fields of shortwave "//&
1419 "radiation are read in", default=.false.)
1420 call get_param(param_file, mdl, "KD_MAX", cs%Kd_max, &
1421 "The maximum permitted increment for the diapycnal "//&
1422 "diffusivity from TKE-based parameterizations, or a "//&
1423 "negative value for no limit.", units="m2 s-1", default=-1.0, scale=gv%m2_s_to_HZ_T)
1424 call get_param(param_file, mdl, "MIN_RESIDUAL_TRANSPORT", cs%min_residual, &
1425 "How much remaining transport before the main offline advection is exited. "//&
1426 "The default value corresponds to about 1 meter of difference in a grid cell", &
1427 default=1.e9, units="m3", scale=gv%m_to_H*us%m_to_L**2)
1428 call get_param(param_file, mdl, "READ_ALL_TS_UVH", cs%read_all_ts_uvh, &
1429 "Reads all time levels of a subset of the fields necessary to run " // &
1430 "the model offline. This can require a large amount of memory "// &
1431 "and will make initialization very slow. However, for offline "// &
1432 "runs spanning more than a year this can reduce total I/O overhead", &
1433 default=.false.)
1434
1435 ! Concatenate offline directory and file names
1436 cs%snap_file = trim(cs%offlinedir)//trim(cs%snap_file)
1437 cs%mean_file = trim(cs%offlinedir)//trim(cs%mean_file)
1438 cs%sum_file = trim(cs%offlinedir)//trim(cs%sum_file)
1439 cs%surf_file = trim(cs%offlinedir)//trim(cs%surf_file)
1440
1441 cs%num_vert_iter = cs%dt_offline / cs%dt_offline_vertical
1442
1443 ! Map redistribute_method onto logicals in CS
1444 select case (redistribute_method)
1445 case ('barotropic')
1446 cs%redistribute_barotropic = .true.
1447 cs%redistribute_upwards = .false.
1448 case ('upwards')
1449 cs%redistribute_barotropic = .false.
1450 cs%redistribute_upwards = .true.
1451 case ('both')
1452 cs%redistribute_barotropic = .true.
1453 cs%redistribute_upwards = .true.
1454 case ('none')
1455 cs%redistribute_barotropic = .false.
1456 cs%redistribute_upwards = .false.
1457 end select
1458
1459 ! Set the accumulated time to zero
1460 cs%accumulated_time = real_to_time(0.0)
1461 cs%vertical_time = cs%accumulated_time
1462 ! Set the starting read index for time-averaged and snapshotted fields
1463 cs%ridx_sum = cs%start_index
1464 if (cs%fields_are_offset) cs%ridx_snap = next_modulo_time(cs%start_index,cs%numtime)
1465 if (.not. cs%fields_are_offset) cs%ridx_snap = cs%start_index
1466
1467 ! Copy members from other modules
1468 call extract_diabatic_member(diabatic_csp, opacity_csp=cs%opacity_CSp, optics_csp=cs%optics, &
1469 diabatic_aux_csp=cs%diabatic_aux_CSp, &
1470 evap_cfl_limit=cs%evap_CFL_limit, &
1471 minimum_forcing_depth=cs%minimum_forcing_depth)
1472
1473 ! Allocate arrays
1474 allocate(cs%uhtr(isdb:iedb,jsd:jed,nz), source=0.0)
1475 allocate(cs%vhtr(isd:ied,jsdb:jedb,nz), source=0.0)
1476 allocate(cs%eatr(isd:ied,jsd:jed,nz), source=0.0)
1477 allocate(cs%ebtr(isd:ied,jsd:jed,nz), source=0.0)
1478 allocate(cs%h_end(isd:ied,jsd:jed,nz), source=0.0)
1479 allocate(cs%Kd(isd:ied,jsd:jed,nz+1), source=0.0)
1480 if (cs%read_mld) allocate(cs%mld(g%isd:g%ied,g%jsd:g%jed), source=0.0)
1481
1482 if (cs%read_all_ts_uvh) then
1483 call read_all_input(cs, g, gv, us)
1484 endif
1485
1486 ! Initialize ids for clocks used in offline routines
1487 cs%id_clock_read_fields = cpu_clock_id('(Offline read fields)',grain=clock_module)
1488 cs%id_clock_offline_diabatic = cpu_clock_id('(Offline diabatic)',grain=clock_module)
1489 cs%id_clock_offline_adv = cpu_clock_id('(Offline transport)',grain=clock_module)
1490 cs%id_clock_redistribute = cpu_clock_id('(Offline redistribute)',grain=clock_module)
1491
1492 call calltree_leave("offline_transport_init")
1493
1494end subroutine offline_transport_init
1495
1496!> Coordinates the allocation and reading in all time levels of uh, vh, hend, temp, and salt from files. Used
1497!! when read_all_ts_uvh
1498subroutine read_all_input(CS, G, GV, US)
1499 type(offline_transport_cs), intent(inout) :: CS !< Control structure for offline module
1500 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1501 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1502 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1503
1504 integer :: isd, ied, jsd, jed, nz, t, ntime
1505 integer :: IsdB, IedB, JsdB, JedB
1506
1507 nz = gv%ke ; ntime = cs%numtime
1508 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1509 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1510
1511 ! Extra safety check that we're not going to overallocate any arrays
1512 if (cs%read_all_ts_uvh) then
1513 if (allocated(cs%uhtr_all)) call mom_error(fatal, "uhtr_all is already allocated")
1514 if (allocated(cs%vhtr_all)) call mom_error(fatal, "vhtr_all is already allocated")
1515 if (allocated(cs%hend_all)) call mom_error(fatal, "hend_all is already allocated")
1516 if (allocated(cs%temp_all)) call mom_error(fatal, "temp_all is already allocated")
1517 if (allocated(cs%salt_all)) call mom_error(fatal, "salt_all is already allocated")
1518
1519 allocate(cs%uhtr_all(isdb:iedb,jsd:jed,nz,ntime), source=0.0)
1520 allocate(cs%vhtr_all(isd:ied,jsdb:jedb,nz,ntime), source=0.0)
1521 allocate(cs%hend_all(isd:ied,jsd:jed,nz,ntime), source=0.0)
1522 allocate(cs%temp_all(isd:ied,jsd:jed,nz,1:ntime), source=0.0)
1523 allocate(cs%salt_all(isd:ied,jsd:jed,nz,1:ntime), source=0.0)
1524
1525 call mom_mesg("Reading in uhtr, vhtr, h_start, h_end, temp, salt")
1526 do t = 1,ntime
1527 call mom_read_vector(cs%snap_file, 'uhtr_sum', 'vhtr_sum', cs%uhtr_all(:,:,1:cs%nk_input,t), &
1528 cs%vhtr_all(:,:,1:cs%nk_input,t), g%Domain, timelevel=t, &
1529 scale=us%m_to_L**2*gv%kg_m2_to_H)
1530 call mom_read_data(cs%snap_file,'h_end', cs%hend_all(:,:,1:cs%nk_input,t), g%Domain, &
1531 timelevel=t, position=center, scale=gv%kg_m2_to_H)
1532 call mom_read_data(cs%mean_file,'temp', cs%temp_all(:,:,1:cs%nk_input,t), g%Domain, &
1533 timelevel=t, position=center, scale=us%degC_to_C)
1534 call mom_read_data(cs%mean_file,'salt', cs%salt_all(:,:,1:cs%nk_input,t), g%Domain, &
1535 timelevel=t, position=center, scale=us%ppt_to_S)
1536 enddo
1537 endif
1538
1539end subroutine read_all_input
1540
1541!> Deallocates (if necessary) arrays within the offline control structure
1542subroutine offline_transport_end(CS)
1543 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
1544
1545 ! Explicitly allocate all allocatable arrays
1546 deallocate(cs%uhtr)
1547 deallocate(cs%vhtr)
1548 deallocate(cs%eatr)
1549 deallocate(cs%ebtr)
1550 deallocate(cs%h_end)
1551 deallocate(cs%Kd)
1552 if (cs%read_mld) deallocate(cs%mld)
1553 if (cs%read_all_ts_uvh) then
1554 deallocate(cs%uhtr_all)
1555 deallocate(cs%vhtr_all)
1556 deallocate(cs%hend_all)
1557 deallocate(cs%temp_all)
1558 deallocate(cs%salt_all)
1559 endif
1560
1561 deallocate(cs)
1562
1563end subroutine offline_transport_end
1564
1565!> \namespace mom_offline_main
1566!! \section offline_overview Offline Tracer Transport in MOM6
1567!! 'Offline tracer modeling' uses physical fields (e.g. mass transports and layer thicknesses) saved
1568!! from a previous integration of the physical model to transport passive tracers. These fields are
1569!! accumulated or averaged over a period of time (in this test case, 1 day) and used to integrate
1570!! portions of the MOM6 code base that handle the 3d advection and diffusion of passive tracers.
1571!!
1572!! The distribution of tracers in the ocean modeled offline should not be expected to match an online
1573!! simulation. Accumulating transports over more than one online model timestep implicitly assumes
1574!! homogeneity over that time period and essentially aliases over processes that occur with higher
1575!! frequency. For example, consider the case of a surface boundary layer with a strong diurnal cycle.
1576!! An offline simulation with a 1 day timestep, captures the net transport into or out of that layer,
1577!! but not the exact cycling. This effective aliasing may also complicate online model configurations
1578!! which strongly-eddying regions. In this case, the offline model timestep must be limited to some
1579!! fraction of the eddy correlation timescale. Lastly, the nonlinear advection scheme which applies
1580!! limited mass-transports over a sequence of iterations means that tracers are not transported along
1581!! exactly the same path as they are in the online model.
1582!!
1583!! This capability has currently targeted the Baltic_ALE_z test case, though some work has also been
1584!! done with the OM4 1/2 degree configuration. Work is ongoing to develop recommendations and best
1585!! practices for investigators seeking to use MOM6 for offline tracer modeling.
1586!!
1587!! \section offline_technical Implementation of offline routine in MOM6
1588!!
1589!! The subroutine step_tracers that coordinates this can be found in MOM.F90 and is only called
1590!! using the solo ocean driver. This is to avoid issues with coupling to other climate components
1591!! that may be relying on fluxes from the ocean to be coupled more often than the offline time step.
1592!! Other routines related to offline tracer modeling can be found in tracers/MOM_offline_control.F90
1593!!
1594!! As can also be seen in the comments for the step_tracers subroutine, an offline time step
1595!! comprises the following steps:
1596!! -# Using the layer thicknesses and tracer concentrations from the previous timestep,
1597!! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to
1598!! tracer_column_fns.
1599!! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline
1600!! -# Half of the accumulated surface freshwater fluxes are applied
1601!! START ITERATION
1602!! -# Accumulated mass fluxes are used to do horizontal transport. The number of iterations
1603!! used in advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are
1604!! stored for later use and resulting layer thicknesses fed into the next step
1605!! -# Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for
1606!! layers which might 'vanish' because of horizontal mass transport to be 'reinflated'
1607!! and essentially allows for the vertical transport of tracers
1608!! -# Check that transport is done if the remaining mass fluxes equals 0 or if the max
1609!! number of iterations has been reached
1610!! END ITERATION
1611!! -# Repeat steps 1 and 2
1612!! -# Redistribute any residual mass fluxes that remain after the advection iterations
1613!! in a barotropic manner, progressively upward through the water column.
1614!! -# Force a remapping to the stored layer thicknesses that correspond to the snapshot of
1615!! the online model at the end of an accumulation interval
1616!! -# Reset T/S and h to their stored snapshotted values to prevent model drift
1617!!
1618!! \section offline_evaluation Evaluating the utility of an offline tracer model
1619!! How well an offline tracer model can be used as an alternative to integrating tracers online
1620!! with the prognostic model must be evaluated for each application. This efficacy may be related
1621!! to the native coordinate of the online model, to the length of the offline timestep, and to the
1622!! behavior of the tracer itself.
1623!!
1624!! A framework for formally regression testing the offline capability still needs to be developed.
1625!! However, as a simple way of testing whether the offline model is nominally behaving as expected,
1626!! the total inventory of the advection test tracers (tr1, tr2, etc.) should be conserved between
1627!! time steps except for the last 4 decimal places. As a general guideline, an offline timestep of
1628!! 5 days or less.
1629!!
1630!! \section offline_parameters Runtime parameters for offline tracers
1631!! - OFFLINEDIR: Input directory where the offline fields can be found
1632!! - OFF_SUM_FILE: Filename where the accumulated fields can be found (e.g. horizontal mass transports)
1633!! - OFF_SNAP_FILE: Filename where snapshot fields can be found (e.g. end of timestep layer thickness)
1634!! - START_INDEX: Which timelevel of the input files to read first
1635!! - NUMTIME: How many timelevels to read before 'looping' back to 1
1636!! - FIELDS_ARE_OFFSET: True if the time-averaged fields and snapshot fields are offset by one
1637!! time level, probably not needed
1638!! -NUM_OFF_ITER: Maximum number of iterations to do for the nonlinear advection scheme
1639!! -REDISTRIBUTE_METHOD: Redistributes any remaining horizontal fluxes throughout the rest of water column.
1640!! Options are 'barotropic' which "evenly distributes flux throughout the entire water
1641!! column,'upwards' which adds the maximum of the remaining flux in each layer above,
1642!! and 'none' which does no redistribution"
1643
1644end module mom_offline_main
1645