ROMS
Loading...
Searching...
No Matches
dotproduct.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE dotproduct_mod
3
4#if defined TLM_CHECK || defined PROPAGATOR
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! These routines compute various dot products between nonlinear, !
14! tangent linear, and adjoint state variables. !
15! !
16!=======================================================================
17!
18 implicit none
19!
20 PRIVATE
21# ifdef TLM_CHECK
22 PUBLIC :: ad_dotproduct
23# endif
24# ifdef TLM_CHECK
25 PUBLIC :: nl_dotproduct, tl_dotproduct
26# endif
27# ifdef PROPAGATOR
28# ifdef ADJOINT
29 PUBLIC :: ad_statenorm
30# endif
31# ifdef TANGENT
32 PUBLIC :: tl_statenorm
33# endif
34# endif
35!
36 CONTAINS
37!
38# ifdef TLM_CHECK
39 SUBROUTINE nl_dotproduct (ng, tile, Linp)
40!
41!=======================================================================
42! !
43! This routine computes dot product function (g1) associated with !
44! the nonlinear solution. !
45! !
46!=======================================================================
47!
48 USE mod_param
49# ifdef MASKING
50 USE mod_grid
51# endif
52 USE mod_ocean
53 USE mod_stepping
54!
55! Imported variable declarations.
56!
57 integer, intent(in) :: ng, tile, linp
58!
59! Local variable declarations.
60!
61 character (len=*), parameter :: myfile = &
62 & __FILE__//", nl_dotproduct"
63!
64# include "tile.h"
65!
66# ifdef PROFILE
67 CALL wclock_on (ng, inlm, 7, __line__, myfile)
68# endif
69 CALL nl_dotproduct_tile (ng, tile, &
70 & lbi, ubi, lbj, ubj, &
71 & imins, imaxs, jmins, jmaxs, &
72 & kout, linp, &
73# ifdef SOLVE3D
74 & nout, &
75# endif
76# ifdef MASKING
77 & grid(ng) % rmask, &
78 & grid(ng) % umask, &
79 & grid(ng) % vmask, &
80# endif
81# ifdef SOLVE3D
82 & ocean(ng) % ad_t, &
83 & ocean(ng) % t, &
84 & ocean(ng) % tG, &
85 & ocean(ng) % ad_u, &
86 & ocean(ng) % u, &
87 & ocean(ng) % uG, &
88 & ocean(ng) % ad_v, &
89 & ocean(ng) % v, &
90 & ocean(ng) % vG, &
91# endif
92 & ocean(ng) % ad_ubar, &
93 & ocean(ng) % ubar, &
94 & ocean(ng) % ubarG, &
95 & ocean(ng) % ad_vbar, &
96 & ocean(ng) % vbar, &
97 & ocean(ng) % vbarG, &
98 & ocean(ng) % ad_zeta, &
99 & ocean(ng) % zeta, &
100 & ocean(ng) % zetaG)
101# ifdef PROFILE
102 CALL wclock_off (ng, inlm, 7, __line__, myfile)
103# endif
104!
105 RETURN
106 END SUBROUTINE nl_dotproduct
107!
108!***********************************************************************
109 SUBROUTINE nl_dotproduct_tile (ng, tile, &
110 & LBi, UBi, LBj, UBj, &
111 & IminS, ImaxS, JminS, JmaxS, &
112 & Kinp, Linp, &
113# ifdef SOLVE3D
114 & Ninp, &
115# endif
116# ifdef MASKING
117 & rmask, umask, vmask, &
118# endif
119# ifdef SOLVE3D
120 & ad_t, t, tG, &
121 & ad_u, u, uG, &
122 & ad_v, v, vG, &
123# endif
124 & ad_ubar, ubar, ubarG, &
125 & ad_vbar, vbar, vbarG, &
126 & ad_zeta, zeta, zetaG)
127!***********************************************************************
128!
129 USE mod_param
130 USE mod_parallel
131 USE mod_fourdvar
132 USE mod_ncparam
133 USE mod_iounits
134 USE mod_scalars
135
136# ifdef DISTRIBUTE
137!
138 USE distribute_mod, ONLY : mp_reduce
139# endif
140!
141! Imported variable declarations.
142!
143 integer, intent(in) :: ng, tile
144 integer, intent(in) :: LBi, UBi, LBj, UBj
145 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
146 integer, intent(in) :: Kinp, Linp
147# ifdef SOLVE3D
148 integer, intent(in) :: Ninp
149# endif
150!
151# ifdef ASSUMED_SHAPE
152# ifdef MASKING
153 real(r8), intent(in) :: rmask(LBi:,LBj:)
154 real(r8), intent(in) :: umask(LBi:,LBj:)
155 real(r8), intent(in) :: vmask(LBi:,LBj:)
156# endif
157# ifdef SOLVE3D
158 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
159 real(r8), intent(in) :: t (LBi:,LBj:,:,:,:)
160 real(r8), intent(in) :: tG(LBi:,LBj:,:,:,:)
161
162 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
163 real(r8), intent(in) :: u (LBi:,LBj:,:,:)
164 real(r8), intent(in) :: uG(LBi:,LBj:,:,:)
165
166 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
167 real(r8), intent(in) :: v (LBi:,LBj:,:,:)
168 real(r8), intent(in) :: vG(LBi:,LBj:,:,:)
169# endif
170 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
171 real(r8), intent(in) :: ubar (LBi:,LBj:,:)
172 real(r8), intent(in) :: ubarG(LBi:,LBj:,:)
173
174 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
175 real(r8), intent(in) :: vbar (LBi:,LBj:,:)
176 real(r8), intent(in) :: vbarG(LBi:,LBj:,:)
177
178 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
179 real(r8), intent(in) :: zeta (LBi:,LBj:,:)
180 real(r8), intent(in) :: zetaG(LBi:,LBj:,:)
181# else
182# ifdef MASKING
183 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
184 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
185 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
186# endif
187# ifdef SOLVE3D
188 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
189 real(r8), intent(in) :: t (LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
190 real(r8), intent(in) :: tG(LBi:UBi,LBj:UBj,N(ng),2,NT(ng))
191
192 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
193 real(r8), intent(in) :: u (LBi:UBi,LBj:UBj,N(ng),2)
194 real(r8), intent(in) :: uG(LBi:UBi,LBj:UBj,N(ng),2)
195
196 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
197 real(r8), intent(in) :: v (LBi:UBi,LBj:UBj,N(ng),2)
198 real(r8), intent(in) :: vG(LBi:UBi,LBj:UBj,N(ng),2)
199# endif
200 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,:)
201 real(r8), intent(in) :: ubar (LBi:UBi,LBj:UBj,:)
202 real(r8), intent(in) :: ubarG(LBi:UBi,LBj:UBj,2)
203
204 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,:)
205 real(r8), intent(in) :: vbar (LBi:UBi,LBj:UBj,:)
206 real(r8), intent(in) :: vbarG(LBi:UBi,LBj:UBj,2)
207
208 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,:)
209 real(r8), intent(in) :: zeta (LBi:UBi,LBj:UBj,:)
210 real(r8), intent(in) :: zetaG(LBi:UBi,LBj:UBj,2)
211# endif
212!
213! Local variable declarations.
214!
215 integer :: NSUB, Tindex, i, j
216# ifdef SOLVE3D
217 integer :: itrc, k
218# endif
219 real(r8) :: my_dot, p, scale
220
221# ifdef DISTRIBUTE
222 character(len=3 ) :: op_handle
223# endif
224 character(len=12) :: text
225
226# include "set_bounds.h"
227!
228!-----------------------------------------------------------------------
229! Add a perturbation to nonlinear state variables: steepest descent
230! direction times the perturbation amplitude "p".
231!-----------------------------------------------------------------------
232!
233 IF ((outer.ge.1).and.(mod(iic(ng)-1,nhis(ng)).eq.0).and. &
234 & ((nrrec(ng).eq.0).or.(iic(ng).ne.ntstart(ng)))) THEN
235!
236! Set perturbation amplitude as a function of the inner loop.
237!
238 p=10.0_r8**real(-inner,r8)
239 scale=1.0_r8/sqrt(addotproduct)
240 my_dot=0.0_r8
241 IF (time(ng).eq.tintrp(1,idfsur,ng)) THEN
242 tindex=1
243 ELSE IF (time(ng).eq.tintrp(2,idfsur,ng)) THEN
244 tindex=2
245 END IF
246 IF (iic(ng).eq.ntstart(ng)) THEN
247 text='------------'
248 ELSE
249 text=' '
250 END IF
251!
252! Free-surface.
253!
254 DO j=jstrr,jendr
255 DO i=istrr,iendr
256 my_dot=my_dot+ &
257 & p*scale*ad_zeta(i,j,linp)* &
258 & (zeta(i,j,kinp)-zetag(i,j,tindex))
259# ifdef MASKING
260 my_dot=my_dot*rmask(i,j)
261# endif
262 END DO
263 END DO
264!
265! 2D u-momentum component.
266!
267 DO j=jstrr,jendr
268 DO i=istr,iendr
269 my_dot=my_dot+ &
270 & p*scale*ad_ubar(i,j,linp)* &
271 & (ubar(i,j,kinp)-ubarg(i,j,tindex))
272# ifdef MASKING
273 my_dot=my_dot*umask(i,j)
274# endif
275 END DO
276 END DO
277!
278! 2D v-momentum component.
279!
280 DO j=jstr,jendr
281 DO i=istrr,iendr
282 my_dot=my_dot+ &
283 & p*scale*ad_vbar(i,j,linp)* &
284 & (vbar(i,j,kinp)-vbarg(i,j,tindex))
285# ifdef MASKING
286 my_dot=my_dot*vmask(i,j)
287# endif
288 END DO
289 END DO
290# ifdef SOLVE3D
291!
292! 3D u-momentum component.
293!
294 DO k=1,n(ng)
295 DO j=jstrr,jendr
296 DO i=istr,iendr
297 my_dot=my_dot+ &
298 & p*scale*ad_u(i,j,k,linp)* &
299 & (u(i,j,k,ninp)-ug(i,j,k,tindex))
300# ifdef MASKING
301 my_dot=my_dot*umask(i,j)
302# endif
303 END DO
304 END DO
305 END DO
306!
307! 3D v-momentum component.
308!
309 DO k=1,n(ng)
310 DO j=jstr,jendr
311 DO i=istrr,iendr
312 my_dot=my_dot+ &
313 & p*scale*ad_v(i,j,k,linp)* &
314 & (v(i,j,k,ninp)-vg(i,j,k,tindex))
315# ifdef MASKING
316 my_dot=my_dot*vmask(i,j)
317# endif
318 END DO
319 END DO
320 END DO
321!
322! Tracers.
323!
324 DO itrc=1,nt(ng)
325 DO k=1,n(ng)
326 DO j=jstrr,jendr
327 DO i=istrr,iendr
328 my_dot=my_dot+ &
329 & p*scale*ad_t(i,j,k,linp,itrc)* &
330 & (t(i,j,k,ninp,itrc)-tg(i,j,k,tindex,itrc))
331# ifdef MASKING
332 my_dot=my_dot*rmask(i,j)
333# endif
334 END DO
335 END DO
336 END DO
337 END DO
338# endif
339!
340! Perform parallel global reduction operation: dot product.
341!
342# ifdef DISTRIBUTE
343 nsub=1 ! distributed-memory
344# else
345 IF (domain(ng)%SouthWest_Corner(tile).and. &
346 & domain(ng)%NorthEast_Corner(tile)) THEN
347 nsub=1 ! non-tiled application
348 ELSE
349 nsub=ntilex(ng)*ntilee(ng) ! tiled application
350 END IF
351# endif
352!$OMP CRITICAL (NL_DOT)
353 IF (tile_count.eq.0) THEN
354 dotproduct=my_dot
355 ELSE
356 dotproduct=dotproduct+my_dot
357 END IF
358 tile_count=tile_count+1
359 IF (tile_count.eq.nsub) THEN
360 tile_count=0
361# ifdef DISTRIBUTE
362 op_handle='SUM'
363 CALL mp_reduce (ng, inlm, 1, dotproduct, op_handle)
364# endif
365 ig1count=ig1count+1
366 g1(ig1count)=dotproduct
367 IF (master) THEN
368 WRITE (stdout,10) outer, inner, p, text, tindex, dotproduct
369 10 FORMAT (5x,'(Outer,Inner) = ','(',i4.4,',',i4.4,')',3x, &
370 & 'TLM Check, p = ',1p,e21.14,0p,/,5x,a,4x, &
371 & '(Index=',i1,')',5x,'TLM Check, g1 = ',1p,e21.14,0p)
372 END IF
373 END IF
374!$OMP END CRITICAL (NL_DOT)
375 END IF
376!
377 RETURN
378 END SUBROUTINE nl_dotproduct_tile
379# endif
380
381# if defined ADJOINT && defined PROPAGATOR
382 SUBROUTINE ad_statenorm (ng, tile, Kinp, Ninp, StateNorm)
383!
384!=======================================================================
385! !
386! This routine computes the dot product norm of requested adjoint !
387! state. !
388! !
389!=======================================================================
390!
391 USE mod_param
392 USE mod_grid
393 USE mod_ocean
394!
395! Imported variable declarations.
396!
397 integer, intent(in) :: ng, tile, kinp, ninp
398
399 real(r8), intent(out) :: statenorm
400!
401! Local variable declarations.
402!
403 character (len=*), parameter :: myfile = &
404 & __FILE__//", ad_statenorm"
405!
406# include "tile.h"
407!
408# ifdef PROFILE
409 CALL wclock_on (ng, iadm, 7, __line__, myfile)
410# endif
411 CALL ad_statenorm_tile (ng, tile, &
412 & lbi, ubi, lbj, ubj, &
413 & imins, imaxs, jmins, jmaxs, &
414 & kinp, ninp, &
415 & statenorm, &
416# ifdef MASKING
417 & grid(ng) % rmask, &
418 & grid(ng) % umask, &
419 & grid(ng) % vmask, &
420# endif
421# ifdef SOLVE3D
422 & grid(ng) % Hz, &
423# else
424 & grid(ng) % h, &
425# endif
426# ifdef SOLVE3D
427 & ocean(ng) % ad_t, &
428 & ocean(ng) % ad_u, &
429 & ocean(ng) % ad_v, &
430# else
431 & ocean(ng) % ad_ubar, &
432 & ocean(ng) % ad_vbar, &
433# endif
434 & ocean(ng) % ad_zeta)
435# ifdef PROFILE
436 CALL wclock_off (ng, iadm, 7, __line__, myfile)
437# endif
438!
439 RETURN
440 END SUBROUTINE ad_statenorm
441!
442!***********************************************************************
443 SUBROUTINE ad_statenorm_tile (ng, tile, &
444 & LBi, UBi, LBj, UBj, &
445 & IminS, ImaxS, JminS, JmaxS, &
446 & Kinp, Ninp, &
447 & StateNorm, &
448# ifdef MASKING
449 & rmask, umask, vmask, &
450# endif
451# ifdef SOLVE3D
452 & Hz, &
453# else
454 & h, &
455# endif
456# ifdef SOLVE3D
457 & ad_t, ad_u, ad_v, &
458# else
459 & ad_ubar, ad_vbar, &
460# endif
461 & ad_zeta)
462!***********************************************************************
463!
464 USE mod_param
465 USE mod_parallel
466 USE mod_scalars
467
468# ifdef DISTRIBUTE
469!
470 USE distribute_mod, ONLY : mp_reduce
471# endif
472!
473! Imported variable declarations.
474!
475 integer, intent(in) :: ng, tile
476 integer, intent(in) :: LBi, UBi, LBj, UBj
477 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
478 integer, intent(in) :: Kinp, Ninp
479
480 real(r8), intent(out) :: StateNorm
481!
482# ifdef ASSUMED_SHAPE
483# ifdef MASKING
484 real(r8), intent(in) :: rmask(LBi:,LBj:)
485 real(r8), intent(in) :: umask(LBi:,LBj:)
486 real(r8), intent(in) :: vmask(LBi:,LBj:)
487# endif
488# ifdef SOLVE3D
489 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
490# else
491 real(r8), intent(in) :: h(LBi:,LBj:)
492# endif
493# ifdef SOLVE3D
494 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
495 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
496 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
497# else
498 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
499 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
500# endif
501 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
502# else
503# ifdef MASKING
504 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
505 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
506 real(r8), intent(in) :: vmask(LBi:Ubi,LBj:UBj)
507# endif
508# ifdef SOLVE3D
509 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
510# else
511 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
512# endif
513# ifdef SOLVE3D
514 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
515 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
516 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
517# else
518 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,:)
519 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,:)
520# endif
521 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,:)
522# endif
523!
524! Local variable declarations.
525!
526 integer :: NSUB, i, j
527# ifdef SOLVE3D
528 integer :: itrc, k
529# endif
530 real(r8) :: cff, cff1, my_norm, scale
531
532# ifdef DISTRIBUTE
533 character(len=3) :: op_handle
534# endif
535
536# include "set_bounds.h"
537!
538!-----------------------------------------------------------------------
539! Compute dot product norm of current adjoint solution.
540!-----------------------------------------------------------------------
541
542# ifdef FULL_GRID
543# define IR_RANGE IstrT,IendT
544# define IU_RANGE IstrP,IendT
545# define JR_RANGE JstrT,JendT
546# define JV_RANGE JstrP,JendT
547# else
548# define IR_RANGE Istr,Iend
549# define IU_RANGE IstrU,Iend
550# define JR_RANGE Jstr,Jend
551# define JV_RANGE JstrV,Jend
552# endif
553!
554! Initialize private norm.
555!
556 my_norm=0.0_r8
557!
558! Free-surface.
559!
560 scale=0.5_r8*g*rho0
561 DO j=jr_range
562 DO i=ir_range
563 cff1=scale*ad_zeta(i,j,kinp)*ad_zeta(i,j,kinp)
564# ifdef MASKING
565 cff1=cff1*rmask(i,j)
566# endif
567 my_norm=my_norm+cff1
568 END DO
569 END DO
570
571# ifndef SOLVE3D
572!
573! 2D u-momentum component.
574!
575 cff=0.25_r8*rho0
576 DO j=jr_range
577 DO i=iu_range
578 scale=cff*(h(i-1,j)+h(i,j))
579 cff1=scale*ad_ubar(i,j,kinp)*ad_ubar(i,j,kinp)
580# ifdef MASKING
581 cff1=cff1*umask(i,j)
582# endif
583 my_norm=my_norm+cff1
584 END DO
585 END DO
586!
587! 2D v-momentum component.
588!
589 cff=0.25_r8*rho0
590 DO j=jv_range
591 DO i=ir_range
592 scale=cff*(h(i,j-1)+h(i,j))
593 cff1=scale*ad_vbar(i,j,kinp)*ad_vbar(i,j,kinp)
594# ifdef MASKING
595 cff1=cff1*vmask(i,j)
596# endif
597 my_norm=my_norm+cff1
598 END DO
599 END DO
600# else
601!
602! 3D u-momentum component.
603!
604 cff=0.25_r8*rho0
605 DO k=1,n(ng)
606 DO j=jr_range
607 DO i=iu_range
608 scale=cff*(hz(i-1,j,k)+hz(i,j,k))
609 cff1=scale*ad_u(i,j,k,ninp)*ad_u(i,j,k,ninp)
610# ifdef MASKING
611 cff1=cff1*umask(i,j)
612# endif
613 my_norm=my_norm+cff1
614 END DO
615 END DO
616 END DO
617!
618! 3D v-momentum component.
619!
620 cff=0.25_r8*rho0
621 DO k=1,n(ng)
622 DO j=jv_range
623 DO i=ir_range
624 scale=cff*(hz(i,j-1,k)+hz(i,j,k))
625 cff1=scale*ad_v(i,j,k,ninp)*ad_v(i,j,k,ninp)
626# ifdef MASKING
627 cff1=cff1*vmask(i,j)
628# endif
629 my_norm=my_norm+cff1
630 END DO
631 END DO
632 END DO
633!
634! Tracers. For now, use salinity scale for passive tracers.
635!
636 DO itrc=1,nt(ng)
637 IF (itrc.eq.itemp) THEN
638 cff=0.5_r8*rho0*tcoef(ng)*tcoef(ng)*g*g/bvf_bak
639 ELSE IF (itrc.eq.isalt) THEN
640 cff=0.5_r8*rho0*scoef(ng)*scoef(ng)*g*g/bvf_bak
641 ELSE
642 cff=0.5_r8*rho0*scoef(ng)*scoef(ng)*g*g/bvf_bak
643 END IF
644 DO k=1,n(ng)
645 DO j=jr_range
646 DO i=ir_range
647 scale=cff*hz(i,j,k)
648 cff1=scale*ad_t(i,j,k,ninp,itrc)* &
649 & ad_t(i,j,k,ninp,itrc)
650# ifdef MASKING
651 cff1=cff1*rmask(i,j)
652# endif
653 my_norm=my_norm+cff1
654 END DO
655 END DO
656 END DO
657 END DO
658# endif
659!
660! Perform parallel global reduction operation: dot product.
661!
662# ifdef DISTRIBUTE
663 nsub=1 ! distributed-memory
664# else
665 IF (domain(ng)%SouthWest_Corner(tile).and. &
666 & domain(ng)%NorthEast_Corner(tile)) THEN
667 nsub=1 ! non-tiled application
668 ELSE
669 nsub=ntilex(ng)*ntilee(ng) ! tiled application
670 END IF
671# endif
672!$OMP CRITICAL (TL_NORM)
673 IF (tile_count.eq.0) THEN
674 statenorm=my_norm
675 ELSE
676 statenorm=statenorm+my_norm
677 END IF
678 tile_count=tile_count+1
679 IF (tile_count.eq.nsub) THEN
680 tile_count=0
681# ifdef DISTRIBUTE
682 op_handle='SUM'
683 CALL mp_reduce (ng, itlm, 1, statenorm, op_handle)
684# endif
685 END IF
686!$OMP END CRITICAL (TL_NORM)
687
688# undef IR_RANGE
689# undef IU_RANGE
690# undef JR_RANGE
691# undef JV_RANGE
692!
693 RETURN
694 END SUBROUTINE ad_statenorm_tile
695# endif
696
697# if defined TANGENT && defined PROPAGATOR
698 SUBROUTINE tl_statenorm (ng, tile, Kinp, Ninp, StateNorm)
699!
700!=======================================================================
701! !
702! This routine computes the dot product norm of requested tangent !
703! linear state. !
704! !
705!=======================================================================
706!
707 USE mod_param
708 USE mod_grid
709 USE mod_ocean
710!
711! Imported variable declarations.
712!
713 integer, intent(in) :: ng, tile, kinp, ninp
714
715 real(r8), intent(out) :: statenorm
716!
717! Local variable declarations.
718!
719 character (len=*), parameter :: myfile = &
720 & __FILE__//", tl_statenorm"
721!
722# include "tile.h"
723!
724# ifdef PROFILE
725 CALL wclock_on (ng, itlm, 7, __line__, myfile)
726# endif
727 CALL tl_statenorm_tile (ng, tile, &
728 & lbi, ubi, lbj, ubj, &
729 & imins, imaxs, jmins, jmaxs, &
730 & kinp, ninp, &
731 & statenorm, &
732# ifdef MASKING
733 & grid(ng) % rmask, &
734 & grid(ng) % umask, &
735 & grid(ng) % vmask, &
736# endif
737# ifdef BNORM
738# ifdef SOLVE3D
739 & ocean(ng) % e_t, &
740 & ocean(ng) % e_u, &
741 & ocean(ng) % e_v, &
742# else
743 & ocean(ng) % e_ubar, &
744 & ocean(ng) % e_vbar, &
745# endif
746 & ocean(ng) % e_zeta, &
747# endif
748# ifdef SOLVE3D
749 & grid(ng) % Hz, &
750# else
751 & grid(ng) % h, &
752# endif
753# ifdef SOLVE3D
754 & ocean(ng) % tl_t, &
755 & ocean(ng) % tl_u, &
756 & ocean(ng) % tl_v, &
757# else
758 & ocean(ng) % tl_ubar, &
759 & ocean(ng) % tl_vbar, &
760# endif
761 & ocean(ng) % tl_zeta)
762# ifdef PROFILE
763 CALL wclock_off (ng, itlm, 7, __line__, myfile)
764# endif
765!
766 RETURN
767 END SUBROUTINE tl_statenorm
768
769!
770!***********************************************************************
771 SUBROUTINE tl_statenorm_tile (ng, tile, &
772 & LBi, UBi, LBj, UBj, &
773 & IminS, ImaxS, JminS, JmaxS, &
774 & Kinp, Ninp, &
775 & StateNorm, &
776# ifdef MASKING
777 & rmask, umask, vmask, &
778# endif
779# ifdef BNORM
780# ifdef SOLVE3D
781 & e_t, e_u, e_v, &
782# else
783 & e_ubar, e_vbar, &
784# endif
785 & e_zeta, &
786# endif
787# ifdef SOLVE3D
788 & Hz, &
789# else
790 & h, &
791# endif
792# ifdef SOLVE3D
793 & tl_t, tl_u, tl_v, &
794# else
795 & tl_ubar, tl_vbar, &
796# endif
797 & tl_zeta)
798!***********************************************************************
799!
800 USE mod_param
801 USE mod_parallel
802 USE mod_scalars
803
804# ifdef DISTRIBUTE
805!
806 USE distribute_mod, ONLY : mp_reduce
807# endif
808!
809! Imported variable declarations.
810!
811 integer, intent(in) :: ng, tile
812 integer, intent(in) :: LBi, UBi, LBj, UBj
813 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
814 integer, intent(in) :: Kinp, Ninp
815
816 real(r8), intent(out) :: StateNorm
817!
818# ifdef ASSUMED_SHAPE
819# ifdef MASKING
820 real(r8), intent(in) :: rmask(LBi:,LBj:)
821 real(r8), intent(in) :: umask(LBi:,LBj:)
822 real(r8), intent(in) :: vmask(LBi:,LBj:)
823# endif
824# ifdef SOLVE3D
825 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
826# else
827 real(r8), intent(in) :: h(LBi:,LBj:)
828# endif
829# ifdef SOLVE3D
830 real(r8), intent(in) :: tl_t(LBi:,LBj:,:,:,:)
831 real(r8), intent(in) :: tl_u(LBi:,LBj:,:,:)
832 real(r8), intent(in) :: tl_v(LBi:,LBj:,:,:)
833# else
834 real(r8), intent(in) :: tl_ubar(LBi:,LBj:,:)
835 real(r8), intent(in) :: tl_vbar(LBi:,LBj:,:)
836# endif
837 real(r8), intent(in) :: tl_zeta(LBi:,LBj:,:)
838# ifdef BNORM
839 real(r8), intent(in) :: e_zeta(LBi:,LBj:,:)
840# ifdef SOLVE3D
841 real(r8), intent(in) :: e_t(LBi:,LBj:,:,:,:)
842 real(r8), intent(in) :: e_u(LBi:,LBj:,:,:)
843 real(r8), intent(in) :: e_v(LBi:,LBj:,:,:)
844# else
845 real(r8), intent(in) :: e_ubar(LBi:,LBj:,:)
846 real(r8), intent(in) :: e_vbar(LBi:,LBj:,:)
847# endif
848# endif
849# else
850# ifdef MASKING
851 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
852 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
853 real(r8), intent(in) :: vmask(LBi:Ubi,LBj:UBj)
854# endif
855# ifdef SOLVE3D
856 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
857# else
858 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
859# endif
860# ifdef SOLVE3D
861 real(r8), intent(in) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
862 real(r8), intent(in) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
863 real(r8), intent(in) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
864# else
865 real(r8), intent(in) :: tl_ubar(LBi:UBi,LBj:UBj,:)
866 real(r8), intent(in) :: tl_vbar(LBi:UBi,LBj:UBj,:)
867# endif
868 real(r8), intent(in) :: tl_zeta(LBi:UBi,LBj:UBj,:)
869# ifdef BNORM
870 real(r8), intent(in) :: e_zeta(LBi:UBi,LBj:UBj,NSA)
871# ifdef SOLVE3D
872 real(r8), intent(in) :: e_t(LBi:UBi,LBj:UBj,N(ng),NSA,NT(ng))
873 real(r8), intent(in) :: e_u(LBi:UBi,LBj:UBj,N(ng),NSA)
874 real(r8), intent(in) :: e_v(LBi:UBi,LBj:UBj,N(ng),NSA)
875# else
876 real(r8), intent(in) :: e_ubar(LBi:UBi,LBj:UBj,NSA)
877 real(r8), intent(in) :: e_vbar(LBi:UBi,LBj:UBj,NSA)
878# endif
879# endif
880# endif
881!
882! Local variable declarations.
883!
884 integer :: NSUB, i, j
885# ifdef SOLVE3D
886 integer :: itrc, k
887# endif
888 real(r8) :: cff, cff1, my_norm, scale
889
890# ifdef DISTRIBUTE
891 character(len=3) :: op_handle
892# endif
893
894# include "set_bounds.h"
895!
896!-----------------------------------------------------------------------
897! Compute dot product norm of current tangent linear solution.
898!-----------------------------------------------------------------------
899
900# ifdef FULL_GRID
901# define IR_RANGE IstrT,IendT
902# define IU_RANGE IstrP,IendT
903# define JR_RANGE JstrT,JendT
904# define JV_RANGE JstrP,JendT
905# else
906# define IR_RANGE Istr,Iend
907# define IU_RANGE IstrU,Iend
908# define JR_RANGE Jstr,Jend
909# define JV_RANGE JstrV,Jend
910# endif
911!
912! Initialize private norm.
913!
914 my_norm=0.0_r8
915!
916! Free-surface.
917!
918 scale=0.5_r8*g*rho0
919 DO j=jr_range
920 DO i=ir_range
921# ifdef BNORM
922 IF (e_zeta(i,j,1).gt.0.0_r8) THEN
923 cff1=1.0_r8/(e_zeta(i,j,1)*e_zeta(i,j,1))
924 ELSE
925 cff1=1.0_r8
926 END IF
927 cff1=cff1*tl_zeta(i,j,kinp)*tl_zeta(i,j,kinp)
928# else
929 cff1=scale*tl_zeta(i,j,kinp)*tl_zeta(i,j,kinp)
930# endif
931# ifdef MASKING
932 cff1=cff1*rmask(i,j)
933# endif
934 my_norm=my_norm+cff1
935 END DO
936 END DO
937
938# ifndef SOLVE3D
939!
940! 2D u-momentum component.
941!
942 cff=0.25_r8*rho0
943 DO j=jr_range
944 DO i=iu_range
945# ifdef BNORM
946 IF (e_ubar(i,j,1).gt.0.0_r8) THEN
947 cff1=1.0_r8/(e_ubar(i,j,1)*e_ubar(i,j,1))
948 ELSE
949 cff1=1.0_r8
950 END IF
951 cff1=cff1*tl_ubar(i,j,kinp)*tl_ubar(i,j,kinp)
952# else
953 scale=cff*(h(i-1,j)+h(i,j))
954 cff1=scale*tl_ubar(i,j,kinp)*tl_ubar(i,j,kinp)
955# endif
956# ifdef MASKING
957 cff1=cff1*umask(i,j)
958# endif
959 my_norm=my_norm+cff1
960 END DO
961 END DO
962!
963! 2D v-momentum component.
964!
965 cff=0.25_r8*rho0
966 DO j=jv_range
967 DO i=ir_range
968# ifdef BNORM
969 IF (e_vbar(i,j,1).gt.0.0_r8) THEN
970 cff1=1.0_r8/(e_vbar(i,j,1)*e_vbar(i,j,1))
971 ELSE
972 cff1=1.0_r8
973 END IF
974 cff1=cff1*tl_vbar(i,j,kinp)*tl_vbar(i,j,kinp)
975# else
976 scale=cff*(h(i,j-1)+h(i,j))
977 cff1=scale*tl_vbar(i,j,kinp)*tl_vbar(i,j,kinp)
978# endif
979# ifdef MASKING
980 cff1=cff1*vmask(i,j)
981# endif
982 my_norm=my_norm+cff1
983 END DO
984 END DO
985# else
986!
987! 3D u-momentum component.
988!
989 cff=0.25_r8*rho0
990 DO k=1,n(ng)
991 DO j=jr_range
992 DO i=iu_range
993# ifdef BNORM
994 IF (e_u(i,j,k,1).gt.0.0_r8) THEN
995 cff1=1.0_r8/(e_u(i,j,k,1)*e_u(i,j,k,1))
996 ELSE
997 cff1=1.0_r8
998 END IF
999 cff1=cff1*tl_u(i,j,k,ninp)*tl_u(i,j,k,ninp)
1000# else
1001 scale=cff*(hz(i-1,j,k)+hz(i,j,k))
1002 cff1=scale*tl_u(i,j,k,ninp)*tl_u(i,j,k,ninp)
1003# endif
1004# ifdef MASKING
1005 cff1=cff1*umask(i,j)
1006# endif
1007 my_norm=my_norm+cff1
1008 END DO
1009 END DO
1010 END DO
1011!
1012! 3D v-momentum component.
1013!
1014 cff=0.25_r8*rho0
1015 DO k=1,n(ng)
1016 DO j=jv_range
1017 DO i=ir_range
1018# ifdef BNORM
1019 IF (e_v(i,j,k,1).gt.0.0_r8) THEN
1020 cff1=1.0_r8/(e_v(i,j,k,1)*e_v(i,j,k,1))
1021 ELSE
1022 cff1=1.0_r8
1023 END IF
1024 cff1=cff1*tl_v(i,j,k,ninp)*tl_v(i,j,k,ninp)
1025# else
1026 scale=cff*(hz(i,j-1,k)+hz(i,j,k))
1027 cff1=scale*tl_v(i,j,k,ninp)*tl_v(i,j,k,ninp)
1028# endif
1029# ifdef MASKING
1030 cff1=cff1*vmask(i,j)
1031# endif
1032 my_norm=my_norm+cff1
1033 END DO
1034 END DO
1035 END DO
1036!
1037! Tracers. For now, use salinity scale for passive tracers.
1038!
1039 DO itrc=1,nt(ng)
1040 IF (itrc.eq.itemp) THEN
1041 cff=0.5_r8*rho0*tcoef(ng)*tcoef(ng)*g*g/bvf_bak
1042 ELSE IF (itrc.eq.isalt) THEN
1043 cff=0.5_r8*rho0*scoef(ng)*scoef(ng)*g*g/bvf_bak
1044 ELSE
1045 cff=0.5_r8*rho0*scoef(ng)*scoef(ng)*g*g/bvf_bak
1046 END IF
1047 DO k=1,n(ng)
1048 DO j=jr_range
1049 DO i=ir_range
1050# ifdef BNORM
1051 IF (e_t(i,j,k,1,itrc).gt.0.0_r8) THEN
1052 cff1=1.0_r8/(e_t(i,j,k,1,itrc)*e_t(i,j,k,1,itrc))
1053 ELSE
1054 cff1=1.0_r8
1055 END IF
1056 cff1=cff1*tl_t(i,j,k,ninp,itrc)* &
1057 & tl_t(i,j,k,ninp,itrc)
1058# else
1059 scale=cff*hz(i,j,k)
1060 cff1=scale*tl_t(i,j,k,ninp,itrc)* &
1061 & tl_t(i,j,k,ninp,itrc)
1062# endif
1063# ifdef MASKING
1064 cff1=cff1*rmask(i,j)
1065# endif
1066 my_norm=my_norm+cff1
1067 END DO
1068 END DO
1069 END DO
1070 END DO
1071# endif
1072!
1073! Perform parallel global reduction operation: dot product.
1074!
1075# ifdef DISTRIBUTE
1076 nsub=1 ! distributed-memory
1077# else
1078 IF (domain(ng)%SouthWest_Corner(tile).and. &
1079 & domain(ng)%NorthEast_Corner(tile)) THEN
1080 nsub=1 ! non-tiled application
1081 ELSE
1082 nsub=ntilex(ng)*ntilee(ng) ! tiled application
1083 END IF
1084# endif
1085!$OMP CRITICAL (TL_NORM)
1086 IF (tile_count.eq.0) THEN
1087 statenorm=my_norm
1088 ELSE
1089 statenorm=statenorm+my_norm
1090 END IF
1091 tile_count=tile_count+1
1092 IF (tile_count.eq.nsub) THEN
1093 tile_count=0
1094# ifdef DISTRIBUTE
1095 op_handle='SUM'
1096 CALL mp_reduce (ng, itlm, 1, statenorm, op_handle)
1097# endif
1098 END IF
1099!$OMP END CRITICAL (TL_NORM)
1100
1101# undef IR_RANGE
1102# undef IU_RANGE
1103# undef JR_RANGE
1104# undef JV_RANGE
1105!
1106 RETURN
1107 END SUBROUTINE tl_statenorm_tile
1108# endif
1109
1110# ifdef TLM_CHECK
1111 SUBROUTINE tl_dotproduct (ng, tile, Linp)
1112!
1113!=======================================================================
1114! !
1115! This routine computes dot product function (g1) associated with !
1116! the tangent linear solution. !
1117! !
1118!=======================================================================
1119!
1120 USE mod_param
1121 USE mod_grid
1122 USE mod_ocean
1123 USE mod_stepping
1124!
1125! Imported variable declarations.
1126!
1127 integer, intent(in) :: ng, tile, linp
1128!
1129! Local variable declarations.
1130!
1131 character (len=*), parameter :: myfile = &
1132 & __FILE__//", tl_dotproduct"
1133!
1134# include "tile.h"
1135!
1136# ifdef PROFILE
1137 CALL wclock_on (ng, itlm, 7, __line__, myfile)
1138# endif
1139 CALL tl_dotproduct_tile (ng, tile, &
1140 & lbi, ubi, lbj, ubj, &
1141 & imins, imaxs, jmins, jmaxs, &
1142 & kout, linp, &
1143# ifdef SOLVE3D
1144 & nout, &
1145# endif
1146# ifdef MASKING
1147 & grid(ng) % rmask, &
1148 & grid(ng) % umask, &
1149 & grid(ng) % vmask, &
1150# endif
1151# ifdef SOLVE3D
1152 & ocean(ng) % ad_t, &
1153 & ocean(ng) % tl_t, &
1154 & ocean(ng) % ad_u, &
1155 & ocean(ng) % tl_u, &
1156 & ocean(ng) % ad_v, &
1157 & ocean(ng) % tl_v, &
1158# endif
1159 & ocean(ng) % ad_ubar, &
1160 & ocean(ng) % tl_ubar, &
1161 & ocean(ng) % ad_vbar, &
1162 & ocean(ng) % tl_vbar, &
1163 & ocean(ng) % ad_zeta, &
1164 & ocean(ng) % tl_zeta)
1165# ifdef PROFILE
1166 CALL wclock_off (ng, itlm, 7, __line__, myfile)
1167# endif
1168!
1169 RETURN
1170 END SUBROUTINE tl_dotproduct
1171
1172!
1173!***********************************************************************
1174 SUBROUTINE tl_dotproduct_tile (ng, tile, &
1175 & LBi, UBi, LBj, UBj, &
1176 & IminS, ImaxS, JminS, JmaxS, &
1177 & Kinp, Linp, &
1178# ifdef SOLVE3D
1179 & Ninp, &
1180# endif
1181# ifdef MASKING
1182 & rmask, umask, vmask, &
1183# endif
1184# ifdef SOLVE3D
1185 & ad_t, tl_t, &
1186 & ad_u, tl_u, &
1187 & ad_v, tl_v, &
1188# endif
1189 & ad_ubar, tl_ubar, &
1190 & ad_vbar, tl_vbar, &
1191 & ad_zeta, tl_zeta)
1192!***********************************************************************
1193!
1194 USE mod_param
1195 USE mod_parallel
1196 USE mod_fourdvar
1197 USE mod_iounits
1198 USE mod_scalars
1199
1200# ifdef DISTRIBUTE
1201!
1202 USE distribute_mod, ONLY : mp_reduce
1203# endif
1204!
1205! Imported variable declarations.
1206!
1207 integer, intent(in) :: ng, tile
1208 integer, intent(in) :: LBi, UBi, LBj, UBj
1209 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1210 integer, intent(in) :: Kinp, Linp
1211# ifdef SOLVE3D
1212 integer, intent(in) :: Ninp
1213# endif
1214!
1215# ifdef ASSUMED_SHAPE
1216# ifdef MASKING
1217 real(r8), intent(in) :: rmask(LBi:,LBj:)
1218 real(r8), intent(in) :: umask(LBi:,LBj:)
1219 real(r8), intent(in) :: vmask(LBi:,LBj:)
1220# endif
1221# ifdef SOLVE3D
1222 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
1223 real(r8), intent(in) :: tl_t(LBi:,LBj:,:,:,:)
1224
1225 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
1226 real(r8), intent(in) :: tl_u(LBi:,LBj:,:,:)
1227
1228 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
1229 real(r8), intent(in) :: tl_v(LBi:,LBj:,:,:)
1230# endif
1231 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
1232 real(r8), intent(in) :: tl_ubar(LBi:,LBj:,:)
1233
1234 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
1235 real(r8), intent(in) :: tl_vbar(LBi:,LBj:,:)
1236
1237 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
1238 real(r8), intent(in) :: tl_zeta(LBi:,LBj:,:)
1239# else
1240# ifdef MASKING
1241 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1242 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1243 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1244# endif
1245# ifdef SOLVE3D
1246 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1247 real(r8), intent(in) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1248
1249 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1250 real(r8), intent(in) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
1251
1252 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1253 real(r8), intent(in) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
1254# endif
1255 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,:)
1256 real(r8), intent(in) :: tl_ubar(LBi:UBi,LBj:UBj,:)
1257
1258 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,:)
1259 real(r8), intent(in) :: tl_vbar(LBi:UBi,LBj:UBj,:)
1260
1261 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,:)
1262 real(r8), intent(in) :: tl_zeta(LBi:UBi,LBj:UBj,:)
1263# endif
1264!
1265! Local variable declarations.
1266!
1267 integer :: NSUB, Tindex, i, j
1268# ifdef SOLVE3D
1269 integer :: itrc, k
1270# endif
1271 real(r8) :: my_dot, p, scale
1272
1273# ifdef DISTRIBUTE
1274 character(len=3 ) :: op_handle
1275# endif
1276 character(len=12) :: text
1277
1278# include "set_bounds.h"
1279!
1280!-----------------------------------------------------------------------
1281! Add a perturbation to nonlinear state variables: steepest descent
1282! direction times the perturbation amplitude "p".
1283!-----------------------------------------------------------------------
1284!
1285 IF ((outer.ge.1).and.(mod(iic(ng)-1,ntlm(ng)).eq.0).and. &
1286 & ((nrrec(ng).eq.0).or.(iic(ng).ne.ntstart(ng)))) THEN
1287!
1288! Set perturbation amplitude as a function of the inner loop.
1289!
1290 p=10.0_r8**real(-inner,r8)
1291 my_dot=0.0_r8
1292 scale=1.0_r8/sqrt(addotproduct)
1293 IF (iic(ng).eq.ntstart(ng)) THEN
1294 text='------------'
1295 ELSE
1296 text=' '
1297 END IF
1298!
1299! Free-surface.
1300!
1301 DO j=jstrr,jendr
1302 DO i=istrr,iendr
1303 my_dot=my_dot+ &
1304 & p*scale*ad_zeta(i,j,linp)*tl_zeta(i,j,kinp)
1305# ifdef MASKING
1306 my_dot=my_dot*rmask(i,j)
1307# endif
1308 END DO
1309 END DO
1310!
1311! 2D u-momentum component.
1312!
1313 DO j=jstrr,jendr
1314 DO i=istr,iendr
1315 my_dot=my_dot+ &
1316 & p*scale*ad_ubar(i,j,linp)*tl_ubar(i,j,kinp)
1317# ifdef MASKING
1318 my_dot=my_dot*umask(i,j)
1319# endif
1320 END DO
1321 END DO
1322!
1323! 2D v-momentum component.
1324!
1325 DO j=jstr,jendr
1326 DO i=istrr,iendr
1327 my_dot=my_dot+ &
1328 & p*scale*ad_vbar(i,j,linp)*tl_vbar(i,j,kinp)
1329# ifdef MASKING
1330 my_dot=my_dot*vmask(i,j)
1331# endif
1332 END DO
1333 END DO
1334# ifdef SOLVE3D
1335!
1336! 3D u-momentum component.
1337!
1338 DO k=1,n(ng)
1339 DO j=jstrr,jendr
1340 DO i=istr,iendr
1341 my_dot=my_dot+ &
1342 & p*scale*ad_u(i,j,k,linp)*tl_u(i,j,k,ninp)
1343# ifdef MASKING
1344 my_dot=my_dot*umask(i,j)
1345# endif
1346 END DO
1347 END DO
1348 END DO
1349!
1350! 3D v-momentum component.
1351!
1352 DO k=1,n(ng)
1353 DO j=jstr,jendr
1354 DO i=istrr,iendr
1355 my_dot=my_dot+ &
1356 & p*scale*ad_v(i,j,k,linp)*tl_v(i,j,k,ninp)
1357# ifdef MASKING
1358 my_dot=my_dot*vmask(i,j)
1359# endif
1360 END DO
1361 END DO
1362 END DO
1363!
1364! Tracers.
1365!
1366 DO itrc=1,nt(ng)
1367 DO k=1,n(ng)
1368 DO j=jstrr,jendr
1369 DO i=istrr,iendr
1370 my_dot=my_dot+ &
1371 & p*scale*ad_t(i,j,k,linp,itrc)* &
1372 & tl_t(i,j,k,ninp,itrc)
1373# ifdef MASKING
1374 my_dot=my_dot*rmask(i,j)
1375# endif
1376 END DO
1377 END DO
1378 END DO
1379 END DO
1380# endif
1381!
1382! Perform parallel global reduction operation: dot product.
1383!
1384# ifdef DISTRIBUTE
1385 nsub=1 ! distributed-memory
1386# else
1387 IF (domain(ng)%SouthWest_Corner(tile).and. &
1388 & domain(ng)%NorthEast_Corner(tile)) THEN
1389 nsub=1 ! non-tiled application
1390 ELSE
1391 nsub=ntilex(ng)*ntilee(ng) ! tiled application
1392 END IF
1393# endif
1394!$OMP CRITICAL (TL_DOT)
1395 IF (tile_count.eq.0) THEN
1396 dotproduct=my_dot
1397 ELSE
1398 dotproduct=dotproduct+my_dot
1399 END IF
1400 tile_count=tile_count+1
1401 IF (tile_count.eq.nsub) THEN
1402 tile_count=0
1403# ifdef DISTRIBUTE
1404 op_handle='SUM'
1405 CALL mp_reduce (ng, itlm, 1, dotproduct, op_handle)
1406# endif
1407 ig2count=ig2count+1
1408 g2(ig2count)=dotproduct
1409 IF (master) THEN
1410 WRITE (stdout,10) outer, inner, p, text, dotproduct
1411 10 FORMAT (5x,'(Outer,Inner) = ','(',i4.4,',',i4.4,')',3x, &
1412 & 'TLM Check, p = ',1p,e21.14,0p,/,5x,a,18x, &
1413 & 'TLM Check, g2 = ',1p,e21.14,0p)
1414 END IF
1415 END IF
1416!$OMP END CRITICAL (TL_DOT)
1417 END IF
1418!
1419 RETURN
1420 END SUBROUTINE tl_dotproduct_tile
1421# endif
1422
1423# ifdef TLM_CHECK
1424 SUBROUTINE ad_dotproduct (ng, tile, Linp)
1425!
1426!=======================================================================
1427! !
1428! This routine computes the adjoint solution dot product to scale !
1429! the adjoint solution. !
1430! !
1431!=======================================================================
1432!
1433 USE mod_param
1434 USE mod_grid
1435 USE mod_ocean
1436!
1437! Imported variable declarations.
1438!
1439 integer, intent(in) :: ng, tile, linp
1440!
1441! Local variable declarations.
1442!
1443 character (len=*), parameter :: myfile = &
1444 & __FILE__//", ad_dotproduct"
1445!
1446# include "tile.h"
1447!
1448# ifdef PROFILE
1449 CALL wclock_on (ng, iadm, 7, __line__, myfile)
1450# endif
1451 CALL ad_dotproduct_tile (ng, tile, &
1452 & lbi, ubi, lbj, ubj, &
1453 & imins, imaxs, jmins, jmaxs, &
1454 & linp, &
1455# ifdef MASKING
1456 & grid(ng) % rmask, &
1457 & grid(ng) % umask, &
1458 & grid(ng) % vmask, &
1459# endif
1460# ifdef SOLVE3D
1461 & ocean(ng) % ad_t, &
1462 & ocean(ng) % ad_u, &
1463 & ocean(ng) % ad_v, &
1464# endif
1465 & ocean(ng) % ad_ubar, &
1466 & ocean(ng) % ad_vbar, &
1467 & ocean(ng) % ad_zeta)
1468# ifdef PROFILE
1469 CALL wclock_off (ng, iadm, 7, __line__, myfile)
1470# endif
1471!
1472 RETURN
1473 END SUBROUTINE ad_dotproduct
1474!
1475!***********************************************************************
1476 SUBROUTINE ad_dotproduct_tile (ng, tile, &
1477 & LBi, UBi, LBj, UBj, &
1478 & IminS, ImaxS, JminS, JmaxS, &
1479 & Linp, &
1480# ifdef MASKING
1481 & rmask, umask, vmask, &
1482# endif
1483# ifdef SOLVE3D
1484 & ad_t, ad_u, ad_v, &
1485# endif
1486 & ad_ubar, ad_vbar, ad_zeta)
1487!***********************************************************************
1488!
1489 USE mod_param
1490 USE mod_parallel
1491 USE mod_fourdvar
1492 USE mod_iounits
1493 USE mod_scalars
1494
1495# ifdef DISTRIBUTE
1496!
1497 USE distribute_mod, ONLY : mp_reduce
1498# endif
1499!
1500! Imported variable declarations.
1501!
1502 integer, intent(in) :: ng, tile
1503 integer, intent(in) :: LBi, UBi, LBj, UBj
1504 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1505 integer, intent(in) :: Linp
1506!
1507# ifdef ASSUMED_SHAPE
1508# ifdef MASKING
1509 real(r8), intent(in) :: rmask(LBi:,LBj:)
1510 real(r8), intent(in) :: umask(LBi:,LBj:)
1511 real(r8), intent(in) :: vmask(LBi:,LBj:)
1512# endif
1513# ifdef SOLVE3D
1514 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
1515 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
1516 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
1517# endif
1518 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
1519 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
1520 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
1521# else
1522# ifdef MASKING
1523 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1524 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1525 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1526# endif
1527# ifdef SOLVE3D
1528 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1529 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1530 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1531# endif
1532 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,:)
1533 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,:)
1534 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,:)
1535# endif
1536!
1537! Local variable declarations.
1538!
1539 integer :: NSUB, i, j
1540# ifdef SOLVE3D
1541 integer :: itrc, k
1542# endif
1543 real(r8) :: my_dot
1544
1545# ifdef DISTRIBUTE
1546 character(len=3) :: op_handle
1547# endif
1548
1549# include "set_bounds.h"
1550!
1551!-----------------------------------------------------------------------
1552! Compute adjoint solution dot product.
1553!-----------------------------------------------------------------------
1554!
1555 my_dot=0.0_r8
1556!
1557! Free-surface.
1558!
1559 DO j=jstrr,jendr
1560 DO i=istrr,iendr
1561 my_dot=my_dot+ &
1562 & ad_zeta(i,j,linp)*ad_zeta(i,j,linp)
1563# ifdef MASKING
1564 my_dot=my_dot*rmask(i,j)
1565# endif
1566 END DO
1567 END DO
1568# ifndef SOLVE3D
1569!
1570! 2D u-momentum component.
1571!
1572 DO j=jstrr,jendr
1573 DO i=istr,iendr
1574 my_dot=my_dot+ &
1575 & ad_ubar(i,j,linp)*ad_ubar(i,j,linp)
1576# ifdef MASKING
1577 my_dot=my_dot*umask(i,j)
1578# endif
1579 END DO
1580 END DO
1581!
1582! 2D v-momentum component.
1583!
1584 DO j=jstr,jendr
1585 DO i=istrr,iendr
1586 my_dot=my_dot+ &
1587 & ad_vbar(i,j,linp)*ad_vbar(i,j,linp)
1588# ifdef MASKING
1589 my_dot=my_dot*vmask(i,j)
1590# endif
1591 END DO
1592 END DO
1593# endif
1594# ifdef SOLVE3D
1595!
1596! 3D u-momentum component.
1597!
1598 DO k=1,n(ng)
1599 DO j=jstrr,jendr
1600 DO i=istr,iendr
1601 my_dot=my_dot+ &
1602 & ad_u(i,j,k,linp)*ad_u(i,j,k,linp)
1603# ifdef MASKING
1604 my_dot=my_dot*umask(i,j)
1605# endif
1606 END DO
1607 END DO
1608 END DO
1609!
1610! 3D v-momentum component.
1611!
1612 DO k=1,n(ng)
1613 DO j=jstr,jendr
1614 DO i=istrr,iendr
1615 my_dot=my_dot+ &
1616 & ad_v(i,j,k,linp)*ad_v(i,j,k,linp)
1617# ifdef MASKING
1618 my_dot=my_dot*vmask(i,j)
1619# endif
1620 END DO
1621 END DO
1622 END DO
1623!
1624! Tracers.
1625!
1626 DO itrc=1,nt(ng)
1627 DO k=1,n(ng)
1628 DO j=jstrr,jendr
1629 DO i=istrr,iendr
1630 my_dot=my_dot+ &
1631 & ad_t(i,j,k,linp,itrc)*ad_t(i,j,k,linp,itrc)
1632# ifdef MASKING
1633 my_dot=my_dot*rmask(i,j)
1634# endif
1635 END DO
1636 END DO
1637 END DO
1638 END DO
1639# endif
1640!
1641! Perform parallel global reduction operation: dot product.
1642!
1643# ifdef DISTRIBUTE
1644 nsub=1 ! distributed-memory
1645# else
1646 IF (domain(ng)%SouthWest_Corner(tile).and. &
1647 & domain(ng)%NorthEast_Corner(tile)) THEN
1648 nsub=1 ! non-tiled application
1649 ELSE
1650 nsub=ntilex(ng)*ntilee(ng) ! tiled application
1651 END IF
1652# endif
1653!$OMP CRITICAL (AD_DOT)
1654 IF (tile_count.eq.0) THEN
1655 addotproduct=my_dot
1656 ELSE
1657 addotproduct=addotproduct+my_dot
1658 END IF
1659 tile_count=tile_count+1
1660 IF (tile_count.eq.nsub) THEN
1661 tile_count=0
1662# ifdef DISTRIBUTE
1663 op_handle='SUM'
1664 CALL mp_reduce (ng, iadm, 1, addotproduct, op_handle)
1665# endif
1666 IF (master) THEN
1667 WRITE (stdout,10) addotproduct
1668 10 FORMAT (/,' Adjoint Solution Dot Product: ',1p,e21.14)
1669 END IF
1670 END IF
1671!$OMP END CRITICAL (AD_DOT)
1672!
1673 RETURN
1674 END SUBROUTINE ad_dotproduct_tile
1675# endif
1676
1677#endif
1678 END MODULE dotproduct_mod
distribute_mod::mp_reduce
Definition distribute.F:113
distribute_mod
Definition distribute.F:3
dotproduct_mod
Definition dotproduct.F:2
dotproduct_mod::tl_statenorm_tile
subroutine tl_statenorm_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, ninp, statenorm, rmask, umask, vmask, e_t, e_u, e_v, e_zeta, hz, tl_t, tl_u, tl_v, tl_zeta)
Definition dotproduct.F:798
dotproduct_mod::tl_dotproduct
subroutine, public tl_dotproduct(ng, tile, linp)
Definition dotproduct.F:1112
dotproduct_mod::ad_dotproduct_tile
subroutine ad_dotproduct_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, linp, rmask, umask, vmask, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta)
Definition dotproduct.F:1487
dotproduct_mod::tl_statenorm
subroutine, public tl_statenorm(ng, tile, kinp, ninp, statenorm)
Definition dotproduct.F:699
dotproduct_mod::nl_dotproduct_tile
subroutine nl_dotproduct_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, linp, ninp, rmask, umask, vmask, ad_t, t, tg, ad_u, u, ug, ad_v, v, vg, ad_ubar, ubar, ubarg, ad_vbar, vbar, vbarg, ad_zeta, zeta, zetag)
Definition dotproduct.F:127
dotproduct_mod::ad_dotproduct
subroutine, public ad_dotproduct(ng, tile, linp)
Definition dotproduct.F:1425
dotproduct_mod::ad_statenorm
subroutine, public ad_statenorm(ng, tile, kinp, ninp, statenorm)
Definition dotproduct.F:383
dotproduct_mod::ad_statenorm_tile
subroutine ad_statenorm_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, ninp, statenorm, rmask, umask, vmask, hz, ad_t, ad_u, ad_v, ad_zeta)
Definition dotproduct.F:462
dotproduct_mod::tl_dotproduct_tile
subroutine tl_dotproduct_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, linp, ninp, rmask, umask, vmask, ad_t, tl_t, ad_u, tl_u, ad_v, tl_v, ad_ubar, tl_ubar, ad_vbar, tl_vbar, ad_zeta, tl_zeta)
Definition dotproduct.F:1192
dotproduct_mod::nl_dotproduct
subroutine, public nl_dotproduct(ng, tile, linp)
Definition dotproduct.F:40
mod_fourdvar
Definition mod_fourdvar.F:2
mod_fourdvar::ig2count
integer ig2count
Definition mod_fourdvar.F:874
mod_fourdvar::g2
real(r8), dimension(1000) g2
Definition mod_fourdvar.F:877
mod_fourdvar::g1
real(r8), dimension(1000) g1
Definition mod_fourdvar.F:876
mod_fourdvar::addotproduct
real(r8) addotproduct
Definition mod_fourdvar.F:869
mod_fourdvar::dotproduct
real(r8) dotproduct
Definition mod_fourdvar.F:868
mod_fourdvar::ig1count
integer ig1count
Definition mod_fourdvar.F:873
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::idfsur
integer idfsur
Definition mod_ncparam.F:206
mod_ncparam::tintrp
real(dp), dimension(:,:,:), allocatable tintrp
Definition mod_ncparam.F:583
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::master
logical master
Definition mod_parallel.F:40
mod_parallel::tile_count
integer tile_count
Definition mod_parallel.F:89
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::ntilex
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
mod_param::iadm
integer, parameter iadm
Definition mod_param.F:665
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::ntilee
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_scalars
Definition mod_scalars.F:2
mod_scalars::nrrec
integer, dimension(:), allocatable nrrec
Definition mod_scalars.F:1171
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::ntlm
integer, dimension(:), allocatable ntlm
Definition mod_scalars.F:1080
mod_scalars::bvf_bak
real(dp) bvf_bak
Definition mod_scalars.F:791
mod_scalars::tcoef
real(r8), dimension(:), allocatable tcoef
Definition mod_scalars.F:732
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::nhis
integer, dimension(:), allocatable nhis
Definition mod_scalars.F:1076
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::time
real(dp), dimension(:), allocatable time
Definition mod_scalars.F:262
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::ntstart
integer, dimension(:), allocatable ntstart
Definition mod_scalars.F:378
mod_scalars::inner
integer inner
Definition mod_scalars.F:360
mod_scalars::scoef
real(r8), dimension(:), allocatable scoef
Definition mod_scalars.F:733
mod_scalars::outer
integer outer
Definition mod_scalars.F:361
mod_stepping
Definition mod_stepping.F:2
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