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 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago

refer to ticket #798

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