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/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynvor_tam.F90 @ 3627

Last change on this file since 3627 was 3627, checked in by rblod, 11 years ago

Correct long lines in TAM 4.3 see ticket #1010

  • Property svn:executable set to *
File size: 83.0 KB
Line 
1MODULE dynvor_tam
2#if defined 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   !! History :  OPA  !  1989-12  (P. Andrich)  vor_ens: Original code
11   !!            5.0  !  1991-11  (G. Madec) vor_ene, vor_mix: Original code
12   !!            6.0  !  1996-01  (G. Madec)  s-coord, suppress work arrays
13   !!            8.5  !  2002-08  (G. Madec)  F90: Free form and module
14   !!   NEMO     1.0  !  2004-02  (G. Madec)  vor_een: Original code
15   !!             -   !  2003-08  (G. Madec)  add vor_ctl
16   !!             -   !  2005-11  (G. Madec)  add dyn_vor (new step architecture)
17   !!            2.0  !  2006-11  (G. Madec)  flux form advection: add metric term
18   !!            3.2  !  2009-04  (R. Benshila)  vvl: correction of een scheme
19   !! History of the TAM module:
20   !!            9.0  ! 2008-06  (A. Vidard) Skeleton
21   !!            9.0  ! 2009-01  (A. Vidard) TAM of the 06-11 version
22   !!            9.0  ! 2010-01  (F. Vigilant) Add een TAM option
23   !!  NEMO      3.2  ! 2010-04 (F. Vigilant) 3.2 version
24   !!  NEMO      3.4  ! 2012-07 (P.-A. Bouttier) 3.4 version
25   !!----------------------------------------------------------------------
26
27   !!----------------------------------------------------------------------
28   !!   dyn_vor     : Update the momentum trend with the vorticity trend
29   !!       vor_ens : enstrophy conserving scheme       (ln_dynvor_ens=T)
30   !!       vor_ene : energy conserving scheme          (ln_dynvor_ene=T)
31   !!       vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
32   !!       vor_een : energy and enstrophy conserving   (ln_dynvor_een=T)
33   !!       vor_ctl : set and control of the different vorticity option
34   !!----------------------------------------------------------------------
35   USE par_oce
36   USE oce
37   USE oce_tam
38   USE divcur
39   USE divcur_tam
40   USE dom_oce
41   USE dynadv
42   USE lbclnk
43   USE lbclnk_tam
44   USE in_out_manager
45   USE gridrandom
46   USE dotprodfld
47   USE tstool_tam
48   USE dommsk
49   USE lib_mpp
50   USE wrk_nemo
51   USE timing
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC   dyn_vor_init_tam
57   PUBLIC   dyn_vor_tan     ! routine called by step_tam.F90
58   PUBLIC   dyn_vor_adj     ! routine called by step_tam.F90
59   PUBLIC   dyn_vor_adj_tst ! routine called by the tst.F90
60
61
62   !!* Namelist nam_dynvor: vorticity term
63   LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
64   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
65   LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
66   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.   !: energy and enstrophy conserving scheme
67
68   INTEGER ::   nvor = 0   ! type of vorticity trend used
69   INTEGER ::   ncor = 1   ! coriolis
70   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
71   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
72
73   !! * Substitutions
74#  include "domzgr_substitute.h90"
75#  include "vectopt_loop_substitute.h90"
76
77CONTAINS
78
79   SUBROUTINE dyn_vor_tan( kt )
80      !!----------------------------------------------------------------------
81      !!
82      !! ** Purpose of the direct routine:
83      !!               compute the lateral ocean tracer physics.
84      !!
85      !! ** Action : - Update (ua,va) with the now vorticity term trend
86      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
87      !!               and planetary vorticity trends) ('key_trddyn')
88      !!----------------------------------------------------------------------
89      !!
90      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
91      !!----------------------------------------------------------------------
92      !
93      IF( nn_timing == 1 )  CALL timing_start('dyn_vor_tan')
94      !
95      !!                                          ! vorticity term
96      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
97      !
98      CASE ( -1 )                                      ! esopa: test all possibility with control print
99         CALL vor_ene_tan( kt, ntot, ua_tl, va_tl )
100         CALL vor_ens_tan( kt, ntot, ua_tl, va_tl )
101         CALL vor_mix_tan( kt )
102         CALL vor_een_tan( kt, ntot, ua_tl, va_tl )
103         !
104      CASE ( 0 )                                       ! energy conserving scheme
105         CALL vor_ene_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
106         !
107      CASE ( 1 )                                       ! enstrophy conserving scheme
108         CALL vor_ens_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
109         !
110      CASE ( 2 )                                       ! mixed ene-ens scheme
111         CALL vor_mix_tan( kt )                               ! total vorticity (mix=ens-ene)
112         !
113      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
114         CALL vor_een_tan( kt, ntot, ua_tl, va_tl )                ! total vorticity
115         !
116      END SELECT
117      !
118      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor_tan')
119      !
120   END SUBROUTINE dyn_vor_tan
121   SUBROUTINE vor_ene_tan( kt, kvor, pua_tl, pva_tl )
122      !!----------------------------------------------------------------------
123      !!                  ***  ROUTINE vor_ene  ***
124      !!
125      !! ** Purpose :   Compute the now total vorticity trend and add it to
126      !!      the general trend of the momentum equation.
127      !!
128      !! ** Method  :   Trend evaluated using now fields (centered in time)
129      !!      and the Sadourny (1975) flux form formulation : conserves the
130      !!      horizontal kinetic energy.
131      !!      The trend of the vorticity term is given by:
132      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
133      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
134      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
135      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
136      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
137      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
138      !!      Add this trend to the general momentum trend (ua,va):
139      !!          (ua,va) = (ua,va) + ( voru , vorv )
140      !!
141      !! ** Action : - Update (ua,va) with the now vorticity term trend
142      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
143      !!               and planetary vorticity trends) ('key_trddyn')
144      !!
145      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
146      !!----------------------------------------------------------------------
147      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
148      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
149      !                                                           ! =nrvm (relative vorticity or metric)
150      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_tl ! total u-trend
151      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_tl ! total v-trend
152      !!
153      INTEGER  ::   ji, jj, jk         ! dummy loop indices
154      REAL(wp) ::   zx1, zy1, zfact2   ! temporary scalars
155      REAL(wp) ::   zx1tl, zy1tl   ! temporary scalars
156      REAL(wp) ::   zx2, zy2           !    "         "
157      REAL(wp) ::   zx2tl, zy2tl           !    "         "
158      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! temporary 2D workspace
159      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxtl, zwytl, zwztl   ! temporary 2D workspace
160      !!----------------------------------------------------------------------
161      !
162      IF( nn_timing == 1 )  CALL timing_start('vor_ene_tan')
163      !
164      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
165      CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl )
166      !
167      IF( kt == nit000 ) THEN
168         IF(lwp) WRITE(numout,*)
169         IF(lwp) WRITE(numout,*) 'dyn:vor_ene_tan : vorticity term: energy conserving scheme'
170         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
171      ENDIF
172
173      zfact2 = 0.5 * 0.5      ! Local constant initialization
174
175!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
176      !                                                ! ===============
177      DO jk = 1, jpkm1                                 ! Horizontal slab
178         !                                             ! ===============
179         !
180         ! Potential vorticity and horizontal fluxes
181         ! -----------------------------------------
182         SELECT CASE( kvor )      ! vorticity considered
183         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)  ; zwztl(:,:) = 0.
184         CASE ( 2 )   ;    zwz(:,:) =   rotn(:,:,jk)  ; zwztl(:,:) =   rotn_tl(:,:,jk)                ! relative  vorticity
185         CASE ( 3 )                                                ! metric term
186            DO jj = 1, jpjm1
187               DO ji = 1, fs_jpim1   ! vector opt.
188                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
189                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
190                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
191                  zwztl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
192                       &         - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
193                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
194               END DO
195            END DO
196         CASE ( 4 )    ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ;   zwztl(:,:) = rotn_tl(:,:,jk)    ! total (relative + planetary vorticity)
197         CASE ( 5 )                                                ! total (coriolis + metric)
198            DO jj = 1, jpjm1
199               DO ji = 1, fs_jpim1   ! vector opt.
200                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
201                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
202                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
203                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
204                       &       )
205                  zwztl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
206                       &           - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
207                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
208               END DO
209            END DO
210         END SELECT
211
212         IF( ln_sco ) THEN
213            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
214            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
215            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
216            zwztl(:,:) = zwztl(:,:) / fse3f(:,:,jk)
217            zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
218            zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
219         ELSE
220            zwx(:,:) = e2u(:,:) * un(:,:,jk)
221            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
222            zwxtl(:,:) = e2u(:,:) * un_tl(:,:,jk)
223            zwytl(:,:) = e1v(:,:) * vn_tl(:,:,jk)
224         ENDIF
225
226         ! Compute and add the vorticity term trend
227         ! ----------------------------------------
228         DO jj = 2, jpjm1
229            DO ji = fs_2, fs_jpim1   ! vector opt.
230               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
231               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
232               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
233               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
234               zy1tl = zwytl(ji,jj-1) + zwytl(ji+1,jj-1)
235               zy2tl = zwytl(ji,jj  ) + zwytl(ji+1,jj  )
236               zx1tl = zwxtl(ji-1,jj) + zwxtl(ji-1,jj+1)
237               zx2tl = zwxtl(ji  ,jj) + zwxtl(ji  ,jj+1)
238               pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwztl(ji  ,jj-1) * zy1 + zwztl(ji,jj) * zy2 )   &
239                                &                  + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1tl + zwz(ji,jj) * zy2tl )
240               pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwztl(ji-1,jj  ) * zx1 + zwztl(ji,jj) * zx2 )   &
241                                &                  - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1tl + zwz(ji,jj) * zx2tl )
242            END DO
243         END DO
244         !                                             ! ===============
245      END DO                                           !   End of slab
246      !                                                ! ===============
247      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
248      CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl )
249      !
250      IF( nn_timing == 1 )  CALL timing_stop('vor_ene_tan')
251      !
252   END SUBROUTINE vor_ene_tan
253
254
255   SUBROUTINE vor_mix_tan( kt )
256      !!----------------------------------------------------------------------
257      !!                 ***  ROUTINE vor_mix  ***
258      !!
259      !! ** Purpose :   Compute the now total vorticity trend and add it to
260      !!      the general trend of the momentum equation.
261      !!
262      !! ** Method  :   Trend evaluated using now fields (centered in time)
263      !!      Mixte formulation : conserves the potential enstrophy of a hori-
264      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
265      !!      ticity term and the horizontal kinetic energy for (f x uh), the
266      !!      coriolis term. the now trend of the vorticity term is given by:
267      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
268      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
269      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
270      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
271      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
272      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
273      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
274      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
275      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
276      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
277      !!      Add this now trend to the general momentum trend (ua,va):
278      !!          (ua,va) = (ua,va) + ( voru , vorv )
279      !!
280      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
281      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
282      !!               and planetary vorticity trends) ('key_trddyn')
283      !!
284      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
285      !!----------------------------------------------------------------------
286      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
287      !!
288      INTEGER  ::   ji, jj, jk   ! dummy loop indices
289      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! temporary scalars
290      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !    "         "
291      REAL(wp) ::   zuatl, zcuatl, zx1tl, zy1tl   ! temporary scalars
292      REAL(wp) ::   zvatl, zcvatl, zx2tl, zy2tl   !    "         "
293      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxtl, zwytl, zwztl, zwwtl   ! temporary 3D workspace
294      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! temporary 3D workspace
295      !!----------------------------------------------------------------------
296      !
297      IF( nn_timing == 1 )  CALL timing_start('vor_mix_tan')
298      !
299      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )
300      CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl, zwwtl )
301      !
302      IF( kt == nit000 ) THEN
303         IF(lwp) WRITE(numout,*)
304         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
305         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
306      ENDIF
307
308      zfact1 = 0.5 * 0.25      ! Local constant initialization
309      zfact2 = 0.5 * 0.5
310
311!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
312      !                                                ! ===============
313      DO jk = 1, jpkm1                                 ! Horizontal slab
314         !                                             ! ===============
315         !
316         ! Relative and planetary potential vorticity and horizontal fluxes
317         ! ----------------------------------------------------------------
318         IF( ln_sco ) THEN
319            IF( ln_dynadv_vec ) THEN
320               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
321               zwwtl(:,:) = rotn_tl(:,:,jk) / fse3f(:,:,jk)
322            ELSE
323               DO jj = 1, jpjm1
324                  DO ji = 1, fs_jpim1   ! vector opt.
325                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
326                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
327                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
328                     zwwtl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
329                        &           - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
330                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
331                  END DO
332               END DO
333            ENDIF
334            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
335            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
336            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
337            zwztl(:,:) = 0._wp
338            zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
339            zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
340         ELSE
341            IF( ln_dynadv_vec ) THEN
342               zww(:,:) = rotn(:,:,jk)
343               zwwtl(:,:) = rotn_tl(:,:,jk)
344            ELSE
345               DO jj = 1, jpjm1
346                  DO ji = 1, fs_jpim1   ! vector opt.
347                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
348                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
349                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
350                     zwwtl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
351                        &           - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
352                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
353                  END DO
354               END DO
355            ENDIF
356            zwz(:,:) = ff (:,:)
357            zwx(:,:) = e2u(:,:) * un(:,:,jk)
358            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
359            zwztl(:,:) = 0.0_wp
360            zwxtl(:,:) = e2u(:,:) * un_tl(:,:,jk)
361            zwytl(:,:) = e1v(:,:) * vn_tl(:,:,jk)
362         ENDIF
363
364         ! Compute and add the vorticity term trend
365         ! ----------------------------------------
366         DO jj = 2, jpjm1
367            DO ji = fs_2, fs_jpim1   ! vector opt.
368               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
369               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
370               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
371               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
372               zy1tl = ( zwytl(ji,jj-1) + zwytl(ji+1,jj-1) ) / e1u(ji,jj)
373               zy2tl = ( zwytl(ji,jj  ) + zwytl(ji+1,jj  ) ) / e1u(ji,jj)
374               zx1tl = ( zwxtl(ji-1,jj) + zwxtl(ji-1,jj+1) ) / e2v(ji,jj)
375               zx2tl = ( zwxtl(ji  ,jj) + zwxtl(ji  ,jj+1) ) / e2v(ji,jj)
376               ! enstrophy conserving formulation for relative vorticity term
377               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
378               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
379               zuatl = zfact1 * ( zwwtl(ji  ,jj-1) + zwwtl(ji,jj) ) * ( zy1 + zy2 )   &
380                 &   + zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1tl + zy2tl )
381               zvatl =-zfact1 * ( zwwtl(ji-1,jj  ) + zwwtl(ji,jj) ) * ( zx1 + zx2 )   &
382                 &   - zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1tl + zx2tl )
383               ! energy conserving formulation for planetary vorticity term
384               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
385               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
386               zcuatl = zfact2 * ( zwztl(ji  ,jj-1) * zy1 + zwztl(ji,jj) * zy2 )   &
387                  &   + zfact2 * ( zwz(ji  ,jj-1) * zy1tl + zwz(ji,jj) * zy2tl )
388               zcvatl =-zfact2 * ( zwztl(ji-1,jj  ) * zx1 + zwztl(ji,jj) * zx2 )   &
389                  &    -zfact2 * ( zwz(ji-1,jj  ) * zx1tl + zwz(ji,jj) * zx2tl )
390               ! mixed vorticity trend added to the momentum trends
391               ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zcuatl + zuatl
392               va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zcvatl + zvatl
393            END DO
394         END DO
395         !                                             ! ===============
396      END DO                                           !   End of slab
397      !                                                ! ===============
398      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )
399      CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl, zwwtl )
400      !
401      IF( nn_timing == 1 )  CALL timing_stop('vor_mix_tan')
402      !
403   END SUBROUTINE vor_mix_tan
404   SUBROUTINE vor_ens_tan( kt, kvor, pua_tl, pva_tl )
405      !!----------------------------------------------------------------------
406      !!                ***  ROUTINE vor_ens_tan  ***
407      !!
408      !! ** Purpose of the direct routine:
409      !!      Compute the now total vorticity trend and add it to
410      !!      the general trend of the momentum equation.
411      !!
412      !! ** Method of the direct routine:
413      !!      Trend evaluated using now fields (centered in time)
414      !!      and the Sadourny (1975) flux FORM formulation : conserves the
415      !!      potential enstrophy of a horizontally non-divergent flow. the
416      !!      trend of the vorticity term is given by:
417      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
418      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
419      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
420      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
421      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
422      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
423      !!      Add this trend to the general momentum trend (ua,va):
424      !!          (ua,va) = (ua,va) + ( voru , vorv )
425      !!
426      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
427      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
428      !!               and planetary vorticity trends) ('key_trddyn')
429      !!
430      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
431      !!----------------------------------------------------------------------
432      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
433      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
434         !                                                        ! =nrvm (relative vorticity or metric)
435      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_tl    ! total u-trend
436      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_tl    ! total v-trend
437      !!
438      INTEGER  ::   ji, jj, jk           ! dummy loop indices
439      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
440      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz         ! temporary 3D workspace
441      REAL(wp) ::   zuavtl, zvautl   ! temporary scalars
442      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxtl, zwytl, zwztl   ! temporary 3D workspace
443      !!----------------------------------------------------------------------
444      !
445      IF( nn_timing == 1 )  CALL timing_start('vor_ens_tan')
446      !
447      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
448      CALL wrk_alloc( jpi, jpj, zwxtl, zwytl, zwztl )
449      !
450      IF( kt == nit000 ) THEN
451         IF(lwp) WRITE(numout,*)
452         IF(lwp) WRITE(numout,*) 'dyn_vor_ens_tan : vorticity term: enstrophy conserving scheme'
453         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
454      ENDIF
455      ! Local constant initialization
456      zfact1 = 0.5 * 0.25
457
458!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz  zwxtl, zwytl, zwztl)
459      !                                                ! ===============
460      DO jk = 1, jpkm1                                 ! Horizontal slab
461         !                                             ! ===============
462         ! Potential vorticity and horizontal fluxes
463         ! -----------------------------------------
464         SELECT CASE( kvor )      ! vorticity considered
465         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
466         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
467         CASE ( 3 )                                                ! metric term
468            DO jj = 1, jpjm1
469               DO ji = 1, fs_jpim1   ! vector opt.
470                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
471                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
472                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
473               END DO
474            END DO
475         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
476         CASE ( 5 )                                                ! total (coriolis + metric)
477            DO jj = 1, jpjm1
478               DO ji = 1, fs_jpim1   ! vector opt.
479                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
480                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
481                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
482                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
483                       &       )
484               END DO
485            END DO
486         END SELECT
487
488         IF( ln_sco ) THEN
489            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
490               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
491                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
492                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
493                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
494               END DO
495            END DO
496         ELSE
497            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
498               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
499                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
500                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
501               END DO
502            END DO
503         ENDIF
504
505         ! ===================
506         ! Tangent counterpart
507         ! ===================
508         ! Potential vorticity and horizontal fluxes
509         ! -----------------------------------------
510         SELECT CASE( kvor )      ! vorticity considered
511         CASE ( 1 )     ;   zwztl(:,:) =  0.0_wp                      ! planetary vorticity (Coriolis)
512         CASE ( 2 ,4)   ;   zwztl(:,:) =  rotn_tl(:,:,jk)             ! relative  vorticity
513         CASE ( 3 ,5 )                                                ! metric term
514            DO jj = 1, jpjm1
515               DO ji = 1, fs_jpim1   ! vector opt.
516                  zwztl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) )  &
517                     &               * ( e2v(ji+1,jj  ) - e2v(ji,jj) )           &
518                     &             - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) )  &
519                     &               * ( e1u(ji  ,jj+1) - e1u(ji,jj) )           &
520                     &           ) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
521               END DO
522            END DO
523         END SELECT
524
525         IF( ln_sco ) THEN
526            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
527               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
528                  zwztl(ji,jj) = zwztl(ji,jj) / fse3f(ji,jj,jk)
529                  zwxtl(ji,jj) = e2u(ji,jj)   * fse3u(ji,jj,jk) * un_tl(ji,jj,jk)
530                  zwytl(ji,jj) = e1v(ji,jj)   * fse3v(ji,jj,jk) * vn_tl(ji,jj,jk)
531               END DO
532            END DO
533         ELSE
534            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
535               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
536                  zwxtl(ji,jj) = e2u(ji,jj) * un_tl(ji,jj,jk)
537                  zwytl(ji,jj) = e1v(ji,jj) * vn_tl(ji,jj,jk)
538               END DO
539            END DO
540         ENDIF
541
542         ! Compute and add the vorticity term trend
543         ! ----------------------------------------
544         DO jj = 2, jpjm1
545            DO ji = fs_2, fs_jpim1   ! vector opt.
546               zuav   = zfact1 / e1u(ji,jj) * ( zwy(  ji  ,jj-1) + zwy(  ji+1,jj-1)   &
547                  &                           + zwy(  ji  ,jj  ) + zwy(  ji+1,jj  ) )
548               zvau   =-zfact1 / e2v(ji,jj) * ( zwx(  ji-1,jj  ) + zwx(  ji-1,jj+1)   &
549                  &                           + zwx(  ji  ,jj  ) + zwx(  ji  ,jj+1) )
550               zuavtl = zfact1 / e1u(ji,jj) * ( zwytl(ji  ,jj-1) + zwytl(ji+1,jj-1)   &
551                  &                           + zwytl(ji  ,jj  ) + zwytl(ji+1,jj  ) )
552               zvautl =-zfact1 / e2v(ji,jj) * ( zwxtl(ji-1,jj  ) + zwxtl(ji-1,jj+1)   &
553                  &                           + zwxtl(ji  ,jj  ) + zwxtl(ji  ,jj+1) )
554               pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk)                   &
555                  &     + zuavtl * ( zwz(  ji,jj-1) + zwz(  ji,jj) ) &
556                  &     + zuav   * ( zwztl(ji,jj-1) + zwztl(ji,jj) )
557               pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk)                   &
558                  &     + zvautl * ( zwz(  ji-1,jj) + zwz(  ji,jj) ) &
559                  &     + zvau   * ( zwztl(ji-1,jj) + zwztl(ji,jj) )
560            END DO
561         END DO
562         !                                             ! ===============
563      END DO                                           !   End of slab
564      !                                                ! ===============
565      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
566      CALL wrk_dealloc( jpi, jpj, zwxtl, zwytl, zwztl )
567      !
568      IF( nn_timing == 1 )  CALL timing_stop('vor_ens_tan')
569      !
570   END SUBROUTINE vor_ens_tan
571
572   SUBROUTINE vor_een_tan( kt, kvor, pua_tl, pva_tl )
573      !!----------------------------------------------------------------------
574      !!                ***  ROUTINE vor_een_tan  ***
575      !!
576      !! ** Purpose :   Compute the now total vorticity trend and add it to
577      !!      the general trend of the momentum equation.
578      !!
579      !! ** Method  :   Trend evaluated using now fields (centered in time)
580      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
581      !!      both the horizontal kinetic energy and the potential enstrophy
582      !!      when horizontal divergence is zero (see the NEMO documentation)
583      !!      Add this trend to the general momentum trend (ua,va).
584      !!
585      !! ** Action : - Update (ua,va) with the now vorticity term trend
586      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
587      !!               and planetary vorticity trends) ('key_trddyn')
588      !!
589      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
590      !!----------------------------------------------------------------------
591      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
592      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
593         !                                                        ! =nrvm (relative vorticity or metric)
594      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_tl ! total u-trend
595      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_tl ! total v-trend
596      !!
597      INTEGER ::   ji, jj, jk          ! dummy loop indices
598      INTEGER :: ierr
599      REAL(wp) ::   zfac12, zua, zva   ! temporary scalars
600      REAL(wp) ::   zuatl, zvatl       ! temporary scalars
601      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz                    ! temporary 2D workspace
602      REAL(wp), POINTER, DIMENSION(:,:) ::   ztnw, ztne, ztsw, ztse           ! temporary 3D workspace
603      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxtl, zwytl, zwztl              ! temporary 2D workspace
604      REAL(wp), POINTER, DIMENSION(:,:) ::   ztnwtl, ztnetl, ztswtl, ztsetl   ! temporary 3D workspace
605#if defined key_vvl
606      REAL(wp), POINTER, DIMENSION(:,:,:)       ::   ze3f_tl
607#else
608      REAL(wp), POINTER, DIMENSION(:,:,:), SAVE ::   ze3f_tl
609#endif
610      !!----------------------------------------------------------------------
611      IF( nn_timing == 1 )  CALL timing_start('vor_een_tan')
612      !
613      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        )
614      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )
615      CALL wrk_alloc( jpi, jpj,      zwxtl , zwytl , zwztl        )
616      CALL wrk_alloc( jpi, jpj,      ztnwtl, ztnetl, ztswtl, ztsetl )
617      zuatl = 0.0_wp ; zvatl = 0.0_wp
618      zwx (:,:) = 0.0_wp ; zwy (:,:) = 0.0_wp ; zwz (:,:) = 0.0_wp
619      ztnw(:,:) = 0.0_wp ; ztne(:,:) = 0.0_wp ; ztsw(:,:) = 0.0_wp ; ztse(:,:) = 0.0_wp
620      zwxtl (:,:) = 0.0_wp ; zwytl (:,:) = 0.0_wp ; zwztl (:,:) = 0.0_wp
621      ztnwtl(:,:) = 0.0_wp ; ztnetl(:,:) = 0.0_wp ; ztswtl(:,:) = 0.0_wp ; ztsetl(:,:) = 0.0_wp
622#if defined key_vvl
623      CALL wrk_alloc( jpi, jpj, jpk, ze3f_tl                   )
624#endif
625      IF( kt == nit000 ) THEN
626         IF(lwp) WRITE(numout,*)
627         IF(lwp) WRITE(numout,*) 'dyn:vor_een_tam : vorticity term: energy and enstrophy conserving scheme'
628         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
629         IF( .NOT.lk_vvl ) THEN
630            ALLOCATE( ze3f_tl(jpi,jpj,jpk) , STAT=ierr )
631            IF( lk_mpp    )   CALL mpp_sum ( ierr )
632            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
633            ze3f_tl = 0._wp
634         ENDIF
635      ENDIF
636
637      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
638         DO jk = 1, jpk
639            DO jj = 1, jpjm1
640               DO ji = 1, jpim1
641                  ze3f_tl(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
642                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25_wp
643                  IF( ze3f_tl(ji,jj,jk) /= 0.0_wp )   ze3f_tl(ji,jj,jk) = 1.0_wp / ze3f_tl(ji,jj,jk)
644               END DO
645            END DO
646         END DO
647         CALL lbc_lnk( ze3f_tl, 'F', 1._wp )
648      ENDIF
649
650      zfac12 = 1.0_wp / 12.0_wp      ! Local constant initialization
651
652
653!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
654      !                                                ! ===============
655      DO jk = 1, jpkm1                                 ! Horizontal slab
656         !                                             ! ===============
657
658         ! Potential vorticity and horizontal fluxes
659         ! -----------------------------------------
660         SELECT CASE( kvor )      ! vorticity considered
661         CASE ( 1 )
662            zwz(:,:) = ff(:,:)      * ze3f_tl(:,:,jk)   ! planetary vorticity (Coriolis)
663            zwztl(:,:) = 0.0_wp
664         CASE ( 2 )
665            zwz(:,:) = rotn(:,:,jk) * ze3f_tl(:,:,jk)   ! relative  vorticity
666            zwztl(:,:) = rotn_tl(:,:,jk) * ze3f_tl(:,:,jk)
667         CASE ( 3 )                                                ! metric term
668            DO jj = 1, jpjm1
669               DO ji = 1, fs_jpim1   ! vector opt.
670                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
671                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
672                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_tl(ji,jj,jk)
673               END DO
674            END DO
675            CALL lbc_lnk( zwz, 'F', 1._wp )
676            DO jj = 1, jpjm1
677               DO ji = 1, fs_jpim1   ! vector opt.
678                  zwztl(ji,jj) = (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
679                       &         - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
680                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_tl(ji,jj,jk)
681               END DO
682            END DO
683            CALL lbc_lnk( zwztl, 'F', 1._wp )
684         CASE ( 4 )
685            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f_tl(:,:,jk) ! total (relative + planetary vorticity)
686            zwztl(:,:) = ( rotn_tl(:,:,jk) ) * ze3f_tl(:,:,jk)
687         CASE ( 5 )                                                ! total (coriolis + metric)
688            DO jj = 1, jpjm1
689               DO ji = 1, fs_jpim1   ! vector opt.
690                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
691                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
692                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
693                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
694                       &       ) * ze3f_tl(ji,jj,jk)
695               END DO
696            END DO
697            CALL lbc_lnk( zwz, 'F', 1._wp )
698            DO jj = 1, jpjm1
699               DO ji = 1, fs_jpim1   ! vector opt.
700                  zwztl(ji,jj) = ( (   ( vn_tl(ji+1,jj  ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
701                       &           - ( un_tl(ji  ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
702                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
703                       &       ) * ze3f_tl(ji,jj,jk)
704               END DO
705            END DO
706            CALL lbc_lnk( zwztl, 'F', 1._wp )
707         END SELECT
708
709         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
710         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
711
712         zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)
713         zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
714
715         ! Compute and add the vorticity term trend
716         ! ----------------------------------------
717         jj=2
718         ztne(1,:)   = 0.0_wp ; ztnw(1,:)   = 0.0_wp ; ztse(1,:)   = 0.0_wp ; ztsw(1,:)   = 0.0_wp
719         ztnetl(1,:) = 0.0_wp ; ztnwtl(1,:) = 0.0_wp ; ztsetl(1,:) = 0.0_wp ; ztswtl(1,:) = 0.0_wp
720         DO ji = 2, jpi
721               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
722               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
723               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
724               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
725
726               ztnetl(ji,jj) = zwztl(ji-1,jj  ) + zwztl(ji  ,jj  ) + zwztl(ji  ,jj-1)
727               ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj  ) + zwztl(ji  ,jj  )
728               ztsetl(ji,jj) = zwztl(ji  ,jj  ) + zwztl(ji  ,jj-1) + zwztl(ji-1,jj-1)
729               ztswtl(ji,jj) = zwztl(ji  ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj  )
730         END DO
731         DO jj = 3, jpj
732            DO ji = fs_2, jpi   ! vector opt.
733               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
734               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
735               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
736               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
737
738               ztnetl(ji,jj) = zwztl(ji-1,jj  ) + zwztl(ji  ,jj  ) + zwztl(ji  ,jj-1)
739               ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj  ) + zwztl(ji  ,jj  )
740               ztsetl(ji,jj) = zwztl(ji  ,jj  ) + zwztl(ji  ,jj-1) + zwztl(ji-1,jj-1)
741               ztswtl(ji,jj) = zwztl(ji  ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj  )
742            END DO
743         END DO
744         DO jj = 2, jpjm1
745            DO ji = fs_2, fs_jpim1   ! vector opt.
746               zuatl = + zfac12 / e1u(ji,jj) * (  ztnetl(ji,jj  ) * zwy(ji  ,jj  ) + ztne(ji,jj  ) * zwytl(ji  ,jj  )  &
747                  &                             + ztnwtl(ji+1,jj) * zwy(ji+1,jj  ) + ztnw(ji+1,jj) * zwytl(ji+1,jj  )  &
748                  &                             + ztsetl(ji,jj  ) * zwy(ji  ,jj-1) + ztse(ji,jj  ) * zwytl(ji  ,jj-1)    &
749                  &                             + ztswtl(ji+1,jj) * zwy(ji+1,jj-1) + ztsw(ji+1,jj) * zwytl(ji+1,jj-1))
750
751               zvatl = - zfac12 / e2v(ji,jj) * (  ztswtl(ji,jj+1) * zwx(ji-1,jj+1) + ztsw(ji,jj+1) * zwxtl(ji-1,jj+1)  &
752                  &                           +   ztsetl(ji,jj+1) * zwx(ji  ,jj+1) + ztse(ji,jj+1) * zwxtl(ji  ,jj+1)   &
753                  &                           +   ztnwtl(ji,jj  ) * zwx(ji-1,jj  ) + ztnw(ji,jj  ) * zwxtl(ji-1,jj  )   &
754                  &                           +   ztnetl(ji,jj  ) * zwx(ji  ,jj  ) + ztne(ji,jj  ) * zwxtl(ji  ,jj  ) )
755               pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zuatl
756               pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) + zvatl
757            END DO
758         END DO
759         !                                             ! ===============
760      END DO                                           !   End of slab
761      !                                                ! ===============
762      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        )
763      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )
764      CALL wrk_dealloc( jpi, jpj,      zwxtl , zwytl , zwztl        )
765      CALL wrk_dealloc( jpi, jpj,      ztnwtl, ztnetl, ztswtl, ztsetl )
766#if defined key_vvl
767      CALL wrk_dealloc( jpi, jpj, jpk, ze3f_tl                   )
768#endif
769      !
770      IF( nn_timing == 1 )  CALL timing_stop('vor_een_tan')
771      !
772   END SUBROUTINE vor_een_tan
773
774
775   SUBROUTINE dyn_vor_adj( kt )
776      !!----------------------------------------------------------------------
777      !!
778      !! ** Purpose of the direct routine:
779      !!               compute the lateral ocean tracer physics.
780      !!
781      !! ** Action : - Update (ua,va) with the now vorticity term trend
782      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
783      !!               and planetary vorticity trends) ('key_trddyn')
784      !!----------------------------------------------------------------------
785      !!
786      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
787      !!----------------------------------------------------------------------
788      !
789      IF( nn_timing == 1 )  CALL timing_start('dyn_vor_adj')
790      !
791      !                                          ! vorticity term
792      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
793      !
794      CASE ( -1 )                                      ! esopa: test all possibility with control print
795         CALL vor_een_adj( kt, ntot, ua_ad, va_ad )
796!         CALL vor_mix_adj( kt )
797         CALL vor_ens_adj( kt, ntot, ua_ad, va_ad )
798!         CALL vor_ene_adj( kt, ntot, ua_ad, va_ad )
799         !
800      CASE ( 0 )                                       ! energy conserving scheme
801         CALL ctl_stop ('vor_ene_adj not available yet')
802!         CALL vor_ene_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
803         !
804      CASE ( 1 )                                       ! enstrophy conserving scheme
805         CALL vor_ens_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
806         !
807      CASE ( 2 )                                       ! mixed ene-ens scheme
808         CALL ctl_stop ('vor_mix_adj not available yet')
809!         CALL vor_mix_adj( kt )                               ! total vorticity (mix=ens-ene)
810         !
811      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
812         CALL vor_een_adj( kt, ntot, ua_ad, va_ad )                ! total vorticity
813         !
814      END SELECT
815      !
816      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor_adj')
817      !
818   END SUBROUTINE dyn_vor_adj
819   SUBROUTINE vor_ens_adj( kt, kvor, pua_ad, pva_ad )
820      !!----------------------------------------------------------------------
821      !!                ***  ROUTINE vor_ens_adj  ***
822      !!
823      !! ** Purpose of the direct routine:
824      !!      Compute the now total vorticity trend and add it to
825      !!      the general trend of the momentum equation.
826      !!
827      !! ** Method of the direct routine:
828      !!      Trend evaluated using now fields (centered in time)
829      !!      and the Sadourny (1975) flux FORM formulation : conserves the
830      !!      potential enstrophy of a horizontally non-divergent flow. the
831      !!      trend of the vorticity term is given by:
832      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
833      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
834      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
835      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
836      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
837      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
838      !!      Add this trend to the general momentum trend (ua,va):
839      !!          (ua,va) = (ua,va) + ( voru , vorv )
840      !!
841      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
842      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
843      !!               and planetary vorticity trends) ('key_trddyn')
844      !!
845      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
846      !!----------------------------------------------------------------------
847      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
848      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
849         !                                                        ! =nrvm (relative vorticity or metric)
850      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_ad    ! total u-trend
851      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_ad    ! total v-trend
852      !!
853      INTEGER  ::   ji, jj, jk           ! dummy loop indices
854      REAL(wp) ::   zfact1               ! temporary scalars
855      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz         ! temporary 3D workspace
856      REAL(wp) ::   zuav, zvau       ! temporary scalars
857      REAL(wp) ::   zuavad, zvauad   ! temporary scalars
858      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxad, zwyad, zwzad   ! temporary 3D workspace
859      !!----------------------------------------------------------------------
860      !
861      IF( nn_timing == 1 )  CALL timing_start('vor_ens_adj')
862      !
863      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
864      CALL wrk_alloc( jpi, jpj, zwxad, zwyad, zwzad )
865      !
866      IF( kt == nitend ) THEN
867         IF(lwp) WRITE(numout,*)
868         IF(lwp) WRITE(numout,*) 'dyn_vor_ens_adj : vorticity term: enstrophy conserving scheme'
869         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
870      ENDIF
871
872      ! Local constant initialization
873      zfact1 = 0.5 * 0.25
874
875!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zwxad, zwyad, zwzad )
876      !                                                ! ===============
877      DO jk = 1, jpkm1                                 ! Horizontal slab
878         !                                             ! ===============
879         ! Potential vorticity and horizontal fluxes
880         ! -----------------------------------------
881         SELECT CASE( kvor )      ! vorticity considered
882         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
883         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
884         CASE ( 3 )                                                ! metric term
885            DO jj = 1, jpjm1
886               DO ji = 1, fs_jpim1   ! vector opt.
887                  zwz(ji,jj) = ( ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )  &
888                       &       - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) ) )&
889                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
890               END DO
891            END DO
892         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
893         CASE ( 5 )                                                ! total (coriolis + metric)
894            DO jj = 1, jpjm1
895               DO ji = 1, fs_jpim1   ! vector opt.
896                  zwz(ji,jj) = ( ff (ji,jj)                         &
897                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
898                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
899                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )  &
900                       &       )
901               END DO
902            END DO
903         END SELECT
904
905         IF( ln_sco ) THEN
906            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
907               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
908                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
909                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
910                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
911               END DO
912            END DO
913         ELSE
914            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
915               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
916                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
917                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
918               END DO
919            END DO
920         ENDIF
921         ! ===================
922         ! Adjoint counterpart
923         ! ===================
924         zuavad     = 0.0_wp
925         zvauad     = 0.0_wp
926         zwxad(:,:) = 0.0_wp
927         zwyad(:,:) = 0.0_wp
928         zwzad(:,:) = 0.0_wp
929         ! Compute and add the vorticity term trend
930         ! ----------------------------------------
931         DO jj = jpjm1, 2, -1
932            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
933               ! Compute and add the vorticity term trend
934               ! ----------------------------------------
935               !- Direct counterpart
936               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
937                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
938               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
939                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
940               !-
941               zvauad = zvauad + pva_ad(ji,jj,jk) * ( zwz(ji-1,jj) + zwz(ji,jj) )
942               zwzad(ji-1,jj) = zwzad(ji-1,jj) + pva_ad(ji,jj,jk) * zvau
943               zwzad(ji  ,jj) = zwzad(ji  ,jj) + pva_ad(ji,jj,jk) * zvau
944
945               zuavad = zuavad + pua_ad(ji,jj,jk) * ( zwz(ji,jj-1) + zwz(ji,jj) )
946               zwzad(ji,jj-1) = zwzad(ji,jj-1) + pua_ad(ji,jj,jk) * zuav
947               zwzad(ji,jj  ) = zwzad(ji,jj  ) + pua_ad(ji,jj,jk) * zuav
948
949               zwxad(ji-1,jj  ) = zwxad(ji-1,jj  ) - zvauad * zfact1 / e2v(ji,jj)
950               zwxad(ji-1,jj+1) = zwxad(ji-1,jj+1) - zvauad * zfact1 / e2v(ji,jj)
951               zwxad(ji  ,jj  ) = zwxad(ji  ,jj  ) - zvauad * zfact1 / e2v(ji,jj)
952               zwxad(ji  ,jj+1) = zwxad(ji  ,jj+1) - zvauad * zfact1 / e2v(ji,jj)
953               zvauad = 0.0_wp
954
955               zwyad(ji  ,jj-1) = zwyad(ji  ,jj-1) + zuavad * zfact1 / e1u(ji,jj)
956               zwyad(ji+1,jj-1) = zwyad(ji+1,jj-1) + zuavad * zfact1 / e1u(ji,jj)
957               zwyad(ji  ,jj  ) = zwyad(ji  ,jj  ) + zuavad * zfact1 / e1u(ji,jj)
958               zwyad(ji+1,jj  ) = zwyad(ji+1,jj  ) + zuavad * zfact1 / e1u(ji,jj)
959               zuavad = 0.0_wp
960            END DO
961         END DO
962         IF( ln_sco ) THEN
963            DO jj = jpj, 1, -1     ! caution: don't use (:,:) for this loop
964               DO ji = jpi, 1, -1  ! it causes optimization problems on NEC in auto-tasking
965                  zwzad(ji,jj) = zwzad(ji,jj) / fse3f(ji,jj,jk)
966                  un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + zwxad(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
967                  vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + zwyad(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,jk)
968                  zwxad(ji,jj) = 0.0_wp
969                  zwyad(ji,jj) = 0.0_wp
970               END DO
971            END DO
972         ELSE
973            DO jj = jpj, 1, -1    ! caution: don't use (:,:) for this loop
974               DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking
975                  un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e2u(ji,jj) * zwxad(ji,jj)
976                  vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e1v(ji,jj) * zwyad(ji,jj)
977                  zwxad(ji,jj) = 0.0_wp
978                  zwyad(ji,jj) = 0.0_wp
979               END DO
980            END DO
981         ENDIF
982         ! Potential vorticity and horizontal fluxes
983         ! -----------------------------------------
984         SELECT CASE( kvor )      ! vorticity considered
985         CASE ( 1 )                     ! planetary vorticity (Coriolis)
986            zwzad(:,:) =  0.0_wp
987         CASE ( 2 ,4)                   ! relative  vorticity
988            rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:)
989            zwzad(:,:) = 0.0_wp
990         CASE ( 3 ,5 )                  ! metric term
991            DO jj = jpjm1, 1, -1
992               DO ji = fs_jpim1, 1, -1  ! vector opt.
993                  vn_ad(ji+1,jj,jk) = vn_ad(ji+1,jj,jk)                   &
994                     &  + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) )    &
995                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
996                  vn_ad(ji  ,jj,jk) = vn_ad(ji  ,jj,jk)                   &
997                     &  + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) )    &
998                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
999                  un_ad(ji,jj+1,jk) = un_ad(ji,jj+1,jk)                   &
1000                     &  - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) )    &
1001                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
1002                  un_ad(ji,jj  ,jk) = un_ad(ji,jj  ,jk)                   &
1003                     &  - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) )    &
1004                     &                 * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
1005                  zwzad(ji,jj) = 0.0_wp
1006               END DO
1007            END DO
1008         END SELECT
1009         !                                             ! ===============
1010      END DO                                           !   End of slab
1011      !                                                ! ===============
1012      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
1013      CALL wrk_dealloc( jpi, jpj, zwxad, zwyad, zwzad )
1014      !
1015      IF( nn_timing == 1 )  CALL timing_stop('vor_ens_adj')
1016      !
1017   END SUBROUTINE vor_ens_adj
1018
1019
1020   SUBROUTINE vor_een_adj( kt, kvor, pua_ad, pva_ad )
1021      !!----------------------------------------------------------------------
1022      !!                ***  ROUTINE vor_een_adj  ***
1023      !!
1024      !! ** Purpose :   Compute the now total vorticity trend and add it to
1025      !!      the general trend of the momentum equation.
1026      !!
1027      !! ** Method  :   Trend evaluated using now fields (centered in time)
1028      !!      and the Arakawa and Lamb (19XX) flux form formulation : conserves
1029      !!      both the horizontal kinetic energy and the potential enstrophy
1030      !!      when horizontal divergence is zero.
1031      !!      The trend of the vorticity term is given by:
1032      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
1033      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
1034      !!      Add this trend to the general momentum trend (ua,va):
1035      !!          (ua,va) = (ua,va) + ( voru , vorv )
1036      !!
1037      !! ** Action : - Update (ua,va) with the now vorticity term trend
1038      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
1039      !!               and planetary vorticity trends) ('key_trddyn')
1040      !!
1041      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
1042      !!----------------------------------------------------------------------
1043      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
1044      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
1045         !                                                        ! =nrvm (relative vorticity or metric)
1046      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua_ad ! total u-trend
1047      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva_ad ! total v-trend
1048      !!
1049      INTEGER ::   ji, jj, jk          ! dummy loop indices
1050      REAL(wp) ::   zfac12             ! temporary scalars
1051      REAL(wp) ::   zuaad, zvaad       ! temporary scalars
1052      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz                    ! temporary 2D workspace
1053      REAL(wp), POINTER, DIMENSION(:,:) ::   ztnw, ztne, ztsw, ztse           ! temporary 3D workspace
1054      REAL(wp), POINTER, DIMENSION(:,:) ::   zwxad, zwyad, zwzad              ! temporary 2D workspace
1055      REAL(wp), POINTER, DIMENSION(:,:) ::   ztnwad, ztnead, ztswad, ztsead   ! temporary 3D workspace
1056      REAL(wp), POINTER, DIMENSION(:,:,:), SAVE ::   ze3f_ad
1057      !!----------------------------------------------------------------------
1058      !
1059      IF( nn_timing == 1 )  CALL timing_start('vor_een_adj')
1060      !
1061      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        )
1062      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )
1063      CALL wrk_alloc( jpi, jpj,      zwxad , zwyad , zwzad        )
1064      CALL wrk_alloc( jpi, jpj,      ztnwad, ztnead, ztswad, ztsead )
1065      CALL wrk_alloc( jpi, jpj, jpk,  ze3f_ad )
1066      ! local adjoint initailization
1067      zuaad = 0.0_wp ; zvaad = 0.0_wp
1068      zwx (:,:) = 0.0_wp ; zwy (:,:) = 0.0_wp ; zwz (:,:) = 0.0_wp
1069      ztnw(:,:) = 0.0_wp ; ztne(:,:) = 0.0_wp ; ztsw(:,:) = 0.0_wp ; ztse(:,:) = 0.0_wp
1070      zwxad (:,:) = 0.0_wp ; zwyad (:,:) = 0.0_wp ; zwzad (:,:) = 0.0_wp
1071      ztnwad(:,:) = 0.0_wp ; ztnead(:,:) = 0.0_wp ; ztswad(:,:) = 0.0_wp ; ztsead(:,:) = 0.0_wp
1072      ze3f_ad = 0._wp
1073
1074      IF( kt == nitend ) THEN
1075         IF(lwp) WRITE(numout,*)
1076         IF(lwp) WRITE(numout,*) 'dyn:vor_een_adj : vorticity term: energy and enstrophy conserving scheme'
1077         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
1078      ENDIF
1079
1080      IF( kt == nit000 ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
1081         DO jk = 1, jpk
1082            DO jj = 1, jpjm1
1083               DO ji = 1, jpim1
1084                  ze3f_ad(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
1085                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25_wp
1086                  IF( ze3f_ad(ji,jj,jk) /= 0.0_wp )   ze3f_ad(ji,jj,jk) = 1.0_wp / ze3f_ad(ji,jj,jk)
1087               END DO
1088            END DO
1089         END DO
1090         CALL lbc_lnk( ze3f_ad, 'F', 1._wp )
1091      ENDIF
1092
1093      ! Local constant initialization
1094      zfac12 = 1.0_wp / 12.0_wp
1095
1096!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
1097      !                                                ! ===============
1098      DO jk = 1, jpkm1                                 ! Horizontal slab
1099         !                                             ! ===============
1100
1101         ! Potential vorticity and horizontal fluxes  (Direct local variables init)
1102         ! -----------------------------------------
1103         SELECT CASE( kvor )      ! vorticity considered
1104         CASE ( 1 )
1105            zwz(:,:) =  ff(:,:) * ze3f_ad(:,:,jk)      ! planetary vorticity (Coriolis)
1106         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk) * ze3f_ad(:,:,jk)                ! relative  vorticity
1107         CASE ( 3 )                                                ! metric term
1108            DO jj = 1, jpjm1
1109               DO ji = 1, fs_jpim1   ! vector opt.
1110                  zwz(ji,jj) = ( ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )  &
1111                       &       - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) ) )&
1112                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1113               END DO
1114            END DO
1115            CALL lbc_lnk( zwz, 'F', 1._wp )
1116         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f_ad(:,:,jk)    ! total (relative + planetary vorticity)
1117         CASE ( 5 )                                                ! total (coriolis + metric)
1118            DO jj = 1, jpjm1
1119               DO ji = 1, fs_jpim1   ! vector opt.
1120                  zwz(ji,jj) = ( ff (ji,jj)                         &
1121                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
1122                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
1123                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )  &
1124                       &       ) * ze3f_ad(ji,jj,jk)
1125               END DO
1126            END DO
1127            CALL lbc_lnk( zwz, 'F', 1._wp )
1128         END SELECT
1129
1130         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
1131         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
1132
1133         ! Compute and add the vorticity term trend
1134         ! ----------------------------------------
1135         jj=2
1136         ztne(1,:)   = 0.0_wp ; ztnw(1,:)   = 0.0_wp ; ztse(1,:)   = 0.0_wp ; ztsw(1,:)   = 0.0_wp
1137         DO ji = 2, jpi
1138               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1139               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1140               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1141               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1142         END DO
1143         DO jj = 3, jpj
1144            DO ji = fs_2, jpi   ! vector opt.
1145               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1146               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1147               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1148               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1149            END DO
1150         END DO
1151
1152      ! ===================
1153      ! Adjoint counterpart
1154      ! ===================
1155
1156         DO jj = jpjm1, 2, -1
1157            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
1158               zuaad = zuaad + pua_ad(ji,jj,jk)
1159               zvaad = zvaad + pva_ad(ji,jj,jk)
1160
1161               zvaad = - zvaad * zfac12 / e2v(ji,jj)
1162               ztswad(ji  ,jj+1) = ztswad(ji  ,jj+1) + zvaad * zwx (ji-1,jj+1)
1163               zwxad (ji-1,jj+1) = zwxad (ji-1,jj+1) + zvaad * ztsw(ji  ,jj+1)
1164               ztsead(ji  ,jj+1) = ztsead(ji  ,jj+1) + zvaad * zwx (ji  ,jj+1)
1165               zwxad (ji  ,jj+1) = zwxad (ji  ,jj+1) + zvaad * ztse(ji  ,jj+1)
1166               ztnwad(ji  ,jj  ) = ztnwad(ji  ,jj  ) + zvaad * zwx (ji-1,jj  )
1167               zwxad (ji-1,jj  ) = zwxad (ji-1,jj  ) + zvaad * ztnw(ji  ,jj  )
1168               ztnead(ji  ,jj  ) = ztnead(ji  ,jj  ) + zvaad * zwx (ji  ,jj  )
1169               zwxad (ji  ,jj  ) = zwxad (ji  ,jj  ) + zvaad * ztne(ji  ,jj  )
1170               zvaad = 0.0_wp
1171
1172               zuaad = zuaad * zfac12 / e1u(ji,jj)
1173               ztnead(ji  ,jj  ) = ztnead(ji  ,jj  ) + zuaad * zwy (ji  ,jj  )
1174               zwyad (ji  ,jj  ) = zwyad (ji  ,jj  ) + zuaad * ztne(ji  ,jj  )
1175               ztnwad(ji+1,jj  ) = ztnwad(ji+1,jj  ) + zuaad * zwy (ji+1,jj  )
1176               zwyad (ji+1,jj  ) = zwyad (ji+1,jj  ) + zuaad * ztnw(ji+1,jj  )
1177               ztsead(ji  ,jj  ) = ztsead(ji  ,jj  ) + zuaad * zwy (ji  ,jj-1)
1178               zwyad (ji  ,jj-1) = zwyad (ji  ,jj-1) + zuaad * ztse(ji  ,jj  )
1179               ztswad(ji+1,jj  ) = ztswad(ji+1,jj  ) + zuaad * zwy (ji+1,jj-1)
1180               zwyad (ji+1,jj-1) = zwyad (ji+1,jj-1) + zuaad * ztsw(ji+1,jj  )
1181               zuaad = 0.0_wp
1182            END DO
1183         END DO
1184         DO jj = jpj, 3, -1
1185            DO ji = jpi, fs_2, -1   ! vector opt.
1186               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztswad(ji,jj)
1187               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj)
1188               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztswad(ji,jj)
1189               ztswad(ji  ,jj  ) = 0.0_wp
1190               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztsead(ji,jj)
1191               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztsead(ji,jj)
1192               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj)
1193               ztsead(ji,jj)     = 0.0_wp
1194               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj)
1195               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztnwad(ji,jj)
1196               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztnwad(ji,jj)
1197               ztnwad(ji  ,jj  ) = 0.0_wp
1198               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztnead(ji,jj)
1199               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztnead(ji,jj)
1200               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztnead(ji,jj)
1201               ztnead(ji,jj)     = 0.0_wp
1202            END DO
1203         END DO
1204         jj=2
1205         DO ji = jpi, 2, -1
1206               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztswad(ji,jj)
1207               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj)
1208               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztswad(ji,jj)
1209               ztswad(ji,jj) = 0.0_wp
1210               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztsead(ji,jj)
1211               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztsead(ji,jj)
1212               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj)
1213               ztsead(ji  ,jj  ) = 0.0_wp
1214               zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj)
1215               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztnwad(ji,jj)
1216               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztnwad(ji,jj)
1217               ztnwad(ji  ,jj ) = 0.0_wp
1218               zwzad (ji-1,jj  ) = zwzad(ji-1,jj  ) + ztnead(ji,jj)
1219               zwzad (ji  ,jj  ) = zwzad(ji  ,jj  ) + ztnead(ji,jj)
1220               zwzad (ji  ,jj-1) = zwzad(ji  ,jj-1) + ztnead(ji,jj)
1221               ztnead(ji  ,jj  ) = 0.0_wp
1222         END DO
1223         ztnead(1,:)   = 0.0_wp ; ztnwad(1,:)   = 0.0_wp
1224         ztsead(1,:)   = 0.0_wp ; ztswad(1,:)   = 0.0_wp
1225
1226         vn_ad(:,:,jk) = vn_ad(:,:,jk) + zwyad(:,:) * e1v(:,:) * fse3v(:,:,jk)
1227         un_ad(:,:,jk) = un_ad(:,:,jk) + zwxad(:,:) * e2u(:,:) * fse3u(:,:,jk)
1228         zwyad(:,:)    = 0.0_wp
1229         zwxad(:,:)    = 0.0_wp
1230
1231         ! Potential vorticity and horizontal fluxes
1232         ! -----------------------------------------
1233         SELECT CASE( kvor )      ! vorticity considered
1234         CASE ( 1 )
1235            zwzad(:,:) = 0.0_wp
1236         CASE ( 2 )
1237            rotn_ad(:,:,jk) = rotn_ad(:,:,jk) +  zwzad(:,:) * ze3f_ad(:,:,jk)
1238            zwzad(:,:)      = 0.0_wp
1239         CASE ( 3 )
1240            CALL lbc_lnk_adj( zwzad, 'F', 1._wp )                                                ! metric term
1241            DO jj = jpjm1, 1, -1
1242               DO ji = fs_jpim1, 1, -1   ! vector opt.
1243                  zwzad(ji  ,jj     ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1244                  vn_ad(ji+1,jj  ,jk) =  vn_ad(ji+1,jj  ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )
1245                  vn_ad(ji  ,jj  ,jk) =  vn_ad(ji  ,jj  ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )
1246                  un_ad(ji  ,jj+1,jk) =  un_ad(ji  ,jj+1,jk) - zwzad(ji,jj) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )
1247                  un_ad(ji  ,jj  ,jk) =  un_ad(ji  ,jj  ,jk) - zwzad(ji,jj) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )
1248                  zwzad(ji  ,jj     ) = 0.0_wp
1249               END DO
1250            END DO
1251         CASE ( 4 )
1252            rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) * ze3f_ad(:,:,jk)
1253            zwzad(:,:) = 0.0_wp
1254         CASE ( 5 )
1255            CALL lbc_lnk_adj( zwzad, 'F', 1._wp )                                               ! total (coriolis + metric)
1256            DO jj = jpjm1, 1, -1
1257               DO ji = fs_jpim1, 1, -1   ! vector opt.
1258                  zwzad(ji  ,jj     ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f_ad(ji,jj,jk)
1259                  vn_ad(ji+1,jj  ,jk) = vn_ad(ji+1,jj  ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )
1260                  vn_tl(ji  ,jj  ,jk) = vn_tl(ji  ,jj  ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )
1261                  un_ad(ji  ,jj+1,jk) = un_ad(ji  ,jj+1,jk) - zwzad(ji,jj) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )
1262                  un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) - zwzad(ji,jj) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )
1263                  zwzad(ji  ,jj     ) = 0.0_wp
1264               END DO
1265            END DO
1266         END SELECT
1267         !                                             ! ===============
1268      END DO                                           !   End of slab
1269      !                                                ! ===============
1270      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        )
1271      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )
1272      CALL wrk_dealloc( jpi, jpj,      zwxad , zwyad , zwzad        )
1273      CALL wrk_dealloc( jpi, jpj,      ztnwad, ztnead, ztswad, ztsead )
1274      CALL wrk_dealloc( jpi, jpj, jpk, ze3f_ad)
1275      !
1276      IF( nn_timing == 1 )  CALL timing_stop('vor_een_adj')
1277      !
1278   END SUBROUTINE vor_een_adj
1279
1280   SUBROUTINE dyn_vor_init_tam
1281      !!---------------------------------------------------------------------
1282      !!                  ***  ROUTINE vor_ctl_tam  ***
1283      !!
1284      !! ** Purpose :   Control the consistency between cpp options for
1285      !!      tracer advection schemes
1286      !!----------------------------------------------------------------------
1287      INTEGER ::   ioptio          ! temporary integer
1288      INTEGER ::   ji, jj, jk      ! dummy loop indices
1289      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
1290      !!----------------------------------------------------------------------
1291
1292      REWIND ( numnam )               ! Read Namelist namdyn_vor : Vorticity scheme options
1293      READ   ( numnam, namdyn_vor )
1294
1295      IF(lwp) THEN                    ! Namelist print
1296         WRITE(numout,*)
1297         WRITE(numout,*) 'dyn:vor_init_tam : vorticity term : read namelist and control the consistency'
1298         WRITE(numout,*) '~~~~~~~~~~~~~~~'
1299         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
1300         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
1301         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
1302         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
1303         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
1304      ENDIF
1305      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
1306      ! at angles with three ocean points and one land point
1307      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
1308         DO jk = 1, jpk
1309            DO jj = 2, jpjm1
1310               DO ji = 2, jpim1
1311                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
1312                      fmask(ji,jj,jk) = 1._wp
1313               END DO
1314            END DO
1315         END DO
1316          !
1317          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
1318          !
1319      ENDIF
1320      !
1321      ioptio = 0                     ! Control of vorticity scheme options
1322      IF( ln_dynvor_ene )   ioptio = ioptio + 1
1323      IF( ln_dynvor_ens )   ioptio = ioptio + 1
1324      IF( ln_dynvor_mix )   ioptio = ioptio + 1
1325      IF( ln_dynvor_een )   ioptio = ioptio + 1
1326      IF( lk_esopa      )   ioptio =          1
1327
1328      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
1329
1330      !                              ! Set nvor (type of scheme for vorticity)
1331      IF( ln_dynvor_ene )   nvor =  0
1332      IF( ln_dynvor_ens )   nvor =  1
1333      IF( ln_dynvor_mix )   nvor =  2
1334      IF( ln_dynvor_een )   nvor =  3
1335      IF( lk_esopa      )   nvor = -1
1336
1337      !                              ! Set ncor, nrvm, ntot (type of vorticity)
1338      IF(lwp) WRITE(numout,*)
1339      ncor = 1
1340      IF( ln_dynadv_vec ) THEN
1341         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
1342         nrvm = 2
1343         ntot = 4
1344      ELSE
1345         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
1346         nrvm = 3
1347         ntot = 5
1348      ENDIF
1349
1350      IF(lwp) THEN                   ! Print the choice
1351         WRITE(numout,*)
1352         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
1353         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
1354         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
1355         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
1356         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
1357      ENDIF
1358      !
1359   END SUBROUTINE dyn_vor_init_tam
1360
1361   SUBROUTINE dyn_vor_adj_tst( kumadt )
1362      !!-----------------------------------------------------------------------
1363      !!
1364      !!                  ***  ROUTINE dyn_adv_adj_tst ***
1365      !!
1366      !! ** Purpose : Test the adjoint routine.
1367      !!
1368      !! ** Method  : Verify the scalar product
1369      !!
1370      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
1371      !!
1372      !!              where  L   = tangent routine
1373      !!                     L^T = adjoint routine
1374      !!                     W   = diagonal matrix of scale factors
1375      !!                     dx  = input perturbation (random field)
1376      !!                     dy  = L dx
1377      !!
1378      !! ** Action  : Separate tests are applied for the following dx and dy:
1379      !!
1380      !!              1) dx = ( SSH ) and dy = ( SSH )
1381      !!
1382      !! History :
1383      !!        ! 08-08 (A. Vidard)
1384      !!-----------------------------------------------------------------------
1385      !! * Modules used
1386
1387      !! * Arguments
1388      INTEGER, INTENT(IN) :: &
1389         & kumadt             ! Output unit
1390
1391      INTEGER :: &
1392         & ji,    &        ! dummy loop indices
1393         & jj,    &
1394         & jk,    &
1395         & jt
1396      INTEGER, DIMENSION(jpi,jpj) :: &
1397         & iseed_2d        ! 2D seed for the random number generator
1398
1399      !! * Local declarations
1400      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1401         & zun_tlin,     & ! Tangent input:  now   u-velocity
1402         & zvn_tlin,     & ! Tangent input:  now   v-velocity
1403         & zrotn_tlin,   & ! Tangent input:  now   rot
1404         & zun_adout,    & ! Adjoint output: now   u-velocity
1405         & zvn_adout,    & ! Adjoint output: now   v-velocity
1406         & zrotn_adout,  & ! Adjoint output: now   rot
1407         & zua_adout,    & ! Tangent output: after u-velocity
1408         & zva_adout,    & ! Tangent output: after v-velocity
1409         & zua_tlin,     & ! Tangent output: after u-velocity
1410         & zva_tlin,     & ! Tangent output: after v-velocity
1411         & zua_tlout,    & ! Tangent output: after u-velocity
1412         & zva_tlout,    & ! Tangent output: after v-velocity
1413         & zua_adin,     & ! Tangent output: after u-velocity
1414         & zva_adin,     & ! Tangent output: after v-velocity
1415         & zau,          & ! 3D random field for rotn
1416         & zav,          & ! 3D random field for rotn
1417         & znu,          & ! 3D random field for u
1418         & znv             ! 3D random field for v
1419      REAL(KIND=wp) :: &
1420         & zsp1,         & ! scalar product involving the tangent routine
1421         & zsp1_1,       & !   scalar product components
1422         & zsp1_2,       &
1423         & zsp2,         & ! scalar product involving the adjoint routine
1424         & zsp2_1,       & !   scalar product components
1425         & zsp2_2,       &
1426         & zsp2_3,       &
1427         & zsp2_4,       &
1428         & zsp2_5
1429      CHARACTER(LEN=14) :: cl_name
1430
1431      ! Allocate memory
1432
1433      ALLOCATE( &
1434         & zun_tlin(jpi,jpj,jpk),     &
1435         & zvn_tlin(jpi,jpj,jpk),     &
1436         & zrotn_tlin(jpi,jpj,jpk),   &
1437         & zun_adout(jpi,jpj,jpk),    &
1438         & zvn_adout(jpi,jpj,jpk),    &
1439         & zrotn_adout(jpi,jpj,jpk),  &
1440         & zua_adout(jpi,jpj,jpk),    &
1441         & zva_adout(jpi,jpj,jpk),    &
1442         & zua_tlin(jpi,jpj,jpk),     &
1443         & zva_tlin(jpi,jpj,jpk),     &
1444         & zua_tlout(jpi,jpj,jpk),    &
1445         & zva_tlout(jpi,jpj,jpk),    &
1446         & zua_adin(jpi,jpj,jpk),     &
1447         & zva_adin(jpi,jpj,jpk),     &
1448         & zau(jpi,jpj,jpk),          &
1449         & zav(jpi,jpj,jpk),          &
1450         & znu(jpi,jpj,jpk),          &
1451         & znv(jpi,jpj,jpk)           &
1452         & )
1453
1454      ! init ntot parameter
1455      CALL dyn_vor_init_tam          ! initialisation & control of options
1456
1457      DO jt = 1, 2
1458         IF (jt == 1) nvor=1 ! enstrophy conserving scheme
1459         IF (jt == 2) nvor=3 ! energy and enstrophy conserving scheme
1460
1461         ! Initialize rotn
1462!AV: it calls cla_init the tries to reallocate already allocated arrays...nasty
1463!AV         CALL div_cur ( nit000 )
1464
1465         !==================================================================
1466         ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
1467         !    dy = ( hdivb_tl, hdivn_tl )
1468         !==================================================================
1469
1470         !--------------------------------------------------------------------
1471         ! Reset the tangent and adjoint variables
1472         !--------------------------------------------------------------------
1473
1474         zun_tlin(:,:,:) = 0.0_wp
1475         zvn_tlin(:,:,:) = 0.0_wp
1476         zrotn_tlin(:,:,:) = 0.0_wp
1477         zun_adout(:,:,:) = 0.0_wp
1478         zvn_adout(:,:,:) = 0.0_wp
1479         zrotn_adout(:,:,:) = 0.0_wp
1480         zua_tlout(:,:,:) = 0.0_wp
1481         zva_tlout(:,:,:) = 0.0_wp
1482         zua_adin(:,:,:) = 0.0_wp
1483         zva_adin(:,:,:) = 0.0_wp
1484         zua_adout(:,:,:) = 0.0_wp
1485         zva_adout(:,:,:) = 0.0_wp
1486         zua_tlin(:,:,:) = 0.0_wp
1487         zva_tlin(:,:,:) = 0.0_wp
1488         znu(:,:,:) = 0.0_wp
1489         znv(:,:,:) = 0.0_wp
1490         zau(:,:,:) = 0.0_wp
1491         zav(:,:,:) = 0.0_wp
1492
1493
1494         un_tl(:,:,:) = 0.0_wp
1495         vn_tl(:,:,:) = 0.0_wp
1496         ua_tl(:,:,:) = 0.0_wp
1497         va_tl(:,:,:) = 0.0_wp
1498         un_ad(:,:,:) = 0.0_wp
1499         vn_ad(:,:,:) = 0.0_wp
1500         ua_ad(:,:,:) = 0.0_wp
1501         va_ad(:,:,:) = 0.0_wp
1502         rotn_tl(:,:,:) = 0.0_wp
1503         rotn_ad(:,:,:) = 0.0_wp
1504
1505         !--------------------------------------------------------------------
1506         ! Initialize the tangent input with random noise: dx
1507         !--------------------------------------------------------------------
1508
1509         CALL grid_random(  znu, 'U', 0.0_wp, stdu )
1510         CALL grid_random(  znv, 'V', 0.0_wp, stdv )
1511         CALL grid_random(  zau, 'U', 0.0_wp, stdu )
1512         CALL grid_random(  zav, 'V', 0.0_wp, stdv )
1513
1514         DO jk = 1, jpk
1515            DO jj = nldj, nlej
1516               DO ji = nldi, nlei
1517                  zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
1518                  zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
1519                  zua_tlin(ji,jj,jk) = zau(ji,jj,jk)
1520                  zva_tlin(ji,jj,jk) = zav(ji,jj,jk)
1521               END DO
1522            END DO
1523         END DO
1524         un_tl(:,:,:) = zun_tlin(:,:,:)
1525         vn_tl(:,:,:) = zvn_tlin(:,:,:)
1526         ua_tl(:,:,:) = zua_tlin(:,:,:)
1527         va_tl(:,:,:) = zva_tlin(:,:,:)
1528
1529         ! initialize rotn_tl with noise
1530         CALL div_cur_tan ( nit000 )
1531
1532         DO jk = 1, jpk
1533            DO jj = nldj, nlej
1534               DO ji = nldi, nlei
1535                  zrotn_tlin(ji,jj,jk) = rotn_tl(ji,jj,jk)
1536               END DO
1537            END DO
1538         END DO
1539         rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1540
1541
1542         IF (nvor == 1 )  CALL vor_ens_tan( nit000, ntot, ua_tl, va_tl )
1543         IF (nvor == 3 )  CALL vor_een_tan( nit000, ntot, ua_tl, va_tl )
1544         zua_tlout(:,:,:) = ua_tl(:,:,:)
1545         zva_tlout(:,:,:) = va_tl(:,:,:)
1546
1547         !--------------------------------------------------------------------
1548         ! Initialize the adjoint variables: dy^* = W dy
1549         !--------------------------------------------------------------------
1550
1551         DO jk = 1, jpk
1552            DO jj = nldj, nlej
1553               DO ji = nldi, nlei
1554                  zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
1555                       &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
1556                       &               * umask(ji,jj,jk)
1557                  zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
1558                       &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
1559                       &               * vmask(ji,jj,jk)
1560               END DO
1561            END DO
1562         END DO
1563         !--------------------------------------------------------------------
1564         ! Compute the scalar product: ( L dx )^T W dy
1565         !--------------------------------------------------------------------
1566
1567         zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin )
1568         zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin )
1569         zsp1   = zsp1_1 + zsp1_2
1570
1571         !--------------------------------------------------------------------
1572         ! Call the adjoint routine: dx^* = L^T dy^*
1573         !--------------------------------------------------------------------
1574
1575         ua_ad(:,:,:) = zua_adin(:,:,:)
1576         va_ad(:,:,:) = zva_adin(:,:,:)
1577
1578
1579         IF (nvor == 1 )  CALL vor_ens_adj( nitend, ntot, ua_ad, va_ad )
1580         IF (nvor == 3 )  CALL vor_een_adj( nitend, ntot, ua_ad, va_ad )
1581         zun_adout(:,:,:)   = un_ad(:,:,:)
1582         zvn_adout(:,:,:)   = vn_ad(:,:,:)
1583         zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1584         zua_adout(:,:,:)   = ua_ad(:,:,:)
1585         zva_adout(:,:,:)   = va_ad(:,:,:)
1586
1587         zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
1588         zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
1589         zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1590         zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout )
1591         zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout )
1592         zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5
1593
1594         ! Compare the scalar products
1595
1596         ! 14 char:'12345678901234'
1597         IF (nvor == 1 )  cl_name = 'dynvor_adj ens'
1598         IF (nvor == 3 )  cl_name = 'dynvor_adj een'
1599
1600         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1601      END DO
1602
1603      DEALLOCATE( &
1604         & zun_tlin,     &
1605         & zvn_tlin,     &
1606         & zrotn_tlin,   &
1607         & zun_adout,    &
1608         & zvn_adout,    &
1609         & zrotn_adout,  &
1610         & zua_adout,    &
1611         & zva_adout,    &
1612         & zua_tlin,     &
1613         & zva_tlin,     &
1614         & zua_tlout,    &
1615         & zva_tlout,    &
1616         & zua_adin,     &
1617         & zva_adin,     &
1618         & zau,          &
1619         & zav,          &
1620         & znu,          &
1621         & znv           &
1622         & )
1623   END SUBROUTINE dyn_vor_adj_tst
1624   !!=============================================================================
1625#endif
1626END MODULE dynvor_tam
Note: See TracBrowser for help on using the repository browser.