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

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

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

add TAM sources

File size: 42.8 KB
Line 
1MODULE dynvor_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  dynvor_tam  ***
5   !! Ocean dynamics: Update the momentum trend with the relative and
6   !!                 planetary vorticity trends
7   !!                 Tangent and Adjoint Module
8   !!======================================================================
9   !! History of the drect module:
10   !!            1.0  !  89-12  (P. Andrich)  vor_ens: Original code
11   !!            5.0  !  91-11  (G. Madec) vor_ene, vor_mix: Original code
12   !!            6.0  !  96-01  (G. Madec)  s-coord, suppress work arrays
13   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
14   !!            8.5  !  04-02  (G. Madec)  vor_een: Original code
15   !!            9.0  !  03-08  (G. Madec)  vor_ctl: Original code
16   !!            9.0  !  05-11  (G. Madec)  dyn_vor: Original code (new step architecture)
17   !!            9.0  !  06-11  (G. Madec)  flux form advection: add metric term
18   !! History of the TAM module:
19   !!            9.0  !  08-06  (A. Vidard) Skeleton
20   !!            9.0  !  09-01  (A. Vidard) TAM of the 06-11 version
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dyn_vor     : Update the momentum trend with the vorticity trend
25   !!       vor_ens : enstrophy conserving scheme       (ln_dynvor_ens=T)
26   !!       vor_ene : energy conserving scheme          (ln_dynvor_ene=T)
27   !!       vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
28   !!       vor_een : energy and enstrophy conserving   (ln_dynvor_een=T)
29   !!       vor_ctl : set and control of the different vorticity option
30   !!----------------------------------------------------------------------
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      & jpim1,               &
38      & jpjm1,               & 
39      & jpkm1,               & 
40      & jpiglo
41   USE oce           , ONLY: & ! ocean dynamics and tracers
42      & un,                  &
43      & vn,                  &
44      & rotn
45   USE oce_tam       , ONLY: &
46      & un_tl,         &
47      & vn_tl,         &
48      & ua_tl,         &
49      & va_tl,         &
50      & rotn_tl,       &
51      & un_ad,         &
52      & vn_ad,         &
53      & ua_ad,         &
54      & va_ad,         &
55      & rotn_ad
56   USE divcur        , ONLY: &
57      & div_cur
58   USE divcur_tam    , ONLY: &
59      & div_cur_tan
60   USE dom_oce       , ONLY: & ! ocean space and time domain
61      & ln_sco,              &
62      & ff,                  &
63      & e1u,                 &
64      & e2u,                 &
65      & e1v,                 &
66      & e2v,                 &
67#if defined key_zco
68      & e3t_0,               &
69#else
70      & e3u,                 &
71      & e3v,                 &
72      & e3f,                 &
73#endif
74      & e1f,                 &
75      & e2f,                 &
76      & mig,                 &
77      & mjg,                 &
78      & nldi,                &
79      & nldj,                &
80      & nlei,                &
81      & nlej,                &
82      & umask,               &
83      & vmask
84   USE dynadv        , ONLY: &
85      & ln_dynadv_vec ! vector form flag
86   USE in_out_manager, ONLY: & ! I/O manager
87      & ctl_stop,            &
88      & lk_esopa,            &
89      & numnam,              &
90      & numout,              &
91      & nit000,              &
92      & nitend,              &
93      & lwp
94   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
95      & grid_random
96   USE dotprodfld,     ONLY: & ! Computes dot product for 3D and 2D fields
97      & dot_product
98   USE tstool_tam    , ONLY: &
99      & prntst_adj,          & !
100                               ! random field standard deviation for:
101      & stdu,                & !   u-velocity
102      & stdv                   !   v-velocity
103
104   IMPLICIT NONE
105   PRIVATE
106
107   PUBLIC   dyn_vor_tan     ! routine called by step_tam.F90
108   PUBLIC   dyn_vor_adj     ! routine called by step_tam.F90
109   PUBLIC   dyn_vor_adj_tst ! routine called by the tst.F90
110
111
112   !!* Namelist nam_dynvor: vorticity term
113   LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
114   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
115   LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
116   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.   !: energy and enstrophy conserving scheme
117
118   INTEGER ::   nvor = 0   ! type of vorticity trend used
119   INTEGER ::   ncor = 1   ! coriolis
120   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
121   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
122
123   !! * Substitutions
124#  include "domzgr_substitute.h90"
125#  include "vectopt_loop_substitute.h90"
126
127CONTAINS
128
129   SUBROUTINE dyn_vor_tan( kt )
130      !!----------------------------------------------------------------------
131      !!
132      !! ** Purpose of the direct routine:
133      !!               compute the lateral ocean tracer physics.
134      !!
135      !! ** Action : - Update (ua,va) with the now vorticity term trend
136      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
137      !!               and planetary vorticity trends) ('key_trddyn')
138      !!----------------------------------------------------------------------
139      !!
140      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
141      !!----------------------------------------------------------------------
142      IF( kt == nit000 )   CALL vor_ctl_tam          ! initialisation & control of options
143
144      !                                          ! vorticity term
145      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
146      !
147      CASE ( -1 )                                      ! esopa: test all possibility with control print
148!         CALL vor_ene_tan( kt, ntot, ua_tl, va_tl )
149         CALL vor_ens_tan( kt, ntot, ua_tl, va_tl )
150!         CALL vor_mix_tan( kt )
151!         CALL vor_een_tan( kt, ntot, ua_tl, va_tl )
152         !
153      CASE ( 0 )                                       ! energy conserving scheme
154         CALL ctl_stop ('vor_ene_tan not available yet')
155!         CALL vor_ene_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
156         !
157      CASE ( 1 )                                       ! enstrophy conserving scheme
158         CALL vor_ens_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
159         !
160      CASE ( 2 )                                       ! mixed ene-ens scheme
161         CALL ctl_stop ('vor_mix_tan not available yet')
162!         CALL vor_mix_tan( kt )                               ! total vorticity (mix=ens-ene)
163         !
164      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
165         CALL ctl_stop ('vor_een_tan not available yet')
166!         CALL vor_een_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
167         !
168      END SELECT
169
170   END SUBROUTINE dyn_vor_tan
171   SUBROUTINE vor_ens_tan( kt, kvor, pua_tl, pva_tl )
172      !!----------------------------------------------------------------------
173      !!                ***  ROUTINE vor_ens_tan  ***
174      !!
175      !! ** Purpose of the direct routine:
176      !!      Compute the now total vorticity trend and add it to
177      !!      the general trend of the momentum equation.
178      !!
179      !! ** Method of the direct routine:
180      !!      Trend evaluated using now fields (centered in time)
181      !!      and the Sadourny (1975) flux FORM formulation : conserves the
182      !!      potential enstrophy of a horizontally non-divergent flow. the
183      !!      trend of the vorticity term is given by:
184      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
185      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
186      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
187      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
188      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
189      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
190      !!      Add this trend to the general momentum trend (ua,va):
191      !!          (ua,va) = (ua,va) + ( voru , vorv )
192      !!
193      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
194      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
195      !!               and planetary vorticity trends) ('key_trddyn')
196      !!
197      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
198      !!----------------------------------------------------------------------
199      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
200      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
201         !                                                        ! =nrvm (relative vorticity or metric)
202      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_tl    ! total u-trend
203      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_tl    ! total v-trend
204      !!
205      INTEGER  ::   ji, jj, jk           ! dummy loop indices
206      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
207      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz         ! temporary 3D workspace
208      REAL(wp) ::   zuavtl, zvautl   ! temporary scalars
209      REAL(wp), DIMENSION(jpi,jpj) ::   zwxtl, zwytl, zwztl   ! temporary 3D workspace
210      !!----------------------------------------------------------------------
211     
212      IF( kt == nit000 ) THEN
213         IF(lwp) WRITE(numout,*)
214         IF(lwp) WRITE(numout,*) 'dyn_vor_ens_tan : vorticity term: enstrophy conserving scheme'
215         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
216      ENDIF
217
218      ! Local constant initialization
219      zfact1 = 0.5 * 0.25
220
221!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
222      !                                                ! ===============
223      DO jk = 1, jpkm1                                 ! Horizontal slab
224         !                                             ! ===============
225         ! Potential vorticity and horizontal fluxes
226         ! -----------------------------------------
227         SELECT CASE( kvor )      ! vorticity considered
228         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
229         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
230         CASE ( 3 )                                                ! metric term
231            DO jj = 1, jpjm1
232               DO ji = 1, fs_jpim1   ! vector opt.
233                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
234                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
235                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
236               END DO
237            END DO
238         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
239         CASE ( 5 )                                                ! total (coriolis + metric)
240            DO jj = 1, jpjm1
241               DO ji = 1, fs_jpim1   ! vector opt.
242                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
243                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
244                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
245                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
246                       &       )
247               END DO
248            END DO
249         END SELECT
250
251         IF( ln_sco ) THEN
252            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
253               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
254                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
255                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
256                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
257               END DO
258            END DO
259         ELSE
260            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
261               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
262                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
263                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
264               END DO
265            END DO
266         ENDIF
267
268         ! Compute and add the vorticity term trend
269         ! ----------------------------------------
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1   ! vector opt.
272               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
273                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
274               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
275                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
276            END DO 
277         END DO 
278         !                                             ! ===============
279      END DO                                           !   End of slab
280      !                                                ! ===============
281
282!CDIR PARALLEL DO PRIVATE( zwxtl, zwytl, zwztl )
283      ! ===================
284      ! Tangent counterpart
285      ! ===================
286      !                                                ! ===============
287      DO jk = 1, jpkm1                                 ! Horizontal slab
288         !                                             ! ===============
289         ! Potential vorticity and horizontal fluxes
290         ! -----------------------------------------
291         SELECT CASE( kvor )      ! vorticity considered
292         CASE ( 1 )   ;   zwztl(:,:) =  0.0_wp      ! planetary vorticity (Coriolis)
293         CASE ( 2 ,4)   ;   zwztl(:,:) =   rotn_tl(:,:,jk)                ! relative  vorticity
294         CASE ( 3 ,5 )                                                ! metric term
295            DO jj = 1, jpjm1
296               DO ji = 1, fs_jpim1   ! vector opt.
297                  zwztl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )   &
298                       &           - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) ) ) &
299                       &        * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
300               END DO
301            END DO
302         END SELECT
303
304         IF( ln_sco ) THEN
305            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
306               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
307                  zwztl(ji,jj) = zwztl(ji,jj) / fse3f(ji,jj,jk)
308                  zwxtl(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un_tl(ji,jj,jk)
309                  zwytl(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn_tl(ji,jj,jk)
310               END DO
311            END DO
312         ELSE
313            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
314               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
315                  zwxtl(ji,jj) = e2u(ji,jj) * un_tl(ji,jj,jk)
316                  zwytl(ji,jj) = e1v(ji,jj) * vn_tl(ji,jj,jk)
317               END DO
318            END DO
319         ENDIF
320
321         ! Compute and add the vorticity term trend
322         ! ----------------------------------------
323         DO jj = 2, jpjm1
324            DO ji = fs_2, fs_jpim1   ! vector opt.
325               zuavtl = zfact1 / e1u(ji,jj) * ( zwytl(ji  ,jj-1) + zwytl(ji+1,jj-1)   &
326                  &                         + zwytl(ji  ,jj  ) + zwytl(ji+1,jj  ) )
327               zvautl =-zfact1 / e2v(ji,jj) * ( zwxtl(ji-1,jj  ) + zwxtl(ji-1,jj+1)   &
328                  &                         + zwxtl(ji  ,jj  ) + zwxtl(ji  ,jj+1) )
329               pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk)                   &
330                  &     + zuavtl * ( zwz(  ji,jj-1) + zwz(  ji,jj) ) &
331                  &     + zuav   * ( zwztl(ji,jj-1) + zwztl(ji,jj) )
332               pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk)                   &
333                  &     + zvautl * ( zwz(  ji-1,jj) + zwz(  ji,jj) ) &
334                  &     + zvau   * ( zwztl(ji-1,jj) + zwztl(ji,jj) )
335            END DO 
336         END DO 
337         !                                             ! ===============
338      END DO                                           !   End of slab
339      !                                                ! ===============
340   END SUBROUTINE vor_ens_tan
341   
342   SUBROUTINE dyn_vor_adj( kt )
343      !!----------------------------------------------------------------------
344      !!
345      !! ** Purpose of the direct routine:
346      !!               compute the lateral ocean tracer physics.
347      !!
348      !! ** Action : - Update (ua,va) with the now vorticity term trend
349      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
350      !!               and planetary vorticity trends) ('key_trddyn')
351      !!----------------------------------------------------------------------
352      !!
353      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
354      !!----------------------------------------------------------------------
355
356      IF( kt == nitend )   CALL vor_ctl_tam          ! initialisation & control of options
357
358      !                                          ! vorticity term
359      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
360      !
361      CASE ( -1 )                                      ! esopa: test all possibility with control print
362!         CALL vor_een_adj( kt, ntot, ua_ad, va_ad )
363!         CALL vor_mix_adj( kt )
364         CALL vor_ens_adj( kt, ntot, ua_ad, va_ad )
365!         CALL vor_ene_adj( kt, ntot, ua_ad, va_ad )
366         !
367      CASE ( 0 )                                       ! energy conserving scheme
368         CALL ctl_stop ('vor_ene_adj not available yet')
369!         CALL vor_ene_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
370         !
371      CASE ( 1 )                                       ! enstrophy conserving scheme
372         CALL vor_ens_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
373         !
374      CASE ( 2 )                                       ! mixed ene-ens scheme
375         CALL ctl_stop ('vor_mix_adj not available yet')
376!         CALL vor_mix_adj( kt )                               ! total vorticity (mix=ens-ene)
377         !
378      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
379         CALL ctl_stop ('vor_een_adj not available yet')
380!         CALL vor_een_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
381         !
382      END SELECT
383   END SUBROUTINE dyn_vor_adj
384   SUBROUTINE vor_ens_adj( kt, kvor, pua_ad, pva_ad )
385      !!----------------------------------------------------------------------
386      !!                ***  ROUTINE vor_ens_adj  ***
387      !!
388      !! ** Purpose of the direct routine:
389      !!      Compute the now total vorticity trend and add it to
390      !!      the general trend of the momentum equation.
391      !!
392      !! ** Method of the direct routine:
393      !!      Trend evaluated using now fields (centered in time)
394      !!      and the Sadourny (1975) flux FORM formulation : conserves the
395      !!      potential enstrophy of a horizontally non-divergent flow. the
396      !!      trend of the vorticity term is given by:
397      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
398      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
399      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
400      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
401      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
402      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
403      !!      Add this trend to the general momentum trend (ua,va):
404      !!          (ua,va) = (ua,va) + ( voru , vorv )
405      !!
406      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
407      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
408      !!               and planetary vorticity trends) ('key_trddyn')
409      !!
410      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
411      !!----------------------------------------------------------------------
412      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
413      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
414         !                                                        ! =nrvm (relative vorticity or metric)
415      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_ad    ! total u-trend
416      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_ad    ! total v-trend
417      !!
418      INTEGER  ::   ji, jj, jk           ! dummy loop indices
419      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
420      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz         ! temporary 3D workspace
421      REAL(wp) ::   zuavad, zvauad   ! temporary scalars
422      REAL(wp), DIMENSION(jpi,jpj) ::   zwxad, zwyad, zwzad   ! temporary 3D workspace
423      !!----------------------------------------------------------------------
424     
425      IF( kt == nitend ) THEN
426         IF(lwp) WRITE(numout,*)
427         IF(lwp) WRITE(numout,*) 'dyn_vor_ens_adj : vorticity term: enstrophy conserving scheme'
428         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
429      ENDIF
430
431      ! Local constant initialization
432      zfact1 = 0.5 * 0.25
433
434!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
435      !                                                ! ===============
436      DO jk = 1, jpkm1                                 ! Horizontal slab
437         !                                             ! ===============
438         ! Potential vorticity and horizontal fluxes
439         ! -----------------------------------------
440         SELECT CASE( kvor )      ! vorticity considered
441         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
442         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
443         CASE ( 3 )                                                ! metric term
444            DO jj = 1, jpjm1
445               DO ji = 1, fs_jpim1   ! vector opt.
446                  zwz(ji,jj) = ( ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )  &
447                       &       - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) ) )&
448                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
449               END DO
450            END DO
451         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
452         CASE ( 5 )                                                ! total (coriolis + metric)
453            DO jj = 1, jpjm1
454               DO ji = 1, fs_jpim1   ! vector opt.
455                  zwz(ji,jj) = ( ff (ji,jj)                         &
456                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
457                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
458                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )  &
459                       &       )
460               END DO
461            END DO
462         END SELECT
463
464         IF( ln_sco ) THEN
465            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
466               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
467                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
468                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
469                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
470               END DO
471            END DO
472         ELSE
473            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
474               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
475                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
476                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
477               END DO
478            END DO
479         ENDIF
480
481         ! Compute and add the vorticity term trend
482         ! ----------------------------------------
483         DO jj = 2, jpjm1
484            DO ji = fs_2, fs_jpim1   ! vector opt.
485               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
486                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
487               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
488                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
489            END DO 
490         END DO 
491         !                                             ! ===============
492      END DO                                           !   End of slab
493      !                                                ! ===============
494!CDIR PARALLEL DO PRIVATE( zwxad, zwyad, zwzad )
495      ! ===================
496      ! Adjoint counterpart
497      ! ===================
498      zuavad     = 0.0_wp
499      zvauad     = 0.0_wp
500      zwxad(:,:) = 0.0_wp
501      zwyad(:,:) = 0.0_wp
502      zwzad(:,:) = 0.0_wp
503      !                                                ! ===============
504      DO jk = jpkm1, 1, -1                             ! Horizontal slab
505         !                                             ! ===============
506         ! Compute and add the vorticity term trend
507         ! ----------------------------------------
508         DO jj = jpjm1, 2, -1
509            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
510               zuavad = zuavad + pua_ad(ji,jj,jk) * ( zwz(ji,jj-1) + zwz(ji,jj) ) 
511               zwzad(ji,jj-1) = zwzad(ji,jj-1) + pua_ad(ji,jj,jk) * zuav
512               zwzad(ji,jj  ) = zwzad(ji,jj  ) + pua_ad(ji,jj,jk) * zuav
513               
514               zvauad = zvauad + pva_ad(ji,jj,jk) * ( zwz(ji-1,jj) + zwz(ji,jj) )
515               zwzad(ji-1,jj) = zwzad(ji-1,jj) + pva_ad(ji,jj,jk) * zvau
516               zwzad(ji  ,jj) = zwzad(ji  ,jj) + pva_ad(ji,jj,jk) * zvau
517               
518               zwyad(ji  ,jj-1) = zwyad(ji  ,jj-1) + zuavad * zfact1 / e1u(ji,jj)
519               zwyad(ji+1,jj-1) = zwyad(ji+1,jj-1) + zuavad * zfact1 / e1u(ji,jj)
520               zwyad(ji  ,jj  ) = zwyad(ji  ,jj  ) + zuavad * zfact1 / e1u(ji,jj)
521               zwyad(ji+1,jj  ) = zwyad(ji+1,jj  ) + zuavad * zfact1 / e1u(ji,jj)
522               zuavad = 0.0_wp
523
524               zwxad(ji-1,jj  ) = zwxad(ji-1,jj  ) - zvauad * zfact1 / e2v(ji,jj)
525               zwxad(ji-1,jj+1) = zwxad(ji-1,jj+1) - zvauad * zfact1 / e2v(ji,jj)
526               zwxad(ji  ,jj  ) = zwxad(ji  ,jj  ) - zvauad * zfact1 / e2v(ji,jj)
527               zwxad(ji  ,jj+1) = zwxad(ji  ,jj+1) - zvauad * zfact1 / e2v(ji,jj)
528               zvauad = 0.0_wp
529            END DO 
530         END DO
531         IF( ln_sco ) THEN
532            DO jj = jpj, 1, -1     ! caution: don't use (:,:) for this loop
533               DO ji = jpi, 1, -1  ! it causes optimization problems on NEC in auto-tasking
534                  zwzad(ji,jj) = zwzad(ji,jj) / fse3f(ji,jj,jk)
535                  un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + zwxad(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
536                  vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + zwyad(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,jk)
537                  zwxad(ji,jj) = 0.0_wp
538                  zwyad(ji,jj) = 0.0_wp
539               END DO
540            END DO
541         ELSE
542            DO jj = jpj, 1, -1    ! caution: don't use (:,:) for this loop
543               DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking
544                  un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e2u(ji,jj) * zwxad(ji,jj)
545                  vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e1v(ji,jj) * zwyad(ji,jj)
546                  zwxad(ji,jj) = 0.0_wp
547                  zwyad(ji,jj) = 0.0_wp
548               END DO
549            END DO
550         ENDIF
551         ! Potential vorticity and horizontal fluxes
552         ! -----------------------------------------
553         SELECT CASE( kvor )      ! vorticity considered
554         CASE ( 1 )                     ! planetary vorticity (Coriolis)
555            zwzad(:,:) =  0.0_wp
556         CASE ( 2 ,4)                   ! relative  vorticity
557            rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:)
558            zwzad(:,:) = 0.0_wp
559         CASE ( 3 ,5 )                  ! metric term
560            DO jj = jpjm1, 1, -1
561               DO ji = fs_jpim1, 1, -1  ! vector opt.
562                  vn_ad(ji+1,jj,jk) = vn_ad(ji+1,jj,jk)                   &
563                     &  + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) )    &
564                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
565                  vn_ad(ji  ,jj,jk) = vn_ad(ji  ,jj,jk)                   &
566                     &  + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) )    &
567                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
568                  un_ad(ji,jj+1,jk) = un_ad(ji,jj+1,jk)                   &
569                     &  - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) )    &
570                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
571                  un_ad(ji,jj  ,jk) = un_ad(ji,jj  ,jk)                   &
572                     &  - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) )    &
573                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
574                  zwzad(ji,jj) = 0.0_wp
575               END DO
576            END DO
577         END SELECT
578         !                                             ! ===============
579      END DO                                           !   End of slab
580      !                                                ! ===============
581   END SUBROUTINE vor_ens_adj
582
583   SUBROUTINE vor_ctl_tam
584      !!---------------------------------------------------------------------
585      !!                  ***  ROUTINE vor_ctl_tam  ***
586      !!
587      !! ** Purpose :   Control the consistency between cpp options for
588      !!      tracer advection schemes
589      !!----------------------------------------------------------------------
590      INTEGER ::   ioptio          ! temporary integer
591      NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
592      !!----------------------------------------------------------------------
593
594      REWIND ( numnam )               ! Read Namelist nam_dynvor : Vorticity scheme options
595      READ   ( numnam, nam_dynvor )
596
597      IF(lwp) THEN                    ! Namelist print
598         WRITE(numout,*)
599         WRITE(numout,*) 'dyn:vor_ctl_tam : vorticity term : read namelist and control the consistency'
600         WRITE(numout,*) '~~~~~~~~~~~~~~~'
601         WRITE(numout,*) '        Namelist nam_dynvor : oice of the vorticity term scheme'
602         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
603         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
604         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
605         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
606      ENDIF
607
608      ioptio = 0                     ! Control of vorticity scheme options
609      IF( ln_dynvor_ene )   ioptio = ioptio + 1
610      IF( ln_dynvor_ens )   ioptio = ioptio + 1
611      IF( ln_dynvor_mix )   ioptio = ioptio + 1
612      IF( ln_dynvor_een )   ioptio = ioptio + 1
613      IF( lk_esopa      )   ioptio =          1
614
615      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
616
617      !                              ! Set nvor (type of scheme for vorticity)
618      IF( ln_dynvor_ene )   nvor =  0
619      IF( ln_dynvor_ens )   nvor =  1
620      IF( ln_dynvor_mix )   nvor =  2
621      IF( ln_dynvor_een )   nvor =  3
622      IF( lk_esopa      )   nvor = -1
623     
624      !                              ! Set ncor, nrvm, ntot (type of vorticity)
625      IF(lwp) WRITE(numout,*)
626      ncor = 1
627      IF( ln_dynadv_vec ) THEN     
628         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
629         nrvm = 2
630         ntot = 4
631      ELSE                       
632         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
633         nrvm = 3
634         ntot = 5
635      ENDIF
636     
637      IF(lwp) THEN                   ! Print the choice
638         WRITE(numout,*)
639         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
640         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
641         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
642         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
643         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
644      ENDIF
645      !
646   END SUBROUTINE vor_ctl_tam
647
648   SUBROUTINE dyn_vor_adj_tst( kumadt )
649      !!-----------------------------------------------------------------------
650      !!
651      !!                  ***  ROUTINE dyn_adv_adj_tst ***
652      !!
653      !! ** Purpose : Test the adjoint routine.
654      !!
655      !! ** Method  : Verify the scalar product
656      !!           
657      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
658      !!
659      !!              where  L   = tangent routine
660      !!                     L^T = adjoint routine
661      !!                     W   = diagonal matrix of scale factors
662      !!                     dx  = input perturbation (random field)
663      !!                     dy  = L dx
664      !!
665      !! ** Action  : Separate tests are applied for the following dx and dy:
666      !!               
667      !!              1) dx = ( SSH ) and dy = ( SSH )
668      !!                   
669      !! History :
670      !!        ! 08-08 (A. Vidard)
671      !!-----------------------------------------------------------------------
672      !! * Modules used
673
674      !! * Arguments
675      INTEGER, INTENT(IN) :: &
676         & kumadt             ! Output unit
677 
678      INTEGER :: &
679         & ji,    &        ! dummy loop indices
680         & jj,    &       
681         & jk     
682      INTEGER, DIMENSION(jpi,jpj) :: &
683         & iseed_2d        ! 2D seed for the random number generator
684
685      !! * Local declarations
686      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
687         & zun_tlin,     & ! Tangent input:  now   u-velocity
688         & zvn_tlin,     & ! Tangent input:  now   v-velocity
689         & zrotn_tlin,   & ! Tangent input:  now   rot
690         & zun_adout,    & ! Adjoint output: now   u-velocity
691         & zvn_adout,    & ! Adjoint output: now   v-velocity
692         & zrotn_adout,  & ! Adjoint output: now   rot
693         & zua_adout,    & ! Tangent output: after u-velocity
694         & zva_adout,    & ! Tangent output: after v-velocity
695         & zua_tlin,     & ! Tangent output: after u-velocity
696         & zva_tlin,     & ! Tangent output: after v-velocity
697         & zua_tlout,    & ! Tangent output: after u-velocity
698         & zva_tlout,    & ! Tangent output: after v-velocity
699         & zua_adin,     & ! Tangent output: after u-velocity
700         & zva_adin,     & ! Tangent output: after v-velocity
701         & zau,          & ! 3D random field for rotn
702         & zav,          & ! 3D random field for rotn
703         & znu,          & ! 3D random field for u
704         & znv             ! 3D random field for v
705      REAL(KIND=wp) :: &
706         & zsp1,         & ! scalar product involving the tangent routine
707         & zsp1_1,       & !   scalar product components
708         & zsp1_2,       & 
709         & zsp2,         & ! scalar product involving the adjoint routine
710         & zsp2_1,       & !   scalar product components
711         & zsp2_2,       & 
712         & zsp2_3,       & 
713         & zsp2_4,       & 
714         & zsp2_5
715      CHARACTER(LEN=14) :: cl_name
716
717      ! Allocate memory
718
719      ALLOCATE( &
720         & zun_tlin(jpi,jpj,jpk),     & 
721         & zvn_tlin(jpi,jpj,jpk),     & 
722         & zrotn_tlin(jpi,jpj,jpk),   & 
723         & zun_adout(jpi,jpj,jpk),    & 
724         & zvn_adout(jpi,jpj,jpk),    & 
725         & zrotn_adout(jpi,jpj,jpk),  & 
726         & zua_adout(jpi,jpj,jpk),    & 
727         & zva_adout(jpi,jpj,jpk),    & 
728         & zua_tlin(jpi,jpj,jpk),     & 
729         & zva_tlin(jpi,jpj,jpk),     & 
730         & zua_tlout(jpi,jpj,jpk),    & 
731         & zva_tlout(jpi,jpj,jpk),    & 
732         & zua_adin(jpi,jpj,jpk),     & 
733         & zva_adin(jpi,jpj,jpk),     & 
734         & zau(jpi,jpj,jpk),          & 
735         & zav(jpi,jpj,jpk),          & 
736         & znu(jpi,jpj,jpk),          & 
737         & znv(jpi,jpj,jpk)           &
738         & )
739
740      ! Initialize rotn
741      CALL div_cur ( nit000 )
742     
743      !==================================================================
744      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
745      !    dy = ( hdivb_tl, hdivn_tl )
746      !==================================================================
747
748      !--------------------------------------------------------------------
749      ! Reset the tangent and adjoint variables
750      !--------------------------------------------------------------------
751     
752      zun_tlin(:,:,:) = 0.0_wp     
753      zvn_tlin(:,:,:) = 0.0_wp     
754      zrotn_tlin(:,:,:) = 0.0_wp   
755      zun_adout(:,:,:) = 0.0_wp     
756      zvn_adout(:,:,:) = 0.0_wp     
757      zrotn_adout(:,:,:) = 0.0_wp   
758      zua_tlout(:,:,:) = 0.0_wp     
759      zva_tlout(:,:,:) = 0.0_wp     
760      zua_adin(:,:,:) = 0.0_wp     
761      zva_adin(:,:,:) = 0.0_wp     
762      zua_adout(:,:,:) = 0.0_wp     
763      zva_adout(:,:,:) = 0.0_wp     
764      zua_tlin(:,:,:) = 0.0_wp     
765      zva_tlin(:,:,:) = 0.0_wp     
766      znu(:,:,:) = 0.0_wp           
767      znv(:,:,:) = 0.0_wp   
768      zau(:,:,:) = 0.0_wp           
769      zav(:,:,:) = 0.0_wp   
770
771
772      un_tl(:,:,:) = 0.0_wp
773      vn_tl(:,:,:) = 0.0_wp
774      ua_tl(:,:,:) = 0.0_wp
775      va_tl(:,:,:) = 0.0_wp
776      un_ad(:,:,:) = 0.0_wp
777      vn_ad(:,:,:) = 0.0_wp
778      ua_ad(:,:,:) = 0.0_wp
779      va_ad(:,:,:) = 0.0_wp
780      rotn_tl(:,:,:) = 0.0_wp
781      rotn_ad(:,:,:) = 0.0_wp
782         
783      !--------------------------------------------------------------------
784      ! Initialize the tangent input with random noise: dx
785      !--------------------------------------------------------------------
786
787      DO jj = 1, jpj
788         DO ji = 1, jpi
789            iseed_2d(ji,jj) = - ( 596035 + &
790               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
791         END DO
792      END DO
793      CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu )
794
795      DO jj = 1, jpj
796         DO ji = 1, jpi
797            iseed_2d(ji,jj) = - ( 523432 + &
798               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
799         END DO
800      END DO
801      CALL grid_random( iseed_2d, znv, 'V', 0.0_wp, stdv )
802
803      DO jj = 1, jpj
804         DO ji = 1, jpi
805            iseed_2d(ji,jj) = - ( 432545 + &
806               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
807         END DO
808      END DO
809      CALL grid_random( iseed_2d, zau, 'U', 0.0_wp, stdu )
810
811      DO jj = 1, jpj
812         DO ji = 1, jpi
813            iseed_2d(ji,jj) = - ( 287503 + &
814               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
815         END DO
816      END DO
817      CALL grid_random( iseed_2d, zav, 'V', 0.0_wp, stdv )
818!zun_tlin(:,:,:) = znu(:,:,:)
819!zvn_tlin(:,:,:) = znv(:,:,:)
820!zua_tlin(:,:,:) = zau(:,:,:)
821!zva_tlin(:,:,:) = zav(:,:,:)
822      DO jk = 1, jpk
823         DO jj = nldj, nlej
824            DO ji = nldi, nlei
825               zun_tlin(ji,jj,jk) = znu(ji,jj,jk) 
826               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk) 
827               zua_tlin(ji,jj,jk) = zau(ji,jj,jk) 
828               zva_tlin(ji,jj,jk) = zav(ji,jj,jk)
829            END DO
830         END DO
831      END DO
832      un_tl(:,:,:) = zun_tlin(:,:,:)
833      vn_tl(:,:,:) = zvn_tlin(:,:,:)
834      ua_tl(:,:,:) = zua_tlin(:,:,:)
835      va_tl(:,:,:) = zva_tlin(:,:,:)
836
837      ! initialize rotn_tl with noise
838      CALL div_cur_tan ( nit000 )
839!zrotn_tlin(:,:,:) = rotn_tl(:,:,:)
840      DO jk = 1, jpk
841        DO jj = nldj, nlej
842           DO ji = nldi, nlei
843              zrotn_tlin(ji,jj,jk) = rotn_tl(ji,jj,jk) 
844            END DO
845         END DO
846      END DO
847      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
848
849      CALL dyn_vor_tan( nit000 )
850      zua_tlout(:,:,:) = ua_tl(:,:,:)
851      zva_tlout(:,:,:) = va_tl(:,:,:)
852
853      !--------------------------------------------------------------------
854      ! Initialize the adjoint variables: dy^* = W dy
855      !--------------------------------------------------------------------
856
857      DO jk = 1, jpk
858        DO jj = nldj, nlej
859           DO ji = nldi, nlei
860              zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
861                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
862                 &               * umask(ji,jj,jk)
863              zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
864                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
865                 &               * vmask(ji,jj,jk)
866            END DO
867         END DO
868      END DO
869      !--------------------------------------------------------------------
870      ! Compute the scalar product: ( L dx )^T W dy
871      !--------------------------------------------------------------------
872
873      zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin )
874      zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin )
875      zsp1   = zsp1_1 + zsp1_2
876
877      !--------------------------------------------------------------------
878      ! Call the adjoint routine: dx^* = L^T dy^*
879      !--------------------------------------------------------------------
880
881      ua_ad(:,:,:) = zua_adin(:,:,:)
882      va_ad(:,:,:) = zva_adin(:,:,:)
883
884      CALL dyn_vor_adj ( nitend )
885
886      zun_adout(:,:,:)   = un_ad(:,:,:)
887      zvn_adout(:,:,:)   = vn_ad(:,:,:)
888      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
889      zua_adout(:,:,:)   = ua_ad(:,:,:)
890      zva_adout(:,:,:)   = va_ad(:,:,:)
891     
892      zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
893      zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
894      zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
895      zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout )
896      zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout )
897      zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5
898
899      ! Compare the scalar products
900
901      ! 14 char:'12345678901234'
902      cl_name = 'dyn_vor_adj   '
903      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
904
905      DEALLOCATE( &
906         & zun_tlin,     &
907         & zvn_tlin,     &
908         & zrotn_tlin,   &
909         & zun_adout,    &
910         & zvn_adout,    &
911         & zrotn_adout,  &
912         & zua_adout,    &
913         & zva_adout,    &
914         & zua_tlin,     &
915         & zva_tlin,     &
916         & zua_tlout,    &
917         & zva_tlout,    &
918         & zua_adin,     &
919         & zva_adin,     &
920         & zau,          &
921         & zav,          &
922         & znu,          &
923         & znv           &
924         & )
925   END SUBROUTINE dyn_vor_adj_tst
926   !!=============================================================================
927#endif
928END MODULE dynvor_tam
Note: See TracBrowser for help on using the repository browser.