Line data Source code
1 : !=======================================================================
2 : ! Transports quantities using the second-order conservative remapping
3 : ! scheme developed by John Dukowicz and John Baumgardner (DB) and modified
4 : ! for sea ice by William Lipscomb and Elizabeth Hunke.
5 : !
6 : ! References:
7 : !
8 : ! Dukowicz, J. K., and J. R. Baumgardner, 2000: Incremental
9 : ! remapping as a transport/advection algorithm, J. Comput. Phys.,
10 : ! 160, 318-335.
11 : !
12 : ! Lipscomb, W. H., and E. C. Hunke, 2004: Modeling sea ice
13 : ! transport using incremental remapping, Mon. Wea. Rev., 132,
14 : ! 1341-1354.
15 : !
16 : ! authors William H. Lipscomb, LANL
17 : ! John Baumgardner, LANL
18 : !
19 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
20 : ! 2004-05: Block structure added (WHL)
21 : ! 2006: Moved remap driver to ice_transport_driver
22 : ! Geometry changes:
23 : ! (1) Reconstruct fields in stretched logically rectangular coordinates
24 : ! (2) Modify geometry so that the area flux across each edge
25 : ! can be specified (following an idea of Mats Bentsen)
26 : ! 2010: ECH removed unnecessary grid arrays and optional arguments from
27 : ! horizontal_remap
28 :
29 : module ice_transport_remap
30 :
31 : use ice_kinds_mod
32 : use ice_blocks, only: nx_block, ny_block
33 : use ice_calendar, only: istep1
34 : use ice_communicate, only: my_task
35 : use ice_constants, only: c0, c1, c2, c12, p333, p4, p5, p6, &
36 : eps13, eps16, & ! LCOV_EXCL_LINE
37 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
38 : field_loc_NEcorner, field_type_vector
39 : use ice_diagnostics, only: diagnostic_abort
40 : use ice_domain_size, only: max_blocks, ncat
41 : use ice_fileunits, only: nu_diag
42 : use ice_exit, only: abort_ice
43 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
44 : use icepack_intfc, only: icepack_query_parameters
45 : use ice_grid, only : grid_ice
46 :
47 : implicit none
48 : private
49 : public :: init_remap, horizontal_remap, make_masks
50 :
51 : integer (kind=int_kind), parameter :: &
52 : ngroups = 6 ,&! number of groups of triangles that ! LCOV_EXCL_LINE
53 : ! contribute transports across each edge
54 : nvert = 3 ! number of vertices in a triangle
55 :
56 : ! for triangle integral formulas
57 : real (kind=dbl_kind), parameter :: &
58 : p5625m = -9._dbl_kind/16._dbl_kind ,& ! LCOV_EXCL_LINE
59 : p52083 = 25._dbl_kind/48._dbl_kind
60 :
61 : logical :: &
62 : l_fixed_area ! if true, prescribe area flux across each edge
63 : ! if false, area flux is determined internally
64 : ! and is passed out
65 :
66 : logical (kind=log_kind), parameter :: bugcheck = .false.
67 :
68 : !=======================================================================
69 : ! Here is some information about how the incremental remapping scheme
70 : ! works in CICE and how it can be adapted for use in other models.
71 : !
72 : ! The remapping routine is designed to transport a generic mass-like
73 : ! field (in CICE, the ice fractional area) along with an arbitrary number
74 : ! of tracers in two dimensions. The velocity components are assumed
75 : ! to lie at grid cell corners and the transported scalars at cell centers.
76 : ! Incremental remapping has the following desirable properties:
77 : !
78 : ! (1) Tracer monotonicity is preserved. That is, no new local
79 : ! extrema are produced in fields like ice thickness or internal
80 : ! energy.
81 : ! (2) The reconstucted mass and tracer fields vary linearly in x and y.
82 : ! This means that remapping is 2nd-order accurate in space,
83 : ! except where horizontal gradients are limited to preserve
84 : ! monotonicity.
85 : ! (3) There are economies of scale. Transporting a single field
86 : ! is rather expensive, but additional fields have a relatively
87 : ! low marginal cost.
88 : !
89 : ! The following generic conservation equations may be solved:
90 : !
91 : ! dm/dt = del*(u*m) (0)
92 : ! d(m*T1)/dt = del*(u*m*T1) (1)
93 : ! d(m*T1*T2)/dt = del*(u*m*T1*T2) (2)
94 : ! d(m*T1*T2*T3)/dt = del*(u*m*T1*T2*T3) (3)
95 : !
96 : ! where d is a partial derivative, del is the 2D divergence operator,
97 : ! u is the horizontal velocity, m is the mass density field, and
98 : ! T1, T2, and T3 are tracers.
99 : !
100 : ! In CICE, these equations have the form
101 : !
102 : ! da/dt = del*(u*a) (4)
103 : ! dv/dt = d(a*h)/dt = del*(u*a*h) (5)
104 : ! de/dt = d(a*h*q)/dt = del*(u*a*h*q) (6)
105 : ! d(aT)/dt = del*(u*a*t) (7)
106 : !
107 : ! where a = fractional ice area, v = ice/snow volume, h = v/a = thickness,
108 : ! e = ice/snow internal energy (J/m^2), q = e/v = internal energy per
109 : ! unit volume (J/m^3), and T is a tracer. These equations express
110 : ! conservation of ice area, volume, internal energy, and area-weighted
111 : ! tracer, respectively.
112 : !
113 : ! (Note: In CICE, a, v and e are prognostic quantities from which
114 : ! h and q are diagnosed. The remapping routine works with tracers,
115 : ! which means that h and q must be derived from a, v, and e before
116 : ! calling the remapping routine.)
117 : !
118 : ! Earlier versions of CICE assumed fixed ice and snow density.
119 : ! Beginning with CICE 4.0, the ice and snow density can be variable.
120 : ! In this case, equations (5) and (6) are replaced by
121 : !
122 : ! dv/dt = d(a*h)/dt = del*(u*a*h) (8)
123 : ! dm/dt = d(a*h*rho)/dt = del*(u*a*h*rho) (9)
124 : ! de/dt = d(a*h*rho*qm)/dt = del*(u*a*h*rho*qm) (10)
125 : !
126 : ! where rho = density and qm = internal energy per unit mass (J/kg).
127 : ! Eq. (9) expresses mass conservation, which in the variable-density
128 : ! case is no longer equivalent to volume conservation (8).
129 : !
130 : ! Tracers satisfying equations of the form (1) are called "type 1."
131 : ! In CICE the paradigmatic type 1 tracers are hi and hs.
132 : !
133 : ! Tracers satisfying equations of the form (2) are called "type 2".
134 : ! The paradigmatic type 2 tracers are qi and qs (or rhoi and rhos
135 : ! in the variable-density case).
136 : !
137 : ! Tracers satisfying equations of the form (3) are called "type 3."
138 : ! The paradigmatic type 3 tracers are qmi and qms in the variable-density
139 : ! case. There are no such tracers in the constant-density case.
140 : !
141 : ! The fields a, T1, and T2 are reconstructed in each grid cell with
142 : ! 2nd-order accuracy. T3 is reconstructed with 1st-order accuracy
143 : ! (i.e., it is transported in upwind fashion) in order to avoid
144 : ! additional mathematical complexity.
145 : !
146 : ! The mass-like field lives in the array "mm" (shorthand for mean
147 : ! mass) and the tracers fields in the array "tm" (mean tracers).
148 : ! In order to transport tracers correctly, the remapping routine
149 : ! needs to know the tracers types and relationships. This is done
150 : ! as follows:
151 : !
152 : ! Each field in the "tm" array is assigned an index, 1:ntrace.
153 : ! (Note: ntrace is not the same as ntrcr, the number of tracers
154 : ! in the trcrn state variable array. For remapping purposes we
155 : ! have additional tracers hi and hs.)
156 : !
157 : ! The tracer types (1,2,3) are contained in the "tracer_type" array.
158 : ! For standard CICE:
159 : !
160 : ! tracer_type = (1 1 1 2 2 2 2 2)
161 : !
162 : ! Type 2 and type 3 tracers are said to depend on type 1 tracers.
163 : ! For instance, qi depends on hi, which is to say that
164 : ! there is a conservation equation of the form (2) or (6).
165 : ! Thus we define a "depend" array. For standard CICE:
166 : !
167 : ! depend = (0 0 0 1 1 1 1 2)
168 : !
169 : ! which implies that elements 1-3 (hi, hs, Ts) are type 1,
170 : ! elements 4-7 (qi) depend on element 1 (hi), and element 8 (qs)
171 : ! depends on element 2 (hs).
172 : !
173 : ! We also define a logical array "has_dependents". In standard CICE:
174 : !
175 : ! has_dependents = (T T F F F F F F),
176 : !
177 : ! which means that only elements 1 and 2 (hi and hs) have dependent
178 : ! tracers.
179 : !
180 : ! For the variable-density case, things are a bit more complicated.
181 : ! Suppose we have 4 variable-density ice layers and one variable-
182 : ! density snow layer. Then the indexing is as follows:
183 : ! 1 = hi
184 : ! 2 = hs
185 : ! 3 = Ts
186 : ! 4-7 = rhoi
187 : ! 8 = rhos
188 : ! 9-12 = qmi
189 : ! 13 = qms
190 : !
191 : ! The key arrays are:
192 : !
193 : ! tracer_type = (1 1 1 2 2 2 2 2 3 3 3 3 3)
194 : !
195 : ! depend = (0 0 0 1 1 1 1 2 4 5 6 7 8)
196 : !
197 : ! has_dependents = (T T F T T T T T F F F F F)
198 : !
199 : ! which imply that hi and hs are type 1 with dependents rhoi and rhos,
200 : ! while rhoi and rhos are type 2 with dependents qmi and qms.
201 : !
202 : ! Tracers added to the ntrcr array are handled automatically
203 : ! by the remapping with little extra coding. It is necessary
204 : ! only to provide the correct type and dependency information.
205 : !
206 : ! When using this routine in other models, most of the tracer dependency
207 : ! apparatus may be irrelevant. In a layered ocean model, for example,
208 : ! the transported fields are the layer thickness h (the mass density
209 : ! field) and two or more tracers (T, S, and various trace species).
210 : ! Suppose there are just two tracers, T and S. Then the tracer arrays
211 : ! have the values:
212 : !
213 : ! tracer_type = (1 1)
214 : ! depend = (0 0)
215 : ! has_dependents = (F F)
216 : !
217 : ! which is to say that all tracer transport equations are of the form (1).
218 : !
219 : ! The tracer dependency arrays are optional input arguments for the
220 : ! main remapping subroutine. If these arrays are not passed in, they
221 : ! take on the default values tracer_type(:) = 1, depend(:) = 0, and
222 : ! has_dependents(:) = F, which are appropriate for most purposes.
223 : !
224 : ! Another optional argument is integral_order. If integral_order = 1,
225 : ! then the triangle integrals are exact for linear functions of x and y.
226 : ! If integral_order = 2, these integrals are exact for both linear and
227 : ! quadratic functions. If integral_order = 3, integrals are exact for
228 : ! cubic functions as well. If all tracers are of type 1, then the
229 : ! integrals of mass*tracer are quadratic, and integral_order = 2 is
230 : ! sufficient. In CICE, where there are type 2 tracers, we integrate
231 : ! functions of the form mass*tracer1*tracer2. Thus integral_order = 3
232 : ! is required for exactness, though integral_order = 2 may be good enough
233 : ! in practice.
234 : !
235 : ! Finally, a few words about the edgearea fields:
236 : !
237 : ! In earlier versions of this scheme, the divergence of the velocity
238 : ! field implied by the remapping was, in general, different from the
239 : ! value of del*u computed in the dynamics. For energetic consistency
240 : ! (in CICE as well as in layered ocean models such as HYPOP),
241 : ! these two values should agree. This can be ensured by setting
242 : ! l_fixed_area = T and specifying the area transported across each grid
243 : ! cell edge in the arrays edgearea_e and edgearea_n. The departure
244 : ! regions are then tweaked, following an idea by Mats Bentsen, such
245 : ! that they have the desired area. If l_fixed_area = F, these regions
246 : ! are not tweaked, and the edgearea arrays are output variables.
247 : !
248 : !=======================================================================
249 :
250 : contains
251 :
252 : !=======================================================================
253 : !
254 : ! Grid quantities used by the remapping transport scheme
255 : !
256 : ! Note: the arrays xyav, xxxav, etc are not needed for rectangular grids
257 : ! but may be needed in the future for other nonuniform grids. They have
258 : ! been commented out here to save memory and flops.
259 : !
260 : ! author William H. Lipscomb, LANL
261 :
262 37 : subroutine init_remap
263 :
264 : use ice_domain, only: nblocks
265 : use ice_grid, only: xav, yav, xxav, yyav
266 : ! dxT, dyT, xyav, &
267 : ! xxxav, xxyav, xyyav, yyyav
268 :
269 : integer (kind=int_kind) :: &
270 : i, j, iblk ! standard indices
271 :
272 : character(len=*), parameter :: subname = '(init_remap)'
273 :
274 : ! Compute grid cell average geometric quantities on the scaled
275 : ! rectangular grid with dx = 1, dy = 1.
276 : !
277 : ! Note: On a rectangular grid, the integral of any odd function
278 : ! of x or y = 0.
279 :
280 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
281 50 : do iblk = 1, nblocks
282 1276 : do j = 1, ny_block
283 50267 : do i = 1, nx_block
284 49028 : xav(i,j,iblk) = c0
285 49028 : yav(i,j,iblk) = c0
286 : !!! These formulas would be used on a rectangular grid
287 : !!! with dimensions (dxT, dyT):
288 : !!! xxav(i,j,iblk) = dxT(i,j,iblk)**2 / c12
289 : !!! yyav(i,j,iblk) = dyT(i,j,iblk)**2 / c12
290 49028 : xxav(i,j,iblk) = c1/c12
291 50234 : yyav(i,j,iblk) = c1/c12
292 : ! xyav(i,j,iblk) = c0
293 : ! xxxav(i,j,iblk) = c0
294 : ! xxyav(i,j,iblk) = c0
295 : ! xyyav(i,j,iblk) = c0
296 : ! yyyav(i,j,iblk) = c0
297 : enddo
298 : enddo
299 : enddo
300 : !$OMP END PARALLEL DO
301 :
302 : !-------------------------------------------------------------------
303 : ! Set logical l_fixed_area depending of the grid type.
304 : !
305 : ! If l_fixed_area is true, the area of each departure region is
306 : ! computed in advance (e.g., by taking the divergence of the
307 : ! velocity field and passed to locate_triangles. The departure
308 : ! regions are adjusted to obtain the desired area.
309 : ! If false, edgearea is computed in locate_triangles and passed out.
310 : !
311 : ! l_fixed_area = .false. has been the default approach in CICE. It is
312 : ! used like this for the B-grid. However, idealized tests with the
313 : ! C-grid have shown that l_fixed_area = .false. leads to a checkerboard
314 : ! pattern in prognostic fields (e.g. aice). Using l_fixed_area = .true.
315 : ! eliminates the checkerboard pattern in C-grid simulations.
316 : !
317 : !-------------------------------------------------------------------
318 :
319 37 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
320 0 : l_fixed_area = .true.
321 : else
322 37 : l_fixed_area = .false.
323 : endif
324 :
325 37 : end subroutine init_remap
326 :
327 : !=======================================================================
328 : !
329 : ! Solve the transport equations for one timestep using the incremental
330 : ! remapping scheme developed by John Dukowicz and John Baumgardner (DB)
331 : ! and modified for sea ice by William Lipscomb and Elizabeth Hunke.
332 : !
333 : ! This scheme preserves monotonicity of ice area and tracers. That is,
334 : ! it does not produce new extrema. It is second-order accurate in space,
335 : ! except where gradients are limited to preserve monotonicity.
336 : !
337 : ! This version of the remapping allows the user to specify the areal
338 : ! flux across each edge, based on an idea developed by Mats Bentsen.
339 : !
340 : ! author William H. Lipscomb, LANL
341 : ! 2006: Moved driver (subroutine transport_remap) into separate module.
342 : ! Geometry changes (logically rectangular coordinates, fixed
343 : ! area fluxes)
344 :
345 5784 : subroutine horizontal_remap (dt, ntrace, &
346 : uvel, vvel, & ! LCOV_EXCL_LINE
347 : mm, tm, & ! LCOV_EXCL_LINE
348 : tracer_type, depend, & ! LCOV_EXCL_LINE
349 : has_dependents, & ! LCOV_EXCL_LINE
350 : integral_order, & ! LCOV_EXCL_LINE
351 : l_dp_midpt, & ! LCOV_EXCL_LINE
352 0 : uvelE, vvelN)
353 :
354 : use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
355 : ice_HaloDestroy
356 : use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_remap
357 : use ice_blocks, only: block, get_block, nghost
358 : use ice_grid, only: HTE, HTN, dxu, dyu, &
359 : earea, narea, tarear, hm, & ! LCOV_EXCL_LINE
360 : xav, yav, xxav, yyav
361 : ! xyav, xxxav, xxyav, xyyav, yyyav
362 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
363 :
364 : real (kind=dbl_kind), intent(in) :: &
365 : dt ! time step
366 :
367 : integer (kind=int_kind), intent(in) :: &
368 : ntrace ! number of tracers in use
369 :
370 : real (kind=dbl_kind), intent(in), dimension(nx_block,ny_block,max_blocks) :: &
371 : uvel, & ! x-component of velocity (m/s) ugrid ! LCOV_EXCL_LINE
372 : vvel ! y-component of velocity (m/s) ugrid
373 :
374 : real (kind=dbl_kind), intent(in), optional, dimension(nx_block,ny_block,max_blocks) :: &
375 : uvelE, & ! x-component of velocity (m/s) egrid ! LCOV_EXCL_LINE
376 : vvelN ! y-component of velocity (m/s) ngrid
377 :
378 : real (kind=dbl_kind), intent(inout), dimension (nx_block,ny_block,0:ncat,max_blocks) :: &
379 : mm ! mean mass values in each grid cell
380 :
381 : real (kind=dbl_kind), intent(inout), dimension (nx_block,ny_block,ntrace,ncat,max_blocks) :: &
382 : tm ! mean tracer values in each grid cell
383 :
384 : integer (kind=int_kind), dimension (ntrace), intent(in) :: &
385 : tracer_type , & ! = 1, 2, or 3 (see comments above) ! LCOV_EXCL_LINE
386 : depend ! tracer dependencies (see above)
387 :
388 : logical (kind=log_kind), dimension (ntrace), intent(in) :: &
389 : has_dependents ! true if a tracer has dependent tracers
390 :
391 : integer (kind=int_kind), intent(in) :: &
392 : integral_order ! polynomial order for triangle integrals
393 :
394 : logical (kind=log_kind), intent(in) :: &
395 : l_dp_midpt ! if true, find departure points using
396 : ! corrected midpoint velocity
397 : ! local variables
398 :
399 : integer (kind=int_kind) :: &
400 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
401 : iblk , & ! block index ! LCOV_EXCL_LINE
402 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
403 : n, m ! ice category, tracer indices
404 :
405 : integer (kind=int_kind), dimension(0:ncat,max_blocks) :: &
406 11568 : icellsnc ! number of cells with ice
407 :
408 : integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat) :: &
409 11568 : indxinc, indxjnc ! compressed i/j indices
410 :
411 : real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
412 : edgearea_e , & ! area of departure regions for east edges ! LCOV_EXCL_LINE
413 1260048 : edgearea_n ! area of departure regions for north edges
414 :
415 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
416 : dpx , & ! x coordinates of departure points at cell corners ! LCOV_EXCL_LINE
417 5015568 : dpy ! y coordinates of departure points at cell corners
418 :
419 : real (kind=dbl_kind), dimension(nx_block,ny_block,0:ncat,max_blocks) :: &
420 : mc , & ! mass at geometric center of cell ! LCOV_EXCL_LINE
421 60086928 : mx, my ! limited derivative of mass wrt x and y
422 :
423 : real (kind=dbl_kind), dimension(nx_block,ny_block,0:ncat) :: &
424 7518288 : mmask ! = 1. if mass is present, = 0. otherwise
425 :
426 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat,max_blocks) :: &
427 : tc , & ! tracer values at geometric center of cell ! LCOV_EXCL_LINE
428 1301493648 : tx, ty ! limited derivative of tracer wrt x and y
429 :
430 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat) :: &
431 162694128 : tmask ! = 1. if tracer is present, = 0. otherwise
432 :
433 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat) :: &
434 15026448 : mflxe, mflxn ! mass transports across E and N cell edges
435 :
436 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace,ncat) :: &
437 325378128 : mtflxe, mtflxn ! mass*tracer transports across E and N cell edges
438 :
439 : real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups) :: &
440 7518288 : triarea ! area of east-edge departure triangle
441 :
442 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups) :: &
443 60092688 : xp, yp ! x and y coordinates of special triangle points
444 : ! (need 4 points for triangle integrals)
445 : integer (kind=int_kind), dimension (nx_block,ny_block,ngroups) :: &
446 : iflux , & ! i index of cell contributing transport ! LCOV_EXCL_LINE
447 11568 : jflux ! j index of cell contributing transport
448 :
449 : integer (kind=int_kind), dimension(ngroups,max_blocks) :: &
450 11568 : icellsng ! number of cells with ice
451 :
452 : integer (kind=int_kind), dimension(nx_block*ny_block,ngroups) :: &
453 11568 : indxing, indxjng ! compressed i/j indices
454 :
455 : integer (kind=int_kind), dimension(nx_block,ny_block,max_blocks) :: &
456 10128 : halomask ! temporary mask for fast halo updates
457 :
458 : logical (kind=log_kind) :: &
459 : l_stop ! if true, abort the model
460 :
461 : integer (kind=int_kind) :: &
462 : istop, jstop ! indices of grid cell where model aborts
463 :
464 : character (len=char_len) :: &
465 : edge ! 'north' or 'east'
466 :
467 : type (ice_halo) :: &
468 : halo_info_tracer ! masked halo
469 :
470 : type (block) :: &
471 : this_block ! block information for current block
472 :
473 : character(len=*), parameter :: subname = '(horizontal_remap)'
474 :
475 : !-------------------------------------------------------------------
476 : ! Remap the ice area and associated tracers.
477 : ! Remap the open water area (without tracers).
478 : !-------------------------------------------------------------------
479 :
480 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,n, &
481 : !$OMP indxinc,indxjnc,mmask,tmask,istop,jstop,l_stop) & ! LCOV_EXCL_LINE
482 2880 : !$OMP SCHEDULE(runtime)
483 8688 : do iblk = 1, nblocks
484 :
485 5784 : l_stop = .false.
486 5784 : istop = 0
487 5784 : jstop = 0
488 :
489 5784 : this_block = get_block(blocks_ice(iblk),iblk)
490 5784 : ilo = this_block%ilo
491 5784 : ihi = this_block%ihi
492 5784 : jlo = this_block%jlo
493 5784 : jhi = this_block%jhi
494 :
495 : !-------------------------------------------------------------------
496 : ! Compute masks and count ice cells.
497 : ! Masks are used to prevent tracer values in cells without ice from
498 : ! being used to compute tracer gradients.
499 : !-------------------------------------------------------------------
500 :
501 : call make_masks (nx_block, ny_block, &
502 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
503 : nghost, ntrace, & ! LCOV_EXCL_LINE
504 : has_dependents, icellsnc(:,iblk), & ! LCOV_EXCL_LINE
505 : indxinc(:,:), indxjnc(:,:), & ! LCOV_EXCL_LINE
506 : mm (:,:,:,iblk), mmask(:,:,:), & ! LCOV_EXCL_LINE
507 5784 : tm(:,:,:,:,iblk), tmask(:,:,:,:))
508 :
509 : !-------------------------------------------------------------------
510 : ! Construct linear fields, limiting gradients to preserve monotonicity.
511 : ! Note: Pass in unit arrays instead of true distances HTE, HTN, etc.
512 : ! The resulting gradients are in scaled coordinates.
513 : !-------------------------------------------------------------------
514 :
515 : ! open water
516 : call construct_fields(nx_block, ny_block, &
517 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
518 : nghost, ntrace, & ! LCOV_EXCL_LINE
519 : tracer_type, depend, & ! LCOV_EXCL_LINE
520 : has_dependents, icellsnc(0,iblk), & ! LCOV_EXCL_LINE
521 : indxinc(:,0), indxjnc(:,0), & ! LCOV_EXCL_LINE
522 : hm (:,:,iblk), xav (:,:,iblk), & ! LCOV_EXCL_LINE
523 : yav (:,:,iblk), xxav (:,:,iblk), & ! LCOV_EXCL_LINE
524 : yyav (:,:,iblk), & ! LCOV_EXCL_LINE
525 : ! xyav (:,:,iblk), & ! LCOV_EXCL_LINE
526 : ! xxxav (:,:,iblk), xxyav (:,:,iblk), & ! LCOV_EXCL_LINE
527 : ! xyyav (:,:,iblk), yyyav (:,:,iblk), & ! LCOV_EXCL_LINE
528 : mm (:,:,0,iblk), mc (:,:,0,iblk), & ! LCOV_EXCL_LINE
529 : mx (:,:,0,iblk), my (:,:,0,iblk), & ! LCOV_EXCL_LINE
530 5784 : mmask(:,:,0) )
531 :
532 : ! ice categories
533 :
534 34704 : do n = 1, ncat
535 :
536 : call construct_fields(nx_block, ny_block, &
537 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
538 : nghost, ntrace, & ! LCOV_EXCL_LINE
539 : tracer_type, depend, & ! LCOV_EXCL_LINE
540 : has_dependents, icellsnc (n,iblk), & ! LCOV_EXCL_LINE
541 : indxinc (:,n), indxjnc(:,n), & ! LCOV_EXCL_LINE
542 : hm (:,:,iblk), xav (:,:,iblk), & ! LCOV_EXCL_LINE
543 : yav (:,:,iblk), xxav (:,:,iblk), & ! LCOV_EXCL_LINE
544 : yyav (:,:,iblk), & ! LCOV_EXCL_LINE
545 : ! xyav (:,:,iblk), & ! LCOV_EXCL_LINE
546 : ! xxxav (:,:,iblk), xxyav (:,:,iblk), & ! LCOV_EXCL_LINE
547 : ! xyyav (:,:,iblk), yyyav (:,:,iblk), & ! LCOV_EXCL_LINE
548 : mm (:,:,n,iblk), mc (:,:,n,iblk), & ! LCOV_EXCL_LINE
549 : mx (:,:,n,iblk), my (:,:,n,iblk), & ! LCOV_EXCL_LINE
550 : mmask (:,:,n), & ! LCOV_EXCL_LINE
551 : tm (:,:,:,n,iblk), tc(:,:,:,n,iblk), & ! LCOV_EXCL_LINE
552 : tx (:,:,:,n,iblk), ty(:,:,:,n,iblk), & ! LCOV_EXCL_LINE
553 34704 : tmask(:,:,:,n) )
554 :
555 : enddo ! n
556 :
557 : !-------------------------------------------------------------------
558 : ! Given velocity field at cell corners, compute departure points
559 : ! of trajectories.
560 : !-------------------------------------------------------------------
561 :
562 : call departure_points(nx_block, ny_block, &
563 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
564 : nghost, dt, & ! LCOV_EXCL_LINE
565 : uvel (:,:,iblk), vvel(:,:,iblk), & ! LCOV_EXCL_LINE
566 : dxu (:,:,iblk), dyu (:,:,iblk), & ! LCOV_EXCL_LINE
567 : HTN (:,:,iblk), HTE (:,:,iblk), & ! LCOV_EXCL_LINE
568 : dpx (:,:,iblk), dpy (:,:,iblk), & ! LCOV_EXCL_LINE
569 : l_dp_midpt, l_stop, & ! LCOV_EXCL_LINE
570 5784 : istop, jstop)
571 :
572 11568 : if (l_stop) then
573 0 : call diagnostic_abort(istop,jstop,iblk,'bad departure points')
574 : endif
575 :
576 : enddo ! iblk
577 : !$OMP END PARALLEL DO
578 :
579 : !-------------------------------------------------------------------
580 : ! Ghost cell updates
581 : ! If nghost >= 2, these calls are not needed
582 : !-------------------------------------------------------------------
583 :
584 : if (nghost==1) then
585 :
586 5784 : call ice_timer_start(timer_bound)
587 :
588 : ! departure points
589 0 : call ice_HaloUpdate (dpx, halo_info, &
590 5784 : field_loc_NEcorner, field_type_vector)
591 0 : call ice_HaloUpdate (dpy, halo_info, &
592 5784 : field_loc_NEcorner, field_type_vector)
593 :
594 : ! mass field
595 0 : call ice_HaloUpdate (mc, halo_info, &
596 5784 : field_loc_center, field_type_scalar)
597 0 : call ice_HaloUpdate (mx, halo_info, &
598 5784 : field_loc_center, field_type_vector)
599 0 : call ice_HaloUpdate (my, halo_info, &
600 5784 : field_loc_center, field_type_vector)
601 :
602 : ! tracer fields
603 5784 : if (maskhalo_remap) then
604 0 : halomask(:,:,:) = 0
605 0 : !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,n,m,j,i) SCHEDULE(runtime)
606 0 : do iblk = 1, nblocks
607 0 : this_block = get_block(blocks_ice(iblk),iblk)
608 0 : ilo = this_block%ilo
609 0 : ihi = this_block%ihi
610 0 : jlo = this_block%jlo
611 0 : jhi = this_block%jhi
612 :
613 0 : do n = 1, ncat
614 0 : do m = 1, ntrace
615 0 : do j = jlo, jhi
616 0 : do i = ilo, ihi
617 0 : if (tc(i,j,m,n,iblk) /= c0) halomask(i,j,iblk) = 1
618 0 : if (tx(i,j,m,n,iblk) /= c0) halomask(i,j,iblk) = 1
619 0 : if (ty(i,j,m,n,iblk) /= c0) halomask(i,j,iblk) = 1
620 : enddo
621 : enddo
622 : enddo
623 : enddo
624 : enddo
625 : !$OMP END PARALLEL DO
626 0 : call ice_HaloUpdate (halomask, halo_info, &
627 0 : field_loc_center, field_type_scalar)
628 0 : call ice_HaloMask(halo_info_tracer, halo_info, halomask)
629 :
630 0 : call ice_HaloUpdate (tc, halo_info_tracer, &
631 0 : field_loc_center, field_type_scalar)
632 0 : call ice_HaloUpdate (tx, halo_info_tracer, &
633 0 : field_loc_center, field_type_vector)
634 0 : call ice_HaloUpdate (ty, halo_info_tracer, &
635 0 : field_loc_center, field_type_vector)
636 0 : call ice_HaloDestroy(halo_info_tracer)
637 : else
638 0 : call ice_HaloUpdate (tc, halo_info, &
639 5784 : field_loc_center, field_type_scalar)
640 0 : call ice_HaloUpdate (tx, halo_info, &
641 5784 : field_loc_center, field_type_vector)
642 0 : call ice_HaloUpdate (ty, halo_info, &
643 5784 : field_loc_center, field_type_vector)
644 : endif
645 5784 : call ice_timer_stop(timer_bound)
646 :
647 : endif ! nghost
648 :
649 : ! tcraig, this OMP loop sometimes fails with cce/14.0.3, compiler bug??
650 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,n, &
651 : !$OMP edgearea_e,edgearea_n,edge,iflux,jflux, & ! LCOV_EXCL_LINE
652 : !$OMP xp,yp,indxing,indxjng,mflxe,mflxn, & ! LCOV_EXCL_LINE
653 : !$OMP mtflxe,mtflxn,triarea,istop,jstop,l_stop) & ! LCOV_EXCL_LINE
654 2880 : !$OMP SCHEDULE(runtime)
655 8688 : do iblk = 1, nblocks
656 :
657 5784 : l_stop = .false.
658 5784 : istop = 0
659 5784 : jstop = 0
660 :
661 5784 : this_block = get_block(blocks_ice(iblk),iblk)
662 5784 : ilo = this_block%ilo
663 5784 : ihi = this_block%ihi
664 5784 : jlo = this_block%jlo
665 5784 : jhi = this_block%jhi
666 :
667 : !-------------------------------------------------------------------
668 : ! If l_fixed_area is true, compute edgearea by taking the divergence
669 : ! of the velocity field. Otherwise, initialize edgearea.
670 : !-------------------------------------------------------------------
671 :
672 204456 : do j = 1, ny_block
673 7151880 : do i = 1, nx_block
674 6947424 : edgearea_e(i,j) = c0
675 7146096 : edgearea_n(i,j) = c0
676 : enddo
677 : enddo
678 :
679 5784 : if (l_fixed_area) then
680 0 : if (grid_ice == 'CD' .or. grid_ice == 'C') then ! velocities are already on the center
681 0 : if (.not.present(uvelE).or..not.present(vvelN)) then
682 0 : call abort_ice (subname//'ERROR: uvelE,vvelN required with C|CD and l_fixed_area')
683 : endif
684 :
685 0 : do j = jlo, jhi
686 0 : do i = ilo-1, ihi
687 0 : edgearea_e(i,j) = uvelE(i,j,iblk) * HTE(i,j,iblk) * dt
688 : enddo
689 : enddo
690 :
691 0 : do j = jlo-1, jhi
692 0 : do i = ilo, ihi
693 0 : edgearea_n(i,j) = vvelN(i,j,iblk)*HTN(i,j,iblk) * dt
694 : enddo
695 : enddo
696 :
697 : else
698 0 : do j = jlo, jhi
699 0 : do i = ilo-1, ihi
700 : edgearea_e(i,j) = (uvel(i,j,iblk) + uvel(i,j-1,iblk)) &
701 0 : * p5 * HTE(i,j,iblk) * dt
702 : enddo
703 : enddo
704 :
705 0 : do j = jlo-1, jhi
706 0 : do i = ilo, ihi
707 : edgearea_n(i,j) = (vvel(i,j,iblk) + vvel(i-1,j,iblk)) &
708 0 : * p5 * HTN(i,j,iblk) * dt
709 : enddo
710 : enddo
711 : endif
712 : endif
713 :
714 : !-------------------------------------------------------------------
715 : ! Transports for east cell edges.
716 : !-------------------------------------------------------------------
717 :
718 : !-------------------------------------------------------------------
719 : ! Compute areas and vertices of departure triangles.
720 : !-------------------------------------------------------------------
721 :
722 5784 : edge = 'east'
723 : call locate_triangles(nx_block, ny_block, &
724 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
725 : nghost, edge, & ! LCOV_EXCL_LINE
726 : icellsng(:,iblk), & ! LCOV_EXCL_LINE
727 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
728 : dpx (:,:,iblk), dpy(:,:,iblk), & ! LCOV_EXCL_LINE
729 : dxu (:,:,iblk), dyu(:,:,iblk), & ! LCOV_EXCL_LINE
730 : earea (:,:,iblk), narea (:,:,iblk), & ! LCOV_EXCL_LINE
731 : xp (:,:,:,:), yp (:,:,:,:), & ! LCOV_EXCL_LINE
732 : iflux, jflux, & ! LCOV_EXCL_LINE
733 5784 : triarea, edgearea_e(:,:))
734 :
735 : !-------------------------------------------------------------------
736 : ! Given triangle vertices, compute coordinates of triangle points
737 : ! needed for transport integrals.
738 : !-------------------------------------------------------------------
739 :
740 : call triangle_coordinates (nx_block, ny_block, &
741 : integral_order, icellsng(:,iblk), & ! LCOV_EXCL_LINE
742 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
743 5784 : xp, yp)
744 :
745 : !-------------------------------------------------------------------
746 : ! Compute the transport across east cell edges by summing contributions
747 : ! from each triangle.
748 : !-------------------------------------------------------------------
749 :
750 : ! open water
751 : call transport_integrals(nx_block, ny_block, &
752 : ntrace, icellsng (:,iblk), & ! LCOV_EXCL_LINE
753 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
754 : tracer_type, depend, & ! LCOV_EXCL_LINE
755 : integral_order, triarea, & ! LCOV_EXCL_LINE
756 : iflux, jflux, & ! LCOV_EXCL_LINE
757 : xp, yp, & ! LCOV_EXCL_LINE
758 : mc(:,:,0,iblk), mx (:,:,0,iblk), & ! LCOV_EXCL_LINE
759 5784 : my(:,:,0,iblk), mflxe(:,:,0))
760 :
761 : ! ice categories
762 34704 : do n = 1, ncat
763 : call transport_integrals &
764 : (nx_block, ny_block, & ! LCOV_EXCL_LINE
765 : ntrace, icellsng (:,iblk), & ! LCOV_EXCL_LINE
766 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
767 : tracer_type, depend, & ! LCOV_EXCL_LINE
768 : integral_order, triarea, & ! LCOV_EXCL_LINE
769 : iflux, jflux, & ! LCOV_EXCL_LINE
770 : xp, yp, & ! LCOV_EXCL_LINE
771 : mc(:,:, n,iblk), mx (:,:, n,iblk), & ! LCOV_EXCL_LINE
772 : my(:,:, n,iblk), mflxe (:,:, n), & ! LCOV_EXCL_LINE
773 : tc(:,:,:,n,iblk), tx (:,:,:,n,iblk), & ! LCOV_EXCL_LINE
774 34704 : ty(:,:,:,n,iblk), mtflxe(:,:,:,n))
775 :
776 : enddo
777 :
778 : !-------------------------------------------------------------------
779 : ! Repeat for north edges
780 : !-------------------------------------------------------------------
781 :
782 5784 : edge = 'north'
783 : call locate_triangles(nx_block, ny_block, &
784 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
785 : nghost, edge, & ! LCOV_EXCL_LINE
786 : icellsng(:,iblk), & ! LCOV_EXCL_LINE
787 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
788 : dpx (:,:,iblk), dpy (:,:,iblk), & ! LCOV_EXCL_LINE
789 : dxu (:,:,iblk), dyu (:,:,iblk), & ! LCOV_EXCL_LINE
790 : earea (:,:,iblk), narea (:,:,iblk), & ! LCOV_EXCL_LINE
791 : xp (:,:,:,:), yp(:,:,:,:), & ! LCOV_EXCL_LINE
792 : iflux, jflux, & ! LCOV_EXCL_LINE
793 5784 : triarea, edgearea_n(:,:))
794 :
795 : call triangle_coordinates (nx_block, ny_block, &
796 : integral_order, icellsng(:,iblk), & ! LCOV_EXCL_LINE
797 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
798 5784 : xp, yp)
799 :
800 : ! open water
801 : call transport_integrals(nx_block, ny_block, &
802 : ntrace, icellsng (:,iblk), & ! LCOV_EXCL_LINE
803 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
804 : tracer_type, depend, & ! LCOV_EXCL_LINE
805 : integral_order, triarea, & ! LCOV_EXCL_LINE
806 : iflux, jflux, & ! LCOV_EXCL_LINE
807 : xp, yp, & ! LCOV_EXCL_LINE
808 : mc(:,:,0,iblk), mx (:,:,0,iblk), & ! LCOV_EXCL_LINE
809 5784 : my(:,:,0,iblk), mflxn(:,:,0))
810 :
811 : ! ice categories
812 34704 : do n = 1, ncat
813 : call transport_integrals &
814 : (nx_block, ny_block, & ! LCOV_EXCL_LINE
815 : ntrace, icellsng (:,iblk), & ! LCOV_EXCL_LINE
816 : indxing(:,:), indxjng(:,:), & ! LCOV_EXCL_LINE
817 : tracer_type, depend, & ! LCOV_EXCL_LINE
818 : integral_order, triarea, & ! LCOV_EXCL_LINE
819 : iflux, jflux, & ! LCOV_EXCL_LINE
820 : xp, yp, & ! LCOV_EXCL_LINE
821 : mc(:,:, n,iblk), mx (:,:, n,iblk), & ! LCOV_EXCL_LINE
822 : my(:,:, n,iblk), mflxn (:,:, n), & ! LCOV_EXCL_LINE
823 : tc(:,:,:,n,iblk), tx (:,:,:,n,iblk), & ! LCOV_EXCL_LINE
824 34704 : ty(:,:,:,n,iblk), mtflxn(:,:,:,n))
825 :
826 : enddo ! n
827 :
828 : !-------------------------------------------------------------------
829 : ! Update the ice area and tracers.
830 : !-------------------------------------------------------------------
831 :
832 : ! open water
833 : call update_fields (nx_block, ny_block, &
834 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
835 : ntrace, & ! LCOV_EXCL_LINE
836 : tracer_type, depend, & ! LCOV_EXCL_LINE
837 : tarear(:,:,iblk), l_stop, & ! LCOV_EXCL_LINE
838 : istop, jstop, & ! LCOV_EXCL_LINE
839 : mflxe(:,:,0), mflxn(:,:,0), & ! LCOV_EXCL_LINE
840 5784 : mm (:,:,0,iblk))
841 :
842 5784 : if (l_stop) then
843 0 : call diagnostic_abort(istop,jstop,iblk,'negative area (open water)')
844 : endif
845 :
846 : ! ice categories
847 40488 : do n = 1, ncat
848 :
849 : call update_fields(nx_block, ny_block, &
850 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
851 : ntrace, & ! LCOV_EXCL_LINE
852 : tracer_type, depend, & ! LCOV_EXCL_LINE
853 : tarear (:,:,iblk), l_stop, & ! LCOV_EXCL_LINE
854 : istop, jstop, & ! LCOV_EXCL_LINE
855 : mflxe (:,:, n), mflxn (:,:, n), & ! LCOV_EXCL_LINE
856 : mm (:,:, n,iblk), & ! LCOV_EXCL_LINE
857 : mtflxe(:,:,:,n), mtflxn(:,:,:,n), & ! LCOV_EXCL_LINE
858 28920 : tm (:,:,:,n,iblk))
859 :
860 34704 : if (l_stop) then
861 0 : write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
862 0 : istep1, my_task, iblk, n
863 0 : call diagnostic_abort(istop,jstop,iblk,'negative area (ice)')
864 : endif
865 : enddo ! n
866 :
867 : enddo ! iblk
868 : !$OMP END PARALLEL DO
869 :
870 5784 : end subroutine horizontal_remap
871 :
872 : !=======================================================================
873 : !
874 : ! Make area and tracer masks.
875 : !
876 : ! If an area is masked out (mm < puny), then the values of tracers
877 : ! in that grid cell are assumed to have no physical meaning.
878 : !
879 : ! Similarly, if a tracer with dependents is masked out
880 : ! (abs(tm) < puny), then the values of its dependent tracers in that
881 : ! grid cell are assumed to have no physical meaning.
882 : ! For example, the enthalpy value has no meaning if the thickness
883 : ! is zero.
884 : !
885 : ! author William H. Lipscomb, LANL
886 :
887 22944 : subroutine make_masks (nx_block, ny_block, &
888 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
889 : nghost, ntrace, & ! LCOV_EXCL_LINE
890 : has_dependents, & ! LCOV_EXCL_LINE
891 : icells, & ! LCOV_EXCL_LINE
892 : indxi, indxj, & ! LCOV_EXCL_LINE
893 : mm, mmask, & ! LCOV_EXCL_LINE
894 45888 : tm, tmask)
895 :
896 : integer (kind=int_kind), intent(in) :: &
897 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
898 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
899 : nghost , & ! number of ghost cells ! LCOV_EXCL_LINE
900 : ntrace ! number of tracers in use
901 :
902 :
903 : logical (kind=log_kind), dimension (ntrace), intent(in) :: &
904 : has_dependents ! true if a tracer has dependent tracers
905 :
906 : integer (kind=int_kind), dimension(0:ncat), intent(out) :: &
907 : icells ! number of cells with ice
908 :
909 : integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat), intent(out) :: &
910 : indxi , & ! compressed i/j indices ! LCOV_EXCL_LINE
911 : indxj
912 :
913 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), intent(in) :: &
914 : mm ! mean ice area in each grid cell
915 :
916 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), intent(out) :: &
917 : mmask ! = 1. if ice is present, else = 0.
918 :
919 : real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace, ncat), intent(in), optional :: &
920 : tm ! mean tracer values in each grid cell
921 :
922 : real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace, ncat), intent(out), optional :: &
923 : tmask ! = 1. if tracer is present, else = 0.
924 :
925 : ! local variables
926 :
927 : integer (kind=int_kind) :: &
928 : i, j, ij , & ! horizontal indices ! LCOV_EXCL_LINE
929 : n , & ! ice category index ! LCOV_EXCL_LINE
930 : nt ! tracer index
931 :
932 : real (kind=dbl_kind) :: &
933 5760 : puny !
934 :
935 : character(len=*), parameter :: subname = '(make_masks)'
936 :
937 22944 : call icepack_query_parameters(puny_out=puny)
938 22944 : call icepack_warnings_flush(nu_diag)
939 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
940 0 : file=__FILE__, line=__LINE__)
941 :
942 160608 : do n = 0, ncat
943 92757072 : do ij = 1, nx_block*ny_block
944 92596464 : indxi(ij,n) = 0
945 92734128 : indxj(ij,n) = 0
946 : enddo
947 : enddo
948 :
949 : !-------------------------------------------------------------------
950 : ! open water mask
951 : !-------------------------------------------------------------------
952 :
953 22944 : icells(0) = 0
954 753576 : do j = 1, ny_block
955 16186320 : do i = 1, nx_block
956 16163376 : if (mm(i,j,0) > puny) then
957 8718039 : mmask(i,j,0) = c1
958 8718039 : icells(0) = icells(0) + 1
959 8718039 : ij = icells(0)
960 8718039 : indxi(ij,0) = i
961 8718039 : indxj(ij,0) = j
962 : else
963 6714705 : mmask(i,j,0) = c0
964 : endif
965 : enddo
966 : enddo
967 :
968 137664 : do n = 1, ncat
969 :
970 : !-------------------------------------------------------------------
971 : ! Find grid cells where ice is present.
972 : !-------------------------------------------------------------------
973 :
974 114720 : icells(n) = 0
975 3767880 : do j = 1, ny_block
976 80931600 : do i = 1, nx_block
977 80816880 : if (mm(i,j,n) > puny) then
978 37655263 : icells(n) = icells(n) + 1
979 37655263 : ij = icells(n)
980 37655263 : indxi(ij,n) = i
981 37655263 : indxj(ij,n) = j
982 : endif ! mm > puny
983 : enddo
984 : enddo
985 :
986 : !-------------------------------------------------------------------
987 : ! ice area mask
988 : !-------------------------------------------------------------------
989 :
990 80931600 : mmask(:,:,n) = c0
991 37769983 : do ij = 1, icells(n)
992 37655263 : i = indxi(ij,n)
993 37655263 : j = indxj(ij,n)
994 37769983 : mmask(i,j,n) = c1
995 : enddo
996 :
997 : !-------------------------------------------------------------------
998 : ! tracer masks
999 : !-------------------------------------------------------------------
1000 :
1001 114720 : if (present(tm)) then
1002 :
1003 2104336320 : tmask(:,:,:,n) = c0
1004 :
1005 3097440 : do nt = 1, ntrace
1006 3097440 : if (has_dependents(nt)) then
1007 151079932 : do ij = 1, icells(n)
1008 150621052 : i = indxi(ij,n)
1009 150621052 : j = indxj(ij,n)
1010 151079932 : if (abs(tm(i,j,nt,n)) > puny) then
1011 105244210 : tmask(i,j,nt,n) = c1
1012 : endif
1013 : enddo
1014 : endif
1015 : enddo
1016 :
1017 : endif ! present(tm)
1018 :
1019 : !-------------------------------------------------------------------
1020 : ! Redefine icells
1021 : ! For nghost = 1, exclude ghost cells
1022 : ! For nghost = 2, include one layer of ghost cells
1023 : !-------------------------------------------------------------------
1024 :
1025 114720 : icells(n) = 0
1026 3561384 : do j = jlo-nghost+1, jhi+nghost-1
1027 69134640 : do i = ilo-nghost+1, ihi+nghost-1
1028 69019920 : if (mm(i,j,n) > puny) then
1029 33882955 : icells(n) = icells(n) + 1
1030 33882955 : ij = icells(n)
1031 33882955 : indxi(ij,n) = i
1032 33882955 : indxj(ij,n) = j
1033 : endif ! mm > puny
1034 : enddo
1035 : enddo
1036 :
1037 : enddo ! n
1038 :
1039 68832 : end subroutine make_masks
1040 :
1041 : !=======================================================================
1042 : !
1043 : ! Construct fields of ice area and tracers.
1044 : !
1045 : ! authors William H. Lipscomb, LANL
1046 : ! John R. Baumgardner, LANL
1047 :
1048 137664 : subroutine construct_fields (nx_block, ny_block, &
1049 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1050 : nghost, ntrace, & ! LCOV_EXCL_LINE
1051 : tracer_type, depend, & ! LCOV_EXCL_LINE
1052 : has_dependents, icells, & ! LCOV_EXCL_LINE
1053 : indxi, indxj, & ! LCOV_EXCL_LINE
1054 : hm, xav, & ! LCOV_EXCL_LINE
1055 : yav, xxav, & ! LCOV_EXCL_LINE
1056 : yyav, & ! LCOV_EXCL_LINE
1057 : ! xyav, & ! LCOV_EXCL_LINE
1058 : ! xxxav, xxyav, & ! LCOV_EXCL_LINE
1059 : ! xyyav, yyyav, & ! LCOV_EXCL_LINE
1060 : mm, mc, & ! LCOV_EXCL_LINE
1061 : mx, my, & ! LCOV_EXCL_LINE
1062 : mmask, & ! LCOV_EXCL_LINE
1063 : tm, tc, & ! LCOV_EXCL_LINE
1064 : tx, ty, & ! LCOV_EXCL_LINE
1065 114720 : tmask)
1066 :
1067 : integer (kind=int_kind), intent(in) :: &
1068 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1069 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
1070 : nghost , & ! number of ghost cells ! LCOV_EXCL_LINE
1071 : ntrace , & ! number of tracers in use ! LCOV_EXCL_LINE
1072 : icells ! number of cells with mass
1073 :
1074 : integer (kind=int_kind), dimension (ntrace), intent(in) :: &
1075 : tracer_type , & ! = 1, 2, or 3 (see comments above) ! LCOV_EXCL_LINE
1076 : depend ! tracer dependencies (see above)
1077 :
1078 : logical (kind=log_kind), dimension (ntrace), intent(in) :: &
1079 : has_dependents ! true if a tracer has dependent tracers
1080 :
1081 : integer (kind=int_kind), dimension(nx_block*ny_block), intent(in) :: &
1082 : indxi , & ! compressed i/j indices ! LCOV_EXCL_LINE
1083 : indxj
1084 :
1085 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1086 : hm , & ! land/boundary mask, thickness (T-cell) ! LCOV_EXCL_LINE
1087 : xav, yav , & ! mean T-cell values of x, y ! LCOV_EXCL_LINE
1088 : xxav, yyav ! mean T-cell values of xx, yy
1089 : ! xyav, , & ! mean T-cell values of xy
1090 : ! xxxav,xxyav,xyyav,yyyav ! mean T-cell values of xxx, xxy, xyy, yyy
1091 :
1092 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1093 : mm , & ! mean value of mass field ! LCOV_EXCL_LINE
1094 : mmask ! = 1. if ice is present, = 0. otherwise
1095 :
1096 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace), intent(in), optional :: &
1097 : tm , & ! mean tracer ! LCOV_EXCL_LINE
1098 : tmask ! = 1. if tracer is present, = 0. otherwise
1099 :
1100 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1101 : mc , & ! mass value at geometric center of cell ! LCOV_EXCL_LINE
1102 : mx, my ! limited derivative of mass wrt x and y
1103 :
1104 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace), intent(out), optional :: &
1105 : tc , & ! tracer at geometric center of cell ! LCOV_EXCL_LINE
1106 : tx, ty ! limited derivative of tracer wrt x and y
1107 :
1108 : ! local variables
1109 :
1110 : integer (kind=int_kind) :: &
1111 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
1112 : nt, nt1 , & ! tracer indices ! LCOV_EXCL_LINE
1113 : ij ! combined i/j horizontal index
1114 :
1115 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1116 : mxav , & ! x coordinate of center of mass ! LCOV_EXCL_LINE
1117 30238848 : myav ! y coordinate of center of mass
1118 :
1119 : real (kind=dbl_kind), dimension (nx_block,ny_block,ntrace) :: &
1120 : mtxav , & ! x coordinate of center of mass*tracer ! LCOV_EXCL_LINE
1121 781089408 : mtyav ! y coordinate of center of mass*tracer
1122 :
1123 : real (kind=dbl_kind) :: &
1124 : puny, & ! LCOV_EXCL_LINE
1125 34560 : w1, w2, w3, w7 ! work variables
1126 :
1127 : character(len=*), parameter :: subname = '(construct_fields)'
1128 :
1129 : !-------------------------------------------------------------------
1130 : ! Compute field values at the geometric center of each grid cell,
1131 : ! and compute limited gradients in the x and y directions.
1132 : !
1133 : ! For second order accuracy, each state variable is approximated as
1134 : ! a field varying linearly over x and y within each cell. For each
1135 : ! category, the integrated value of m(x,y) over the cell must
1136 : ! equal mm(i,j,n)*tarea(i,j), where tarea(i,j) is the cell area.
1137 : ! Similarly, the integrated value of m(x,y)*t(x,y) must equal
1138 : ! the total mass*tracer, mm(i,j,n)*tm(i,j,n)*tarea(i,j).
1139 : !
1140 : ! These integral conditions are satisfied for linear fields if we
1141 : ! stipulate the following:
1142 : ! (1) The mean mass, mm, is equal to the mass at the cell centroid.
1143 : ! (2) The mean value tm1 of type 1 tracers is equal to the value
1144 : ! at the center of mass.
1145 : ! (3) The mean value tm2 of type 2 tracers is equal to the value
1146 : ! at the center of mass*tm1, where tm2 depends on tm1.
1147 : ! (See comments at the top of the module.)
1148 : !
1149 : ! We want to find the value of each state variable at a standard
1150 : ! reference point, which we choose to be the geometric center of
1151 : ! the cell. The geometric center is located at the intersection
1152 : ! of the line joining the midpoints of the north and south edges
1153 : ! with the line joining the midpoints of the east and west edges.
1154 : ! To find the value at the geometric center, we must know the
1155 : ! location of the cell centroid/center of mass, along with the
1156 : ! mean value and the gradients with respect to x and y.
1157 : !
1158 : ! The cell gradients are first computed from the difference between
1159 : ! values in the neighboring cells, then limited by requiring that
1160 : ! no new extrema are created within the cell.
1161 : !
1162 : ! For rectangular coordinates the centroid and the geometric
1163 : ! center coincide, which means that some of the equations in this
1164 : ! subroutine could be simplified. However, the full equations
1165 : ! are retained for generality.
1166 : !-------------------------------------------------------------------
1167 :
1168 : !-------------------------------------------------------------------
1169 : ! Initialize
1170 : !-------------------------------------------------------------------
1171 :
1172 137664 : call icepack_query_parameters(puny_out=puny)
1173 137664 : call icepack_warnings_flush(nu_diag)
1174 137664 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1175 0 : file=__FILE__, line=__LINE__)
1176 :
1177 4521456 : do j = 1, ny_block
1178 97117920 : do i = 1, nx_block
1179 92596464 : mc(i,j) = c0
1180 92596464 : mx(i,j) = c0
1181 92596464 : my(i,j) = c0
1182 92596464 : mxav(i,j) = c0
1183 96980256 : myav(i,j) = c0
1184 : enddo
1185 : enddo
1186 :
1187 137664 : if (present(tm)) then
1188 3097440 : do nt = 1, ntrace
1189 98079600 : do j = 1, ny_block
1190 2104221600 : do i = 1, nx_block
1191 2006256720 : tc(i,j,nt) = c0
1192 2006256720 : tx(i,j,nt) = c0
1193 2101238880 : ty(i,j,nt) = c0
1194 : enddo
1195 : enddo
1196 : enddo
1197 : endif
1198 :
1199 : ! limited gradient of mass field in each cell (except masked cells)
1200 : ! Note: The gradient is computed in scaled coordinates with
1201 : ! dxT = dyT = hte = htn = 1.
1202 :
1203 : call limited_gradient (nx_block, ny_block, &
1204 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1205 : nghost, & ! LCOV_EXCL_LINE
1206 : mm, hm, & ! LCOV_EXCL_LINE
1207 : xav, yav, & ! LCOV_EXCL_LINE
1208 137664 : mx, my)
1209 :
1210 42738658 : do ij = 1,icells ! ice is present
1211 42600994 : i = indxi(ij)
1212 42600994 : j = indxj(ij)
1213 :
1214 : ! mass field at geometric center
1215 : ! echmod: xav = yav = 0
1216 42738658 : mc(i,j) = mm(i,j)
1217 :
1218 : ! mc(i,j) = mm(i,j) - xav(i,j)*mx(i,j) &
1219 : ! - yav(i,j)*my(i,j)
1220 :
1221 : enddo ! ij
1222 :
1223 : ! tracers
1224 :
1225 137664 : if (present(tm)) then
1226 :
1227 33997675 : do ij = 1,icells ! cells with mass
1228 33882955 : i = indxi(ij)
1229 33882955 : j = indxj(ij)
1230 :
1231 : ! center of mass (mxav,myav) for each cell
1232 : ! echmod: xyav = 0
1233 3085762 : mxav(i,j) = (mx(i,j)*xxav(i,j) &
1234 33882955 : + mc(i,j)*xav (i,j)) / mm(i,j)
1235 3085762 : myav(i,j) = (my(i,j)*yyav(i,j) &
1236 33997675 : + mc(i,j)*yav(i,j)) / mm(i,j)
1237 :
1238 : ! mxav(i,j) = (mx(i,j)*xxav(i,j) &
1239 : ! + my(i,j)*xyav(i,j) & ! LCOV_EXCL_LINE
1240 : ! + mc(i,j)*xav (i,j)) / mm(i,j)
1241 : ! myav(i,j) = (mx(i,j)*xyav(i,j) &
1242 : ! + my(i,j)*yyav(i,j) & ! LCOV_EXCL_LINE
1243 : ! + mc(i,j)*yav(i,j)) / mm(i,j)
1244 : enddo
1245 :
1246 3097440 : do nt = 1, ntrace
1247 :
1248 3097440 : if (tracer_type(nt)==1) then ! independent of other tracers
1249 :
1250 : call limited_gradient(nx_block, ny_block, &
1251 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1252 : nghost, & ! LCOV_EXCL_LINE
1253 : tm(:,:,nt), mmask, & ! LCOV_EXCL_LINE
1254 : mxav, myav, & ! LCOV_EXCL_LINE
1255 688320 : tx(:,:,nt), ty(:,:,nt))
1256 :
1257 688320 : if (has_dependents(nt)) then ! need center of area*tracer
1258 :
1259 11303640 : do j = 1, ny_block
1260 242794800 : do i = 1, nx_block
1261 231491160 : mtxav(i,j,nt) = c0
1262 242450640 : mtyav(i,j,nt) = c0
1263 : enddo
1264 : enddo
1265 :
1266 101993025 : do ij = 1, icells ! Note: no tx or ty in ghost cells
1267 : ! (bound calls are later)
1268 101648865 : i = indxi(ij)
1269 101648865 : j = indxj(ij)
1270 :
1271 : ! tracer value at geometric center
1272 9257286 : tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
1273 101648865 : - ty(i,j,nt)*myav(i,j)
1274 :
1275 101993025 : if (tmask(i,j,nt) > puny) then
1276 :
1277 : ! center of area*tracer
1278 82794204 : w1 = mc(i,j)*tc(i,j,nt)
1279 9056294 : w2 = mc(i,j)*tx(i,j,nt) &
1280 82794204 : + mx(i,j)*tc(i,j,nt)
1281 9056294 : w3 = mc(i,j)*ty(i,j,nt) &
1282 82794204 : + my(i,j)*tc(i,j,nt)
1283 : ! w4 = mx(i,j)*tx(i,j,nt)
1284 : ! w5 = mx(i,j)*ty(i,j,nt) &
1285 : ! + my(i,j)*tx(i,j,nt)
1286 : ! w6 = my(i,j)*ty(i,j,nt)
1287 82794204 : w7 = c1 / (mm(i,j)*tm(i,j,nt))
1288 : ! echmod: grid arrays = 0
1289 9056294 : mtxav(i,j,nt) = (w1*xav (i,j) + w2*xxav (i,j)) &
1290 82794204 : * w7
1291 9056294 : mtyav(i,j,nt) = (w1*yav(i,j) + w3*yyav(i,j)) &
1292 82794204 : * w7
1293 :
1294 : ! mtxav(i,j,nt) = (w1*xav (i,j) + w2*xxav (i,j) &
1295 : ! + w3*xyav (i,j) + w4*xxxav(i,j) & ! LCOV_EXCL_LINE
1296 : ! + w5*xxyav(i,j) + w6*xyyav(i,j)) & ! LCOV_EXCL_LINE
1297 : ! * w7
1298 : ! mtyav(i,j,nt) = (w1*yav(i,j) + w2*xyav (i,j) &
1299 : ! + w3*yyav(i,j) + w4*xxyav(i,j) & ! LCOV_EXCL_LINE
1300 : ! + w5*xyyav(i,j) + w6*yyyav(i,j)) & ! LCOV_EXCL_LINE
1301 : ! * w7
1302 : endif ! tmask
1303 :
1304 : enddo ! ij
1305 :
1306 : else ! no dependents
1307 :
1308 101993025 : do ij = 1, icells ! mass is present
1309 101648865 : i = indxi(ij)
1310 101648865 : j = indxj(ij)
1311 :
1312 : ! tracer value at geometric center
1313 9257286 : tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
1314 101993025 : - ty(i,j,nt)*myav(i,j)
1315 : enddo ! ij
1316 :
1317 : endif ! has_dependents
1318 :
1319 2294400 : elseif (tracer_type(nt)==2) then ! tracer nt depends on nt1
1320 2064960 : nt1 = depend(nt)
1321 :
1322 : call limited_gradient(nx_block, ny_block, &
1323 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1324 : nghost, & ! LCOV_EXCL_LINE
1325 : tm (:,:,nt), tmask(:,:,nt1), & ! LCOV_EXCL_LINE
1326 : mtxav(:,:,nt1), mtyav(:,:,nt1), & ! LCOV_EXCL_LINE
1327 2064960 : tx (:,:,nt), ty (:,:,nt))
1328 :
1329 611958150 : do ij = 1, icells ! ice is present
1330 609893190 : i = indxi(ij)
1331 609893190 : j = indxj(ij)
1332 55543716 : tc(i,j,nt) = tm(i,j,nt) &
1333 : - tx(i,j,nt) * mtxav(i,j,nt1) & ! LCOV_EXCL_LINE
1334 611958150 : - ty(i,j,nt) * mtyav(i,j,nt1)
1335 : enddo ! ij
1336 :
1337 229440 : elseif (tracer_type(nt)==3) then ! upwind approx; gradient = 0
1338 :
1339 67995350 : do ij = 1, icells
1340 67765910 : i = indxi(ij)
1341 67765910 : j = indxj(ij)
1342 :
1343 67995350 : tc(i,j,nt) = tm(i,j,nt)
1344 : ! tx(i,j,nt) = c0 ! already initialized to 0.
1345 : ! ty(i,j,nt) = c0
1346 : enddo ! ij
1347 :
1348 : endif ! tracer_type
1349 :
1350 : enddo ! ntrace
1351 :
1352 : endif ! present (tm)
1353 :
1354 596544 : end subroutine construct_fields
1355 :
1356 : !=======================================================================
1357 : !
1358 : ! Compute a limited gradient of the scalar field phi in scaled coordinates.
1359 : ! "Limited" means that we do not create new extrema in phi. For
1360 : ! instance, field values at the cell corners can neither exceed the
1361 : ! maximum of phi(i,j) in the cell and its eight neighbors, nor fall
1362 : ! below the minimum.
1363 : !
1364 : ! authors William H. Lipscomb, LANL
1365 : ! John R. Baumgardner, LANL
1366 :
1367 2890944 : subroutine limited_gradient (nx_block, ny_block, &
1368 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1369 : nghost, & ! LCOV_EXCL_LINE
1370 : phi, phimask, & ! LCOV_EXCL_LINE
1371 : cnx, cny, & ! LCOV_EXCL_LINE
1372 2890944 : gx, gy)
1373 :
1374 : integer (kind=int_kind), intent(in) :: &
1375 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1376 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
1377 : nghost ! number of ghost cells
1378 :
1379 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent (in) :: &
1380 : phi , & ! input tracer field (mean values in each grid cell) ! LCOV_EXCL_LINE
1381 : cnx , & ! x-coordinate of phi relative to geometric center of cell ! LCOV_EXCL_LINE
1382 : cny , & ! y-coordinate of phi relative to geometric center of cell ! LCOV_EXCL_LINE
1383 : phimask ! phimask(i,j) = 1 if phi(i,j) has physical meaning, = 0 otherwise.
1384 : ! For instance, aice has no physical meaning in land cells,
1385 : ! and hice no physical meaning where aice = 0.
1386 :
1387 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1388 : gx , & ! limited x-direction gradient ! LCOV_EXCL_LINE
1389 : gy ! limited y-direction gradient
1390 :
1391 : ! local variables
1392 :
1393 : integer (kind=int_kind) :: &
1394 : i, j, ij , & ! standard indices ! LCOV_EXCL_LINE
1395 : icells ! number of cells to limit
1396 :
1397 : integer (kind=int_kind), dimension(nx_block*ny_block) :: &
1398 5781888 : indxi, indxj ! combined i/j horizontal indices
1399 :
1400 : real (kind=dbl_kind) :: &
1401 : phi_nw, phi_n, phi_ne , & ! values of phi in 8 neighbor cells ! LCOV_EXCL_LINE
1402 : phi_w, phi_e , & ! LCOV_EXCL_LINE
1403 : phi_sw, phi_s, phi_se , & ! LCOV_EXCL_LINE
1404 : qmn, qmx , & ! min and max value of phi within grid cell ! LCOV_EXCL_LINE
1405 : pmn, pmx , & ! min and max value of phi among neighbor cells ! LCOV_EXCL_LINE
1406 725760 : w1, w2, w3, w4 ! work variables
1407 :
1408 : real (kind=dbl_kind) :: &
1409 : puny , & ! ! LCOV_EXCL_LINE
1410 725760 : gxtmp, gytmp ! temporary term for x- and y- limited gradient
1411 :
1412 : character(len=*), parameter :: subname = '(limited_gradient)'
1413 :
1414 2890944 : call icepack_query_parameters(puny_out=puny)
1415 2890944 : call icepack_warnings_flush(nu_diag)
1416 2890944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1417 0 : file=__FILE__, line=__LINE__)
1418 :
1419 2039476320 : gx(:,:) = c0
1420 2039476320 : gy(:,:) = c0
1421 :
1422 : ! For nghost = 1, loop over physical cells and update ghost cells later
1423 : ! For nghost = 2, loop over a layer of ghost cells and skip the update
1424 :
1425 2890944 : icells = 0
1426 89168688 : do j = jlo-nghost+1, jhi+nghost-1
1427 1742192928 : do i = ilo-nghost+1, ihi+nghost-1
1428 1739301984 : if (phimask(i,j) > puny) then
1429 858654723 : icells = icells + 1
1430 858654723 : indxi(icells) = i
1431 858654723 : indxj(icells) = j
1432 : endif ! phimask > puny
1433 : enddo
1434 : enddo
1435 :
1436 861545667 : do ij = 1, icells
1437 858654723 : i = indxi(ij)
1438 858654723 : j = indxj(ij)
1439 :
1440 : ! Store values of phi in the 8 neighbor cells.
1441 : ! Note: phimask = 1. or 0. If phimask = 1., use the true value;
1442 : ! if phimask = 0., use the home cell value so that non-physical
1443 : ! values of phi do not contribute to the gradient.
1444 364601024 : phi_nw = phimask(i-1,j+1) * phi(i-1,j+1) &
1445 1223255747 : + (c1-phimask(i-1,j+1))* phi(i ,j )
1446 182300512 : phi_n = phimask(i ,j+1) * phi(i ,j+1) &
1447 1040955235 : + (c1-phimask(i ,j+1))* phi(i ,j )
1448 364601024 : phi_ne = phimask(i+1,j+1) * phi(i+1,j+1) &
1449 1223255747 : + (c1-phimask(i+1,j+1))* phi(i ,j )
1450 273450768 : phi_w = phimask(i-1,j ) * phi(i-1,j ) &
1451 1040955235 : + (c1-phimask(i-1,j ))* phi(i ,j )
1452 273450768 : phi_e = phimask(i+1,j ) * phi(i+1,j ) &
1453 1040955235 : + (c1-phimask(i+1,j ))* phi(i ,j )
1454 364601024 : phi_sw = phimask(i-1,j-1) * phi(i-1,j-1) &
1455 1223255747 : + (c1-phimask(i-1,j-1))* phi(i ,j )
1456 182300512 : phi_s = phimask(i ,j-1) * phi(i ,j-1) &
1457 1040955235 : + (c1-phimask(i ,j-1))* phi(i ,j )
1458 364601024 : phi_se = phimask(i+1,j-1) * phi(i+1,j-1) &
1459 1223255747 : + (c1-phimask(i+1,j-1))* phi(i ,j )
1460 :
1461 : ! unlimited gradient components
1462 : ! (factors of two cancel out)
1463 :
1464 858654723 : gxtmp = (phi_e - phi_w) * p5
1465 858654723 : gytmp = (phi_n - phi_s) * p5
1466 :
1467 : ! minimum and maximum among the nine local cells
1468 91150256 : pmn = min (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
1469 858654723 : phi_e, phi_sw, phi_s, phi_se)
1470 91150256 : pmx = max (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
1471 858654723 : phi_e, phi_sw, phi_s, phi_se)
1472 :
1473 858654723 : pmn = pmn - phi(i,j)
1474 858654723 : pmx = pmx - phi(i,j)
1475 :
1476 : ! minimum and maximum deviation of phi within the cell
1477 91150256 : w1 = (p5 - cnx(i,j)) * gxtmp &
1478 858654723 : + (p5 - cny(i,j)) * gytmp
1479 91150256 : w2 = (p5 - cnx(i,j)) * gxtmp &
1480 858654723 : - (p5 + cny(i,j)) * gytmp
1481 91150256 : w3 = -(p5 + cnx(i,j)) * gxtmp &
1482 858654723 : - (p5 + cny(i,j)) * gytmp
1483 91150256 : w4 = (p5 - cny(i,j)) * gytmp &
1484 858654723 : - (p5 + cnx(i,j)) * gxtmp
1485 :
1486 858654723 : qmn = min (w1, w2, w3, w4)
1487 858654723 : qmx = max (w1, w2, w3, w4)
1488 :
1489 : ! the limiting coefficient
1490 858654723 : if ( abs(qmn) > abs(pmn) ) then ! 'abs(qmn) > puny' not sufficient
1491 86773338 : w1 = max(c0, pmn/qmn)
1492 : else
1493 771881385 : w1 = c1
1494 : endif
1495 :
1496 858654723 : if ( abs(qmx) > abs(pmx) ) then
1497 65479305 : w2 = max(c0, pmx/qmx)
1498 : else
1499 793175418 : w2 = c1
1500 : endif
1501 :
1502 858654723 : w1 = min(w1, w2)
1503 :
1504 : ! Limit the gradient components
1505 858654723 : gx(i,j) = w1 * gxtmp
1506 861545667 : gy(i,j) = w1 * gytmp
1507 :
1508 : enddo ! ij
1509 :
1510 2890944 : end subroutine limited_gradient
1511 :
1512 : !=======================================================================
1513 : !
1514 : ! Given velocity fields on cell corners, compute departure points
1515 : ! of back trajectories in nondimensional coordinates.
1516 : !
1517 : ! author William H. Lipscomb, LANL
1518 :
1519 22944 : subroutine departure_points (nx_block, ny_block, &
1520 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1521 : nghost, dt, & ! LCOV_EXCL_LINE
1522 : uvel, vvel, & ! LCOV_EXCL_LINE
1523 : dxu, dyu, & ! LCOV_EXCL_LINE
1524 : HTN, HTE, & ! LCOV_EXCL_LINE
1525 : dpx, dpy, & ! LCOV_EXCL_LINE
1526 : l_dp_midpt, l_stop, & ! LCOV_EXCL_LINE
1527 : istop, jstop)
1528 :
1529 : integer (kind=int_kind), intent(in) :: &
1530 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1531 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
1532 : nghost ! number of ghost cells
1533 :
1534 : real (kind=dbl_kind), intent(in) :: &
1535 : dt ! time step (s)
1536 :
1537 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1538 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1539 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1540 : dxu , & ! E-W dimensions of U-cell (m) ! LCOV_EXCL_LINE
1541 : dyu , & ! N-S dimensions of U-cell (m) ! LCOV_EXCL_LINE
1542 : HTN , & ! length of north face of T-cell (m) ! LCOV_EXCL_LINE
1543 : HTE ! length of east face of T-cell (m)
1544 :
1545 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1546 : dpx , & ! coordinates of departure points (m) ! LCOV_EXCL_LINE
1547 : dpy ! coordinates of departure points (m)
1548 :
1549 : logical (kind=log_kind), intent(in) :: &
1550 : l_dp_midpt ! if true, find departure points using
1551 : ! corrected midpoint velocity
1552 :
1553 : logical (kind=log_kind), intent(inout) :: &
1554 : l_stop ! if true, abort on return
1555 :
1556 : integer (kind=int_kind), intent(inout) :: &
1557 : istop, jstop ! indices of grid cell where model aborts
1558 :
1559 : ! local variables
1560 :
1561 : integer (kind=int_kind) :: &
1562 : i, j, i2, j2 ! horizontal indices
1563 :
1564 : real (kind=dbl_kind) :: &
1565 : mpx, mpy , & ! coordinates of midpoint of back trajectory, ! LCOV_EXCL_LINE
1566 : ! relative to cell corner
1567 5760 : mpxt, mpyt , & ! midpoint coordinates relative to cell center
1568 5760 : ump, vmp ! corrected velocity at midpoint
1569 :
1570 : character(len=*), parameter :: subname = '(departure_points)'
1571 :
1572 : !-------------------------------------------------------------------
1573 : ! Estimate departure points.
1574 : ! This estimate is 1st-order accurate in time; improve accuracy by
1575 : ! using midpoint approximation (to add later).
1576 : ! For nghost = 1, loop over physical cells and update ghost cells later.
1577 : ! For nghost = 2, loop over a layer of ghost cells and skip update.
1578 : !-------------------------------------------------------------------
1579 :
1580 16186320 : dpx(:,:) = c0
1581 16186320 : dpy(:,:) = c0
1582 :
1583 707688 : do j = jlo-nghost+1, jhi+nghost-1
1584 13826928 : do i = ilo-nghost+1, ihi+nghost-1
1585 :
1586 13119240 : dpx(i,j) = -dt*uvel(i,j)
1587 13119240 : dpy(i,j) = -dt*vvel(i,j)
1588 :
1589 : ! Check for values out of bounds (more than one grid cell away)
1590 8352000 : if (dpx(i,j) < -HTN(i,j) .or. dpx(i,j) > HTN(i+1,j) .or. &
1591 17979984 : dpy(i,j) < -HTE(i,j) .or. dpy(i,j) > HTE(i,j+1)) then
1592 0 : l_stop = .true.
1593 0 : istop = i
1594 0 : jstop = j
1595 : endif
1596 :
1597 : enddo
1598 : enddo
1599 :
1600 22944 : if (l_stop) then
1601 0 : i = istop
1602 0 : j = jstop
1603 0 : write (nu_diag,*) ' '
1604 : write (nu_diag,*) &
1605 0 : 'Warning: Departure points out of bounds in remap'
1606 0 : write (nu_diag,*) 'my_task, i, j =', my_task, i, j
1607 0 : write (nu_diag,*) 'dpx, dpy =', dpx(i,j), dpy(i,j)
1608 0 : write (nu_diag,*) 'HTN(i,j), HTN(i+1,j) =', HTN(i,j), HTN(i+1,j)
1609 0 : write (nu_diag,*) 'HTE(i,j), HTE(i,j+1) =', HTE(i,j), HTE(i,j+1)
1610 0 : return
1611 : endif
1612 :
1613 22944 : if (l_dp_midpt) then ! find dep pts using corrected midpt velocity
1614 :
1615 707688 : do j = jlo-nghost+1, jhi+nghost-1
1616 13826928 : do i = ilo-nghost+1, ihi+nghost-1
1617 13803984 : if (uvel(i,j)/=c0 .or. vvel(i,j)/=c0) then
1618 :
1619 : !-------------------------------------------------------------------
1620 : ! Scale departure points to coordinate system in which grid cells
1621 : ! have sides of unit length.
1622 : !-------------------------------------------------------------------
1623 :
1624 6452519 : dpx(i,j) = dpx(i,j) / dxu(i,j)
1625 6452519 : dpy(i,j) = dpy(i,j) / dyu(i,j)
1626 :
1627 : !-------------------------------------------------------------------
1628 : ! Estimate midpoint of backward trajectory relative to corner (i,j).
1629 : !-------------------------------------------------------------------
1630 :
1631 6452519 : mpx = p5 * dpx(i,j)
1632 6452519 : mpy = p5 * dpy(i,j)
1633 :
1634 : !-------------------------------------------------------------------
1635 : ! Determine the indices (i2,j2) of the cell where the trajectory lies.
1636 : ! Compute the coordinates of the midpoint of the backward trajectory
1637 : ! relative to the cell center in a stretch coordinate system
1638 : ! with vertices at (1/2, 1/2), (1/2, -1/2), etc.
1639 : !-------------------------------------------------------------------
1640 :
1641 6452519 : if (mpx >= c0 .and. mpy >= c0) then ! cell (i+1,j+1)
1642 214554 : i2 = i+1
1643 214554 : j2 = j+1
1644 214554 : mpxt = mpx - p5
1645 214554 : mpyt = mpy - p5
1646 6237965 : elseif (mpx < c0 .and. mpy < c0) then ! cell (i,j)
1647 5198378 : i2 = i
1648 5198378 : j2 = j
1649 5198378 : mpxt = mpx + p5
1650 5198378 : mpyt = mpy + p5
1651 1039587 : elseif (mpx >= c0 .and. mpy < c0) then ! cell (i+1,j)
1652 226767 : i2 = i+1
1653 226767 : j2 = j
1654 226767 : mpxt = mpx - p5
1655 226767 : mpyt = mpy + p5
1656 812820 : elseif (mpx < c0 .and. mpy >= c0) then ! cell (i,j+1)
1657 812820 : i2 = i
1658 812820 : j2 = j+1
1659 812820 : mpxt = mpx + p5
1660 812820 : mpyt = mpy - p5
1661 : endif
1662 :
1663 : !-------------------------------------------------------------------
1664 : ! Using a bilinear approximation, estimate the velocity at the
1665 : ! trajectory midpoint in the (i2,j2) reference frame.
1666 : !-------------------------------------------------------------------
1667 :
1668 944754 : ump = uvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
1669 : - uvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) & ! LCOV_EXCL_LINE
1670 : + uvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) & ! LCOV_EXCL_LINE
1671 7869650 : - uvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
1672 :
1673 944754 : vmp = vvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
1674 : - vvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) & ! LCOV_EXCL_LINE
1675 : + vvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) & ! LCOV_EXCL_LINE
1676 7869650 : - vvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
1677 :
1678 : !-------------------------------------------------------------------
1679 : ! Use the midpoint velocity to estimate the coordinates of the
1680 : ! departure point relative to corner (i,j).
1681 : !-------------------------------------------------------------------
1682 :
1683 6452519 : dpx(i,j) = -dt * ump
1684 6452519 : dpy(i,j) = -dt * vmp
1685 :
1686 : endif ! nonzero velocity
1687 :
1688 : enddo ! i
1689 : enddo ! j
1690 :
1691 : endif ! l_dp_midpt
1692 :
1693 : end subroutine departure_points
1694 :
1695 : !=======================================================================
1696 : !
1697 : ! Compute areas and vertices of transport triangles for north or
1698 : ! east cell edges.
1699 : !
1700 : ! authors William H. Lipscomb, LANL
1701 : ! John R. Baumgardner, LANL
1702 :
1703 45888 : subroutine locate_triangles (nx_block, ny_block, &
1704 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1705 : nghost, edge, & ! LCOV_EXCL_LINE
1706 : icells, & ! LCOV_EXCL_LINE
1707 : indxi, indxj, & ! LCOV_EXCL_LINE
1708 : dpx, dpy, & ! LCOV_EXCL_LINE
1709 : dxu, dyu, & ! LCOV_EXCL_LINE
1710 : earea, narea, & ! LCOV_EXCL_LINE
1711 : xp, yp, & ! LCOV_EXCL_LINE
1712 : iflux, jflux, & ! LCOV_EXCL_LINE
1713 45888 : triarea, edgearea)
1714 :
1715 : integer (kind=int_kind), intent(in) :: &
1716 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1717 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
1718 : nghost ! number of ghost cells
1719 :
1720 : character (len=char_len), intent(in) :: &
1721 : edge ! 'north' or 'east'
1722 :
1723 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
1724 : dpx , & ! x coordinates of departure points at cell corners ! LCOV_EXCL_LINE
1725 : dpy , & ! y coordinates of departure points at cell corners ! LCOV_EXCL_LINE
1726 : dxu , & ! E-W dimension of U-cell (m) ! LCOV_EXCL_LINE
1727 : dyu , & ! N-S dimension of U-cell (m) ! LCOV_EXCL_LINE
1728 : earea , & ! area of E-cell ! LCOV_EXCL_LINE
1729 : narea ! area of N-cell
1730 :
1731 : real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups), intent(out) :: &
1732 : xp, yp ! coordinates of triangle vertices
1733 :
1734 : real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups), intent(out) :: &
1735 : triarea ! area of departure triangle
1736 :
1737 : integer (kind=int_kind), dimension (nx_block,ny_block,ngroups), intent(out) :: &
1738 : iflux , & ! i index of cell contributing transport ! LCOV_EXCL_LINE
1739 : jflux ! j index of cell contributing transport
1740 :
1741 : integer (kind=int_kind), dimension (ngroups), intent(out) :: &
1742 : icells ! number of cells where triarea > puny
1743 :
1744 : integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), intent(out) :: &
1745 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1746 : indxj ! compressed index in j-direction
1747 :
1748 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
1749 : edgearea ! area of departure region for each edge
1750 : ! edgearea > 0 for eastward/northward flow
1751 :
1752 : ! local variables
1753 :
1754 : integer (kind=int_kind) :: &
1755 : i, j, ij, ic , & ! horizontal indices ! LCOV_EXCL_LINE
1756 : ib, jb , & ! limits for loops for bugcheck ! LCOV_EXCL_LINE
1757 : ng, nv , & ! triangle indices ! LCOV_EXCL_LINE
1758 : ishift , jshift , & ! differences between neighbor cells ! LCOV_EXCL_LINE
1759 : ishift_tl, jshift_tl , & ! i,j indices of TL cell relative to edge ! LCOV_EXCL_LINE
1760 : ishift_bl, jshift_bl , & ! i,j indices of BL cell relative to edge ! LCOV_EXCL_LINE
1761 : ishift_tr, jshift_tr , & ! i,j indices of TR cell relative to edge ! LCOV_EXCL_LINE
1762 : ishift_br, jshift_br , & ! i,j indices of BR cell relative to edge ! LCOV_EXCL_LINE
1763 : ishift_tc, jshift_tc , & ! i,j indices of TC cell relative to edge ! LCOV_EXCL_LINE
1764 : ishift_bc, jshift_bc , & ! i,j indices of BC cell relative to edge ! LCOV_EXCL_LINE
1765 : is_l, js_l , & ! i,j shifts for TL1, BL2 for area consistency ! LCOV_EXCL_LINE
1766 : is_r, js_r , & ! i,j shifts for TR1, BR2 for area consistency ! LCOV_EXCL_LINE
1767 : ise_tl, jse_tl , & ! i,j of TL other edge relative to edge ! LCOV_EXCL_LINE
1768 : ise_bl, jse_bl , & ! i,j of BL other edge relative to edge ! LCOV_EXCL_LINE
1769 : ise_tr, jse_tr , & ! i,j of TR other edge relative to edge ! LCOV_EXCL_LINE
1770 : ise_br, jse_br ! i,j of BR other edge relative to edge
1771 :
1772 : integer (kind=int_kind) :: &
1773 : icellsd ! number of cells where departure area > 0.
1774 :
1775 : integer (kind=int_kind), dimension (nx_block*ny_block) :: &
1776 : indxid , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1777 91776 : indxjd ! compressed index in j-direction
1778 :
1779 : real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
1780 : dx, dy , & ! scaled departure points ! LCOV_EXCL_LINE
1781 : areafac_c , & ! earea or narea ! LCOV_EXCL_LINE
1782 10079616 : areafac_ce ! areafac_c on other edge (narea or earea)
1783 :
1784 : real (kind=dbl_kind) :: &
1785 : xcl, ycl , & ! coordinates of left corner point ! LCOV_EXCL_LINE
1786 : ! (relative to midpoint of edge)
1787 11520 : xdl, ydl , & ! left departure point
1788 : xil, yil , & ! left intersection point ! LCOV_EXCL_LINE
1789 : xcr, ycr , & ! right corner point ! LCOV_EXCL_LINE
1790 : xdr, ydr , & ! right departure point ! LCOV_EXCL_LINE
1791 : xir, yir , & ! right intersection point ! LCOV_EXCL_LINE
1792 : xic, yic , & ! x-axis intersection point ! LCOV_EXCL_LINE
1793 : xicl, yicl , & ! left-hand x-axis intersection point ! LCOV_EXCL_LINE
1794 : xicr, yicr , & ! right-hand x-axis intersection point ! LCOV_EXCL_LINE
1795 : xdm, ydm , & ! midpoint of segment connecting DL and DR; ! LCOV_EXCL_LINE
1796 : ! shifted if l_fixed_area = T
1797 11520 : md , & ! slope of line connecting DL and DR
1798 : mdl , & ! slope of line connecting DL and DM ! LCOV_EXCL_LINE
1799 : mdr , & ! slope of line connecting DR and DM ! LCOV_EXCL_LINE
1800 : area1, area2 , & ! temporary triangle areas ! LCOV_EXCL_LINE
1801 : area3, area4 , & ! ! LCOV_EXCL_LINE
1802 : area_c , & ! center polygon area ! LCOV_EXCL_LINE
1803 : puny , & ! ! LCOV_EXCL_LINE
1804 11520 : w1, w2 ! work variables
1805 :
1806 : real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups) :: &
1807 60145536 : areafact ! = 1 for positive flux, -1 for negative
1808 :
1809 : real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
1810 10091136 : areasum ! sum of triangle areas for a given edge
1811 :
1812 : character(len=*), parameter :: subname = '(locate_triangles)'
1813 :
1814 : !-------------------------------------------------------------------
1815 : ! Triangle notation:
1816 : ! For each edge, there are 20 triangles that can contribute,
1817 : ! but many of these are mutually exclusive. It turns out that
1818 : ! at most 5 triangles can contribute to transport integrals at once.
1819 : !
1820 : ! See Figure 3 in DB for pictures of these triangles.
1821 : ! See Table 1 in DB for logical conditions.
1822 : !
1823 : ! For the north edge, DB refer to these triangles as:
1824 : ! (1) NW, NW1, W, W2
1825 : ! (2) NE, NE1, E, E2
1826 : ! (3) NW2, W1, NE2, E1
1827 : ! (4) H1a, H1b, N1a, N1b
1828 : ! (5) H2a, H2b, N2a, N2b
1829 : !
1830 : ! For the east edge, DB refer to these triangles as:
1831 : ! (1) NE, NE1, N, N2
1832 : ! (2) SE, SE1, S, S2
1833 : ! (3) NE2, N1, SE2, S1
1834 : ! (4) H1a, H1b, E1a, E2b
1835 : ! (5) H2a, H2b, E2a, E2b
1836 : !
1837 : ! The code below works for either north or east edges.
1838 : ! The respective triangle labels are:
1839 : ! (1) TL, TL1, BL, BL2
1840 : ! (2) TR, TR1, BR, BR2
1841 : ! (3) TL2, BL1, TR2, BR1
1842 : ! (4) BC1a, BC1b, TC1a, TC2b
1843 : ! (5) BC2a, BC2b, TC2a, TC2b
1844 : !
1845 : ! where the cell labels are:
1846 : !
1847 : ! | |
1848 : ! TL | TC | TR (top left, center, right)
1849 : ! | |
1850 : ! ------------------------
1851 : ! | |
1852 : ! BL | BC | BR (bottom left, center, right)
1853 : ! | |
1854 : !
1855 : ! and the transport is across the edge between cells TC and BC.
1856 : !
1857 : ! Departure points are scaled to a local coordinate system
1858 : ! whose origin is at the midpoint of the edge.
1859 : ! In this coordinate system, the lefthand corner CL = (-0.5,0)
1860 : ! and the righthand corner CR = (0.5, 0).
1861 : !-------------------------------------------------------------------
1862 :
1863 : !-------------------------------------------------------------------
1864 : ! Initialize
1865 : !-------------------------------------------------------------------
1866 :
1867 45888 : call icepack_query_parameters(puny_out=puny)
1868 45888 : call icepack_warnings_flush(nu_diag)
1869 45888 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1870 0 : file=__FILE__, line=__LINE__)
1871 :
1872 32372640 : areafac_c(:,:) = c0
1873 32372640 : areafac_ce(:,:) = c0
1874 :
1875 321216 : do ng = 1, ngroups
1876 9042912 : do j = 1, ny_block
1877 194235840 : do i = 1, nx_block
1878 185192928 : triarea (i,j,ng) = c0
1879 185192928 : areafact(i,j,ng) = c0
1880 185192928 : iflux (i,j,ng) = i
1881 193960512 : jflux (i,j,ng) = j
1882 : enddo
1883 : enddo
1884 1422528 : do nv = 0, nvert
1885 36446976 : do j = 1, ny_block
1886 776943360 : do i = 1, nx_block
1887 740771712 : xp(i,j,nv,ng) = c0
1888 775842048 : yp(i,j,nv,ng) = c0
1889 : enddo
1890 : enddo
1891 : enddo
1892 : enddo
1893 :
1894 45888 : if (trim(edge) == 'north') then
1895 :
1896 : ! index shifts for neighbor cells
1897 :
1898 22944 : ishift_tl = -1
1899 22944 : jshift_tl = 1
1900 22944 : ishift_bl = -1
1901 22944 : jshift_bl = 0
1902 22944 : ishift_tr = 1
1903 22944 : jshift_tr = 1
1904 22944 : ishift_br = 1
1905 22944 : jshift_br = 0
1906 22944 : ishift_tc = 0
1907 22944 : jshift_tc = 1
1908 22944 : ishift_bc = 0
1909 22944 : jshift_bc = 0
1910 :
1911 : ! index shifts for TL1, BL2, TR1 and BR2 for area consistency
1912 :
1913 22944 : is_l = -1
1914 22944 : js_l = 0
1915 22944 : is_r = 1
1916 22944 : js_r = 0
1917 :
1918 : ! index shifts for neighbor east edges
1919 :
1920 22944 : ise_tl = -1
1921 22944 : jse_tl = 1
1922 22944 : ise_bl = -1
1923 22944 : jse_bl = 0
1924 22944 : ise_tr = 0
1925 22944 : jse_tr = 1
1926 22944 : ise_br = 0
1927 22944 : jse_br = 0
1928 :
1929 : ! area scale factor
1930 : ! earea, narea valid on halo
1931 :
1932 753576 : do j = 1, ny_block
1933 16186320 : do i = 1, nx_block
1934 16163376 : areafac_c(i,j) = narea(i,j)
1935 : enddo
1936 : enddo
1937 :
1938 : ! area scale factor for other edge (east)
1939 :
1940 753576 : do j = 1, ny_block
1941 16186320 : do i = 1, nx_block
1942 16163376 : areafac_ce(i,j) = earea(i,j)
1943 : enddo
1944 : enddo
1945 :
1946 : else ! east edge
1947 :
1948 : ! index shifts for neighbor cells
1949 :
1950 22944 : ishift_tl = 1
1951 22944 : jshift_tl = 1
1952 22944 : ishift_bl = 0
1953 22944 : jshift_bl = 1
1954 22944 : ishift_tr = 1
1955 22944 : jshift_tr = -1
1956 22944 : ishift_br = 0
1957 22944 : jshift_br = -1
1958 22944 : ishift_tc = 1
1959 22944 : jshift_tc = 0
1960 22944 : ishift_bc = 0
1961 22944 : jshift_bc = 0
1962 :
1963 : ! index shifts for TL1, BL2, TR1 and BR2 for area consistency
1964 :
1965 22944 : is_l = 0
1966 22944 : js_l = 1
1967 22944 : is_r = 0
1968 22944 : js_r = -1
1969 :
1970 : ! index shifts for neighbor north edges
1971 :
1972 22944 : ise_tl = 1
1973 22944 : jse_tl = 0
1974 22944 : ise_bl = 0
1975 22944 : jse_bl = 0
1976 22944 : ise_tr = 1
1977 22944 : jse_tr = -1
1978 22944 : ise_br = 0
1979 22944 : jse_br = -1
1980 :
1981 : ! area scale factors
1982 : ! earea, narea valid on halo
1983 :
1984 753576 : do j = 1, ny_block
1985 16186320 : do i = 1, nx_block
1986 16163376 : areafac_c(i,j) = earea(i,j)
1987 : enddo
1988 : enddo
1989 :
1990 : ! area scale factor for other edge (north)
1991 :
1992 753576 : do j = 1, ny_block
1993 16186320 : do i = 1, nx_block
1994 16163376 : areafac_ce(i,j) = narea(i,j)
1995 : enddo
1996 : enddo
1997 :
1998 : endif
1999 :
2000 : !-------------------------------------------------------------------
2001 : ! Compute mask for edges with nonzero departure areas and for
2002 : ! one grid-cell wide channels
2003 : !-------------------------------------------------------------------
2004 :
2005 45888 : icellsd = 0
2006 45888 : if (trim(edge) == 'north') then
2007 730632 : do j = jlo-1, jhi
2008 14275992 : do i = ilo, ihi
2009 12960000 : if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 &
2010 : .or. & ! LCOV_EXCL_LINE
2011 22893048 : dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
2012 6691810 : icellsd = icellsd + 1
2013 6691810 : indxid(icellsd) = i
2014 6691810 : indxjd(icellsd) = j
2015 : else
2016 6853550 : if ( abs(edgearea(i,j)) > c0 ) then ! 1 grid-cell wide channel: dpx,y = 0, edgearea /= 0
2017 0 : icellsd = icellsd + 1
2018 0 : indxid(icellsd) = i
2019 0 : indxjd(icellsd) = j
2020 : endif
2021 : endif
2022 : enddo
2023 : enddo
2024 : else ! east edge
2025 707688 : do j = jlo, jhi
2026 14511672 : do i = ilo-1, ihi
2027 8686080 : if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 &
2028 : .or. & ! LCOV_EXCL_LINE
2029 23174808 : dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
2030 6844823 : icellsd = icellsd + 1
2031 6844823 : indxid(icellsd) = i
2032 6844823 : indxjd(icellsd) = j
2033 : else
2034 6959161 : if ( abs(edgearea(i,j)) > c0 ) then ! 1 grid-cell wide channel: dpx,y = 0, edgearea /= 0
2035 0 : icellsd = icellsd + 1
2036 0 : indxid(icellsd) = i
2037 0 : indxjd(icellsd) = j
2038 : endif
2039 : endif
2040 : enddo
2041 : enddo
2042 : endif ! edge = north/east
2043 :
2044 : !-------------------------------------------------------------------
2045 : ! Scale the departure points
2046 : !-------------------------------------------------------------------
2047 :
2048 1461264 : do j = 1, jhi
2049 29967360 : do i = 1, ihi
2050 28506096 : dx(i,j) = dpx(i,j) / dxu(i,j)
2051 29921472 : dy(i,j) = dpy(i,j) / dyu(i,j)
2052 : enddo
2053 : enddo
2054 :
2055 : !-------------------------------------------------------------------
2056 : ! Compute departure regions, divide into triangles, and locate
2057 : ! vertices of each triangle.
2058 : ! Work in a nondimensional coordinate system in which lengths are
2059 : ! scaled by the local metric coefficients (dxu and dyu).
2060 : ! Note: The do loop includes north faces of the j = 1 ghost cells
2061 : ! when edge = 'north'. The loop includes east faces of i = 1
2062 : ! ghost cells when edge = 'east'.
2063 : !-------------------------------------------------------------------
2064 :
2065 13582521 : do ij = 1, icellsd
2066 13536633 : i = indxid(ij)
2067 13536633 : j = indxjd(ij)
2068 :
2069 13536633 : xcl = -p5
2070 13536633 : ycl = c0
2071 :
2072 13536633 : xcr = p5
2073 13536633 : ycr = c0
2074 :
2075 : ! Departure points
2076 :
2077 13536633 : if (trim(edge) == 'north') then ! north edge
2078 6691810 : xdl = xcl + dx(i-1,j )
2079 6691810 : ydl = ycl + dy(i-1,j )
2080 6691810 : xdr = xcr + dx(i ,j )
2081 6691810 : ydr = ycr + dy(i ,j )
2082 : else ! east edge; rotate trajectory by pi/2
2083 6844823 : xdl = xcl - dy(i ,j )
2084 6844823 : ydl = ycl + dx(i ,j )
2085 6844823 : xdr = xcr - dy(i ,j-1)
2086 6844823 : ydr = ycr + dx(i ,j-1)
2087 : endif
2088 :
2089 13536633 : xdm = p5 * (xdr + xdl)
2090 13536633 : ydm = p5 * (ydr + ydl)
2091 :
2092 : ! Intersection points
2093 :
2094 13536633 : xil = xcl
2095 13536633 : yil = (xcl*(ydm-ydl) + xdm*ydl - xdl*ydm) / (xdm - xdl)
2096 :
2097 13536633 : xir = xcr
2098 13536633 : yir = (xcr*(ydr-ydm) - xdm*ydr + xdr*ydm) / (xdr - xdm)
2099 :
2100 13536633 : md = (ydr - ydl) / (xdr - xdl)
2101 :
2102 13536633 : if (abs(md) > puny) then
2103 13536631 : xic = xdl - ydl/md
2104 : else
2105 2 : xic = c0
2106 : endif
2107 13536633 : yic = c0
2108 :
2109 13536633 : xicl = xic
2110 13536633 : yicl = yic
2111 13536633 : xicr = xic
2112 13536633 : yicr = yic
2113 :
2114 : !-------------------------------------------------------------------
2115 : ! Locate triangles in TL cell (NW for north edge, NE for east edge)
2116 : ! and BL cell (W for north edge, N for east edge).
2117 : !
2118 : ! areafact_c or areafac_ce (areafact_c for the other edge) are used
2119 : ! (with shifted indices) to make sure that a flux area on one edge
2120 : ! is consistent with the analogous area on the other edge and to
2121 : ! ensure that areas add up when using l_fixed_area = T. See PR #849
2122 : ! for details.
2123 : !
2124 : !-------------------------------------------------------------------
2125 :
2126 13536633 : if (yil > c0 .and. xdl < xcl .and. ydl >= c0) then
2127 :
2128 : ! TL (group 1)
2129 :
2130 1055992 : ng = 1
2131 1055992 : xp (i,j,1,ng) = xcl
2132 1055992 : yp (i,j,1,ng) = ycl
2133 1055992 : xp (i,j,2,ng) = xil
2134 1055992 : yp (i,j,2,ng) = yil
2135 1055992 : xp (i,j,3,ng) = xdl
2136 1055992 : yp (i,j,3,ng) = ydl
2137 1055992 : iflux (i,j,ng) = i + ishift_tl
2138 1055992 : jflux (i,j,ng) = j + jshift_tl
2139 1055992 : areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl)
2140 :
2141 12480641 : elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then
2142 :
2143 : ! BL (group 1)
2144 :
2145 6163929 : ng = 1
2146 6163929 : xp (i,j,1,ng) = xcl
2147 6163929 : yp (i,j,1,ng) = ycl
2148 6163929 : xp (i,j,2,ng) = xdl
2149 6163929 : yp (i,j,2,ng) = ydl
2150 6163929 : xp (i,j,3,ng) = xil
2151 6163929 : yp (i,j,3,ng) = yil
2152 6163929 : iflux (i,j,ng) = i + ishift_bl
2153 6163929 : jflux (i,j,ng) = j + jshift_bl
2154 6163929 : areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl)
2155 :
2156 6316712 : elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then
2157 :
2158 : ! TL1 (group 1)
2159 :
2160 172 : ng = 1
2161 172 : xp (i,j,1,ng) = xcl
2162 172 : yp (i,j,1,ng) = ycl
2163 172 : xp (i,j,2,ng) = xdl
2164 172 : yp (i,j,2,ng) = ydl
2165 172 : xp (i,j,3,ng) = xic
2166 172 : yp (i,j,3,ng) = yic
2167 172 : iflux (i,j,ng) = i + ishift_tl
2168 172 : jflux (i,j,ng) = j + jshift_tl
2169 172 : areafact(i,j,ng) = areafac_c(i+is_l,j+js_l)
2170 :
2171 : ! BL1 (group 3)
2172 :
2173 172 : ng = 3
2174 172 : xp (i,j,1,ng) = xcl
2175 172 : yp (i,j,1,ng) = ycl
2176 172 : xp (i,j,2,ng) = xic
2177 172 : yp (i,j,2,ng) = yic
2178 172 : xp (i,j,3,ng) = xil
2179 172 : yp (i,j,3,ng) = yil
2180 172 : iflux (i,j,ng) = i + ishift_bl
2181 172 : jflux (i,j,ng) = j + jshift_bl
2182 172 : areafact(i,j,ng) = areafac_ce(i+ise_bl,j+jse_bl)
2183 :
2184 6316540 : elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then
2185 :
2186 : ! TL2 (group 3)
2187 :
2188 217 : ng = 3
2189 217 : xp (i,j,1,ng) = xcl
2190 217 : yp (i,j,1,ng) = ycl
2191 217 : xp (i,j,2,ng) = xil
2192 217 : yp (i,j,2,ng) = yil
2193 217 : xp (i,j,3,ng) = xic
2194 217 : yp (i,j,3,ng) = yic
2195 217 : iflux (i,j,ng) = i + ishift_tl
2196 217 : jflux (i,j,ng) = j + jshift_tl
2197 217 : areafact(i,j,ng) = -areafac_ce(i+ise_tl,j+jse_tl)
2198 :
2199 : ! BL2 (group 1)
2200 :
2201 217 : ng = 1
2202 217 : xp (i,j,1,ng) = xcl
2203 217 : yp (i,j,1,ng) = ycl
2204 217 : xp (i,j,2,ng) = xic
2205 217 : yp (i,j,2,ng) = yic
2206 217 : xp (i,j,3,ng) = xdl
2207 217 : yp (i,j,3,ng) = ydl
2208 217 : iflux (i,j,ng) = i + ishift_bl
2209 217 : jflux (i,j,ng) = j + jshift_bl
2210 217 : areafact(i,j,ng) = -areafac_c(i+is_l,j+js_l)
2211 :
2212 : endif ! TL and BL triangles
2213 :
2214 : !-------------------------------------------------------------------
2215 : ! Locate triangles in TR cell (NE for north edge, SE for east edge)
2216 : ! and in BR cell (E for north edge, S for east edge).
2217 : !-------------------------------------------------------------------
2218 :
2219 13536633 : if (yir > c0 .and. xdr >= xcr .and. ydr >= c0) then
2220 :
2221 : ! TR (group 2)
2222 :
2223 455946 : ng = 2
2224 455946 : xp (i,j,1,ng) = xcr
2225 455946 : yp (i,j,1,ng) = ycr
2226 455946 : xp (i,j,2,ng) = xdr
2227 455946 : yp (i,j,2,ng) = ydr
2228 455946 : xp (i,j,3,ng) = xir
2229 455946 : yp (i,j,3,ng) = yir
2230 455946 : iflux (i,j,ng) = i + ishift_tr
2231 455946 : jflux (i,j,ng) = j + jshift_tr
2232 455946 : areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr)
2233 :
2234 13080687 : elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then
2235 :
2236 : ! BR (group 2)
2237 :
2238 5555742 : ng = 2
2239 5555742 : xp (i,j,1,ng) = xcr
2240 5555742 : yp (i,j,1,ng) = ycr
2241 5555742 : xp (i,j,2,ng) = xir
2242 5555742 : yp (i,j,2,ng) = yir
2243 5555742 : xp (i,j,3,ng) = xdr
2244 5555742 : yp (i,j,3,ng) = ydr
2245 5555742 : iflux (i,j,ng) = i + ishift_br
2246 5555742 : jflux (i,j,ng) = j + jshift_br
2247 5555742 : areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br)
2248 :
2249 7524945 : elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then
2250 :
2251 : ! TR1 (group 2)
2252 :
2253 212 : ng = 2
2254 212 : xp (i,j,1,ng) = xcr
2255 212 : yp (i,j,1,ng) = ycr
2256 212 : xp (i,j,2,ng) = xic
2257 212 : yp (i,j,2,ng) = yic
2258 212 : xp (i,j,3,ng) = xdr
2259 212 : yp (i,j,3,ng) = ydr
2260 212 : iflux (i,j,ng) = i + ishift_tr
2261 212 : jflux (i,j,ng) = j + jshift_tr
2262 212 : areafact(i,j,ng) = areafac_c(i+is_r,j+js_r)
2263 :
2264 : ! BR1 (group 3)
2265 :
2266 212 : ng = 3
2267 212 : xp (i,j,1,ng) = xcr
2268 212 : yp (i,j,1,ng) = ycr
2269 212 : xp (i,j,2,ng) = xir
2270 212 : yp (i,j,2,ng) = yir
2271 212 : xp (i,j,3,ng) = xic
2272 212 : yp (i,j,3,ng) = yic
2273 212 : iflux (i,j,ng) = i + ishift_br
2274 212 : jflux (i,j,ng) = j + jshift_br
2275 212 : areafact(i,j,ng) = areafac_ce(i+ise_br,j+jse_br)
2276 :
2277 7524733 : elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then
2278 :
2279 : ! TR2 (group 3)
2280 :
2281 85 : ng = 3
2282 85 : xp (i,j,1,ng) = xcr
2283 85 : yp (i,j,1,ng) = ycr
2284 85 : xp (i,j,2,ng) = xic
2285 85 : yp (i,j,2,ng) = yic
2286 85 : xp (i,j,3,ng) = xir
2287 85 : yp (i,j,3,ng) = yir
2288 85 : iflux (i,j,ng) = i + ishift_tr
2289 85 : jflux (i,j,ng) = j + jshift_tr
2290 85 : areafact(i,j,ng) = -areafac_ce(i+ise_tr,j+jse_tr)
2291 :
2292 : ! BR2 (group 2)
2293 :
2294 85 : ng = 2
2295 85 : xp (i,j,1,ng) = xcr
2296 85 : yp (i,j,1,ng) = ycr
2297 85 : xp (i,j,2,ng) = xdr
2298 85 : yp (i,j,2,ng) = ydr
2299 85 : xp (i,j,3,ng) = xic
2300 85 : yp (i,j,3,ng) = yic
2301 85 : iflux (i,j,ng) = i + ishift_br
2302 85 : jflux (i,j,ng) = j + jshift_br
2303 85 : areafact(i,j,ng) = -areafac_c(i+is_r,j+js_r)
2304 :
2305 : endif ! TR and BR triangles
2306 :
2307 : !-------------------------------------------------------------------
2308 : ! Redefine departure points if not located in central cells (TC or BC)
2309 : !-------------------------------------------------------------------
2310 :
2311 13536633 : if (xdl < xcl) then
2312 7220310 : xdl = xil
2313 7220310 : ydl = yil
2314 : endif
2315 :
2316 13536633 : if (xdr > xcr) then
2317 6011985 : xdr = xir
2318 6011985 : ydr = yir
2319 : endif
2320 :
2321 : !-------------------------------------------------------------------
2322 : ! For l_fixed_area = T, shift the midpoint so that the departure
2323 : ! region has the prescribed area
2324 : !-------------------------------------------------------------------
2325 :
2326 13536633 : if (l_fixed_area) then
2327 :
2328 : ! Sum the areas of the left and right triangles.
2329 : ! Note that yp(i,j,1,ng) = 0 for all triangles, so we can
2330 : ! drop those terms from the area formula.
2331 :
2332 0 : ng = 1
2333 0 : area1 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
2334 : yp(i,j,3,ng) & ! LCOV_EXCL_LINE
2335 : - yp(i,j,2,ng) * & ! LCOV_EXCL_LINE
2336 : (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & ! LCOV_EXCL_LINE
2337 0 : * areafact(i,j,ng)
2338 :
2339 0 : ng = 2
2340 0 : area2 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
2341 : yp(i,j,3,ng) & ! LCOV_EXCL_LINE
2342 : - yp(i,j,2,ng) * & ! LCOV_EXCL_LINE
2343 : (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & ! LCOV_EXCL_LINE
2344 0 : * areafact(i,j,ng)
2345 :
2346 0 : ng = 3
2347 0 : area3 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
2348 : yp(i,j,3,ng) & ! LCOV_EXCL_LINE
2349 : - yp(i,j,2,ng) * & ! LCOV_EXCL_LINE
2350 : (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & ! LCOV_EXCL_LINE
2351 0 : * areafact(i,j,ng)
2352 :
2353 : !-----------------------------------------------------------
2354 : ! Check whether the central triangles lie in one grid cell or two.
2355 : ! If all are in one grid cell, then adjust the area of the central
2356 : ! region so that the sum of all triangle areas is equal to the
2357 : ! prescribed value.
2358 : ! If two triangles are in one grid cell and one is in the other,
2359 : ! then compute the area of the lone triangle. Then adjust
2360 : ! the area of the remaining two-triangle region so that the sum of
2361 : ! all triangle areas has the prescribed value.
2362 : !-----------------------------------------------------------
2363 :
2364 0 : if (ydl*ydr >= c0) then ! Both DPs lie on same side of x-axis
2365 :
2366 : ! compute required area of central departure region
2367 0 : area_c = edgearea(i,j) - area1 - area2 - area3
2368 :
2369 : ! shift midpoint so that the area of remaining triangles = area_c
2370 0 : w1 = c2*area_c/areafac_c(i,j) &
2371 0 : + (xdr-xcl)*ydl + (xcr-xdl)*ydr
2372 0 : w2 = (xdr-xdl)**2 + (ydr-ydl)**2
2373 0 : w1 = w1/w2
2374 0 : xdm = xdm + (ydr - ydl) * w1
2375 0 : ydm = ydm - (xdr - xdl) * w1
2376 :
2377 : ! compute left and right intersection points
2378 0 : mdl = (ydm - ydl) / (xdm - xdl)
2379 0 : mdr = (ydr - ydm) / (xdr - xdm)
2380 :
2381 0 : if (abs(mdl) > puny) then
2382 0 : xicl = xdl - ydl/mdl
2383 : else
2384 0 : xicl = c0
2385 : endif
2386 0 : yicl = c0
2387 :
2388 0 : if (abs(mdr) > puny) then
2389 0 : xicr = xdr - ydr/mdr
2390 : else
2391 0 : xicr = c0
2392 : endif
2393 0 : yicr = c0
2394 :
2395 0 : elseif (xic < c0 .and. xic > xcl) then ! fix ICL = IC
2396 :
2397 0 : xicl = xic
2398 0 : yicl = yic
2399 :
2400 : ! compute midpoint between ICL and DR
2401 0 : xdm = p5 * (xdr + xicl)
2402 0 : ydm = p5 * ydr
2403 :
2404 : ! compute area of (lone) triangle adjacent to left corner
2405 0 : area4 = p5 * (xcl - xic) * ydl * areafac_c(i,j)
2406 0 : area_c = edgearea(i,j) - area1 - area2 - area3 - area4
2407 :
2408 : ! shift midpoint so that area of remaining triangles = area_c
2409 0 : w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr
2410 0 : w2 = (xdr-xic)**2 + ydr**2
2411 0 : w1 = w1/w2
2412 0 : xdm = xdm + ydr*w1
2413 0 : ydm = ydm - (xdr - xic) * w1
2414 :
2415 : ! compute ICR
2416 0 : mdr = (ydr - ydm) / (xdr - xdm)
2417 0 : if (abs(mdr) > puny) then
2418 0 : xicr = xdr - ydr/mdr
2419 : else
2420 0 : xicr = c0
2421 : endif
2422 0 : yicr = c0
2423 :
2424 0 : elseif (xic >= c0 .and. xic < xcr) then ! fix ICR = IR
2425 :
2426 0 : xicr = xic
2427 0 : yicr = yic
2428 :
2429 : ! compute midpoint between ICR and DL
2430 0 : xdm = p5 * (xicr + xdl)
2431 0 : ydm = p5 * ydl
2432 :
2433 : ! compute area of (lone) triangle adjacent to right corner
2434 0 : area4 = p5 * (xic - xcr) * ydr * areafac_c(i,j)
2435 0 : area_c = edgearea(i,j) - area1 - area2 - area3 - area4
2436 :
2437 : ! shift midpoint so that area of remaining triangles = area_c
2438 0 : w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl
2439 0 : w2 = (xic-xdl)**2 + ydl**2
2440 0 : w1 = w1/w2
2441 0 : xdm = xdm - ydl*w1
2442 0 : ydm = ydm - (xic - xdl) * w1
2443 :
2444 : ! compute ICL
2445 :
2446 0 : mdl = (ydm - ydl) / (xdm - xdl)
2447 0 : if (abs(mdl) > puny) then
2448 0 : xicl = xdl - ydl/mdl
2449 : else
2450 0 : xicl = c0
2451 : endif
2452 0 : yicl = c0
2453 :
2454 : endif ! ydl*ydr >= c0
2455 :
2456 : endif ! l_fixed_area
2457 :
2458 : !-------------------------------------------------------------------
2459 : ! Locate triangles in BC cell (H for both north and east edges)
2460 : ! and TC cell (N for north edge and E for east edge).
2461 : !-------------------------------------------------------------------
2462 :
2463 : ! Start with cases where both DPs lie in the same grid cell
2464 :
2465 13582521 : if (ydl >= c0 .and. ydr >= c0 .and. ydm >= c0) then
2466 :
2467 : ! TC1a (group 4)
2468 :
2469 1497825 : ng = 4
2470 1497825 : xp (i,j,1,ng) = xcl
2471 1497825 : yp (i,j,1,ng) = ycl
2472 1497825 : xp (i,j,2,ng) = xcr
2473 1497825 : yp (i,j,2,ng) = ycr
2474 1497825 : xp (i,j,3,ng) = xdl
2475 1497825 : yp (i,j,3,ng) = ydl
2476 1497825 : iflux (i,j,ng) = i + ishift_tc
2477 1497825 : jflux (i,j,ng) = j + jshift_tc
2478 1497825 : areafact(i,j,ng) = -areafac_c(i,j)
2479 :
2480 : ! TC2a (group 5)
2481 :
2482 1497825 : ng = 5
2483 1497825 : xp (i,j,1,ng) = xcr
2484 1497825 : yp (i,j,1,ng) = ycr
2485 1497825 : xp (i,j,2,ng) = xdr
2486 1497825 : yp (i,j,2,ng) = ydr
2487 1497825 : xp (i,j,3,ng) = xdl
2488 1497825 : yp (i,j,3,ng) = ydl
2489 1497825 : iflux (i,j,ng) = i + ishift_tc
2490 1497825 : jflux (i,j,ng) = j + jshift_tc
2491 1497825 : areafact(i,j,ng) = -areafac_c(i,j)
2492 :
2493 : ! TC3a (group 6)
2494 :
2495 1497825 : ng = 6
2496 1497825 : xp (i,j,1,ng) = xdl
2497 1497825 : yp (i,j,1,ng) = ydl
2498 1497825 : xp (i,j,2,ng) = xdr
2499 1497825 : yp (i,j,2,ng) = ydr
2500 1497825 : xp (i,j,3,ng) = xdm
2501 1497825 : yp (i,j,3,ng) = ydm
2502 1497825 : iflux (i,j,ng) = i + ishift_tc
2503 1497825 : jflux (i,j,ng) = j + jshift_tc
2504 1497825 : areafact(i,j,ng) = -areafac_c(i,j)
2505 :
2506 12038808 : elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then ! rare
2507 :
2508 : ! TC1b (group 4)
2509 :
2510 0 : ng = 4
2511 0 : xp (i,j,1,ng) = xcl
2512 0 : yp (i,j,1,ng) = ycl
2513 0 : xp (i,j,2,ng) = xicl
2514 0 : yp (i,j,2,ng) = yicl
2515 0 : xp (i,j,3,ng) = xdl
2516 0 : yp (i,j,3,ng) = ydl
2517 0 : iflux (i,j,ng) = i + ishift_tc
2518 0 : jflux (i,j,ng) = j + jshift_tc
2519 0 : areafact(i,j,ng) = -areafac_c(i,j)
2520 :
2521 : ! TC2b (group 5)
2522 :
2523 0 : ng = 5
2524 0 : xp (i,j,1,ng) = xcr
2525 0 : yp (i,j,1,ng) = ycr
2526 0 : xp (i,j,2,ng) = xdr
2527 0 : yp (i,j,2,ng) = ydr
2528 0 : xp (i,j,3,ng) = xicr
2529 0 : yp (i,j,3,ng) = yicr
2530 0 : iflux (i,j,ng) = i + ishift_tc
2531 0 : jflux (i,j,ng) = j + jshift_tc
2532 0 : areafact(i,j,ng) = -areafac_c(i,j)
2533 :
2534 : ! BC3b (group 6)
2535 :
2536 0 : ng = 6
2537 0 : xp (i,j,1,ng) = xicr
2538 0 : yp (i,j,1,ng) = yicr
2539 0 : xp (i,j,2,ng) = xicl
2540 0 : yp (i,j,2,ng) = yicl
2541 0 : xp (i,j,3,ng) = xdm
2542 0 : yp (i,j,3,ng) = ydm
2543 0 : iflux (i,j,ng) = i + ishift_bc
2544 0 : jflux (i,j,ng) = j + jshift_bc
2545 0 : areafact(i,j,ng) = areafac_c(i,j)
2546 :
2547 12038808 : elseif (ydl <= c0 .and. ydr <= c0 .and. ydm <= c0) then
2548 :
2549 : ! BC1a (group 4)
2550 :
2551 11777113 : ng = 4
2552 11777113 : xp (i,j,1,ng) = xcl
2553 11777113 : yp (i,j,1,ng) = ycl
2554 11777113 : xp (i,j,2,ng) = xdl
2555 11777113 : yp (i,j,2,ng) = ydl
2556 11777113 : xp (i,j,3,ng) = xcr
2557 11777113 : yp (i,j,3,ng) = ycr
2558 11777113 : iflux (i,j,ng) = i + ishift_bc
2559 11777113 : jflux (i,j,ng) = j + jshift_bc
2560 11777113 : areafact(i,j,ng) = areafac_c(i,j)
2561 :
2562 : ! BC2a (group 5)
2563 :
2564 11777113 : ng = 5
2565 11777113 : xp (i,j,1,ng) = xcr
2566 11777113 : yp (i,j,1,ng) = ycr
2567 11777113 : xp (i,j,2,ng) = xdl
2568 11777113 : yp (i,j,2,ng) = ydl
2569 11777113 : xp (i,j,3,ng) = xdr
2570 11777113 : yp (i,j,3,ng) = ydr
2571 11777113 : iflux (i,j,ng) = i + ishift_bc
2572 11777113 : jflux (i,j,ng) = j + jshift_bc
2573 11777113 : areafact(i,j,ng) = areafac_c(i,j)
2574 :
2575 : ! BC3a (group 6)
2576 :
2577 11777113 : ng = 6
2578 11777113 : xp (i,j,1,ng) = xdl
2579 11777113 : yp (i,j,1,ng) = ydl
2580 11777113 : xp (i,j,2,ng) = xdm
2581 11777113 : yp (i,j,2,ng) = ydm
2582 11777113 : xp (i,j,3,ng) = xdr
2583 11777113 : yp (i,j,3,ng) = ydr
2584 11777113 : iflux (i,j,ng) = i + ishift_bc
2585 11777113 : jflux (i,j,ng) = j + jshift_bc
2586 11777113 : areafact(i,j,ng) = areafac_c(i,j)
2587 :
2588 261695 : elseif (ydl <= c0 .and. ydr <= c0 .and. ydm > c0) then ! rare
2589 :
2590 : ! BC1b (group 4)
2591 :
2592 0 : ng = 4
2593 0 : xp (i,j,1,ng) = xcl
2594 0 : yp (i,j,1,ng) = ycl
2595 0 : xp (i,j,2,ng) = xdl
2596 0 : yp (i,j,2,ng) = ydl
2597 0 : xp (i,j,3,ng) = xicl
2598 0 : yp (i,j,3,ng) = yicl
2599 0 : iflux (i,j,ng) = i + ishift_bc
2600 0 : jflux (i,j,ng) = j + jshift_bc
2601 0 : areafact(i,j,ng) = areafac_c(i,j)
2602 :
2603 : ! BC2b (group 5)
2604 :
2605 0 : ng = 5
2606 0 : xp (i,j,1,ng) = xcr
2607 0 : yp (i,j,1,ng) = ycr
2608 0 : xp (i,j,2,ng) = xicr
2609 0 : yp (i,j,2,ng) = yicr
2610 0 : xp (i,j,3,ng) = xdr
2611 0 : yp (i,j,3,ng) = ydr
2612 0 : iflux (i,j,ng) = i + ishift_bc
2613 0 : jflux (i,j,ng) = j + jshift_bc
2614 0 : areafact(i,j,ng) = areafac_c(i,j)
2615 :
2616 : ! TC3b (group 6)
2617 :
2618 0 : ng = 6
2619 0 : xp (i,j,1,ng) = xicl
2620 0 : yp (i,j,1,ng) = yicl
2621 0 : xp (i,j,2,ng) = xicr
2622 0 : yp (i,j,2,ng) = yicr
2623 0 : xp (i,j,3,ng) = xdm
2624 0 : yp (i,j,3,ng) = ydm
2625 0 : iflux (i,j,ng) = i + ishift_tc
2626 0 : jflux (i,j,ng) = j + jshift_tc
2627 0 : areafact(i,j,ng) = -areafac_c(i,j)
2628 :
2629 : ! Now consider cases where the two DPs lie in different grid cells
2630 :
2631 : elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 &
2632 261695 : .and. ydm >= c0) then
2633 :
2634 : ! TC1b (group 4)
2635 :
2636 57424 : ng = 4
2637 57424 : xp (i,j,1,ng) = xcl
2638 57424 : yp (i,j,1,ng) = ycl
2639 57424 : xp (i,j,2,ng) = xicr
2640 57424 : yp (i,j,2,ng) = yicr
2641 57424 : xp (i,j,3,ng) = xdl
2642 57424 : yp (i,j,3,ng) = ydl
2643 57424 : iflux (i,j,ng) = i + ishift_tc
2644 57424 : jflux (i,j,ng) = j + jshift_tc
2645 57424 : areafact(i,j,ng) = -areafac_c(i,j)
2646 :
2647 : ! BC2b (group 5) lone triangle
2648 :
2649 57424 : ng = 5
2650 57424 : xp (i,j,1,ng) = xcr
2651 57424 : yp (i,j,1,ng) = ycr
2652 57424 : xp (i,j,2,ng) = xicr
2653 57424 : yp (i,j,2,ng) = yicr
2654 57424 : xp (i,j,3,ng) = xdr
2655 57424 : yp (i,j,3,ng) = ydr
2656 57424 : iflux (i,j,ng) = i + ishift_bc
2657 57424 : jflux (i,j,ng) = j + jshift_bc
2658 57424 : areafact(i,j,ng) = areafac_c(i,j)
2659 :
2660 : ! TC3b (group 6)
2661 :
2662 57424 : ng = 6
2663 57424 : xp (i,j,1,ng) = xdl
2664 57424 : yp (i,j,1,ng) = ydl
2665 57424 : xp (i,j,2,ng) = xicr
2666 57424 : yp (i,j,2,ng) = yicr
2667 57424 : xp (i,j,3,ng) = xdm
2668 57424 : yp (i,j,3,ng) = ydm
2669 57424 : iflux (i,j,ng) = i + ishift_tc
2670 57424 : jflux (i,j,ng) = j + jshift_tc
2671 57424 : areafact(i,j,ng) = -areafac_c(i,j)
2672 :
2673 : elseif (ydl > c0 .and. ydr < c0 .and. xic >= c0 &
2674 204271 : .and. ydm < c0 ) then ! less common
2675 :
2676 : ! TC1b (group 4)
2677 :
2678 94 : ng = 4
2679 94 : xp (i,j,1,ng) = xcl
2680 94 : yp (i,j,1,ng) = ycl
2681 94 : xp (i,j,2,ng) = xicl
2682 94 : yp (i,j,2,ng) = yicl
2683 94 : xp (i,j,3,ng) = xdl
2684 94 : yp (i,j,3,ng) = ydl
2685 94 : iflux (i,j,ng) = i + ishift_tc
2686 94 : jflux (i,j,ng) = j + jshift_tc
2687 94 : areafact(i,j,ng) = -areafac_c(i,j)
2688 :
2689 : ! BC2b (group 5) lone triangle
2690 :
2691 94 : ng = 5
2692 94 : xp (i,j,1,ng) = xcr
2693 94 : yp (i,j,1,ng) = ycr
2694 94 : xp (i,j,2,ng) = xicr
2695 94 : yp (i,j,2,ng) = yicr
2696 94 : xp (i,j,3,ng) = xdr
2697 94 : yp (i,j,3,ng) = ydr
2698 94 : iflux (i,j,ng) = i + ishift_bc
2699 94 : jflux (i,j,ng) = j + jshift_bc
2700 94 : areafact(i,j,ng) = areafac_c(i,j)
2701 :
2702 : ! BC3b (group 6)
2703 :
2704 94 : ng = 6
2705 94 : xp (i,j,1,ng) = xicr
2706 94 : yp (i,j,1,ng) = yicr
2707 94 : xp (i,j,2,ng) = xicl
2708 94 : yp (i,j,2,ng) = yicl
2709 94 : xp (i,j,3,ng) = xdm
2710 94 : yp (i,j,3,ng) = ydm
2711 94 : iflux (i,j,ng) = i + ishift_bc
2712 94 : jflux (i,j,ng) = j + jshift_bc
2713 94 : areafact(i,j,ng) = areafac_c(i,j)
2714 :
2715 : elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 &
2716 204177 : .and. ydm < c0) then
2717 :
2718 : ! TC1b (group 4) lone triangle
2719 :
2720 51291 : ng = 4
2721 51291 : xp (i,j,1,ng) = xcl
2722 51291 : yp (i,j,1,ng) = ycl
2723 51291 : xp (i,j,2,ng) = xicl
2724 51291 : yp (i,j,2,ng) = yicl
2725 51291 : xp (i,j,3,ng) = xdl
2726 51291 : yp (i,j,3,ng) = ydl
2727 51291 : iflux (i,j,ng) = i + ishift_tc
2728 51291 : jflux (i,j,ng) = j + jshift_tc
2729 51291 : areafact(i,j,ng) = -areafac_c(i,j)
2730 :
2731 : ! BC2b (group 5)
2732 :
2733 51291 : ng = 5
2734 51291 : xp (i,j,1,ng) = xcr
2735 51291 : yp (i,j,1,ng) = ycr
2736 51291 : xp (i,j,2,ng) = xicl
2737 51291 : yp (i,j,2,ng) = yicl
2738 51291 : xp (i,j,3,ng) = xdr
2739 51291 : yp (i,j,3,ng) = ydr
2740 51291 : iflux (i,j,ng) = i + ishift_bc
2741 51291 : jflux (i,j,ng) = j + jshift_bc
2742 51291 : areafact(i,j,ng) = areafac_c(i,j)
2743 :
2744 : ! BC3b (group 6)
2745 :
2746 51291 : ng = 6
2747 51291 : xp (i,j,1,ng) = xdr
2748 51291 : yp (i,j,1,ng) = ydr
2749 51291 : xp (i,j,2,ng) = xicl
2750 51291 : yp (i,j,2,ng) = yicl
2751 51291 : xp (i,j,3,ng) = xdm
2752 51291 : yp (i,j,3,ng) = ydm
2753 51291 : iflux (i,j,ng) = i + ishift_bc
2754 51291 : jflux (i,j,ng) = j + jshift_bc
2755 51291 : areafact(i,j,ng) = areafac_c(i,j)
2756 :
2757 : elseif (ydl > c0 .and. ydr < c0 .and. xic < c0 &
2758 152886 : .and. ydm >= c0) then ! less common
2759 :
2760 : ! TC1b (group 4) lone triangle
2761 :
2762 138 : ng = 4
2763 138 : xp (i,j,1,ng) = xcl
2764 138 : yp (i,j,1,ng) = ycl
2765 138 : xp (i,j,2,ng) = xicl
2766 138 : yp (i,j,2,ng) = yicl
2767 138 : xp (i,j,3,ng) = xdl
2768 138 : yp (i,j,3,ng) = ydl
2769 138 : iflux (i,j,ng) = i + ishift_tc
2770 138 : jflux (i,j,ng) = j + jshift_tc
2771 138 : areafact(i,j,ng) = -areafac_c(i,j)
2772 :
2773 : ! BC2b (group 5)
2774 :
2775 138 : ng = 5
2776 138 : xp (i,j,1,ng) = xcr
2777 138 : yp (i,j,1,ng) = ycr
2778 138 : xp (i,j,2,ng) = xicr
2779 138 : yp (i,j,2,ng) = yicr
2780 138 : xp (i,j,3,ng) = xdr
2781 138 : yp (i,j,3,ng) = ydr
2782 138 : iflux (i,j,ng) = i + ishift_bc
2783 138 : jflux (i,j,ng) = j + jshift_bc
2784 138 : areafact(i,j,ng) = areafac_c(i,j)
2785 :
2786 : ! TC3b (group 6)
2787 :
2788 138 : ng = 6
2789 138 : xp (i,j,1,ng) = xicl
2790 138 : yp (i,j,1,ng) = yicl
2791 138 : xp (i,j,2,ng) = xicr
2792 138 : yp (i,j,2,ng) = yicr
2793 138 : xp (i,j,3,ng) = xdm
2794 138 : yp (i,j,3,ng) = ydm
2795 138 : iflux (i,j,ng) = i + ishift_tc
2796 138 : jflux (i,j,ng) = j + jshift_tc
2797 138 : areafact(i,j,ng) = -areafac_c(i,j)
2798 :
2799 : elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 &
2800 152748 : .and. ydm >= c0) then
2801 :
2802 : ! BC1b (group 4) lone triangle
2803 :
2804 74220 : ng = 4
2805 74220 : xp (i,j,1,ng) = xcl
2806 74220 : yp (i,j,1,ng) = ycl
2807 74220 : xp (i,j,2,ng) = xdl
2808 74220 : yp (i,j,2,ng) = ydl
2809 74220 : xp (i,j,3,ng) = xicl
2810 74220 : yp (i,j,3,ng) = yicl
2811 74220 : iflux (i,j,ng) = i + ishift_bc
2812 74220 : jflux (i,j,ng) = j + jshift_bc
2813 74220 : areafact(i,j,ng) = areafac_c(i,j)
2814 :
2815 : ! TC2b (group 5)
2816 :
2817 74220 : ng = 5
2818 74220 : xp (i,j,1,ng) = xcr
2819 74220 : yp (i,j,1,ng) = ycr
2820 74220 : xp (i,j,2,ng) = xdr
2821 74220 : yp (i,j,2,ng) = ydr
2822 74220 : xp (i,j,3,ng) = xicl
2823 74220 : yp (i,j,3,ng) = yicl
2824 74220 : iflux (i,j,ng) = i + ishift_tc
2825 74220 : jflux (i,j,ng) = j + jshift_tc
2826 74220 : areafact(i,j,ng) = -areafac_c(i,j)
2827 :
2828 : ! TC3b (group 6)
2829 :
2830 74220 : ng = 6
2831 74220 : xp (i,j,1,ng) = xicl
2832 74220 : yp (i,j,1,ng) = yicl
2833 74220 : xp (i,j,2,ng) = xdr
2834 74220 : yp (i,j,2,ng) = ydr
2835 74220 : xp (i,j,3,ng) = xdm
2836 74220 : yp (i,j,3,ng) = ydm
2837 74220 : iflux (i,j,ng) = i + ishift_tc
2838 74220 : jflux (i,j,ng) = j + jshift_tc
2839 74220 : areafact(i,j,ng) = -areafac_c(i,j)
2840 :
2841 : elseif (ydl < c0 .and. ydr > c0 .and. xic < c0 &
2842 78528 : .and. ydm < c0) then ! less common
2843 :
2844 : ! BC1b (group 4) lone triangle
2845 :
2846 271 : ng = 4
2847 271 : xp (i,j,1,ng) = xcl
2848 271 : yp (i,j,1,ng) = ycl
2849 271 : xp (i,j,2,ng) = xdl
2850 271 : yp (i,j,2,ng) = ydl
2851 271 : xp (i,j,3,ng) = xicl
2852 271 : yp (i,j,3,ng) = yicl
2853 271 : iflux (i,j,ng) = i + ishift_bc
2854 271 : jflux (i,j,ng) = j + jshift_bc
2855 271 : areafact(i,j,ng) = areafac_c(i,j)
2856 :
2857 : ! TC2b (group 5)
2858 :
2859 271 : ng = 5
2860 271 : xp (i,j,1,ng) = xcr
2861 271 : yp (i,j,1,ng) = ycr
2862 271 : xp (i,j,2,ng) = xdr
2863 271 : yp (i,j,2,ng) = ydr
2864 271 : xp (i,j,3,ng) = xicr
2865 271 : yp (i,j,3,ng) = yicr
2866 271 : iflux (i,j,ng) = i + ishift_tc
2867 271 : jflux (i,j,ng) = j + jshift_tc
2868 271 : areafact(i,j,ng) = -areafac_c(i,j)
2869 :
2870 : ! BC3b (group 6)
2871 :
2872 271 : ng = 6
2873 271 : xp (i,j,1,ng) = xicr
2874 271 : yp (i,j,1,ng) = yicr
2875 271 : xp (i,j,2,ng) = xicl
2876 271 : yp (i,j,2,ng) = yicl
2877 271 : xp (i,j,3,ng) = xdm
2878 271 : yp (i,j,3,ng) = ydm
2879 271 : iflux (i,j,ng) = i + ishift_bc
2880 271 : jflux (i,j,ng) = j + jshift_bc
2881 271 : areafact(i,j,ng) = areafac_c(i,j)
2882 :
2883 : elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 &
2884 78257 : .and. ydm < c0) then
2885 :
2886 : ! BC1b (group 4)
2887 :
2888 78115 : ng = 4
2889 78115 : xp (i,j,1,ng) = xcl
2890 78115 : yp (i,j,1,ng) = ycl
2891 78115 : xp (i,j,2,ng) = xdl
2892 78115 : yp (i,j,2,ng) = ydl
2893 78115 : xp (i,j,3,ng) = xicr
2894 78115 : yp (i,j,3,ng) = yicr
2895 78115 : iflux (i,j,ng) = i + ishift_bc
2896 78115 : jflux (i,j,ng) = j + jshift_bc
2897 78115 : areafact(i,j,ng) = areafac_c(i,j)
2898 :
2899 : ! TC2b (group 5) lone triangle
2900 :
2901 78115 : ng = 5
2902 78115 : xp (i,j,1,ng) = xcr
2903 78115 : yp (i,j,1,ng) = ycr
2904 78115 : xp (i,j,2,ng) = xdr
2905 78115 : yp (i,j,2,ng) = ydr
2906 78115 : xp (i,j,3,ng) = xicr
2907 78115 : yp (i,j,3,ng) = yicr
2908 78115 : iflux (i,j,ng) = i + ishift_tc
2909 78115 : jflux (i,j,ng) = j + jshift_tc
2910 78115 : areafact(i,j,ng) = -areafac_c(i,j)
2911 :
2912 : ! BC3b (group 6)
2913 :
2914 78115 : ng = 6
2915 78115 : xp (i,j,1,ng) = xicr
2916 78115 : yp (i,j,1,ng) = yicr
2917 78115 : xp (i,j,2,ng) = xdl
2918 78115 : yp (i,j,2,ng) = ydl
2919 78115 : xp (i,j,3,ng) = xdm
2920 78115 : yp (i,j,3,ng) = ydm
2921 78115 : iflux (i,j,ng) = i + ishift_bc
2922 78115 : jflux (i,j,ng) = j + jshift_bc
2923 78115 : areafact(i,j,ng) = areafac_c(i,j)
2924 :
2925 : elseif (ydl < c0 .and. ydr > c0 .and. xic >= c0 &
2926 142 : .and. ydm >= c0) then ! less common
2927 :
2928 : ! BC1b (group 4)
2929 :
2930 142 : ng = 4
2931 142 : xp (i,j,1,ng) = xcl
2932 142 : yp (i,j,1,ng) = ycl
2933 142 : xp (i,j,2,ng) = xdl
2934 142 : yp (i,j,2,ng) = ydl
2935 142 : xp (i,j,3,ng) = xicl
2936 142 : yp (i,j,3,ng) = yicl
2937 142 : iflux (i,j,ng) = i + ishift_bc
2938 142 : jflux (i,j,ng) = j + jshift_bc
2939 142 : areafact(i,j,ng) = areafac_c(i,j)
2940 :
2941 : ! TC2b (group 5) lone triangle
2942 :
2943 142 : ng = 5
2944 142 : xp (i,j,1,ng) = xcr
2945 142 : yp (i,j,1,ng) = ycr
2946 142 : xp (i,j,2,ng) = xdr
2947 142 : yp (i,j,2,ng) = ydr
2948 142 : xp (i,j,3,ng) = xicr
2949 142 : yp (i,j,3,ng) = yicr
2950 142 : iflux (i,j,ng) = i + ishift_tc
2951 142 : jflux (i,j,ng) = j + jshift_tc
2952 142 : areafact(i,j,ng) = -areafac_c(i,j)
2953 :
2954 : ! TC3b (group 6)
2955 :
2956 142 : ng = 6
2957 142 : xp (i,j,1,ng) = xicl
2958 142 : yp (i,j,1,ng) = yicl
2959 142 : xp (i,j,2,ng) = xicr
2960 142 : yp (i,j,2,ng) = yicr
2961 142 : xp (i,j,3,ng) = xdm
2962 142 : yp (i,j,3,ng) = ydm
2963 142 : iflux (i,j,ng) = i + ishift_tc
2964 142 : jflux (i,j,ng) = j + jshift_tc
2965 142 : areafact(i,j,ng) = -areafac_c(i,j)
2966 :
2967 : endif ! TC and BC triangles
2968 :
2969 : enddo ! ij
2970 :
2971 : !-------------------------------------------------------------------
2972 : ! Compute triangle areas with appropriate sign.
2973 : ! These are found by computing the area in scaled coordinates and
2974 : ! multiplying by a scale factor (areafact).
2975 : ! Note that the scale factor is positive for fluxes out of the cell
2976 : ! and negative for fluxes into the cell.
2977 : !
2978 : ! Note: The triangle area formula below gives A >=0 iff the triangle
2979 : ! points x1, x2, and x3 are taken in counterclockwise order.
2980 : ! These points are defined above in such a way that the
2981 : ! order is nearly always CCW.
2982 : ! In rare cases, we may compute A < 0. In this case,
2983 : ! the quadrilateral departure area is equal to the
2984 : ! difference of two triangle areas instead of the sum.
2985 : ! The fluxes work out correctly in the end.
2986 : !
2987 : ! Also compute the cumulative area transported across each edge.
2988 : ! If l_fixed_area = T, this area is compared to edgearea as a bug check.
2989 : ! If l_fixed_area = F, this area is passed as an output array.
2990 : !-------------------------------------------------------------------
2991 :
2992 32372640 : areasum(:,:) = c0
2993 :
2994 321216 : do ng = 1, ngroups
2995 275328 : icells(ng) = 0
2996 :
2997 81541014 : do ij = 1, icellsd
2998 81219798 : i = indxid(ij)
2999 81219798 : j = indxjd(ij)
3000 :
3001 6488196 : triarea(i,j,ng) = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
3002 : (yp(i,j,3,ng)-yp(i,j,1,ng)) & ! LCOV_EXCL_LINE
3003 : - (yp(i,j,2,ng)-yp(i,j,1,ng)) * & ! LCOV_EXCL_LINE
3004 : (xp(i,j,3,ng)-xp(i,j,1,ng)) ) & ! LCOV_EXCL_LINE
3005 81219798 : * areafact(i,j,ng)
3006 :
3007 81219798 : if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then
3008 41522255 : triarea(i,j,ng) = c0
3009 : else
3010 39697543 : icells(ng) = icells(ng) + 1
3011 39697543 : ic = icells(ng)
3012 39697543 : indxi(ic,ng) = i
3013 39697543 : indxj(ic,ng) = j
3014 : endif
3015 :
3016 81495126 : areasum(i,j) = areasum(i,j) + triarea(i,j,ng)
3017 :
3018 : enddo ! ij
3019 : enddo ! ng
3020 :
3021 45888 : if (l_fixed_area) then
3022 : if (bugcheck) then ! set bugcheck = F to speed up code
3023 : do ij = 1, icellsd
3024 : i = indxid(ij)
3025 : j = indxjd(ij)
3026 : if ( abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j) .and. abs(edgearea(i,j)) > c0 ) then
3027 : write(nu_diag,*) ''
3028 : write(nu_diag,*) 'Areas do not add up: m, i, j, edge =', &
3029 : my_task, i, j, trim(edge)
3030 : write(nu_diag,*) 'edgearea =', edgearea(i,j)
3031 : write(nu_diag,*) 'areasum =', areasum(i,j)
3032 : write(nu_diag,*) 'areafac_c =', areafac_c(i,j)
3033 : write(nu_diag,*) ''
3034 : write(nu_diag,*) 'Triangle areas:'
3035 : do ng = 1, ngroups ! not vector friendly
3036 : if (abs(triarea(i,j,ng)) > eps16*abs(areafact(i,j,ng))) then
3037 : write(nu_diag,*) ng, triarea(i,j,ng)
3038 : endif
3039 : enddo
3040 : endif
3041 : enddo
3042 : endif ! bugcheck
3043 :
3044 : else ! l_fixed_area = F
3045 13582521 : do ij = 1, icellsd
3046 13536633 : i = indxid(ij)
3047 13536633 : j = indxjd(ij)
3048 13582521 : edgearea(i,j) = areasum(i,j)
3049 : enddo
3050 : endif ! l_fixed_area
3051 :
3052 : !-------------------------------------------------------------------
3053 : ! Transform triangle vertices to a scaled coordinate system centered
3054 : ! in the cell containing the triangle.
3055 : !-------------------------------------------------------------------
3056 :
3057 45888 : if (trim(edge) == 'north') then
3058 160608 : do ng = 1, ngroups
3059 573600 : do nv = 1, nvert
3060 59849262 : do ij = 1, icells(ng)
3061 59298606 : i = indxi(ij,ng)
3062 59298606 : j = indxj(ij,ng)
3063 59298606 : ishift = iflux(i,j,ng) - i
3064 59298606 : jshift = jflux(i,j,ng) - j
3065 59298606 : xp(i,j,nv,ng) = xp(i,j,nv,ng) - c1*ishift
3066 59711598 : yp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*jshift
3067 : enddo ! ij
3068 : enddo ! nv
3069 : enddo ! ng
3070 : else ! east edge
3071 160608 : do ng = 1, ngroups
3072 573600 : do nv = 1, nvert
3073 60344679 : do ij = 1, icells(ng)
3074 59794023 : i = indxi(ij,ng)
3075 59794023 : j = indxj(ij,ng)
3076 59794023 : ishift = iflux(i,j,ng) - i
3077 59794023 : jshift = jflux(i,j,ng) - j
3078 : ! Note rotation of pi/2 here
3079 59794023 : w1 = xp(i,j,nv,ng)
3080 59794023 : xp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*ishift
3081 60207015 : yp(i,j,nv,ng) = -w1 - c1*jshift
3082 : enddo ! ij
3083 : enddo ! nv
3084 : enddo ! ng
3085 : endif
3086 :
3087 : if (bugcheck) then
3088 : if (trim(edge) == 'north') then
3089 : ib = ilo
3090 : jb = jlo-1
3091 : else ! east edge
3092 : ib = ilo-1
3093 : jb = jlo
3094 : endif
3095 : do ng = 1, ngroups
3096 : do nv = 1, nvert
3097 : do j = jb, jhi
3098 : do i = ib, ihi
3099 : if (abs(triarea(i,j,ng)) > puny) then
3100 : if (abs(xp(i,j,nv,ng)) > p5+puny) then
3101 : write(nu_diag,*) ''
3102 : write(nu_diag,*) 'WARNING: xp =', xp(i,j,nv,ng)
3103 : write(nu_diag,*) 'm, i, j, ng, nv =', my_task, i, j, ng, nv
3104 : ! write(nu_diag,*) 'yil,xdl,xcl,ydl=',yil,xdl,xcl,ydl
3105 : ! write(nu_diag,*) 'yir,xdr,xcr,ydr=',yir,xdr,xcr,ydr
3106 : ! write(nu_diag,*) 'ydm=',ydm
3107 : ! stop
3108 : endif
3109 : if (abs(yp(i,j,nv,ng)) > p5+puny) then
3110 : write(nu_diag,*) ''
3111 : write(nu_diag,*) 'WARNING: yp =', yp(i,j,nv,ng)
3112 : write(nu_diag,*) 'm, i, j, ng, nv =', my_task, i, j, ng, nv
3113 : endif
3114 : endif ! triarea
3115 : enddo
3116 : enddo
3117 : enddo
3118 : enddo
3119 : endif ! bugcheck
3120 :
3121 45888 : end subroutine locate_triangles
3122 :
3123 : !=======================================================================
3124 : !
3125 : ! For each triangle, find the coordinates of the quadrature points needed
3126 : ! to compute integrals of linear, quadratic, or cubic polynomials,
3127 : ! using formulas from A.H. Stroud, Approximate Calculation of Multiple
3128 : ! Integrals, Prentice-Hall, 1971. (Section 8.8, formula 3.1.)
3129 : ! Linear functions can be integrated exactly by evaluating the function
3130 : ! at just one point (the midpoint). Quadratic functions require
3131 : ! 3 points, and cubics require 4 points.
3132 : ! The default is cubic, but the code can be sped up slightly using
3133 : ! linear or quadratic integrals, usually with little loss of accuracy.
3134 : !
3135 : ! The formulas are as follows:
3136 : !
3137 : ! I1 = integral of f(x,y)*dA
3138 : ! = A * f(x0,y0)
3139 : ! where A is the traingle area and (x0,y0) is the midpoint.
3140 : !
3141 : ! I2 = A * (f(x1,y1) + f(x2,y2) + f(x3,y3))
3142 : ! where these three points are located halfway between the midpoint
3143 : ! and the three vertics of the triangle.
3144 : !
3145 : ! I3 = A * [ -9/16 * f(x0,y0)
3146 : ! + 25/48 * (f(x1,y1) + f(x2,y2) + f(x3,y3))]
3147 : ! where (x0,y0) is the midpoint, and the other three points are
3148 : ! located 2/5 of the way from the midpoint to the three vertices.
3149 : !
3150 : ! author William H. Lipscomb, LANL
3151 :
3152 45888 : subroutine triangle_coordinates (nx_block, ny_block, &
3153 : integral_order, icells, & ! LCOV_EXCL_LINE
3154 : indxi, indxj, & ! LCOV_EXCL_LINE
3155 45888 : xp, yp)
3156 :
3157 : integer (kind=int_kind), intent(in) :: &
3158 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
3159 : integral_order ! polynomial order for quadrature integrals
3160 :
3161 : integer (kind=int_kind), dimension (ngroups), intent(in) :: &
3162 : icells ! number of cells where triarea > puny
3163 :
3164 : integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), intent(in) :: &
3165 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
3166 : indxj ! compressed index in j-direction
3167 :
3168 : real (kind=dbl_kind), intent(inout), dimension (nx_block, ny_block, 0:nvert, ngroups) :: &
3169 : xp, yp ! coordinates of triangle points
3170 :
3171 : ! local variables
3172 :
3173 : integer (kind=int_kind) :: &
3174 : i, j, ij , & ! horizontal indices ! LCOV_EXCL_LINE
3175 : ng ! triangle index
3176 :
3177 : character(len=*), parameter :: subname = '(triangle_coordinates)'
3178 :
3179 45888 : if (integral_order == 1) then ! linear (1-point formula)
3180 :
3181 0 : do ng = 1, ngroups
3182 0 : do ij = 1, icells(ng)
3183 0 : i = indxi(ij,ng)
3184 0 : j = indxj(ij,ng)
3185 :
3186 : ! coordinates of midpoint
3187 0 : xp(i,j,0,ng) = p333 &
3188 0 : * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
3189 0 : yp(i,j,0,ng) = p333 &
3190 0 : * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
3191 :
3192 : enddo ! ij
3193 : enddo ! ng
3194 :
3195 45888 : elseif (integral_order == 2) then ! quadratic (3-point formula)
3196 :
3197 0 : do ng = 1, ngroups
3198 0 : do ij = 1, icells(ng)
3199 0 : i = indxi(ij,ng)
3200 0 : j = indxj(ij,ng)
3201 :
3202 : ! coordinates of midpoint
3203 0 : xp(i,j,0,ng) = p333 &
3204 0 : * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
3205 0 : yp(i,j,0,ng) = p333 &
3206 0 : * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
3207 :
3208 : ! coordinates of the 3 points needed for integrals
3209 :
3210 0 : xp(i,j,1,ng) = p5*xp(i,j,1,ng) + p5*xp(i,j,0,ng)
3211 0 : yp(i,j,1,ng) = p5*yp(i,j,1,ng) + p5*yp(i,j,0,ng)
3212 :
3213 0 : xp(i,j,2,ng) = p5*xp(i,j,2,ng) + p5*xp(i,j,0,ng)
3214 0 : yp(i,j,2,ng) = p5*yp(i,j,2,ng) + p5*yp(i,j,0,ng)
3215 :
3216 0 : xp(i,j,3,ng) = p5*xp(i,j,3,ng) + p5*xp(i,j,0,ng)
3217 0 : yp(i,j,3,ng) = p5*yp(i,j,3,ng) + p5*yp(i,j,0,ng)
3218 :
3219 : enddo ! ij
3220 : enddo ! ng
3221 :
3222 : else ! cubic (4-point formula)
3223 :
3224 321216 : do ng = 1, ngroups
3225 40018759 : do ij = 1, icells(ng)
3226 39697543 : i = indxi(ij,ng)
3227 39697543 : j = indxj(ij,ng)
3228 :
3229 : ! coordinates of midpoint
3230 2884556 : xp(i,j,0,ng) = p333 &
3231 39697543 : * (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
3232 2884556 : yp(i,j,0,ng) = p333 &
3233 39697543 : * (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
3234 :
3235 : ! coordinates of the other 3 points needed for integrals
3236 :
3237 39697543 : xp(i,j,1,ng) = p4*xp(i,j,1,ng) + p6*xp(i,j,0,ng)
3238 39697543 : yp(i,j,1,ng) = p4*yp(i,j,1,ng) + p6*yp(i,j,0,ng)
3239 :
3240 39697543 : xp(i,j,2,ng) = p4*xp(i,j,2,ng) + p6*xp(i,j,0,ng)
3241 39697543 : yp(i,j,2,ng) = p4*yp(i,j,2,ng) + p6*yp(i,j,0,ng)
3242 :
3243 39697543 : xp(i,j,3,ng) = p4*xp(i,j,3,ng) + p6*xp(i,j,0,ng)
3244 39972871 : yp(i,j,3,ng) = p4*yp(i,j,3,ng) + p6*yp(i,j,0,ng)
3245 :
3246 : enddo ! ij
3247 : enddo ! ng
3248 :
3249 : endif
3250 :
3251 45888 : end subroutine triangle_coordinates
3252 :
3253 : !=======================================================================
3254 : !
3255 : ! Compute the transports across each edge by integrating the mass
3256 : ! and tracers over each departure triangle.
3257 : ! Input variables have the same meanings as in the main subroutine.
3258 : ! Repeated use of certain sums makes the calculation more efficient.
3259 : ! Integral formulas are described in triangle_coordinates subroutine.
3260 : !
3261 : ! author William H. Lipscomb, LANL
3262 :
3263 275328 : subroutine transport_integrals (nx_block, ny_block, &
3264 : ntrace, icells, & ! LCOV_EXCL_LINE
3265 : indxi, indxj, & ! LCOV_EXCL_LINE
3266 : tracer_type, depend, & ! LCOV_EXCL_LINE
3267 : integral_order, triarea, & ! LCOV_EXCL_LINE
3268 : iflux, jflux, & ! LCOV_EXCL_LINE
3269 : xp, yp, & ! LCOV_EXCL_LINE
3270 : mc, mx, & ! LCOV_EXCL_LINE
3271 : my, mflx, & ! LCOV_EXCL_LINE
3272 : tc, tx, & ! LCOV_EXCL_LINE
3273 458880 : ty, mtflx)
3274 :
3275 : integer (kind=int_kind), intent(in) :: &
3276 : nx_block, ny_block , & ! block dimensions ! LCOV_EXCL_LINE
3277 : ntrace , & ! number of tracers in use ! LCOV_EXCL_LINE
3278 : integral_order ! polynomial order for quadrature integrals
3279 :
3280 : integer (kind=int_kind), dimension (ntrace), intent(in) :: &
3281 : tracer_type , & ! = 1, 2, or 3 (see comments above) ! LCOV_EXCL_LINE
3282 : depend ! tracer dependencies (see above)
3283 :
3284 : integer (kind=int_kind), dimension (ngroups), intent(in) :: &
3285 : icells ! number of cells where triarea > puny
3286 :
3287 : integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), intent(in) :: &
3288 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
3289 : indxj ! compressed index in j-direction
3290 :
3291 : real (kind=dbl_kind), intent(in), dimension (nx_block, ny_block, 0:nvert, ngroups) :: &
3292 : xp, yp ! coordinates of triangle points
3293 :
3294 : real (kind=dbl_kind), intent(in), dimension (nx_block, ny_block, ngroups) :: &
3295 : triarea ! triangle area
3296 :
3297 : integer (kind=int_kind), intent(in), dimension (nx_block, ny_block, ngroups) :: &
3298 : iflux ,& ! LCOV_EXCL_LINE
3299 : jflux
3300 :
3301 : real (kind=dbl_kind), intent(in), dimension (nx_block, ny_block) :: &
3302 : mc, mx, my
3303 :
3304 : real (kind=dbl_kind), intent(out), dimension (nx_block, ny_block) :: &
3305 : mflx
3306 :
3307 : real (kind=dbl_kind), intent(in), dimension (nx_block, ny_block, ntrace), optional :: &
3308 : tc, tx, ty
3309 :
3310 : real (kind=dbl_kind), intent(out), dimension (nx_block, ny_block, ntrace), optional :: &
3311 : mtflx
3312 :
3313 : ! local variables
3314 :
3315 : integer (kind=int_kind) :: &
3316 : i, j, ij , & ! horizontal indices of edge ! LCOV_EXCL_LINE
3317 : i2, j2 , & ! horizontal indices of cell contributing transport ! LCOV_EXCL_LINE
3318 : ng , & ! triangle index ! LCOV_EXCL_LINE
3319 : nt, nt1 ! tracer indices
3320 :
3321 : real (kind=dbl_kind) :: &
3322 : m0, m1, m2, m3 , & ! mass field at internal points ! LCOV_EXCL_LINE
3323 69120 : w0, w1, w2, w3 ! work variables
3324 :
3325 : real (kind=dbl_kind), dimension (nx_block, ny_block) :: &
3326 : msum, mxsum, mysum , & ! sum of mass, mass*x, and mass*y ! LCOV_EXCL_LINE
3327 180539136 : mxxsum, mxysum, myysum ! sum of mass*x*x, mass*x*y, mass*y*y
3328 :
3329 : real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace) :: &
3330 : mtsum , & ! sum of mass*tracer ! LCOV_EXCL_LINE
3331 : mtxsum , & ! sum of mass*tracer*x ! LCOV_EXCL_LINE
3332 1562247936 : mtysum ! sum of mass*tracer*y
3333 :
3334 : character(len=*), parameter :: subname = '(transport_integrals)'
3335 :
3336 : !-------------------------------------------------------------------
3337 : ! Initialize
3338 : !-------------------------------------------------------------------
3339 :
3340 194235840 : mflx(:,:) = c0
3341 275328 : if (present(mtflx)) then
3342 6194880 : do nt = 1, ntrace
3343 4208672640 : mtflx(:,:,nt) = c0
3344 : enddo
3345 : endif
3346 :
3347 : !-------------------------------------------------------------------
3348 : ! Main loop
3349 : !-------------------------------------------------------------------
3350 :
3351 1927296 : do ng = 1, ngroups
3352 :
3353 1651968 : if (integral_order == 1) then ! linear (1-point formula)
3354 :
3355 0 : do ij = 1, icells(ng)
3356 0 : i = indxi(ij,ng)
3357 0 : j = indxj(ij,ng)
3358 :
3359 0 : i2 = iflux(i,j,ng)
3360 0 : j2 = jflux(i,j,ng)
3361 :
3362 : ! mass transports
3363 :
3364 0 : m0 = mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
3365 0 : + yp(i,j,0,ng)*my(i2,j2)
3366 0 : msum(i,j) = m0
3367 :
3368 0 : mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
3369 :
3370 : ! quantities needed for tracer transports
3371 0 : mxsum(i,j) = m0*xp(i,j,0,ng)
3372 0 : mxxsum(i,j) = mxsum(i,j)*xp(i,j,0,ng)
3373 0 : mxysum(i,j) = mxsum(i,j)*yp(i,j,0,ng)
3374 0 : mysum(i,j) = m0*yp(i,j,0,ng)
3375 0 : myysum(i,j) = mysum(i,j)*yp(i,j,0,ng)
3376 : enddo ! ij
3377 :
3378 1651968 : elseif (integral_order == 2) then ! quadratic (3-point formula)
3379 :
3380 0 : do ij = 1, icells(ng)
3381 0 : i = indxi(ij,ng)
3382 0 : j = indxj(ij,ng)
3383 :
3384 0 : i2 = iflux(i,j,ng)
3385 0 : j2 = jflux(i,j,ng)
3386 :
3387 : ! mass transports
3388 : ! Weighting factor of 1/3 is incorporated into the ice
3389 : ! area terms m1, m2, and m3.
3390 0 : m1 = p333 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
3391 0 : + yp(i,j,1,ng)*my(i2,j2))
3392 0 : m2 = p333 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
3393 0 : + yp(i,j,2,ng)*my(i2,j2))
3394 0 : m3 = p333 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
3395 0 : + yp(i,j,3,ng)*my(i2,j2))
3396 0 : msum(i,j) = m1 + m2 + m3
3397 0 : mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
3398 :
3399 : ! quantities needed for mass_tracer transports
3400 0 : w1 = m1 * xp(i,j,1,ng)
3401 0 : w2 = m2 * xp(i,j,2,ng)
3402 0 : w3 = m3 * xp(i,j,3,ng)
3403 :
3404 0 : mxsum(i,j) = w1 + w2 + w3
3405 :
3406 0 : mxxsum(i,j) = w1*xp(i,j,1,ng) + w2*xp(i,j,2,ng) &
3407 0 : + w3*xp(i,j,3,ng)
3408 :
3409 0 : mxysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
3410 0 : + w3*yp(i,j,3,ng)
3411 :
3412 0 : w1 = m1 * yp(i,j,1,ng)
3413 0 : w2 = m2 * yp(i,j,2,ng)
3414 0 : w3 = m3 * yp(i,j,3,ng)
3415 :
3416 0 : mysum(i,j) = w1 + w2 + w3
3417 :
3418 0 : myysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
3419 0 : + w3*yp(i,j,3,ng)
3420 : enddo ! ij
3421 :
3422 : else ! cubic (4-point formula)
3423 :
3424 239837226 : do ij = 1, icells(ng)
3425 238185258 : i = indxi(ij,ng)
3426 238185258 : j = indxj(ij,ng)
3427 :
3428 238185258 : i2 = iflux(i,j,ng)
3429 238185258 : j2 = jflux(i,j,ng)
3430 :
3431 : ! mass transports
3432 :
3433 : ! Weighting factors are incorporated into the
3434 : ! terms m0, m1, m2, and m3.
3435 17307336 : m0 = p5625m * (mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
3436 238185258 : + yp(i,j,0,ng)*my(i2,j2))
3437 17307336 : m1 = p52083 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
3438 238185258 : + yp(i,j,1,ng)*my(i2,j2))
3439 17307336 : m2 = p52083 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
3440 238185258 : + yp(i,j,2,ng)*my(i2,j2))
3441 17307336 : m3 = p52083 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
3442 238185258 : + yp(i,j,3,ng)*my(i2,j2))
3443 238185258 : msum(i,j) = m0 + m1 + m2 + m3
3444 238185258 : mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
3445 :
3446 : ! quantities needed for tracer transports
3447 238185258 : w0 = m0 * xp(i,j,0,ng)
3448 238185258 : w1 = m1 * xp(i,j,1,ng)
3449 238185258 : w2 = m2 * xp(i,j,2,ng)
3450 238185258 : w3 = m3 * xp(i,j,3,ng)
3451 :
3452 238185258 : mxsum(i,j) = w0 + w1 + w2 + w3
3453 :
3454 17307336 : mxxsum(i,j) = w0*xp(i,j,0,ng) + w1*xp(i,j,1,ng) &
3455 238185258 : + w2*xp(i,j,2,ng) + w3*xp(i,j,3,ng)
3456 :
3457 17307336 : mxysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
3458 238185258 : + w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
3459 :
3460 238185258 : w0 = m0 * yp(i,j,0,ng)
3461 238185258 : w1 = m1 * yp(i,j,1,ng)
3462 238185258 : w2 = m2 * yp(i,j,2,ng)
3463 238185258 : w3 = m3 * yp(i,j,3,ng)
3464 :
3465 238185258 : mysum(i,j) = w0 + w1 + w2 + w3
3466 :
3467 17307336 : myysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
3468 239837226 : + w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
3469 :
3470 : enddo ! ij
3471 :
3472 : endif ! integral_order
3473 :
3474 : ! mass * tracer transports
3475 :
3476 1927296 : if (present(mtflx)) then
3477 :
3478 37169280 : do nt = 1, ntrace
3479 37169280 : if (tracer_type(nt)==1) then ! does not depend on another tracer
3480 :
3481 1199186130 : do ij = 1, icells(ng)
3482 1190926290 : i = indxi(ij,ng)
3483 1190926290 : j = indxj(ij,ng)
3484 :
3485 1190926290 : i2 = iflux(i,j,ng)
3486 1190926290 : j2 = jflux(i,j,ng)
3487 :
3488 86536680 : mtsum(i,j,nt) = msum(i,j) * tc(i2,j2,nt) &
3489 : + mxsum(i,j) * tx(i2,j2,nt) & ! LCOV_EXCL_LINE
3490 1190926290 : + mysum(i,j) * ty(i2,j2,nt)
3491 :
3492 86536680 : mtflx(i,j,nt) = mtflx(i,j,nt) &
3493 1190926290 : + triarea(i,j,ng) * mtsum(i,j,nt)
3494 :
3495 : ! quantities needed for dependent tracers
3496 :
3497 86536680 : mtxsum(i,j,nt) = mxsum(i,j) * tc(i2,j2,nt) &
3498 : + mxxsum(i,j) * tx(i2,j2,nt) & ! LCOV_EXCL_LINE
3499 1190926290 : + mxysum(i,j) * ty(i2,j2,nt)
3500 :
3501 86536680 : mtysum(i,j,nt) = mysum(i,j) * tc(i2,j2,nt) &
3502 : + mxysum(i,j) * tx(i2,j2,nt) & ! LCOV_EXCL_LINE
3503 1199186130 : + myysum(i,j) * ty(i2,j2,nt)
3504 : enddo ! ij
3505 :
3506 27532800 : elseif (tracer_type(nt)==2) then ! depends on another tracer
3507 24779520 : nt1 = depend(nt)
3508 :
3509 3597558390 : do ij = 1, icells(ng)
3510 3572778870 : i = indxi(ij,ng)
3511 3572778870 : j = indxj(ij,ng)
3512 :
3513 3572778870 : i2 = iflux(i,j,ng)
3514 3572778870 : j2 = jflux(i,j,ng)
3515 :
3516 259610040 : mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt) &
3517 : + mtxsum(i,j,nt1) * tx(i2,j2,nt) & ! LCOV_EXCL_LINE
3518 3572778870 : + mtysum(i,j,nt1) * ty(i2,j2,nt)
3519 :
3520 259610040 : mtflx(i,j,nt) = mtflx(i,j,nt) &
3521 3597558390 : + triarea(i,j,ng) * mtsum(i,j,nt)
3522 : enddo ! ij
3523 :
3524 :
3525 2753280 : elseif (tracer_type(nt)==3) then ! depends on two tracers
3526 2753280 : nt1 = depend(nt)
3527 :
3528 399728710 : do ij = 1, icells(ng)
3529 396975430 : i = indxi(ij,ng)
3530 396975430 : j = indxj(ij,ng)
3531 :
3532 396975430 : i2 = iflux(i,j,ng)
3533 396975430 : j2 = jflux(i,j,ng)
3534 :
3535 : ! upwind approx (tx=ty=0) for type 3 tracers
3536 396975430 : mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt)
3537 :
3538 28845560 : mtflx(i,j,nt) = mtflx(i,j,nt) &
3539 399728710 : + triarea(i,j,ng) * mtsum(i,j,nt)
3540 : enddo ! ij
3541 :
3542 : endif ! tracer type
3543 : enddo ! ntrace
3544 : endif ! present(mtflx)
3545 : enddo ! ng
3546 :
3547 963648 : end subroutine transport_integrals
3548 :
3549 : !=======================================================================
3550 : !
3551 : ! Given transports through cell edges, compute new area and tracers.
3552 : !
3553 : ! author William H. Lipscomb, LANL
3554 :
3555 137664 : subroutine update_fields (nx_block, ny_block, &
3556 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
3557 : ntrace, & ! LCOV_EXCL_LINE
3558 : tracer_type, depend, & ! LCOV_EXCL_LINE
3559 : tarear, l_stop, & ! LCOV_EXCL_LINE
3560 : istop, jstop, & ! LCOV_EXCL_LINE
3561 : mflxe, mflxn, & ! LCOV_EXCL_LINE
3562 : mm, & ! LCOV_EXCL_LINE
3563 : mtflxe, mtflxn, & ! LCOV_EXCL_LINE
3564 114720 : tm)
3565 :
3566 : integer (kind=int_kind), intent(in) :: &
3567 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
3568 : ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ! LCOV_EXCL_LINE
3569 : ntrace ! number of tracers in use
3570 :
3571 : integer (kind=int_kind), dimension (ntrace), intent(in) :: &
3572 : tracer_type , & ! = 1, 2, or 3 (see comments above) ! LCOV_EXCL_LINE
3573 : depend ! tracer dependencies (see above)
3574 :
3575 : real (kind=dbl_kind), dimension (nx_block, ny_block), intent(in) :: &
3576 : mflxe, mflxn , & ! mass transport across east and north cell edges ! LCOV_EXCL_LINE
3577 : tarear ! 1/tarea
3578 :
3579 : real (kind=dbl_kind), dimension (nx_block, ny_block), intent(inout) :: &
3580 : mm ! mass field (mean)
3581 :
3582 : real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace), intent(in), optional :: &
3583 : mtflxe, mtflxn ! mass*tracer transport across E and N cell edges
3584 :
3585 : real (kind=dbl_kind), dimension (nx_block, ny_block, ntrace), intent(inout), optional :: &
3586 : tm ! tracer fields
3587 :
3588 : logical (kind=log_kind), intent(inout) :: &
3589 : l_stop ! if true, abort on return
3590 :
3591 : integer (kind=int_kind), intent(inout) :: &
3592 : istop, jstop ! indices of grid cell where model aborts
3593 :
3594 : ! local variables
3595 :
3596 : integer (kind=int_kind) :: &
3597 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
3598 : nt, nt1, nt2 ! tracer indices
3599 :
3600 : real (kind=dbl_kind), dimension(nx_block,ny_block,ntrace) :: &
3601 781089408 : mtold ! old mass*tracer
3602 :
3603 : real (kind=dbl_kind) :: &
3604 : puny, & ! ! LCOV_EXCL_LINE
3605 34560 : w1 ! work variable
3606 :
3607 : integer (kind=int_kind), dimension(nx_block*ny_block) :: &
3608 : indxi , & ! compressed indices in i and j directions ! LCOV_EXCL_LINE
3609 172224 : indxj
3610 :
3611 : integer (kind=int_kind) :: &
3612 : icells , & ! number of cells with mm > 0. ! LCOV_EXCL_LINE
3613 : ij ! combined i/j horizontal index
3614 :
3615 : character(len=*), parameter :: subname = '(update_fields)'
3616 :
3617 : !-------------------------------------------------------------------
3618 : ! Save starting values of mass*tracer
3619 : !-------------------------------------------------------------------
3620 :
3621 137664 : call icepack_query_parameters(puny_out=puny)
3622 137664 : call icepack_warnings_flush(nu_diag)
3623 137664 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
3624 0 : file=__FILE__, line=__LINE__)
3625 :
3626 137664 : if (present(tm)) then
3627 3097440 : do nt = 1, ntrace
3628 3097440 : if (tracer_type(nt)==1) then ! does not depend on other tracers
3629 21230640 : do j = jlo, jhi
3630 414807840 : do i = ilo, ihi
3631 414119520 : mtold(i,j,nt) = mm(i,j) * tm(i,j,nt)
3632 : enddo ! i
3633 : enddo ! j
3634 2294400 : elseif (tracer_type(nt)==2) then ! depends on another tracer
3635 2064960 : nt1 = depend(nt)
3636 63691920 : do j = jlo, jhi
3637 1244423520 : do i = ilo, ihi
3638 1242358560 : mtold(i,j,nt) = mm(i,j) * tm(i,j,nt1) * tm(i,j,nt)
3639 : enddo ! i
3640 : enddo ! j
3641 229440 : elseif (tracer_type(nt)==3) then ! depends on two tracers
3642 229440 : nt1 = depend(nt)
3643 229440 : nt2 = depend(nt1)
3644 7076880 : do j = jlo, jhi
3645 138269280 : do i = ilo, ihi
3646 41760000 : mtold(i,j,nt) = mm(i,j) &
3647 138039840 : * tm(i,j,nt2) * tm(i,j,nt1) * tm(i,j,nt)
3648 : enddo ! i
3649 : enddo ! j
3650 : endif ! depend(nt) = 0
3651 : enddo ! nt
3652 : endif ! present(tm)
3653 :
3654 : !-------------------------------------------------------------------
3655 : ! Update mass field
3656 : !-------------------------------------------------------------------
3657 :
3658 4246128 : do j = jlo, jhi
3659 82961568 : do i = ilo, ihi
3660 :
3661 50112000 : w1 = mflxe(i,j) - mflxe(i-1,j ) &
3662 103771440 : + mflxn(i,j) - mflxn(i ,j-1)
3663 78715440 : mm(i,j) = mm(i,j) - w1*tarear(i,j)
3664 :
3665 82823904 : if (mm(i,j) < -puny) then ! abort with negative value
3666 0 : l_stop = .true.
3667 0 : istop = i
3668 0 : jstop = j
3669 78715440 : elseif (mm(i,j) < c0) then ! set to zero
3670 86196 : mm(i,j) = c0
3671 : endif
3672 :
3673 : enddo
3674 : enddo
3675 :
3676 137664 : if (l_stop) then
3677 0 : i = istop
3678 0 : j = jstop
3679 0 : w1 = mflxe(i,j) - mflxe(i-1,j ) &
3680 0 : + mflxn(i,j) - mflxn(i ,j-1)
3681 0 : write (nu_diag,*) ' '
3682 0 : write (nu_diag,*) 'New mass < 0, i, j =', i, j
3683 0 : write (nu_diag,*) 'Old mass =', mm(i,j) + w1*tarear(i,j)
3684 0 : write (nu_diag,*) 'New mass =', mm(i,j)
3685 0 : write (nu_diag,*) 'Net transport =', -w1*tarear(i,j)
3686 0 : return
3687 : endif
3688 :
3689 : !-------------------------------------------------------------------
3690 : ! Update tracers
3691 : !-------------------------------------------------------------------
3692 :
3693 137664 : if (present(tm)) then
3694 :
3695 114720 : icells = 0
3696 3538440 : do j = jlo, jhi
3697 69134640 : do i = ilo, ihi
3698 69019920 : if (mm(i,j) > c0) then ! grid cells with positive areas
3699 33892021 : icells = icells + 1
3700 33892021 : indxi(icells) = i
3701 33892021 : indxj(icells) = j
3702 : endif
3703 : enddo ! i
3704 : enddo ! j
3705 :
3706 3097440 : do nt = 1, ntrace
3707 :
3708 91999440 : do j = jlo, jhi
3709 1797500640 : do i = ilo, ihi
3710 1794517920 : tm(i,j,nt) = c0
3711 : enddo
3712 : enddo
3713 :
3714 3097440 : if (tracer_type(nt)==1) then ! does not depend on other tracers
3715 :
3716 204040446 : do ij = 1, icells
3717 203352126 : i = indxi(ij)
3718 203352126 : j = indxj(ij)
3719 :
3720 37080468 : w1 = mtflxe(i,j,nt) - mtflxe(i-1,j ,nt) &
3721 221892360 : + mtflxn(i,j,nt) - mtflxn(i ,j-1,nt)
3722 18540234 : tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
3723 204040446 : / mm(i,j)
3724 : enddo ! ij
3725 :
3726 2294400 : elseif (tracer_type(nt)==2) then ! depends on another tracer
3727 2064960 : nt1 = depend(nt)
3728 :
3729 612121338 : do ij = 1, icells
3730 610056378 : i = indxi(ij)
3731 610056378 : j = indxj(ij)
3732 :
3733 612121338 : if (abs(tm(i,j,nt1)) > c0) then
3734 110928632 : w1 = mtflxe(i,j,nt) - mtflxe(i-1,j ,nt) &
3735 648459748 : + mtflxn(i,j,nt) - mtflxn(i ,j-1,nt)
3736 55464316 : tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
3737 592995432 : / (mm(i,j) * tm(i,j,nt1))
3738 : endif
3739 :
3740 : enddo ! ij
3741 :
3742 229440 : elseif (tracer_type(nt)==3) then ! depends on two tracers
3743 229440 : nt1 = depend(nt)
3744 229440 : nt2 = depend(nt1)
3745 :
3746 68013482 : do ij = 1, icells
3747 67784042 : i = indxi(ij)
3748 67784042 : j = indxj(ij)
3749 :
3750 73964120 : if (abs(tm(i,j,nt1)) > c0 .and. &
3751 6409518 : abs(tm(i,j,nt2)) > c0) then
3752 9344144 : w1 = mtflxe(i,j,nt) - mtflxe(i-1,j ,nt) &
3753 46485216 : + mtflxn(i,j,nt) - mtflxn(i ,j-1,nt)
3754 4672072 : tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
3755 41813144 : / (mm(i,j) * tm(i,j,nt2) * tm(i,j,nt1))
3756 : endif
3757 : enddo ! ij
3758 :
3759 : endif ! tracer_type
3760 : enddo ! nt
3761 : endif ! present(tm)
3762 :
3763 367104 : end subroutine update_fields
3764 :
3765 : !=======================================================================
3766 :
3767 : end module ice_transport_remap
3768 :
3769 : !=======================================================================
|
