New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
lbcnfd_tam.F90 in branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/LBC – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/LBC/lbcnfd_tam.F90 @ 4584

Last change on this file since 4584 was 4584, checked in by pabouttier, 10 years ago

Fix wrong loop in lbcnfd_tam, see Ticket #1278. Add calls to mpp_sum in adjoint tests
  • Property svn:executable set to *
File size: 35.7 KB
Line 
1MODULE lbcnfd_tam
2   !!======================================================================
3   !!                       ***  MODULE  lbcnfd_tam  ***
4   !! Ocean        : TAM of north fold  boundary conditions
5   !!======================================================================
6   !! History :  3.2  ! 2009-03  (R. Benshila)  Original code
7   !! History of TAM : 3.2 ! 2010-04 (F. Vigilant) Original Code
8   !!                  3.4 ! 2012-03 (P.-A. Bouttier) Update
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   lbc_nfd_adj       : generic interface for lbc_nfd_3d_adj and lbc_nfd_2d_adj routines
13   !!   lbc_nfd_3d_adj    : lateral boundary condition: North fold treatment for a 3D arrays   (lbc_nfd_adj)
14   !!   lbc_nfd_2d_adj    : lateral boundary condition: North fold treatment for a 2D arrays   (lbc_nfd_adj)
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE in_out_manager ! I/O manager
18
19   IMPLICIT NONE
20   PRIVATE
21
22   INTERFACE lbc_nfd_adj
23      MODULE PROCEDURE   lbc_nfd_3d_adj, lbc_nfd_2d_adj
24   END INTERFACE
25
26   PUBLIC   lbc_nfd_adj      ! north fold conditions
27#if defined key_tam
28   PUBLIC   lbc_nfd_adj_tst  ! Adjoint test routine
29#endif
30   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
32   !! $Id$
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE lbc_nfd_3d_adj( pt3d, cd_type, psgn )
38      !!----------------------------------------------------------------------
39      !!                  ***  routine lbc_nfd_3d_adj  ***
40      !!
41      !! ** Purpose :   Adjoint of 3D lateral boundary condition : North fold treatment
42      !!              without processor exchanges.
43      !!
44      !! ** Method  :
45      !!
46      !! ** Action  :   pt3d with updated values along the north fold
47      !!----------------------------------------------------------------------
48      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type   ! define the nature of ptab array grid-points
49      !                                                        !   = T , U , V , F , W points
50      REAL(wp)                  , INTENT(in   ) ::   psgn      ! control of the sign change
51      !                                                        !   = -1. , the sign is changed if north fold boundary
52      !                                                        !   =  1. , the sign is kept  if north fold boundary
53      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pt3d      ! 3D array on which the boundary condition is applied
54      !
55      INTEGER  ::   ji, jk
56      INTEGER  ::   ijt, iju, ijpj, ijpjm1
57      !
58      REAL(wp) :: ztmp
59      !!----------------------------------------------------------------------
60
61      SELECT CASE ( jpni )
62      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
63      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
64      END SELECT
65      ijpjm1 = ijpj-1
66
67      DO jk = 1, jpk
68         !
69         SELECT CASE ( npolj )
70         !
71         CASE ( 3 , 4 )                        ! *  North fold  T-point pivot
72            !
73            SELECT CASE ( cd_type )
74         CASE ( 'T' , 'W' )                         ! T-, W-point
75            DO ji = jpiglo, jpiglo/2+1, -1
76               ijt = jpiglo-ji+2
77               ztmp = psgn * pt3d(ji,ijpjm1,jk)
78               pt3d(ji ,ijpjm1,jk) = 0.0_wp
79               pt3d(ijt,ijpjm1,jk) = pt3d(ijt,ijpjm1,jk) + ztmp
80            END DO
81
82            DO ji = jpiglo, 2, -1
83               ijt = jpiglo-ji+2
84               pt3d(ijt,ijpj-2,jk) = pt3d(ijt,ijpj-2,jk) + psgn * pt3d(ji,ijpj,jk)
85               pt3d(ji ,ijpj  ,jk) = 0.0_wp
86            END DO
87         CASE ( 'U' )                               ! U-point
88            DO ji = jpiglo-1, jpiglo/2, -1
89               iju = jpiglo-ji+1
90               ztmp = psgn * pt3d(ji,ijpjm1,jk)
91               pt3d(ji ,ijpjm1,jk) = 0.0_wp
92               pt3d(iju,ijpjm1,jk) = pt3d(iju,ijpjm1,jk) + ztmp
93            END DO
94
95            DO ji = jpiglo-1, 1, -1
96               iju = jpiglo-ji+1
97               pt3d(iju,ijpj-2,jk) = pt3d(iju,ijpj-2,jk) + psgn * pt3d(ji,ijpj,jk)
98               pt3d(ji ,ijpj  ,jk) = 0.0_wp
99            END DO
100         CASE ( 'V' )                               ! V-point
101            DO ji = jpiglo, 2, -1
102               ijt = jpiglo-ji+2
103               pt3d(ijt,ijpj-3,jk) = pt3d(ijt,ijpj-3,jk) + psgn * pt3d(ji,ijpj  ,jk)
104               pt3d(ji,ijpj  ,jk) = 0.0_wp
105               pt3d(ijt,ijpj-2,jk) = pt3d(ijt,ijpj-2,jk) + psgn * pt3d(ji,ijpj-1,jk)
106               pt3d(ji,ijpj-1,jk) = 0.0_wp
107            END DO
108         CASE ( 'F' )                               ! F-point
109            DO ji = jpiglo-1, 1, -1
110               iju = jpiglo-ji+1
111               pt3d(iju,ijpj-3,jk) = pt3d(iju,ijpj-3,jk) + psgn * pt3d(ji,ijpj  ,jk)
112               pt3d(ji ,ijpj  ,jk) = 0.0_wp
113               pt3d(iju,ijpj-2,jk) = pt3d(iju,ijpj-2,jk) + psgn * pt3d(ji,ijpj-1,jk)
114               pt3d(ji ,ijpj-1,jk) = 0.0_wp
115            END DO
116         END SELECT
117            !
118         CASE ( 5 , 6 )                        ! *  North fold  F-point pivot
119            !
120            SELECT CASE ( cd_type )
121            CASE ( 'T' , 'W' )                         ! T-, W-point
122            DO ji = jpiglo, 1, -1
123               ijt = jpiglo-ji+1
124               pt3d(ijt,ijpj-1,jk) = pt3d(ijt,ijpj-1,jk) + psgn * pt3d(ji,ijpj,jk)
125               pt3d(ji ,ijpj  ,jk) = 0.0_wp
126            END DO
127         CASE ( 'U' )                               ! U-point
128            DO ji = jpiglo-1, 1, -1
129               iju = jpiglo-ji
130               pt3d(iju,ijpj-1,jk) = pt3d(iju,ijpj-1,jk) + psgn * pt3d(ji,ijpj,jk)
131               pt3d(ji ,ijpj  ,jk) = 0.0_wp
132            END DO
133         CASE ( 'V' )                               ! V-point
134            DO ji = jpiglo, jpiglo/2+1, -1
135               ijt = jpiglo-ji+1
136               ztmp = psgn * pt3d(ji,ijpjm1,jk)
137               pt3d(ji ,ijpjm1,jk) = 0.0_wp
138               pt3d(ijt,ijpjm1,jk) = pt3d(ijt,ijpjm1,jk) + ztmp
139            END DO
140            DO ji = jpiglo, 1, -1
141               ijt = jpiglo-ji+1
142               pt3d(ijt,ijpj-2,jk) = pt3d(ijt,ijpj-2,jk) + psgn * pt3d(ji,ijpj,jk)
143               pt3d(ji ,ijpj  ,jk) = 0.0_wp
144            END DO
145         CASE ( 'F' )                               ! F-point
146            DO ji = jpiglo-1, jpiglo/2+1, -1
147               iju = jpiglo-ji
148               ztmp = psgn * pt3d(ji,ijpjm1,jk)
149               pt3d(ji ,ijpjm1,jk) = 0.0_wp
150               pt3d(iju,ijpjm1,jk) = pt3d(iju,ijpjm1,jk) + ztmp
151            END DO
152            DO ji = jpiglo-1, 1, -1
153               iju = jpiglo-ji
154               pt3d(iju,ijpj-2,jk) = pt3d(iju,ijpj-2,jk) + psgn * pt3d(ji,ijpj  ,jk)
155               pt3d(ji ,ijpj  ,jk) = 0.0_wp
156            END DO
157         END SELECT
158            !
159         CASE DEFAULT                           ! *  closed : the code probably never go through
160            !
161            SELECT CASE ( cd_type)
162            CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
163               pt3d(:, 1  ,jk) = 0.e0
164               pt3d(:,ijpj,jk) = 0.e0
165            CASE ( 'F' )                               ! F-point
166               pt3d(:,ijpj,jk) = 0.e0
167            END SELECT
168            !
169         END SELECT     !  npolj
170         !
171      END DO
172      !
173   END SUBROUTINE lbc_nfd_3d_adj
174
175
176   SUBROUTINE lbc_nfd_2d_adj( pt2d, cd_type, psgn, pr2dj )
177      !!----------------------------------------------------------------------
178      !!                  ***  routine lbc_nfd_2d_adj  ***
179      !!
180      !! ** Purpose :   Adjoint of 2D lateral boundary condition : North fold treatment
181      !!       without processor exchanges.
182      !!
183      !! ** Method  :
184      !!
185      !! ** Action  :   pt2d with updated values along the north fold
186      !!----------------------------------------------------------------------
187      CHARACTER(len=1)        , INTENT(in   ) ::   cd_type   ! define the nature of ptab array grid-points
188      !                                                      ! = T , U , V , F , W points
189      REAL(wp)                , INTENT(in   ) ::   psgn      ! control of the sign change
190      !                                                      !   = -1. , the sign is changed if north fold boundary
191      !                                                      !   =  1. , the sign is kept  if north fold boundary
192      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d      ! 2D array on which the boundary condition is applied
193      INTEGER , OPTIONAL      , INTENT(in   ) ::   pr2dj     ! number of additional halos
194      !
195      INTEGER  ::   ji, jl, ipr2dj
196      INTEGER  ::   ijt, iju, ijpj, ijpjm1
197      !
198      REAL(wp) :: ztmp
199      !!----------------------------------------------------------------------
200
201      SELECT CASE ( jpni )
202      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
203      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
204      END SELECT
205      !
206      IF( PRESENT(pr2dj) ) THEN           ! use of additional halos
207         ipr2dj = pr2dj
208         IF( jpni > 1 )   ijpj = ijpj + ipr2dj
209      ELSE
210         ipr2dj = 0
211      ENDIF
212      !
213      ijpjm1 = ijpj-1
214      SELECT CASE ( npolj )
215      !
216      CASE ( 3, 4 )                       ! *  North fold  T-point pivot
217         !
218         SELECT CASE ( cd_type )
219         !
220         CASE ( 'T', 'W' )
221            DO ji = jpiglo, jpiglo/2+1, -1
222               ijt=jpiglo-ji+2
223               ztmp = psgn * pt2d(ji,ijpj-1)
224               pt2d(ji,ijpj-1) = 0.0_wp
225               pt2d(ijt,ijpj-1) = pt2d(ijt,ijpj-1) + ztmp
226            END DO
227            DO jl = ipr2dj, 0, -1
228               DO ji = jpiglo, 2, -1
229                  ijt=jpiglo-ji+2
230                  ztmp = psgn * pt2d(ji,ijpj+jl)
231                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
232                  pt2d(ijt,ijpj-2-jl) = pt2d(ijt,ijpj-2-jl) + ztmp
233               END DO
234            END DO
235         CASE ( 'U' )                                     ! U-point
236            DO ji = jpiglo-1, jpiglo/2, -1
237               iju = jpiglo-ji+1
238               ztmp = psgn * pt2d(ji,ijpjm1)
239               pt2d(ji,ijpjm1) = 0.0_wp
240               pt2d(iju,ijpjm1) = pt2d(iju,ijpjm1) + ztmp
241            END DO
242            DO jl = ipr2dj, 0, -1
243               DO ji = jpiglo-1, 1, -1
244                  iju = jpiglo-ji+1
245                  ztmp = psgn * pt2d(ji,ijpj+jl)
246                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
247                  pt2d(iju,ijpj-2-jl) = pt2d(iju,ijpj-2-jl) + ztmp
248               END DO
249            END DO
250         CASE ( 'V' )                                     ! V-point
251            DO jl = ipr2dj, -1, -1
252               DO ji = jpiglo, 2, -1
253                  ijt = jpiglo-ji+2
254                  ztmp = psgn * pt2d(ji,ijpj+jl)
255                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
256                  pt2d(ijt,ijpj-3-jl) = pt2d(ijt,ijpj-3-jl) + ztmp
257               END DO
258            END DO
259         CASE ( 'F' )                               ! F-point
260            DO jl = ipr2dj, -1, -1
261               DO ji = jpiglo-1, 1, -1
262                  iju = jpiglo-ji+1
263                  ztmp = psgn * pt2d(ji,ijpj+jl)
264                  pt2d(ji,ijpj+jl) = 0.0_wp
265                  pt2d(iju,ijpj-3-jl) = pt2d(iju,ijpj-3-jl) + ztmp
266               END DO
267            END DO
268         CASE ( 'I' )                                     ! ice U-V point
269            DO jl = ipr2dj, 0, -1
270               DO ji = jpiglo, 3, -1
271                  iju = jpiglo - ji + 3
272                  ztmp = psgn * pt2d(ji,ijpj+jl)
273                  pt2d(ji,ijpj+jl) = 0.0_wp
274                  pt2d(iju,ijpj-1-jl) = pt2d(iju,ijpj-1-jl) + ztmp
275               END DO
276               ztmp = psgn * pt2d(2,ijpj+jl)
277               pt2d(2,ijpj+jl) = 0.0_wp
278               pt2d(3,ijpj-1+jl) = pt2d(3,ijpj-1+jl) + ztmp
279            END DO
280         CASE ( 'J' )                                     ! first ice U-V point
281            DO jl =ipr2dj,0,-1
282               DO ji = 3, jpiglo
283                  iju = jpiglo - ji + 3
284                  ztmp = psgn * pt2d(ji,ijpj+jl)
285                  pt2d(ji,ijpj+jl) = 0.0_wp
286                  pt2d(iju,ijpj-1-jl) = pt2d(iju,ijpj-1-jl) + ztmp
287               END DO
288               pt2d(3,ijpj-1+jl) = pt2d(3,ijpj-1+jl) + psgn * pt2d(2,ijpj+jl)
289               pt2d(2,ijpj+jl) = 0.0_wp
290            END DO
291         CASE ( 'K' )                                     ! second ice U-V point
292            DO jl =0, ipr2dj
293               DO ji = 3, jpiglo
294                  iju = jpiglo - ji + 3
295                  pt2d(3,ijpj-1+jl) = pt2d(3,ijpj-1+jl) + psgn * pt2d(ji,ijpj+jl)
296                  pt2d(ji,ijpj+jl) = 0.0_wp
297               END DO
298               ztmp = psgn * pt2d(2,ijpj+jl)
299               pt2d(2,ijpj+jl) = 0.0_wp
300               pt2d(3,ijpj-1+jl) = pt2d(3,ijpj-1+jl) + ztmp
301            END DO
302         END SELECT
303         !
304      CASE ( 5, 6 )                        ! *  North fold  F-point pivot
305         !
306         SELECT CASE ( cd_type )
307         CASE ( 'T' , 'W' )                          ! T-, W-point
308            DO jl = ipr2dj, 0, -1
309               DO ji = jpiglo, 1, -1
310                  ijt = jpiglo-ji+1
311                  ztmp = psgn * pt2d(ji,ijpj+jl)
312                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
313                  pt2d(ijt,ijpj-1-jl) = pt2d(ijt,ijpj-1-jl) + ztmp
314               END DO
315            END DO
316         CASE ( 'U' )                                     ! U-point
317            DO jl = ipr2dj, 0, -1
318               DO ji = jpiglo-1, 1, -1
319                  iju = jpiglo-ji
320                  ztmp = psgn * pt2d(ji,ijpj+jl)
321                  pt2d(ji,ijpj+jl) = 0.0_wp
322                  pt2d(iju,ijpj-1-jl) = pt2d(iju,ijpj-1-jl) + ztmp
323               END DO
324            END DO
325         CASE ( 'V' )                                     ! V-point
326            DO ji = jpiglo, jpiglo/2+1, -1
327               ijt = jpiglo-ji+1
328               ztmp = psgn * pt2d(ji,ijpjm1)
329               pt2d(ji ,ijpjm1) = 0.0_wp
330               pt2d(ijt,ijpjm1) = pt2d(ijt,ijpjm1) + ztmp
331            END DO
332            DO jl = ipr2dj, 0, -1
333               DO ji = jpiglo, 1, -1
334                  ijt = jpiglo-ji+1
335                  ztmp = psgn * pt2d(ji,ijpj+jl)
336                  pt2d(ji,ijpj+jl) = 0.0_wp
337                  pt2d(ijt,ijpj-2-jl) = pt2d(ijt,ijpj-2-jl) + ztmp
338               END DO
339            END DO
340         CASE ( 'F' )                               ! F-point
341            DO ji = jpiglo-1, jpiglo/2+1, -1
342               iju = jpiglo-ji
343               ztmp = psgn * pt2d(ji,ijpjm1)
344               pt2d(ji ,ijpjm1) = 0.0_wp
345               pt2d(iju,ijpjm1) = pt2d(iju,ijpjm1) + ztmp
346            END DO
347            DO jl = ipr2dj, 0, -1
348               DO ji = jpiglo-1, 1, -1
349                  iju = jpiglo-ji
350                  ztmp = psgn * pt2d(ji,ijpj+jl)
351                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
352                  pt2d(iju,ijpj-2-jl) = pt2d(iju,ijpj-2-jl) + ztmp
353               END DO
354            END DO
355         CASE ( 'I' )                                  ! ice U-V point
356            DO jl = ipr2dj, 0, -1
357               DO ji = jpiglo-1, 2, -1
358                  ijt = jpiglo - ji + 2
359                  pt2d(ji ,ijpj-1-jl) = pt2d(ji ,ijpj-1-jl) + 0.5 * pt2d(ji,ijpj+jl)
360                  pt2d(ijt,ijpj-1-jl) = pt2d(ijt,ijpj-1-jl) + 0.5 * psgn * pt2d(ji,ijpj+jl)
361                  pt2d(ji ,ijpj+jl  ) = 0.0_wp
362               END DO
363            END DO
364            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.0_wp
365         CASE ( 'J' )                                  ! first ice U-V point
366            DO jl = 0, ipr2dj
367               DO ji = 2 , jpiglo-1
368                  ijt = jpiglo - ji + 2
369                  pt2d(ji,ijpj-1-jl) = pt2d(ji,ijpj-1-jl) + pt2d(ji,ijpj+jl)
370                  pt2d(ji,ijpj+jl) = 0.0_wp
371               END DO
372            END DO
373            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0
374         CASE ( 'K' )                                  ! second ice U-V point
375            DO jl = 0, ipr2dj
376               DO ji = 2 , jpiglo-1
377                  ijt = jpiglo - ji + 2
378                  pt2d(ijt,ijpj-1-jl) = pt2d(ijt,ijpj-1-jl) + pt2d(ji,ijpj+jl)
379                  pt2d(ji,ijpj+jl) = 0.0_wp
380               END DO
381            END DO
382            pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0
383         END SELECT
384         !
385      CASE DEFAULT                           ! *  closed : the code probably never go through
386         !
387         SELECT CASE ( cd_type)
388         CASE ( 'T' , 'U' , 'V' , 'W' )                 ! T-, U-, V-, W-points
389            pt2d(:, 1:1-ipr2dj     ) = 0.e0
390            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
391         CASE ( 'F' )                                   ! F-point
392            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
393         CASE ( 'I' )                                   ! ice U-V point
394            pt2d(:, 1:1-ipr2dj     ) = 0.e0
395            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
396         CASE ( 'J' )                                   ! first ice U-V point
397            pt2d(:, 1:1-ipr2dj     ) = 0.e0
398            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
399         CASE ( 'K' )                                   ! second ice U-V point
400            pt2d(:, 1:1-ipr2dj     ) = 0.e0
401            pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0
402         END SELECT
403         !
404      END SELECT
405      !
406   END SUBROUTINE lbc_nfd_2d_adj
407
408   SUBROUTINE lbc_nfd_3d_adj_tst( kumadt )
409      !!-----------------------------------------------------------------------
410      !!
411      !!                  ***  ROUTINE lbc_nfd_3d_adj_tst ***
412      !!
413      !! ** Purpose : Test the adjoint routine.
414      !!
415      !! ** Method  : Verify the scalar product
416      !!
417      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
418      !!
419      !!              where  L   = tangent routine
420      !!                     L^T = adjoint routine
421      !!                     W   = diagonal matrix of scale factors
422      !!                     dx  = input perturbation (random field)
423      !!                     dy  = L dx
424      !!
425      !! ** Action  :
426      !!
427      !! History :
428      !!        ! 2010-09 (F. Vigilant)
429      !!-----------------------------------------------------------------------
430      !! * Modules used
431   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
432      & grid_random
433   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
434      & dot_product
435   USE lbcnfd,         ONLY: &
436      & lbc_nfd
437   USE tstool_tam    , ONLY: &
438      & stdu,                &
439      & stdv,                &
440      & stdt,                &
441      & prntst_adj
442   USE dom_oce       , ONLY: & ! Ocean space and time domain
443      & tmask,               &
444      & umask,               &
445      & vmask,               &
446      & mig,                 &
447      & mjg,                 &
448      & nldi,                &
449      & nldj,                &
450      & nlei,                &
451      & nlej
452      !! * Arguments
453      INTEGER, INTENT(IN) :: &
454         & kumadt        ! Output unit
455
456      !! * Local declarations
457      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
458         & zu_tlin,    & ! Tangent input: ua_tl
459         & zv_tlin,    & ! Tangent input: va_tl
460         & zt_tlin,    & ! Tangent input: ub_tl
461         & zu_tlout,   & ! Tangent output: ua_tl
462         & zv_tlout,   & ! Tangent output: va_tl
463         & zt_tlout,   & ! Tangent output: ua_tl
464         & zu_adin,    & ! Tangent output: ua_ad
465         & zv_adin,    & ! Tangent output: va_ad
466         & zt_adin,    & ! Adjoint input: ua_ad
467         & z_adin,    & ! Adjoint input: va_ad
468         & zu_adout,   & ! Adjoint output: ua_ad
469         & zv_adout,   & ! Adjoint output: va_ad
470         & zt_adout,   & ! Adjoint oputput: ub_ad
471         & znu            ! 3D random field for u
472
473      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
474         & zu_tl,    & ! Tangent input: ua_tl
475         & zv_tl,    & ! Tangent input: va_tl
476         & zt_tl,    & ! Tangent input: ub_tl
477         & zu_ad,    & ! Tangent output: ua_ad
478         & zv_ad,    & ! Tangent output: va_ad
479         & zt_ad
480      REAL(wp) :: &
481         & zsp1,    &   ! scalar product involving the tangent routine
482         & zsp2         ! scalar product involving the adjoint routine
483      CHARACTER(len=1) :: &
484         & cd_type
485      INTEGER :: &
486         & indic, &
487         & istp
488      INTEGER :: &
489         & ji, &
490         & jj, &
491         & jk, &
492         & kmod, &
493         & jstp
494      CHARACTER (LEN=14) :: &
495         & cl_name
496      INTEGER ::             &
497         & zijpj, ziglo, zip, zdiff
498
499      ! Allocate memory
500
501      zijpj = 4
502
503      SELECT CASE ( jpni )
504      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
505      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
506      END SELECT
507
508      ALLOCATE( &
509         & zu_tlin(jpiglo,zijpj,jpk),  &
510         & zv_tlin(jpiglo,zijpj,jpk),  &
511         & zt_tlin(jpiglo,zijpj,jpk),  &
512         & zu_tlout(jpiglo,zijpj,jpk), &
513         & zv_tlout(jpiglo,zijpj,jpk), &
514         & zt_tlout(jpiglo,zijpj,jpk), &
515         & zu_adin(jpiglo,zijpj,jpk),  &
516         & zv_adin(jpiglo,zijpj,jpk),  &
517         & zt_adin(jpiglo,zijpj,jpk),  &
518         & zu_adout(jpiglo,zijpj,jpk), &
519         & zv_adout(jpiglo,zijpj,jpk), &
520         & zt_adout(jpiglo,zijpj,jpk), &
521         & znu(jpi,jpj,jpk)         &
522         & )
523
524      ALLOCATE( &
525         & zu_tl(jpiglo,zijpj,jpk),  &
526         & zv_tl(jpiglo,zijpj,jpk),  &
527         & zt_tl(jpiglo,zijpj,jpk),  &
528         & zu_ad(jpiglo,zijpj,jpk),  &
529         & zv_ad(jpiglo,zijpj,jpk),  &
530         & zt_ad(jpiglo,zijpj,jpk)   &
531         & )
532
533        zu_tlin (:,:,:) = 0.0_wp
534        zv_tlin (:,:,:) = 0.0_wp
535        zt_tlin (:,:,:) = 0.0_wp
536        zu_tlout(:,:,:) = 0.0_wp
537        zv_tlout(:,:,:) = 0.0_wp
538        zt_tlout(:,:,:) = 0.0_wp
539        zu_adin (:,:,:) = 0.0_wp
540        zv_adin (:,:,:) = 0.0_wp
541        zt_adin (:,:,:) = 0.0_wp
542        zu_adout(:,:,:) = 0.0_wp
543        zv_adout(:,:,:) = 0.0_wp
544        zt_adout(:,:,:) = 0.0_wp
545
546        zu_tl (:,:,:) = 0.0_wp
547        zv_tl (:,:,:) = 0.0_wp
548        zt_tl (:,:,:) = 0.0_wp
549        zu_ad (:,:,:) = 0.0_wp
550        zv_ad (:,:,:) = 0.0_wp
551        zt_ad (:,:,:) = 0.0_wp
552
553        ziglo = INT(jpiglo / (nlei - nldi) )
554        zdiff = nlei - nldi
555        !--------------------------------------------------------------------
556        ! Initialize the tangent input with random noise: dx
557        !--------------------------------------------------------------------
558        CALL grid_random( znu, 'U', 0.0_wp, stdu )
559        DO jk = 1, jpk
560           DO jj = 1, 4
561              DO ji = nldi, nlei
562                 DO zip = 1, ziglo
563                    zu_tlin(ji+(zip-1)*zdiff,jj,jk) = znu(ji,jj,jk)
564                 ENDDO
565              END DO
566           END DO
567        END DO
568        CALL grid_random( znu, 'V', 0.0_wp, stdu )
569        DO jk = 1, jpk
570           DO jj = 1, 4
571              DO ji = nldi, nlei
572                 DO zip = 1, ziglo
573                    zv_tlin(ji+(zip-1)*zdiff,jj,jk) = znu(ji,jj,jk)
574                 ENDDO
575              END DO
576           END DO
577        END DO
578        CALL grid_random( znu, 'T', 0.0_wp, stdt )
579        DO jk = 1, jpk
580           DO jj = 1, 4
581              DO ji = nldi, nlei
582                 DO zip = 1, ziglo
583                    zt_tlin(ji+(zip-1)*zdiff,jj,jk) = znu(ji,jj,jk)
584                 ENDDO
585              END DO
586           END DO
587        END DO
588
589        !--------------------------------------------------------------------
590        ! Call the tangent routine: dy = L dx
591        !--------------------------------------------------------------------
592        zu_tl(:,:,:) = zu_tlin(:,:,:)
593        zv_tl(:,:,:) = zv_tlin(:,:,:)
594        zt_tl(:,:,:) = zt_tlin(:,:,:)
595        CALL lbc_nfd( zu_tl, 'U',  -1.0_wp )
596        CALL lbc_nfd( zv_tl, 'V',  -1.0_wp )
597        CALL lbc_nfd( zt_tl, 'T',   1.0_wp )
598        zu_tlout(:,:,:) = zu_tl(:,:,:)
599        zv_tlout(:,:,:) = zv_tl(:,:,:)
600        zt_tlout(:,:,:) = zt_tl(:,:,:)
601
602        DO jk = 1, jpk
603           DO jj = 1,4
604              DO ji = nldi, nlei
605                 DO zip = 1, ziglo
606                    zu_adin(ji+(zip-1)*zdiff,jj,jk) = zu_tlout(ji+(zip-1)*ziglo,jj,jk) &
607                       &               * e1u(ji,jj) * e2u(ji,jj) &!* fse3u(ji,jj,jk) &
608                       &               * umask(ji,jj,jk)
609                    zv_adin(ji+(zip-1)*zdiff,jj,jk) = zv_tlout(ji+(zip-1)*ziglo,jj,jk) &
610                       &               * e1v(ji,jj) * e2v(ji,jj) &!* fse3v(ji,jj,jk) &
611                       &               * vmask(ji,jj,jk)
612                    zt_adin(ji+(zip-1)*zdiff,jj,jk) = zt_tlout(ji+(zip-1)*ziglo,jj,jk) &
613                       &               * e1t(ji,jj) * e2t(ji,jj) &!* fse3t(ji,jj,jk) &
614                       &               * tmask(ji,jj,jk)
615                 ENDDO
616              END DO
617           END DO
618        END DO
619
620        ! DOT_PRODUCT
621        zsp1 = sum( PACK(zu_tlout(:,:,:),.TRUE.) * &
622         &                       PACK( zu_adin(:,:,:),.TRUE.) )
623
624        zu_ad(:,:,:) = zu_adin(:,:,:)
625        zv_ad(:,:,:) = zv_adin(:,:,:)
626        zt_ad(:,:,:) = zt_adin(:,:,:)
627
628        CALL lbc_nfd_adj( zu_ad, 'U',  -1.0_wp )
629        CALL lbc_nfd_adj( zv_ad, 'V',  -1.0_wp )
630        CALL lbc_nfd_adj( zt_ad, 'T',   1.0_wp )
631
632        zu_adout(:,:,:) = zu_ad(:,:,:)
633        zv_adout(:,:,:) = zv_ad(:,:,:)
634        zt_adout(:,:,:) = zt_ad(:,:,:)
635
636        zsp2 = sum( PACK(zu_tlin(:,:,:),.TRUE.) * &
637         &                       PACK( zu_adout(:,:,:),.TRUE.) )
638
639        CALL mpp_sum( zsp1 )
640        CALL mpp_sum( zsp2 )
641
642        cl_name = 'lbc_nfd  U  3d'
643        CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
644
645        zsp1 = sum( PACK(zv_tlout(:,:,:),.TRUE.) * &
646         &                       PACK( zv_adin(:,:,:),.TRUE.) )
647
648        zsp2 = sum( PACK(zv_tlin(:,:,:),.TRUE.) * &
649         &                       PACK( zv_adout(:,:,:),.TRUE.) )
650
651        CALL mpp_sum( zsp1 )
652        CALL mpp_sum( zsp2 )
653        cl_name = 'lbc_nfd  V  3d'
654        CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
655
656        zsp1 = sum( PACK(zt_tlout(:,:,:),.TRUE.) * &
657         &                       PACK( zt_adin(:,:,:),.TRUE.) )
658
659        zsp2 = sum( PACK(zt_tlin(:,:,:),.TRUE.) * &
660         &                       PACK( zt_adout(:,:,:),.TRUE.) )
661
662        CALL mpp_sum( zsp1 )
663        CALL mpp_sum( zsp2 )
664
665      cl_name = 'lbc_nfd  T  3d'
666      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
667
668   END SUBROUTINE lbc_nfd_3d_adj_tst
669
670   SUBROUTINE lbc_nfd_2d_adj_tst( kumadt )
671      !!-----------------------------------------------------------------------
672      !!
673      !!                  ***  ROUTINE lbc_nfd_2d_adj_tst ***
674      !!
675      !! ** Purpose : Test the adjoint routine.
676      !!
677      !! ** Method  : Verify the scalar product
678      !!
679      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
680      !!
681      !!              where  L   = tangent routine
682      !!                     L^T = adjoint routine
683      !!                     W   = diagonal matrix of scale factors
684      !!                     dx  = input perturbation (random field)
685      !!                     dy  = L dx
686      !!
687      !! ** Action  :
688      !!
689      !! History :
690      !!        ! 2010-09 (F. Vigilant)
691      !!-----------------------------------------------------------------------
692      !! * Modules used
693   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
694      & grid_random
695   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
696      & dot_product
697   USE lbcnfd,         ONLY: &
698      & lbc_nfd
699   USE tstool_tam    , ONLY: &
700      & stdu,                &
701      & stdv,                &
702      & stdt,                &
703      & prntst_adj
704   USE dom_oce       , ONLY: & ! Ocean space and time domain
705      & tmask,               &
706      & umask,               &
707      & vmask,               &
708      & mig,                 &
709      & mjg,                 &
710      & nldi,                &
711      & nldj,                &
712      & nlei,                &
713      & nlej
714      !! * Arguments
715      INTEGER, INTENT(IN) :: &
716         & kumadt        ! Output unit
717
718      !! * Local declarations
719      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
720         & zu_tlin,    & ! Tangent input: ua_tl
721         & zv_tlin,    & ! Tangent input: va_tl
722         & zt_tlin,    & ! Tangent input: ub_tl
723         & zu_tlout,   & ! Tangent output: ua_tl
724         & zv_tlout,   & ! Tangent output: va_tl
725         & zt_tlout,   & ! Tangent output: ua_tl
726         & zu_adin,    & ! Tangent output: ua_ad
727         & zv_adin,    & ! Tangent output: va_ad
728         & zt_adin,    & ! Adjoint input: ua_ad
729         & z_adin,    & ! Adjoint input: va_ad
730         & zu_adout,   & ! Adjoint output: ua_ad
731         & zv_adout,   & ! Adjoint output: va_ad
732         & zt_adout,   & ! Adjoint oputput: ub_ad
733         & znu            ! 3D random field for u
734
735      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
736         & zu_tl,    & ! Tangent input: ua_tl
737         & zv_tl,    & ! Tangent input: va_tl
738         & zt_tl,    & ! Tangent input: ub_tl
739         & zu_ad,    & ! Tangent output: ua_ad
740         & zv_ad,    & ! Tangent output: va_ad
741         & zt_ad
742      REAL(wp) :: &
743         & zsp1,    &   ! scalar product involving the tangent routine
744         & zsp2         ! scalar product involving the adjoint routine
745      INTEGER, DIMENSION(jpi,jpj) :: &
746         & iseed_2d    ! 2D seed for the random number generator
747      CHARACTER(len=1) :: &
748         & cd_type
749      INTEGER :: &
750         & indic, &
751         & istp
752      INTEGER :: &
753         & ji, &
754         & jj, &
755         & jk, &
756         & kmod, &
757         & jstp
758      CHARACTER (LEN=14) :: &
759         & cl_name
760      INTEGER ::             &
761         & zijpj, ziglo, zip, zdiff
762
763      ! Allocate memory
764
765      SELECT CASE ( jpni )
766      CASE ( 1 )     ;   ijpj = nlcj      ! 1 proc only  along the i-direction
767      CASE DEFAULT   ;   ijpj = 4         ! several proc along the i-direction
768      END SELECT
769
770      ALLOCATE( &
771         & zu_tlin(jpiglo,zijpj),  &
772         & zv_tlin(jpiglo,zijpj),  &
773         & zt_tlin(jpiglo,zijpj),  &
774         & zu_tlout(jpiglo,zijpj), &
775         & zv_tlout(jpiglo,zijpj), &
776         & zt_tlout(jpiglo,zijpj), &
777         & zu_adin(jpiglo,zijpj),  &
778         & zv_adin(jpiglo,zijpj),  &
779         & zt_adin(jpiglo,zijpj),  &
780         & zu_adout(jpiglo,zijpj), &
781         & zv_adout(jpiglo,zijpj), &
782         & zt_adout(jpiglo,zijpj), &
783         & znu(jpi,jpj)         &
784         & )
785
786      ALLOCATE( &
787         & zu_tl(jpiglo,zijpj),  &
788         & zv_tl(jpiglo,zijpj),  &
789         & zt_tl(jpiglo,zijpj),  &
790         & zu_ad(jpiglo,zijpj),  &
791         & zv_ad(jpiglo,zijpj),  &
792         & zt_ad(jpiglo,zijpj)   &
793         & )
794
795        zu_tlin (:,:) = 0.0_wp
796        zv_tlin (:,:) = 0.0_wp
797        zt_tlin (:,:) = 0.0_wp
798        zu_tlout(:,:) = 0.0_wp
799        zv_tlout(:,:) = 0.0_wp
800        zt_tlout(:,:) = 0.0_wp
801        zu_adin (:,:) = 0.0_wp
802        zv_adin (:,:) = 0.0_wp
803        zt_adin (:,:) = 0.0_wp
804        zu_adout(:,:) = 0.0_wp
805        zv_adout(:,:) = 0.0_wp
806        zt_adout(:,:) = 0.0_wp
807
808        zu_tl (:,:) = 0.0_wp
809        zv_tl (:,:) = 0.0_wp
810        zt_tl (:,:) = 0.0_wp
811        zu_ad (:,:) = 0.0_wp
812        zv_ad (:,:) = 0.0_wp
813        zt_ad (:,:) = 0.0_wp
814
815        ziglo = INT(jpiglo / (nlei - nldi) )
816        zdiff = nlei - nldi
817        !--------------------------------------------------------------------
818        ! Initialize the tangent input with random noise: dx
819        !--------------------------------------------------------------------
820        CALL grid_random( znu, 'U', 0.0_wp, stdu )
821        DO jj = 1, 4
822           DO ji = nldi, nlei
823              DO zip = 1, ziglo
824                 zu_tlin(ji+(zip-1)*zdiff,jj) = znu(ji,jj)
825              END DO
826           END DO
827        END DO
828        CALL grid_random( znu, 'V', 0.0_wp, stdu )
829        DO jj = 1, 4
830           DO ji = nldi, nlei
831              DO zip = 1, ziglo
832                 zv_tlin(ji+(zip-1)*zdiff,jj) = znu(ji,jj)
833              END DO
834           END DO
835        END DO
836        CALL grid_random(  znu, 'T', 0.0_wp, stdt )
837        DO jj = 1, 4
838           DO ji = nldi, nlei
839              DO zip = 1, ziglo
840                 zt_tlin(ji+(zip-1)*zdiff,jj) = znu(ji,jj)
841              END DO
842           END DO
843        END DO
844
845        !--------------------------------------------------------------------
846        ! Call the tangent routine: dy = L dx
847        !--------------------------------------------------------------------
848        zu_tl(:,:) = zu_tlin(:,:)
849        zv_tl(:,:) = zv_tlin(:,:)
850        zt_tl(:,:) = zt_tlin(:,:)
851
852        CALL lbc_nfd( zu_tl, 'U',  -1.0_wp )
853        CALL lbc_nfd( zv_tl, 'V',  -1.0_wp )
854        CALL lbc_nfd( zt_tl, 'T',   1.0_wp )
855        zu_tlout(:,:) = zu_tl(:,:)
856        zv_tlout(:,:) = zv_tl(:,:)
857        zt_tlout(:,:) = zt_tl(:,:)
858
859        DO jj = 1,4
860           DO ji = nldi, nlei
861              DO zip = 1, ziglo
862              zu_adin(ji+(zip-1)*zdiff,jj) = zu_tlout(ji+(zip-1)*ziglo,jj) &
863                 &               * e1u(ji,jj) * e2u(ji,jj) &!* fse3u(ji,jj) &
864                 &               * umask(ji,jj,1)
865              zv_adin(ji+(zip-1)*zdiff,jj) = zv_tlout(ji+(zip-1)*ziglo,jj) &
866                 &               * e1v(ji,jj) * e2v(ji,jj) &!* fse3v(ji,jj) &
867                 &               * vmask(ji,jj,1)
868              zt_adin(ji+(zip-1)*zdiff,jj) = zt_tlout(ji+(zip-1)*ziglo,jj) &
869                 &               * e1t(ji,jj) * e2t(ji,jj) &!* fse3t(ji,jj) &
870                 &               * tmask(ji,jj,1)
871              END DO
872           END DO
873        END DO
874
875        zsp1 = sum( PACK(zu_tlout(:,:),.TRUE.) * &
876         &                       PACK( zu_adin(:,:),.TRUE.) )
877        zu_ad(:,:) = zu_adin(:,:)
878        zv_ad(:,:) = zv_adin(:,:)
879        zt_ad(:,:) = zt_adin(:,:)
880
881        CALL lbc_nfd_adj( zu_ad, 'U',  -1.0_wp )
882        CALL lbc_nfd_adj( zv_ad, 'V',  -1.0_wp )
883        CALL lbc_nfd_adj( zt_ad, 'T',   1.0_wp )
884
885        zu_adout(:,:) = zu_ad(:,:)
886        zv_adout(:,:) = zv_ad(:,:)
887        zt_adout(:,:) = zt_ad(:,:)
888
889        zsp2 = sum( PACK(zu_tlin(:,:),.TRUE.) * &
890         &                       PACK( zu_adout(:,:),.TRUE.) )
891
892        CALL mpp_sum( zsp1 )
893        CALL mpp_sum( zsp2 )
894
895        cl_name = 'lbc_nfd  U  2d'
896        CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
897
898        zsp1 = sum( PACK(zv_tlout(:,:),.TRUE.) * &
899         &                       PACK( zv_adin(:,:),.TRUE.) )
900
901        zsp2 = sum( PACK(zv_tlin(:,:),.TRUE.) * &
902         &                       PACK( zv_adout(:,:),.TRUE.) )
903       
904        CALL mpp_sum( zsp1 )
905        CALL mpp_sum( zsp2 )
906
907        cl_name = 'lbc_nfd  V  2d'
908        CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
909
910        zsp1 = sum( PACK(zt_tlout(:,:),.TRUE.) * &
911         &                       PACK( zt_adin(:,:),.TRUE.) )
912
913        zsp2 = sum( PACK(zt_tlin(:,:),.TRUE.) * &
914         &                       PACK( zt_adout(:,:),.TRUE.) )
915
916        CALL mpp_sum( zsp1 )
917        CALL mpp_sum( zsp2 )
918
919        cl_name = 'lbc_nfd  T  2d'
920        CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
921
922   END SUBROUTINE lbc_nfd_2d_adj_tst
923
924   SUBROUTINE lbc_nfd_adj_tst( kumadt )
925      !!-----------------------------------------------------------------------
926      !!
927      !!                  ***  ROUTINE lbc_nfd_adj_tst ***
928      !!
929      !! ** Purpose : Test the adjoint routine.
930      !!
931      !! ** Method  : Verify the scalar product
932      !!
933      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
934      !!
935      !!              where  L   = tangent routine
936      !!                     L^T = adjoint routine
937      !!                     W   = diagonal matrix of scale factors
938      !!                     dx  = input perturbation (random field)
939      !!                     dy  = L dx
940      !!
941      !! ** Action  :
942      !!
943      !! History :
944      !!        ! 2010-09 (F. Vigilant)
945      !!-----------------------------------------------------------------------
946      !! * Modules used
947      !! * Arguments
948      INTEGER, INTENT(IN) :: &
949         & kumadt        ! Output unit
950
951      CALL lbc_nfd_3d_adj_tst( kumadt )
952      CALL lbc_nfd_2d_adj_tst( kumadt )
953
954   END SUBROUTINE lbc_nfd_adj_tst
955   !!======================================================================
956END MODULE lbcnfd_tam
Note: See TracBrowser for help on using the repository browser.