ROMS
Loading...
Searching...
No Matches
tl_nesting.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE tl_nesting_mod
3
4#if defined NESTING && defined TANGENT
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! This module contains several tangent linear routines to process the !
14! connectivity between nested grids. It process the contact region !
15! points between data donor and data receiver grids. !
16! !
17! The locations of the linear interpolation weights in the donor !
18! grid with respect the receiver grid contact region at contact !
19! point x(Irg,Jrg,Krg) are: !
20! !
21! 8___________7 (Idg+1,Jdg+1,Kdg) !
22! /. /| !
23! / . / | !
24! (Idg,Jdg+1,Kdg) 5/___________/6 | !
25! | . | | !
26! | . x | | !
27! | 4.........|..|3 (Idg+1,Jdg+1,Kdg-1) !
28! | . | / !
29! |. | / !
30! |___________|/ !
31! (Idg,Jdg,Kdg-1) 1 2 !
32! !
33! Suffix: dg = donor grid !
34! rg = receiver grid !
35! !
36! Routines: !
37! ======== !
38! !
39! tl_nesting Public interface to time-stepping kernel !
40! !
41! tl_get_composite Composite grid, extracts contact points donor !
42! data !
43! tl_get_refine Refinement grid, extracts contact points donor !
44! data !
45! tl_put_composite Composite grid, fills contact points (spatial !
46! interpolation) !
47! tl_put_refine Refinement grid, fills contact points (spatial !
48! and temporal interpolation) !
49! !
50# ifdef NESTING_DEBUG
51! tl_check_massflux If refinement, checks mass fluxes between coarse !
52! and fine grids for volume conservation. It is !
53! use only for debugging and diagnostics. !
54# endif
55! tl_correct_tracer Corrects coarse grid tracer at the refinement !
56! grid boundary with the refined accumulated !
57! fluxes !
58! tl_fine2coarse Replaces coarse grid state variables with the !
59! averaged fine grid values (two-way nesting) !
60! !
61! tl_put_refine2d Interpolates (space-time) 2D state variables !
62! tl_put_refine3d Interpolates (space-time) 3D state variables !
63! !
64! tl_z_weights Sets donor grid vertical indices (cell holding !
65! contact point) and vertical interpolation !
66! weights !
67! !
68! WARNINGS: !
69! ======== !
70! !
71! All the routines contained in this module are inside of a parallel !
72! region, except the main driver routine "nesting", which is called !
73! serially several times from main2d or main3d to perform different !
74! tasks. Notice that the calls to private "get_***" and "put_***" !
75! routines need to be in separated parallel loops because of serial !
76! with partitions and shared-memory rules. Furthermore, the donor !
77! and receiver grids may have different tile partitions. There is no !
78! I/O management inside the nesting routines. !
79! !
80! The connectivity between donor and receiver grids can be complex. !
81! The horizontal mapping between grids is static and done outside of !
82! ROMS. Only the time-dependent vertical interpolation weights are !
83! computed here. The contact region points I- and J-cell indices !
84! between donor and receiver grids, and the horizontal interpolation !
85! weights are read from the input nesting connectivity NetCDF file. !
86! It makes the nesting efficient and greatly simplifies parallelism. !
87! !
88!=======================================================================
89!
90 implicit none
91!
92 PUBLIC :: tl_nesting
93# ifdef NESTING_DEBUG
94 PUBLIC :: tl_check_massflux
95# endif
96# ifdef SOLVE3D
97 PRIVATE :: tl_correct_tracer
98 PRIVATE :: tl_correct_tracer_tile
99# endif
100 PRIVATE :: tl_fine2coarse
101 PRIVATE :: tl_get_composite
102 PRIVATE :: tl_get_refine
103 PRIVATE :: tl_put_composite
104 PRIVATE :: tl_put_refine
105 PRIVATE :: tl_put_refine2d
106# ifdef SOLVE3D
107 PRIVATE :: tl_put_refine3d
108 PRIVATE :: tl_z_weights
109# endif
110!
111 CONTAINS
112!
113 SUBROUTINE tl_nesting (ng, model, isection)
114!
115!=======================================================================
116! !
117! This routine process the contact region points between composite !
118! grids. In composite grids, it is possible to have more than one !
119! contact region. !
120! !
121! On Input: !
122! !
123! ng Data receiver grid number (integer) !
124! model Calling model identifier (integer) !
125! isection Governing equations time-stepping section in !
126! main2d or main3d indicating which state !
127! variables to process (integer) !
128! !
129!=======================================================================
130!
131 USE mod_param
132 USE mod_parallel
133 USE mod_ncparam
134 USE mod_nesting
135 USE mod_scalars
136!
137# ifdef SOLVE3D
138!
139 USE set_depth_mod, ONLY : set_depth
140 USE tl_set_depth_mod, ONLY : tl_set_depth
141# endif
142 USE nesting_mod, ONLY : get_metrics
143 USE nesting_mod, ONLY : mask_hweights
144 USE strings_mod, ONLY : founderror
145!
146! Imported variable declarations.
147!
148 integer, intent(in) :: ng, model, isection
149!
150! Local variable declarations.
151!
152 logical :: lputfsur
153!
154 integer :: subs, tile, thread
155 integer :: ngc
156!
157 character (len=*), parameter :: myfile = &
158 & __FILE__
159
160# ifdef PROFILE
161!
162!-----------------------------------------------------------------------
163! Turn on time clocks.
164!-----------------------------------------------------------------------
165!
166 CALL wclock_on (ng, model, 36, __line__, myfile)
167# endif
168# ifdef SOLVE3D
169!
170!-----------------------------------------------------------------------
171! Process vertical indices and interpolation weigths associated with
172! depth. Currently, vertical weights are used in composite grids
173! because their grids are not coincident. They are not needed in
174! refinement grids because the donor and receiver grids have the same
175! number of vertical levels and matching bathymetry. However, in
176! the future, it is possible to have configurations that require
177! vertical weights in refinement. The switch "get_Vweights" controls
178! if such weights are computed or not. If false, it will accelerate
179! computations because of less distributed-memory communications.
180!-----------------------------------------------------------------------
181!
182 IF ((isection.eq.nzwgt).and.get_vweights) THEN
183 DO tile=last_tile(ng),first_tile(ng),-1
184!^ CALL z_weights (ng, model, tile)
185!^
186 CALL tl_z_weights (ng, model, tile)
187 END DO
188 RETURN
189 END IF
190# endif
191
192# if defined MASKING || defined WET_DRY
193!
194!-----------------------------------------------------------------------
195! If Land/Sea masking, scale horizontal interpolation weights to
196! account for land contact points. If wetting and drying, the scaling
197! is done at every time-step because masking is time dependent.
198!-----------------------------------------------------------------------
199!
200 IF (isection.eq.nmask) THEN
201 DO tile=last_tile(ng),first_tile(ng),-1
202 CALL mask_hweights (ng, model, tile)
203 END DO
204 RETURN
205 END IF
206# endif
207!
208!-----------------------------------------------------------------------
209! If refinement grid, process contact points.
210!-----------------------------------------------------------------------
211!
212 IF (refinedgrid(ng)) THEN
213!
214! Extract grid spacing metrics (on_u and om_v) and load then to
215! REFINE(:) structure. These metrics are needed to impose mass
216! flux at the finer grid physical boundaries. It need to be done
217! separately because parallelism partions between all nested grid.
218!
219 IF (isection.eq.ndxdy) THEN
220 DO tile=first_tile(ng),last_tile(ng),+1
221 CALL get_metrics (ng, model, tile)
222 END DO
223!
224! Extract and store donor grid data at contact points.
225!
226 ELSE IF (isection.eq.ngetd) THEN
227 DO tile=first_tile(ng),last_tile(ng),+1
228!^ CALL get_refine (ng, model, tile)
229!^
230 CALL tl_get_refine (ng, model, tile)
231 END DO
232!
233! Fill refinement grid contact points variables by interpolating
234! (space, time) from extracted donor grid data. The free-surface
235! needs to be processed first and in a separate parallel region
236! because of shared-memory applications.
237!
238 ELSE IF (isection.eq.nputd) THEN
239 lputfsur=.true.
240 DO tile=first_tile(ng),last_tile(ng),+1
241!^ CALL put_refine (ng, model, tile, LputFsur)
242!^
243 CALL tl_put_refine (ng, model, tile, lputfsur)
244 END DO
245!
246 lputfsur=.false.
247 DO tile=first_tile(ng),last_tile(ng),+1
248!^ CALL put_refine (ng, model, tile, LputFsur)
249!^
250 CALL tl_put_refine (ng, model, tile, lputfsur)
251 END DO
252
253
254# ifdef NESTING_DEBUG
255!
256! If refinement, check mass flux conservation between coarser and
257! finer grids.
258!
259 ELSE IF (isection.eq.nmflx) THEN
260 DO tile=first_tile(ng),last_tile(ng),+1
261!^ CALL check_massflux (ng, model, tile)
262!^
263 CALL tl_check_massflux (ng, model, tile)
264 END DO
265# endif
266
267# ifndef ONE_WAY
268!
269! Fine to coarse coupling: two-way nesting.
270!
271 ELSE IF (isection.eq.n2way) THEN
272
273 ngc=coarserdonor(ng) ! coarse grid number
274
275# if defined SOLVE3D && !defined NO_CORRECT_TRACER
276!
277! Correct coarse grid tracer values at the refinement grid, ng,
278! boundary with the refined accumulated fluxes (Hz*u*T/n, Hz*v*T/m).
279!
280 DO tile=first_tile(ngc),last_tile(ngc),+1
281!^ CALL correct_tracer (ngc, ng, model, tile)
282!^
283 CALL tl_correct_tracer (ngc, ng, model, tile)
284 END DO
285# endif
286!
287! Replace coarse grid 2D state variables with the averaged fine grid
288! values (two-way coupling).
289!
290 DO tile=last_tile(ngc),first_tile(ngc),-1
291!^ CALL fine2coarse (ng, model, r2dvar, tile)
292!^
293 CALL tl_fine2coarse (ng, model, r2dvar, tile)
294 END DO
295 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
296
297# ifdef SOLVE3D
298!
299! Update coarse grid depth variables. We have a new coarse grid
300! adjusted free-surface, Zt_avg1.
301!
302 DO tile=first_tile(ngc),last_tile(ngc),+1
303!^ CALL set_depth (ngc, tile, model)
304!^
305 CALL tl_set_depth (ngc, tile, model)
306 END DO
307!
308! Replace coarse grid 3D state variables with the averaged fine grid
309! values (two-way coupling).
310!
311 DO tile=last_tile(ngc),first_tile(ngc),-1
312!^ CALL fine2coarse (ng, model, r3dvar, tile)
313!^
314 CALL tl_fine2coarse (ng, model, r3dvar, tile)
315 END DO
316 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
317# endif
318# else
319!
320! Fine to coarse coupling (two-way nesting) is not activated!
321!
322 ELSE IF (isection.eq.n2way) THEN
323# endif
324 END IF
325!
326!-----------------------------------------------------------------------
327! Otherwise, process contact points in composite grid.
328!-----------------------------------------------------------------------
329!
330 ELSE
331!
332! Get composite grid contact points data from donor grid. It extracts
333! the donor grid cell data necessary to interpolate state variables
334! at each contact point.
335!
336 DO tile=first_tile(ng),last_tile(ng),+1
337!^ CALL get_composite (ng, model, isection, tile)
338!^
339 CALL tl_get_composite (ng, model, isection, tile)
340 END DO
341!
342! Fill composite grid contact points variables by interpolating from
343! extracted donor grid data.
344!
345 DO tile=last_tile(ng),first_tile(ng),-1
346!^ CALL put_composite (ng, model, isection, tile)
347!^
348 CALL tl_put_composite (ng, model, isection, tile)
349 END DO
350
351 END IF
352
353# ifdef PROFILE
354!
355!-----------------------------------------------------------------------
356! Turn off time clocks.
357!-----------------------------------------------------------------------
358!
359 CALL wclock_off (ng, model, 36, __line__, myfile)
360# endif
361!
362 RETURN
363 END SUBROUTINE tl_nesting
364!
365 SUBROUTINE tl_get_composite (ng, model, isection, tile)
366!
367!=======================================================================
368! !
369! This routine gets the donor grid data required to process the !
370! contact points of the current composite grid. It extracts the !
371! donor cell points containing each contact point. In composite !
372! grids, it is possible to have more than one contact region. !
373! !
374! The interpolation of composite grid contact points from donor !
375! grid data is carried out in a different parallel region using !
376! 'put_composite'. !
377! !
378! On Input: !
379! !
380! ng Composite grid number (integer) !
381! model Calling model identifier (integer) !
382! isection Governing equations time-stepping section in !
383! main2d or main3d indicating which state !
384! variables to process (integer) !
385! tile Domain tile partition (integer) !
386! !
387! On Output: (mod_nesting) !
388! !
389! COMPOSITE Updated contact points structure. !
390! !
391!=======================================================================
392!
393 USE mod_param
394 USE mod_coupling
395 USE mod_forces
396 USE mod_grid
397 USE mod_ncparam
398 USE mod_nesting
399 USE mod_ocean
400 USE mod_scalars
401 USE mod_stepping
402!
403 USE nesting_mod, ONLY : get_contact2d
404# ifdef SOLVE3D
405 USE nesting_mod, ONLY : get_contact3d
406# endif
407!
408! Imported variable declarations.
409!
410 integer, intent(in) :: ng, model, isection, tile
411!
412! Local variable declarations.
413!
414 integer :: cr, dg, rg, nrec, rec
415# ifdef SOLVE3D
416 integer :: itrc
417# endif
418 integer :: lbi, ubi, lbj, ubj
419 integer :: tindex
420!
421!-----------------------------------------------------------------------
422! Get donor grid data needed to process composite grid contact points.
423! Only process those variables associated with the governing equation
424! time-stepping section.
425!-----------------------------------------------------------------------
426!
427 DO cr=1,ncontact
428!
429! Get data donor and data receiver grid numbers.
430!
431 dg=rcontact(cr)%donor_grid
432 rg=rcontact(cr)%receiver_grid
433!
434! Process only contact region data for requested nested grid "ng".
435!
436 IF (rg.eq.ng) THEN
437!
438! Set donor grid lower and upper array indices.
439!
440 lbi=bounds(dg)%LBi(tile)
441 ubi=bounds(dg)%UBi(tile)
442 lbj=bounds(dg)%LBj(tile)
443 ubj=bounds(dg)%UBj(tile)
444!
445! Process bottom stress (bustr, bvstr).
446!
447 IF (isection.eq.nbstr) THEN
448!^ CALL get_contact2d (dg, model, tile, &
449!^ & u2dvar, Vname(1,idUbms), &
450!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
451!^ & LBi, UBi, LBj, UBj, &
452!^ & FORCES(dg) % bustr, &
453!^ & COMPOSITE(cr) % bustr)
454!^
455 CALL get_contact2d (dg, model, tile, &
456 & u2dvar, vname(1,idubms), &
457 & cr, ucontact(cr)%Npoints, ucontact, &
458 & lbi, ubi, lbj, ubj, &
459 & forces(dg) % tl_bustr, &
460 & composite(cr) % tl_bustr)
461!^ CALL get_contact2d (dg, model, tile, &
462!^ & v2dvar, Vname(1,idVbms), &
463!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
464!^ & LBi, UBi, LBj, UBj, &
465!^ & FORCES(dg) % bvstr, &
466!^ & COMPOSITE(cr) % bvstr)
467!^
468 CALL get_contact2d (dg, model, tile, &
469 & v2dvar, vname(1,idvbms), &
470 & cr, vcontact(cr)%Npoints, vcontact, &
471 & lbi, ubi, lbj, ubj, &
472 & forces(dg) % tl_bvstr, &
473 & composite(cr) % tl_bvstr)
474 END IF
475!
476! Process free-surface (zeta) at the appropriate time index.
477!
478 IF ((isection.eq.nfsic).or. &
479 & (isection.eq.nzeta).or. &
480 & (isection.eq.n2dps).or. &
481 & (isection.eq.n2dcs)) THEN
482 IF (isection.eq.nzeta) THEN
483 nrec=2 ! process time records 1 and 2
484 ELSE
485 nrec=1 ! process knew record
486 END IF
487 DO rec=1,nrec
488 IF (isection.eq.nzeta) THEN
489 tindex=rec
490 ELSE
491 tindex=knew(dg)
492 END IF
493!^ CALL get_contact2d (dg, model, tile, &
494!^ & r2dvar, Vname(1,idFsur), &
495!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
496!^ & LBi, UBi, LBj, UBj, &
497!^ & OCEAN(dg) % zeta(:,:,Tindex), &
498!^ & COMPOSITE(cr) % zeta(:,:,rec))
499!^
500 CALL get_contact2d (dg, model, tile, &
501 & r2dvar, vname(1,idfsur), &
502 & cr, rcontact(cr)%Npoints, rcontact, &
503 & lbi, ubi, lbj, ubj, &
504 & ocean(dg) % tl_zeta(:,:,tindex), &
505 & composite(cr) % tl_zeta(:,:,rec))
506 END DO
507 END IF
508!
509! Process free-surface equation rigth-hand-side (rzeta) term.
510!
511 IF (isection.eq.n2dps) THEN
512 tindex=1
513!^ CALL get_contact2d (dg, model, tile, &
514!^ & r2dvar, Vname(1,idRzet), &
515!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
516!^ & LBi, UBi, LBj, UBj, &
517!^ & OCEAN(dg) % rzeta(:,:,Tindex), &
518!^ & COMPOSITE(cr) % rzeta)
519!^
520 CALL get_contact2d (dg, model, tile, &
521 & r2dvar, vname(1,idrzet), &
522 & cr, rcontact(cr)%Npoints, rcontact, &
523 & lbi, ubi, lbj, ubj, &
524 & ocean(dg) % tl_rzeta(:,:,tindex), &
525 & composite(cr) % tl_rzeta)
526 END IF
527!
528! Process 2D momentum components (ubar,vbar) at the appropriate time
529! index.
530!
531 IF ((isection.eq.n2dic).or. &
532 & (isection.eq.n2dps).or. &
533 & (isection.eq.n2dcs).or. &
534 & (isection.eq.n3duv)) THEN
535 IF (isection.eq.n3duv) THEN
536 nrec=2 ! process time records 1 and 2
537 ELSE
538 nrec=1 ! process knew record
539 END IF
540 DO rec=1,nrec
541 IF (isection.eq.n3duv) THEN
542 tindex=rec
543 ELSE
544 tindex=knew(dg)
545 END IF
546!^ CALL get_contact2d (dg, model, tile, &
547!^ & u2dvar, Vname(1,idUbar), &
548!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
549!^ & LBi, UBi, LBj, UBj, &
550!^ & OCEAN(dg) % ubar(:,:,Tindex), &
551!^ & COMPOSITE(cr) % ubar(:,:,rec))
552!^
553 CALL get_contact2d (dg, model, tile, &
554 & u2dvar, vname(1,idubar), &
555 & cr, ucontact(cr)%Npoints, ucontact, &
556 & lbi, ubi, lbj, ubj, &
557 & ocean(dg) % tl_ubar(:,:,tindex), &
558 & composite(cr) % tl_ubar(:,:,rec))
559!^ CALL get_contact2d (dg, model, tile, &
560!^ & v2dvar, Vname(1,idVbar), &
561!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
562!^ & LBi, UBi, LBj, UBj, &
563!^ & OCEAN(dg) % vbar(:,:,Tindex), &
564!^ & COMPOSITE(cr) % vbar(:,:,rec))
565!^
566 CALL get_contact2d (dg, model, tile, &
567 & v2dvar, vname(1,idvbar), &
568 & cr, vcontact(cr)%Npoints, vcontact, &
569 & lbi, ubi, lbj, ubj, &
570 & ocean(dg) % tl_vbar(:,:,tindex), &
571 & composite(cr) % tl_vbar(:,:,rec))
572 END DO
573 END IF
574
575# ifdef SOLVE3D
576!
577! Process time averaged free-surface (Zt_avg1) and 2D momentum fluxes
578! (DU_avg1, DV_avg1).
579!
580 IF (isection.eq.n2dfx) THEN
581!^ CALL get_contact2d (dg, model, tile, &
582!^ & r2dvar, 'Zt_avg1', &
583!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
584!^ & LBi, UBi, LBj, UBj, &
585!^ & COUPLING(dg) % Zt_avg1, &
586!^ & COMPOSITE(cr) % Zt_avg1)
587!^
588 CALL get_contact2d (dg, model, tile, &
589 & r2dvar, 'Zt_avg1', &
590 & cr, rcontact(cr)%Npoints, rcontact, &
591 & lbi, ubi, lbj, ubj, &
592 & coupling(dg) % tl_Zt_avg1, &
593 & composite(cr) % tl_Zt_avg1)
594
595 CALL get_contact2d (dg, model, tile, &
596 & u2dvar, 'DU_avg1', &
597 & cr, ucontact(cr)%Npoints, ucontact, &
598 & lbi, ubi, lbj, ubj, &
599 & coupling(dg) % DU_avg1, &
600 & composite(cr) % DU_avg1)
601 CALL get_contact2d (dg, model, tile, &
602 & u2dvar, 'DU_avg1', &
603 & cr, ucontact(cr)%Npoints, ucontact, &
604 & lbi, ubi, lbj, ubj, &
605 & coupling(dg) % tl_DU_avg1, &
606 & composite(cr) % tl_DU_avg1)
607
608 CALL get_contact2d (dg, model, tile, &
609 & v2dvar, 'DV_avg1', &
610 & cr, vcontact(cr)%Npoints, vcontact, &
611 & lbi, ubi, lbj, ubj, &
612 & coupling(dg) % DV_avg1, &
613 & composite(cr) % DV_avg1)
614 CALL get_contact2d (dg, model, tile, &
615 & v2dvar, 'DV_avg1', &
616 & cr, vcontact(cr)%Npoints, vcontact, &
617 & lbi, ubi, lbj, ubj, &
618 & coupling(dg) % tl_DV_avg1, &
619 & composite(cr) % tl_DV_avg1)
620 END IF
621
622# if !defined TS_FIXED
623!
624! Process tracer variables (t) at the appropriate time index.
625!
626 IF ((isection.eq.ntvic).or. &
627 & (isection.eq.nrhst).or. &
628 & (isection.eq.n3dtv)) THEN
629 DO itrc=1,nt(ng)
630 IF (isection.eq.nrhst) THEN
631 tindex=3
632 ELSE
633 tindex=nnew(dg)
634 END IF
635!^ CALL get_contact3d (dg, model, tile, &
636!^ & r3dvar, Vname(1,idTvar(itrc)), &
637!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
638!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
639!^ & OCEAN(dg) % t(:,:,:,Tindex,itrc), &
640!^ & COMPOSITE(cr) % t(:,:,:,itrc))
641!^
642 CALL get_contact3d (dg, model, tile, &
643 & r3dvar, vname(1,idtvar(itrc)), &
644 & cr, rcontact(cr)%Npoints, rcontact, &
645 & lbi, ubi, lbj, ubj, 1, n(dg), &
646 & ocean(dg) % tl_t(:,:,:,tindex,itrc), &
647 & composite(cr) % tl_t(:,:,:,itrc))
648 END DO
649 END IF
650# endif
651!
652! Process 3D momentum (u, v) at the appropriate time-index.
653!
654 IF ((isection.eq.n3dic).or. &
655 & (isection.eq.n3duv)) THEN
656 tindex=nnew(dg)
657!^ CALL get_contact3d (dg, model, tile, &
658!^ & u3dvar, Vname(1,idUvel), &
659!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
660!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
661!^ & OCEAN(dg) % u(:,:,:,Tindex), &
662!^ & COMPOSITE(cr) % u)
663!^
664 CALL get_contact3d (dg, model, tile, &
665 & u3dvar, vname(1,iduvel), &
666 & cr, ucontact(cr)%Npoints, ucontact, &
667 & lbi, ubi, lbj, ubj, 1, n(dg), &
668 & ocean(dg) % tl_u(:,:,:,tindex), &
669 & composite(cr) % tl_u)
670!^ CALL get_contact3d (dg, model, tile, &
671!^ & v3dvar, Vname(1,idVvel), &
672!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
673!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
674!^ & OCEAN(dg) % v(:,:,:,Tindex), &
675!^ & COMPOSITE(cr) % v)
676!^
677 CALL get_contact3d (dg, model, tile, &
678 & v3dvar, vname(1,idvvel), &
679 & cr, vcontact(cr)%Npoints, vcontact, &
680 & lbi, ubi, lbj, ubj, 1, n(dg), &
681 & ocean(dg) % tl_v(:,:,:,tindex), &
682 & composite(cr) % tl_v)
683 END IF
684!
685! Process 3D momentum fluxes (Huon, Hvom).
686!
687 IF (isection.eq.n3duv) THEN
688!^ CALL get_contact3d (dg, model, tile, &
689!^ & u3dvar, 'Huon', &
690!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
691!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
692!^ & GRID(dg) % Huon, &
693!^ & COMPOSITE(cr) % Huon)
694!^
695 CALL get_contact3d (dg, model, tile, &
696 & u3dvar, 'Huon', &
697 & cr, ucontact(cr)%Npoints, ucontact, &
698 & lbi, ubi, lbj, ubj, 1, n(dg), &
699 & grid(dg) % tl_Huon, &
700 & composite(cr) % tl_Huon)
701
702!^ CALL get_contact3d (dg, model, tile, &
703!^ & v3dvar, 'Hvom', &
704!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
705!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
706!^ & GRID(dg) % Hvom, &
707!^ & COMPOSITE(cr) % Hvom)
708!^
709 CALL get_contact3d (dg, model, tile, &
710 & v3dvar, 'Hvom', &
711 & cr, vcontact(cr)%Npoints, vcontact, &
712 & lbi, ubi, lbj, ubj, 1, n(dg), &
713 & grid(dg) % tl_Hvom, &
714 & composite(cr) % tl_Hvom)
715 END IF
716# endif
717
718 END IF
719 END DO
720!
721 RETURN
722 END SUBROUTINE tl_get_composite
723!
724 SUBROUTINE tl_get_refine (ng, model, tile)
725!
726!=======================================================================
727! !
728! This routine gets the donor grid data required to process the !
729! contact points of the current refinement grid. It extracts !
730! the donor cell points containing each contact point. !
731! !
732! The extracted data is stored in two-time rolling records which !
733! are needed for the space and time interpolation in 'put_refine'. !
734! !
735! Except for initialization, this routine is called at the bottom !
736! of the donor grid time step so all the values are updated for the !
737! time(dg)+dt(dg). That is, in 2D applications it is called after !
738! "step2d" corrector step and in 3D applications it is called after !
739! "step3d_t". This is done to have the coarser grid snapshots at !
740! time(dg) and time(dg)+dt(dg) to bound the interpolation of the !
741! finer grid contact points. !
742! !
743! On Input: !
744! !
745! ng Refinement grid number (integer) !
746! model Calling model identifier (integer) !
747! tile Domain tile partition (integer) !
748! !
749! On Output: (mod_nesting) !
750! !
751! REFINED Updated contact points structure. !
752! !
753!=======================================================================
754!
755 USE mod_param
756 USE mod_parallel
757 USE mod_coupling
758 USE mod_ncparam
759 USE mod_nesting
760 USE mod_ocean
761 USE mod_scalars
762 USE mod_stepping
763!
764 USE nesting_mod, ONLY : get_contact2d
765# ifdef SOLVE3D
766 USE nesting_mod, ONLY : get_contact3d
767# endif
768 USE nesting_mod, ONLY : get_persisted2d
769!
770! Imported variable declarations.
771!
772 integer, intent(in) :: ng, model, tile
773!
774! Local variable declarations.
775!
776# ifdef NESTING_DEBUG
777 logical, save :: first = .true.
778# endif
779 integer :: tindex2d, cr, dg, ir, rg, tnew
780# ifdef SOLVE3D
781 integer :: tindex3d, itrc
782# endif
783 integer :: lbi, ubi, lbj, ubj
784!
785!-----------------------------------------------------------------------
786! Get donor grid data needed to process refinement grid contact points.
787! The extracted contact point data is stored in two time records to
788! facilitate the space-time interpolation elsewhere.
789!-----------------------------------------------------------------------
790!
791 DO cr=1,ncontact
792!
793! Get data donor and data receiver grid numbers.
794!
795 dg=rcontact(cr)%donor_grid
796 rg=rcontact(cr)%receiver_grid
797!
798! Process only contact region data for requested nested grid "ng".
799!
800 IF ((dg.eq.coarserdonor(rg)).and.(dg.eq.ng)) THEN
801!
802! Set donor grid lower and upper array indices.
803!
804 lbi=bounds(dg)%LBi(tile)
805 ubi=bounds(dg)%UBi(tile)
806 lbj=bounds(dg)%LBj(tile)
807 ubj=bounds(dg)%UBj(tile)
808!
809! Update rolling time indices. The contact data is stored in two time
810! levels. We need a special case for ROMS initialization in "main2d"
811! or "main3d" after the processing "ini_fields". Notice that a dt(dg)
812! is added because this routine is called after the end of the time
813! step.
814!
815 IF (rollingindex(cr).eq.0) THEN
816 tnew=1 ! ROMS
817 rollingindex(cr)=tnew ! initialization
818 rollingtime(tnew,cr)=time(dg) ! step
819 ELSE
820 tnew=3-rollingindex(cr)
821 rollingindex(cr)=tnew
822 rollingtime(tnew,cr)=time(dg)+dt(dg)
823 END IF
824!
825! Set donor grid time index to process. In 3D applications, the 2D
826! record index to use can be either 1 or 2 since both ubar(:,:,1:2)
827! and vbar(:,:,1:2) are set to its time-averaged values in "step3d_uv".
828! That is, we can use Tindex2d=kstp(dg) or Tindex2d=knew(dg). However,
829! in 2D applications we need to use Tindex2d=knew(dg).
830!
831 tindex2d=knew(dg)
832# ifdef SOLVE3D
833 tindex3d=nnew(dg)
834# endif
835
836# ifdef NESTING_DEBUG
837!
838! If debugging, write information into Fortran unit 101 to check the
839! logic of processing donor grid data.
840!
841 IF (domain(ng)%SouthWest_Test(tile)) THEN
842 IF (master) THEN
843 IF (first) THEN
844 first=.false.
845 WRITE (101,10)
846 END IF
847 WRITE (101,20) ng, cr, dg, rg, iic(dg), iic(rg), &
848 & 3-tnew, tnew, tindex2d, tindex3d, &
849 & int(time(rg)), &
850 & int(rollingtime(3-tnew,cr)), &
851 & int(time(ng)), &
852 & int(rollingtime(tnew,cr))
853 10 FORMAT (2x,'ng',2x,'cr',2x,'dg',2x,'rg',5x,'iic',5x,'iic',&
854 & 2x,'told',2x,'tnew',2x,'Tindex',2x,'Tindex', &
855 & 9x,'time',8x,'time',8x,'time',8x,'time',/, &
856 & 20x,'(dg)',4x,'(rg)',18x,'2D',6x,'3D',9x,'(rg)', &
857 & 8x,'told',8x,'(ng)',8x,'tnew',/)
858 20 FORMAT (4(1x,i3),2(1x,i7),2(2x,i4),2(4x,i4),1x,4(2x,i10))
859 FLUSH (101)
860 END IF
861 END IF
862# endif
863!
864! Extract free-surface.
865!
866# ifdef SOLVE3D
867!^ CALL get_contact2d (dg, model, tile, &
868!^ & r2dvar, 'Zt_avg1', &
869!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
870!^ & LBi, UBi, LBj, UBj, &
871!^ & COUPLING(dg) % Zt_avg1, &
872!^ & REFINED(cr) % zeta(:,:,tnew))
873!^
874 CALL get_contact2d (dg, model, tile, &
875 & r2dvar, 'Zt_avg1', &
876 & cr, rcontact(cr)%Npoints, rcontact, &
877 & lbi, ubi, lbj, ubj, &
878 & coupling(dg) % tl_Zt_avg1, &
879 & refined(cr) % tl_zeta(:,:,tnew))
880# else
881!^ CALL get_contact2d (dg, model, tile, &
882!^ & r2dvar, 'zeta', &
883!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
884!^ & LBi, UBi, LBj, UBj, &
885!^ & OCEAN(dg) % zeta(:,:,Tindex2d), &
886!^ & REFINED(cr) % zeta(:,:,tnew))
887!^
888 CALL get_contact2d (dg, model, tile, &
889 & r2dvar, 'zeta', &
890 & cr, rcontact(cr)%Npoints, rcontact, &
891 & lbi, ubi, lbj, ubj, &
892 & ocean(dg) % tl_zeta(:,:,tindex2d), &
893 & refined(cr) % tl_zeta(:,:,tnew))
894# endif
895!
896! Extract 2D momentum components (ubar, vbar).
897!
898!^ CALL get_contact2d (dg, model, tile, &
899!^ & u2dvar, Vname(1,idUbar), &
900!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
901!^ & LBi, UBi, LBj, UBj, &
902!^ & OCEAN(dg) % ubar(:,:,Tindex2d), &
903!^ & REFINED(cr) % ubar(:,:,tnew))
904!^
905 CALL get_contact2d (dg, model, tile, &
906 & u2dvar, vname(1,idubar), &
907 & cr, ucontact(cr)%Npoints, ucontact, &
908 & lbi, ubi, lbj, ubj, &
909 & ocean(dg) % tl_ubar(:,:,tindex2d), &
910 & refined(cr) % tl_ubar(:,:,tnew))
911
912!^ CALL get_contact2d (dg, model, tile, &
913!^ & v2dvar, Vname(1,idVbar), &
914!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
915!^ & LBi, UBi, LBj, UBj, &
916!^ & OCEAN(dg) % vbar(:,:,Tindex2d), &
917!^ & REFINED(cr) % vbar(:,:,tnew))
918!^
919 CALL get_contact2d (dg, model, tile, &
920 & v2dvar, vname(1,idvbar), &
921 & cr, vcontact(cr)%Npoints, vcontact, &
922 & lbi, ubi, lbj, ubj, &
923 & ocean(dg) % tl_vbar(:,:,tindex2d), &
924 & refined(cr) % tl_vbar(:,:,tnew))
925
926# ifdef SOLVE3D
927!
928! Extract time-averaged fluxes (DU_avg2, DV_avg2). We will use latter
929! only the values at the finer grid physical boundary to impose mass
930! flux conservation in routine "put_refine2d".
931!
932 CALL get_persisted2d (dg, rg, model, tile, &
933 & u2dvar, 'DU_avg2', &
934 & cr, ucontact(cr)%Npoints, ucontact, &
935 & lbi, ubi, lbj, ubj, &
936 & coupling(dg) % DU_avg2, &
937 & refined(cr) % DU_avg2(:,:,tnew))
938 CALL get_persisted2d (dg, rg, model, tile, &
939 & u2dvar, 'DU_avg2', &
940 & cr, ucontact(cr)%Npoints, ucontact, &
941 & lbi, ubi, lbj, ubj, &
942 & coupling(dg) % tl_DU_avg2, &
943 & refined(cr) % tl_DU_avg2(:,:,tnew))
944
945 CALL get_persisted2d (dg, rg, model, tile, &
946 & v2dvar, 'DV_avg2', &
947 & cr, vcontact(cr)%Npoints, vcontact, &
948 & lbi, ubi, lbj, ubj, &
949 & coupling(dg) % DV_avg2, &
950 & refined(cr) % DV_avg2(:,:,tnew))
951 CALL get_persisted2d (dg, rg, model, tile, &
952 & v2dvar, 'DV_avg2', &
953 & cr, vcontact(cr)%Npoints, vcontact, &
954 & lbi, ubi, lbj, ubj, &
955 & coupling(dg) % tl_DV_avg2, &
956 & refined(cr) % tl_DV_avg2(:,:,tnew))
957!
958! Tracer-type variables.
959!
960 DO itrc=1,nt(dg)
961!^ CALL get_contact3d (dg, model, tile, &
962!^ & r3dvar, Vname(1,idTvar(itrc)), &
963!^ & cr, Rcontact(cr)%Npoints, Rcontact, &
964!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
965!^ & OCEAN(dg) % t(:,:,:,Tindex3d,itrc), &
966!^ & REFINED(cr) % t(:,:,:,tnew,itrc))
967!^
968 CALL get_contact3d (dg, model, tile, &
969 & r3dvar, vname(1,idtvar(itrc)), &
970 & cr, rcontact(cr)%Npoints, rcontact, &
971 & lbi, ubi, lbj, ubj, 1, n(dg), &
972 & ocean(dg) % tl_t(:,:,:,tindex3d,itrc), &
973 & refined(cr) % tl_t(:,:,:,tnew,itrc))
974 END DO
975!
976! Extract 3D momentum components (u, v).
977!
978!^ CALL get_contact3d (dg, model, tile, &
979!^ & u3dvar, Vname(1,idUvel), &
980!^ & cr, Ucontact(cr)%Npoints, Ucontact, &
981!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
982!^ & OCEAN(dg) % u(:,:,:,Tindex3d), &
983!^ & REFINED(cr) % u(:,:,:,tnew))
984!^
985 CALL get_contact3d (dg, model, tile, &
986 & u3dvar, vname(1,iduvel), &
987 & cr, ucontact(cr)%Npoints, ucontact, &
988 & lbi, ubi, lbj, ubj, 1, n(dg), &
989 & ocean(dg) % tl_u(:,:,:,tindex3d), &
990 & refined(cr) % tl_u(:,:,:,tnew))
991
992!^ CALL get_contact3d (dg, model, tile, &
993!^ & v3dvar, Vname(1,idVvel), &
994!^ & cr, Vcontact(cr)%Npoints, Vcontact, &
995!^ & LBi, UBi, LBj, UBj, 1, N(dg), &
996!^ & OCEAN(dg) % v(:,:,:,Tindex3d), &
997!^ & REFINED(cr) % v(:,:,:,tnew))
998!^
999 CALL get_contact3d (dg, model, tile, &
1000 & v3dvar, vname(1,idvvel), &
1001 & cr, vcontact(cr)%Npoints, vcontact, &
1002 & lbi, ubi, lbj, ubj, 1, n(dg), &
1003 & ocean(dg) % tl_v(:,:,:,tindex3d), &
1004 & refined(cr) % tl_v(:,:,:,tnew))
1005# endif
1006 END IF
1007 END DO
1008!
1009 RETURN
1010 END SUBROUTINE tl_get_refine
1011!
1012 SUBROUTINE tl_put_composite (ng, model, isection, tile)
1013!
1014!=======================================================================
1015! !
1016! This routine interpolates composite grid contact points from donor !
1017! grid data extracted in routine 'get_composite'. !
1018! !
1019! On Input: !
1020! !
1021! ng Composite grid number (integer) !
1022! model Calling model identifier (integer) !
1023! isection Governing equations time-stepping section in !
1024! main2d or main3d indicating which state !
1025! variables to process (integer) !
1026! tile Domain tile partition (integer) !
1027! !
1028!=======================================================================
1029!
1030 USE mod_param
1031 USE mod_coupling
1032 USE mod_forces
1033 USE mod_grid
1034 USE mod_ncparam
1035 USE mod_nesting
1036 USE mod_ocean
1037 USE mod_scalars
1038 USE mod_stepping
1039!
1040 USE nesting_mod, ONLY : put_contact2d
1041
1042# ifdef DISTRIBUTE
1043!
1044 USE mp_exchange_mod, ONLY : mp_exchange2d
1045# ifdef SOLVE3D
1046 USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
1047# endif
1048# endif
1049!
1050! Imported variable declarations.
1051!
1052 integer, intent(in) :: ng, model, isection, tile
1053!
1054! Local variable declarations.
1055!
1056 integer :: dg, rg, cr, nrec, rec
1057# ifdef SOLVE3D
1058 integer :: itrc
1059# endif
1060 integer :: lbi, ubi, lbj, ubj
1061 integer :: tindex
1062!
1063!-----------------------------------------------------------------------
1064! Interpolate composite grid contact points from donor grid data.
1065! Only process those variables associated with the governing equation
1066! time-stepping section.
1067!-----------------------------------------------------------------------
1068!
1069 cr_loop : DO cr=1,ncontact
1070!
1071! Get data donor and data receiver grid numbers.
1072!
1073 dg=rcontact(cr)%donor_grid
1074 rg=rcontact(cr)%receiver_grid
1075!
1076! Process only contact region data for requested nested grid "ng".
1077!
1078 IF (rg.eq.ng) THEN
1079!
1080! Set receiver grid lower and upper array indices.
1081!
1082 lbi=bounds(rg)%LBi(tile)
1083 ubi=bounds(rg)%UBi(tile)
1084 lbj=bounds(rg)%LBj(tile)
1085 ubj=bounds(rg)%UBj(tile)
1086!
1087! Process bottom stress (bustr, bvstr).
1088!
1089 IF (isection.eq.nbstr) THEN
1090 CALL put_contact2d (rg, model, tile, &
1091 & u2dvar, vname(1,idubms), &
1092 & cr, ucontact(cr)%Npoints, ucontact, &
1093 & lbi, ubi, lbj, ubj, &
1094# ifdef MASKING
1095 & grid(rg) % umask, &
1096# endif
1097 & composite(cr) % tl_bustr, &
1098 & forces(rg) % tl_bustr)
1099 CALL put_contact2d (rg, model, tile, &
1100 & v2dvar, vname(1,idvbms), &
1101 & cr, vcontact(cr)%Npoints, vcontact, &
1102 & lbi, ubi, lbj, ubj, &
1103# ifdef MASKING
1104 & grid(rg) % vmask, &
1105# endif
1106 & composite(cr) % tl_bvstr, &
1107 & forces(rg) % tl_bvstr)
1108# ifdef DISTRIBUTE
1109 CALL mp_exchange2d (rg, tile, model, 2, &
1110 & lbi, ubi, lbj, ubj, &
1111 & nghostpoints, &
1112 & ewperiodic(rg), nsperiodic(rg), &
1113 & forces(rg) % tl_bustr, &
1114 & forces(rg) % tl_bvstr)
1115# endif
1116 END IF
1117!
1118! Process free-surface (zeta) at the appropriate time index.
1119!
1120 IF ((isection.eq.nfsic).or. &
1121 & (isection.eq.nzeta).or. &
1122 & (isection.eq.n2dps).or. &
1123 & (isection.eq.n2dcs)) THEN
1124 IF (isection.eq.nzeta) THEN
1125 nrec=2 ! process time records 1 and 2
1126 ELSE
1127 nrec=1 ! process knew record
1128 END IF
1129 DO rec=1,nrec
1130 IF (isection.eq.nzeta) THEN
1131 tindex=rec
1132 ELSE
1133 tindex=knew(rg)
1134 END IF
1135 CALL put_contact2d (rg, model, tile, &
1136 & r2dvar, vname(1,idfsur), &
1137 & cr, rcontact(cr)%Npoints, rcontact, &
1138 & lbi, ubi, lbj, ubj, &
1139# ifdef MASKING
1140 & grid(rg) % rmask, &
1141# endif
1142 & composite(cr) % tl_zeta(:,:,rec), &
1143 & ocean(rg) % tl_zeta(:,:,tindex))
1144# ifdef DISTRIBUTE
1145 CALL mp_exchange2d (rg, tile, model, 1, &
1146 & lbi, ubi, lbj, ubj, &
1147 & nghostpoints, &
1148 & ewperiodic(rg), nsperiodic(rg), &
1149 & ocean(rg) % tl_zeta(:,:,tindex))
1150# endif
1151 END DO
1152 END IF
1153!
1154! Process free-surface equation rigth-hand-side (rzeta) term.
1155!
1156 IF (isection.eq.n2dps) THEN
1157 tindex=1
1158 CALL put_contact2d (rg, model, tile, &
1159 & r2dvar, vname(1,idrzet), &
1160 & cr, rcontact(cr)%Npoints, rcontact, &
1161 & lbi, ubi, lbj, ubj, &
1162# ifdef MASKING
1163 & grid(rg) % rmask, &
1164# endif
1165 & composite(cr) % tl_rzeta, &
1166 & ocean(rg) % tl_rzeta(:,:,tindex))
1167# ifdef DISTRIBUTE
1168 CALL mp_exchange2d (rg, tile, model, 1, &
1169 & lbi, ubi, lbj, ubj, &
1170 & nghostpoints, &
1171 & ewperiodic(rg), nsperiodic(rg), &
1172 & ocean(rg) % tl_rzeta(:,:,tindex))
1173# endif
1174 END IF
1175!
1176! Process 2D momentum components (ubar,vbar) at the appropriate time
1177! index.
1178!
1179 IF ((isection.eq.n2dic).or. &
1180 & (isection.eq.n2dps).or. &
1181 & (isection.eq.n2dcs).or. &
1182 & (isection.eq.n3duv)) THEN
1183 IF (isection.eq.n3duv) THEN
1184 nrec=2 ! process time records 1 and 2
1185 ELSE
1186 nrec=1 ! process knew record
1187 END IF
1188 DO rec=1,nrec
1189 IF (isection.eq.n3duv) THEN
1190 tindex=rec
1191 ELSE
1192 tindex=knew(rg)
1193 END IF
1194 CALL put_contact2d (rg, model, tile, &
1195 & u2dvar, vname(1,idubar), &
1196 & cr, ucontact(cr)%Npoints, ucontact, &
1197 & lbi, ubi, lbj, ubj, &
1198# ifdef MASKING
1199 & grid(rg) % umask, &
1200# endif
1201 & composite(cr) % tl_ubar(:,:,rec), &
1202 & ocean(rg) % tl_ubar(:,:,tindex))
1203 CALL put_contact2d (rg, model, tile, &
1204 & v2dvar, vname(1,idvbar), &
1205 & cr, vcontact(cr)%Npoints, vcontact, &
1206 & lbi, ubi, lbj, ubj, &
1207# ifdef MASKING
1208 & grid(rg) % vmask, &
1209# endif
1210 & composite(cr) % tl_vbar(:,:,rec), &
1211 & ocean(rg) % tl_vbar(:,:,tindex))
1212# ifdef DISTRIBUTE
1213 CALL mp_exchange2d (rg, tile, model, 2, &
1214 & lbi, ubi, lbj, ubj, &
1215 & nghostpoints, &
1216 & ewperiodic(rg), nsperiodic(rg), &
1217 & ocean(rg) % tl_ubar(:,:,tindex), &
1218 & ocean(rg) % tl_vbar(:,:,tindex))
1219# endif
1220 END DO
1221 END IF
1222
1223# ifdef SOLVE3D
1224!
1225! Process time averaged free-surface (Zt_avg1) and 2D momentum fluxes
1226! (DU_avg1, DV_avg1).
1227!
1228 IF (isection.eq.n2dfx) THEN
1229 CALL put_contact2d (rg, model, tile, &
1230 & r2dvar, 'Zt_avg1', &
1231 & cr, rcontact(cr)%Npoints, rcontact, &
1232 & lbi, ubi, lbj, ubj, &
1233# ifdef MASKING
1234 & grid(rg) % rmask, &
1235# endif
1236 & composite(cr) % tl_Zt_avg1, &
1237 & coupling(rg) % tl_Zt_avg1)
1238 CALL put_contact2d (rg, model, tile, &
1239 & u2dvar, vname(1,idufx1), &
1240 & cr, ucontact(cr)%Npoints, ucontact, &
1241 & lbi, ubi, lbj, ubj, &
1242# ifdef MASKING
1243 & grid(rg) % umask, &
1244# endif
1245 & composite(cr) % DU_avg1, &
1246 & coupling(rg) % DU_avg1)
1247 CALL put_contact2d (rg, model, tile, &
1248 & u2dvar, vname(1,idufx1), &
1249 & cr, ucontact(cr)%Npoints, ucontact, &
1250 & lbi, ubi, lbj, ubj, &
1251# ifdef MASKING
1252 & grid(rg) % umask, &
1253# endif
1254 & composite(cr) % tl_DU_avg1, &
1255 & coupling(rg) % tl_DU_avg1)
1256 CALL put_contact2d (rg, model, tile, &
1257 & v2dvar, vname(1,idvfx1), &
1258 & cr, vcontact(cr)%Npoints, vcontact, &
1259 & lbi, ubi, lbj, ubj, &
1260# ifdef MASKING
1261 & grid(rg) % vmask, &
1262# endif
1263 & composite(cr) % DV_avg1, &
1264 & coupling(rg) % DV_avg1)
1265 CALL put_contact2d (rg, model, tile, &
1266 & v2dvar, vname(1,idvfx1), &
1267 & cr, vcontact(cr)%Npoints, vcontact, &
1268 & lbi, ubi, lbj, ubj, &
1269# ifdef MASKING
1270 & grid(rg) % vmask, &
1271# endif
1272 & composite(cr) % tl_DV_avg1, &
1273 & coupling(rg) % tl_DV_avg1)
1274# ifdef DISTRIBUTE
1275 CALL mp_exchange2d (rg, tile, model, 2, &
1276 & lbi, ubi, lbj, ubj, &
1277 & nghostpoints, &
1278 & ewperiodic(rg), nsperiodic(rg), &
1279 & coupling(rg) % DU_avg1, &
1280 & coupling(rg) % DV_avg1)
1281 CALL mp_exchange2d (rg, tile, model, 3, &
1282 & lbi, ubi, lbj, ubj, &
1283 & nghostpoints, &
1284 & ewperiodic(rg), nsperiodic(rg), &
1285 & coupling(rg) % tl_Zt_avg1, &
1286 & coupling(rg) % tl_DU_avg1, &
1287 & coupling(rg) % tl_DV_avg1)
1288# endif
1289 END IF
1290
1291# if !defined TS_FIXED
1292!
1293! Process tracer variables (t) at the appropriate time index.
1294!
1295 IF ((isection.eq.ntvic).or. &
1296 & (isection.eq.nrhst).or. &
1297 & (isection.eq.n3dtv)) THEN
1298 DO itrc=1,nt(ng)
1299 IF (isection.eq.nrhst) THEN
1300 tindex=3
1301 ELSE
1302 tindex=nnew(rg)
1303 END IF
1304 CALL tl_put_contact3d (rg, model, tile, &
1305 & r3dvar, vname(1,idtvar(itrc)), &
1306 & cr, rcontact(cr)%Npoints, rcontact,&
1307 & lbi, ubi, lbj, ubj, 1, n(rg), &
1308# ifdef MASKING
1309 & grid(rg) % rmask, &
1310# endif
1311 & composite(cr) % t(:,:,:,itrc), &
1312 & composite(cr) % tl_t(:,:,:,itrc), &
1313 & ocean(rg) % tl_t(:,:,:,tindex, &
1314 & itrc))
1315 END DO
1316# ifdef DISTRIBUTE
1317 CALL mp_exchange4d (rg, tile, model, 1, &
1318 & lbi, ubi, lbj, ubj, 1, n(rg), 1, nt(rg),&
1319 & nghostpoints, &
1320 & ewperiodic(rg), nsperiodic(rg), &
1321 & ocean(rg) % tl_t(:,:,:,tindex,:))
1322# endif
1323 END IF
1324# endif
1325!
1326! Process 3D momentum (u, v) at the appropriate time-index.
1327!
1328 IF ((isection.eq.n3dic).or. &
1329 & (isection.eq.n3duv)) THEN
1330 tindex=nnew(rg)
1331 CALL tl_put_contact3d (rg, model, tile, &
1332 & u3dvar, vname(1,iduvel), &
1333 & cr, ucontact(cr)%Npoints, ucontact, &
1334 & lbi, ubi, lbj, ubj, 1, n(rg), &
1335# ifdef MASKING
1336 & grid(rg) % umask, &
1337# endif
1338 & composite(cr) % u, &
1339 & composite(cr) % tl_u, &
1340 & ocean(rg) % tl_u(:,:,:,tindex))
1341 CALL tl_put_contact3d (rg, model, tile, &
1342 & v3dvar, vname(1,idvvel), &
1343 & cr, vcontact(cr)%Npoints, vcontact, &
1344 & lbi, ubi, lbj, ubj, 1, n(rg), &
1345# ifdef MASKING
1346 & grid(rg) % vmask, &
1347# endif
1348 & composite(cr) % v, &
1349 & composite(cr) % tl_v, &
1350 & ocean(rg) % tl_v(:,:,:,tindex))
1351# ifdef DISTRIBUTE
1352 CALL mp_exchange3d (rg, tile, model, 2, &
1353 & lbi, ubi, lbj, ubj, 1, n(rg), &
1354 & nghostpoints, &
1355 & ewperiodic(rg), nsperiodic(rg), &
1356 & ocean(rg) % tl_u(:,:,:,tindex), &
1357 & ocean(rg) % tl_v(:,:,:,tindex))
1358# endif
1359 END IF
1360!
1361! Process 3D momentum fluxes (Huon, Hvom).
1362!
1363 IF (isection.eq.n3duv) THEN
1364 CALL tl_put_contact3d (rg, model, tile, &
1365 & u3dvar, 'Huon', &
1366 & cr, ucontact(cr)%Npoints, ucontact, &
1367 & lbi, ubi, lbj, ubj, 1, n(rg), &
1368# ifdef MASKING
1369 & grid(rg) % umask, &
1370# endif
1371 & composite(cr) % Huon, &
1372 & composite(cr) % tl_Huon, &
1373 & grid(rg) % tl_Huon)
1374 CALL tl_put_contact3d (rg, model, tile, &
1375 & v3dvar, 'Hvom', &
1376 & cr, vcontact(cr)%Npoints, vcontact, &
1377 & lbi, ubi, lbj, ubj, 1, n(rg), &
1378# ifdef MASKING
1379 & grid(rg) % vmask, &
1380# endif
1381 & composite(cr) % Hvom, &
1382 & composite(cr) % tl_Hvom, &
1383 & grid(rg) % tl_Hvom)
1384# ifdef DISTRIBUTE
1385 CALL mp_exchange3d (rg, tile, model, 2, &
1386 & lbi, ubi, lbj, ubj, 1, n(rg), &
1387 & nghostpoints, &
1388 & ewperiodic(rg), nsperiodic(rg), &
1389 & grid(rg) % tl_Huon, &
1390 & grid(rg) % tl_Hvom)
1391# endif
1392 END IF
1393# endif
1394
1395 END IF
1396 END DO cr_loop
1397!
1398 RETURN
1399 END SUBROUTINE tl_put_composite
1400!
1401 SUBROUTINE tl_put_refine (ng, model, tile, LputFsur)
1402!
1403!=======================================================================
1404! !
1405! This routine interpolates refinement grid contact points from donor !
1406! grid data extracted in routine 'get_refine'. Notice that because of !
1407! shared-memory parallelism, the free-surface is processed first and !
1408! in a different parallel region.
1409! !
1410! On Input: !
1411! !
1412! ng Refinement grid number (integer) !
1413! model Calling model identifier (integer) !
1414! tile Domain tile partition (integer) !
1415! LputFsur Switch to process or not free-surface (logical) !
1416! !
1417!=======================================================================
1418!
1419 USE mod_param
1420 USE mod_coupling
1421 USE mod_forces
1422 USE mod_grid
1423 USE mod_ncparam
1424 USE mod_nesting
1425 USE mod_ocean
1426 USE mod_scalars
1427 USE mod_stepping
1428!
1429! Imported variable declarations.
1430!
1431 logical, intent(in) :: lputfsur
1432 integer, intent(in) :: ng, model, tile
1433!
1434! Local variable declarations.
1435!
1436 integer :: dg, rg, cr, nrec, rec
1437# ifdef SOLVE3D
1438 integer :: itrc
1439# endif
1440 integer :: lbi, ubi, lbj, ubj
1441 integer :: tindex
1442!
1443!-----------------------------------------------------------------------
1444! Interpolate refinement grid contact points from donor grid data
1445! (space-time interpolation)
1446!-----------------------------------------------------------------------
1447!
1448 DO cr=1,ncontact
1449!
1450! Get data donor and data receiver grid numbers.
1451!
1452 dg=rcontact(cr)%donor_grid
1453 rg=rcontact(cr)%receiver_grid
1454!
1455! Process only contact region data for requested nested grid "ng", if
1456! donor grid is coarser than receiver grid. That is, we are only
1457! processing external contact points areas.
1458!
1459 IF ((rg.eq.ng).and.(dxmax(dg).gt.dxmax(rg))) THEN
1460!
1461! Set receiver grid lower and upper array indices.
1462!
1463 lbi=bounds(rg)%LBi(tile)
1464 ubi=bounds(rg)%UBi(tile)
1465 lbj=bounds(rg)%LBj(tile)
1466 ubj=bounds(rg)%UBj(tile)
1467!
1468! Fill free-surface separatelly.
1469!
1470 IF (lputfsur) THEN
1471 CALL tl_put_refine2d (ng, dg, cr, model, tile, lputfsur, &
1472 & lbi, ubi, lbj, ubj)
1473 ELSE
1474!
1475! Fill other 2D state variables (like momentum) contact points.
1476!
1477 CALL tl_put_refine2d (ng, dg, cr, model, tile, lputfsur, &
1478 & lbi, ubi, lbj, ubj)
1479
1480# ifdef SOLVE3D
1481!
1482! Fill 3D state variables contact points.
1483!
1484 CALL tl_put_refine3d (ng, dg, cr, model, tile, &
1485 & lbi, ubi, lbj, ubj)
1486# endif
1487 END IF
1488 END IF
1489 END DO
1490!
1491 RETURN
1492 END SUBROUTINE tl_put_refine
1493
1494# ifdef SOLVE3D
1495!
1496 SUBROUTINE tl_correct_tracer (ng, ngf, model, tile)
1497!
1498!=======================================================================
1499! !
1500! This routine corrects the tracer values in the coarser grid at the !
1501! location of the finer grid physical domain perimeter by comparing !
1502! vertically accumulated horizontal tracer flux (Hz*u*T/n, Hz*v*T/m) !
1503! in two-way nesting refinement: !
1504! !
1505! coarse grid, t(:,jb,:,nstp,:) = t(:,jb,:,nstp,:) - FacJ (west, !
1506! east) !
1507! t(ib,:,:,nstp,:) = t(ib,:,:,nstp,:) - FacI (south, !
1508! north) !
1509! where !
1510! !
1511! FacJ = (TFF(jb,itrc) - TFC(jb,itrc)) * !
1512! pm(:,jb) * pn(:,jb) / D(:,jb) !
1513! !
1514! TFF(ib,itrc) = SUM[SUM[Tflux(ib,k,itrc)]] finer !
1515! grid !
1516! for k=1:N, 1:RefineScale flux !
1517! !
1518! TFC(ib,itrc) = SUM[Tflux(ib,k,itrc)] coarser !
1519! grid !
1520! for k=1:N flux !
1521! !
1522! Similarly, for the southern and northern tracer fluxes. !
1523! !
1524! !
1525! On Input: !
1526! !
1527! ngc Coarser grid number (integer) !
1528! ngf Finer grid number (integer) !
1529! model Calling model identifier (integer) !
1530! tile Domain tile partition (integer) !
1531! !
1532! On Output: (mod_ocean) !
1533! !
1534! t Updated coarse grid tracer values at finer grid !
1535! perimeter !
1536! !
1537!=======================================================================
1538!
1539 USE mod_param
1540!
1541! Imported variable declarations.
1542!
1543 integer, intent(in) :: ng, ngf, model, tile
1544!
1545! Local variable declarations.
1546!
1547# include "tile.h"
1548!
1549 CALL tl_correct_tracer_tile (ng, ngf, model, tile, &
1550 lbi, ubi, lbj, ubj, &
1551 & imins, imaxs, jmins, jmaxs)
1552!
1553 RETURN
1554 END SUBROUTINE tl_correct_tracer
1555!
1556!***********************************************************************
1557 SUBROUTINE tl_correct_tracer_tile (ngc, ngf, model, tile, &
1558 & LBi, UBi, LBj, UBj, &
1559 & IminS, ImaxS, JminS, JmaxS)
1560!***********************************************************************
1561!
1562 USE mod_param
1563 USE mod_clima
1564 USE mod_grid
1565 USE mod_ocean
1566 USE mod_nesting
1567 USE mod_scalars
1568 USE mod_stepping
1569
1570# ifdef DISTRIBUTE
1571!
1572 USE mp_exchange_mod, ONLY : mp_exchange4d
1573# endif
1574!
1575! Imported variable declarations.
1576!
1577 integer, intent(in) :: ngc, ngf, model, tile
1578 integer, intent(in) :: lbi, ubi, lbj, ubj
1579 integer, intent(in) :: imins, imaxs, jmins, jmaxs
1580!
1581! Local variable declarations.
1582!
1583 integer :: iedge, ibc, ibc_min, ibc_max, ibf, io
1584 integer :: jedge, jbc, jbc_min, jbc_max, jbf, jo
1585 integer :: istr, iend, jstr, jend
1586 integer :: istrm2, iendp2, jstrm2, jendp2
1587 integer :: tindex, i, ic, isum, itrc, j, jsum, k, half
1588 integer :: cr, dg, dgcr, rg, rgcr
1589
1590 real(r8) :: tfc, tff, tvalue, cff
1591 real(r8) :: tl_tfc, tl_tff, tl_tvalue, tl_cff
1592
1593 real(r8) :: dinv(imins:imaxs,jmins:jmaxs)
1594 real(r8) :: tl_dinv(imins:imaxs,jmins:jmaxs)
1595!
1596!-----------------------------------------------------------------------
1597! Correct coarser grid tracer values at finer grid perimeter.
1598!-----------------------------------------------------------------------
1599!
1600! Determine contact regions where coarse grid is the donor and coarse
1601! grid is the receiver..
1602!
1603 DO cr=1,ncontact
1604 dg=donor_grid(cr)
1605 rg=receiver_grid(cr)
1606 IF ((ngc.eq.dg).and.(ngf.eq.rg)) THEN
1607 dgcr=cr ! coarse is donor
1608 ELSE IF ((ngc.eq.rg).and.(ngf.eq.dg)) THEN
1609 rgcr=cr ! coarse is receiver
1610 END IF
1611 END DO
1612!
1613! Set tile starting and ending indices for coarser grid.
1614!
1615 istr =bounds(ngc)%Istr (tile)
1616 iend =bounds(ngc)%Iend (tile)
1617 jstr =bounds(ngc)%Jstr (tile)
1618 jend =bounds(ngc)%Jend (tile)
1619!
1620 istrm2=bounds(ngc)%Istrm2(tile)
1621 iendp2=bounds(ngc)%Iendp2(tile)
1622 jstrm2=bounds(ngc)%Jstrm2(tile)
1623 jendp2=bounds(ngc)%Jendp2(tile)
1624!
1625! Compute coarser grid inverse water colunm thickness.
1626!
1627 DO j=jstrm2,jendp2
1628 DO i=istrm2,iendp2
1629 cff=grid(ngc)%Hz(i,j,1)
1630 tl_cff=grid(ngc)%tl_Hz(i,j,1)
1631 DO k=2,n(rg)
1632 cff=cff+grid(ngc)%Hz(i,j,k)
1633 tl_cff=tl_cff+grid(ngc)%tl_Hz(i,j,k)
1634 END DO
1635 dinv(i,j)=1.0_r8/cff
1636 tl_dinv(i,j)=-tl_cff*dinv(i,j)/cff
1637 END DO
1638 END DO
1639!
1640! Set finer grid center (half) and offset indices (Io and Jo) for
1641! coarser grid (I,J) coordinates.
1642!
1643 half=(refinescale(ngf)-1)/2
1644 io=half+1
1645 jo=half+1
1646!
1647! Set coarse grid tracer index to correct. Since the exchange of data
1648! is done at the bottom of main3d, we need to use the newest time
1649! index, I think.
1650!
1651 tindex=nstp(ngc) ! HGA: Why this index is stable?
1652!! Tindex=nnew(ngc) ! Gets a lot of noise at boundary
1653!
1654!=======================================================================
1655! Compute vertically integrated horizontal advective tracer flux for
1656! coarser at the finer grid physical boundary. Then, correct coarser
1657! grid tracer values at that boundary.
1658!=======================================================================
1659!
1660! Initialize tracer counter index. The "tclm" array is only allocated
1661! to the NTCLM fields that need to be processed. This is done to
1662! reduce memory.
1663!
1664 ic=0
1665!
1666 t_loop : DO itrc=1,nt(ngc)
1667 ic=ic+1
1668!
1669!-----------------------------------------------------------------------
1670! Finer grid western boundary.
1671!-----------------------------------------------------------------------
1672!
1673 ibc=i_left(ngf)
1674 jbc_min=j_bottom(ngf)
1675 jbc_max=j_top(ngf)-1 ! interior points, no top
1676! left corner
1677 DO jbc=jstr,jend
1678 IF (((istr.le.ibc-1).and.(ibc-1.le.iend)).and. &
1679 & ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
1680!
1681! Sum vertically coarse grid horizontal advective tracer flux,
1682! Hz*u*T/n, from last time-step.
1683!
1684 tfc=0.0_r8
1685 tl_tfc=0.0_r8
1686 DO k=1,n(ngc)
1687 tfc=tfc+bry_contact(iwest,rgcr)%Tflux(jbc,k,itrc)
1688 tl_tfc=tl_tfc+bry_contact(iwest,rgcr)%tl_Tflux(jbc,k,itrc)
1689 END DO
1690!
1691! Sum vertically and horizontally finer grid advective tracer flux.
1692! This is a vertical and horizontal J-integral because "RefineScale"
1693! sub-divisions are done in the finer grid in each single coarse grid
1694! at the J-edge.
1695!
1696 tff=0.0_r8
1697 tl_tff=0.0_r8
1698 jedge=jo+(jbc-jbc_min)*refinescale(ngf)
1699 DO jsum=-half,half
1700 jbf=jedge+jsum
1701 DO k=1,n(ngf)
1702 tff=tff+bry_contact(iwest,dgcr)%Tflux(jbf,k,itrc)
1703 tl_tff=tl_tff+ &
1704 & bry_contact(iwest,dgcr)%tl_Tflux(jbf,k,itrc)
1705 END DO
1706 END DO
1707!
1708! Zeroth order correction to fine grid time integral (RIL, 2016).
1709!
1710 tff=tff*dt(ngc)/dt(ngf)
1711 tl_tff=tl_tff*dt(ngc)/dt(ngf)
1712!
1713! Correct coarse grid tracer at the finer grid western boundary.
1714!
1715 cff=grid(ngc)%pm(ibc-1,jbc)* &
1716 & grid(ngc)%pn(ibc-1,jbc)* &
1717 & dinv(ibc-1,jbc)
1718 tl_cff=grid(ngc)%pm(ibc-1,jbc)* &
1719 & grid(ngc)%pn(ibc-1,jbc)* &
1720 & tl_dinv(ibc-1,jbc)
1721 DO k=1,n(ngc)
1722 tvalue=max(0.0_r8, &
1723 & ocean(ngc)%t(ibc-1,jbc,k,tindex,itrc)- &
1724 & cff*(tff-tfc))
1725 tl_tvalue=(0.5_r8- &
1726 & sign(0.5_r8,-(ocean(ngc)%t(ibc-1,jbc,k,tindex,itrc)- &
1727 & cff*(tff-tfc))))* &
1728 & (ocean(ngc)%tl_t(ibc-1,jbc,k,tindex,itrc)- &
1729 & tl_cff*(tff-tfc)-cff*(tl_tff-tl_tfc))
1730 IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
1731 tvalue=tvalue+ &
1732 & dt(ngc)*clima(ngc)%Tnudgcof(ibc-1,jbc,k,itrc)* &
1733 & (clima(ngc)%tclm(ibc-1,jbc,k,ic)-tvalue)
1734 tl_tvalue=tl_tvalue- &
1735 & dt(ngc)*clima(ngc)%Tnudgcof(ibc-1,jbc,k,itrc)* &
1736 & tl_tvalue
1737 END IF
1738# ifdef MASKING
1739 tvalue=tvalue*grid(ngc)%rmask(ibc-1,jbc)
1740 tl_tvalue=tl_tvalue*grid(ngc)%rmask(ibc-1,jbc)
1741# endif
1742!^ OCEAN(ngc)%t(Ibc-1,Jbc,k,Tindex,itrc)=Tvalue
1743!^
1744 ocean(ngc)%tl_t(ibc-1,jbc,k,tindex,itrc)=tl_tvalue
1745 END DO
1746 END IF
1747 END DO
1748!
1749!-----------------------------------------------------------------------
1750! Finer grid eastern boundary.
1751!-----------------------------------------------------------------------
1752!
1753 ibc=i_right(ngf)
1754 jbc_min=j_bottom(ngf)
1755 jbc_max=j_top(ngf)-1 ! interior points, no top
1756! right corner
1757 DO jbc=jstr,jend
1758 IF (((istr.le.ibc).and.(ibc.le.iend)).and. &
1759 & ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
1760!
1761! Sum vertically coarse grid horizontal advective tracer flux,
1762! Hz*u*T/n, from last time-step.
1763!
1764 tfc=0.0_r8
1765 tl_tfc=0.0_r8
1766 DO k=1,n(ngc)
1767 tfc=tfc+bry_contact(ieast,rgcr)%Tflux(jbc,k,itrc)
1768 tl_tfc=tl_tfc+bry_contact(ieast,rgcr)%tl_Tflux(jbc,k,itrc)
1769 END DO
1770!
1771! Sum vertically and horizontally finer grid advective tracer flux.
1772! This is a vertical and horizontal J-integral because "RefineScale"
1773! sub-divisions are done in the finer grid in each single coarse grid
1774! at the J-edge.
1775!
1776 tff=0.0_r8
1777 tl_tff=0.0_r8
1778 jedge=jo+(jbc-jbc_min)*refinescale(ngf)
1779 DO jsum=-half,half
1780 jbf=jedge+jsum
1781 DO k=1,n(ngf)
1782 tff=tff+bry_contact(ieast,dgcr)%Tflux(jbf,k,itrc)
1783 tl_tff=tl_tff+ &
1784 & bry_contact(ieast,dgcr)%tl_Tflux(jbf,k,itrc)
1785 END DO
1786 END DO
1787!
1788! Zeroth order correction to fine grid time integral (RIL, 2016).
1789!
1790 tff=tff*dt(ngc)/dt(ngf)
1791 tl_tff=tl_tff*dt(ngc)/dt(ngf)
1792!
1793! Correct coarse grid tracer at the finer grid eastern boundary.
1794!
1795 cff=grid(ngc)%pm(ibc,jbc)* &
1796 & grid(ngc)%pn(ibc,jbc)* &
1797 & dinv(ibc,jbc)
1798 tl_cff=grid(ngc)%pm(ibc,jbc)* &
1799 & grid(ngc)%pn(ibc,jbc)* &
1800 & tl_dinv(ibc,jbc)
1801 DO k=1,n(ngc)
1802 tvalue=max(0.0_r8, &
1803 & ocean(ngc)%t(ibc,jbc,k,tindex,itrc)- &
1804 & cff*(tff-tfc))
1805 tl_tvalue=(0.5_r8- &
1806 & sign(0.5_r8,-(ocean(ngc)%t(ibc,jbc,k,tindex,itrc)- &
1807 & cff*(tff-tfc))))* &
1808 & (ocean(ngc)%tl_t(ibc,jbc,k,tindex,itrc)- &
1809 & tl_cff*(tff-tfc)-cff*(tl_tff-tl_tfc))
1810 IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
1811 tvalue=tvalue+ &
1812 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)* &
1813 & (clima(ngc)%tclm(ibc,jbc,k,ic)-tvalue)
1814 tl_tvalue=tl_tvalue- &
1815 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)* &
1816 & tl_tvalue
1817 END IF
1818# ifdef MASKING
1819 tvalue=tvalue*grid(ngc)%rmask(ibc,jbc)
1820 tl_tvalue=tl_tvalue*grid(ngc)%rmask(ibc,jbc)
1821# endif
1822!^ OCEAN(ngc)%t(Ibc,Jbc,k,Tindex,itrc)=Tvalue
1823!^
1824 ocean(ngc)%tl_t(ibc,jbc,k,tindex,itrc)=tl_tvalue
1825 END DO
1826 END IF
1827 END DO
1828!
1829!-----------------------------------------------------------------------
1830! Finer grid southern boundary.
1831!-----------------------------------------------------------------------
1832!
1833 jbc=j_bottom(ngf)
1834 ibc_min=i_left(ngf)
1835 ibc_max=i_right(ngf)-1 ! interior points, no bottom
1836! right corner
1837 DO ibc=istr,iend
1838 IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and. &
1839 & ((jstr.le.jbc-1).and.(jbc-1.le.jend))) THEN
1840!
1841! Sum vertically coarse grid horizontal advective tracer flux,
1842! Hz*v*T/m, from last time-step.
1843!
1844 tfc=0.0_r8
1845 tl_tfc=0.0_r8
1846 DO k=1,n(ngc)
1847 tfc=tfc+bry_contact(isouth,rgcr)%Tflux(ibc,k,itrc)
1848 tl_tfc=tl_tfc+bry_contact(isouth,rgcr)%tl_Tflux(ibc,k,itrc)
1849 END DO
1850!
1851! Sum vertically and horizontally finer grid advective tracer flux.
1852! This is a vertical and horizontal I-integral because "RefineScale"
1853! sub-divisions are done in the finer grid in each single coarse grid
1854! at the I-edge.
1855!
1856 tff=0.0_r8
1857 tl_tff=0.0_r8
1858 iedge=io+(ibc-ibc_min)*refinescale(ngf)
1859 DO isum=-half,half
1860 ibf=iedge+isum
1861 DO k=1,n(ngf)
1862 tff=tff+bry_contact(isouth,dgcr)%Tflux(ibf,k,itrc)
1863 tl_tff=tl_tff+ &
1864 & bry_contact(isouth,dgcr)%tl_Tflux(ibf,k,itrc)
1865 END DO
1866 END DO
1867!
1868! Zeroth order correction to fine grid time integral (RIL, 2016).
1869!
1870 tff=tff*dt(ngc)/dt(ngf)
1871 tl_tff=tl_tff*dt(ngc)/dt(ngf)
1872!
1873! Correct coarse grid tracer at the finer grid southern boundary.
1874!
1875 cff=grid(ngc)%pm(ibc,jbc-1)* &
1876 & grid(ngc)%pn(ibc,jbc-1)* &
1877 & dinv(ibc,jbc-1)
1878 tl_cff=grid(ngc)%pm(ibc,jbc-1)* &
1879 & grid(ngc)%pn(ibc,jbc-1)* &
1880 & tl_dinv(ibc,jbc-1)
1881 DO k=1,n(ngc)
1882 tvalue=max(0.0_r8, &
1883 & ocean(ngc)%t(ibc,jbc-1,k,tindex,itrc)- &
1884 & cff*(tff-tfc))
1885 tl_tvalue=(0.5_r8- &
1886 & sign(0.5_r8,-(ocean(ngc)%t(ibc,jbc-1,k,tindex,itrc)- &
1887 & cff*(tff-tfc))))* &
1888 & (ocean(ngc)%tl_t(ibc,jbc-1,k,tindex,itrc)- &
1889 & tl_cff*(tff-tfc)-cff*(tl_tff-tl_tfc))
1890 IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
1891 tvalue=tvalue+ &
1892 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc-1,k,itrc)* &
1893 & (clima(ngc)%tclm(ibc,jbc-1,k,ic)-tvalue)
1894 tl_tvalue=tl_tvalue- &
1895 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc-1,k,itrc)* &
1896 & tl_tvalue
1897 END IF
1898# ifdef MASKING
1899 tvalue=tvalue*grid(ngc)%rmask(ibc,jbc-1)
1900 tl_tvalue=tl_tvalue*grid(ngc)%rmask(ibc,jbc-1)
1901# endif
1902!^ OCEAN(ngc)%t(Ibc,Jbc-1,k,Tindex,itrc)=Tvalue
1903!^
1904 ocean(ngc)%tl_t(ibc,jbc-1,k,tindex,itrc)=tl_tvalue
1905 END DO
1906 END IF
1907 END DO
1908!
1909!-----------------------------------------------------------------------
1910! Finer grid northern boundary.
1911!-----------------------------------------------------------------------
1912!
1913 jbc=j_top(ngf)
1914 ibc_min=i_left(ngf)
1915 ibc_max=i_right(ngf)-1 ! interior points, no top
1916! right corner
1917 DO ibc=istr,iend
1918 IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and. &
1919 & ((jstr.le.jbc).and.(jbc.le.jend))) THEN
1920!
1921! Sum vertically coarse grid horizontal advective tracer flux,
1922! Hz*v*T/m, from last time-step.
1923!
1924 tfc=0.0_r8
1925 tl_tfc=0.0_r8
1926 DO k=1,n(ngc)
1927 tfc=tfc+bry_contact(inorth,rgcr)%Tflux(ibc,k,itrc)
1928 tl_tfc=tl_tfc+ &
1929 & bry_contact(inorth,rgcr)%tl_Tflux(ibc,k,itrc)
1930 END DO
1931!
1932! Sum vertically and horizontally finer grid advective tracer flux.
1933! This is a vertical and horizontal I-integral because "RefineScale"
1934! sub-divisions are done in the finer grid in each single coarse grid
1935! at the I-edge.
1936!
1937 tff=0.0_r8
1938 tl_tff=0.0_r8
1939 iedge=io+(ibc-ibc_min)*refinescale(ngf)
1940 DO isum=-half,half
1941 ibf=iedge+isum
1942 DO k=1,n(ngf)
1943 tff=tff+bry_contact(inorth,dgcr)%Tflux(ibf,k,itrc)
1944 tl_tff=tl_tff+ &
1945 & bry_contact(inorth,dgcr)%tl_Tflux(ibf,k,itrc)
1946 END DO
1947 END DO
1948!
1949! Zeroth order correction to fine grid time integral.
1950!
1951 tff=tff*dt(ngc)/dt(ngf)
1952 tl_tff=tl_tff*dt(ngc)/dt(ngf)
1953!
1954! Correct coarse grid tracer at the finer grid northern boundary.
1955!
1956 cff=grid(ngc)%pm(ibc,jbc)* &
1957 & grid(ngc)%pn(ibc,jbc)* &
1958 & dinv(ibc,jbc)
1959 tl_cff=grid(ngc)%pm(ibc,jbc)* &
1960 & grid(ngc)%pn(ibc,jbc)* &
1961 & tl_dinv(ibc,jbc)
1962 DO k=1,n(ngc)
1963 tvalue=max(0.0_r8, &
1964 & ocean(ngc)%t(ibc,jbc,k,tindex,itrc)- &
1965 & cff*(tff-tfc))
1966 tl_tvalue=(0.5_r8- &
1967 & sign(0.5_r8,-(ocean(ngc)%t(ibc,jbc,k,tindex,itrc)- &
1968 & cff*(tff-tfc))))* &
1969 & (ocean(ngc)%tl_t(ibc,jbc,k,tindex,itrc)- &
1970 & tl_cff*(tff-tfc)-cff*(tl_tff-tl_tfc))
1971 IF (ltracerclm(itrc,ngc).and.lnudgetclm(itrc,ngc)) THEN
1972 tvalue=tvalue+ &
1973 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)* &
1974 & (clima(ngc)%tclm(ibc,jbc,k,ic)-tvalue)
1975 tl_tvalue=tl_tvalue- &
1976 & dt(ngc)*clima(ngc)%Tnudgcof(ibc,jbc,k,itrc)* &
1977 & tl_tvalue
1978 END IF
1979# ifdef MASKING
1980 tvalue=tvalue*grid(ngc)%rmask(ibc,jbc)
1981 tl_tvalue=tl_tvalue*grid(ngc)%rmask(ibc,jbc)
1982# endif
1983!^ OCEAN(ngc)%t(Ibc,Jbc,k,Tindex,itrc)=Tvalue
1984!^
1985 ocean(ngc)%tl_t(ibc,jbc,k,tindex,itrc)=tl_tvalue
1986 END DO
1987 END IF
1988 END DO
1989 END DO t_loop
1990
1991# ifdef DISTRIBUTE
1992!
1993!-----------------------------------------------------------------------
1994! Exchange boundary data.
1995!-----------------------------------------------------------------------
1996!
1997!^ CALL mp_exchange4d (ngc, tile, model, 1, &
1998!^ & LBi, UBi, LBj, UBj, 1, N(ngc), &
1999!^ & 1, NT(ngc), &
2000!^ & NghostPoints, &
2001!^ & EWperiodic(ngc), NSperiodic(ngc), &
2002!^ & OCEAN(ngc)%t(:,:,:,Tindex,:))
2003!^
2004 CALL mp_exchange4d (ngc, tile, model, 1, &
2005 & lbi, ubi, lbj, ubj, 1, n(ngc), &
2006 & 1, nt(ngc), &
2007 & nghostpoints, &
2008 & ewperiodic(ngc), nsperiodic(ngc), &
2009 & ocean(ngc)%tl_t(:,:,:,tindex,:))
2010# endif
2011!
2012 RETURN
2013 END SUBROUTINE tl_correct_tracer_tile
2014# endif
2015!
2016 SUBROUTINE tl_fine2coarse (ng, model, vtype, tile)
2017!
2018!=======================================================================
2019! !
2020! This routine replaces interior coarse grid data with the refined !
2021! averaged values: two-way nesting. !
2022! !
2023! On Input: !
2024! !
2025! ng Refinement grid number (integer) !
2026! model Calling model identifier (integer) !
2027! vtype State variables to process (integer): !
2028! vtype = r2dvar 2D state variables !
2029! vtype = r3dvar 3D state variables !
2030! tile Domain tile partition (integer) !
2031! !
2032! On Output: (mod_coupling, mod_ocean) !
2033! !
2034! Updated state variable with average refined grid !
2035! solution !
2036! !
2037!=======================================================================
2038!
2039 USE mod_param
2040 USE mod_parallel
2041 USE mod_coupling
2042 USE mod_forces
2043 USE mod_grid
2044 USE mod_iounits
2045 USE mod_ncparam
2046 USE mod_nesting
2047 USE mod_ocean
2048 USE mod_scalars
2049 USE mod_stepping
2050
2051 USE nesting_mod, ONLY : fine2coarse2d
2052# ifdef SOLVE3D
2053 USE nesting_mod, ONLY : fine2coarse3d
2054# endif
2055!
2056 USE exchange_2d_mod
2057# ifdef SOLVE3D
2058 USE exchange_3d_mod
2059# endif
2060# ifdef DISTRIBUTE
2061 USE mp_exchange_mod, ONLY : mp_exchange2d
2062# ifdef SOLVE3D
2063 USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
2064# endif
2065# endif
2066 USE strings_mod, ONLY : founderror
2067!
2068! Imported variable declarations.
2069!
2070 integer, intent(in) :: ng, model, vtype, tile
2071!
2072! Local variable declarations.
2073!
2074 logical :: areaavg
2075!
2076 integer :: lbid, ubid, lbjd, ubjd
2077 integer :: lbir, ubir, lbjr, ubjr
2078 integer :: dindex2d, rindex2d
2079# ifdef SOLVE3D
2080 integer :: dindex3d, rindex3d
2081# endif
2082 integer :: cr, dg, k, rg, nrec, rec
2083# ifdef SOLVE3D
2084 integer :: itrc
2085# endif
2086!
2087 character (len=*), parameter :: myfile = &
2088 & __FILE__//", tl_fine2coarse"
2089!
2090!-----------------------------------------------------------------------
2091! Average interior fine grid state variable data to the coarse grid
2092! location. Then, replace coarse grid values with averaged data.
2093!-----------------------------------------------------------------------
2094!
2095 DO cr=1,ncontact
2096!
2097! Get data donor and data receiver grid numbers.
2098!
2099 dg=rcontact(cr)%donor_grid
2100 rg=rcontact(cr)%receiver_grid
2101!
2102! Process contact region if the current refinement grid "ng" is the
2103! donor grid. The coarse grid "rg" is the receiver grid and the
2104! contact structure has all the information necessary for fine to
2105! coarse coupling. The donor grid size is always smaller than the
2106! receiver coarser grid.
2107!
2108 IF ((ng.eq.dg).and.(dxmax(dg).lt.dxmax(rg))) THEN
2109!
2110! Set donor and receiver grids lower and upper array indices.
2111!
2112 lbid=bounds(dg)%LBi(tile)
2113 ubid=bounds(dg)%UBi(tile)
2114 lbjd=bounds(dg)%LBj(tile)
2115 ubjd=bounds(dg)%UBj(tile)
2116!
2117 lbir=bounds(rg)%LBi(tile)
2118 ubir=bounds(rg)%UBi(tile)
2119 lbjr=bounds(rg)%LBj(tile)
2120 ubjr=bounds(rg)%UBj(tile)
2121!
2122! Report.
2123!
2124 IF (domain(ng)%SouthWest_Test(tile)) THEN
2125 IF (master.and.(vtype.eq.r2dvar)) THEN
2126 WRITE (stdout,10) dg, rg, cr
2127 10 FORMAT (6x,'TL_FINE2COARSE - exchanging data between ', &
2128 & 'grids: dg = ',i2.2,' and rg = ',i2.2, &
2129 & ' at cr = ',i2.2)
2130 END IF
2131 END IF
2132!
2133! Set state variable indices to process for donor and receiver grids.
2134! Since the exchange of data is done at the bottom of main2d/main3d,
2135! we need to use the newest time indices.
2136!
2137 dindex2d=knew(dg) ! Donor 2D variables index
2138 rindex2d=knew(rg) ! Receiver 3D variables index
2139# ifdef SOLVE3D
2140 dindex3d=nnew(dg) ! Donor 3D variables index
2141 rindex3d=nnew(rg) ! Receiver 3D variables index
2142# endif
2143!
2144!-----------------------------------------------------------------------
2145! Process 2D state variables.
2146!-----------------------------------------------------------------------
2147!
2148 IF (vtype.eq.r2dvar) THEN
2149!
2150! Free-surface.
2151!
2152 areaavg=.false.
2153# ifdef SOLVE3D
2154 CALL fine2coarse2d (rg, dg, model, tile, &
2155 & r2dvar, 'Zt_avg1', &
2156 & areaavg, refinescale(dg), &
2157 & cr, rcontact(cr)%Npoints, rcontact, &
2158 & lbid, ubid, lbjd, ubjd, &
2159 & lbir, ubir, lbjr, ubjr, &
2160 & grid(dg)%om_r, &
2161 & grid(dg)%on_r, &
2162 & grid(rg)%pm, &
2163 & grid(rg)%pn, &
2164# ifdef MASKING
2165 & grid(dg)%rmask_full, &
2166 & grid(rg)%rmask, &
2167# endif
2168 & coupling(dg)%tl_Zt_avg1, &
2169 & coupling(rg)%tl_Zt_avg1)
2170# else
2171 CALL fine2coarse2d (rg, dg, model, tile, &
2172 & r2dvar, vname(1,idfsur), &
2173 & areaavg, refinescale(dg), &
2174 & cr, rcontact(cr)%Npoints, rcontact, &
2175 & lbid, ubid, lbjd, ubjd, &
2176 & lbir, ubir, lbjr, ubjr, &
2177 & grid(dg)%om_r, &
2178 & grid(dg)%on_r, &
2179 & grid(rg)%pm, &
2180 & grid(rg)%pn, &
2181# ifdef MASKING
2182 & grid(dg)%rmask, &
2183 & grid(rg)%rmask, &
2184# endif
2185 & ocean(dg)%tl_zeta(:,:,dindex2d), &
2186 & ocean(rg)%tl_zeta(:,:,rindex2d))
2187# endif
2188 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2189!
2190! Process 2D momentum components (ubar,vbar).
2191!
2192 areaavg=.false.
2193 CALL fine2coarse2d (rg, dg, model, tile, &
2194 & u2dvar, vname(1,idubar), &
2195 & areaavg, refinescale(dg), &
2196 & cr, ucontact(cr)%Npoints, ucontact, &
2197 & lbid, ubid, lbjd, ubjd, &
2198 & lbir, ubir, lbjr, ubjr, &
2199 & grid(dg)%om_u, &
2200 & grid(dg)%on_u, &
2201 & grid(rg)%pm, &
2202 & grid(rg)%pn, &
2203# ifdef MASKING
2204 & grid(dg)%umask_full, &
2205 & grid(rg)%umask_full, &
2206# endif
2207 & ocean(dg)%tl_ubar(:,:,dindex2d), &
2208# ifdef SOLVE3D
2209 & ocean(rg)%tl_ubar(:,:,1), &
2210 & ocean(rg)%tl_ubar(:,:,2))
2211# else
2212 & ocean(rg)%tl_ubar(:,:,rindex2d))
2213# endif
2214 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2215!
2216 CALL fine2coarse2d (rg, dg, model, tile, &
2217 & v2dvar, vname(1,idvbar), &
2218 & areaavg, refinescale(dg), &
2219 & cr, vcontact(cr)%Npoints, vcontact, &
2220 & lbid, ubid, lbjd, ubjd, &
2221 & lbir, ubir, lbjr, ubjr, &
2222 & grid(dg)%om_v, &
2223 & grid(dg)%on_v, &
2224 & grid(rg)%pm, &
2225 & grid(rg)%pn, &
2226# ifdef MASKING
2227 & grid(dg)%vmask_full, &
2228 & grid(rg)%vmask_full, &
2229# endif
2230 & ocean(dg)%tl_vbar(:,:,dindex2d), &
2231# ifdef SOLVE3D
2232 & ocean(rg)%tl_vbar(:,:,1), &
2233 & ocean(rg)%tl_vbar(:,:,2))
2234# else
2235 & ocean(rg)%tl_vbar(:,:,rindex2d))
2236# endif
2237 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2238
2239# ifdef SOLVE3D
2240!
2241!-----------------------------------------------------------------------
2242! Process 3D state variables.
2243!-----------------------------------------------------------------------
2244!
2245 ELSE IF (vtype.eq.r3dvar) THEN
2246!
2247! Tracer type-variables.
2248!
2249 areaavg=.false.
2250 DO itrc=1,nt(rg)
2251 CALL fine2coarse3d (rg, dg, model, tile, &
2252 & r3dvar, vname(1,idtvar(itrc)), &
2253 & areaavg, refinescale(dg), &
2254 & cr, rcontact(cr)%Npoints, rcontact, &
2255 & lbid, ubid, lbjd, ubjd, 1, n(dg), &
2256 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2257 & grid(dg)%om_r, &
2258 & grid(dg)%on_r, &
2259 & grid(rg)%pm, &
2260 & grid(rg)%pn, &
2261# ifdef MASKING
2262 & grid(dg)%rmask, &
2263 & grid(rg)%rmask, &
2264# endif
2265 & ocean(dg)%tl_t(:,:,:,dindex3d,itrc), &
2266 & ocean(rg)%tl_t(:,:,:,rindex3d,itrc))
2267 IF (founderror(exit_flag, noerror, &
2268 & __line__, myfile)) RETURN
2269 END DO
2270!
2271! Process 3D momentum components (u, v).
2272!
2273 areaavg=.false.
2274 CALL fine2coarse3d (rg, dg, model, tile, &
2275 & u3dvar, vname(1,iduvel), &
2276 & areaavg, refinescale(dg), &
2277 & cr, ucontact(cr)%Npoints, ucontact, &
2278 & lbid, ubid, lbjd, ubjd, 1, n(dg), &
2279 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2280 & grid(dg)%om_u, &
2281 & grid(dg)%on_u, &
2282 & grid(rg)%pm, &
2283 & grid(rg)%pn, &
2284# ifdef MASKING
2285 & grid(dg)%umask_full, &
2286 & grid(rg)%umask_full, &
2287# endif
2288 & ocean(dg)%tl_u(:,:,:,dindex3d), &
2289 & ocean(rg)%tl_u(:,:,:,rindex3d))
2290 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2291!
2292 CALL fine2coarse3d (rg, dg, model, tile, &
2293 & v3dvar, vname(1,idvvel), &
2294 & areaavg, refinescale(dg), &
2295 & cr, vcontact(cr)%Npoints, vcontact, &
2296 & lbid, ubid, lbjd, ubjd, 1, n(dg), &
2297 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2298 & grid(dg)%om_v, &
2299 & grid(dg)%on_v, &
2300 & grid(rg)%pm, &
2301 & grid(rg)%pn, &
2302# ifdef MASKING
2303 & grid(dg)%vmask_full, &
2304 & grid(rg)%vmask_full, &
2305# endif
2306 & ocean(dg)%tl_v(:,:,:,dindex3d), &
2307 & ocean(rg)%tl_v(:,:,:,rindex3d))
2308 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2309# endif
2310 END IF
2311!
2312!-----------------------------------------------------------------------
2313! Exchange boundary data.
2314!-----------------------------------------------------------------------
2315!
2316 IF (ewperiodic(rg).or.nsperiodic(rg)) THEN
2317 IF (vtype.eq.r2dvar) THEN
2318# ifdef SOLVE3D
2319!^ CALL exchange_r2d_tile (rg, tile, &
2320!^ & LBiR, UBiR, LBjR, UBjR, &
2321!^ & COUPLING(rg)%Zt_avg1)
2322!^
2323 CALL exchange_r2d_tile (rg, tile, &
2324 & lbir, ubir, lbjr, ubjr, &
2325 & coupling(rg)%tl_Zt_avg1)
2326 DO k=1,2
2327!^ CALL exchange_u2d_tile (rg, tile, &
2328!^ & LBiR, UBiR, LBjR, UBjR, &
2329!^ & OCEAN(rg)%ubar(:,:,k))
2330!^
2331 CALL exchange_u2d_tile (rg, tile, &
2332 & lbir, ubir, lbjr, ubjr, &
2333 & ocean(rg)%tl_ubar(:,:,k))
2334!^ CALL exchange_v2d_tile (rg, tile, &
2335!^ & LBiR, UBiR, LBjR, UBjR, &
2336!^ & OCEAN(rg)%vbar(:,:,k))
2337!^
2338 CALL exchange_v2d_tile (rg, tile, &
2339 & lbir, ubir, lbjr, ubjr, &
2340 & ocean(rg)%tl_vbar(:,:,k))
2341 END DO
2342# else
2343> CALL exchange_r2d_tile (rg, tile, &
2344!^ & LBiR, UBiR, LBjR, UBjR, &
2345!^ & OCEAN(rg)%zeta(:,:,Rindex2d))
2346!^
2347 CALL exchange_r2d_tile (rg, tile, &
2348 & lbir, ubir, lbjr, ubjr, &
2349 & ocean(rg)%tl_zeta(:,:,rindex2d))
2350!^ CALL exchange_u2d_tile (rg, tile, &
2351!^ & LBiR, UBiR, LBjR, UBjR, &
2352!^ & OCEAN(rg)%ubar(:,:,Rindex2d))
2353!^
2354 CALL exchange_u2d_tile (rg, tile, &
2355 & lbir, ubir, lbjr, ubjr, &
2356 & ocean(rg)%tl_ubar(:,:,rindex2d))
2357!^ CALL exchange_v2d_tile (rg, tile, &
2358!^ & LBiR, UBiR, LBjR, UBjR, &
2359!^ & OCEAN(rg)%vbar(:,:,Rindex2d))
2360!^
2361 CALL exchange_v2d_tile (rg, tile, &
2362 & lbir, ubir, lbjr, ubjr, &
2363 & ocean(rg)%tl_vbar(:,:,rindex2d))
2364
2365# endif
2366# ifdef SOLVE3D
2367 ELSE IF (vtype.eq.r3dvar) THEN
2368!^ CALL exchange_u3d_tile (rg, tile, &
2369!^ & LBiR, UBiR, LBjR, UBjR, 1, N(rg), &
2370!^ & OCEAN(rg)%u(:,:,:,Rindex3d))
2371!^
2372 CALL exchange_u3d_tile (rg, tile, &
2373 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2374 & ocean(rg)%tl_u(:,:,:,rindex3d))
2375!^ CALL exchange_v3d_tile (rg, tile, &
2376!^ & LBiR, UBiR, LBjR, UBjR, 1, N(rg), &
2377!^ & OCEAN(rg)%v(:,:,:,Rindex3d))
2378!^
2379 CALL exchange_v3d_tile (rg, tile, &
2380 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2381 & ocean(rg)%tl_v(:,:,:,rindex3d))
2382 DO itrc=1,nt(rg)
2383!^ CALL exchange_r3d_tile (rg, tile, &
2384!^ & LBiR, UBiR, LBjR, UBjR, &
2385!^ 1, N(rg), &
2386!^ & OCEAN(rg)%t(:,:,:,Rindex3d, &
2387!^ & itrc))
2388!^
2389 CALL exchange_r3d_tile (rg, tile, &
2390 & lbir, ubir, lbjr, ubjr, &
2391 & 1, n(rg), &
2392 & ocean(rg)%tl_t(:,:,:,rindex3d, &
2393 & itrc))
2394 END DO
2395# endif
2396 END IF
2397 END IF
2398
2399# ifdef DISTRIBUTE
2400!
2401 IF (vtype.eq.r2dvar) THEN
2402# ifdef SOLVE3D
2403!^ CALL mp_exchange2d (rg, tile, model, 1, &
2404!^ & LBiR, UBiR, LBjR, UBjR, &
2405!^ & NghostPoints, &
2406!^ & EWperiodic(rg), NSperiodic(rg), &
2407!^ & COUPLING(rg)%Zt_avg1)
2408!^
2409 CALL mp_exchange2d (rg, tile, model, 1, &
2410 & lbir, ubir, lbjr, ubjr, &
2411 & nghostpoints, &
2412 & ewperiodic(rg), nsperiodic(rg), &
2413 & coupling(rg)%tl_Zt_avg1)
2414!^ CALL mp_exchange2d (rg, tile, model, 4, &
2415!^ & LBiR, UBiR, LBjR, UBjR, &
2416!^ & NghostPoints, &
2417!^ & EWperiodic(rg), NSperiodic(rg), &
2418!^ & OCEAN(rg)%ubar(:,:,1), &
2419!^ & OCEAN(rg)%vbar(:,:,1), &
2420!^ & OCEAN(rg)%ubar(:,:,2), &
2421!^ & OCEAN(rg)%vbar(:,:,2))
2422!^
2423 CALL mp_exchange2d (rg, tile, model, 4, &
2424 & lbir, ubir, lbjr, ubjr, &
2425 & nghostpoints, &
2426 & ewperiodic(rg), nsperiodic(rg), &
2427 & ocean(rg)%tl_ubar(:,:,1), &
2428 & ocean(rg)%tl_vbar(:,:,1), &
2429 & ocean(rg)%tl_ubar(:,:,2), &
2430 & ocean(rg)%tl_vbar(:,:,2))
2431# else
2432!^ CALL mp_exchange2d (rg, tile, model, 3, &
2433!^ & LBiR, UBiR, LBjR, UBjR, &
2434!^ & NghostPoints, &
2435!^ & EWperiodic(rg), NSperiodic(rg), &
2436!^ & OCEAN(rg)%zeta(:,:,Rindex2d), &
2437!^ & OCEAN(rg)%ubar(:,:,Rindex2d), &
2438!^ & OCEAN(rg)%vbar(:,:,Rindex2d))
2439!^
2440 CALL mp_exchange2d (rg, tile, model, 3, &
2441 & lbir, ubir, lbjr, ubjr, &
2442 & nghostpoints, &
2443 & ewperiodic(rg), nsperiodic(rg), &
2444 & ocean(rg)%tl_zeta(:,:,rindex2d), &
2445 & ocean(rg)%tl_ubar(:,:,rindex2d), &
2446 & ocean(rg)%tl_vbar(:,:,rindex2d))
2447# endif
2448# ifdef SOLVE3D
2449 ELSE IF (vtype.eq.r3dvar) THEN
2450!^ CALL mp_exchange3d (rg, tile, model, 2, &
2451!^ & LBiR, UBiR, LBjR, UBjR, 1, N(rg), &
2452!^ & NghostPoints, &
2453!^ & EWperiodic(rg), NSperiodic(rg), &
2454!^ & OCEAN(rg)%u(:,:,:,Rindex3d), &
2455!^ & OCEAN(rg)%v(:,:,:,Rindex3d))
2456!^
2457 CALL mp_exchange3d (rg, tile, model, 2, &
2458 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2459 & nghostpoints, &
2460 & ewperiodic(rg), nsperiodic(rg), &
2461 & ocean(rg)%tl_u(:,:,:,rindex3d), &
2462 & ocean(rg)%tl_v(:,:,:,rindex3d))
2463!^ CALL mp_exchange4d (rg, tile, model, 1, &
2464!^ & LBiR, UBiR, LBjR, UBjR, 1, N(rg), &
2465!^ & 1, NT(rg), &
2466!^ & NghostPoints, &
2467!^ & EWperiodic(rg), NSperiodic(rg), &
2468!^ & OCEAN(rg)%t(:,:,:,Rindex3d,:))
2469!^
2470 CALL mp_exchange4d (rg, tile, model, 1, &
2471 & lbir, ubir, lbjr, ubjr, 1, n(rg), &
2472 & 1, nt(rg), &
2473 & nghostpoints, &
2474 & ewperiodic(rg), nsperiodic(rg), &
2475 & ocean(rg)%tl_t(:,:,:,rindex3d,:))
2476# endif
2477 END IF
2478# endif
2479 END IF
2480 END DO
2481!
2482 RETURN
2483 END SUBROUTINE tl_fine2coarse
2484!
2485 SUBROUTINE tl_put_refine2d (ng, dg, cr, model, tile, LputFsur, &
2486 & LBi, UBi, LBj, UBj)
2487!
2488!=======================================================================
2489! !
2490! This routine interpolates (space, time) refinement grid 2D state !
2491! variables contact points using data from the donor grid. Notice !
2492! that because of shared-memory parallelism, the free-surface is !
2493! processed first and in a different parallel region. !
2494! !
2495! On Input: !
2496! !
2497! ng Refinement (receiver) grid number (integer) !
2498! dg Donor grid number (integer) !
2499! cr Contact region number to process (integer) !
2500! model Calling model identifier (integer) !
2501! tile Domain tile partition (integer) !
2502! LputFsur Switch to process or not free-surface (logical) !
2503! LBi Receiver grid, I-dimension Lower bound (integer) !
2504! UBi Receiver grid, I-dimension Upper bound (integer) !
2505! LBj Receiver grid, J-dimension Lower bound (integer) !
2506! UBj Receiver grid, J-dimension Upper bound (integer) !
2507! !
2508! On Output: OCEAN(ng) structure !
2509! !
2510! zeta Updated free-surface !
2511! ubar Updated 2D momentum in the XI-direction !
2512! vbar Updated 2D momentum in the ETA-direction !
2513! !
2514!=======================================================================
2515!
2516 USE mod_param
2517 USE mod_parallel
2518 USE mod_coupling
2519 USE mod_grid
2520 USE mod_nesting
2521 USE mod_ocean
2522 USE mod_scalars
2523 USE mod_stepping
2524 USE mod_iounits
2525
2526# ifdef DISTRIBUTE
2527!
2528 USE distribute_mod, ONLY : mp_assemble
2529 USE mp_exchange_mod, ONLY : mp_exchange2d
2530# endif
2531 USE strings_mod, ONLY : founderror
2532!
2533! Imported variable declarations.
2534!
2535 logical, intent(in) :: lputfsur
2536 integer, intent(in) :: ng, dg, cr, model, tile
2537 integer, intent(in) :: lbi, ubi, lbj, ubj
2538!
2539! Local variable declarations.
2540!
2541 logical :: uboundary, vboundary
2542!
2543# ifdef DISTRIBUTE
2544 integer :: ilb, iub, jlb, jub, nptssn, nptswe, my_tile
2545# endif
2546 integer :: nsub, i, irec, j, m, tnew, told
2547 integer :: idg, jdg
2548!
2549# ifdef DISTRIBUTE
2550 real(r8), parameter :: spv = 0.0_r8
2551# endif
2552 real(dp) :: wnew, wold, secscale, fac
2553 real(r8) :: cff, cff1, tl_cff
2554 real(r8) :: my_value, tl_my_value
2555!
2556 character (len=*), parameter :: myfile = &
2557 & __FILE__//", tl_put_refined2d"
2558
2559# include "set_bounds.h"
2560!
2561!-----------------------------------------------------------------------
2562! Interpolate (space, time) refinement grid contact points for 2D state
2563! variables from donor grid.
2564!-----------------------------------------------------------------------
2565
2566# ifdef DISTRIBUTE
2567!
2568! Set global size of boundary edges.
2569!
2570 IF (.not.lputfsur) THEN
2571 my_tile=-1
2572 ilb=bounds(ng)%LBi(my_tile)
2573 iub=bounds(ng)%UBi(my_tile)
2574 jlb=bounds(ng)%LBj(my_tile)
2575 jub=bounds(ng)%UBj(my_tile)
2576 nptswe=jub-jlb+1
2577 nptssn=iub-ilb+1
2578
2579# ifdef NESTING_DEBUG
2580!
2581! If distributed-memory, initialize arrays used to check mass flux
2582! conservation with special value (zero) to facilitate the global
2583! reduction when collecting data between all nodes.
2584!
2585 bry_contact(iwest ,cr)%Mflux=spv
2586 bry_contact(ieast ,cr)%Mflux=spv
2587 bry_contact(isouth,cr)%Mflux=spv
2588 bry_contact(inorth,cr)%Mflux=spv
2589# endif
2590 END IF
2591# endif
2592!
2593! Set time snapshot indices for the donor grid data.
2594!
2595 told=3-rollingindex(cr)
2596 tnew=rollingindex(cr)
2597!
2598! Set linear time interpolation weights. Fractional seconds are
2599! rounded to the nearest milliseconds integer towards zero in the
2600! time interpolation weights.
2601!
2602 secscale=1000.0_dp ! seconds to milliseconds
2603!
2604 wold=anint((rollingtime(tnew,cr)-time(ng))*secscale,dp)
2605 wnew=anint((time(ng)-rollingtime(told,cr))*secscale,dp)
2606 fac=1.0_dp/(wold+wnew)
2607 wold=fac*wold
2608 wnew=fac*wnew
2609!
2610 IF (((wold*wnew).lt.0.0_dp).or.((wold+wnew).le.0.0_dp)) THEN
2611 IF (domain(ng)%SouthWest_Test(tile)) THEN
2612 IF (master) THEN
2613 WRITE (stdout,10) cr, dg, ng, &
2614 & iic(dg), told, tnew, &
2615 & iic(ng), wold, wnew, &
2616 & int(time(ng)), &
2617 & int(rollingtime(told,cr)), &
2618 & int(rollingtime(tnew,cr))
2619 END IF
2620 exit_flag=8
2621 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2622 END IF
2623 END IF
2624!
2625!-----------------------------------------------------------------------
2626! Process free-surface.
2627!-----------------------------------------------------------------------
2628!
2629 free_surface : IF (lputfsur) THEN
2630 DO m=1,rcontact(cr)%Npoints
2631 i=rcontact(cr)%Irg(m)
2632 j=rcontact(cr)%Jrg(m)
2633 IF (((istrt.le.i).and.(i.le.iendt)).and. &
2634 & ((jstrt.le.j).and.(j.le.jendt))) THEN
2635!^ my_value=Wold* &
2636!^ & (Rcontact(cr)%Lweight(1,m)* &
2637!^ & REFINED(cr)%zeta(1,m,told)+ &
2638!^ & Rcontact(cr)%Lweight(2,m)* &
2639!^ & REFINED(cr)%zeta(2,m,told)+ &
2640!^ & Rcontact(cr)%Lweight(3,m)* &
2641!^ & REFINED(cr)%zeta(3,m,told)+ &
2642!^ & Rcontact(cr)%Lweight(4,m)* &
2643!^ & REFINED(cr)%zeta(4,m,told))+ &
2644!^ & Wnew* &
2645!^ & (Rcontact(cr)%Lweight(1,m)* &
2646!^ & REFINED(cr)%zeta(1,m,tnew)+ &
2647!^ & Rcontact(cr)%Lweight(2,m)* &
2648!^ & REFINED(cr)%zeta(2,m,tnew)+ &
2649!^ & Rcontact(cr)%Lweight(3,m)* &
2650!^ & REFINED(cr)%zeta(3,m,tnew)+ &
2651!^ & Rcontact(cr)%Lweight(4,m)* &
2652!^ & REFINED(cr)%zeta(4,m,tnew))
2653!^
2654 tl_my_value=wold* &
2655 & (rcontact(cr)%Lweight(1,m)* &
2656 & refined(cr)%tl_zeta(1,m,told)+ &
2657 & rcontact(cr)%Lweight(2,m)* &
2658 & refined(cr)%tl_zeta(2,m,told)+ &
2659 & rcontact(cr)%Lweight(3,m)* &
2660 & refined(cr)%tl_zeta(3,m,told)+ &
2661 & rcontact(cr)%Lweight(4,m)* &
2662 & refined(cr)%tl_zeta(4,m,told))+ &
2663 & wnew* &
2664 & (rcontact(cr)%Lweight(1,m)* &
2665 & refined(cr)%tl_zeta(1,m,tnew)+ &
2666 & rcontact(cr)%Lweight(2,m)* &
2667 & refined(cr)%tl_zeta(2,m,tnew)+ &
2668 & rcontact(cr)%Lweight(3,m)* &
2669 & refined(cr)%tl_zeta(3,m,tnew)+ &
2670 & rcontact(cr)%Lweight(4,m)* &
2671 & refined(cr)%tl_zeta(4,m,tnew))
2672# ifdef MASKING
2673!^ my_value=my_value*GRID(ng)%rmask(i,j)
2674!^
2675 tl_my_value=tl_my_value*grid(ng)%rmask(i,j)
2676# endif
2677# ifdef WET_DRY
2678 IF (my_value.le.(dcrit(ng)-grid(ng)%h(i,j))) THEN
2679!^ my_value=Dcrit(ng)-GRID(ng)%h(i,j)
2680!^
2681 tl_my_value=-grid(ng)%tl_h(i,j)
2682 END IF
2683# endif
2684!^ OCEAN(ng)%zeta(i,j,1)=my_value
2685!^
2686 ocean(ng)%tl_zeta(i,j,1)=tl_my_value
2687!^ OCEAN(ng)%zeta(i,j,2)=my_value
2688!^
2689 ocean(ng)%tl_zeta(i,j,2)=tl_my_value
2690!^ OCEAN(ng)%zeta(i,j,3)=my_value
2691!^
2692 ocean(ng)%tl_zeta(i,j,3)=tl_my_value
2693# ifdef SOLVE3D
2694!^ COUPLING(ng)%Zt_avg1(i,j)=my_value
2695!^
2696 coupling(ng)%tl_Zt_avg1(i,j)=tl_my_value
2697# endif
2698 END IF
2699 END DO
2700
2701 ELSE
2702!
2703!-----------------------------------------------------------------------
2704! Process 2D momentum.
2705!-----------------------------------------------------------------------
2706!
2707! 2D momentum in the XI-direction.
2708# ifdef SOLVE3D
2709!
2710! Notice that contact points at the domain western and eastern
2711! boundaries are avoided for indx1(ng) time record. They are be
2712! assigned in the mass flux computations below. This exception is
2713! done for adjoint correctness.
2714# endif
2715!
2716 DO m=1,ucontact(cr)%Npoints
2717 i=ucontact(cr)%Irg(m)
2718 j=ucontact(cr)%Jrg(m)
2719 IF (((istrp.le.i).and.(i.le.iendt)).and. &
2720 & ((jstrt.le.j).and.(j.le.jendt))) THEN
2721!^ my_value=Wold* &
2722!^ & (Ucontact(cr)%Lweight(1,m)* &
2723!^ & REFINED(cr)%ubar(1,m,told)+ &
2724!^ & Ucontact(cr)%Lweight(2,m)* &
2725!^ & REFINED(cr)%ubar(2,m,told)+ &
2726!^ & Ucontact(cr)%Lweight(3,m)* &
2727!^ & REFINED(cr)%ubar(3,m,told)+ &
2728!^ & Ucontact(cr)%Lweight(4,m)* &
2729!^ & REFINED(cr)%ubar(4,m,told))+ &
2730!^ & Wnew* &
2731!^ & (Ucontact(cr)%Lweight(1,m)* &
2732!^ & REFINED(cr)%ubar(1,m,tnew)+ &
2733!^ & Ucontact(cr)%Lweight(2,m)* &
2734!^ & REFINED(cr)%ubar(2,m,tnew)+ &
2735!^ & Ucontact(cr)%Lweight(3,m)* &
2736!^ & REFINED(cr)%ubar(3,m,tnew)+ &
2737!^ & Ucontact(cr)%Lweight(4,m)* &
2738!^ & REFINED(cr)%ubar(4,m,tnew))
2739!^
2740 tl_my_value=wold* &
2741 & (ucontact(cr)%Lweight(1,m)* &
2742 & refined(cr)%tl_ubar(1,m,told)+ &
2743 & ucontact(cr)%Lweight(2,m)* &
2744 & refined(cr)%tl_ubar(2,m,told)+ &
2745 & ucontact(cr)%Lweight(3,m)* &
2746 & refined(cr)%tl_ubar(3,m,told)+ &
2747 & ucontact(cr)%Lweight(4,m)* &
2748 & refined(cr)%tl_ubar(4,m,told))+ &
2749 & wnew* &
2750 & (ucontact(cr)%Lweight(1,m)* &
2751 & refined(cr)%tl_ubar(1,m,tnew)+ &
2752 & ucontact(cr)%Lweight(2,m)* &
2753 & refined(cr)%tl_ubar(2,m,tnew)+ &
2754 & ucontact(cr)%Lweight(3,m)* &
2755 & refined(cr)%tl_ubar(3,m,tnew)+ &
2756 & ucontact(cr)%Lweight(4,m)* &
2757 & refined(cr)%tl_ubar(4,m,tnew))
2758# ifdef MASKING
2759!^ my_value=my_value*GRID(ng)%umask(i,j)
2760!^
2761 tl_my_value=tl_my_value*grid(ng)%umask(i,j)
2762# endif
2763# ifdef WET_DRY
2764!^ my_value=my_value*GRID(ng)%umask_wet(i,j)
2765!^
2766 tl_my_value=tl_my_value*grid(ng)%umask_wet(i,j)
2767# endif
2768 DO irec=1,3
2769# ifdef SOLVE3D
2770 uboundary=(m.eq.bry_contact(iwest,cr)%C2Bindex(j)).or. &
2771 & (m.eq.bry_contact(ieast,cr)%C2Bindex(j))
2772 IF(.not.(uboundary.and.(irec.eq.indx1(ng)))) THEN
2773!^ OCEAN(ng)%ubar(i,j,irec)=my_value
2774!^
2775 ocean(ng)%tl_ubar(i,j,irec)=tl_my_value
2776!! ELSE ! for debugging
2777!! OCEAN(ng)%ubar(i,j,irec)=0.0_r8 ! purposes
2778 END IF
2779# else
2780!^ OCEAN(ng)%ubar(i,j,irec)=my_value
2781!^
2782 ocean(ng)%tl_ubar(i,j,irec)=tl_my_value
2783# endif
2784 END DO
2785 END IF
2786 END DO
2787!
2788! 2D momentum in the ETA-direction.
2789# ifdef SOLVE3D
2790!
2791! Notice that contact points at the domain southern and northern
2792! boundaries are avoided for indx1(ng) time record. They are be
2793! assigned in the mass flux computations below. This exception is
2794! done for adjoint correctness.
2795# endif
2796!
2797 DO m=1,vcontact(cr)%Npoints
2798 i=vcontact(cr)%Irg(m)
2799 j=vcontact(cr)%Jrg(m)
2800 IF (((istrt.le.i).and.(i.le.iendt)).and. &
2801 & ((jstrp.le.j).and.(j.le.jendt))) THEN
2802!^ my_value=Wold* &
2803!^ & (Vcontact(cr)%Lweight(1,m)* &
2804!^ & REFINED(cr)%vbar(1,m,told)+ &
2805!^ & Vcontact(cr)%Lweight(2,m)* &
2806!^ & REFINED(cr)%vbar(2,m,told)+ &
2807!^ & Vcontact(cr)%Lweight(3,m)* &
2808!^ & REFINED(cr)%vbar(3,m,told)+ &
2809!^ & Vcontact(cr)%Lweight(4,m)* &
2810!^ & REFINED(cr)%vbar(4,m,told))+ &
2811!^ & Wnew* &
2812!^ & (Vcontact(cr)%Lweight(1,m)* &
2813!^ & REFINED(cr)%vbar(1,m,tnew)+ &
2814!^ & Vcontact(cr)%Lweight(2,m)* &
2815!^ & REFINED(cr)%vbar(2,m,tnew)+ &
2816!^ & Vcontact(cr)%Lweight(3,m)* &
2817!^ & REFINED(cr)%vbar(3,m,tnew)+ &
2818!^ & Vcontact(cr)%Lweight(4,m)* &
2819!^ & REFINED(cr)%vbar(4,m,tnew))
2820!^
2821 tl_my_value=wold* &
2822 & (vcontact(cr)%Lweight(1,m)* &
2823 & refined(cr)%tl_vbar(1,m,told)+ &
2824 & vcontact(cr)%Lweight(2,m)* &
2825 & refined(cr)%tl_vbar(2,m,told)+ &
2826 & vcontact(cr)%Lweight(3,m)* &
2827 & refined(cr)%tl_vbar(3,m,told)+ &
2828 & vcontact(cr)%Lweight(4,m)* &
2829 & refined(cr)%tl_vbar(4,m,told))+ &
2830 & wnew* &
2831 & (vcontact(cr)%Lweight(1,m)* &
2832 & refined(cr)%tl_vbar(1,m,tnew)+ &
2833 & vcontact(cr)%Lweight(2,m)* &
2834 & refined(cr)%tl_vbar(2,m,tnew)+ &
2835 & vcontact(cr)%Lweight(3,m)* &
2836 & refined(cr)%tl_vbar(3,m,tnew)+ &
2837 & vcontact(cr)%Lweight(4,m)* &
2838 & refined(cr)%tl_vbar(4,m,tnew))
2839# ifdef MASKING
2840!^ my_value=my_value*GRID(ng)%vmask(i,j)
2841!^
2842 tl_my_value=tl_my_value*grid(ng)%vmask(i,j)
2843# endif
2844# ifdef WET_DRY
2845!^ my_value=my_value*GRID(ng)%vmask_wet(i,j)
2846 tl_my_value=tl_my_value*grid(ng)%vmask_wet(i,j)
2847# endif
2848 DO irec=1,3
2849# ifdef SOLVE3D
2850 vboundary=(m.eq.bry_contact(isouth,cr)%C2Bindex(i)).or. &
2851 & (m.eq.bry_contact(inorth,cr)%C2Bindex(i))
2852 IF(.not.(vboundary.and.(irec.eq.indx1(ng)))) THEN
2853!^ OCEAN(ng)%vbar(i,j,irec)=my_value
2854!^
2855 ocean(ng)%tl_vbar(i,j,irec)=tl_my_value
2856!! ELSE ! for debugging
2857!! OCEAN(ng)%vbar(i,j,irec)=0.0_r8 ! purposes
2858 END IF
2859# else
2860!^ OCEAN(ng)%vbar(i,j,irec)=my_value
2861!^
2862 ocean(ng)%tl_vbar(i,j,irec)=tl_my_value
2863# endif
2864 END DO
2865 END IF
2866 END DO
2867
2868# ifdef SOLVE3D
2869!
2870!-----------------------------------------------------------------------
2871! Impose mass flux at the finer grid physical boundaries. This is only
2872! done for indx1(ng) time record.
2873!
2874! Western/Eastern boundary:
2875!
2876! ubar(Ibry,:,indx1) = DU_avg2(Ibry,:) * pn(Ibry,:) / D(Ibry,:)
2877!
2878! Southern/Northern boundary:
2879!
2880! vbar(:,Jbry,indx1) = DV_avg2(:,Jbry) * pm(:,Jbry) / D(:,Jbry)
2881!
2882! We use the latest coarse grid mass flux REFINED(cr)%DU_avg(1,:,tnew)
2883! with a linear variation (cff1) to ensure that the sum of the refined
2884! grid fluxes equals the coarse grid flux.
2885!-----------------------------------------------------------------------
2886!
2887! Western edge.
2888!
2889 IF (domain(ng)%Western_Edge(tile)) THEN
2890 DO j=jstr,jend
2891 m=bry_contact(iwest,cr)%C2Bindex(j)
2892 idg=ucontact(cr)%Idg(m) ! for debugging
2893 jdg=ucontact(cr)%Jdg(m) ! purposes
2894 cff=0.5_r8*grid(ng)%on_u(istr,j)* &
2895 (grid(ng)%h(istr-1,j)+ &
2896 & ocean(ng)%zeta(istr-1,j,indx1(ng))+ &
2897 & grid(ng)%h(istr ,j)+ &
2898 & ocean(ng)%zeta(istr ,j,indx1(ng)))
2899 tl_cff=0.5_r8*grid(ng)%on_u(istr,j)* &
2900 (grid(ng)%tl_h(istr-1,j)+ &
2901 & ocean(ng)%tl_zeta(istr-1,j,indx1(ng))+ &
2902 & grid(ng)%tl_h(istr ,j)+ &
2903 & ocean(ng)%tl_zeta(istr ,j,indx1(ng)))
2904 cff1=grid(ng)%on_u(istr,j)/refined(cr)%on_u(m)
2905# ifdef TIME_INTERP_FLUX
2906 my_value=cff1*(wold*refined(cr)%DU_avg2(1,m,told)+ &
2907 & wnew*refined(cr)%DU_avg2(1,m,tnew))/cff
2908 tl_my_value=cff1*(wold*refined(cr)%tl_DU_avg2(1,m,told)+ &
2909 & wnew*refined(cr)%tl_DU_avg2(1,m,tnew))/cff- &
2910 & tl_cff*my_value/cff
2911# else
2912 my_value=cff1*refined(cr)%DU_avg2(1,m,tnew)/cff
2913 tl_my_value=cff1*refined(cr)%tl_DU_avg2(1,m,tnew)/cff- &
2914 & tl_cff*my_value/cff
2915# endif
2916# ifdef MASKING
2917 my_value=my_value*grid(ng)%umask(istr,j)
2918 tl_my_value=tl_my_value*grid(ng)%umask(istr,j)
2919# endif
2920# ifdef WET_DRY
2921 my_value=my_value*grid(ng)%umask_wet(istr,j)
2922 tl_my_value=tl_my_value*grid(ng)%umask_wet(istr,j)
2923# endif
2924# ifdef NESTING_DEBUG
2925!^ BRY_CONTACT(iwest,cr)%Mflux(j)=cff*my_value
2926!^
2927 bry_contact(iwest,cr)%tl_Mflux(j)=cff*tl_my_value+ &
2928 & tl_cff*my_value
2929# endif
2930!^ OCEAN(ng)%ubar(Istr,j,indx1(ng))=my_value
2931!^
2932 ocean(ng)%tl_ubar(istr,j,indx1(ng))=tl_my_value
2933 END DO
2934 END IF
2935!
2936! Eastern edge.
2937!
2938 IF (domain(ng)%Eastern_Edge(tile)) THEN
2939 DO j=jstr,jend
2940 m=bry_contact(ieast,cr)%C2Bindex(j)
2941 idg=ucontact(cr)%Idg(m) ! for debugging
2942 jdg=ucontact(cr)%Jdg(m) ! purposes
2943 cff=0.5_r8*grid(ng)%on_u(iend+1,j)* &
2944 & (grid(ng)%h(iend+1,j)+ &
2945 & ocean(ng)%zeta(iend+1,j,indx1(ng))+ &
2946 & grid(ng)%h(iend ,j)+ &
2947 & ocean(ng)%zeta(iend ,j,indx1(ng)))
2948 tl_cff=0.5_r8*grid(ng)%on_u(iend+1,j)* &
2949 & (grid(ng)%tl_h(iend+1,j)+ &
2950 & ocean(ng)%tl_zeta(iend+1,j,indx1(ng))+ &
2951 & grid(ng)%tl_h(iend ,j)+ &
2952 & ocean(ng)%tl_zeta(iend ,j,indx1(ng)))
2953 cff1=grid(ng)%on_u(iend+1,j)/refined(cr)%on_u(m)
2954# ifdef TIME_INTERP_FLUX
2955 my_value=cff1*(wold*refined(cr)%DU_avg2(1,m,told)+ &
2956 & wnew*refined(cr)%DU_avg2(1,m,tnew))/cff
2957 tl_my_value=cff1*(wold*refined(cr)%tl_DU_avg2(1,m,told)+ &
2958 & wnew*refined(cr)%tl_DU_avg2(1,m,tnew))/cff- &
2959 & tl_cff*my_value/cff
2960# else
2961 my_value=cff1*refined(cr)%DU_avg2(1,m,tnew)/cff
2962 tl_my_value=cff1*refined(cr)%tl_DU_avg2(1,m,tnew)/cff- &
2963 & tl_cff*my_value/cff
2964# endif
2965# ifdef MASKING
2966 my_value=my_value*grid(ng)%umask(iend+1,j)
2967 tl_my_value=tl_my_value*grid(ng)%umask(iend+1,j)
2968# endif
2969# ifdef WET_DRY
2970 my_value=my_value*grid(ng)%umask_wet(iend+1,j)
2971 tl_my_value=tl_my_value*grid(ng)%umask_wet(iend+1,j)
2972# endif
2973# ifdef NESTING_DEBUG
2974!^ BRY_CONTACT(ieast,cr)%Mflux(j)=cff*my_value
2975!^
2976 bry_contact(ieast,cr)%tl_Mflux(j)=tl_cff*my_value+ &
2977 & cff*tl_my_value
2978# endif
2979!^ OCEAN(ng)%ubar(Iend+1,j,indx1(ng))=my_value
2980!^
2981 ocean(ng)%tl_ubar(iend+1,j,indx1(ng))=tl_my_value
2982 END DO
2983 END IF
2984!
2985! Southern edge.
2986!
2987 IF (domain(ng)%Southern_Edge(tile)) THEN
2988 DO i=istr,iend
2989 m=bry_contact(isouth,cr)%C2Bindex(i)
2990 idg=vcontact(cr)%Idg(m) ! for debugging
2991 jdg=vcontact(cr)%Jdg(m) ! purposes
2992 cff=0.5_r8*grid(ng)%om_v(i,jstr)* &
2993 & (grid(ng)%h(i,jstr-1)+ &
2994 & ocean(ng)%zeta(i,jstr-1,indx1(ng))+ &
2995 & grid(ng)%h(i,jstr )+ &
2996 & ocean(ng)%zeta(i,jstr ,indx1(ng)))
2997 tl_cff=0.5_r8*grid(ng)%om_v(i,jstr)* &
2998 & (grid(ng)%tl_h(i,jstr-1)+ &
2999 & ocean(ng)%tl_zeta(i,jstr-1,indx1(ng))+ &
3000 & grid(ng)%tl_h(i,jstr )+ &
3001 & ocean(ng)%tl_zeta(i,jstr ,indx1(ng)))
3002 cff1=grid(ng)%om_v(i,jstr)/refined(cr)%om_v(m)
3003# ifdef TIME_INTERP_FLUX
3004 my_value=cff1*(wold*refined(cr)%DV_avg2(1,m,told)+ &
3005 & wnew*refined(cr)%DV_avg2(1,m,tnew))/cff
3006 tl_my_value=cff1*(wold*refined(cr)%tl_DV_avg2(1,m,told)+ &
3007 & wnew*refined(cr)%tl_DV_avg2(1,m,tnew))/cff- &
3008 & tl_cff*my_value/cff
3009# else
3010 my_value=cff1*refined(cr)%DV_avg2(1,m,tnew)/cff
3011 tl_my_value=cff1*refined(cr)%tl_DV_avg2(1,m,tnew)/cff- &
3012 & tl_cff*my_value/cff
3013# endif
3014# ifdef MASKING
3015 my_value=my_value*grid(ng)%vmask(i,jstr)
3016 tl_my_value=tl_my_value*grid(ng)%vmask(i,jstr)
3017# endif
3018# ifdef WET_DRY
3019 my_value=my_value*grid(ng)%vmask_wet(i,jstr)
3020 tl_my_value=tl_my_value*grid(ng)%vmask_wet(i,jstr)
3021# endif
3022# ifdef NESTING_DEBUG
3023!^ BRY_CONTACT(isouth,cr)%Mflux(i)=cff*my_value
3024!^
3025 bry_contact(isouth,cr)%tl_Mflux(i)=tl_cff*my_value+ &
3026 & cff*tl_my_value
3027# endif
3028!^ OCEAN(ng)%vbar(i,Jstr,indx1(ng))=my_value
3029!^
3030 ocean(ng)%tl_vbar(i,jstr,indx1(ng))=tl_my_value
3031 END DO
3032 END IF
3033!
3034! Northern edge.
3035!
3036 IF (domain(ng)%Northern_Edge(tile)) THEN
3037 DO i=istr,iend
3038 m=bry_contact(inorth,cr)%C2Bindex(i)
3039 idg=vcontact(cr)%Idg(m) ! for debugging
3040 jdg=vcontact(cr)%Jdg(m) ! purposes
3041 cff=0.5_r8*grid(ng)%om_v(i,jend+1)* &
3042 & (grid(ng)%h(i,jend+1)+ &
3043 & ocean(ng)%zeta(i,jend+1,indx1(ng))+ &
3044 & grid(ng)%h(i,jend )+ &
3045 & ocean(ng)%zeta(i,jend ,indx1(ng)))
3046 tl_cff=0.5_r8*grid(ng)%om_v(i,jend+1)* &
3047 & (grid(ng)%tl_h(i,jend+1)+ &
3048 & ocean(ng)%tl_zeta(i,jend+1,indx1(ng))+ &
3049 & grid(ng)%tl_h(i,jend )+ &
3050 & ocean(ng)%tl_zeta(i,jend ,indx1(ng)))
3051 cff1=grid(ng)%om_v(i,jend+1)/refined(cr)%om_v(m)
3052# ifdef TIME_INTERP_FLUX
3053 my_value=cff1*(wold*refined(cr)%DV_avg2(1,m,told)+ &
3054 & wnew*refined(cr)%DV_avg2(1,m,tnew))/cff
3055 tl_my_value=cff1*(wold*refined(cr)%tl_DV_avg2(1,m,told)+ &
3056 & wnew*refined(cr)%tl_DV_avg2(1,m,tnew))/cff- &
3057 & tl_cff*my_value/cff
3058# else
3059 my_value=cff1*refined(cr)%DV_avg2(1,m,tnew)/cff
3060 tl_my_value=cff1*refined(cr)%tl_DV_avg2(1,m,tnew)/cff- &
3061 & tl_cff*my_value/cff
3062# endif
3063# ifdef MASKING
3064 my_value=my_value*grid(ng)%vmask(i,jend+1)
3065 tl_my_value=tl_my_value*grid(ng)%vmask(i,jend+1)
3066# endif
3067# ifdef WET_DRY
3068 my_value=my_value*grid(ng)%vmask_wet(i,jend+1)
3069 tl_my_value=tl_my_value*grid(ng)%vmask_wet(i,jend+1)
3070# endif
3071# ifdef NESTING_DEBUG
3072!^ BRY_CONTACT(inorth,cr)%Mflux(i)=cff*my_value
3073!^
3074 bry_contact(inorth,cr)%tl_Mflux(i)=tl_cff*my_value+ &
3075 & cff*tl_my_value
3076# endif
3077!^ OCEAN(ng)%vbar(i,Jend+1,indx1(ng))=my_value
3078!^
3079 ocean(ng)%tl_vbar(i,jend+1,indx1(ng))=tl_my_value
3080 END DO
3081 END IF
3082# endif
3083 END IF free_surface
3084
3085# ifdef DISTRIBUTE
3086!
3087!-----------------------------------------------------------------------
3088! Exchange tile information.
3089!-----------------------------------------------------------------------
3090!
3091! Free-surface.
3092!
3093 IF (lputfsur) THEN
3094!^ CALL mp_exchange2d (ng, tile, model, &
3095# ifdef SOLVE3D
3096!^ & 4, &
3097# else
3098!^ & 3, &
3099# endif
3100!^ & LBi, UBi, LBj, UBj, &
3101!^ & NghostPoints, &
3102!^ & EWperiodic(ng), NSperiodic(ng), &
3103# ifdef SOLVE3D
3104!^ & COUPLING(ng)%Zt_avg1, &
3105# endif
3106!^ & OCEAN(ng)%zeta(:,:,1), &
3107!^ & OCEAN(ng)%zeta(:,:,2), &
3108!^ & OCEAN(ng)%zeta(:,:,3))
3109!^
3110 CALL mp_exchange2d (ng, tile, model, &
3111# ifdef SOLVE3D
3112 & 4, &
3113# else
3114 & 3, &
3115# endif
3116 & lbi, ubi, lbj, ubj, &
3117 & nghostpoints, &
3118 & ewperiodic(ng), nsperiodic(ng), &
3119# ifdef SOLVE3D
3120 & coupling(ng)%tl_Zt_avg1, &
3121# endif
3122 & ocean(ng)%tl_zeta(:,:,1), &
3123 & ocean(ng)%tl_zeta(:,:,2), &
3124 & ocean(ng)%tl_zeta(:,:,3))
3125!
3126! 2D momentum.
3127!
3128 ELSE
3129!^ CALL mp_exchange2d (ng, tile, model, 3, &
3130!^ & LBi, UBi, LBj, UBj, &
3131!^ & NghostPoints, &
3132!^ & EWperiodic(ng), NSperiodic(ng), &
3133!^ & OCEAN(ng)%ubar(:,:,1), &
3134!^ & OCEAN(ng)%ubar(:,:,2), &
3135!^ & OCEAN(ng)%ubar(:,:,3))
3136!^
3137 CALL mp_exchange2d (ng, tile, model, 3, &
3138 & lbi, ubi, lbj, ubj, &
3139 & nghostpoints, &
3140 & ewperiodic(ng), nsperiodic(ng), &
3141 & ocean(ng)%tl_ubar(:,:,1), &
3142 & ocean(ng)%tl_ubar(:,:,2), &
3143 & ocean(ng)%tl_ubar(:,:,3))
3144
3145!^ CALL mp_exchange2d (ng, tile, model, 3, &
3146!^ & LBi, UBi, LBj, UBj, &
3147!^ & NghostPoints, &
3148!^ & EWperiodic(ng), NSperiodic(ng), &
3149!^ & OCEAN(ng)%vbar(:,:,1), &
3150!^ & OCEAN(ng)%vbar(:,:,2), &
3151!^ & OCEAN(ng)%vbar(:,:,3))
3152!^
3153 CALL mp_exchange2d (ng, tile, model, 3, &
3154 & lbi, ubi, lbj, ubj, &
3155 & nghostpoints, &
3156 & ewperiodic(ng), nsperiodic(ng), &
3157 & ocean(ng)%tl_vbar(:,:,1), &
3158 & ocean(ng)%tl_vbar(:,:,2), &
3159 & ocean(ng)%tl_vbar(:,:,3))
3160
3161# ifdef NESTING_DEBUG
3162!
3163!^ CALL mp_assemble (ng, model, NptsWE, spv, &
3164!^ & BRY_CONTACT(iwest ,cr)%Mflux(JLB:))
3165!^
3166 CALL mp_assemble (ng, model, nptswe, spv, &
3167 & bry_contact(iwest ,cr)%tl_Mflux(jlb:))
3168!^ CALL mp_assemble (ng, model, NptsWE, spv, &
3169!^ & BRY_CONTACT(ieast ,cr)%Mflux(JLB:))
3170!^
3171 CALL mp_assemble (ng, model, nptswe, spv, &
3172 & bry_contact(ieast ,cr)%tl_Mflux(jlb:))
3173!^ CALL mp_assemble (ng, model, NptsSN, spv, &
3174!^ & BRY_CONTACT(isouth,cr)%Mflux(ILB:))
3175!^
3176 CALL mp_assemble (ng, model, nptssn, spv, &
3177 & bry_contact(isouth,cr)%tl_Mflux(ilb:))
3178!^ CALL mp_assemble (ng, model, NptsSN, spv, &
3179!^ & BRY_CONTACT(inorth,cr)%Mflux(ILB:))
3180!^
3181 CALL mp_assemble (ng, model, nptssn, spv, &
3182 & bry_contact(inorth,cr)%tl_Mflux(ilb:))
3183# endif
3184 END IF
3185# endif
3186!
3187 10 FORMAT (/,' PUT_REFINE2D - unbounded contact points temporal: ', &
3188 & ' interpolation:', &
3189 & /,2x, 'cr = ',i2.2, &
3190 & 8x,'dg = ',i2.2, &
3191 & 8x,'ng = ',i2.2, &
3192 & /,2x, 'iic(dg) = ',i7.7, &
3193 & 3x,'told = ',i1, &
3194 & 9x,'tnew = ',i1, &
3195 & /,2x, 'iic(ng) = ',i7.7, &
3196 & 3x,'Wold = ',f8.5, &
3197 & 2x,'Wnew = ',f8.5, &
3198 & /,2x, 'time(ng) = ',i7.7, &
3199 & 3x,'time(told) = ',i7.7, &
3200 & 3x,'time(tnew) = ',i7.7)
3201!
3202 RETURN
3203 END SUBROUTINE tl_put_refine2d
3204
3205# ifdef SOLVE3D
3206!
3207 SUBROUTINE tl_put_refine3d (ng, dg, cr, model, tile, &
3208 & LBi, UBi, LBj, UBj)
3209!
3210!=======================================================================
3211! !
3212! This routine interpolates (space, time) refinement grid 3D state !
3213! variables contact points using data from the donor grid. !
3214! !
3215! On Input: !
3216! !
3217! ng Refinement (receiver) grid number (integer) !
3218! dg Donor grid number (integer) !
3219! cr Contact region number to process (integer) !
3220! model Calling model identifier (integer) !
3221! tile Domain tile partition (integer) !
3222! LBi Receiver grid, I-dimension Lower bound (integer) !
3223! UBi Receiver grid, I-dimension Upper bound (integer) !
3224! LBj Receiver grid, J-dimension Lower bound (integer) !
3225! UBj Receiver grid, J-dimension Upper bound (integer) !
3226! !
3227! On Output: OCEAN(ng) structure !
3228! !
3229! t Updated tracer-type variables !
3230! u Updated 3D momentum in the XI-direction !
3231! v Updated 3D momentum in the ETA-direction !
3232! !
3233!=======================================================================
3234!
3235 USE mod_param
3236 USE mod_parallel
3237 USE mod_grid
3238 USE mod_nesting
3239 USE mod_ocean
3240 USE mod_scalars
3241 USE mod_stepping
3242 USE mod_iounits
3243!
3244# ifdef DISTRIBUTE
3245 USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
3246# endif
3247 USE strings_mod, ONLY : founderror
3248!
3249! Imported variable declarations.
3250!
3251 integer, intent(in) :: ng, dg, cr, model, tile
3252 integer, intent(in) :: lbi, ubi, lbj, ubj
3253!
3254! Local variable declarations.
3255!
3256# ifdef NESTING_DEBUG
3257 logical, save :: first = .true.
3258!
3259# endif
3260 integer :: i, itrc, j, k, m, tnew, told
3261!
3262 real(dp) :: wnew, wold, secscale, fac
3263 real(r8) :: my_value, tl_my_value
3264!
3265 character (len=*), parameter :: myfile = &
3266 & __FILE__//", tl_put_refine3d"
3267
3268# include "set_bounds.h"
3269!
3270!-----------------------------------------------------------------------
3271! Interpolate (space, time) refinement grid contact points for 2D state
3272! variables from donor grid.
3273!-----------------------------------------------------------------------
3274!
3275! Set time snapshot indices for the donor grid data.
3276!
3277 told=3-rollingindex(cr)
3278 tnew=rollingindex(cr)
3279!
3280! Set linear time interpolation weights. Fractional seconds are
3281! rounded to the nearest milliseconds integer towards zero in the
3282! time interpolation weights.
3283!
3284 secscale=1000.0_dp ! seconds to milliseconds
3285!
3286 wold=anint((rollingtime(tnew,cr)-time(ng))*secscale,dp)
3287 wnew=anint((time(ng)-rollingtime(told,cr))*secscale,dp)
3288 fac=1.0_dp/(wold+wnew)
3289 wold=fac*wold
3290 wnew=fac*wnew
3291!
3292 IF (((wold*wnew).lt.0.0_dp).or.((wold+wnew).le.0.0_dp)) THEN
3293 IF (domain(ng)%SouthWest_Test(tile)) THEN
3294 IF (master) THEN
3295 WRITE (stdout,10) cr, dg, ng, &
3296 & iic(dg), told, tnew, &
3297 & iic(ng), wold, wnew, &
3298 & int(time(ng)), &
3299 & int(rollingtime(told,cr)), &
3300 & int(rollingtime(tnew,cr))
3301 END IF
3302 exit_flag=8
3303 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3304 END IF
3305 END IF
3306
3307# ifdef NESTING_DEBUG
3308!
3309! If debugging, write information into Fortran unit 201 to check the
3310! logic of interpolating from donor grid data.
3311!
3312 IF (domain(ng)%SouthWest_Test(tile)) THEN
3313 IF (master) THEN
3314 IF (first) THEN
3315 first=.false.
3316 WRITE (201,20)
3317 END IF
3318 WRITE (201,30) cr, dg, ng, iic(dg), iic(ng), told, tnew, &
3319 & int(time(dg)), &
3320 & int(rollingtime(told,cr)), &
3321 & int(time(ng)), &
3322 & int(rollingtime(tnew,cr)), &
3323 & wold, wnew
3324 20 FORMAT (3x,'cr',3x,'dg',3x,'ng',4x,'iic',4x,'iic',2x,'told', &
3325 & 2x,'tnew',7x,'time',7x,'time',7x,'time',7x,'time', &
3326 & 7x,'Wold',7x,'Wnew',/,18x,'(dg)',3x,'(ng)', &
3327 & 19x,'(dg)',7x,'told',7x,'(ng)',7x,'tnew',/)
3328 30 FORMAT (3i5,2i7,2i6,4(2x,i9),2f11.4)
3329 FLUSH (201)
3330 END IF
3331 END IF
3332# endif
3333!
3334! Tracer-type variables.
3335!
3336 DO m=1,rcontact(cr)%Npoints
3337 i=rcontact(cr)%Irg(m)
3338 j=rcontact(cr)%Jrg(m)
3339 IF (((istrt.le.i).and.(i.le.iendt)).and. &
3340 & ((jstrt.le.j).and.(j.le.jendt))) THEN
3341 DO itrc=1,nt(ng)
3342 DO k=1,n(ng)
3343!^ my_value=Wold* &
3344!^ & (Rcontact(cr)%Lweight(1,m)* &
3345!^ & REFINED(cr)%t(1,k,m,told,itrc)+ &
3346!^ & Rcontact(cr)%Lweight(2,m)* &
3347!^ & REFINED(cr)%t(2,k,m,told,itrc)+ &
3348!^ & Rcontact(cr)%Lweight(3,m)* &
3349!^ & REFINED(cr)%t(3,k,m,told,itrc)+ &
3350!^ & Rcontact(cr)%Lweight(4,m)* &
3351!^ & REFINED(cr)%t(4,k,m,told,itrc))+ &
3352!^ & Wnew* &
3353!^ & (Rcontact(cr)%Lweight(1,m)* &
3354!^ & REFINED(cr)%t(1,k,m,tnew,itrc)+ &
3355!^ & Rcontact(cr)%Lweight(2,m)* &
3356!^ & REFINED(cr)%t(2,k,m,tnew,itrc)+ &
3357!^ & Rcontact(cr)%Lweight(3,m)* &
3358!^ & REFINED(cr)%t(3,k,m,tnew,itrc)+ &
3359!^ & Rcontact(cr)%Lweight(4,m)* &
3360!^ & REFINED(cr)%t(4,k,m,tnew,itrc))
3361!^
3362 tl_my_value=wold* &
3363 & (rcontact(cr)%Lweight(1,m)* &
3364 & refined(cr)%tl_t(1,k,m,told,itrc)+ &
3365 & rcontact(cr)%Lweight(2,m)* &
3366 & refined(cr)%tl_t(2,k,m,told,itrc)+ &
3367 & rcontact(cr)%Lweight(3,m)* &
3368 & refined(cr)%tl_t(3,k,m,told,itrc)+ &
3369 & rcontact(cr)%Lweight(4,m)* &
3370 & refined(cr)%tl_t(4,k,m,told,itrc))+ &
3371 & wnew* &
3372 & (rcontact(cr)%Lweight(1,m)* &
3373 & refined(cr)%tl_t(1,k,m,tnew,itrc)+ &
3374 & rcontact(cr)%Lweight(2,m)* &
3375 & refined(cr)%tl_t(2,k,m,tnew,itrc)+ &
3376 & rcontact(cr)%Lweight(3,m)* &
3377 & refined(cr)%tl_t(3,k,m,tnew,itrc)+ &
3378 & rcontact(cr)%Lweight(4,m)* &
3379 & refined(cr)%tl_t(4,k,m,tnew,itrc))
3380# ifdef MASKING
3381!^ my_value=my_value*GRID(ng)%rmask(i,j)
3382!^
3383 tl_my_value=tl_my_value*grid(ng)%rmask(i,j)
3384# endif
3385!^ OCEAN(ng)%t(i,j,k,1,itrc)=my_value
3386!^
3387 ocean(ng)%tl_t(i,j,k,1,itrc)=tl_my_value
3388!^ OCEAN(ng)%t(i,j,k,2,itrc)=my_value
3389!^
3390 ocean(ng)%tl_t(i,j,k,2,itrc)=tl_my_value
3391!^ OCEAN(ng)%t(i,j,k,3,itrc)=my_value
3392!^
3393 ocean(ng)%tl_t(i,j,k,3,itrc)=tl_my_value
3394 END DO
3395 END DO
3396 END IF
3397 END DO
3398!
3399! 3D momentum in the XI-direction.
3400!
3401 DO m=1,ucontact(cr)%Npoints
3402 i=ucontact(cr)%Irg(m)
3403 j=ucontact(cr)%Jrg(m)
3404 IF (((istrp.le.i).and.(i.le.iendt)).and. &
3405 & ((jstrt.le.j).and.(j.le.jendt))) THEN
3406 DO k=1,n(ng)
3407!^ my_value=Wold* &
3408!^ & (Ucontact(cr)%Lweight(1,m)* &
3409!^ & REFINED(cr)%u(1,k,m,told)+ &
3410!^ & Ucontact(cr)%Lweight(2,m)* &
3411!^ & REFINED(cr)%u(2,k,m,told)+ &
3412!^ & Ucontact(cr)%Lweight(3,m)* &
3413!^ & REFINED(cr)%u(3,k,m,told)+ &
3414!^ & Ucontact(cr)%Lweight(4,m)* &
3415!^ & REFINED(cr)%u(4,k,m,told))+ &
3416!^ & Wnew* &
3417!^ & (Ucontact(cr)%Lweight(1,m)* &
3418!^ & REFINED(cr)%u(1,k,m,tnew)+ &
3419!^ & Ucontact(cr)%Lweight(2,m)* &
3420!^ & REFINED(cr)%u(2,k,m,tnew)+ &
3421!^ & Ucontact(cr)%Lweight(3,m)* &
3422!^ & REFINED(cr)%u(3,k,m,tnew)+ &
3423!^ & Ucontact(cr)%Lweight(4,m)* &
3424!^ & REFINED(cr)%u(4,k,m,tnew))
3425!^
3426 tl_my_value=wold* &
3427 & (ucontact(cr)%Lweight(1,m)* &
3428 & refined(cr)%tl_u(1,k,m,told)+ &
3429 & ucontact(cr)%Lweight(2,m)* &
3430 & refined(cr)%tl_u(2,k,m,told)+ &
3431 & ucontact(cr)%Lweight(3,m)* &
3432 & refined(cr)%tl_u(3,k,m,told)+ &
3433 & ucontact(cr)%Lweight(4,m)* &
3434 & refined(cr)%tl_u(4,k,m,told))+ &
3435 & wnew* &
3436 & (ucontact(cr)%Lweight(1,m)* &
3437 & refined(cr)%tl_u(1,k,m,tnew)+ &
3438 & ucontact(cr)%Lweight(2,m)* &
3439 & refined(cr)%tl_u(2,k,m,tnew)+ &
3440 & ucontact(cr)%Lweight(3,m)* &
3441 & refined(cr)%tl_u(3,k,m,tnew)+ &
3442 & ucontact(cr)%Lweight(4,m)* &
3443 & refined(cr)%tl_u(4,k,m,tnew))
3444# ifdef MASKING
3445!^ my_value=my_value*GRID(ng)%umask(i,j)
3446!^
3447 tl_my_value=tl_my_value*grid(ng)%umask(i,j)
3448# endif
3449!^ OCEAN(ng)%u(i,j,k,1)=my_value
3450!^
3451 ocean(ng)%tl_u(i,j,k,1)=tl_my_value
3452!^ OCEAN(ng)%u(i,j,k,2)=my_value
3453!^
3454 ocean(ng)%tl_u(i,j,k,2)=tl_my_value
3455 END DO
3456 END IF
3457 END DO
3458!
3459! 3D momentum in the ETA-direction.
3460!
3461 DO m=1,vcontact(cr)%Npoints
3462 i=vcontact(cr)%Irg(m)
3463 j=vcontact(cr)%Jrg(m)
3464 IF (((istrt.le.i).and.(i.le.iendt)).and. &
3465 & ((jstrp.le.j).and.(j.le.jendt))) THEN
3466 DO k=1,n(ng)
3467!^ my_value=Wold* &
3468!^ & (Vcontact(cr)%Lweight(1,m)* &
3469!^ & REFINED(cr)%v(1,k,m,told)+ &
3470!^ & Vcontact(cr)%Lweight(2,m)* &
3471!^ & REFINED(cr)%v(2,k,m,told)+ &
3472!^ & Vcontact(cr)%Lweight(3,m)* &
3473!^ & REFINED(cr)%v(3,k,m,told)+ &
3474!^ & Vcontact(cr)%Lweight(4,m)* &
3475!^ & REFINED(cr)%v(4,k,m,told))+ &
3476!^ & Wnew* &
3477!^ & (Vcontact(cr)%Lweight(1,m)* &
3478!^ & REFINED(cr)%v(1,k,m,tnew)+ &
3479!^ & Vcontact(cr)%Lweight(2,m)* &
3480!^ & REFINED(cr)%v(2,k,m,tnew)+ &
3481!^ & Vcontact(cr)%Lweight(3,m)* &
3482!^ & REFINED(cr)%v(3,k,m,tnew)+ &
3483!^ & Vcontact(cr)%Lweight(4,m)* &
3484!^ & REFINED(cr)%v(4,k,m,tnew))
3485!^
3486 tl_my_value=wold* &
3487 & (vcontact(cr)%Lweight(1,m)* &
3488 & refined(cr)%tl_v(1,k,m,told)+ &
3489 & vcontact(cr)%Lweight(2,m)* &
3490 & refined(cr)%tl_v(2,k,m,told)+ &
3491 & vcontact(cr)%Lweight(3,m)* &
3492 & refined(cr)%tl_v(3,k,m,told)+ &
3493 & vcontact(cr)%Lweight(4,m)* &
3494 & refined(cr)%tl_v(4,k,m,told))+ &
3495 & wnew* &
3496 & (vcontact(cr)%Lweight(1,m)* &
3497 & refined(cr)%tl_v(1,k,m,tnew)+ &
3498 & vcontact(cr)%Lweight(2,m)* &
3499 & refined(cr)%tl_v(2,k,m,tnew)+ &
3500 & vcontact(cr)%Lweight(3,m)* &
3501 & refined(cr)%tl_v(3,k,m,tnew)+ &
3502 & vcontact(cr)%Lweight(4,m)* &
3503 & refined(cr)%tl_v(4,k,m,tnew))
3504# ifdef MASKING
3505!^ my_value=my_value*GRID(ng)%vmask(i,j)
3506!^
3507 tl_my_value=tl_my_value*grid(ng)%vmask(i,j)
3508# endif
3509!^ OCEAN(ng)%v(i,j,k,1)=my_value
3510!^
3511 ocean(ng)%tl_v(i,j,k,1)=tl_my_value
3512!^ OCEAN(ng)%v(i,j,k,2)=my_value
3513!^
3514 ocean(ng)%tl_v(i,j,k,2)=tl_my_value
3515 END DO
3516 END IF
3517 END DO
3518
3519# ifdef DISTRIBUTE
3520!
3521!-----------------------------------------------------------------------
3522! Exchange tile information.
3523!-----------------------------------------------------------------------
3524!
3525!^ CALL mp_exchange4d (ng, tile, model, 3, &
3526!^ & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), &
3527!^ & NghostPoints, &
3528!^ & EWperiodic(ng), NSperiodic(ng), &
3529!^ & OCEAN(ng)%t(:,:,:,1,:), &
3530!^ & OCEAN(ng)%t(:,:,:,2,:), &
3531!^ & OCEAN(ng)%t(:,:,:,3,:))
3532!^
3533 CALL mp_exchange4d (ng, tile, model, 3, &
3534 & lbi, ubi, lbj, ubj, 1, n(ng), 1, nt(ng), &
3535 & nghostpoints, &
3536 & ewperiodic(ng), nsperiodic(ng), &
3537 & ocean(ng)%tl_t(:,:,:,1,:), &
3538 & ocean(ng)%tl_t(:,:,:,2,:), &
3539 & ocean(ng)%tl_t(:,:,:,3,:))
3540!^ CALL mp_exchange3d (ng, tile, model, 4, &
3541!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
3542!^ & NghostPoints, &
3543!^ & EWperiodic(ng), NSperiodic(ng), &
3544!^ & OCEAN(ng)%u(:,:,:,1), &
3545!^ & OCEAN(ng)%u(:,:,:,2), &
3546!^ & OCEAN(ng)%v(:,:,:,1), &
3547!^ & OCEAN(ng)%v(:,:,:,2))
3548!^
3549 CALL mp_exchange3d (ng, tile, model, 4, &
3550 & lbi, ubi, lbj, ubj, 1, n(ng), &
3551 & nghostpoints, &
3552 & ewperiodic(ng), nsperiodic(ng), &
3553 & ocean(ng)%tl_u(:,:,:,1), &
3554 & ocean(ng)%tl_u(:,:,:,2), &
3555 & ocean(ng)%tl_v(:,:,:,1), &
3556 & ocean(ng)%tl_v(:,:,:,2))
3557# endif
3558!
3559 10 FORMAT (/,' PUT_REFINE3D - unbounded contact points temporal: ', &
3560 & ' interpolation:', &
3561 & /,2x, 'cr = ',i2.2, &
3562 & 8x,'dg = ',i2.2, &
3563 & 8x,'ng = ',i2.2, &
3564 & /,2x, 'iic(dg) = ',i7.7, &
3565 & 3x,'told = ',i1, &
3566 & 9x,'tnew = ',i1, &
3567 & /,2x, 'iic(ng) = ',i7.7, &
3568 & 3x,'Wold = ',f8.5, &
3569 & 2x,'Wnew = ',f8.5, &
3570 & /,2x, 'time(ng) = ',i7.7, &
3571 & 3x,'time(told) = ',i7.7, &
3572 & 3x,'time(tnew) = ',i7.7)
3573!
3574 RETURN
3575 END SUBROUTINE tl_put_refine3d
3576# endif
3577
3578# ifdef SOLVE3D
3579!
3580 SUBROUTINE tl_z_weights (ng, model, tile)
3581!
3582!=======================================================================
3583! !
3584! This routine determines the vertical indices and interpolation !
3585! weights associated with depth, which are needed to process 3D !
3586! fields in the contact region. !
3587! !
3588! On Input: !
3589! !
3590! model Calling model identifier (integer) !
3591! tile Domain partition for composite grid ng (integer) !
3592! !
3593! On Output: Updated T_NGC type structures in mod_param: !
3594! !
3595! Rcontact Updated values for Kdg(:,:) and Vweigths (:,:,:) !
3596! Ucontact Updated values for Kdg(:,:) and Vweigths (:,:,:) !
3597! Vcontact Updated values for Kdg(:,:) and Vweigths (:,:,:) !
3598! !
3599!=======================================================================
3600!
3601 USE mod_param
3602 USE mod_grid
3603 USE mod_nesting
3604 USE mod_scalars
3605!
3606# ifdef DISTRIBUTE
3607 USE distribute_mod, ONLY : mp_assemble
3608# endif
3609 USE strings_mod, ONLY : founderror
3610!
3611! Imported variable declarations.
3612!
3613 integer, intent(in) :: ng, model, tile
3614!
3615! Local variable declarations.
3616!
3617 integer :: cr, dg, rg, i, j, k, m
3618 integer :: idg, jdg, kdg, imind, imaxd, jmind, jmaxd
3619 integer :: irg, jrg, krg, iminr, imaxr, jminr, jmaxr
3620 integer :: idgm1, idgp1, jdgm1, jdgp1
3621 integer :: npoints
3622# ifdef DISTRIBUTE
3623 integer :: nkpts, nwpts, nzpts
3624
3625 integer, parameter :: ispv = 0
3626# endif
3627!
3628 real(r8), parameter :: spv = 0.0_r8
3629
3630 real(r8) :: zbot, zr, ztop, dz, r1, r2
3631 real(r8) :: tl_zbot, tl_zr, tl_ztop, tl_dz, tl_r1, tl_r2
3632
3633 real(r8), allocatable :: zd(:,:,:)
3634 real(r8), allocatable :: tl_zd(:,:,:)
3635!
3636 character (len=*), parameter :: myfile = &
3637 & __FILE__//", tl_z_weights"
3638!
3639!=======================================================================
3640! Compute vertical indices and weights for each contact region.
3641!=======================================================================
3642!
3643 DO cr=1,ncontact
3644!
3645! Get donor and receiver grid numbers.
3646!
3647 dg=rcontact(cr)%donor_grid
3648 rg=rcontact(cr)%receiver_grid
3649!
3650! Process only contact region data for requested nested grid "ng".
3651!
3652 IF (rg.eq.ng) THEN
3653!
3654!-----------------------------------------------------------------------
3655! Process variables in structure Rcontact(cr).
3656!-----------------------------------------------------------------------
3657!
3658! Get number of contact points to process.
3659!
3660 npoints=rcontact(cr)%Npoints
3661!
3662! Set starting and ending tile indices for the donor and receiver
3663! grids.
3664!
3665 imind=bounds(dg) % IstrT(tile)
3666 imaxd=bounds(dg) % IendT(tile)
3667 jmind=bounds(dg) % JstrT(tile)
3668 jmaxd=bounds(dg) % JendT(tile)
3669!
3670 iminr=bounds(rg) % IstrT(tile)
3671 imaxr=bounds(rg) % IendT(tile)
3672 jminr=bounds(rg) % JstrT(tile)
3673 jmaxr=bounds(rg) % JendT(tile)
3674
3675# ifdef DISTRIBUTE
3676!
3677! If distributed-memory, initialize with special value (zero) to
3678! facilitate the global reduction when collecting data between all
3679! nodes.
3680!
3681 nkpts=n(rg)*npoints
3682 nwpts=2*nkpts
3683 nzpts=4*nkpts
3684
3685 rcontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
3686 rcontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
3687 rcontact(cr)%tl_Vweight(1:2,1:n(rg),1:npoints)=0.0_r8
3688# endif
3689!
3690! If coincident grids and requested, avoid vertical interpolation.
3691!
3692 r_contact : IF (.not.rcontact(cr)%interpolate.and. &
3693 & rcontact(cr)%coincident) THEN
3694 DO krg=1,n(rg)
3695 DO m=1,npoints
3696 irg=rcontact(cr)%Irg(m)
3697 jrg=rcontact(cr)%Jrg(m)
3698 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
3699 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
3700 rcontact(cr)%Kdg(krg,m)=krg
3701 rcontact(cr)%Vweight(1,krg,m)=1.0_r8
3702 rcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
3703 rcontact(cr)%Vweight(2,krg,m)=0.0_r8
3704 rcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
3705 END IF
3706 END DO
3707 END DO
3708!
3709! Otherwise, vertically interpolate because donor and receiver grids
3710! are not coincident.
3711!
3712 ELSE
3713!
3714! Allocate and initialize local working arrays.
3715!
3716 IF (.not.allocated(zd)) THEN
3717 allocate ( zd(4,n(dg),npoints) )
3718 END IF
3719 zd=spv
3720 IF (.not.allocated(tl_zd)) THEN
3721 allocate ( tl_zd(4,n(dg),npoints) )
3722 END IF
3723 tl_zd=0.0_r8
3724!
3725! Extract donor grid depths for each cell containing the receiver grid
3726! contact point. Notice that indices i+1 and j+1 are bounded to the
3727! maximum possible values in contact points at the edge of the contact
3728! region. In such cases, Lweight(1,m)=1 and Lweight(2:3,m)=0. This is
3729! done to avoid out of range errors. We need to take care of this in
3730! the adjoint code.
3731!
3732 DO kdg=1,n(dg)
3733 DO m=1,npoints
3734 idg =rcontact(cr)%Idg(m)
3735 idgp1=min(idg+1, bounds(dg)%UBi(-1))
3736 jdg =rcontact(cr)%Jdg(m)
3737 jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
3738 IF (((imind.le.idg).and.(idg.le.imaxd)).and. &
3739 & ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
3740 zd(1,kdg,m)=grid(dg)%z_r(idg ,jdg ,kdg)
3741 tl_zd(1,kdg,m)=grid(dg)%tl_z_r(idg ,jdg ,kdg)
3742 zd(2,kdg,m)=grid(dg)%z_r(idgp1,jdg ,kdg)
3743 tl_zd(2,kdg,m)=grid(dg)%tl_z_r(idgp1,jdg ,kdg)
3744 zd(3,kdg,m)=grid(dg)%z_r(idgp1,jdgp1,kdg)
3745 tl_zd(3,kdg,m)=grid(dg)%tl_z_r(idgp1,jdgp1,kdg)
3746 zd(4,kdg,m)=grid(dg)%z_r(idg ,jdgp1,kdg)
3747 tl_zd(4,kdg,m)=grid(dg)%tl_z_r(idg ,jdgp1,kdg)
3748 END IF
3749 END DO
3750 END DO
3751
3752# ifdef DISTRIBUTE
3753!
3754! Exchange data between all parallel nodes.
3755!
3756 CALL mp_assemble (dg, model, nzpts, spv, zd)
3757 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3758
3759 CALL mp_assemble (dg, model, nzpts, spv, tl_zd)
3760 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3761# endif
3762!
3763! Determine donor grid vertical indices (Kdg) and weights (Vweight)
3764! needed for the interpolation of data at the receiver grid contact
3765! points.
3766!
3767 DO krg=1,n(rg)
3768 DO m=1,npoints
3769 irg=rcontact(cr)%Irg(m)
3770 jrg=rcontact(cr)%Jrg(m)
3771 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
3772 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
3773 ztop=rcontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+ &
3774 & rcontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+ &
3775 & rcontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+ &
3776 & rcontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
3777 tl_ztop=rcontact(cr)%Lweight(1,m)*tl_zd(1,n(dg),m)+ &
3778 & rcontact(cr)%Lweight(2,m)*tl_zd(2,n(dg),m)+ &
3779 & rcontact(cr)%Lweight(3,m)*tl_zd(3,n(dg),m)+ &
3780 & rcontact(cr)%Lweight(4,m)*tl_zd(4,n(dg),m)
3781 zbot=rcontact(cr)%Lweight(1,m)*zd(1,1 ,m)+ &
3782 & rcontact(cr)%Lweight(2,m)*zd(2,1 ,m)+ &
3783 & rcontact(cr)%Lweight(3,m)*zd(3,1 ,m)+ &
3784 & rcontact(cr)%Lweight(4,m)*zd(4,1 ,m)
3785 tl_zbot=rcontact(cr)%Lweight(1,m)*tl_zd(1,1 ,m)+ &
3786 & rcontact(cr)%Lweight(2,m)*tl_zd(2,1 ,m)+ &
3787 & rcontact(cr)%Lweight(3,m)*tl_zd(3,1 ,m)+ &
3788 & rcontact(cr)%Lweight(4,m)*tl_zd(4,1 ,m)
3789 zr=grid(rg)%z_r(irg,jrg,krg)
3790 tl_zr=grid(rg)%tl_z_r(irg,jrg,krg)
3791 IF (zr.ge.ztop) THEN ! If shallower, use top
3792 rcontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
3793 rcontact(cr)%Vweight(1,krg,m)=0.0_r8
3794 rcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
3795 rcontact(cr)%Vweight(2,krg,m)=1.0_r8
3796 rcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
3797 ELSE IF (zbot.ge.zr) THEN ! If deeper, use bottom
3798 rcontact(cr)%Kdg(krg,m)=1 ! donor grid cell value
3799 rcontact(cr)%Vweight(1,krg,m)=0.0_r8
3800 rcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
3801 rcontact(cr)%Vweight(2,krg,m)=1.0_r8
3802 rcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
3803 ELSE ! bounded, interpolate
3804 DO kdg=n(dg),2,-1
3805 ztop=rcontact(cr)%Lweight(1,m)*zd(1,kdg ,m)+ &
3806 & rcontact(cr)%Lweight(2,m)*zd(2,kdg ,m)+ &
3807 & rcontact(cr)%Lweight(3,m)*zd(3,kdg ,m)+ &
3808 & rcontact(cr)%Lweight(4,m)*zd(4,kdg ,m)
3809 tl_ztop=rcontact(cr)%Lweight(1,m)* &
3810 & tl_zd(1,kdg ,m)+ &
3811 & rcontact(cr)%Lweight(2,m)* &
3812 tl_zd(2,kdg ,m)+ &
3813 & rcontact(cr)%Lweight(3,m)* &
3814 & tl_zd(3,kdg ,m)+ &
3815 & rcontact(cr)%Lweight(4,m)*tl_zd(4,kdg ,m)
3816 zbot=rcontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+ &
3817 & rcontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+ &
3818 & rcontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+ &
3819 & rcontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
3820 tl_zbot=rcontact(cr)%Lweight(1,m)* &
3821 & tl_zd(1,kdg-1,m)+ &
3822 & rcontact(cr)%Lweight(2,m)* &
3823 & tl_zd(2,kdg-1,m)+ &
3824 & rcontact(cr)%Lweight(3,m)* &
3825 & tl_zd(3,kdg-1,m)+ &
3826 & rcontact(cr)%Lweight(4,m)*tl_zd(4,kdg-1,m)
3827 IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
3828 dz=ztop-zbot
3829 tl_dz=tl_ztop-tl_zbot
3830 r2=(zr-zbot)/dz
3831 tl_r2=(tl_zr-tl_zbot)/dz-tl_dz*r2/dz
3832 r1=1.0_r8-r2
3833 tl_r1=-tl_r2
3834 rcontact(cr)%Kdg(krg,m)=kdg
3835 rcontact(cr)%Vweight(1,krg,m)=r1
3836 rcontact(cr)%tl_Vweight(1,krg,m)=tl_r1
3837 rcontact(cr)%Vweight(2,krg,m)=r2
3838 rcontact(cr)%tl_Vweight(2,krg,m)=tl_r2
3839 END IF
3840 END DO
3841 END IF
3842 END IF
3843 END DO
3844 END DO
3845 END IF r_contact
3846
3847# ifdef DISTRIBUTE
3848!
3849! Exchange data between all parallel nodes.
3850!
3851 CALL mp_assemble (rg, model, nkpts, ispv, rcontact(cr)%Kdg)
3852 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3853
3854 CALL mp_assemble (rg, model, nwpts, spv, rcontact(cr)%Vweight)
3855 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3856
3857 CALL mp_assemble (rg, model, nwpts, spv, &
3858 & rcontact(cr)%tl_Vweight)
3859 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3860# endif
3861!
3862! Deallocate local work arrays.
3863!
3864 IF (allocated(zd)) THEN
3865 deallocate (zd)
3866 END IF
3867 IF (allocated(tl_zd)) THEN
3868 deallocate (tl_zd)
3869 END IF
3870!
3871!-----------------------------------------------------------------------
3872! Process variables in structure Ucontact(cr).
3873!-----------------------------------------------------------------------
3874!
3875! Get number of contact points to process.
3876!
3877 npoints=ucontact(cr)%Npoints
3878!
3879! Set starting and ending tile indices for the donor and receiver
3880! grids.
3881!
3882 imind=bounds(dg) % IstrP(tile)
3883 imaxd=bounds(dg) % IendT(tile)
3884 jmind=bounds(dg) % JstrT(tile)
3885 jmaxd=bounds(dg) % JendT(tile)
3886!
3887 iminr=bounds(rg) % IstrP(tile)
3888 imaxr=bounds(rg) % IendT(tile)
3889 jminr=bounds(rg) % JstrT(tile)
3890 jmaxr=bounds(rg) % JendT(tile)
3891
3892# ifdef DISTRIBUTE
3893!
3894! If distributed-memory, initialize with special value (zero) to
3895! facilitate the global reduction when collecting data between all
3896! nodes.
3897!
3898 nkpts=n(rg)*npoints
3899 nwpts=2*nkpts
3900 nzpts=4*nkpts
3901
3902 ucontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
3903 ucontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
3904 ucontact(cr)%tl_Vweight(1:2,1:n(rg),1:npoints)=0.0_r8
3905# endif
3906!
3907! If coincident grids and requested, avoid vertical interpolation.
3908!
3909 u_contact : IF (.not.ucontact(cr)%interpolate.and. &
3910 & ucontact(cr)%coincident) THEN
3911 DO krg=1,n(rg)
3912 DO m=1,npoints
3913 irg=ucontact(cr)%Irg(m)
3914 jrg=ucontact(cr)%Jrg(m)
3915 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
3916 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
3917 ucontact(cr)%Kdg(krg,m)=krg
3918 ucontact(cr)%Vweight(1,krg,m)=1.0_r8
3919 ucontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
3920 ucontact(cr)%Vweight(2,krg,m)=0.0_r8
3921 ucontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
3922 END IF
3923 END DO
3924 END DO
3925!
3926! Otherwise, vertically interpolate because donor and receiver grids
3927! are not coincident.
3928!
3929 ELSE
3930!
3931! Allocate and initialize local working arrays.
3932!
3933 IF (.not.allocated(zd)) THEN
3934 allocate (zd(4,n(dg),npoints))
3935 END IF
3936 zd=spv
3937 IF (.not.allocated(tl_zd)) THEN
3938 allocate (tl_zd(4,n(dg),npoints))
3939 END IF
3940 tl_zd=0.0_r8
3941!
3942! Extract donor grid depths for each cell containing the receiver grid
3943! contact point. Notice that indices i-1, i+1 and j-1, j+1 are bounded
3944! the minimum/maximum possible values in contact points at the edge of
3945! the contact region. In such cases, the interpolation weights
3946! Lweight(1,m)=1 and Lweight(2:3,m)=0. This is done to avoid out of
3947! range errors. We need to take care of this in the adjoint code.
3948!
3949 DO kdg=1,n(dg)
3950 DO m=1,npoints
3951 idg =ucontact(cr)%Idg(m)
3952 idgm1=max(idg-1, bounds(dg)%LBi(-1))
3953 idgp1=min(idg+1, bounds(dg)%UBi(-1))
3954 jdg =ucontact(cr)%Jdg(m)
3955 jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
3956 IF (((imind.le.idg).and.(idg.le.imaxd)).and. &
3957 & ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
3958 zd(1,kdg,m)=0.5_r8*(grid(dg)%z_r(idgm1,jdg ,kdg)+ &
3959 & grid(dg)%z_r(idg ,jdg ,kdg))
3960 tl_zd(1,kdg,m)=0.5_r8* &
3961 & (grid(dg)%tl_z_r(idgm1,jdg ,kdg)+ &
3962 & grid(dg)%tl_z_r(idg ,jdg ,kdg))
3963 zd(2,kdg,m)=0.5_r8*(grid(dg)%z_r(idg ,jdg ,kdg)+ &
3964 & grid(dg)%z_r(idgp1,jdg ,kdg))
3965 tl_zd(2,kdg,m)=0.5_r8* &
3966 & (grid(dg)%tl_z_r(idg ,jdg ,kdg)+ &
3967 & grid(dg)%tl_z_r(idgp1,jdg ,kdg))
3968 zd(3,kdg,m)=0.5_r8*(grid(dg)%z_r(idg ,jdgp1,kdg)+ &
3969 & grid(dg)%z_r(idgp1,jdgp1,kdg))
3970 tl_zd(3,kdg,m)=0.5_r8* &
3971 & (grid(dg)%tl_z_r(idg ,jdgp1,kdg)+ &
3972 & grid(dg)%tl_z_r(idgp1,jdgp1,kdg))
3973 zd(4,kdg,m)=0.5_r8*(grid(dg)%z_r(idgm1,jdgp1,kdg)+ &
3974 & grid(dg)%z_r(idg ,jdgp1,kdg))
3975 tl_zd(4,kdg,m)=0.5_r8* &
3976 & (grid(dg)%tl_z_r(idgm1,jdgp1,kdg)+ &
3977 & grid(dg)%tl_z_r(idg ,jdgp1,kdg))
3978 END IF
3979 END DO
3980 END DO
3981
3982# ifdef DISTRIBUTE
3983!
3984! Exchange data between all parallel nodes.
3985!
3986 CALL mp_assemble (dg, model, nzpts, spv, zd)
3987 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3988
3989 CALL mp_assemble (dg, model, nzpts, spv, tl_zd)
3990 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3991# endif
3992!
3993! Determine donor grid vertical indices (Kdg) and weights (Vweight)
3994! needed for the interpolation of data at the receiver grid contact
3995! points.
3996!
3997 DO krg=1,n(rg)
3998 DO m=1,npoints
3999 irg=ucontact(cr)%Irg(m)
4000 jrg=ucontact(cr)%Jrg(m)
4001 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
4002 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
4003 ztop=ucontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+ &
4004 & ucontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+ &
4005 & ucontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+ &
4006 & ucontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
4007 tl_ztop=ucontact(cr)%Lweight(1,m)*tl_zd(1,n(dg),m)+ &
4008 & ucontact(cr)%Lweight(2,m)*tl_zd(2,n(dg),m)+ &
4009 & ucontact(cr)%Lweight(3,m)*tl_zd(3,n(dg),m)+ &
4010 & ucontact(cr)%Lweight(4,m)*tl_zd(4,n(dg),m)
4011 zbot=ucontact(cr)%Lweight(1,m)*zd(1,1 ,m)+ &
4012 & ucontact(cr)%Lweight(2,m)*zd(2,1 ,m)+ &
4013 & ucontact(cr)%Lweight(3,m)*zd(3,1 ,m)+ &
4014 & ucontact(cr)%Lweight(4,m)*zd(4,1 ,m)
4015 tl_zbot=ucontact(cr)%Lweight(1,m)*tl_zd(1,1 ,m)+ &
4016 & ucontact(cr)%Lweight(2,m)*tl_zd(2,1 ,m)+ &
4017 & ucontact(cr)%Lweight(3,m)*tl_zd(3,1 ,m)+ &
4018 & ucontact(cr)%Lweight(4,m)*tl_zd(4,1 ,m)
4019 zr=0.5_r8*(grid(rg)%z_r(irg ,jrg,krg)+ &
4020 & grid(rg)%z_r(irg-1,jrg,krg))
4021 tl_zr=0.5_r8*(grid(rg)%tl_z_r(irg ,jrg,krg)+ &
4022 & grid(rg)%tl_z_r(irg-1,jrg,krg))
4023 IF (zr.ge.ztop) THEN ! If shallower, use top
4024 ucontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
4025 ucontact(cr)%Vweight(1,krg,m)=0.0_r8
4026 ucontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
4027 ucontact(cr)%Vweight(2,krg,m)=1.0_r8
4028 ucontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
4029 ELSE IF (zbot.ge.zr) THEN ! If deeper, use bottom
4030 ucontact(cr)%Kdg(krg,m)=1 ! donor grid cell value
4031 ucontact(cr)%Vweight(1,krg,m)=0.0_r8
4032 ucontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
4033 ucontact(cr)%Vweight(2,krg,m)=1.0_r8
4034 ucontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
4035 ELSE ! bounded, interpolate
4036 DO kdg=n(dg),2,-1
4037 ztop=ucontact(cr)%Lweight(1,m)*zd(1,kdg ,m)+ &
4038 & ucontact(cr)%Lweight(2,m)*zd(2,kdg ,m)+ &
4039 & ucontact(cr)%Lweight(3,m)*zd(3,kdg ,m)+ &
4040 & ucontact(cr)%Lweight(4,m)*zd(4,kdg ,m)
4041 tl_ztop=ucontact(cr)%Lweight(1,m)* &
4042 & tl_zd(1,kdg ,m)+ &
4043 & ucontact(cr)%Lweight(2,m)* &
4044 & tl_zd(2,kdg ,m)+ &
4045 & ucontact(cr)%Lweight(3,m)* &
4046 & tl_zd(3,kdg ,m)+ &
4047 & ucontact(cr)%Lweight(4,m)*tl_zd(4,kdg ,m)
4048 zbot=ucontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+ &
4049 & ucontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+ &
4050 & ucontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+ &
4051 & ucontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
4052 tl_zbot=ucontact(cr)%Lweight(1,m)* &
4053 & tl_zd(1,kdg-1,m)+ &
4054 & ucontact(cr)%Lweight(2,m)* &
4055 & tl_zd(2,kdg-1,m)+ &
4056 & ucontact(cr)%Lweight(3,m)* &
4057 & tl_zd(3,kdg-1,m)+ &
4058 & ucontact(cr)%Lweight(4,m)*tl_zd(4,kdg-1,m)
4059 IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
4060 dz=ztop-zbot
4061 tl_dz=tl_ztop-tl_zbot
4062 r2=(zr-zbot)/dz
4063 tl_r2=(tl_zr-tl_zbot)/dz-tl_dz*r2/dz
4064 r1=1.0_r8-r2
4065 tl_r1=-tl_r2
4066 ucontact(cr)%Kdg(krg,m)=kdg
4067 ucontact(cr)%Vweight(1,krg,m)=r1
4068 ucontact(cr)%tl_Vweight(1,krg,m)=tl_r1
4069 ucontact(cr)%Vweight(2,krg,m)=r2
4070 ucontact(cr)%tl_Vweight(2,krg,m)=tl_r2
4071 END IF
4072 END DO
4073 END IF
4074 END IF
4075 END DO
4076 END DO
4077 END IF u_contact
4078
4079# ifdef DISTRIBUTE
4080!
4081! Exchange data between all parallel nodes.
4082!
4083 CALL mp_assemble (rg, model, nkpts, ispv, ucontact(cr)%Kdg)
4084 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4085
4086 CALL mp_assemble (rg, model, nwpts, spv, ucontact(cr)%Vweight)
4087 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4088
4089 CALL mp_assemble (rg, model, nwpts, spv, &
4090 & ucontact(cr)%tl_Vweight)
4091 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4092# endif
4093!
4094! Deallocate local work arrays.
4095!
4096 IF (allocated(zd)) THEN
4097 deallocate (zd)
4098 END IF
4099 IF (allocated(tl_zd)) THEN
4100 deallocate (tl_zd)
4101 END IF
4102!
4103!-----------------------------------------------------------------------
4104! Process variables in structure Vcontact(cr).
4105!-----------------------------------------------------------------------
4106!
4107! Get number of contact points to process.
4108!
4109 npoints=vcontact(cr)%Npoints
4110!
4111! Set starting and ending tile indices for the donor and receiver
4112! grids.
4113!
4114 imind=bounds(dg) % IstrT(tile)
4115 imaxd=bounds(dg) % IendT(tile)
4116 jmind=bounds(dg) % JstrP(tile)
4117 jmaxd=bounds(dg) % JendT(tile)
4118!
4119 iminr=bounds(rg) % IstrT(tile)
4120 imaxr=bounds(rg) % IendT(tile)
4121 jminr=bounds(rg) % JstrP(tile)
4122 jmaxr=bounds(rg) % JendT(tile)
4123
4124# ifdef DISTRIBUTE
4125!
4126! If distributed-memory, initialize with special value (zero) to
4127! facilitate the global reduction when collecting data between all
4128! nodes.
4129!
4130 nkpts=n(rg)*npoints
4131 nwpts=2*nkpts
4132 nzpts=4*nkpts
4133
4134 vcontact(cr)%Kdg(1:n(rg),1:npoints)=ispv
4135 vcontact(cr)%Vweight(1:2,1:n(rg),1:npoints)=spv
4136 vcontact(cr)%tl_Vweight(1:2,1:n(rg),1:npoints)=0.0_r8
4137# endif
4138!
4139! If coincident grids and requested, avoid vertical interpolation.
4140!
4141 v_contact : IF (.not.vcontact(cr)%interpolate.and. &
4142 & vcontact(cr)%coincident) THEN
4143 DO krg=1,n(rg)
4144 DO m=1,npoints
4145 irg=vcontact(cr)%Irg(m)
4146 jrg=vcontact(cr)%Jrg(m)
4147 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
4148 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
4149 vcontact(cr)%Kdg(krg,m)=krg
4150 vcontact(cr)%Vweight(1,krg,m)=1.0_r8
4151 vcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
4152 vcontact(cr)%Vweight(2,krg,m)=0.0_r8
4153 vcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
4154 END IF
4155 END DO
4156 END DO
4157!
4158! Otherwise, vertically interpolate because donor and receiver grids
4159! are not coincident.
4160!
4161 ELSE
4162!
4163! Allocate and initialize local working arrays.
4164!
4165 IF (.not.allocated(zd)) THEN
4166 allocate (zd(4,n(dg),npoints))
4167 END IF
4168 zd=spv
4169 IF (.not.allocated(tl_zd)) THEN
4170 allocate (tl_zd(4,n(dg),npoints))
4171 END IF
4172 tl_zd=0.0_r8
4173!
4174! Extract donor grid depths for each cell containing the receiver grid
4175! contact point.
4176!
4177 DO kdg=1,n(dg)
4178 DO m=1,npoints
4179 idg=vcontact(cr)%Idg(m)
4180 idgp1=min(idg+1, bounds(dg)%UBi(-1))
4181 jdg=vcontact(cr)%Jdg(m)
4182 jdgm1=max(jdg-1, bounds(dg)%LBj(-1))
4183 jdgp1=min(jdg+1, bounds(dg)%UBj(-1))
4184 IF (((imind.le.idg).and.(idg.le.imaxd)).and. &
4185 & ((jmind.le.jdg).and.(jdg.le.jmaxd))) THEN
4186 zd(1,kdg,m)=0.5_r8*(grid(dg)%z_r(idg ,jdgm1,kdg)+ &
4187 & grid(dg)%z_r(idg ,jdg ,kdg))
4188 tl_zd(1,kdg,m)=0.5_r8* &
4189 & (grid(dg)%tl_z_r(idg ,jdgm1,kdg)+ &
4190 & grid(dg)%tl_z_r(idg ,jdg ,kdg))
4191 zd(2,kdg,m)=0.5_r8*(grid(dg)%z_r(idgp1,jdgm1,kdg)+ &
4192 & grid(dg)%z_r(idgp1,jdg ,kdg))
4193 tl_zd(2,kdg,m)=0.5_r8* &
4194 & (grid(dg)%tl_z_r(idgp1,jdgm1,kdg)+ &
4195 & grid(dg)%tl_z_r(idgp1,jdg ,kdg))
4196 zd(3,kdg,m)=0.5_r8*(grid(dg)%z_r(idgp1,jdg ,kdg)+ &
4197 & grid(dg)%z_r(idgp1,jdgp1,kdg))
4198 tl_zd(3,kdg,m)=0.5_r8* &
4199 & (grid(dg)%tl_z_r(idgp1,jdg ,kdg)+ &
4200 & grid(dg)%tl_z_r(idgp1,jdgp1,kdg))
4201 zd(4,kdg,m)=0.5_r8*(grid(dg)%z_r(idg ,jdg ,kdg)+ &
4202 & grid(dg)%z_r(idg ,jdgp1,kdg))
4203 tl_zd(4,kdg,m)=0.5_r8* &
4204 & (grid(dg)%tl_z_r(idg ,jdg ,kdg)+ &
4205 & grid(dg)%tl_z_r(idg ,jdgp1,kdg))
4206 END IF
4207 END DO
4208 END DO
4209
4210# ifdef DISTRIBUTE
4211!
4212! Exchange data between all parallel nodes.
4213!
4214 CALL mp_assemble (dg, model, nzpts, spv, zd)
4215 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4216
4217 CALL mp_assemble (dg, model, nzpts, spv, tl_zd)
4218 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4219# endif
4220!
4221! Determine donor grid vertical indices (Kdg) and weights (Vweight)
4222! needed for the interpolation of data at the receiver grid contact
4223! points.
4224!
4225 DO krg=1,n(rg)
4226 DO m=1,npoints
4227 irg=vcontact(cr)%Irg(m)
4228 jrg=vcontact(cr)%Jrg(m)
4229 IF (((iminr.le.irg).and.(irg.le.imaxr)).and. &
4230 & ((jminr.le.jrg).and.(jrg.le.jmaxr))) THEN
4231 ztop=vcontact(cr)%Lweight(1,m)*zd(1,n(dg),m)+ &
4232 & vcontact(cr)%Lweight(2,m)*zd(2,n(dg),m)+ &
4233 & vcontact(cr)%Lweight(3,m)*zd(3,n(dg),m)+ &
4234 & vcontact(cr)%Lweight(4,m)*zd(4,n(dg),m)
4235 tl_ztop=vcontact(cr)%Lweight(1,m)*tl_zd(1,n(dg),m)+ &
4236 & vcontact(cr)%Lweight(2,m)*tl_zd(2,n(dg),m)+ &
4237 & vcontact(cr)%Lweight(3,m)*tl_zd(3,n(dg),m)+ &
4238 & vcontact(cr)%Lweight(4,m)*tl_zd(4,n(dg),m)
4239 zbot=vcontact(cr)%Lweight(1,m)*zd(1,1 ,m)+ &
4240 & vcontact(cr)%Lweight(2,m)*zd(2,1 ,m)+ &
4241 & vcontact(cr)%Lweight(3,m)*zd(3,1 ,m)+ &
4242 & vcontact(cr)%Lweight(4,m)*zd(4,1 ,m)
4243 tl_zbot=vcontact(cr)%Lweight(1,m)*tl_zd(1,1 ,m)+ &
4244 & vcontact(cr)%Lweight(2,m)*tl_zd(2,1 ,m)+ &
4245 & vcontact(cr)%Lweight(3,m)*tl_zd(3,1 ,m)+ &
4246 & vcontact(cr)%Lweight(4,m)*tl_zd(4,1 ,m)
4247 zr=0.5_r8*(grid(rg)%z_r(irg,jrg ,krg)+ &
4248 & grid(rg)%z_r(irg,jrg-1,krg))
4249 tl_zr=0.5_r8*(grid(rg)%tl_z_r(irg,jrg ,krg)+ &
4250 & grid(rg)%tl_z_r(irg,jrg-1,krg))
4251 IF (zr.ge.ztop) THEN ! If shallower, use top
4252 vcontact(cr)%Kdg(krg,m)=n(dg)! donor grid cell value
4253 vcontact(cr)%Vweight(1,krg,m)=0.0_r8
4254 vcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
4255 vcontact(cr)%Vweight(2,krg,m)=1.0_r8
4256 vcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
4257 ELSE IF (zbot.ge.zr) THEN ! If deeper, use bottom
4258 vcontact(cr)%Kdg(krg,m)=1 ! donor grid cell value
4259 vcontact(cr)%Vweight(1,krg,m)=0.0_r8
4260 vcontact(cr)%tl_Vweight(1,krg,m)=0.0_r8
4261 vcontact(cr)%Vweight(2,krg,m)=1.0_r8
4262 vcontact(cr)%tl_Vweight(2,krg,m)=0.0_r8
4263 ELSE ! bounded, interpolate
4264 DO kdg=n(dg),2,-1
4265 ztop=vcontact(cr)%Lweight(1,m)*zd(1,kdg ,m)+ &
4266 & vcontact(cr)%Lweight(2,m)*zd(2,kdg ,m)+ &
4267 & vcontact(cr)%Lweight(3,m)*zd(3,kdg ,m)+ &
4268 & vcontact(cr)%Lweight(4,m)*zd(4,kdg ,m)
4269 tl_ztop=vcontact(cr)%Lweight(1,m)* &
4270 & tl_zd(1,kdg ,m)+ &
4271 & vcontact(cr)%Lweight(2,m)* &
4272 & tl_zd(2,kdg ,m)+ &
4273 & vcontact(cr)%Lweight(3,m)* &
4274 & tl_zd(3,kdg ,m)+ &
4275 & vcontact(cr)%Lweight(4,m)*tl_zd(4,kdg ,m)
4276 zbot=vcontact(cr)%Lweight(1,m)*zd(1,kdg-1,m)+ &
4277 & vcontact(cr)%Lweight(2,m)*zd(2,kdg-1,m)+ &
4278 & vcontact(cr)%Lweight(3,m)*zd(3,kdg-1,m)+ &
4279 & vcontact(cr)%Lweight(4,m)*zd(4,kdg-1,m)
4280 tl_zbot=vcontact(cr)%Lweight(1,m)* &
4281 & tl_zd(1,kdg-1,m)+ &
4282 & vcontact(cr)%Lweight(2,m)* &
4283 & tl_zd(2,kdg-1,m)+ &
4284 & vcontact(cr)%Lweight(3,m)* &
4285 & tl_zd(3,kdg-1,m)+ &
4286 & vcontact(cr)%Lweight(4,m)*tl_zd(4,kdg-1,m)
4287 IF ((ztop.gt.zr).and.(zr.ge.zbot)) THEN
4288 dz=ztop-zbot
4289 tl_dz=tl_ztop-tl_zbot
4290 r2=(zr-zbot)/dz
4291 tl_r2=(tl_zr-tl_zbot)/dz-tl_dz*r2/dz
4292 r1=1.0_r8-r2
4293 tl_r1=-tl_r2
4294 vcontact(cr)%Kdg(krg,m)=kdg
4295 vcontact(cr)%Vweight(1,krg,m)=r1
4296 vcontact(cr)%tl_Vweight(1,krg,m)=tl_r1
4297 vcontact(cr)%Vweight(2,krg,m)=r2
4298 vcontact(cr)%tl_Vweight(2,krg,m)=tl_r2
4299 END IF
4300 END DO
4301 END IF
4302 END IF
4303 END DO
4304 END DO
4305 END IF v_contact
4306
4307# ifdef DISTRIBUTE
4308!
4309! Exchange data between all parallel nodes.
4310!
4311 CALL mp_assemble (rg, model, nkpts, ispv, vcontact(cr)%Kdg)
4312 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4313
4314 CALL mp_assemble (rg, model, nwpts, spv, vcontact(cr)%Vweight)
4315 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4316
4317 CALL mp_assemble (rg, model, nwpts, spv, &
4318 & vcontact(cr)%tl_Vweight)
4319 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4320# endif
4321!
4322! Deallocate local work arrays.
4323!
4324 IF (allocated(zd)) THEN
4325 deallocate (zd)
4326 END IF
4327 IF (allocated(tl_zd)) THEN
4328 deallocate (tl_zd)
4329 END IF
4330
4331 END IF
4332 END DO
4333!
4334 RETURN
4335 END SUBROUTINE tl_z_weights
4336# endif
4337
4338# ifdef SOLVE3D
4339!
4340 SUBROUTINE tl_put_contact3d (rg, model, tile, &
4341 & gtype, svname, &
4342 & cr, Npoints, contact, &
4343 & LBi, UBi, LBj, UBj, LBk, UBk, &
4344# ifdef MASKING
4345 & Amask, &
4346# endif
4347 & Ac, tl_Ac, tl_Ar)
4348!
4349!=======================================================================
4350! !
4351! This routine uses extracted donor grid data (Ac) to spatially !
4352! interpolate a 3D state variable at the receiver grid contact !
4353! points. If the donor and receiver grids are concident, the !
4354! Lweight(1,:) is unity and Lweight(2:4,:) are zero. !
4355! !
4356! On Input: !
4357! !
4358! rg Receiver grid number (integer) !
4359! model Calling model identifier (integer) !
4360! tile Domain tile partition (integer) !
4361! gtype C-grid variable type (integer) !
4362! svname State variable name (string) !
4363! cr Contact region number to process (integer) !
4364! Npoints Number of points in the contact region (integer) !
4365! contact Contact region information variables (T_NGC structure)!
4366! LBi Receiver grid, I-dimension Lower bound (integer) !
4367! UBi Receiver grid, I-dimension Upper bound (integer) !
4368! LBj Receiver grid, J-dimension Lower bound (integer) !
4369! UBj Receiver grid, J-dimension Upper bound (integer) !
4370! LBk Receiver grid, K-dimension Lower bound (integer) !
4371! UBk Receiver grid, K-dimension Upper bound (integer) !
4372# ifdef MASKING
4373! Amask Receiver grid land/sea masking !
4374# endif
4375! Ac Contact point data extracted from donor grid !
4376! !
4377! On Output: !
4378! !
4379! Ar Updated receiver grid 3D state array !
4380! !
4381!=======================================================================
4382!
4383 USE mod_param
4384 USE mod_ncparam
4385 USE mod_nesting
4386!
4387! Imported variable declarations.
4388!
4389 integer, intent(in) :: rg, model, tile
4390 integer, intent(in) :: gtype, cr, Npoints
4391 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
4392!
4393 character(len=*), intent(in) :: svname
4394!
4395 TYPE (T_NGC), intent(in) :: contact(:)
4396!
4397# ifdef ASSUMED_SHAPE
4398 real(r8), intent(in) :: Ac(:,:,:)
4399 real(r8), intent(in) :: tl_Ac(:,:,:)
4400# ifdef MASKING
4401 real(r8), intent(in) :: Amask(LBi:,LBj:)
4402# endif
4403 real(r8), intent(inout) :: tl_Ar(LBi:,LBj:,LBk:)
4404# else
4405 real(r8), intent(in) :: Ac(Npoints,LBk:UBk,4)
4406 real(r8), intent(in) :: tl_Ac(Npoints,LBk:UBk,4)
4407# ifdef MASKING
4408 real(r8), intent(in) :: Amask(LBi:UBi,LBj:UBj)
4409# endif
4410 real(r8), intent(inout) :: tl_Ar(LBi:UBi,LBj:UBj,LBk:UBk)
4411# endif
4412!
4413! Local variable declarations.
4414!
4415 integer :: i, j, k, kdg, kdgm1, m
4416 integer :: Istr, Iend, Jstr, Jend, Kmin
4417
4418 real(r8), dimension(8) :: cff
4419 real(r8), dimension(8) :: tl_cff
4420!
4421!-----------------------------------------------------------------------
4422! Interpolate 3D data from donor grid to receiver grid contact points.
4423!-----------------------------------------------------------------------
4424!
4425! Set starting and ending tile indices for the receiver grid.
4426!
4427 SELECT CASE (gtype)
4428 CASE (r3dvar)
4429 istr=bounds(rg) % IstrT(tile)
4430 iend=bounds(rg) % IendT(tile)
4431 jstr=bounds(rg) % JstrT(tile)
4432 jend=bounds(rg) % JendT(tile)
4433 kmin=1
4434 CASE (u3dvar)
4435 istr=bounds(rg) % IstrP(tile)
4436 iend=bounds(rg) % IendT(tile)
4437 jstr=bounds(rg) % JstrT(tile)
4438 jend=bounds(rg) % JendT(tile)
4439 kmin=1
4440 CASE (v3dvar)
4441 istr=bounds(rg) % IstrT(tile)
4442 iend=bounds(rg) % IendT(tile)
4443 jstr=bounds(rg) % JstrP(tile)
4444 jend=bounds(rg) % JendT(tile)
4445 kmin=1
4446 CASE (w3dvar)
4447 istr=bounds(rg) % IstrT(tile)
4448 iend=bounds(rg) % IendT(tile)
4449 jstr=bounds(rg) % JstrT(tile)
4450 jend=bounds(rg) % JendT(tile)
4451 kmin=0
4452 END SELECT
4453!
4454! Interpolate.
4455!
4456 DO k=lbk,ubk
4457 DO m=1,npoints
4458 i=contact(cr)%Irg(m)
4459 j=contact(cr)%Jrg(m)
4460 kdg=contact(cr)%Kdg(k,m)
4461 kdgm1=max(kdg-1,kmin)
4462 IF (((istr.le.i).and.(i.le.iend)).and. &
4463 & ((jstr.le.j).and.(j.le.jend))) THEN
4464 cff(1)=contact(cr)%Lweight(1,m)*contact(cr)%Vweight(1,k,m)
4465 tl_cff(1)=contact(cr)%Lweight(1,m)* &
4466 & contact(cr)%tl_Vweight(1,k,m)
4467 cff(2)=contact(cr)%Lweight(2,m)*contact(cr)%Vweight(1,k,m)
4468 tl_cff(2)=contact(cr)%Lweight(2,m)* &
4469 & contact(cr)%tl_Vweight(1,k,m)
4470 cff(3)=contact(cr)%Lweight(3,m)*contact(cr)%Vweight(1,k,m)
4471 tl_cff(3)=contact(cr)%Lweight(3,m)* &
4472 & contact(cr)%tl_Vweight(1,k,m)
4473 cff(4)=contact(cr)%Lweight(4,m)*contact(cr)%Vweight(1,k,m)
4474 tl_cff(4)=contact(cr)%Lweight(4,m)* &
4475 & contact(cr)%tl_Vweight(1,k,m)
4476 cff(5)=contact(cr)%Lweight(1,m)*contact(cr)%Vweight(2,k,m)
4477 tl_cff(5)=contact(cr)%Lweight(1,m)* &
4478 & contact(cr)%tl_Vweight(2,k,m)
4479 cff(6)=contact(cr)%Lweight(2,m)*contact(cr)%Vweight(2,k,m)
4480 tl_cff(6)=contact(cr)%Lweight(2,m)* &
4481 & contact(cr)%tl_Vweight(2,k,m)
4482 cff(7)=contact(cr)%Lweight(3,m)*contact(cr)%Vweight(2,k,m)
4483 tl_cff(7)=contact(cr)%Lweight(3,m)* &
4484 & contact(cr)%tl_Vweight(2,k,m)
4485 cff(8)=contact(cr)%Lweight(4,m)*contact(cr)%Vweight(2,k,m)
4486 tl_cff(8)=contact(cr)%Lweight(4,m)* &
4487 & contact(cr)%tl_Vweight(2,k,m)
4488!^ Ar(i,j,k)=cff(1)*Ac(1,kdgm1,m)+ &
4489!^ & cff(2)*Ac(2,kdgm1,m)+ &
4490!^ & cff(3)*Ac(3,kdgm1,m)+ &
4491!^ & cff(4)*Ac(4,kdgm1,m)+ &
4492!^ & cff(5)*Ac(1,kdg ,m)+ &
4493!^ & cff(6)*Ac(2,kdg ,m)+ &
4494!^ & cff(7)*Ac(3,kdg ,m)+ &
4495!^ & cff(8)*Ac(4,kdg ,m)
4496 tl_ar(i,j,k)=tl_cff(1)*ac(1,kdgm1,m)+ &
4497 & cff(1)*tl_ac(1,kdgm1,m)+ &
4498 & tl_cff(2)*ac(2,kdgm1,m)+ &
4499 & cff(2)*tl_ac(2,kdgm1,m)+ &
4500 & tl_cff(3)*ac(3,kdgm1,m)+ &
4501 & cff(3)*tl_ac(3,kdgm1,m)+ &
4502 & tl_cff(4)*ac(4,kdgm1,m)+ &
4503 & cff(4)*tl_ac(4,kdgm1,m)+ &
4504 & tl_cff(5)*ac(1,kdg ,m)+ &
4505 & cff(5)*tl_ac(1,kdg ,m)+ &
4506 & tl_cff(6)*ac(2,kdg ,m)+ &
4507 & cff(6)*tl_ac(2,kdg ,m)+ &
4508 & tl_cff(7)*ac(3,kdg ,m)+ &
4509 & cff(7)*tl_ac(3,kdg ,m)+ &
4510 & tl_cff(8)*ac(4,kdg ,m)+ &
4511 & cff(8)*tl_ac(4,kdg ,m)
4512# ifdef MASKING
4513!^ Ar(i,j,k)=Ar(i,j,k)*Amask(i,j)
4514 tl_ar(i,j,k)=tl_ar(i,j,k)*amask(i,j)
4515# endif
4516 END IF
4517 END DO
4518 END DO
4519!
4520 RETURN
4521 END SUBROUTINE tl_put_contact3d
4522# endif
4523 SUBROUTINE tl_check_massflux (ngf, model, tile)
4524!
4525!=======================================================================
4526! !
4527! If refinement, this routine check mass fluxes between coarse and !
4528! fine grids for mass and volume conservation. It is only used for !
4529! diagnostic purposes. !
4530! !
4531! On Input: !
4532! !
4533! ngf Finer grid number (integer) !
4534! model Calling model identifier (integer) !
4535! tile Domain tile partition (integer) !
4536! !
4537! On Output: (mod_nesting) !
4538! !
4539! BRY_CONTACT Updated Mflux in structure. !
4540! !
4541!=======================================================================
4542!
4543 USE mod_param
4544 USE mod_parallel
4545 USE mod_nesting
4546 USE mod_scalars
4547
4548# ifdef DISTRIBUTE
4549!
4550 USE distribute_mod, ONLY : mp_assemble
4551# endif
4552!
4553! Imported variable declarations.
4554!
4555 integer, intent(in) :: ngf, model, tile
4556!
4557! Local variable declarations.
4558!
4559# ifdef DISTRIBUTE
4560 integer :: ilb, iub, jlb, jub, nptssn, nptswe, my_tile
4561# endif
4562 integer :: iedge, ibc, ibc_min, ibc_max, ibf, io
4563 integer :: jedge, jbc, jbc_min, jbc_max, jbf, jo
4564 integer :: istr, iend, jstr, jend
4565 integer :: cjcr, cr, dg, half, icr, isum, jsum, m, rg
4566 integer :: tnew, told
4567
4568# ifdef DISTRIBUTE
4569 real(r8), parameter :: spv = 0.0_r8
4570# endif
4571 real(r8) :: eastsum, northsum, southsum, westsum
4572 real(r8) :: tl_eastsum, tl_northsum, tl_southsum, tl_westsum
4573# ifdef NESTING_DEBUG
4574 real(r8) :: mfratio
4575# endif
4576!
4577!-----------------------------------------------------------------------
4578! Check mass and volume conservation during refinement between coarse
4579! and fine grids.
4580!-----------------------------------------------------------------------
4581!
4582 DO cr=1,ncontact
4583!
4584! Get data donor and data receiver grid numbers.
4585!
4586 dg=rcontact(cr)%donor_grid
4587 rg=rcontact(cr)%receiver_grid
4588!
4589! Process only contact region data for requested nested finer grid
4590! "ngf". Notice that the donor grid is coarser than receiver grid.
4591!
4592 IF ((rg.eq.ngf).and.(dxmax(dg).gt.dxmax(rg))) THEN
4593!
4594! Set tile starting and ending indices for donor coarser grid.
4595!
4596 istr=bounds(dg)%Istr(tile)
4597 iend=bounds(dg)%Iend(tile)
4598 jstr=bounds(dg)%Jstr(tile)
4599 jend=bounds(dg)%Jend(tile)
4600!
4601! Set time rolling indices and conjugate region where the coarser
4602! donor grid becomes the receiver grid.
4603!
4604 told=3-rollingindex(cr)
4605 tnew=rollingindex(cr)
4606 DO icr=1,ncontact
4607 IF ((rg.eq.rcontact(icr)%donor_grid).and. &
4608 & (dg.eq.rcontact(icr)%receiver_grid)) THEN
4609 cjcr=icr
4610 EXIT
4611 END IF
4612 END DO
4613
4614# ifdef DISTRIBUTE
4615!
4616! Set global size of boundary edges for coarse grid (donor index).
4617!
4618 my_tile=-1
4619 ilb=bounds(dg)%LBi(my_tile)
4620 iub=bounds(dg)%UBi(my_tile)
4621 jlb=bounds(dg)%LBj(my_tile)
4622 jub=bounds(dg)%UBj(my_tile)
4623 nptswe=jub-jlb+1
4624 nptssn=iub-ilb+1
4625!
4626! If distributed-memory, initialize arrays used to check mass flux
4627! conservation with special value (zero) to facilitate the global
4628! reduction when collecting data between all nodes.
4629!
4630 bry_contact(iwest ,cjcr)%Mflux=spv
4631 bry_contact(ieast ,cjcr)%Mflux=spv
4632 bry_contact(isouth,cjcr)%Mflux=spv
4633 bry_contact(inorth,cjcr)%Mflux=spv
4634 bry_contact(iwest ,cjcr)%tl_Mflux=0.0_r8
4635 bry_contact(ieast ,cjcr)%tl_Mflux=0.0_r8
4636 bry_contact(isouth,cjcr)%tl_Mflux=0.0_r8
4637 bry_contact(inorth,cjcr)%tl_Mflux=0.0_r8
4638# endif
4639!
4640! Set finer grid center (half) and offset indices (Io and Jo) for
4641! coarser grid (I,J) coordinates.
4642!
4643 half=(refinescale(ngf)-1)/2
4644 io=half+1
4645 jo=half+1
4646!
4647!-----------------------------------------------------------------------
4648! Average finer grid western boundary mass fluxes and load them to the
4649! BRY_CONTACT structure.
4650!-----------------------------------------------------------------------
4651!
4652 ibc=i_left(ngf)
4653 jbc_min=j_bottom(ngf)
4654 jbc_max=j_top(ngf)-1 ! interior points, no top
4655! left corner
4656# ifdef NESTING_DEBUG
4657 IF (domain(ngf)%SouthWest_Test(tile)) THEN
4658 IF (master) THEN
4659 WRITE (301,10) 'Western Boundary Mass Fluxes: ', &
4660 & cr, dg, rg, iif(rg), iic(rg), int(time(rg))
4661 FLUSH (301)
4662 END IF
4663 END IF
4664!
4665# endif
4666 DO jbc=jstr,jend
4667 IF (((istr.le.ibc).and.(ibc.le.iend)).and. &
4668 & ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
4669!
4670! Sum finer grid western boundary mass fluxes within coarser grid cell.
4671!
4672 westsum=0.0_r8
4673 tl_westsum=0.0_r8
4674 jedge=jo+(jbc-jbc_min)*refinescale(ngf)
4675 DO jsum=-half,half
4676 jbf=jedge+jsum
4677 westsum=westsum+bry_contact(iwest,cr)%Mflux(jbf)
4678 tl_westsum=tl_westsum+ &
4679 & bry_contact(iwest,cr)%tl_Mflux(jbf)
4680 END DO
4681 m=bry_contact(iwest,cr)%C2Bindex(jbf) ! pick last one
4682!
4683! Load coarser grid western boundary mass flux that have been averaged
4684! from finer grid. These values can be compared with the coarser grid
4685! values REFINED(cr)%DU_avg2 to check if the mass flux between coarser
4686! and finer grid is conserved.
4687!
4688 bry_contact(iwest,cjcr)%Mflux(jbc)=westsum
4689 bry_contact(iwest,cjcr)%tl_Mflux(jbc)=tl_westsum
4690
4691# ifdef NESTING_DEBUG
4692 IF (westsum.ne.0) THEN
4693 mfratio=refined(cr)%DU_avg2(1,m,tnew)/westsum
4694 ELSE
4695 mfratio=1.0_r8
4696 END IF
4697 WRITE (301,30) jbc, refined(cr)%DU_avg2(1,m,tnew), &
4698 & westsum, mfratio
4699 FLUSH (301)
4700# endif
4701 END IF
4702 END DO
4703!
4704!-----------------------------------------------------------------------
4705! Average finer grid eastern boundary mass fluxes and load them to the
4706! BRY_CONTACT structure.
4707!-----------------------------------------------------------------------
4708!
4709 ibc=i_right(ngf)
4710 jbc_min=j_bottom(ngf)
4711 jbc_max=j_top(ngf)-1 ! interior points, no top
4712! right corner
4713# ifdef NESTING_DEBUG
4714 IF (domain(ngf)%SouthWest_Test(tile)) THEN
4715 IF (master) THEN
4716 WRITE (301,10) 'Eastern Boundary Mass Fluxes: ', &
4717 & cr, dg, rg, iif(rg), iic(rg), int(time(rg))
4718 FLUSH (301)
4719 END IF
4720 END IF
4721!
4722# endif
4723 DO jbc=jstr,jend
4724 IF (((istr.le.ibc).and.(ibc.le.iend)).and. &
4725 & ((jbc_min.le.jbc).and.(jbc.le.jbc_max))) THEN
4726!
4727! Sum finer grid eastern boundary mass fluxes within coarser grid cell.
4728!
4729 eastsum=0.0_r8
4730 tl_eastsum=0.0_r8
4731 jedge=jo+(jbc-jbc_min)*refinescale(ngf)
4732 DO jsum=-half,half
4733 jbf=jedge+jsum
4734 eastsum=eastsum+bry_contact(ieast,cr)%Mflux(jbf)
4735 tl_eastsum=tl_eastsum+ &
4736 & bry_contact(ieast,cr)%tl_Mflux(jbf)
4737 END DO
4738 m=bry_contact(ieast,cr)%C2Bindex(jbf) ! pick last one
4739!
4740! Load coarser grid eastern boundary mass flux that have been averaged
4741! from finer grid. These values can be compared with the coarser grid
4742! values REFINED(cr)%DU_avg2 to check if the mass flux between coarser
4743! and finer grid is conserved.
4744!
4745 bry_contact(ieast,cjcr)%Mflux(jbc)=eastsum
4746 bry_contact(ieast,cjcr)%tl_Mflux(jbc)=tl_eastsum
4747
4748# ifdef NESTING_DEBUG
4749 IF (eastsum.ne.0) THEN
4750 mfratio=refined(cr)%DU_avg2(1,m,tnew)/eastsum
4751 ELSE
4752 mfratio=1.0_r8
4753 END IF
4754 WRITE (301,30) jbc, refined(cr)%DU_avg2(1,m,tnew), &
4755 & eastsum, mfratio
4756 FLUSH (301)
4757# endif
4758 END IF
4759 END DO
4760!
4761!-----------------------------------------------------------------------
4762! Average finer grid southern boundary mass fluxes and load them to the
4763! BRY_CONTACT structure.
4764!-----------------------------------------------------------------------
4765!
4766 jbc=j_bottom(ngf)
4767 ibc_min=i_left(ngf)
4768 ibc_max=i_right(ngf)-1 ! interior points, no bottom
4769! right corner
4770# ifdef NESTING_DEBUG
4771 IF (domain(ngf)%SouthWest_Test(tile)) THEN
4772 IF (master) THEN
4773 WRITE (301,20) 'Southern Boundary Mass Fluxes: ', &
4774 & cr, dg, rg, iif(rg), iic(rg), int(time(rg))
4775 FLUSH (301)
4776 END IF
4777 END IF
4778!
4779# endif
4780 DO ibc=istr,iend
4781 IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and. &
4782 & ((jstr.le.jbc).and.(jbc.le.jend))) THEN
4783!
4784! Sum finer grid southern boundary mass fluxes within coarser grid
4785! cell.
4786!
4787 southsum=0.0_r8
4788 tl_southsum=0.0_r8
4789 iedge=io+(ibc-ibc_min)*refinescale(ngf)
4790 DO isum=-half,half
4791 ibf=iedge+isum
4792 southsum=southsum+bry_contact(isouth,cr)%Mflux(ibf)
4793 tl_southsum=tl_southsum+ &
4794 & bry_contact(isouth,cr)%tl_Mflux(ibf)
4795 END DO
4796 m=bry_contact(isouth,cr)%C2Bindex(ibf) ! pick last one
4797!
4798! Load coarser grid southern boundary mass flux that have been averaged
4799! from finer grid. These values can be compared with the coarser grid
4800! values REFINED(cr)%DU_avg2 to check if the mass flux between coarser
4801! and finer grid is conserved.
4802!
4803 bry_contact(isouth,cjcr)%Mflux(ibc)=southsum
4804 bry_contact(isouth,cjcr)%tl_Mflux(ibc)=tl_southsum
4805
4806# ifdef NESTING_DEBUG
4807 IF (southsum.ne.0) THEN
4808 mfratio=refined(cr)%DV_avg2(1,m,tnew)/southsum
4809 ELSE
4810 mfratio=1.0_r8
4811 END IF
4812 WRITE (301,30) ibc, refined(cr)%DV_avg2(1,m,tnew), &
4813 & southsum, mfratio
4814 FLUSH (301)
4815# endif
4816 END IF
4817 END DO
4818!
4819!-----------------------------------------------------------------------
4820! Average finer grid northern boundary mass fluxes and load them to the
4821! BRY_CONTACT structure.
4822!-----------------------------------------------------------------------
4823!
4824 jbc=j_top(ngf)
4825 ibc_min=i_left(ngf)
4826 ibc_max=i_right(ngf)-1 ! interior points, no top
4827! right corner
4828# ifdef NESTING_DEBUG
4829 IF (domain(ngf)%SouthWest_Test(tile)) THEN
4830 IF (master) THEN
4831 WRITE (301,20) 'Northern Boundary Mass Fluxes: ', &
4832 & cr, dg, rg, iif(rg), iic(rg), int(time(rg))
4833 FLUSH (301)
4834 END IF
4835 END IF
4836!
4837# endif
4838 DO ibc=istr,iend
4839 IF (((ibc_min.le.ibc).and.(ibc.le.ibc_max)).and. &
4840 & ((jstr.le.jbc).and.(jbc.le.jend))) THEN
4841!
4842! Sum finer grid northern boundary mass fluxes within coarser grid
4843! cell.
4844!
4845 northsum=0.0_r8
4846 tl_northsum=0.0_r8
4847 iedge=io+(ibc-ibc_min)*refinescale(ngf)
4848 DO isum=-half,half
4849 ibf=iedge+isum
4850 northsum=northsum+bry_contact(inorth,cr)%Mflux(ibf)
4851 tl_northsum=tl_northsum+ &
4852 & bry_contact(inorth,cr)%tl_Mflux(ibf)
4853 END DO
4854 m=bry_contact(inorth,cr)%C2Bindex(ibf) ! pick last one
4855!
4856! Load coarser grid northern boundary mass flux that have been averaged
4857! from finer grid. These values can be compared with the coarser grid
4858! values REFINED(cr)%DU_avg2 to check if the mass flux between coarser
4859! and finer grid is conserved.
4860!
4861 bry_contact(inorth,cjcr)%Mflux(ibc)=northsum
4862 bry_contact(inorth,cjcr)%tl_Mflux(ibc)=tl_northsum
4863
4864# ifdef NESTING_DEBUG
4865 IF (northsum.ne.0) THEN
4866 mfratio=refined(cr)%DV_avg2(1,m,tnew)/northsum
4867 ELSE
4868 mfratio=1.0_r8
4869 END IF
4870 WRITE (301,30) ibc, refined(cr)%DV_avg2(1,m,tnew), &
4871 & northsum, mfratio
4872# endif
4873 END IF
4874 END DO
4875
4876# ifdef DISTRIBUTE
4877!
4878! Collect data from all nodes.
4879!
4880 CALL mp_assemble (dg, model, nptswe, spv, &
4881 & bry_contact(iwest ,cjcr)%Mflux(jlb:))
4882 CALL mp_assemble (dg, model, nptswe, spv, &
4883 & bry_contact(ieast ,cjcr)%Mflux(jlb:))
4884 CALL mp_assemble (dg, model, nptssn, spv, &
4885 & bry_contact(isouth,cjcr)%Mflux(ilb:))
4886 CALL mp_assemble (dg, model, nptssn, spv, &
4887 & bry_contact(inorth,cjcr)%Mflux(ilb:))
4888 CALL mp_assemble (dg, model, nptswe, spv, &
4889 & bry_contact(iwest ,cjcr)%tl_Mflux(jlb:))
4890 CALL mp_assemble (dg, model, nptswe, spv, &
4891 & bry_contact(ieast ,cjcr)%tl_Mflux(jlb:))
4892 CALL mp_assemble (dg, model, nptssn, spv, &
4893 & bry_contact(isouth,cjcr)%tl_Mflux(ilb:))
4894 CALL mp_assemble (dg, model, nptssn, spv, &
4895 & bry_contact(inorth,cjcr)%tl_Mflux(ilb:))
4896# endif
4897 END IF
4898 END DO
4899
4900# ifdef NESTING_DEBUG
4901!
4902 FLUSH (301)
4903!
4904 10 FORMAT (/,1x,a,/,/,4x,'cr = ',i2.2,4x,'dg = ',i2.2,4x,'rg = ', &
4905 & i2.2,4x,'iif(rg) = ',i3.3,4x,'iic(rg) = ',i10.10,4x, &
4906 & 'time(rg) = ',i10.10,/,/,2x,'Coarse',6x,'Coarse Grid',8x, &
4907 & 'Fine Grid',11x,'Ratio',/,4x,'Jb',9x,'DU_avg2',9x, &
4908 & 'SUM(DU_avg2)',/)
4909 20 FORMAT (/,1x,a,/,/,4x,'cr = ',i2.2,4x,'dg = ',i2.2,4x,'rg = ', &
4910 & i2.2,4x,'iif(rg) = ',i3.3,4x,'iic(rg) = ',i10.10,4x, &
4911 & 'time(rg) = ',i10.10,/,/,2x,'Coarse',6x,'Coarse Grid',8x, &
4912 & 'Fine Grid',11x,'Ratio',/,4x,'Ib',9x,'DV_avg2',9x, &
4913 & 'SUM(DV_avg2)',/)
4914 30 FORMAT (4x,i4.4,3(3x,1p,e15.8))
4915# endif
4916!
4917 RETURN
4918 END SUBROUTINE tl_check_massflux
4919#endif
4920 END MODULE tl_nesting_mod
distribute_mod::mp_assemble
Definition distribute.F:67
distribute_mod
Definition distribute.F:3
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_u3d_tile
subroutine exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:495
exchange_3d_mod::exchange_r3d_tile
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:283
exchange_3d_mod::exchange_v3d_tile
subroutine exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:707
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_iounits
Definition mod_iounits.F:2
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::idubar
integer idubar
Definition mod_ncparam.F:325
mod_ncparam::idvvel
integer idvvel
Definition mod_ncparam.F:381
mod_ncparam::idfsur
integer idfsur
Definition mod_ncparam.F:206
mod_ncparam::idtvar
integer, dimension(:), allocatable idtvar
Definition mod_ncparam.F:568
mod_ncparam::idvfx1
integer idvfx1
Definition mod_ncparam.F:372
mod_ncparam::idvbms
integer idvbms
Definition mod_ncparam.F:365
mod_ncparam::iduvel
integer iduvel
Definition mod_ncparam.F:343
mod_ncparam::vname
character(len=maxlen), dimension(6, 0:nv) vname
Definition mod_ncparam.F:591
mod_ncparam::idrzet
integer idrzet
Definition mod_ncparam.F:286
mod_ncparam::idufx1
integer idufx1
Definition mod_ncparam.F:337
mod_ncparam::idubms
integer idubms
Definition mod_ncparam.F:330
mod_ncparam::idvbar
integer idvbar
Definition mod_ncparam.F:360
mod_nesting
Definition mod_nesting.F:2
mod_nesting::nrhst
integer, parameter nrhst
Definition mod_nesting.F:87
mod_nesting::ngetd
integer, parameter ngetd
Definition mod_nesting.F:77
mod_nesting::n2dfx
integer, parameter n2dfx
Definition mod_nesting.F:92
mod_nesting::nputd
integer, parameter nputd
Definition mod_nesting.F:79
mod_nesting::nzeta
integer, parameter nzeta
Definition mod_nesting.F:88
mod_nesting::vcontact
type(t_ngc), dimension(:), allocatable vcontact
Definition mod_nesting.F:288
mod_nesting::bry_contact
type(t_bcp), dimension(:,:), allocatable bry_contact
Definition mod_nesting.F:335
mod_nesting::ntvic
integer, parameter ntvic
Definition mod_nesting.F:85
mod_nesting::rollingindex
integer, dimension(:), allocatable rollingindex
Definition mod_nesting.F:178
mod_nesting::nzwgt
integer, parameter nzwgt
Definition mod_nesting.F:89
mod_nesting::ndxdy
integer, parameter ndxdy
Definition mod_nesting.F:76
mod_nesting::i_right
integer, dimension(:), allocatable i_right
Definition mod_nesting.F:194
mod_nesting::nfsic
integer, parameter nfsic
Definition mod_nesting.F:82
mod_nesting::i_left
integer, dimension(:), allocatable i_left
Definition mod_nesting.F:193
mod_nesting::n3dic
integer, parameter n3dic
Definition mod_nesting.F:84
mod_nesting::nmask
integer, parameter nmask
Definition mod_nesting.F:78
mod_nesting::n2way
integer, parameter n2way
Definition mod_nesting.F:80
mod_nesting::receiver_grid
integer, dimension(:), allocatable receiver_grid
Definition mod_nesting.F:173
mod_nesting::coarserdonor
integer, dimension(:), allocatable coarserdonor
Definition mod_nesting.F:138
mod_nesting::n3duv
integer, parameter n3duv
Definition mod_nesting.F:93
mod_nesting::donor_grid
integer, dimension(:), allocatable donor_grid
Definition mod_nesting.F:172
mod_nesting::j_bottom
integer, dimension(:), allocatable j_bottom
Definition mod_nesting.F:195
mod_nesting::n2dps
integer, parameter n2dps
Definition mod_nesting.F:90
mod_nesting::refined
type(t_refined), dimension(:), allocatable refined
Definition mod_nesting.F:437
mod_nesting::j_top
integer, dimension(:), allocatable j_top
Definition mod_nesting.F:196
mod_nesting::composite
type(t_composite), dimension(:), allocatable composite
Definition mod_nesting.F:409
mod_nesting::n2dic
integer, parameter n2dic
Definition mod_nesting.F:83
mod_nesting::n3dtv
integer, parameter n3dtv
Definition mod_nesting.F:94
mod_nesting::rcontact
type(t_ngc), dimension(:), allocatable rcontact
Definition mod_nesting.F:286
mod_nesting::nbstr
integer, parameter nbstr
Definition mod_nesting.F:86
mod_nesting::rollingtime
real(dp), dimension(:,:), allocatable rollingtime
Definition mod_nesting.F:179
mod_nesting::get_vweights
logical get_vweights
Definition mod_nesting.F:130
mod_nesting::ucontact
type(t_ngc), dimension(:), allocatable ucontact
Definition mod_nesting.F:287
mod_nesting::nmflx
integer, parameter nmflx
Definition mod_nesting.F:75
mod_nesting::ncontact
integer ncontact
Definition mod_nesting.F:256
mod_nesting::n2dcs
integer, parameter n2dcs
Definition mod_nesting.F:91
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_parallel::first_tile
integer, dimension(:), allocatable first_tile
Definition mod_parallel.F:62
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_parallel::last_tile
integer, dimension(:), allocatable last_tile
Definition mod_parallel.F:63
mod_param
Definition mod_param.F:2
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::bounds
type(t_bounds), dimension(:), allocatable bounds
Definition mod_param.F:232
mod_param::r3dvar
integer, parameter r3dvar
Definition mod_param.F:721
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::u3dvar
integer, parameter u3dvar
Definition mod_param.F:722
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::u2dvar
integer, parameter u2dvar
Definition mod_param.F:718
mod_param::w3dvar
integer, parameter w3dvar
Definition mod_param.F:724
mod_param::nt
integer, dimension(:), allocatable nt
Definition mod_param.F:489
mod_param::r2dvar
integer, parameter r2dvar
Definition mod_param.F:717
mod_param::v2dvar
integer, parameter v2dvar
Definition mod_param.F:719
mod_param::v3dvar
integer, parameter v3dvar
Definition mod_param.F:723
mod_scalars
Definition mod_scalars.F:2
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::dcrit
real(r8), dimension(:), allocatable dcrit
Definition mod_scalars.F:782
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::exit_flag
integer exit_flag
Definition mod_scalars.F:559
mod_scalars::indx1
integer, dimension(:), allocatable indx1
Definition mod_scalars.F:248
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::dxmax
real(dp), dimension(:), allocatable dxmax
Definition mod_scalars.F:675
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::time
real(dp), dimension(:), allocatable time
Definition mod_scalars.F:262
mod_scalars::refinedgrid
logical, dimension(:), allocatable refinedgrid
Definition mod_scalars.F:494
mod_scalars::ltracerclm
logical, dimension(:,:), allocatable ltracerclm
Definition mod_scalars.F:528
mod_scalars::refinescale
integer, dimension(:), allocatable refinescale
Definition mod_scalars.F:498
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::lnudgetclm
logical, dimension(:,:), allocatable lnudgetclm
Definition mod_scalars.F:535
mod_scalars::iif
integer, dimension(:), allocatable iif
Definition mod_scalars.F:250
mod_scalars::noerror
integer noerror
Definition mod_scalars.F:561
mod_stepping
Definition mod_stepping.F:2
mod_stepping::knew
integer, dimension(:), allocatable knew
Definition mod_stepping.F:64
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange4d
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
Definition mp_exchange.F:3377
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
nesting_mod
Definition nesting.F:2
nesting_mod::get_persisted2d
subroutine, public get_persisted2d(dg, rg, model, tile, gtype, svname, cr, npoints, contact, lbi, ubi, lbj, ubj, ad, ac)
Definition nesting.F:5150
nesting_mod::get_contact3d
subroutine, public get_contact3d(dg, model, tile, gtype, svname, cr, npoints, contact, lbi, ubi, lbj, ubj, lbk, ubk, ad, ac)
Definition nesting.F:4988
nesting_mod::get_metrics
subroutine, public get_metrics(ng, model, tile)
Definition nesting.F:2245
nesting_mod::put_contact2d
subroutine put_contact2d(rg, model, tile, gtype, svname, cr, npoints, contact, lbi, ubi, lbj, ubj, amask, ac, ar)
Definition nesting.F:5356
nesting_mod::get_contact2d
subroutine, public get_contact2d(dg, model, tile, gtype, svname, cr, npoints, contact, lbi, ubi, lbj, ubj, ad, ac)
Definition nesting.F:4829
nesting_mod::mask_hweights
subroutine, public mask_hweights(ng, model, tile)
Definition nesting.F:1200
nesting_mod::fine2coarse2d
subroutine, public fine2coarse2d(ng, dg, model, tile, gtype, svname, areaavg, rscale, cr, npoints, contact, lbif, ubif, lbjf, ubjf, lbic, ubic, lbjc, ubjc, adx, ady, pmc, pnc, hhc, amsk, cmsk, a, c1, c2)
Definition nesting.F:3879
nesting_mod::fine2coarse3d
subroutine, public fine2coarse3d(ng, dg, model, tile, gtype, svname, areaavg, rscale, cr, npoints, contact, lbif, ubif, lbjf, ubjf, lbkf, ubkf, lbic, ubic, lbjc, ubjc, lbkc, ubkc, adx, ady, pmc, pnc, amsk, cmsk, a, c)
Definition nesting.F:4365
set_depth_mod
Definition set_depth.F:2
set_depth_mod::set_depth
subroutine, public set_depth(ng, tile, model)
Definition set_depth.F:34
strings_mod
Definition strings.F:2
strings_mod::founderror
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
tl_nesting_mod
Definition tl_nesting.F:2
tl_nesting_mod::tl_get_composite
subroutine, private tl_get_composite(ng, model, isection, tile)
Definition tl_nesting.F:366
tl_nesting_mod::tl_put_refine2d
subroutine, private tl_put_refine2d(ng, dg, cr, model, tile, lputfsur, lbi, ubi, lbj, ubj)
Definition tl_nesting.F:2487
tl_nesting_mod::tl_put_refine3d
subroutine, private tl_put_refine3d(ng, dg, cr, model, tile, lbi, ubi, lbj, ubj)
Definition tl_nesting.F:3209
tl_nesting_mod::tl_put_contact3d
subroutine tl_put_contact3d(rg, model, tile, gtype, svname, cr, npoints, contact, lbi, ubi, lbj, ubj, lbk, ubk, amask, ac, tl_ac, tl_ar)
Definition tl_nesting.F:4348
tl_nesting_mod::tl_check_massflux
subroutine, public tl_check_massflux(ngf, model, tile)
Definition tl_nesting.F:4524
tl_nesting_mod::tl_nesting
subroutine, public tl_nesting(ng, model, isection)
Definition tl_nesting.F:114
tl_nesting_mod::tl_put_composite
subroutine, private tl_put_composite(ng, model, isection, tile)
Definition tl_nesting.F:1013
tl_nesting_mod::tl_correct_tracer
subroutine, private tl_correct_tracer(ng, ngf, model, tile)
Definition tl_nesting.F:1497
tl_nesting_mod::tl_get_refine
subroutine, private tl_get_refine(ng, model, tile)
Definition tl_nesting.F:725
tl_nesting_mod::tl_put_refine
subroutine, private tl_put_refine(ng, model, tile, lputfsur)
Definition tl_nesting.F:1402
tl_nesting_mod::tl_fine2coarse
subroutine, private tl_fine2coarse(ng, model, vtype, tile)
Definition tl_nesting.F:2017
tl_nesting_mod::tl_z_weights
subroutine, private tl_z_weights(ng, model, tile)
Definition tl_nesting.F:3581
tl_nesting_mod::tl_correct_tracer_tile
subroutine, private tl_correct_tracer_tile(ngc, ngf, model, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs)
Definition tl_nesting.F:1560
tl_set_depth_mod
Definition tl_set_depth.F:2
tl_set_depth_mod::tl_set_depth
subroutine, public tl_set_depth(ng, tile, model)
Definition tl_set_depth.F:41
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3