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/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/DYN – NEMO

source: branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/DYN/dynnxt_tam.F90 @ 4555

Last change on this file since 4555 was 2578, checked in by rblod, 13 years ago

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