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_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/LBC – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/LBC/lbcnfd_tam.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

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