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 real :: hmix_fixed !< A fixed mixed layer depth to use when read_mld is false [Z ~> m]
106 logical :: diurnal_sw !< Adds a synthetic diurnal cycle on shortwave radiation
107 logical :: debug !< If true, write verbose debugging messages
108 logical :: redistribute_barotropic !< Redistributes column-summed residual transports throughout
109 !! a column weighted by thickness
110 logical :: redistribute_upwards !< Redistributes remaining fluxes only in layers above
111 !! the current one based as the max allowable transport
112 !! in that cell
113 logical :: read_all_ts_uvh !< If true, then all timelevels of temperature, salinity, mass transports, and
114 !! Layer thicknesses are read during initialization
115 !! Variables controlling some of the numerical considerations of offline transport
116 integer :: num_off_iter !< Number of advection iterations per offline step
117 integer :: num_vert_iter !< Number of vertical iterations per offline step
118 integer :: off_ale_mod !< Sets how frequently the ALE step is done during the advection
119 real :: dt_offline !< Timestep used for offline tracers [T ~> s]
120 real :: dt_offline_vertical !< Timestep used for calls to tracer vertical physics [T ~> s]
121 real :: evap_cfl_limit !< Limit on the fraction of the water that can be fluxed out of the top
122 !! layer in a timestep [nondim]. This is Copied from diabatic_CS controlling
123 !! how tracers follow freshwater fluxes
124 real :: minimum_forcing_depth !< The smallest depth over which fluxes can be applied [H ~> m or kg m-2].
125 !! This is copied from diabatic_CS controlling how tracers follow freshwater fluxes
126
127 real :: kd_max !< Runtime parameter specifying the maximum value of vertical diffusivity
128 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
129 real :: min_residual !< The minimum amount of total mass flux before exiting the main advection
130 !! routine [H L2 ~> m3 or kg]
131 !>@{ Diagnostic manager IDs for some fields that may be of interest when doing offline transport
132 integer :: &
133 id_uhr = -1, id_vhr = -1, &
134 ! Unused: id_ear = -1, id_ebr = -1, &
135 id_hr = -1, &
136 id_hdiff = -1, &
137 id_uhr_redist = -1, &
138 id_vhr_redist = -1, &
139 id_uhr_end = -1, &
140 id_vhr_end = -1, &
141 id_eta_pre_distribute = -1, &
142 id_eta_post_distribute = -1, &
143 id_h_redist = -1, &
144 id_eta_diff_end = -1
145
146 ! Diagnostic IDs for the regridded/remapped input fields
147 integer :: &
148 id_uhtr_regrid = -1, &
149 id_vhtr_regrid = -1, &
150 id_temp_regrid = -1, &
151 id_salt_regrid = -1, &
152 id_h_regrid = -1
153 !>@}
154
155 ! IDs for timings of various offline components
156 integer :: id_clock_read_fields = -1 !< A CPU time clock
157 integer :: id_clock_offline_diabatic = -1 !< A CPU time clock
158 integer :: id_clock_offline_adv = -1 !< A CPU time clock
159 integer :: id_clock_redistribute = -1 !< A CPU time clock
160
161 !> Zonal transport that may need to be stored between calls to step_MOM [H L2 ~> m3 or kg]
162 real, allocatable, dimension(:,:,:) :: uhtr
163 !> Meridional transport that may need to be stored between calls to step_MOM [H L2 ~> m3 or kg]
164 real, allocatable, dimension(:,:,:) :: vhtr
165
166 ! Fields at T-point
167 real, allocatable, dimension(:,:,:) :: eatr
168 !< Amount of fluid entrained from the layer above within
169 !! one time step [H ~> m or kg m-2]
170 real, allocatable, dimension(:,:,:) :: ebtr
171 !< Amount of fluid entrained from the layer below within
172 !! one time step [H ~> m or kg m-2]
173 ! Fields at T-points on interfaces
174 real, allocatable, dimension(:,:,:) :: kd !< Vertical diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
175 real, allocatable, dimension(:,:,:) :: h_end !< Thicknesses at the end of offline timestep [H ~> m or kg m-2]
176
177 real, allocatable, dimension(:,:) :: mld !< Mixed layer depths at thickness points [Z ~> m]
178
179 ! Allocatable arrays to read in entire fields during initialization
180 real, allocatable, dimension(:,:,:,:) :: uhtr_all !< Entire field of zonal transport [H L2 ~> m3 or kg]
181 real, allocatable, dimension(:,:,:,:) :: vhtr_all !< Entire field of meridional transport [H L2 ~> m3 or kg]
182 real, allocatable, dimension(:,:,:,:) :: hend_all !< Entire field of layer thicknesses [H ~> m or kg m-2]
183 real, allocatable, dimension(:,:,:,:) :: temp_all !< Entire field of temperatures [C ~> degC]
184 real, allocatable, dimension(:,:,:,:) :: salt_all !< Entire field of salinities [S ~> ppt]
185
187
201
202contains
203
204!> 3D advection is done by doing flux-limited nonlinear horizontal advection interspersed with an ALE
205!! regridding/remapping step. The loop in this routine is exited if remaining residual transports are below
206!! a runtime-specified value or a maximum number of iterations is reached.
207subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, CS, id_clock_ale, &
208 h_pre, uhtr, vhtr, converged)
209 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
210 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
211 real, intent(in) :: time_interval !< time interval covered by this call [T ~> s]
212 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
213 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
214 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
215 type(offline_transport_cs), pointer :: cs !< control structure for offline module
216 integer, intent(in) :: id_clock_ale !< Clock for ALE routines
217 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
218 intent(inout) :: h_pre !< layer thicknesses before advection
219 !! [H ~> m or kg m-2]
220 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
221 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
222 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
223 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
224 logical, intent( out) :: converged !< True if the iterations have converged
225
226 ! Local variables
227 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhtr_sub ! Substep zonal mass transports [H L2 ~> m3 or kg]
228 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub ! Substep meridional mass transports [H L2 ~> m3 or kg]
229
230 real :: prev_tot_residual, tot_residual ! Used to keep track of how close to convergence we are [H L2 ~> m3 or kg]
231
232 ! Variables used to keep track of layer thicknesses at various points in the code
233 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
234 h_new, & ! Updated layer thicknesses [H ~> m or kg m-2]
235 h_post_remap, & ! Layer thicknesses after remapping [H ~> m or kg m-2]
236 h_vol ! Layer volumes [H L2 ~> m3 or kg]
237 real :: dzregrid(szi_(g),szj_(g),szk_(gv)+1) ! The change in grid interface positions due to regridding,
238 ! in the same units as thicknesses [H ~> m or kg m-2]
239 integer :: niter, iter
240 real :: inum_iter ! The inverse of the number of iterations [nondim]
241 character(len=256) :: mesg ! The text of an error message
242 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
243 integer :: isdb, iedb, jsdb, jedb
244 logical :: x_before_y
245 real :: evap_cfl_limit ! Limit on the fraction of the water that can be fluxed out of the
246 ! top layer in a timestep [nondim]
247 real :: minimum_forcing_depth ! The smallest depth over which fluxes can be applied [H ~> m or kg m-2]
248 real :: dt_iter ! The timestep to use for each iteration [T ~> s]
249 real :: hl2_to_kg_scale ! Unit conversion factors to cell mass [kg H-1 L-2 ~> kg m-3 or 1]
250 character(len=20) :: debug_msg
251 call cpu_clock_begin(cs%id_clock_offline_adv)
252
253 ! Grid-related pointer assignments
254
255 x_before_y = cs%x_before_y
256
257 ! Initialize some shorthand variables from other structures
258 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
259 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
260 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
261
262 evap_cfl_limit = cs%evap_CFL_limit
263 minimum_forcing_depth = cs%minimum_forcing_depth
264
265 hl2_to_kg_scale = us%L_to_m**2*gv%H_to_kg_m2
266 niter = cs%num_off_iter
267 inum_iter = 1./real(niter)
268 dt_iter = cs%dt_offline*inum_iter
269
270 ! Initialize working arrays
271 h_new(:,:,:) = 0.0
272 h_vol(:,:,:) = 0.0
273 uhtr_sub(:,:,:) = 0.0
274 vhtr_sub(:,:,:) = 0.0
275
276 ! converged should only be true if there are no remaining mass fluxes
277 converged = .false.
278
279 ! Tracers are transported using the stored mass fluxes. Where possible, operators are Strang-split around
280 ! the call to
281 ! 1) Using the layer thicknesses and tracer concentrations from the previous timestep,
282 ! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to tracer_column_fns.
283 ! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline
284 ! 2) Half of the accumulated surface freshwater fluxes are applied
285 !! START ITERATION
286 ! 3) Accumulated mass fluxes are used to do horizontal transport. The number of iterations used in
287 ! advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are stored for later use
288 ! and resulting layer thicknesses fed into the next step
289 ! 4) Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for layers which might
290 ! 'vanish' because of horizontal mass transport to be 'reinflated'
291 ! 5) Check that transport is done if the remaining mass fluxes equals 0 or if the max number of iterations
292 ! has been reached
293 !! END ITERATION
294 ! 6) Repeat steps 1 and 2
295 ! 7) Force a remapping to the stored layer thicknesses that correspond to the snapshot of the online model
296 ! 8) Reset T/S and h to their stored snapshotted values to prevent model drift
297
298 ! Copy over the horizontal mass fluxes from the total mass fluxes
299 do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
300 uhtr_sub(i,j,k) = uhtr(i,j,k)
301 enddo ; enddo ; enddo
302 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
303 vhtr_sub(i,j,k) = vhtr(i,j,k)
304 enddo ; enddo ; enddo
305 do k=1,nz ; do j=js,je ; do i=is,ie
306 h_new(i,j,k) = h_pre(i,j,k)
307 enddo ; enddo ; enddo
308
309 if (cs%debug) then
310 call hchksum(h_pre, "h_pre before transport", g%HI, unscale=gv%H_to_MKS)
311 call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, g%HI, unscale=hl2_to_kg_scale)
312 endif
313 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
314 if (cs%print_adv_offline) then
315 write(mesg,'(A,ES24.16)') "Main advection starting transport: ", tot_residual*hl2_to_kg_scale
316 call mom_mesg(mesg)
317 endif
318
319 ! This loop does essentially a flux-limited, nonlinear advection scheme until all mass fluxes
320 ! are used. ALE is done after the horizontal advection.
321 do iter=1,cs%num_off_iter
322
323 do k=1,nz ; do j=js,je ; do i=is,ie
324 h_vol(i,j,k) = h_new(i,j,k) * g%areaT(i,j)
325 h_pre(i,j,k) = h_new(i,j,k)
326 enddo ; enddo ; enddo
327
328 if (cs%debug) then
329 call hchksum(h_vol, "h_vol before advect", g%HI, unscale=hl2_to_kg_scale)
330 call uvchksum("[uv]htr_sub before advect", uhtr_sub, vhtr_sub, g%HI, unscale=hl2_to_kg_scale)
331 write(debug_msg, '(A,I4.4)') 'Before advect ', iter
332 call mom_tracer_chkinv(debug_msg, g, gv, h_pre, cs%tracer_reg)
333 endif
334
335 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, cs%dt_offline, g, gv, us, &
336 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
337 max_iter_in=1, update_vol_prev=.true., uhr_out=uhtr, vhr_out=vhtr)
338
339 ! Switch the direction every iteration
340 x_before_y = .not. x_before_y
341
342 ! Update the new layer thicknesses after one round of advection has happened
343 do k=1,nz ; do j=js,je ; do i=is,ie
344 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
345 enddo ; enddo ; enddo
346
347 if (modulo(iter,cs%off_ale_mod)==0) then
348 ! Do ALE remapping/regridding to allow for more advection to occur in the next iteration
349 call pass_var(h_new,g%Domain)
350 if (cs%debug) then
351 call hchksum(h_new,"h_new before ALE", g%HI, unscale=gv%H_to_MKS)
352 write(debug_msg, '(A,I4.4)') 'Before ALE ', iter
353 call mom_tracer_chkinv(debug_msg, g, gv, h_new, cs%tracer_reg)
354 endif
355 call cpu_clock_begin(id_clock_ale)
356
357 call ale_update_regrid_weights(cs%dt_offline, cs%ALE_CSp)
358 call pre_ale_adjustments(g, gv, us, h_new, cs%tv, cs%tracer_Reg, cs%ALE_CSp)
359 ! Uncomment this to adjust the target grids for diagnostics, if there have been thickness
360 ! adjustments, but the offline tracer code does not yet have the other corresponding calls
361 ! that would be needed to support remapping its output.
362 ! call diag_update_remap_grids(CS%diag, alt_h=h_new)
363
364 call ale_regrid(g, gv, us, h_new, h_post_remap, dzregrid, cs%tv, cs%ALE_CSp)
365
366 ! Remap all variables from the old grid h_new onto the new grid h_post_remap
367 call ale_remap_tracers(cs%ALE_CSp, g, gv, h_new, h_post_remap, cs%tracer_Reg, &
368 cs%debug, dt=cs%dt_offline)
369 if (allocated(cs%tv%SpV_avg)) cs%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
370
371 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
372 h_new(i,j,k) = h_post_remap(i,j,k)
373 enddo ; enddo ; enddo
374 call cpu_clock_end(id_clock_ale)
375
376 if (cs%debug) then
377 call hchksum(h_new, "h_new after ALE", g%HI, unscale=gv%H_to_MKS)
378 write(debug_msg, '(A,I4.4)') 'After ALE ', iter
379 call mom_tracer_chkinv(debug_msg, g, gv, h_new, cs%tracer_reg)
380 endif
381 endif
382
383 do k=1,nz ; do j=js,je ; do i=is,ie
384 uhtr_sub(i,j,k) = uhtr(i,j,k)
385 vhtr_sub(i,j,k) = vhtr(i,j,k)
386 enddo ; enddo ; enddo
387 call pass_var(h_new, g%Domain)
388 call pass_vector(uhtr_sub, vhtr_sub, g%Domain)
389
390 ! Check for whether we've used up all the advection, or if we need to move on because
391 ! advection has stalled
392 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
393 if (cs%print_adv_offline) then
394 write(mesg,'(A,ES24.16)') "Main advection remaining transport: ", tot_residual*hl2_to_kg_scale
395 call mom_mesg(mesg)
396 endif
397 ! If all the mass transports have been used u, then quit
398 if (tot_residual == 0.0) then
399 write(mesg,*) "Converged after iteration ", iter
400 call mom_mesg(mesg)
401 converged = .true.
402 exit
403 endif
404 ! If advection has stalled or the remaining residual is less than a specified amount, quit
405 if ( (tot_residual == prev_tot_residual) .or. (tot_residual<cs%min_residual) ) then
406 converged = .false.
407 exit
408 endif
409
410 prev_tot_residual = tot_residual
411
412 enddo
413
414 ! Make sure that uhtr and vhtr halos are updated
415 h_pre(:,:,:) = h_new(:,:,:)
416 call pass_vector(uhtr, vhtr, g%Domain)
417
418 if (cs%debug) then
419 call hchksum(h_pre, "h after offline_advection_ale", g%HI, unscale=gv%H_to_MKS)
420 call uvchksum("[uv]htr after offline_advection_ale", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
421 call mom_tracer_chkinv("After offline_advection_ale", g, gv, h_pre, cs%tracer_reg)
422 endif
423
424 call cpu_clock_end(cs%id_clock_offline_adv)
425
426end subroutine offline_advection_ale
427
428!> In the case where the main advection routine did not converge, something needs to be done with the remaining
429!! transport. Two different ways are offered, 'barotropic' means that the residual is distributed equally
430!! throughout the water column. 'upwards' attempts to redistribute the transport in the layers above and will
431!! eventually work down the entire water column
432subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, converged)
433 type(offline_transport_cs), pointer :: cs !< control structure from initialize_MOM
434 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
435 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
436 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
437 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
438 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
439 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
440 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
441 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
442 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
443 logical, intent(in ) :: converged !< True if the iterations have converged
444
445 logical :: x_before_y
446 ! Variables used to keep track of layer thicknesses at various points in the code
447 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
448 h_new, & ! New layer thicknesses [H ~> m or kg m-2]
449 h_vol ! Cell volume [H L2 ~> m3 or kg]
450
451 ! Used to calculate the eta diagnostics
452 real, dimension(SZI_(G),SZJ_(G)) :: eta_work ! The total column thickness [H ~> m or kg m-2]
453 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhr !< Remaining zonal mass transport [H L2 ~> m3 or kg]
454 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhr !< Remaining meridional mass transport [H L2 ~> m3 or kg]
455
456 character(len=256) :: mesg ! The text of an error message
457 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, iter
458 real :: hl2_to_kg_scale ! Unit conversion factors to cell mass [kg H-1 L-2 ~> kg m-3 or 1]
459 real :: prev_tot_residual, tot_residual ! The absolute value of the remaining transports [H L2 ~> m3 or kg]
460
461 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
462 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
463
464 x_before_y = cs%x_before_y
465 hl2_to_kg_scale = us%L_to_m**2*gv%H_to_kg_m2
466
467 if (cs%id_eta_pre_distribute>0) then
468 eta_work(:,:) = 0.0
469 do k=1,nz ; do j=js,je ; do i=is,ie
470 if (h_pre(i,j,k) > gv%Angstrom_H) then
471 eta_work(i,j) = eta_work(i,j) + h_pre(i,j,k)
472 endif
473 enddo ; enddo ; enddo
474 call post_data(cs%id_eta_pre_distribute, eta_work, cs%diag)
475 endif
476
477 ! These are used to find out how much will be redistributed in this routine
478 if (cs%id_h_redist>0) call post_data(cs%id_h_redist, h_pre, cs%diag)
479 if (cs%id_uhr_redist>0) call post_data(cs%id_uhr_redist, uhtr, cs%diag)
480 if (cs%id_vhr_redist>0) call post_data(cs%id_vhr_redist, vhtr, cs%diag)
481
482 if (converged) return
483
484 if (cs%debug) then
485 call mom_tracer_chkinv("Before redistribute ", g, gv, h_pre, cs%tracer_reg)
486 endif
487
488 call cpu_clock_begin(cs%id_clock_redistribute)
489
490 if (cs%redistribute_upwards .or. cs%redistribute_barotropic) then
491 do iter = 1, cs%num_off_iter
492
493 ! Perform upwards redistribution
494 if (cs%redistribute_upwards) then
495
496 ! Calculate the layer volumes at beginning of redistribute
497 do k=1,nz ; do j=js,je ; do i=is,ie
498 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
499 enddo ; enddo ; enddo
500 call pass_var(h_vol,g%Domain)
501 call pass_vector(uhtr, vhtr, g%Domain)
502
503 if (cs%debug) then
504 call mom_tracer_chksum("Before upwards redistribute ", cs%tracer_Reg, g)
505 call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
506 endif
507
508 if (x_before_y) then
509 call distribute_residual_uh_upwards(g, gv, h_vol, uhtr)
510 call distribute_residual_vh_upwards(g, gv, h_vol, vhtr)
511 else
512 call distribute_residual_vh_upwards(g, gv, h_vol, vhtr)
513 call distribute_residual_uh_upwards(g, gv, h_vol, uhtr)
514 endif
515
516 call advect_tracer(h_pre, uhtr, vhtr, cs%OBC, cs%dt_offline, g, gv, us, &
517 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
518 max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr)
519
520 if (cs%debug) then
521 call mom_tracer_chksum("After upwards redistribute ", cs%tracer_Reg, g)
522 endif
523
524 ! Convert h_new back to layer thickness for ALE remapping
525 do k=1,nz ; do j=js,je ; do i=is,ie
526 uhtr(i,j,k) = uhr(i,j,k)
527 vhtr(i,j,k) = vhr(i,j,k)
528 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
529 h_pre(i,j,k) = h_new(i,j,k)
530 enddo ; enddo ; enddo
531
532 endif ! redistribute upwards
533
534 ! Perform barotropic redistribution
535 if (cs%redistribute_barotropic) then
536
537 ! Calculate the layer volumes at beginning of redistribute
538 do k=1,nz ; do j=js,je ; do i=is,ie
539 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
540 enddo ; enddo ; enddo
541 call pass_var(h_vol, g%Domain)
542 call pass_vector(uhtr, vhtr, g%Domain)
543
544 if (cs%debug) then
545 call mom_tracer_chksum("Before barotropic redistribute ", cs%tracer_Reg, g)
546 call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
547 endif
548
549 if (x_before_y) then
550 call distribute_residual_uh_barotropic(g, gv, h_vol, uhtr)
551 call distribute_residual_vh_barotropic(g, gv, h_vol, vhtr)
552 else
553 call distribute_residual_vh_barotropic(g, gv, h_vol, vhtr)
554 call distribute_residual_uh_barotropic(g, gv, h_vol, uhtr)
555 endif
556
557 call advect_tracer(h_pre, uhtr, vhtr, cs%OBC, cs%dt_offline, g, gv, us, &
558 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, &
559 max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr)
560
561 if (cs%debug) then
562 call mom_tracer_chksum("After barotropic redistribute ", cs%tracer_Reg, g)
563 endif
564
565 ! Convert h_new back to layer thickness for ALE remapping
566 do k=1,nz ; do j=js,je ; do i=is,ie
567 uhtr(i,j,k) = uhr(i,j,k)
568 vhtr(i,j,k) = vhr(i,j,k)
569 h_new(i,j,k) = h_vol(i,j,k) * g%IareaT(i,j)
570 h_pre(i,j,k) = h_new(i,j,k)
571 enddo ; enddo ; enddo
572
573 endif ! redistribute barotropic
574
575 ! Check to see if all transport has been exhausted
576 tot_residual = remaining_transport_sum(g, gv, us, uhtr, vhtr, h_new)
577 if (cs%print_adv_offline) then
578 write(mesg,'(A,ES24.16)') "Residual advection remaining transport: ", tot_residual*hl2_to_kg_scale
579 call mom_mesg(mesg)
580 endif
581 ! If the remaining residual is 0, then this return is done
582 if (tot_residual==0.0 ) then
583 exit
584 endif
585
586 if ( (tot_residual == prev_tot_residual) .or. (tot_residual<cs%min_residual) ) exit
587 prev_tot_residual = tot_residual
588
589 enddo ! Redistribution iteration
590 endif ! If one of the redistribution routines is requested
591
592 if (cs%id_eta_post_distribute>0) then
593 eta_work(:,:) = 0.0
594 do k=1,nz ; do j=js,je ; do i=is,ie
595 if (h_pre(i,j,k)>gv%Angstrom_H) then
596 eta_work(i,j) = eta_work(i,j) + h_pre(i,j,k)
597 endif
598 enddo ; enddo ; enddo
599 call post_data(cs%id_eta_post_distribute, eta_work, cs%diag)
600 endif
601
602 if (cs%id_uhr>0) call post_data(cs%id_uhr, uhtr, cs%diag)
603 if (cs%id_vhr>0) call post_data(cs%id_vhr, vhtr, cs%diag)
604
605 if (cs%debug) then
606 call hchksum(h_pre, "h_pre after redistribute", g%HI, unscale=gv%H_to_MKS)
607 call uvchksum("uhtr after redistribute", uhtr, vhtr, g%HI, unscale=hl2_to_kg_scale)
608 call mom_tracer_chkinv("after redistribute ", g, gv, h_new, cs%tracer_Reg)
609 endif
610
611 call cpu_clock_end(cs%id_clock_redistribute)
612
614
615!> Returns the sums of any non-negligible remaining transport [H L2 ~> m3 or kg] to check for advection convergence
616real function remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new)
617 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
618 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
619 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
620 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
621 intent(in ) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
622 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
623 intent(in ) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
624 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
625 intent(in ) :: h_new !< Layer thicknesses [H ~> m or kg m-2]
626
627 ! Local variables
628 real, dimension(SZI_(G),SZJ_(G)) :: trans_rem_col !< The vertical sum of the absolute value of
629 !! transports through the faces of a column [R Z L2 ~> kg].
630 real :: trans_cell !< The sum of the absolute value of the remaining transports through the faces
631 !! of a tracer cell [H L2 ~> m3 or kg]
632 integer :: i, j, k, is, ie, js, je, nz
633
634 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
635
636 trans_rem_col(:,:) = 0.0
637 do k=1,nz ; do j=js,je ; do i=is,ie
638 trans_cell = (abs(uhtr(i-1,j,k)) + abs(uhtr(i,j,k))) + &
639 (abs(vhtr(i,j-1,k)) + abs(vhtr(i,j,k)))
640 if (trans_cell > max(1.0e-16*h_new(i,j,k), gv%H_subroundoff) * g%areaT(i,j)) &
641 trans_rem_col(i,j) = trans_rem_col(i,j) + gv%H_to_RZ * trans_cell
642 enddo ; enddo ; enddo
643
644 ! The factor of 0.5 here is to avoid double-counting because two cells share a face.
645 remaining_transport_sum = 0.5 * gv%RZ_to_H * reproducing_sum(trans_rem_col, &
646 is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd), unscale=us%RZL2_to_kg)
647
648end function remaining_transport_sum
649
650!> The vertical/diabatic driver for offline tracers. First the eatr/ebtr associated with the interpolated
651!! vertical diffusivities are calculated and then any tracer column functions are done which can include
652!! vertical diffuvities and source/sink terms.
653subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, G, GV, US, CS, h_pre, tv, eatr, ebtr)
654
655 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
656 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
657 type(time_type), intent(in) :: time_end !< ending time of a segment, as a time type
658 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
659 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
660 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
661 type(offline_transport_cs), pointer :: cs !< control structure from initialize_MOM
662 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
663 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
664 type(thermo_var_ptrs), intent(in ) :: tv !< A structure pointing to various thermodynamic variables
665 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
666 intent(inout) :: eatr !< Entrainment from layer above [H ~> m or kg m-2]
667 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
668 intent(inout) :: ebtr !< Entrainment from layer below [H ~> m or kg m-2]
669
670 ! Local variables
671 real, dimension(SZI_(G),SZJ_(G)) :: &
672 sw, sw_vis, sw_nir !< Save old values of shortwave radiation [Q R Z T-1 ~> W m-2]
673 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Vertical distance across layers [Z ~> m]
674 real :: i_dzval ! An inverse distance between layer centers [Z-1 ~> m-1]
675 integer :: i, j, k, is, ie, js, je, nz
676 integer :: k_nonzero
677 real :: kd_bot ! Near-bottom diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
678 nz = gv%ke
679 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
680
681 call cpu_clock_begin(cs%id_clock_offline_diabatic)
682
683 call mom_mesg("Applying tracer source, sinks, and vertical mixing")
684
685 if (cs%debug) then
686 call hchksum(h_pre, "h_pre before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
687 call hchksum(eatr, "eatr before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
688 call hchksum(ebtr, "ebtr before offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
689 call mom_tracer_chkinv("Before offline_diabatic_ale", g, gv, h_pre, cs%tracer_reg)
690 endif
691
692 call thickness_to_dz(h_pre, tv, dz, g, gv, us)
693
694 eatr(:,:,:) = 0.
695 ebtr(:,:,:) = 0.
696 ! Calculate eatr and ebtr if vertical diffusivity is read
697 ! Because the saved remapped diagnostics from the online run assume a zero minimum thickness
698 ! but ALE may have a minimum thickness. Flood the diffusivities for all layers with the value
699 ! of Kd closest to the bottom which is non-zero
700 do j=js,je ; do i=is,ie
701 k_nonzero = nz+1
702 ! Find the nonzero bottom Kd
703 do k=nz+1,1,-1
704 if (cs%Kd(i,j,k)>0.) then
705 kd_bot = cs%Kd(i,j,k)
706 k_nonzero = k
707 exit
708 endif
709 enddo
710 ! Flood the bottom interfaces
711 do k=k_nonzero,nz+1
712 cs%Kd(i,j,k) = kd_bot
713 enddo
714 enddo ; enddo
715
716 do j=js,je ; do i=is,ie
717 eatr(i,j,1) = 0.
718 enddo ; enddo
719 do k=2,nz ; do j=js,je ; do i=is,ie
720 i_dzval = 1.0 / (gv%dZ_subroundoff + 0.5*(dz(i,j,k-1) + dz(i,j,k)))
721 eatr(i,j,k) = cs%dt_offline_vertical * i_dzval * cs%Kd(i,j,k)
722 ebtr(i,j,k-1) = eatr(i,j,k)
723 enddo ; enddo ; enddo
724 do j=js,je ; do i=is,ie
725 ebtr(i,j,nz) = 0.
726 enddo ; enddo
727
728 ! Add diurnal cycle for shortwave radiation (only used if run in ocean-only mode)
729 if (cs%diurnal_SW .and. cs%read_sw) then
730 sw(:,:) = fluxes%sw(:,:)
731 sw_vis(:,:) = fluxes%sw_vis_dir(:,:)
732 sw_nir(:,:) = fluxes%sw_nir_dir(:,:)
733 call offline_add_diurnal_sw(fluxes, g, time_start, time_end)
734 endif
735
736 if (associated(cs%optics)) &
737 call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, &
738 cs%opacity_CSp, cs%tracer_flow_CSp)
739
740 ! Note that tracerBoundaryFluxesInOut within this subroutine should NOT be called
741 ! as the freshwater fluxes have already been accounted for
742 call call_tracer_column_fns(h_pre, h_pre, eatr, ebtr, fluxes, cs%MLD, cs%dt_offline_vertical, &
743 g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
744
745 if (cs%diurnal_SW .and. cs%read_sw) then
746 fluxes%sw(:,:) = sw(:,:)
747 fluxes%sw_vis_dir(:,:) = sw_vis(:,:)
748 fluxes%sw_nir_dir(:,:) = sw_nir(:,:)
749 endif
750
751 if (cs%debug) then
752 call hchksum(h_pre, "h_pre after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
753 call hchksum(eatr, "eatr after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
754 call hchksum(ebtr, "ebtr after offline_diabatic_ale", g%HI, unscale=gv%H_to_MKS)
755 call mom_tracer_chkinv("After offline_diabatic_ale", g, gv, h_pre, cs%tracer_reg)
756 endif
757
758 call cpu_clock_end(cs%id_clock_offline_diabatic)
759
760end subroutine offline_diabatic_ale
761
762!> Apply positive freshwater fluxes (into the ocean) and update netMassOut with only the negative
763!! (out of the ocean) fluxes
764subroutine offline_fw_fluxes_into_ocean(G, GV, CS, fluxes, h, in_flux_optional)
765 type(offline_transport_cs), intent(inout) :: cs !< Offline control structure
766 type(ocean_grid_type), intent(in) :: g !< Grid structure
767 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
768 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
769 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
770 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
771 real, dimension(SZI_(G),SZJ_(G)), &
772 optional, intent(in) :: in_flux_optional !< The total time-integrated amount
773 !! of tracer that leaves with freshwater
774 !! [CU H ~> Conc m or Conc kg m-2]
775
776 integer :: i, j, m
777 real, dimension(SZI_(G),SZJ_(G)) :: negative_fw !< store all negative fluxes [H ~> m or kg m-2]
778 logical :: update_h !< Flag for whether h should be updated
779
780 if ( present(in_flux_optional) ) &
781 call mom_error(warning, "Positive freshwater fluxes with non-zero tracer concentration not supported yet")
782
783 ! Set all fluxes to 0
784 negative_fw(:,:) = 0.
785
786 ! Sort fluxes into positive and negative
787 do j=g%jsc,g%jec ; do i=g%isc,g%iec
788 if (fluxes%netMassOut(i,j)<0.0) then
789 negative_fw(i,j) = fluxes%netMassOut(i,j)
790 fluxes%netMassOut(i,j) = 0.
791 endif
792 enddo ; enddo
793
794 if (cs%debug) then
795 call hchksum(h, "h before fluxes into ocean", g%HI, unscale=gv%H_to_MKS)
796 call mom_tracer_chkinv("Before fluxes into ocean", g, gv, h, cs%tracer_reg)
797 endif
798 do m = 1,cs%tracer_reg%ntr
799 ! Layer thicknesses should only be updated after the last tracer is finished
800 update_h = ( m == cs%tracer_reg%ntr )
801 call applytracerboundaryfluxesinout(g, gv, cs%tracer_reg%tr(m)%t, cs%dt_offline, fluxes, h, &
802 cs%evap_CFL_limit, cs%minimum_forcing_depth, update_h_opt=update_h)
803 enddo
804 if (cs%debug) then
805 call hchksum(h, "h after fluxes into ocean", g%HI, unscale=gv%H_to_MKS)
806 call mom_tracer_chkinv("After fluxes into ocean", g, gv, h, cs%tracer_reg)
807 endif
808
809 ! Now that fluxes into the ocean are done, save the negative fluxes for later
810 fluxes%netMassOut(:,:) = negative_fw(:,:)
811
812end subroutine offline_fw_fluxes_into_ocean
813
814!> Apply negative freshwater fluxes (out of the ocean)
815subroutine offline_fw_fluxes_out_ocean(G, GV, CS, fluxes, h, out_flux_optional)
816 type(offline_transport_cs), intent(inout) :: cs !< Offline control structure
817 type(ocean_grid_type), intent(in) :: g !< Grid structure
818 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
819 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
820 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
821 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
822 real, dimension(SZI_(G),SZJ_(G)), &
823 optional, intent(in) :: out_flux_optional !< The total time-integrated amount
824 !! of tracer that leaves with freshwater
825 !! [CU H ~> Conc m or Conc kg m-2]
826
827 integer :: m
828 logical :: update_h !< Flag for whether h should be updated
829
830 if ( present(out_flux_optional) ) &
831 call mom_error(warning, "Negative freshwater fluxes with non-zero tracer concentration not supported yet")
832
833 if (cs%debug) then
834 call hchksum(h, "h before fluxes out of ocean", g%HI, unscale=gv%H_to_MKS)
835 call mom_tracer_chkinv("Before fluxes out of ocean", g, gv, h, cs%tracer_reg)
836 endif
837 do m = 1, cs%tracer_reg%ntr
838 ! Layer thicknesses should only be updated after the last tracer is finished
839 update_h = ( m == cs%tracer_reg%ntr )
840 call applytracerboundaryfluxesinout(g, gv, cs%tracer_reg%tr(m)%t, cs%dt_offline, fluxes, h, &
841 cs%evap_CFL_limit, cs%minimum_forcing_depth, update_h_opt = update_h)
842 enddo
843 if (cs%debug) then
844 call hchksum(h, "h after fluxes out of ocean", g%HI, unscale=gv%H_to_MKS)
845 call mom_tracer_chkinv("Before fluxes out of ocean", g, gv, h, cs%tracer_reg)
846 endif
847
848end subroutine offline_fw_fluxes_out_ocean
849
850!> When in layer mode, 3D horizontal advection using stored mass fluxes must be used. Horizontal advection is
851!! done via tracer_advect, whereas the vertical component is actually handled by vertdiff in tracer_column_fns
852subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, CS, h_pre, eatr, ebtr, uhtr, vhtr)
853 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
854 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
855 real, intent(in) :: time_interval !< Offline transport time interval [T ~> s]
856 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
857 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
858 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
859 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
860 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
861 intent(inout) :: h_pre !< layer thicknesses before advection [H ~> m or kg m-2]
862 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
863 intent(inout) :: eatr !< Entrainment from layer above [H ~> m or kg m-2]
864 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
865 intent(inout) :: ebtr !< Entrainment from layer below [H ~> m or kg m-2]
866 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
867 intent(inout) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg]
868 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
869 intent(inout) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg]
870
871 ! Local variables
872
873 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uhtr_sub ! Remaining zonal mass transports [H L2 ~> m3 or kg]
874 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub ! Remaining meridional mass transports [H L2 ~> m3 or kg]
875
876 real, dimension(SZI_(G),SZJB_(G)) :: rem_col_flux ! The summed absolute value of the remaining
877 ! mass fluxes through the faces of a column or within a column [R Z L2 ~> kg]
878 real :: sum_flux ! Globally summed absolute value of fluxes [R Z L2 ~> kg], which is
879 ! used to keep track of how close to convergence we are.
880
881 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
882 eatr_sub, & ! Layer entrainment rate from above for this sub-cycle [H ~> m or kg m-2]
883 ebtr_sub ! Layer entrainment rate from below for this sub-cycle [H ~> m or kg m-2]
884 ! Variables used to keep track of layer thicknesses at various points in the code
885 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
886 h_new, & ! Updated thicknesses [H ~> m or kg m-2]
887 h_vol ! Cell volumes [H L2 ~> m3 or kg]
888 ! Work arrays for temperature and salinity
889 integer :: iter
890 real :: dt_iter ! The timestep of each iteration [T ~> s]
891 character(len=160) :: mesg ! The text of an error message
892 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
893 integer :: isdb, iedb, jsdb, jedb
894 logical :: z_first, x_before_y
895
896 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
897 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
898 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
899
900 dt_iter = time_interval / real(max(1, cs%num_off_iter))
901 x_before_y = cs%x_before_y
902
903 do iter=1,cs%num_off_iter
904
905 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
906 eatr_sub(i,j,k) = eatr(i,j,k)
907 ebtr_sub(i,j,k) = ebtr(i,j,k)
908 enddo ; enddo ; enddo
909
910 do k=1,nz ; do j=js-1,je+1 ; do i=is-2,ie+1
911 uhtr_sub(i,j,k) = uhtr(i,j,k)
912 enddo ; enddo ; enddo
913
914 do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
915 vhtr_sub(i,j,k) = vhtr(i,j,k)
916 enddo ; enddo ; enddo
917
918 ! Calculate 3d mass transports to be used in this iteration
919 call limit_mass_flux_3d(g, gv, uhtr_sub, vhtr_sub, eatr_sub, ebtr_sub, h_pre)
920
921 if (z_first) then
922 ! First do vertical advection
923 call update_h_vertical_flux(g, gv, eatr_sub, ebtr_sub, h_pre, h_new)
924 call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, &
925 fluxes, cs%mld, dt_iter, g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
926 ! We are now done with the vertical mass transports, so now h_new is h_sub
927 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
928 h_pre(i,j,k) = h_new(i,j,k)
929 enddo ; enddo ; enddo
930 call pass_var(h_pre,g%Domain)
931
932 ! Second zonal and meridional advection
933 call update_h_horizontal_flux(g, gv, uhtr_sub, vhtr_sub, h_pre, h_new)
934 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
935 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
936 enddo ; enddo ; enddo
937 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, dt_iter, g, gv, us, &
938 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30)
939
940 ! Done with horizontal so now h_pre should be h_new
941 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
942 h_pre(i,j,k) = h_new(i,j,k)
943 enddo ; enddo ; enddo
944
945 endif
946
947 if (.not. z_first) then
948
949 ! First zonal and meridional advection
950 call update_h_horizontal_flux(g, gv, uhtr_sub, vhtr_sub, h_pre, h_new)
951 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
952 h_vol(i,j,k) = h_pre(i,j,k) * g%areaT(i,j)
953 enddo ; enddo ; enddo
954 call advect_tracer(h_pre, uhtr_sub, vhtr_sub, cs%OBC, dt_iter, g, gv, us, cs%tracer_adv_CSp, &
955 cs%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30)
956
957 ! Done with horizontal so now h_pre should be h_new
958 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
959 h_pre(i,j,k) = h_new(i,j,k)
960 enddo ; enddo ; enddo
961
962 ! Second vertical advection
963 call update_h_vertical_flux(g, gv, eatr_sub, ebtr_sub, h_pre, h_new)
964 call call_tracer_column_fns(h_pre, h_new, eatr_sub, ebtr_sub, &
965 fluxes, cs%mld, dt_iter, g, gv, us, cs%tv, cs%optics, cs%tracer_flow_CSp, cs%debug)
966 ! We are now done with the vertical mass transports, so now h_new is h_sub
967 do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1
968 h_pre(i,j,k) = h_new(i,j,k)
969 enddo ; enddo ; enddo
970
971 endif
972
973 ! Update remaining transports
974 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
975 eatr(i,j,k) = eatr(i,j,k) - eatr_sub(i,j,k)
976 ebtr(i,j,k) = ebtr(i,j,k) - ebtr_sub(i,j,k)
977 enddo ; enddo ; enddo
978
979 do k=1,nz ; do j=js-1,je+1 ; do i=is-2,ie+1
980 uhtr(i,j,k) = uhtr(i,j,k) - uhtr_sub(i,j,k)
981 enddo ; enddo ; enddo
982
983 do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
984 vhtr(i,j,k) = vhtr(i,j,k) - vhtr_sub(i,j,k)
985 enddo ; enddo ; enddo
986
987 call pass_var(eatr,g%Domain)
988 call pass_var(ebtr,g%Domain)
989 call pass_var(h_pre,g%Domain)
990 call pass_vector(uhtr,vhtr,g%Domain)
991
992 ! Calculate how close we are to converging by summing the remaining fluxes at each point
993 rem_col_flux(:,:) = 0.0
994 do k=1,nz ; do j=js,je ; do i=is,ie
995 rem_col_flux(i,j) = rem_col_flux(i,j) + gv%H_to_RZ * &
996 ( (abs(eatr(i,j,k)) + abs(ebtr(i,j,k))) + &
997 ((abs(uhtr(i-1,j,k)) + abs(uhtr(i,j,k))) + &
998 (abs(vhtr(i,j-1,k)) + abs(vhtr(i,j,k))) ) )
999 enddo ; enddo ; enddo
1000 sum_flux = reproducing_sum(rem_col_flux, is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd), &
1001 unscale=us%RZL2_to_kg)
1002
1003 if (sum_flux==0) then
1004 write(mesg,*) 'offline_advection_layer: Converged after iteration', iter
1005 call mom_mesg(mesg)
1006 exit
1007 else
1008 write(mesg,*) "offline_advection_layer: Iteration ", iter, " remaining total fluxes: ", sum_flux*us%RZL2_to_kg
1009 call mom_mesg(mesg)
1010 endif
1011
1012 ! Switch order of Strang split every iteration
1013 z_first = .not. z_first
1014 x_before_y = .not. x_before_y
1015 enddo
1016
1017end subroutine offline_advection_layer
1018
1019!> Update fields used in this round of offline transport. First fields are updated from files or from arrays
1020!! read during initialization. Then if in an ALE-dependent coordinate, regrid/remap fields.
1021subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
1022 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
1023 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
1024 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1025 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1026 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1027 intent(inout) :: h !< The regridded layer thicknesses [H ~> m or kg m-2]
1028 type(forcing), intent(inout) :: fluxes !< Pointers to forcing fields
1029 logical, intent(in ) :: do_ale !< True if using ALE
1030 ! Local variables
1031 integer :: stencil
1032 integer :: i, j, k, is, ie, js, je, nz
1033 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_start ! Initial thicknesses [H ~> m or kg m-2]
1034 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1035
1036 call cpu_clock_begin(cs%id_clock_read_fields)
1037 call calltree_enter("update_offline_fields, MOM_offline_main.F90")
1038
1039 if (cs%debug) then
1040 call uvchksum("[uv]htr before update_offline_fields", cs%uhtr, cs%vhtr, g%HI, &
1041 unscale=us%L_to_m**2*gv%H_to_kg_m2)
1042 call hchksum(cs%h_end, "h_end before update_offline_fields", g%HI, unscale=gv%H_to_MKS)
1043 call hchksum(cs%tv%T, "Temp before update_offline_fields", g%HI, unscale=us%C_to_degC)
1044 call hchksum(cs%tv%S, "Salt before update_offline_fields", g%HI, unscale=us%S_to_ppt)
1045 endif
1046
1047 ! Store a copy of the layer thicknesses before ALE regrid/remap
1048 h_start(:,:,:) = h(:,:,:)
1049
1050 ! Most fields will be read in from files
1051 call update_offline_from_files( g, gv, us, cs%nk_input, cs%mean_file, cs%sum_file, cs%snap_file, &
1052 cs%surf_file, cs%h_end, cs%uhtr, cs%vhtr, cs%tv%T, cs%tv%S, &
1053 cs%mld, cs%Kd, fluxes, cs%ridx_sum, cs%ridx_snap, &
1054 cs%read_mld, cs%mld_var_name, cs%read_sw, &
1055 .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 ! Unused: CS%id_ear = register_diag_field('ocean_model', 'ear', diag%axesTL, Time, &
1181 ! 'Remaining thickness entrained from above', 'm')
1182 ! Unused: 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 if (cs%read_mld) then
1417 cs%Hmix_fixed = 0.0
1418 else
1419 call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix_fixed, &
1420 "The prescribed depth over which the near-surface viscosity and "//&
1421 "diffusivity are elevated when the bulk mixed layer is not used.", &
1422 units="m", scale=us%m_to_Z, fail_if_missing=.true.)
1423 endif
1424
1425 call get_param(param_file, mdl, "OFFLINE_ADD_DIURNAL_SW", cs%diurnal_sw, &
1426 "Adds a synthetic diurnal cycle in the same way that the ice "//&
1427 "model would have when time-averaged fields of shortwave "//&
1428 "radiation are read in", default=.false.)
1429 call get_param(param_file, mdl, "KD_MAX", cs%Kd_max, &
1430 "The maximum permitted increment for the diapycnal "//&
1431 "diffusivity from TKE-based parameterizations, or a "//&
1432 "negative value for no limit.", units="m2 s-1", default=-1.0, scale=gv%m2_s_to_HZ_T)
1433 call get_param(param_file, mdl, "MIN_RESIDUAL_TRANSPORT", cs%min_residual, &
1434 "How much remaining transport before the main offline advection is exited. "//&
1435 "The default value corresponds to about 1 meter of difference in a grid cell", &
1436 default=1.e9, units="m3", scale=gv%m_to_H*us%m_to_L**2)
1437 call get_param(param_file, mdl, "READ_ALL_TS_UVH", cs%read_all_ts_uvh, &
1438 "Reads all time levels of a subset of the fields necessary to run " // &
1439 "the model offline. This can require a large amount of memory "// &
1440 "and will make initialization very slow. However, for offline "// &
1441 "runs spanning more than a year this can reduce total I/O overhead", &
1442 default=.false.)
1443
1444 ! Concatenate offline directory and file names
1445 cs%snap_file = trim(cs%offlinedir)//trim(cs%snap_file)
1446 cs%mean_file = trim(cs%offlinedir)//trim(cs%mean_file)
1447 cs%sum_file = trim(cs%offlinedir)//trim(cs%sum_file)
1448 cs%surf_file = trim(cs%offlinedir)//trim(cs%surf_file)
1449
1450 cs%num_vert_iter = cs%dt_offline / cs%dt_offline_vertical
1451
1452 ! Map redistribute_method onto logicals in CS
1453 select case (redistribute_method)
1454 case ('barotropic')
1455 cs%redistribute_barotropic = .true.
1456 cs%redistribute_upwards = .false.
1457 case ('upwards')
1458 cs%redistribute_barotropic = .false.
1459 cs%redistribute_upwards = .true.
1460 case ('both')
1461 cs%redistribute_barotropic = .true.
1462 cs%redistribute_upwards = .true.
1463 case ('none')
1464 cs%redistribute_barotropic = .false.
1465 cs%redistribute_upwards = .false.
1466 end select
1467
1468 ! Set the accumulated time to zero
1469 cs%accumulated_time = real_to_time(0.0)
1470 cs%vertical_time = cs%accumulated_time
1471 ! Set the starting read index for time-averaged and snapshotted fields
1472 cs%ridx_sum = cs%start_index
1473 if (cs%fields_are_offset) cs%ridx_snap = next_modulo_time(cs%start_index,cs%numtime)
1474 if (.not. cs%fields_are_offset) cs%ridx_snap = cs%start_index
1475
1476 ! Copy members from other modules
1477 call extract_diabatic_member(diabatic_csp, opacity_csp=cs%opacity_CSp, optics_csp=cs%optics, &
1478 diabatic_aux_csp=cs%diabatic_aux_CSp, &
1479 evap_cfl_limit=cs%evap_CFL_limit, &
1480 minimum_forcing_depth=cs%minimum_forcing_depth)
1481
1482 ! Allocate arrays
1483 allocate(cs%uhtr(isdb:iedb,jsd:jed,nz), source=0.0)
1484 allocate(cs%vhtr(isd:ied,jsdb:jedb,nz), source=0.0)
1485 allocate(cs%eatr(isd:ied,jsd:jed,nz), source=0.0)
1486 allocate(cs%ebtr(isd:ied,jsd:jed,nz), source=0.0)
1487 allocate(cs%h_end(isd:ied,jsd:jed,nz), source=0.0)
1488 allocate(cs%Kd(isd:ied,jsd:jed,nz+1), source=0.0)
1489 allocate(cs%mld(g%isd:g%ied,g%jsd:g%jed), source=cs%Hmix_fixed)
1490
1491 if (cs%read_all_ts_uvh) then
1492 call read_all_input(cs, g, gv, us)
1493 endif
1494
1495 ! Initialize ids for clocks used in offline routines
1496 cs%id_clock_read_fields = cpu_clock_id('(Offline read fields)',grain=clock_module)
1497 cs%id_clock_offline_diabatic = cpu_clock_id('(Offline diabatic)',grain=clock_module)
1498 cs%id_clock_offline_adv = cpu_clock_id('(Offline transport)',grain=clock_module)
1499 cs%id_clock_redistribute = cpu_clock_id('(Offline redistribute)',grain=clock_module)
1500
1501 call calltree_leave("offline_transport_init")
1502
1503end subroutine offline_transport_init
1504
1505!> Coordinates the allocation and reading in all time levels of uh, vh, hend, temp, and salt from files. Used
1506!! when read_all_ts_uvh
1507subroutine read_all_input(CS, G, GV, US)
1508 type(offline_transport_cs), intent(inout) :: CS !< Control structure for offline module
1509 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1510 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1511 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1512
1513 integer :: isd, ied, jsd, jed, nz, t, ntime
1514 integer :: IsdB, IedB, JsdB, JedB
1515
1516 nz = gv%ke ; ntime = cs%numtime
1517 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1518 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1519
1520 ! Extra safety check that we're not going to overallocate any arrays
1521 if (cs%read_all_ts_uvh) then
1522 if (allocated(cs%uhtr_all)) call mom_error(fatal, "uhtr_all is already allocated")
1523 if (allocated(cs%vhtr_all)) call mom_error(fatal, "vhtr_all is already allocated")
1524 if (allocated(cs%hend_all)) call mom_error(fatal, "hend_all is already allocated")
1525 if (allocated(cs%temp_all)) call mom_error(fatal, "temp_all is already allocated")
1526 if (allocated(cs%salt_all)) call mom_error(fatal, "salt_all is already allocated")
1527
1528 allocate(cs%uhtr_all(isdb:iedb,jsd:jed,nz,ntime), source=0.0)
1529 allocate(cs%vhtr_all(isd:ied,jsdb:jedb,nz,ntime), source=0.0)
1530 allocate(cs%hend_all(isd:ied,jsd:jed,nz,ntime), source=0.0)
1531 allocate(cs%temp_all(isd:ied,jsd:jed,nz,1:ntime), source=0.0)
1532 allocate(cs%salt_all(isd:ied,jsd:jed,nz,1:ntime), source=0.0)
1533
1534 call mom_mesg("Reading in uhtr, vhtr, h_start, h_end, temp, salt")
1535 do t = 1,ntime
1536 call mom_read_vector(cs%snap_file, 'uhtr_sum', 'vhtr_sum', cs%uhtr_all(:,:,1:cs%nk_input,t), &
1537 cs%vhtr_all(:,:,1:cs%nk_input,t), g%Domain, timelevel=t, &
1538 scale=us%m_to_L**2*gv%kg_m2_to_H)
1539 call mom_read_data(cs%snap_file,'h_end', cs%hend_all(:,:,1:cs%nk_input,t), g%Domain, &
1540 timelevel=t, position=center, scale=gv%kg_m2_to_H)
1541 call mom_read_data(cs%mean_file,'temp', cs%temp_all(:,:,1:cs%nk_input,t), g%Domain, &
1542 timelevel=t, position=center, scale=us%degC_to_C)
1543 call mom_read_data(cs%mean_file,'salt', cs%salt_all(:,:,1:cs%nk_input,t), g%Domain, &
1544 timelevel=t, position=center, scale=us%ppt_to_S)
1545 enddo
1546 endif
1547
1548end subroutine read_all_input
1549
1550!> Deallocates (if necessary) arrays within the offline control structure
1551subroutine offline_transport_end(CS)
1552 type(offline_transport_cs), pointer :: cs !< Control structure for offline module
1553
1554 ! Explicitly allocate all allocatable arrays
1555 deallocate(cs%uhtr)
1556 deallocate(cs%vhtr)
1557 deallocate(cs%eatr)
1558 deallocate(cs%ebtr)
1559 deallocate(cs%h_end)
1560 deallocate(cs%Kd)
1561 if (allocated(cs%mld)) deallocate(cs%mld)
1562 if (cs%read_all_ts_uvh) then
1563 deallocate(cs%uhtr_all)
1564 deallocate(cs%vhtr_all)
1565 deallocate(cs%hend_all)
1566 deallocate(cs%temp_all)
1567 deallocate(cs%salt_all)
1568 endif
1569
1570 deallocate(cs)
1571
1572end subroutine offline_transport_end
1573
1574!> \namespace mom_offline_main
1575!! \section offline_overview Offline Tracer Transport in MOM6
1576!! 'Offline tracer modeling' uses physical fields (e.g. mass transports and layer thicknesses) saved
1577!! from a previous integration of the physical model to transport passive tracers. These fields are
1578!! accumulated or averaged over a period of time (in this test case, 1 day) and used to integrate
1579!! portions of the MOM6 code base that handle the 3d advection and diffusion of passive tracers.
1580!!
1581!! The distribution of tracers in the ocean modeled offline should not be expected to match an online
1582!! simulation. Accumulating transports over more than one online model timestep implicitly assumes
1583!! homogeneity over that time period and essentially aliases over processes that occur with higher
1584!! frequency. For example, consider the case of a surface boundary layer with a strong diurnal cycle.
1585!! An offline simulation with a 1 day timestep, captures the net transport into or out of that layer,
1586!! but not the exact cycling. This effective aliasing may also complicate online model configurations
1587!! which strongly-eddying regions. In this case, the offline model timestep must be limited to some
1588!! fraction of the eddy correlation timescale. Lastly, the nonlinear advection scheme which applies
1589!! limited mass-transports over a sequence of iterations means that tracers are not transported along
1590!! exactly the same path as they are in the online model.
1591!!
1592!! This capability has currently targeted the Baltic_ALE_z test case, though some work has also been
1593!! done with the OM4 1/2 degree configuration. Work is ongoing to develop recommendations and best
1594!! practices for investigators seeking to use MOM6 for offline tracer modeling.
1595!!
1596!! \section offline_technical Implementation of offline routine in MOM6
1597!!
1598!! The subroutine step_tracers that coordinates this can be found in MOM.F90 and is only called
1599!! using the solo ocean driver. This is to avoid issues with coupling to other climate components
1600!! that may be relying on fluxes from the ocean to be coupled more often than the offline time step.
1601!! Other routines related to offline tracer modeling can be found in tracers/MOM_offline_control.F90
1602!!
1603!! As can also be seen in the comments for the step_tracers subroutine, an offline time step
1604!! comprises the following steps:
1605!! -# Using the layer thicknesses and tracer concentrations from the previous timestep,
1606!! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to
1607!! tracer_column_fns.
1608!! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline
1609!! -# Half of the accumulated surface freshwater fluxes are applied
1610!! START ITERATION
1611!! -# Accumulated mass fluxes are used to do horizontal transport. The number of iterations
1612!! used in advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are
1613!! stored for later use and resulting layer thicknesses fed into the next step
1614!! -# Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for
1615!! layers which might 'vanish' because of horizontal mass transport to be 'reinflated'
1616!! and essentially allows for the vertical transport of tracers
1617!! -# Check that transport is done if the remaining mass fluxes equals 0 or if the max
1618!! number of iterations has been reached
1619!! END ITERATION
1620!! -# Repeat steps 1 and 2
1621!! -# Redistribute any residual mass fluxes that remain after the advection iterations
1622!! in a barotropic manner, progressively upward through the water column.
1623!! -# Force a remapping to the stored layer thicknesses that correspond to the snapshot of
1624!! the online model at the end of an accumulation interval
1625!! -# Reset T/S and h to their stored snapshotted values to prevent model drift
1626!!
1627!! \section offline_evaluation Evaluating the utility of an offline tracer model
1628!! How well an offline tracer model can be used as an alternative to integrating tracers online
1629!! with the prognostic model must be evaluated for each application. This efficacy may be related
1630!! to the native coordinate of the online model, to the length of the offline timestep, and to the
1631!! behavior of the tracer itself.
1632!!
1633!! A framework for formally regression testing the offline capability still needs to be developed.
1634!! However, as a simple way of testing whether the offline model is nominally behaving as expected,
1635!! the total inventory of the advection test tracers (tr1, tr2, etc.) should be conserved between
1636!! time steps except for the last 4 decimal places. As a general guideline, an offline timestep of
1637!! 5 days or less.
1638!!
1639!! \section offline_parameters Runtime parameters for offline tracers
1640!! - OFFLINEDIR: Input directory where the offline fields can be found
1641!! - OFF_SUM_FILE: Filename where the accumulated fields can be found (e.g. horizontal mass transports)
1642!! - OFF_SNAP_FILE: Filename where snapshot fields can be found (e.g. end of timestep layer thickness)
1643!! - START_INDEX: Which timelevel of the input files to read first
1644!! - NUMTIME: How many timelevels to read before 'looping' back to 1
1645!! - FIELDS_ARE_OFFSET: True if the time-averaged fields and snapshot fields are offset by one
1646!! time level, probably not needed
1647!! -NUM_OFF_ITER: Maximum number of iterations to do for the nonlinear advection scheme
1648!! -REDISTRIBUTE_METHOD: Redistributes any remaining horizontal fluxes throughout the rest of water column.
1649!! Options are 'barotropic' which "evenly distributes flux throughout the entire water
1650!! column,'upwards' which adds the maximum of the remaining flux in each layer above,
1651!! and 'none' which does no redistribution"
1652
1653end module mom_offline_main
1654