ROMS
Loading...
Searching...
No Matches
metrics.F
Go to the documentation of this file.
1#include "cppdefs.h"
3!
4!git $Id$
5!=======================================================================
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md Hernan G. Arango !
9!========================================== Alexander F. Shchepetkin ===
10! !
11! This routine computes various horizontal metric terms. !
12! !
13!=======================================================================
14!
15 implicit none
16
17 PRIVATE
18 PUBLIC :: metrics
19
20 CONTAINS
21!
22!***********************************************************************
23 SUBROUTINE metrics (ng, tile, model)
24!***********************************************************************
25!
26 USE mod_param
27 USE mod_grid
28#if defined DIFF_3DCOEF || defined VISC_3DCOEF
29 USE mod_mixing
30#endif
31 USE mod_ocean
32 USE mod_stepping
33!
34! Imported variable declarations.
35!
36 integer, intent(in) :: ng, tile, model
37!
38! Local variable declarations.
39!
40#include "tile.h"
41!
42 CALL metrics_tile (ng, tile, model, &
43 & lbi, ubi, lbj, ubj, &
44 & imins, imaxs, jmins, jmaxs, &
45 & nstp(ng), nnew(ng), &
46 & grid(ng) % f, &
47 & grid(ng) % h, &
48 & grid(ng) % pm, &
49 & grid(ng) % pn, &
50#ifdef MASKING
51 & grid(ng) % pmask, &
52 & grid(ng) % rmask, &
53#endif
54 & grid(ng) % angler, &
55 & grid(ng) % CosAngler, &
56 & grid(ng) % SinAngler, &
57#if defined TIDE_GENERATING_FORCES
58 & grid(ng) % Cos2Lat, &
59 & grid(ng) % SinLat2, &
60 & grid(ng) % latr, &
61#endif
62#ifdef SOLVE3D
63# ifdef ICESHELF
64 & grid(ng) % zice, &
65# endif
66 & grid(ng) % Hz, &
67 & grid(ng) % z_r, &
68 & grid(ng) % z_w, &
69#endif
70#if defined DIFF_GRID || defined DIFF_3DCOEF || \
71 defined visc_grid || defined visc_3dcoef
72 & grid(ng) % grdscl, &
73#endif
74#if defined DIFF_3DCOEF
75 & mixing(ng) % Hdiffusion, &
76#endif
77#if defined VISC_3DCOEF
78 & mixing(ng) % Hviscosity, &
79#endif
80 & grid(ng) % om_p, &
81 & grid(ng) % om_r, &
82 & grid(ng) % om_u, &
83 & grid(ng) % om_v, &
84 & grid(ng) % on_p, &
85 & grid(ng) % on_r, &
86 & grid(ng) % on_u, &
87 & grid(ng) % on_v, &
88 & grid(ng) % fomn, &
89 & grid(ng) % omn, &
90 & grid(ng) % pnom_p, &
91 & grid(ng) % pnom_r, &
92 & grid(ng) % pnom_u, &
93 & grid(ng) % pnom_v, &
94 & grid(ng) % pmon_p, &
95 & grid(ng) % pmon_r, &
96 & grid(ng) % pmon_u, &
97 & grid(ng) % pmon_v)
98 RETURN
99 END SUBROUTINE metrics
100!
101!***********************************************************************
102 SUBROUTINE metrics_tile (ng, tile, model, &
103 & LBi, UBi, LBj, UBj, &
104 & IminS, ImaxS, JminS, JmaxS, &
105 & nstp, nnew, &
106 & f, h, pm, pn, &
107#ifdef MASKING
108 & pmask, rmask, &
109#endif
110 & angler, CosAngler, SinAngler, &
111#if defined TIDE_GENERATING_FORCES
112 & Cos2Lat, SinLat2, latr, &
113#endif
114#ifdef SOLVE3D
115# ifdef ICESHELF
116 & zice, &
117# endif
118 & Hz, z_r, z_w, &
119#endif
120#if defined DIFF_GRID || defined DIFF_3DCOEF || \
121 defined VISC_GRID || defined VISC_3DCOEF
122 & grdscl, &
123#endif
124#if defined DIFF_3DCOEF
125 & Hdiffusion, &
126#endif
127#if defined VISC_3DCOEF
128 & Hviscosity, &
129#endif
130 & om_p, om_r, om_u, om_v, &
131 & on_p, on_r, on_u, on_v, &
132 & fomn, omn, &
133 & pnom_p, pnom_r, pnom_u, pnom_v, &
134 & pmon_p, pmon_r, pmon_u, pmon_v)
135!***********************************************************************
136!
137 USE mod_param
138 USE mod_parallel
139#ifdef FOUR_DVAR
140 USE mod_fourdvar
141#endif
142 USE mod_iounits
143 USE mod_ncparam
144#ifdef NESTING
145 USE mod_nesting
146#endif
147 USE mod_scalars
148 USE mod_iounits
149!
150#ifdef DISTRIBUTE
151 USE distribute_mod, ONLY : mp_reduce
152#endif
154#ifdef DISTRIBUTE
156#endif
157#ifdef SOLVE3D
158 USE set_depth_mod, ONLY : set_depth_tile
159#endif
160 implicit none
161!
162! Imported variable declarations.
163!
164 integer, intent(in) :: ng, tile, model
165 integer, intent(in) :: LBi, UBi, LBj, UBj
166 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
167 integer, intent(in) :: nstp, nnew
168
169#ifdef ASSUMED_SHAPE
170 real(r8), intent(in) :: f(LBi:,LBj:)
171 real(r8), intent(in) :: pm(LBi:,LBj:)
172 real(r8), intent(in) :: pn(LBi:,LBj:)
173# if defined TIDE_GENERATING_FORCES
174 real(r8), intent(in) :: latr(LBi:,LBj:)
175# endif
176 real(r8), intent(inout) :: h(LBi:,LBj:)
177# ifdef MASKING
178 real(r8), intent(inout) :: pmask(LBi:,LBj:)
179 real(r8), intent(inout) :: rmask(LBi:,LBj:)
180# endif
181 real(r8), intent(inout) :: angler(LBi:,LBj:)
182# ifdef SOLVE3D
183# ifdef ICESHELF
184 real(r8), intent(in) :: zice(LBi:,LBj:)
185# endif
186# endif
187#if defined DIFF_GRID || defined DIFF_3DCOEF || \
188 defined visc_grid || defined visc_3dcoef
189 real(r8), intent(out) :: grdscl(LBi:,LBj:)
190# endif
191# if defined DIFF_3DCOEF
192 real(r8), intent(out) :: Hdiffusion(LBi:,LBj:)
193# endif
194# if defined VISC_3DCOEF
195 real(r8), intent(out) :: Hviscosity(LBi:,LBj:)
196# endif
197 real(r8), intent(out) :: om_p(LBi:,LBj:)
198 real(r8), intent(out) :: om_r(LBi:,LBj:)
199 real(r8), intent(out) :: om_u(LBi:,LBj:)
200 real(r8), intent(out) :: om_v(LBi:,LBj:)
201 real(r8), intent(out) :: on_p(LBi:,LBj:)
202 real(r8), intent(out) :: on_r(LBi:,LBj:)
203 real(r8), intent(out) :: on_u(LBi:,LBj:)
204 real(r8), intent(out) :: on_v(LBi:,LBj:)
205 real(r8), intent(out) :: fomn(LBi:,LBj:)
206 real(r8), intent(out) :: omn(LBi:,LBj:)
207 real(r8), intent(out) :: pnom_p(LBi:,LBj:)
208 real(r8), intent(out) :: pnom_r(LBi:,LBj:)
209 real(r8), intent(out) :: pnom_u(LBi:,LBj:)
210 real(r8), intent(out) :: pnom_v(LBi:,LBj:)
211 real(r8), intent(out) :: pmon_p(LBi:,LBj:)
212 real(r8), intent(out) :: pmon_r(LBi:,LBj:)
213 real(r8), intent(out) :: pmon_u(LBi:,LBj:)
214 real(r8), intent(out) :: pmon_v(LBi:,LBj:)
215 real(r8), intent(out) :: CosAngler(LBi:,LBj:)
216 real(r8), intent(out) :: SinAngler(LBi:,LBj:)
217# if defined TIDE_GENERATING_FORCES
218 real(r8), intent(out) :: Cos2Lat(LBi:,LBj:)
219 real(r8), intent(out) :: SinLat2(LBi:,LBj:)
220# endif
221# ifdef SOLVE3D
222 real(r8), intent(out) :: Hz(LBi:,LBj:,:)
223 real(r8), intent(out) :: z_r(LBi:,LBj:,:)
224 real(r8), intent(out) :: z_w(LBi:,LBj:,0:)
225# endif
226#else
227 real(r8), intent(in) :: f(LBi:UBi,LBj:UBj)
228 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
229 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
230# if defined TIDE_GENERATING_FORCES
231 real(r8), intent(in) :: latr(LBi:UBi,LBj:UBj)
232# endif
233 real(r8), intent(inout) :: h(LBi:UBi,LBj:UBj)
234# ifdef MASKING
235 real(r8), intent(inout) :: pmask(LBi:UBi,LBj:UBj)
236 real(r8), intent(inout) :: rmask(LBi:UBi,LBj:UBj)
237# endif
238 real(r8), intent(inout) :: angler(LBi:UBi,LBj:UBj)
239# ifdef SOLVE3D
240# ifdef ICESHELF
241 real(r8), intent(in) :: zice(LBi:UBi,LBj:UBj)
242# endif
243# endif
244#if defined DIFF_GRID || defined DIFF_3DCOEF || \
245 defined visc_grid || defined visc_3dcoef
246 real(r8), intent(out) :: grdscl(LBi:UBi,LBj:UBj)
247# endif
248# if defined DIFF_3DCOEF
249 real(r8), intent(out) :: Hdiffusion(LBi:UBi,LBj:UBj)
250# endif
251# if defined VISC_3DCOEF
252 real(r8), intent(out) :: Hviscosity(LBi:UBi,LBj:UBj)
253# endif
254 real(r8), intent(out) :: om_p(LBi:UBi,LBj:UBj)
255 real(r8), intent(out) :: om_r(LBi:UBi,LBj:UBj)
256 real(r8), intent(out) :: om_u(LBi:UBi,LBj:UBj)
257 real(r8), intent(out) :: om_v(LBi:UBi,LBj:UBj)
258 real(r8), intent(out) :: on_p(LBi:UBi,LBj:UBj)
259 real(r8), intent(out) :: on_r(LBi:UBi,LBj:UBj)
260 real(r8), intent(out) :: on_u(LBi:UBi,LBj:UBj)
261 real(r8), intent(out) :: on_v(LBi:UBi,LBj:UBj)
262 real(r8), intent(out) :: fomn(LBi:UBi,LBj:UBj)
263 real(r8), intent(out) :: omn(LBi:UBi,LBj:UBj)
264 real(r8), intent(out) :: pnom_p(LBi:UBi,LBj:UBj)
265 real(r8), intent(out) :: pnom_r(LBi:UBi,LBj:UBj)
266 real(r8), intent(out) :: pnom_u(LBi:UBi,LBj:UBj)
267 real(r8), intent(out) :: pnom_v(LBi:UBi,LBj:UBj)
268 real(r8), intent(out) :: pmon_p(LBi:UBi,LBj:UBj)
269 real(r8), intent(out) :: pmon_r(LBi:UBi,LBj:UBj)
270 real(r8), intent(out) :: pmon_u(LBi:UBi,LBj:UBj)
271 real(r8), intent(out) :: pmon_v(LBi:UBi,LBj:UBj)
272 real(r8), intent(out) :: CosAngler(LBi:UBi,LBj:UBj)
273 real(r8), intent(out) :: SinAngler(LBi:UBi,LBj:UBj)
274# if defined TIDE_GENERATING_FORCES
275 real(r8), intent(out) :: Cos2Lat(LBi:UBi,LBj:UBj)
276 real(r8), intent(out) :: SinLat2(LBi:UBi,LBj:UBj)
277# endif
278# ifdef SOLVE3D
279 real(r8), intent(out) :: Hz(LBi:UBi,LBj:UBj,N(ng))
280 real(r8), intent(out) :: z_r(LBi:UBi,LBj:UBj,N(ng))
281 real(r8), intent(out) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
282# endif
283#endif
284!
285! Local variable declarations.
286!
287 integer :: NSUB, i, ibry, is, j, k, rec
288#ifdef NESTING
289 integer :: cr, dg, ig, rg
290#endif
291#if defined DIFF_3DCOEF || defined VISC_3DCOEF
292 real(r8), parameter :: PecletCoef = 1.0_r8 / 12.0_r8
293 real(r8), parameter :: Uscale = 0.1_r8
294#endif
295 real(r8) :: cff, cff1, cff2
296#if defined TIDE_GENERATING_FORCES
297 real(r8) :: cosphi, phi
298#endif
299#ifdef NESTING
300 real(r8) :: SecScale
301#endif
302 real(dp) :: my_DXmax, my_DXmin, my_DYmax, my_DYmin
303#ifdef MASKING
304 real(dp) :: my_DXmaxW, my_DXminW, my_DYmaxW, my_DYminW
305#endif
306
307#ifdef SOLVE3D
308 real(dp) :: my_DZmax, my_DZmin
309# ifdef MASKING
310 real(dp) :: my_DZmaxW, my_DZminW
311# endif
312#endif
313 real(dp) :: my_Cg_Cor, my_Cg_max, my_Cg_min, my_grdmax
314#if defined DIFF_3DCOEF
315 real(dp) :: my_DiffMax, my_DiffMin
316#endif
317#if defined VISC_3DCOEF
318 real(dp) :: my_ViscMax, my_ViscMin
319#endif
320#if defined DIFF_3DCOEF || defined VISC_3DCOEF
321 character (len=4) :: units
322#endif
323#if defined FOUR_DVAR
324 character (len=50) :: Text
325#endif
326
327#ifdef DISTRIBUTE
328# ifdef MASKING
329 real(dp), dimension(20) :: rbuffer
330 character (len=3), dimension(20) :: op_handle
331# else
332 real(dp), dimension(14) :: rbuffer
333 character (len=3), dimension(14) :: op_handle
334# endif
335#endif
336
337#ifdef SOLVE3D
338 real(r8), dimension(LBi:UBi,LBj:UBj) :: A2d
339#endif
340
341#include "set_bounds.h"
342!
343!-----------------------------------------------------------------------
344! Compute 1/m, 1/n, 1/mn, and f/mn at horizontal RHO-points.
345!-----------------------------------------------------------------------
346!
347 DO j=jstrt,jendt
348 DO i=istrt,iendt
349 om_r(i,j)=1.0_r8/pm(i,j)
350 on_r(i,j)=1.0_r8/pn(i,j)
351 omn(i,j)=1.0_r8/(pm(i,j)*pn(i,j))
352 fomn(i,j)=f(i,j)*omn(i,j)
353 END DO
354 END DO
355!
356! Exchange boundary information.
357!
358 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
359 CALL exchange_r2d_tile (ng, tile, &
360 & lbi, ubi, lbj, ubj, &
361 & om_r)
362 CALL exchange_r2d_tile (ng, tile, &
363 & lbi, ubi, lbj, ubj, &
364 & on_r)
365 CALL exchange_r2d_tile (ng, tile, &
366 & lbi, ubi, lbj, ubj, &
367 & omn)
368 CALL exchange_r2d_tile (ng, tile, &
369 & lbi, ubi, lbj, ubj, &
370 & fomn)
371 END IF
372
373#ifdef DISTRIBUTE
374 CALL mp_exchange2d (ng, tile, model, 4, &
375 & lbi, ubi, lbj, ubj, &
376 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
377 & om_r, on_r, omn, fomn)
378#endif
379!
380!-----------------------------------------------------------------------
381! Compute n/m, and m/n at horizontal RHO-points.
382!-----------------------------------------------------------------------
383!
384 DO j=jstrt,jendt
385 DO i=istrt,iendt
386 pnom_r(i,j)=pn(i,j)/pm(i,j)
387 pmon_r(i,j)=pm(i,j)/pn(i,j)
388 END DO
389 END DO
390!
391! Exchange boundary information.
392!
393 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
394 CALL exchange_r2d_tile (ng, tile, &
395 & lbi, ubi, lbj, ubj, &
396 & pnom_r)
397 CALL exchange_r2d_tile (ng, tile, &
398 & lbi, ubi, lbj, ubj, &
399 & pmon_r)
400 END IF
401
402#ifdef DISTRIBUTE
403 CALL mp_exchange2d (ng, tile, model, 2, &
404 & lbi, ubi, lbj, ubj, &
405 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
406 & pnom_r, pmon_r)
407#endif
408!
409!-----------------------------------------------------------------------
410! Compute m/n, 1/m, and 1/n at horizontal U-points.
411!-----------------------------------------------------------------------
412!
413 DO j=jstrt,jendt
414 DO i=istrp,iendt
415 pmon_u(i,j)=(pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
416 pnom_u(i,j)=(pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
417 om_u(i,j)=2.0_r8/(pm(i-1,j)+pm(i,j))
418 on_u(i,j)=2.0_r8/(pn(i-1,j)+pn(i,j))
419 END DO
420 END DO
421!
422! Exchange boundary information.
423!
424 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
425 CALL exchange_u2d_tile (ng, tile, &
426 & lbi, ubi, lbj, ubj, &
427 & pmon_u)
428 CALL exchange_u2d_tile (ng, tile, &
429 & lbi, ubi, lbj, ubj, &
430 & pnom_u)
431 CALL exchange_u2d_tile (ng, tile, &
432 & lbi, ubi, lbj, ubj, &
433 & om_u)
434 CALL exchange_u2d_tile (ng, tile, &
435 & lbi, ubi, lbj, ubj, &
436 & on_u)
437 END IF
438
439#ifdef DISTRIBUTE
440 CALL mp_exchange2d (ng, tile, model, 4, &
441 & lbi, ubi, lbj, ubj, &
442 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
443 & pmon_u, pnom_u, om_u, on_u)
444#endif
445!
446!-----------------------------------------------------------------------
447! Compute n/m, 1/m, and 1/m at horizontal V-points.
448!-----------------------------------------------------------------------
449!
450 DO j=jstrp,jendt
451 DO i=istrt,iendt
452 pmon_v(i,j)=(pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j))
453 pnom_v(i,j)=(pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j))
454 om_v(i,j)=2.0_r8/(pm(i,j-1)+pm(i,j))
455 on_v(i,j)=2.0_r8/(pn(i,j-1)+pn(i,j))
456 END DO
457 END DO
458!
459! Exchange boundary information.
460!
461 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
462 CALL exchange_v2d_tile (ng, tile, &
463 & lbi, ubi, lbj, ubj, &
464 & pmon_v)
465 CALL exchange_v2d_tile (ng, tile, &
466 & lbi, ubi, lbj, ubj, &
467 & pnom_v)
468 CALL exchange_v2d_tile (ng, tile, &
469 & lbi, ubi, lbj, ubj, &
470 & om_v)
471 CALL exchange_v2d_tile (ng, tile, &
472 & lbi, ubi, lbj, ubj, &
473 & on_v)
474 END IF
475
476#ifdef DISTRIBUTE
477 CALL mp_exchange2d (ng, tile, model, 4, &
478 & lbi, ubi, lbj, ubj, &
479 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
480 & pmon_v, pnom_v, om_v, on_v)
481#endif
482!
483!-----------------------------------------------------------------------
484! Compute n/m and m/n at horizontal PSI-points.
485!-----------------------------------------------------------------------
486!
487 DO j=jstrp,jendt
488 DO i=istrp,iendt
489 pnom_p(i,j)=(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))/ &
490 & (pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
491 pmon_p(i,j)=(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))/ &
492 & (pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
493 om_p(i,j)=4.0_r8/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
494 on_p(i,j)=4.0_r8/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
495 END DO
496 END DO
497!
498! Exchange boundary information.
499!
500 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
501 CALL exchange_p2d_tile (ng, tile, &
502 & lbi, ubi, lbj, ubj, &
503 & pnom_p)
504 CALL exchange_p2d_tile (ng, tile, &
505 & lbi, ubi, lbj, ubj, &
506 & pmon_p)
507 CALL exchange_p2d_tile (ng, tile, &
508 & lbi, ubi, lbj, ubj, &
509 & om_p)
510 CALL exchange_p2d_tile (ng, tile, &
511 & lbi, ubi, lbj, ubj, &
512 & on_p)
513 END IF
514
515#ifdef DISTRIBUTE
516 CALL mp_exchange2d (ng, tile, model, 4, &
517 & lbi, ubi, lbj, ubj, &
518 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
519 & pnom_p, pmon_p, om_p, on_p)
520#endif
521
522#ifdef MASKING
523!
524!-----------------------------------------------------------------------
525! Set slipperiness (no-slip) mask at PSI-points.
526!-----------------------------------------------------------------------
527!
528! Set no-slip boundary conditions on land-mask boundaries regardless of
529! supplied value of gamma2.
530!
531 cff1=1.0_r8 ! computation of off-diagonal nonlinear terms
532 cff2=2.0_r8
533 DO j=jstrp,jendp
534 DO i=istrp,iendp
535 IF ((rmask(i-1,j ).gt.0.5_r8).and. &
536 & (rmask(i ,j ).gt.0.5_r8).and. &
537 & (rmask(i-1,j-1).gt.0.5_r8).and. &
538 & (rmask(i ,j-1).gt.0.5_r8)) THEN
539 pmask(i,j)=1.0_r8
540 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. &
541 & (rmask(i ,j ).gt.0.5_r8).and. &
542 & (rmask(i-1,j-1).gt.0.5_r8).and. &
543 & (rmask(i ,j-1).gt.0.5_r8)) THEN
544 pmask(i,j)=cff1
545 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. &
546 & (rmask(i ,j ).lt.0.5_r8).and. &
547 & (rmask(i-1,j-1).gt.0.5_r8).and. &
548 & (rmask(i ,j-1).gt.0.5_r8)) THEN
549 pmask(i,j)=cff1
550 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. &
551 & (rmask(i ,j ).gt.0.5_r8).and. &
552 & (rmask(i-1,j-1).lt.0.5_r8).and. &
553 & (rmask(i ,j-1).gt.0.5_r8)) THEN
554 pmask(i,j)=cff1
555 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. &
556 & (rmask(i ,j ).gt.0.5_r8).and. &
557 & (rmask(i-1,j-1).gt.0.5_r8).and. &
558 & (rmask(i ,j-1).lt.0.5_r8)) THEN
559 pmask(i,j)=cff1
560 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. &
561 & (rmask(i ,j ).lt.0.5_r8).and. &
562 & (rmask(i-1,j-1).gt.0.5_r8).and. &
563 & (rmask(i ,j-1).lt.0.5_r8)) THEN
564 pmask(i,j)=cff2
565 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. &
566 & (rmask(i ,j ).gt.0.5_r8).and. &
567 & (rmask(i-1,j-1).lt.0.5_r8).and. &
568 & (rmask(i ,j-1).gt.0.5_r8)) THEN
569 pmask(i,j)=cff2
570 ELSE IF ((rmask(i-1,j ).gt.0.5_r8).and. &
571 & (rmask(i ,j ).gt.0.5_r8).and. &
572 & (rmask(i-1,j-1).lt.0.5_r8).and. &
573 & (rmask(i ,j-1).lt.0.5_r8)) THEN
574 pmask(i,j)=cff2
575 ELSE IF ((rmask(i-1,j ).lt.0.5_r8).and. &
576 & (rmask(i ,j ).lt.0.5_r8).and. &
577 & (rmask(i-1,j-1).gt.0.5_r8).and. &
578 & (rmask(i ,j-1).gt.0.5_r8)) THEN
579 pmask(i,j)=cff2
580 ELSE
581 pmask(i,j)=0.0_r8
582 END IF
583 END DO
584 END DO
585!
586! Exchange boundary information.
587!
588 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
589 CALL exchange_p2d_tile (ng, tile, &
590 & lbi, ubi, lbj, ubj, &
591 & pmask)
592 END IF
593
594# ifdef DISTRIBUTE
595 CALL mp_exchange2d (ng, tile, model, 1, &
596 & lbi, ubi, lbj, ubj, &
597 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
598 & pmask)
599# endif
600#endif
601!
602!-----------------------------------------------------------------------
603! Compute cosine and sine of grid rotation angle.
604!-----------------------------------------------------------------------
605!
606 DO j=jstrt,jendt
607 DO i=istrt,iendt
608 cosangler(i,j)=cos(angler(i,j))
609 sinangler(i,j)=sin(angler(i,j))
610 END DO
611 END DO
612!
613! Exchange boundary information.
614!
615 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
616 CALL exchange_r2d_tile (ng, tile, &
617 & lbi, ubi, lbj, ubj, &
618 & cosangler)
619 CALL exchange_r2d_tile (ng, tile, &
620 & lbi, ubi, lbj, ubj, &
621 & sinangler)
622 END IF
623
624#ifdef DISTRIBUTE
625 CALL mp_exchange2d (ng, tile, model, 2, &
626 & lbi, ubi, lbj, ubj, &
627 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
628 & cosangler, sinangler)
629#endif
630
631#if defined TIDE_GENERATING_FORCES
632!
633!-----------------------------------------------------------------------
634! Compute squared cosine of latitude and sine of twice latitude.
635!-----------------------------------------------------------------------
636!
637 DO j=jstrt,jendt
638 DO i=istrt,iendt
639 phi=latr(i,j)*deg2rad
640 cosphi=cos(phi)
641 cos2lat(i,j)=cosphi*cosphi ! COS2(latr)
642 sinlat2(i,j)=sin(2.0_r8*phi) ! SIN(2*latr)
643 END DO
644 END DO
645!
646! Exchange boundary information.
647!
648 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
649 CALL exchange_r2d_tile (ng, tile, &
650 & lbi, ubi, lbj, ubj, &
651 & cos2lat)
652 CALL exchange_r2d_tile (ng, tile, &
653 & lbi, ubi, lbj, ubj, &
654 & sinlat2)
655 END IF
656
657# ifdef DISTRIBUTE
658 CALL mp_exchange2d (ng, tile, model, 2, &
659 & lbi, ubi, lbj, ubj, &
660 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
661 & cos2lat, sinlat2)
662# endif
663#endif
664
665#if defined DIFF_GRID || defined DIFF_3DCOEF || \
666 defined visc_grid || defined visc_3dcoef
667!
668!-----------------------------------------------------------------------
669! Determine the squared root of the area of each grid cell used to
670! rescale horizontal mixing by the grid size.
671!-----------------------------------------------------------------------
672!
673 cff=0.0_r8
674 DO j=jstrt,jendt
675 DO i=istrt,iendt
676 grdscl(i,j)=sqrt(om_r(i,j)*on_r(i,j))
677 END DO
678 END DO
679!
680! Exchange boundary information.
681!
682 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
683 CALL exchange_r2d_tile (ng, tile, &
684 & lbi, ubi, lbj, ubj, &
685 & grdscl)
686 END IF
687
688# ifdef DISTRIBUTE
689 CALL mp_exchange2d (ng, tile, model, 1, &
690 & lbi, ubi, lbj, ubj, &
691 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
692 & grdscl)
693# endif
694#endif
695
696#if defined DIFF_3DCOEF || defined VISC_3DCOEF
697!
698!-----------------------------------------------------------------------
699! Compute time invariant, horizontal mixing coefficient using grid
700! scale. Following Holland et (1998), Webb et al. (1998), Griffies
701! and Hallberg (2000), and Lee et al. (2002), the horizontal mixing
702! coefficient can be estimated as:
703!
704! Hmixing = 1/12 * Uscale * grdscl (Harmonic)
705! Bmixing = 1/12 * Uscale * grdscl**3 (Biharmonic)
706!-----------------------------------------------------------------------
707!
708# if defined DIFF_3DCOEF
709 my_diffmin= large
710 my_diffmax=-large
711# endif
712# if defined VISC_3DCOEF
713 my_viscmin= large
714 my_viscmax=-large
715# endif
716 DO j=jstrt,jendt
717 DO i=istrt,iendt
718# if defined DIFF_3DCOEF
719# if defined TS_DIF2
720 hdiffusion(i,j)=pecletcoef*uscale*grdscl(i,j)
721# elif defined TS_DIF4
722 hdiffusion(i,j)=pecletcoef*uscale*grdscl(i,j)**3
723# endif
724 my_diffmin=min(my_diffmin, hdiffusion(i,j))
725 my_diffmax=max(my_diffmax, hdiffusion(i,j))
726# endif
727# if defined VISC_3DCOEF
728# if defined UV_VIS2
729 hviscosity(i,j)=pecletcoef*uscale*grdscl(i,j)
730# elif defined UV_VIS4
731 hviscosity(i,j)=pecletcoef*uscale*grdscl(i,j)**3
732# endif
733 my_viscmin=min(my_viscmin, hviscosity(i,j))
734 my_viscmax=max(my_viscmax, hviscosity(i,j))
735# endif
736 END DO
737 END DO
738!
739! Exchange boundary information.
740!
741 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
742# ifdef DIFF_3DCOEF
743 CALL exchange_r2d_tile (ng, tile, &
744 & lbi, ubi, lbj, ubj, &
745 & hdiffusion)
746# endif
747# ifdef VISC_3DCOEF
748 CALL exchange_r2d_tile (ng, tile, &
749 & lbi, ubi, lbj, ubj, &
750 & hviscosity)
751# endif
752 END IF
753
754# ifdef DISTRIBUTE
755# ifdef DIFF_3DCOEF
756 CALL mp_exchange2d (ng, tile, model, 1, &
757 & lbi, ubi, lbj, ubj, &
758 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
759 & hdiffusion)
760# endif
761# ifdef VISC_3DCOEF
762 CALL mp_exchange2d (ng, tile, model, 1, &
763 & lbi, ubi, lbj, ubj, &
764 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
765 & hviscosity)
766# endif
767# endif
768#endif
769!
770!-----------------------------------------------------------------------
771! Compute minimum and maximum grid spacing.
772!-----------------------------------------------------------------------
773#ifdef SOLVE3D
774!
775! Compute time invariant depths (use zero free-surface).
776!
777 DO i=lbi,ubi
778 DO j=lbj,ubj
779 a2d(i,j)=0.0_r8
780 END DO
781 END DO
782
783 CALL set_depth_tile (ng, tile, inlm, &
784 & lbi, ubi, lbj, ubj, &
785 & imins, imaxs, jmins, jmaxs, &
786 & nstp, nnew, &
787 & h, &
788# ifdef ICESHELF
789 & zice, &
790# endif
791 & a2d, &
792 & hz, z_r, z_w)
793#endif
794!
795! Compute grid spacing range.
796!
797 my_dxmin= large
798 my_dxmax=-large
799 my_dymin= large
800 my_dymax=-large
801#ifdef MASKING
802 my_dxminw= large
803 my_dxmaxw=-large
804 my_dyminw= large
805 my_dymaxw=-large
806#endif
807#ifdef SOLVE3D
808 my_dzmin= large
809 my_dzmax=-large
810# ifdef MASKING
811 my_dzminw= large
812 my_dzmaxw=-large
813# endif
814#endif
815 my_grdmax=-large
816 DO j=jstrt,jendt
817 DO i=istrt,iendt
818#if defined VISC_GRID || defined DIFF_GRID
819 cff=grdscl(i,j)
820#else
821 cff=sqrt(om_r(i,j)*on_r(i,j))
822#endif
823 my_dxmin=min(my_dxmin,om_r(i,j))
824 my_dxmax=max(my_dxmax,om_r(i,j))
825 my_dymin=min(my_dymin,on_r(i,j))
826 my_dymax=max(my_dymax,on_r(i,j))
827#ifdef MASKING
828 IF (rmask(i,j).gt.0.0_r8) THEN
829 my_grdmax=max(my_grdmax,cff)
830 my_dxminw=min(my_dxminw,om_r(i,j))
831 my_dxmaxw=max(my_dxmaxw,om_r(i,j))
832 my_dyminw=min(my_dyminw,on_r(i,j))
833 my_dymaxw=max(my_dymaxw,on_r(i,j))
834 END IF
835#else
836 my_grdmax=max(my_grdmax,cff)
837#endif
838#ifdef SOLVE3D
839 DO k=1,n(ng)
840 my_dzmin=min(my_dzmin,hz(i,j,k))
841 my_dzmax=max(my_dzmax,hz(i,j,k))
842# ifdef MASKING
843 IF (rmask(i,j).gt.0.0_r8) THEN
844 my_dzminw=min(my_dzminw,hz(i,j,k))
845 my_dzmaxw=max(my_dzmaxw,hz(i,j,k))
846 END IF
847# endif
848 END DO
849#endif
850 END DO
851 END DO
852!
853!-----------------------------------------------------------------------
854! Compute gravity waves Courant number.
855!-----------------------------------------------------------------------
856!
857! The 2D Courant number is defined as:
858!
859! Cg = c * dt * SQRT (1/dx^2 + 1/dy^2)
860!
861! where c=SQRT(g*h) is gravity wave speed, and dx, dy are grid spacing
862! in each direction.
863!
864 my_cg_min= large
865 my_cg_max=-large
866 my_cg_cor=-large
867
868 DO j=jstrt,jendt
869 DO i=istrt,iendt
870#ifdef MASKING
871 IF (rmask(i,j).gt.0.0_r8) THEN
872#endif
873#ifdef ICESHELF
874 cff=dtfast(ng)* &
875 & sqrt(g*abs(h(i,j)-abs(zice(i,j)))* &
876 & (pm(i,j)*pm(i,j)+pn(i,j)*pn(i,j)))
877#else
878 cff=dtfast(ng)* &
879 & sqrt(g*abs(h(i,j))*(pm(i,j)*pm(i,j)+pn(i,j)*pn(i,j)))
880#endif
881 my_cg_min=min(my_cg_min,cff)
882 my_cg_max=max(my_cg_max,cff)
883
884 cff=dt(ng)*abs(f(i,j))
885 my_cg_cor=max(my_cg_cor,cff)
886
887#ifdef MASKING
888 END IF
889#endif
890 END DO
891 END DO
892!
893!-----------------------------------------------------------------------
894! Perform global reductions.
895!-----------------------------------------------------------------------
896!
897#ifdef DISTRIBUTE
898 nsub=1 ! distributed-memory
899#else
900 IF (domain(ng)%SouthWest_Corner(tile).and. &
901 & domain(ng)%NorthEast_Corner(tile)) THEN
902 nsub=1 ! non-tiled application
903 ELSE
904 nsub=ntilex(ng)*ntilee(ng) ! tiled application
905 END IF
906#endif
907!$OMP CRITICAL (REDUCTIONS)
908 cg_min(ng)=min(cg_min(ng),my_cg_min)
909 cg_max(ng)=max(cg_max(ng),my_cg_max)
910 cg_cor(ng)=max(cg_cor(ng),my_cg_cor)
911 grdmax(ng)=max(grdmax(ng),my_grdmax)
912 dxmin(ng)=min(dxmin(ng),my_dxmin)
913 dxmax(ng)=max(dxmax(ng),my_dxmax)
914 dymin(ng)=min(dymin(ng),my_dymin)
915 dymax(ng)=max(dymax(ng),my_dymax)
916#ifdef MASKING
917 dxminw(ng)=min(dxminw(ng),my_dxminw)
918 dxmaxw(ng)=max(dxmaxw(ng),my_dxmaxw)
919 dyminw(ng)=min(dyminw(ng),my_dyminw)
920 dymaxw(ng)=max(dymaxw(ng),my_dymaxw)
921#endif
922#ifdef SOLVE3D
923 dzmin(ng)=min(dzmin(ng),my_dzmin)
924 dzmax(ng)=max(dzmax(ng),my_dzmax)
925# ifdef MASKING
926 dzminw(ng)=min(dzminw(ng),my_dzminw)
927 dzmaxw(ng)=max(dzmaxw(ng),my_dzmaxw)
928# endif
929#endif
930#ifdef DIFF_3DCOEF
931 diffmin(ng)=min(diffmin(ng),my_diffmin)
932 diffmax(ng)=max(diffmax(ng),my_diffmax)
933#endif
934#ifdef VISC_3DCOEF
935 viscmin(ng)=min(viscmin(ng),my_viscmin)
936 viscmax(ng)=max(viscmax(ng),my_viscmax)
937#endif
939 IF (tile_count.eq.nsub) THEN
940 tile_count=0
941#ifdef DISTRIBUTE
942 rbuffer(1)=cg_min(ng)
943 op_handle(1)='MIN'
944 rbuffer(2)=cg_max(ng)
945 op_handle(2)='MAX'
946 rbuffer(3)=cg_cor(ng)
947 op_handle(3)='MAX'
948 rbuffer(4)=grdmax(ng)
949 op_handle(4)='MAX'
950 rbuffer(5)=dxmin(ng)
951 op_handle(5)='MIN'
952 rbuffer(6)=dxmax(ng)
953 op_handle(6)='MAX'
954 rbuffer(7)=dymin(ng)
955 op_handle(7)='MIN'
956 rbuffer(8)=dymax(ng)
957 op_handle(8)='MAX'
958# ifdef SOLVE3D
959 rbuffer(9)=dzmin(ng)
960 op_handle(9)='MIN'
961 rbuffer(10)=dzmax(ng)
962 op_handle(10)='MAX'
963# else
964 rbuffer(9)=0.0_r8
965 op_handle(9)='MIN'
966 rbuffer(10)=0.0_r8
967 op_handle(10)='MAX'
968# endif
969# ifdef DIFF_3DCOEF
970 rbuffer(11)=diffmin(ng)
971 op_handle(11)='MIN'
972 rbuffer(12)=diffmax(ng)
973 op_handle(12)='MAX'
974# else
975 rbuffer(11)=0.0_dp
976 op_handle(11)='MIN'
977 rbuffer(12)=0.0_dp
978 op_handle(12)='MAX'
979# endif
980# ifdef VISC_3DCOEF
981 rbuffer(13)=viscmin(ng)
982 op_handle(13)='MIN'
983 rbuffer(14)=viscmax(ng)
984 op_handle(14)='MAX'
985# else
986 rbuffer(13)=0.0_dp
987 op_handle(13)='MIN'
988 rbuffer(14)=0.0_dp
989 op_handle(14)='MAX'
990# endif
991# ifdef MASKING
992 rbuffer(15)=dxminw(ng)
993 op_handle(15)='MIN'
994 rbuffer(16)=dxmaxw(ng)
995 op_handle(16)='MAX'
996 rbuffer(17)=dyminw(ng)
997 op_handle(17)='MIN'
998 rbuffer(18)=dymaxw(ng)
999 op_handle(18)='MAX'
1000# ifdef SOLVE3D
1001 rbuffer(19)=dzminw(ng)
1002 op_handle(19)='MIN'
1003 rbuffer(20)=dzmaxw(ng)
1004 op_handle(20)='MAX'
1005# else
1006 rbuffer(19)=0.0_dp
1007 op_handle(19)='MIN'
1008 rbuffer(20)=0.0_dp
1009 op_handle(20)='MAX'
1010# endif
1011 CALL mp_reduce (ng, model, 20, rbuffer, op_handle)
1012# else
1013 CALL mp_reduce (ng, model, 14, rbuffer, op_handle)
1014# endif
1015 cg_min(ng)=rbuffer(1)
1016 cg_max(ng)=rbuffer(2)
1017 cg_cor(ng)=rbuffer(3)
1018 grdmax(ng)=rbuffer(4)
1019 dxmin(ng)=rbuffer(5)
1020 dxmax(ng)=rbuffer(6)
1021 dymin(ng)=rbuffer(7)
1022 dymax(ng)=rbuffer(8)
1023# ifdef SOLVE3D
1024 dzmin(ng)=rbuffer(9)
1025 dzmax(ng)=rbuffer(10)
1026# endif
1027# ifdef DIFF_3DCOEF
1028 diffmin(ng)=rbuffer(11)
1029 diffmax(ng)=rbuffer(12)
1030# endif
1031# ifdef VISC_3DCOEF
1032 viscmin(ng)=rbuffer(13)
1033 viscmax(ng)=rbuffer(14)
1034# endif
1035# ifdef MASKING
1036 dxminw(ng)=rbuffer(15)
1037 dxmaxw(ng)=rbuffer(16)
1038 dyminw(ng)=rbuffer(17)
1039 dymaxw(ng)=rbuffer(18)
1040# ifdef SOLVE3D
1041 dzminw(ng)=rbuffer(19)
1042 dzmaxw(ng)=rbuffer(20)
1043# endif
1044# endif
1045#endif
1046 IF (master.and.lwrtinfo(ng)) THEN
1047 WRITE (stdout,10) ng, &
1048#ifdef MASKING
1049 & dxmin(ng)*0.001_dp, dxminw(ng)*0.001_dp, &
1050 & dxmax(ng)*0.001_dp, dxmaxw(ng)*0.001_dp, &
1051 & dymin(ng)*0.001_dp, dyminw(ng)*0.001_dp, &
1052 & dymax(ng)*0.001_dp, dymaxw(ng)*0.001_dp
1053 10 FORMAT (/,' Metrics information for Grid ',i2.2,':', &
1054 & /,' ===============================',/, &
1055 & /,' Minimum X-grid spacing, DXmin = ',1pe15.8,' km', &
1056 & 4x,'Water points = ',1pe15.8,' km', &
1057 & /,' Maximum X-grid spacing, DXmax = ',1pe15.8,' km', &
1058 & 4x,'Water points = ',1pe15.8,' km', &
1059 & /,' Minimum Y-grid spacing, DYmin = ',1pe15.8,' km', &
1060 & 4x,'Water points = ',1pe15.8,' km', &
1061 & /,' Maximum Y-grid spacing, DYmax = ',1pe15.8,' km', &
1062 & 4x,'Water points = ',1pe15.8,' km')
1063#else
1064 & dxmin(ng)*0.001_dp, &
1065 & dxmax(ng)*0.001_dp, &
1066 & dymin(ng)*0.001_dp, &
1067 & dymax(ng)*0.001_dp
1068 10 FORMAT (/,' Metrics information for Grid ',i2.2,':', &
1069 & /,' ===============================',/, &
1070 & /,' Minimum X-grid spacing, DXmin = ',1pe15.8,' km', &
1071 & /,' Maximum X-grid spacing, DXmax = ',1pe15.8,' km', &
1072 & /,' Minimum Y-grid spacing, DYmin = ',1pe15.8,' km', &
1073 & /,' Maximum Y-grid spacing, DYmax = ',1pe15.8,' km')
1074#endif
1075#ifdef SOLVE3D
1076# ifdef MASKING
1077 WRITE (stdout,20) dzmin(ng), dzminw(ng), dzmax(ng), dzmaxw(ng)
1078 20 FORMAT (' Minimum Z-grid spacing, DZmin = ',1pe15.8,' m', &
1079 & 5x,'Water points = ',1pe15.8,' m',/, &
1080 & ' Maximum Z-grid spacing, DZmax = ',1pe15.8,' m', &
1081 & 5x,'Water points = ',1pe15.8,' m')
1082# else
1083 WRITE (stdout,20) dzmin(ng), dzmax(ng)
1084 20 FORMAT (' Minimum Z-grid spacing, DZmin = ',1pe15.8,' m',/, &
1085 & ' Maximum Z-grid spacing, DZmax = ',1pe15.8,' m')
1086# endif
1087#endif
1088 WRITE (stdout,30) cg_min(ng), cg_max(ng), cg_cor(ng)
1089 30 FORMAT (/,' Minimum barotropic Courant Number = ', 1pe15.8,/, &
1090 & ' Maximum barotropic Courant Number = ', 1pe15.8,/, &
1091 & ' Maximum Coriolis Courant Number = ', 1pe15.8,/)
1092# if defined VISC_GRID || defined DIFF_GRID
1093 WRITE (stdout,40) grdmax(ng)/1000.0_r8
1094 40 FORMAT (' Horizontal mixing scaled by grid area squared root',&
1095# ifdef MASKING
1096 & ', MAXVAL(grdscl) = ',1pe15.8,' km', &
1097 & 2x,'(Water points)')
1098# else
1099 & ', MAXVAL(grdscl) = ',1pe15.8,' km')
1100# endif
1101# endif
1102# ifdef DIFF_3DCOEF
1103# ifdef TS_DIF2
1104 units='m2/s'
1105# elif defined TS_DIF4
1106 units='m4/s'
1107# endif
1108 WRITE (stdout,50) diffmin(ng), trim(units), &
1109 & diffmax(ng), trim(units)
1110 50 FORMAT (/,' Minimum horizontal diffusion coefficient = ', &
1111 & 1pe15.8,1x,a, &
1112 & /,' Maximum horizontal diffusion coefficient = ', &
1113 & 1pe15.8,1x,a)
1114# endif
1115# ifdef VISC_3DCOEF
1116# ifdef UV_VIS2
1117 units='m2/s'
1118# elif defined UV_VIS4
1119 units='m4/s'
1120# endif
1121 WRITE (stdout,60) viscmin(ng), trim(units), &
1122 & viscmax(ng), trim(units)
1123 60 FORMAT (/,' Minimum horizontal viscosity coefficient = ', &
1124 & 1pe15.8,1x,a, &
1125 & /,' Maximum horizontal viscosity coefficient = ', &
1126 & 1pe15.8,1x,a)
1127# endif
1128 END IF
1129 END IF
1130!$OMP END CRITICAL (REDUCTIONS)
1131
1132#ifdef NESTING
1133!
1134!-----------------------------------------------------------------------
1135! If grid refinement, set various important parameters.
1136!-----------------------------------------------------------------------
1137!
1138 IF (ng.eq.ngrids) THEN
1139!
1140! Set coarser donor grid number to finer grid external contact points.
1141! The reciver grid is always finer than the donor grid, DXmax(dg) >
1142! DXmax(rg). It is done in terms of the reciver grid.
1143!
1144 DO cr=1,ncontact
1145 dg=donor_grid(cr)
1146 rg=receiver_grid(cr)
1147 IF (refinedgrid(rg).and. &
1148 & (refinescale(rg).gt.0).and. &
1149 & dxmax(dg).gt.dxmax(rg)) THEN
1150 coarserdonor(rg)=dg
1151 END IF
1152 END DO
1153!
1154! Set donor finer grid number to receiver coarser grid. This is use for
1155! two-way feeback between grids. The donor grid is always finer than the
1156! receiver grid, DXmax(dg) < DXmax(rg). It is done in terms of the donor
1157! grid.
1158!
1159 DO cr=1,ncontact
1160 dg=donor_grid(cr)
1161 rg=receiver_grid(cr)
1162 IF (refinedgrid(dg).and. &
1163 & (refinescale(dg).gt.0).and. &
1164 & dxmax(dg).gt.dxmax(rg)) THEN
1165 finerdonor(dg)=rg
1166 END IF
1167 END DO
1168!
1169! Set logical indicating which coarser grid is a donor to a finer
1170! grid external contact points.
1171!
1172 DO dg=1,ngrids
1173 DO ig=1,ngrids
1174 IF (coarserdonor(ig).eq.dg) THEN
1175 donortofiner(dg)=.true.
1176 END IF
1177 END DO
1178 END DO
1179!
1180! Set switch indicating which refined grid(s), with RefineScale(dg) > 0,
1181! include finer refined grids inside (telescoping refinement).
1182!
1183 DO cr=1,ncontact
1184 dg=donor_grid(cr)
1185 rg=receiver_grid(cr)
1186 IF ((refinescale(dg).gt.0).and.(coarserdonor(rg).eq.dg)) THEN
1187 telescoping(dg)=.true.
1188 END IF
1189 END DO
1190!
1191! Set Number of refined grid time-steps using the ratio between donor
1192! and receiver grids. Optimally, thes values should be based on the
1193! spatial refinement ratio for stability. However, the user is allowed
1194! to take larger divisible time-step with respect to the donor grid.
1195! The user is responsible to set the appropriate refined grid time-step
1196! for stability.
1197!
1198 DO cr=1,ncontact
1199 dg=donor_grid(cr)
1200 rg=receiver_grid(cr)
1201 IF (refinedgrid(rg).and.(coarserdonor(rg).eq.dg)) THEN
1202 refinesteps(rg)=nint(dt(dg)/dt(rg))
1203 END IF
1204 END DO
1205!
1206! Report information.
1207!
1208 IF (domain(ng)%NorthEast_Test(tile)) THEN
1209 IF (master.and.any(refinedgrid)) THEN
1210 WRITE (stdout,70)
1211 70 FORMAT (/,' Refined Nested Grid(s) Information: ', &
1212 & /,' ==================================',/,/, &
1213 & ' Refined Donor Refined Timestep Refined', &
1214 & /,' Grid Grid Scale Ratio ', &
1215 & 'Timesteps',/)
1216 DO ig=1,ngrids
1217 IF (refinescale(ig).gt.0) THEN
1218 dg=coarserdonor(ig)
1219 WRITE (stdout,80) ig, dg, refinescale(ig), &
1220 & dt(dg)/dt(ig), refinesteps(ig)
1221 80 FORMAT(4x,i2.2,8x,i2.2,7x,i2.2,4x,f9.5,6x,i2.2)
1222 END IF
1223 END DO
1224 WRITE (stdout,90)
1225 90 FORMAT (/,' WARNING: Usually the number of Refined ', &
1226 & 'Timesteps must be the same',/,10x, &
1227 & 'as the Refined Scale for numerical stability.',/)
1228 END IF
1229 END IF
1230!
1231! Check refined grid time-step. The time-step size for refined grids
1232! needs to be an exact mode multiple of its coarser donor grid. In
1233! principle, it can be a "RefineScale" factor smaller. However, other
1234! integer smaller or larger factor is allowed such that:
1235!
1236! MOD(dt(dg)*SecScale, dt(rg)*SecScale) = 0 dg: donor coarse grid
1237! rg: receiver fine grid
1238!
1239! Notice that SecScale is used to avoid roundoff when the timestep
1240! between donor and receiver grids are small and less than one.
1241!
1242 secscale=1000.0_dp ! seconds to milliseconds
1243 DO ig=1,ngrids
1244 IF (refinedgrid(ig).and.(refinescale(ig).gt.0)) THEN
1245 dg=coarserdonor(ig)
1246 IF (mod(dt(dg)*secscale,dt(ig)*secscale).ne.0.0_dp) THEN
1247 IF (domain(ng)%SouthWest_Test(tile)) THEN
1248 IF (master) THEN
1249 WRITE (stdout,100) ig, dt(ig), dg, dt(dg), &
1250 & mod(dt(dg),dt(ig))
1251 100 FORMAT (/,' METRICS - illegal timestep size for ', &
1252 & ' finer reciever grid (rg),',/,11x, &
1253 & 'It must be an exact divisible factor from ', &
1254 & 'donor grid (dg).',/,/,11x,'rg = ',i2.2, &
1255 & ', dt(rg) = ',f12.4,/,11x,'dg = ',i2.2, &
1256 & ', dt(dg) = ',f12.4,4x, &
1257 & 'MOD[dt(dg), dt(rg)] = ',f10.4,/)
1258 END IF
1259 END IF
1260 exit_flag=5
1261 RETURN
1262 END IF
1263 END IF
1264 END DO
1265 END IF
1266#endif
1267
1268#ifdef FOUR_DVAR
1269!
1270!-----------------------------------------------------------------------
1271! Spatial convolution parameters.
1272!-----------------------------------------------------------------------
1273!
1274! Determine diffusion operator time-step size using FTCS stability
1275! criteria:
1276!
1277! Kh DTsizeH / MIN(DXmin,DYmin)^2 < Hgamma / 4
1278!
1279! Kv DTsizeH / DZmin^2 < Vgamma / 2
1280!
1281! where a Hgamma and Vgamma are used to scale the time-step below
1282! its theoretical limit for stability and accurary.
1283!
1284 cff=min(dxmin(ng),dymin(ng))
1285 DO rec=1,2
1286 DO is=1,nstatevar(ng)
1287# ifdef SOLVE3D
1288 IF (is.le.(5+nt(ng))) THEN
1289# else
1290 IF (is.le.3) THEN
1291# endif
1292 dtsizeh(rec,is)=hgamma(rec)*cff*cff/(4.0_r8*khmax(ng))
1293 ELSE ! use stability factor for surface forcing
1294 dtsizeh(rec,is)=hgamma(4)*cff*cff/(4.0_r8*khmax(ng))
1295 END IF
1296# ifdef SOLVE3D
1297# ifdef IMPLICIT_VCONV
1298 dtsizev(rec,is)=vgamma(rec)*dzmax(ng)*dzmax(ng)/ &
1299 & (2.0_r8*kvmax(ng))
1300# else
1301 dtsizev(rec,is)=vgamma(rec)*dzmin(ng)*dzmin(ng)/ &
1302 & (2.0_r8*kvmax(ng))
1303# endif
1304# endif
1305 END DO
1306 END DO
1307
1308# ifdef ADJUST_BOUNDARY
1309!
1310! Open boundaries convolutions. Recall that the horizontal convolutions
1311! are one-dimensional so a factor of 2 is used instead.
1312!
1313 DO ibry=1,4
1314 DO is=1,nstatevar(ng)
1315 IF ((ibry.eq.iwest).or.(ibry.eq.ieast)) THEN
1316 dtsizehb(ibry,is)=hgamma(3)*dymin(ng)*dymin(ng)/ &
1317 & (2.0_r8*khmax(ng))
1318 ELSE IF ((ibry.eq.isouth).or.(ibry.eq.inorth)) THEN
1319 dtsizehb(ibry,is)=hgamma(3)*dxmin(ng)*dxmin(ng)/ &
1320 & (2.0_r8*khmax(ng))
1321 END IF
1322# ifdef SOLVE3D
1323# ifdef IMPLICIT_VCONV
1324 dtsizevb(ibry,is)=vgamma(3)*dzmax(ng)*dzmax(ng)/ &
1325 & (2.0_r8*kvmax(ng))
1326# else
1327 dtsizevb(ibry,is)=vgamma(3)*dzmin(ng)*dzmin(ng)/ &
1328 & (2.0_r8*kvmax(ng))
1329# endif
1330# endif
1331 END DO
1332 END DO
1333# endif
1334!
1335!-----------------------------------------------------------------------
1336! Determine FULL number of integeration steps as function of the
1337! spatial decorrelation scale (Hdecay, Vdecay). Notice that in
1338! the squared-root filter only half of number of step is used.
1339! Thefore, the number of steps are forced to be even.
1340!-----------------------------------------------------------------------
1341!
1342! Initial conditions, surface forcing and model error covariance.
1343!
1344 DO rec=1,2
1345 DO is=1,nstatevar(ng)
1346 IF (hdecay(rec,is,ng).gt.0.0_r8) THEN
1347 nhsteps(rec,is)=nint(hdecay(rec,is,ng)* &
1348 & hdecay(rec,is,ng)/ &
1349 & (2.0_r8*khmax(ng)*dtsizeh(rec,is)))
1350 IF (mod(nhsteps(rec,is),2).ne.0) THEN
1351 nhsteps(rec,is)=nhsteps(rec,is)+1
1352 END IF
1353 END IF
1354# ifdef SOLVE3D
1355 IF (vdecay(rec,is,ng).gt.0.0_r8) THEN
1356 nvsteps(rec,is)=nint(vdecay(rec,is,ng)* &
1357 & vdecay(rec,is,ng)/ &
1358 & (2.0_r8*kvmax(ng)*dtsizev(rec,is)))
1359# ifdef IMPLICIT_VCONV
1360 nvsteps(rec,is)=max(1,nvsteps(rec,is))
1361# endif
1362 IF (mod(nvsteps(rec,is),2).ne.0) THEN
1363 nvsteps(rec,is)=nvsteps(rec,is)+1
1364 END IF
1365 END IF
1366# endif
1367 END DO
1368 END DO
1369
1370# ifdef ADJUST_BOUNDARY
1371!
1372! Boundary conditions model error covariance.
1373!
1374 DO ibry=1,4
1375 DO is=1,nstatevar(ng)
1376 IF (hdecayb(is,ibry,ng).gt.0.0_r8) THEN
1377 nhstepsb(ibry,is)=nint(hdecayb(is,ibry,ng)* &
1378 & hdecayb(is,ibry,ng)/ &
1379 & (2.0_r8*khmax(ng)*dtsizehb(ibry,is)))
1380 IF (mod(nhstepsb(ibry,is),2).ne.0) THEN
1381 nhstepsb(ibry,is)=nhstepsb(ibry,is)+1
1382 END IF
1383 END IF
1384# ifdef SOLVE3D
1385 IF (vdecayb(is,ibry,ng).gt.0.0_r8) THEN
1386 nvstepsb(ibry,is)=nint(vdecayb(is,ibry,ng)* &
1387 & vdecayb(is,ibry,ng)/ &
1388 & (2.0_r8*kvmax(ng)*dtsizevb(ibry,is)))
1389# ifdef IMPLICIT_VCONV
1390 nvstepsb(ibry,is)=max(1,nvstepsb(ibry,is))
1391# endif
1392 IF (mod(nvstepsb(ibry,is),2).ne.0) THEN
1393 nvstepsb(ibry,is)=nvstepsb(ibry,is)+1
1394 END IF
1395# endif
1396 END IF
1397 END DO
1398 END DO
1399# endif
1400!
1401! Report convolution parameters.
1402!
1403 IF (master.and.lwrtinfo(ng)) THEN
1404 DO rec=1,nsa
1405 DO is=1,nstatevar(ng)
1406 IF (rec.eq.1) THEN
1407 text='Horizontal convolution, NHstepsI, DTsizeH ='
1408 ELSE IF (rec.eq.2) THEN
1409 text='Horizontal convolution, NHstepsM, DTsizeH ='
1410# ifdef SOLVE3D
1411 ELSE IF (is.gt.(5+nt(ng))) THEN
1412# else
1413 ELSE IF (is.gt.3) THEN
1414# endif
1415 text='Horizontal convolution, NHstepsF, DTsizeH ='
1416 END IF
1417 IF (hdecay(rec,is,ng).gt.0.0_r8) THEN
1418 WRITE (stdout,110) trim(text), nhsteps(rec,is), &
1419 & dtsizeh(rec,is), &
1420 & trim(vname(1,idsvar(is)))
1421 END IF
1422 END DO
1423 WRITE (stdout,120)
1424 END DO
1425# if defined SOLVE3D && defined VCONVOLUTION
1426 DO rec=1,nsa
1427 DO is=1,nstatevar(ng)
1428 IF (rec.eq.1) THEN
1429 text='Vertical convolution, NVstepsI, DTsizeV ='
1430 ELSE IF (rec.eq.2) THEN
1431 text='Vertical convolution, NVstepsM, DTsizeV ='
1432# ifdef SOLVE3D
1433 ELSE IF (is.gt.(5+nt(ng))) THEN
1434# else
1435 ELSE IF (is.gt.3) THEN
1436# endif
1437 text='Vertical convolution, NVstepsF, DTsizeV ='
1438 END IF
1439 IF (vdecay(rec,is,ng).gt.0.0_r8) THEN
1440 WRITE (stdout,110) trim(text), nvsteps(rec,is), &
1441 & dtsizev(rec,is), &
1442 & trim(vname(1,idsvar(is)))
1443 END IF
1444 END DO
1445 WRITE (stdout,120)
1446 END DO
1447# endif
1448# ifdef ADJUST_BOUNDARY
1449 DO is=1,nstatevar(ng)
1450 DO ibry=1,4
1451 IF (lobc(ibry,is,ng)) THEN
1452 IF (ibry.eq.iwest) THEN
1453 text='West bry Hconvolution, NHstepsB, DTsizeHB ='
1454 ELSE IF (ibry.eq.isouth) THEN
1455 text='South bry Hconvolution, NHstepsB, DTsizeHB ='
1456 ELSE IF (ibry.eq.ieast) THEN
1457 text='East bry Hconvolution, NHstepsB, DTsizeHB ='
1458 ELSE IF (ibry.eq.inorth) THEN
1459 text='North bry Hconvolution, NHstepsB, DTsizeHB ='
1460 END IF
1461 IF (hdecayb(is,ibry,ng).gt.0.0_r8) THEN
1462 WRITE (stdout,110) trim(text), nhstepsb(ibry,is), &
1463 & dtsizehb(ibry,is), &
1464 & trim(vname(1,idsbry(is)))
1465 END IF
1466 END IF
1467 END DO
1468 END DO
1469 WRITE (stdout,120)
1470# if defined SOLVE3D && defined VCONVOLUTION
1471 DO is=1,nstatevar(ng)
1472 DO ibry=1,4
1473 IF (lobc(ibry,is,ng)) THEN
1474 IF (ibry.eq.iwest) THEN
1475 text='West bry Vconvolution, NVstepsB, DTsizeVB ='
1476 ELSE IF (ibry.eq.isouth) THEN
1477 text='South bry Vconvolution, NVstepsB, DTsizeVB ='
1478 ELSE IF (ibry.eq.ieast) THEN
1479 text='East bry Vconvolution, NVstepsB, DTsizeVB ='
1480 ELSE IF (ibry.eq.inorth) THEN
1481 text='North bry Vconvolution, NVstepsB, DTsizeVB ='
1482 END IF
1483 IF (vdecayb(is,ibry,ng).gt.0.0_r8) THEN
1484 WRITE (stdout,110) trim(text), nvstepsb(ibry,is), &
1485 & dtsizevb(ibry,is), &
1486 & trim(vname(1,idsbry(is)))
1487 END IF
1488 END IF
1489 END DO
1490 END DO
1491# endif
1492# endif
1493 END IF
1494
1495 110 FORMAT (1x,a,i5,1x,1pe15.8,' s',2x,a)
1496 120 FORMAT (1x)
1497#endif
1498
1499 RETURN
1500 END SUBROUTINE metrics_tile
1501
1502 END MODULE metrics_mod
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_p2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:66
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine, public metrics(ng, tile, model)
Definition metrics.F:24
subroutine metrics_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, f, h, pm, pn, pmask, rmask, angler, cosangler, sinangler, zice, hz, z_r, z_w, grdscl, om_p, om_r, om_u, om_v, on_p, on_r, on_u, on_v, fomn, omn, pnom_p, pnom_r, pnom_u, pnom_v, pmon_p, pmon_r, pmon_u, pmon_v)
Definition metrics.F:135
integer, dimension(:,:), allocatable nhsteps
real(r8), dimension(:), allocatable kvmax
integer, dimension(:,:), allocatable nvsteps
real(r8), dimension(:,:), allocatable dtsizehb
real(r8), dimension(:,:), allocatable dtsizeh
integer, dimension(:,:), allocatable nvstepsb
real(r8), dimension(:,:), allocatable dtsizevb
real(r8), dimension(:), allocatable khmax
real(r8), dimension(:,:), allocatable dtsizev
integer, dimension(:), allocatable nstatevar
integer, dimension(:,:), allocatable nhstepsb
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
integer stdout
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
integer, dimension(:), allocatable idsbry
character(len=maxlen), dimension(6, 0:nv) vname
integer, dimension(:), allocatable idsvar
integer, dimension(:), allocatable refinesteps
logical, dimension(:), allocatable telescoping
logical, dimension(:), allocatable donortofiner
integer, dimension(:), allocatable receiver_grid
integer, dimension(:), allocatable coarserdonor
integer, dimension(:), allocatable finerdonor
integer, dimension(:), allocatable donor_grid
integer ncontact
logical master
integer tile_count
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
integer nghostpoints
Definition mod_param.F:710
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer ngrids
Definition mod_param.F:113
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
integer, dimension(:), allocatable nt
Definition mod_param.F:489
integer nsa
Definition mod_param.F:612
real(dp), dimension(:), allocatable cg_min
real(dp), dimension(:), allocatable dzmaxw
real(dp), dimension(:), allocatable dxmaxw
real(dp), dimension(:), allocatable dzmin
real(dp), dimension(:), allocatable dzminw
real(dp), dimension(:), allocatable dt
logical, dimension(:,:,:), allocatable lobc
logical, dimension(:), allocatable ewperiodic
integer, parameter iwest
logical, dimension(:), allocatable nsperiodic
real(dp), parameter large
real(dp), dimension(:), allocatable dxmin
real(dp), dimension(:), allocatable dymax
real(dp), dimension(:), allocatable dymaxw
real(dp), parameter deg2rad
real(dp), dimension(:), allocatable cg_max
real(dp), dimension(:), allocatable dxminw
real(r8), dimension(:,:,:), allocatable vdecay
integer exit_flag
logical, dimension(:), allocatable lwrtinfo
real(dp), dimension(:), allocatable grdmax
real(r8), dimension(4) hgamma
real(dp), dimension(:), allocatable dzmax
integer, parameter isouth
real(r8), dimension(:,:,:), allocatable hdecayb
real(dp), dimension(:), allocatable dxmax
real(dp), dimension(:), allocatable dtfast
integer, parameter ieast
real(dp) g
logical, dimension(:), allocatable refinedgrid
integer, dimension(:), allocatable refinescale
integer, parameter inorth
real(dp), dimension(:), allocatable dymin
real(dp), dimension(:), allocatable cg_cor
real(r8), dimension(:,:,:), allocatable vdecayb
real(r8), dimension(:,:,:), allocatable hdecay
real(r8), dimension(4) vgamma
real(dp), dimension(:), allocatable dyminw
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable nstp
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public set_depth_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, h, zice, zt_avg1, hz, z_r, z_w)
Definition set_depth.F:86