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.
dynnxt_tam.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN/dynnxt_tam.F90 @ 1885

Last change on this file since 1885 was 1885, checked in by rblod, 14 years ago

add TAM sources

File size: 33.1 KB
Line 
1MODULE dynnxt_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  dynnxt_tam  ***
5   !! Ocean dynamics: time stepping
6   !!                 Tangent and Adjoint Module
7   !!======================================================================
8   
9   !!----------------------------------------------------------------------
10   !!   dyn_nxt_tan  : update the horizontal velocity from the momentum trend
11   !!   dyn_nxt_adj  : update the horizontal velocity from the momentum trend
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE par_kind      , ONLY: & ! Precision variables
15      & wp
16   USE par_oce       , ONLY: & ! Ocean space and time domain variables
17      & jpi,                 &
18      & jpj,                 & 
19      & jpk,                 & 
20      & jpkm1,               & 
21      & jpiglo
22   USE oce_tam       , ONLY: & ! ocean dynamics and tracers
23      & un_tl,               &
24      & vn_tl,               &
25      & ub_tl,               &
26      & vb_tl,               &
27      & ua_tl,               &
28      & va_tl,               &
29      & un_ad,               &
30      & vn_ad,               &
31      & ub_ad,               &
32      & vb_ad,               &
33      & ua_ad,               &
34      & va_ad
35   USE dom_oce       , ONLY: & ! ocean space and time domain
36      & umask,               &
37      & vmask,               &
38      & lk_vvl,              &
39      & neuler,              &
40      & rdt,                 & 
41      & atfp,                & 
42      & atfp1,               & 
43      & e1u,                 &
44      & e2u,                 &
45      & e1v,                 &
46      & e2v,                 &
47#if defined key_zco
48      & e3t_0,               &
49#else
50      & e3v,                 &
51      & e3u,                 &
52#endif
53      & mig,                 &
54      & mjg,                 &
55      & nldi,                &
56      & nldj,                &
57      & nlei,                &
58      & nlej   
59   USE in_out_manager, ONLY: & ! I/O manager
60      & lwp,                 &
61      & nit000,              &
62      & nitend,              &
63      & numout,              &
64      & ctl_stop
65   USE dynspg_oce    , ONLY: & ! type of surface pressure gradient
66      & lk_dynspg_flt
67   USE lbclnk        , ONLY: & ! lateral boundary condition (or mpp link)
68      & lbc_lnk
69   USE lbclnk_tam    , ONLY: & ! lateral boundary condition (or mpp link)
70      & lbc_lnk_adj
71   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
72      & grid_random
73   USE dotprodfld,     ONLY: & ! Computes dot product for 3D and 2D fields
74      & dot_product
75   USE tstool_tam    , ONLY: &
76      & prntst_adj,          & !
77                               ! random field standard deviation for:
78      & stdu,                & !   u-velocity
79      & stdv                   !   v-velocity
80
81!VERIF
82!   USE bdy_oce         ! unstructured open boundary conditions
83!   USE bdydta          ! unstructured open boundary conditions
84!   USE bdydyn          ! unstructured open boundary conditions
85!!!!!
86
87   IMPLICIT NONE
88   PRIVATE
89
90   !! * Accessibility
91   PUBLIC dyn_nxt_tan            ! routine called by step.F90
92   PUBLIC dyn_nxt_adj            ! routine called by step.F90
93   PUBLIC dyn_nxt_adj_tst        ! routine called by step.F90
94   !! * Substitutions
95#  include "domzgr_substitute.h90"
96   !!----------------------------------------------------------------------
97
98CONTAINS
99
100   SUBROUTINE dyn_nxt_tan ( kt )
101      !!----------------------------------------------------------------------
102      !!                  ***  ROUTINE dyn_nxt_tan  ***
103      !!                   
104      !! ** Purpose of the direct routine:
105      !!      Compute the after horizontal velocity from the
106      !!      momentum trend.
107      !!
108      !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va)
109      !!      through calls to routine lbc_lnk.
110      !!      After velocity is compute using a leap-frog scheme environment:
111      !!         (ua,va) = (ub,vb) + 2 rdt (ua,va)
112      !!      Note that if lk_dynspg_flt=T, the time stepping has already been
113      !!      performed in dynspg module
114      !!      Time filter applied on now horizontal velocity to avoid the
115      !!      divergence of two consecutive time-steps and swap of dynamics
116      !!      arrays to start the next time step:
117      !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
118      !!         (un,vn) = (ua,va)
119      !!
120      !! ** Action : - Update ub,vb arrays, the before horizontal velocity
121      !!             - Update un,vn arrays, the now horizontal velocity
122      !!
123      !! History of the direct routine:
124      !!        !  87-02  (P. Andrich, D. L Hostis)  Original code
125      !!        !  90-10  (C. Levy, G. Madec)
126      !!        !  93-03  (M. Guyon)  symetrical conditions
127      !!        !  97-02  (G. Madec & M. Imbard)  opa, release 8.0
128      !!        !  97-04  (A. Weaver)  Euler forward step
129      !!        !  97-06  (G. Madec)  lateral boudary cond., lbc routine
130      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
131      !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
132      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
133      !!        !  07-07  (D. Storkey) Calls to BDY routines.
134      !! History of the TAM routine:
135      !!   9.0  !  08-06  (A. Vidard) Skeleton
136      !!        !  08-08  (A. Vidard) tangent of the 05-11 version
137      !!        !  08-08  (A. Vidard) tangent of the 07-07 version
138      !!----------------------------------------------------------------------
139      !! * Arguments
140      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
141
142      !! * Local declarations
143      INTEGER  ::   ji, jj, jk   ! dummy loop indices
144      REAL(wp) ::   z2dt         ! temporary scalar
145      REAL(wp) ::   zsshun1, zsshvn1         ! temporary scalar
146      !!----------------------------------------------------------------------
147      !!  OPA 9.0 , LOCEAN-IPSL (2005)
148      !! $Header$
149      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
150      !!----------------------------------------------------------------------
151
152      IF( kt == nit000 ) THEN
153         IF(lwp) WRITE(numout,*)
154         IF(lwp) WRITE(numout,*) 'dyn_nxt_tan : time stepping'
155         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
156      ENDIF
157
158      ! Local constant initialization
159      z2dt = 2. * rdt
160      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
161
162      !! Explicit physics with thickness weighted updates
163      IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
164         CALL ctl_stop ('vvl not available for tangent')
165      ENDIF
166
167      ! Lateral boundary conditions on ( ua, va )
168      CALL lbc_lnk( ua_tl, 'U', -1. )
169      CALL lbc_lnk( va_tl, 'V', -1. )
170
171      !                                                ! ===============
172      DO jk = 1, jpkm1                                 ! Horizontal slab
173         !                                             ! ===============
174         ! Next velocity
175         ! -------------
176#if defined key_dynspg_flt
177         ! Leap-frog time stepping already done in dynspg.F routine
178#else
179         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
180            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
181               ! Leap-frog time stepping
182               ua_tl(ji,jj,jk) = ( ub_tl(ji,jj,jk) + z2dt * ua_tl(ji,jj,jk) ) * umask(ji,jj,jk)
183               va_tl(ji,jj,jk) = ( vb_tl(ji,jj,jk) + z2dt * va_tl(ji,jj,jk) ) * vmask(ji,jj,jk)
184            END DO
185         END DO
186
187         IF( lk_vvl ) THEN
188            CALL ctl_stop ('key_vvl not available yet for tangent')
189         ENDIF
190# if defined key_obc
191         !                                             ! ===============
192      END DO                                           !   End of slab
193      !                                                ! ===============
194      CALL ctl_stop ('key_obc not available yet for tangent')
195
196
197      !                                                ! ===============
198      DO jk = 1, jpkm1                                 ! Horizontal slab
199         !                                             ! ===============
200# elif defined key_bdy 
201         !                                             ! ===============
202      END DO                                           !   End of slab
203      !                                                ! ===============
204      CALL ctl_stop ('key_bdy not available yet for tangent')
205
206      !                                                ! ===============
207      DO jk = 1, jpkm1                                 ! Horizontal slab
208         !                                             ! ===============
209# endif
210# if defined key_agrif
211         !                                             ! ===============
212      END DO                                           !   End of slab
213      !                                                ! ===============
214      CALL ctl_stop ('key_agrif not available yet for tangent')
215      !                                                ! ===============
216      DO jk = 1, jpkm1                                 ! Horizontal slab
217         !                                             ! ===============
218# endif
219#endif
220
221         ! Time filter and swap of dynamics arrays
222         ! ------------------------------------------
223         IF( neuler == 0 .AND. kt == nit000 ) THEN
224            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
225               CALL ctl_stop ('vvl not available for tangent')
226            ELSE                                               ! Fixed levels
227               DO jj = 1, jpj
228                  DO ji = 1, jpi
229                     ! Euler (forward) time stepping
230                     ub_tl(ji,jj,jk) = un_tl(ji,jj,jk)
231                     vb_tl(ji,jj,jk) = vn_tl(ji,jj,jk)
232                     un_tl(ji,jj,jk) = ua_tl(ji,jj,jk)
233                     vn_tl(ji,jj,jk) = va_tl(ji,jj,jk)
234                  END DO
235               END DO
236            ENDIF
237         ELSE
238            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
239               CALL ctl_stop ('vvl not available for tangent')
240            ELSE                                               ! Fixed levels
241               DO jj = 1, jpj
242                  DO ji = 1, jpi
243                     ! Leap-frog time stepping
244                     ub_tl(ji,jj,jk) = atfp  * ( ub_tl(ji,jj,jk) + ua_tl(ji,jj,jk) ) + atfp1 * un_tl(ji,jj,jk)
245                     vb_tl(ji,jj,jk) = atfp  * ( vb_tl(ji,jj,jk) + va_tl(ji,jj,jk) ) + atfp1 * vn_tl(ji,jj,jk)
246                     un_tl(ji,jj,jk) = ua_tl(ji,jj,jk)
247                     vn_tl(ji,jj,jk) = va_tl(ji,jj,jk)
248                  END DO
249               END DO
250            ENDIF
251         ENDIF
252         !                                             ! ===============
253      END DO                                           !   End of slab
254      !                                                ! ===============
255
256#if defined key_agrif
257      CALL ctl_stop ('key_agrif not available yet for tangent')
258     
259!      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
260#endif     
261
262   END SUBROUTINE dyn_nxt_tan
263   SUBROUTINE dyn_nxt_adj ( kt )
264      !!----------------------------------------------------------------------
265      !!                  ***  ROUTINE dyn_nxt_adj  ***
266      !!                   
267      !! ** Purpose of the direct routine:
268      !!      Compute the after horizontal velocity from the
269      !!      momentum trend.
270      !!
271      !! ** Method  :   Apply lateral boundary conditions on the trends (ua,va)
272      !!      through calls to routine lbc_lnk.
273      !!      After velocity is compute using a leap-frog scheme environment:
274      !!         (ua,va) = (ub,vb) + 2 rdt (ua,va)
275      !!      Note that if lk_dynspg_flt=T, the time stepping has already been
276      !!      performed in dynspg module
277      !!      Time filter applied on now horizontal velocity to avoid the
278      !!      divergence of two consecutive time-steps and swap of dynamics
279      !!      arrays to start the next time step:
280      !!         (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
281      !!         (un,vn) = (ua,va)
282      !!
283      !! ** Action : - Update ub,vb arrays, the before horizontal velocity
284      !!             - Update un,vn arrays, the now horizontal velocity
285      !!
286      !! History of the direct routine:
287      !!        !  87-02  (P. Andrich, D. L Hostis)  Original code
288      !!        !  90-10  (C. Levy, G. Madec)
289      !!        !  93-03  (M. Guyon)  symetrical conditions
290      !!        !  97-02  (G. Madec & M. Imbard)  opa, release 8.0
291      !!        !  97-04  (A. Weaver)  Euler forward step
292      !!        !  97-06  (G. Madec)  lateral boudary cond., lbc routine
293      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
294      !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
295      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
296      !! History of the TAM routine:
297      !!   9.0  !  08-06  (A. Vidard) Skeleton
298      !!        !  08-08  (A. Vidard) Adjoint of the 02-10 version
299      !!        !  08-11  (A. Vidard) Adjoint of the 05-12 version
300      !!----------------------------------------------------------------------
301      !! * Arguments
302      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
303
304      !! * Local declarations
305      !! * Local declarations
306      INTEGER  ::   ji, jj, jk   ! dummy loop indices
307      REAL(wp) ::   z2dt         ! temporary scalar
308
309      IF( kt == nit000 ) THEN
310         IF(lwp) WRITE(numout,*)
311         IF(lwp) WRITE(numout,*) 'dyn_nxt_adj : time stepping'
312         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
313      ENDIF
314
315      ! Local constant initialization
316      z2dt = 2. * rdt
317      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
318
319#if defined key_agrif
320      CALL ctl_stop ('key_agrif not available yet for tangent')
321     
322!      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
323#endif     
324
325
326      !                                                ! ===============
327      DO jk = jpkm1, 1, -1                             ! Horizontal slab
328         !                                             ! ===============
329         ! Next velocity
330         ! -------------
331
332         ! Time filter and swap of dynamics arrays
333         ! ------------------------------------------
334         IF( neuler == 0 .AND. kt == nit000 ) THEN
335            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
336               CALL ctl_stop ('vvl not available for tangent')
337            ELSE                                               ! Fixed levels
338               DO jj = 1, jpj
339                  DO ji = 1, jpi
340                     ! Euler (forward) time stepping
341                     ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + un_ad(ji,jj,jk)
342                     va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + vn_ad(ji,jj,jk)
343                     un_ad(ji,jj,jk) = 0.0_wp
344                     vn_ad(ji,jj,jk) = 0.0_wp
345                     un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + ub_ad(ji,jj,jk)
346                     vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + vb_ad(ji,jj,jk)
347                     ub_ad(ji,jj,jk) = 0.0_wp
348                     vb_ad(ji,jj,jk) = 0.0_wp
349                  END DO
350               END DO
351            ENDIF
352         ELSE
353            IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN      ! Varying levels
354               CALL ctl_stop ('vvl not available for tangent')
355            ELSE                                               ! Fixed levels
356               DO jj = 1, jpj
357                  DO ji = 1, jpi
358                     ! Leap-frog time stepping
359                     ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + un_ad(ji,jj,jk)
360                     un_ad(ji,jj,jk) = 0.0_wp
361                     va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + vn_ad(ji,jj,jk)
362                     vn_ad(ji,jj,jk) = 0.0_wp
363                     ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + atfp  * ub_ad(ji,jj,jk)
364                     un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + atfp1 * ub_ad(ji,jj,jk)
365                     ub_ad(ji,jj,jk) =                   atfp  * ub_ad(ji,jj,jk)
366
367                     va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + atfp  * vb_ad(ji,jj,jk)
368                     vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + atfp1 * vb_ad(ji,jj,jk)
369                     vb_ad(ji,jj,jk) =                   atfp  * vb_ad(ji,jj,jk)
370
371                  END DO
372               END DO
373            ENDIF
374         ENDIF
375
376#if defined key_dynspg_flt
377         ! Leap-frog time stepping already done in dynspg.F routine
378#else
379# if defined key_agrif
380         !                                             ! ===============
381      END DO                                           !   End of slab
382      !                                                ! ===============
383      CALL ctl_stop ('key_agrif not available yet for adjoint')
384      !                                                ! ===============
385      DO jk = 1, jpkm1                                 ! Horizontal slab
386         !                                             ! ===============
387# endif
388# if defined key_bdy 
389         !                                             ! ===============
390      END DO                                           !   End of slab
391      !                                                ! ===============
392      CALL ctl_stop ('key_bdy not available yet for adjoint')
393
394      !                                                ! ===============
395      DO jk = 1, jpkm1                                 ! Horizontal slab
396         !                                             ! ===============
397# elif defined key_obc
398         !                                             ! ===============
399      END DO                                           !   End of slab
400      !                                                ! ===============
401      CALL ctl_stop ('key_obc not available yet for adjoint')
402
403
404      !                                                ! ===============
405      DO jk = 1, jpkm1                                 ! Horizontal slab
406         !                                             ! ===============
407# endif
408         IF( lk_vvl ) THEN
409            CALL ctl_stop ('key_vvl not available yet for adjoint')
410         ENDIF
411
412         DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
413            DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
414               ! Leap-frog time stepping
415               ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk) * umask(ji,jj,jk)
416               ua_ad(ji,jj,jk) =            z2dt * ua_ad(ji,jj,jk) * umask(ji,jj,jk)
417               vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk) * vmask(ji,jj,jk)
418               va_ad(ji,jj,jk) =            z2dt * va_ad(ji,jj,jk) * vmask(ji,jj,jk)
419            END DO
420         END DO 
421#endif 
422!                                             ! ===============
423      END DO                                           !   End of slab
424      !                                                ! ===============
425      ! Lateral boundary conditions on ( ua, va )
426      CALL lbc_lnk_adj( ua_ad, 'U', -1. )
427      CALL lbc_lnk_adj( va_ad, 'V', -1. )
428
429      !! Explicit physics with thickness weighted updates
430      IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
431         CALL ctl_stop ('vvl not available for tangent')
432      ENDIF
433
434   END SUBROUTINE dyn_nxt_adj
435   SUBROUTINE dyn_nxt_adj_tst( kumadt )
436      !!-----------------------------------------------------------------------
437      !!
438      !!                  ***  ROUTINE dyn_nxt_adj_tst ***
439      !!
440      !! ** Purpose : Test the adjoint routine.
441      !!
442      !! ** Method  : Verify the scalar product
443      !!           
444      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
445      !!
446      !!              where  L   = tangent routine
447      !!                     L^T = adjoint routine
448      !!                     W   = diagonal matrix of scale factors
449      !!                     dx  = input perturbation (random field)
450      !!                     dy  = L dx
451      !!
452      !! ** Action  : Separate tests are applied for the following dx and dy:
453      !!               
454      !!              1) dx = ( SSH ) and dy = ( SSH )
455      !!                   
456      !! History :
457      !!        ! 08-08 (A. Vidard)
458      !!-----------------------------------------------------------------------
459      !! * Modules used
460
461      !! * Arguments
462      INTEGER, INTENT(IN) :: &
463         & kumadt             ! Output unit
464 
465      INTEGER :: &
466         & ji,    &        ! dummy loop indices
467         & jj,    &       
468         & jk     
469      INTEGER, DIMENSION(jpi,jpj) :: &
470         & iseed_2d        ! 2D seed for the random number generator
471
472      !! * Local declarations
473      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
474         & zun_tlin,     & ! Tangent input:   now   u-velocity
475         & zvn_tlin,     & ! Tangent input:   now   v-velocity
476         & zua_tlin,     & ! Tangent input:  after  u-velocity
477         & zva_tlin,     & ! Tangent input:  after  u-velocity         
478         & zub_tlin,     & ! Tangent input:  before u-velocity
479         & zvb_tlin,     & ! Tangent input:  before u-velocity         
480         & zun_adin,     & ! Adjoint input:   now   u-velocity
481         & zvn_adin,     & ! Adjoint input:   now   v-velocity
482         & zua_adin,     & ! Adjoint input:  after  u-velocity
483         & zva_adin,     & ! Adjoint input:  after  u-velocity         
484         & zub_adin,     & ! Adjoint input:  before u-velocity
485         & zvb_adin,     & ! Adjoint input:  before u-velocity         
486         & zun_tlout,    & ! Tangent output:  now   u-velocity
487         & zvn_tlout,    & ! Tangent output:  now   v-velocity
488         & zua_tlout,    & ! Tangent output: after  u-velocity
489         & zva_tlout,    & ! Tangent output: after  u-velocity         
490         & zub_tlout,    & ! Tangent output: before u-velocity
491         & zvb_tlout,    & ! Tangent output: before u-velocity         
492         & zun_adout,    & ! Adjoint output:  now   u-velocity
493         & zvn_adout,    & ! Adjoint output:  now   v-velocity
494         & zua_adout,    & ! Adjoint output: after  u-velocity
495         & zva_adout,    & ! Adjoint output: after  u-velocity         
496         & zub_adout,    & ! Adjoint output: before u-velocity
497         & zvb_adout,    & ! Adjoint output: before u-velocity         
498         & znu,          & ! 3D random field for u
499         & znv,          & ! 3D random field for v
500         & zbu,          & ! 3D random field for u
501         & zbv,          & ! 3D random field for v
502         & zau,          & ! 3D random field for u
503         & zav             ! 3D random field for v
504         
505      REAL(KIND=wp) :: &
506         & zsp1,         & ! scalar product involving the tangent routine
507         & zsp1_1,       & !   scalar product components
508         & zsp1_2,       & 
509         & zsp1_3,       & 
510         & zsp1_4,       & 
511         & zsp1_5,       & 
512         & zsp1_6,       & 
513         & zsp2,         & ! scalar product involving the adjoint routine
514         & zsp2_1,       & !   scalar product components
515         & zsp2_2,       & 
516         & zsp2_3,       & 
517         & zsp2_4,       & 
518         & zsp2_5,       & 
519         & zsp2_6
520      CHARACTER(LEN=14) :: cl_name
521
522      ! Allocate memory
523
524      ALLOCATE( &
525         & zun_tlin(jpi,jpj,jpk),     & 
526         & zvn_tlin(jpi,jpj,jpk),     & 
527         & zua_tlin(jpi,jpj,jpk),     & 
528         & zva_tlin(jpi,jpj,jpk),     & 
529         & zub_tlin(jpi,jpj,jpk),     & 
530         & zvb_tlin(jpi,jpj,jpk),     & 
531         & zun_adin(jpi,jpj,jpk),     & 
532         & zvn_adin(jpi,jpj,jpk),     & 
533         & zua_adin(jpi,jpj,jpk),     & 
534         & zva_adin(jpi,jpj,jpk),     & 
535         & zub_adin(jpi,jpj,jpk),     & 
536         & zvb_adin(jpi,jpj,jpk),     & 
537         & zun_tlout(jpi,jpj,jpk),    & 
538         & zvn_tlout(jpi,jpj,jpk),    & 
539         & zua_tlout(jpi,jpj,jpk),    & 
540         & zva_tlout(jpi,jpj,jpk),    & 
541         & zub_tlout(jpi,jpj,jpk),    & 
542         & zvb_tlout(jpi,jpj,jpk),    & 
543         & zun_adout(jpi,jpj,jpk),    & 
544         & zvn_adout(jpi,jpj,jpk),    & 
545         & zua_adout(jpi,jpj,jpk),    & 
546         & zva_adout(jpi,jpj,jpk),    & 
547         & zub_adout(jpi,jpj,jpk),    & 
548         & zvb_adout(jpi,jpj,jpk),    & 
549         & znu(jpi,jpj,jpk),          & 
550         & znv(jpi,jpj,jpk),          & 
551         & zbu(jpi,jpj,jpk),          & 
552         & zbv(jpi,jpj,jpk),          & 
553         & zau(jpi,jpj,jpk),          & 
554         & zav(jpi,jpj,jpk)           & 
555         & )
556
557
558      !==================================================================
559      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
560      !    dy = ( hdivb_tl, hdivn_tl )
561      !==================================================================
562
563      !--------------------------------------------------------------------
564      ! Reset the tangent and adjoint variables
565      !--------------------------------------------------------------------
566
567          zun_tlin(:,:,:) = 0.0_wp     
568          zvn_tlin(:,:,:) = 0.0_wp     
569          zua_tlin(:,:,:) = 0.0_wp     
570          zva_tlin(:,:,:) = 0.0_wp     
571          zub_tlin(:,:,:) = 0.0_wp     
572          zvb_tlin(:,:,:) = 0.0_wp     
573          zun_adin(:,:,:) = 0.0_wp     
574          zvn_adin(:,:,:) = 0.0_wp     
575          zua_adin(:,:,:) = 0.0_wp     
576          zva_adin(:,:,:) = 0.0_wp     
577          zub_adin(:,:,:) = 0.0_wp     
578          zvb_adin(:,:,:) = 0.0_wp     
579          zun_tlout(:,:,:) = 0.0_wp     
580          zvn_tlout(:,:,:) = 0.0_wp     
581          zua_tlout(:,:,:) = 0.0_wp     
582          zva_tlout(:,:,:) = 0.0_wp     
583          zub_tlout(:,:,:) = 0.0_wp     
584          zvb_tlout(:,:,:) = 0.0_wp     
585          zun_adout(:,:,:) = 0.0_wp     
586          zvn_adout(:,:,:) = 0.0_wp     
587          zua_adout(:,:,:) = 0.0_wp     
588          zva_adout(:,:,:) = 0.0_wp     
589          zub_adout(:,:,:) = 0.0_wp     
590          zvb_adout(:,:,:) = 0.0_wp     
591          znu(:,:,:) = 0.0_wp           
592          znv(:,:,:) = 0.0_wp           
593          zbu(:,:,:) = 0.0_wp           
594          zbv(:,:,:) = 0.0_wp           
595          zau(:,:,:) = 0.0_wp           
596          zav(:,:,:) = 0.0_wp           
597
598          un_tl(:,:,:) = 0.0_wp     
599          vn_tl(:,:,:) = 0.0_wp     
600          ua_tl(:,:,:) = 0.0_wp     
601          va_tl(:,:,:) = 0.0_wp     
602          ub_tl(:,:,:) = 0.0_wp     
603          vb_tl(:,:,:) = 0.0_wp     
604          un_ad(:,:,:) = 0.0_wp     
605          vn_ad(:,:,:) = 0.0_wp     
606          ua_ad(:,:,:) = 0.0_wp     
607          va_ad(:,:,:) = 0.0_wp     
608          ub_ad(:,:,:) = 0.0_wp     
609          vb_ad(:,:,:) = 0.0_wp     
610
611
612      !--------------------------------------------------------------------
613      ! Initialize the tangent input with random noise: dx
614      !--------------------------------------------------------------------
615
616      DO jj = 1, jpj
617         DO ji = 1, jpi
618            iseed_2d(ji,jj) = - ( 596035 + &
619               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
620         END DO
621      END DO
622      CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu )
623
624      DO jj = 1, jpj
625         DO ji = 1, jpi
626            iseed_2d(ji,jj) = - ( 523432 + &
627               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
628         END DO
629      END DO
630      CALL grid_random( iseed_2d, znv, 'V', 0.0_wp, stdv )
631
632      DO jj = 1, jpj
633         DO ji = 1, jpi
634            iseed_2d(ji,jj) = - ( 456953 + &
635               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
636         END DO
637      END DO
638      CALL grid_random( iseed_2d, zbu, 'U', 0.0_wp, stdu )
639
640      DO jj = 1, jpj
641         DO ji = 1, jpi
642            iseed_2d(ji,jj) = - ( 267406 + &
643               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
644         END DO
645      END DO
646      CALL grid_random( iseed_2d, zbv, 'V', 0.0_wp, stdv )
647
648      DO jj = 1, jpj
649         DO ji = 1, jpi
650            iseed_2d(ji,jj) = - ( 432545 + &
651               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
652         END DO
653      END DO
654      CALL grid_random( iseed_2d, zau, 'U', 0.0_wp, stdu )
655
656      DO jj = 1, jpj
657         DO ji = 1, jpi
658            iseed_2d(ji,jj) = - ( 287503 + &
659               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
660         END DO
661      END DO
662      CALL grid_random( iseed_2d, zav, 'V', 0.0_wp, stdv )
663
664      DO jk = 1, jpk
665         DO jj = nldj, nlej
666            DO ji = nldi, nlei
667               zun_tlin(ji,jj,jk) = znu(ji,jj,jk) 
668               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk) 
669               zub_tlin(ji,jj,jk) = zbu(ji,jj,jk)
670               zvb_tlin(ji,jj,jk) = zbv(ji,jj,jk)
671               zua_tlin(ji,jj,jk) = zau(ji,jj,jk)
672               zva_tlin(ji,jj,jk) = zav(ji,jj,jk)
673            END DO
674         END DO
675      END DO
676
677      un_tl(:,:,:) = zun_tlin(:,:,:)
678      vn_tl(:,:,:) = zvn_tlin(:,:,:)
679      ub_tl(:,:,:) = zub_tlin(:,:,:)
680      vb_tl(:,:,:) = zvb_tlin(:,:,:)
681      ua_tl(:,:,:) = zua_tlin(:,:,:)
682      va_tl(:,:,:) = zva_tlin(:,:,:)
683
684      call dyn_nxt_tan ( nit000 )
685
686      zun_tlout(:,:,:) = un_tl(:,:,:)
687      zvn_tlout(:,:,:) = vn_tl(:,:,:)
688      zub_tlout(:,:,:) = ub_tl(:,:,:)
689      zvb_tlout(:,:,:) = vb_tl(:,:,:)
690      zua_tlout(:,:,:) = ua_tl(:,:,:)
691      zva_tlout(:,:,:) = va_tl(:,:,:)
692
693      !--------------------------------------------------------------------
694      ! Initialize the adjoint variables: dy^* = W dy
695      !--------------------------------------------------------------------
696
697      DO jk = 1, jpk
698        DO jj = nldj, nlej
699           DO ji = nldi, nlei
700              zun_adin(ji,jj,jk) = zun_tlout(ji,jj,jk) &
701                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
702                 &               * umask(ji,jj,jk)
703              zvn_adin(ji,jj,jk) = zvn_tlout(ji,jj,jk) &
704                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
705                 &               * vmask(ji,jj,jk)
706              zub_adin(ji,jj,jk) = zub_tlout(ji,jj,jk) &
707                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
708                 &               * umask(ji,jj,jk)
709              zvb_adin(ji,jj,jk) = zvb_tlout(ji,jj,jk) &
710                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
711                 &               * vmask(ji,jj,jk)
712              zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
713                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
714                 &               * umask(ji,jj,jk)
715              zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
716                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
717                 &               * vmask(ji,jj,jk)
718            END DO
719         END DO
720      END DO
721      !--------------------------------------------------------------------
722      ! Compute the scalar product: ( L dx )^T W dy
723      !--------------------------------------------------------------------
724
725      zsp1_1 = DOT_PRODUCT( zun_tlout, zun_adin )
726      zsp1_2 = DOT_PRODUCT( zvn_tlout, zvn_adin )
727      zsp1_3 = DOT_PRODUCT( zub_tlout, zub_adin )
728      zsp1_4 = DOT_PRODUCT( zvb_tlout, zvb_adin )
729      zsp1_5 = DOT_PRODUCT( zua_tlout, zua_adin )
730      zsp1_6 = DOT_PRODUCT( zva_tlout, zva_adin )
731      zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6 
732
733      !--------------------------------------------------------------------
734      ! Call the adjoint routine: dx^* = L^T dy^*
735      !--------------------------------------------------------------------
736
737      un_ad(:,:,:) = zun_adin(:,:,:)
738      vn_ad(:,:,:) = zvn_adin(:,:,:)
739      ub_ad(:,:,:) = zub_adin(:,:,:)
740      vb_ad(:,:,:) = zvb_adin(:,:,:)
741      ua_ad(:,:,:) = zua_adin(:,:,:)
742      va_ad(:,:,:) = zva_adin(:,:,:)
743
744      CALL dyn_nxt_adj ( nit000 )
745
746      zun_adout(:,:,:) = un_ad(:,:,:)
747      zvn_adout(:,:,:) = vn_ad(:,:,:)
748      zub_adout(:,:,:) = ub_ad(:,:,:)
749      zvb_adout(:,:,:) = vb_ad(:,:,:)
750      zua_adout(:,:,:) = ua_ad(:,:,:)
751      zva_adout(:,:,:) = va_ad(:,:,:)
752
753      zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
754      zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
755      zsp2_3 = DOT_PRODUCT( zub_tlin, zub_adout )
756      zsp2_4 = DOT_PRODUCT( zvb_tlin, zvb_adout )
757      zsp2_5 = DOT_PRODUCT( zua_tlin, zua_adout )
758      zsp2_6 = DOT_PRODUCT( zva_tlin, zva_adout )
759      zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6
760
761      ! Compare the scalar products
762      ! 14 char:'12345678901234'
763      cl_name = 'dyn_nxt_adj   '
764      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
765
766      DEALLOCATE( &
767         & zun_tlin,     & 
768         & zvn_tlin,     & 
769         & zua_tlin,     & 
770         & zva_tlin,     & 
771         & zub_tlin,     & 
772         & zvb_tlin,     & 
773         & zun_adin,     & 
774         & zvn_adin,     & 
775         & zua_adin,     & 
776         & zva_adin,     & 
777         & zub_adin,     & 
778         & zvb_adin,     & 
779         & zun_tlout,    & 
780         & zvn_tlout,    & 
781         & zua_tlout,    & 
782         & zva_tlout,    & 
783         & zub_tlout,    & 
784         & zvb_tlout,    & 
785         & zun_adout,    & 
786         & zvn_adout,    & 
787         & zua_adout,    & 
788         & zva_adout,    & 
789         & zub_adout,    & 
790         & zvb_adout,    & 
791         & znu,          & 
792         & znv,          & 
793         & zbu,          & 
794         & zbv,          & 
795         & zau,          & 
796         & zav           & 
797         & )
798
799
800   END SUBROUTINE dyn_nxt_adj_tst
801   !!======================================================================
802#endif
803END MODULE dynnxt_tam
Note: See TracBrowser for help on using the repository browser.