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.F90 in branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3098

Last change on this file since 3098 was 3098, checked in by cetlod, 12 years ago

dev_LOCEAN_CMCC_2011:minor bug correction
  • Property svn:keywords set to Id
File size: 40.2 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec) vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dyn_vor      : Update the momentum trend with the vorticity trend
21   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
22   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
23   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
24   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
25   !!   dyn_vor_init : set and control of the different vorticity option
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and tracers
28   USE dom_oce        ! ocean space and time domain
29   USE dommsk         ! ocean mask
30   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
31   USE trdmod         ! ocean dynamics trends
32   USE trdmod_oce     ! ocean variables trends
33   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
34   USE prtctl         ! Print control
35   USE in_out_manager ! I/O manager
36   USE lib_mpp
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   dyn_vor        ! routine called by step.F90
42   PUBLIC   dyn_vor_init   ! routine called by opa.F90
43
44   !                                             !!* Namelist namdyn_vor: vorticity term
45   LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
46   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
47   LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
48   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.   !: energy and enstrophy conserving scheme
49
50   INTEGER ::   nvor = 0   ! type of vorticity trend used
51   INTEGER ::   ncor = 1   ! coriolis
52   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
53   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE dyn_vor( kt )
66      !!----------------------------------------------------------------------
67      !!
68      !! ** Purpose :   compute the lateral ocean tracer physics.
69      !!
70      !! ** Action : - Update (ua,va) with the now vorticity term trend
71      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
72      !!               and planetary vorticity trends) ('key_trddyn')
73      !!----------------------------------------------------------------------
74      USE oce, ONLY:   tsa            ! tsa used as 2 3D workspace
75      !!
76      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
77      !
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
79      !!----------------------------------------------------------------------
80      !
81      IF( l_trddyn )   THEN
82         ztrdu => tsa(:,:,:,1) 
83         ztrdv => tsa(:,:,:,2) 
84      END IF
85      !
86      !                                          ! vorticity term
87      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
88      !
89      CASE ( -1 )                                      ! esopa: test all possibility with control print
90         CALL vor_ene( kt, ntot, ua, va )
91         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
92            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
93         CALL vor_ens( kt, ntot, ua, va )
94         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
95            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
96         CALL vor_mix( kt )
97         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
98            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
99         CALL vor_een( kt, ntot, ua, va )
100         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
101            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
102         !
103      CASE ( 0 )                                       ! energy conserving scheme
104         IF( l_trddyn )   THEN
105            ztrdu(:,:,:) = ua(:,:,:)
106            ztrdv(:,:,:) = va(:,:,:)
107            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
108            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
109            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
110            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
111            ztrdu(:,:,:) = ua(:,:,:)
112            ztrdv(:,:,:) = va(:,:,:)
113            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
114            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
115            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
116            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
117            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
118         ELSE
119            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
120         ENDIF
121         !
122      CASE ( 1 )                                       ! enstrophy conserving scheme
123         IF( l_trddyn )   THEN   
124            ztrdu(:,:,:) = ua(:,:,:)
125            ztrdv(:,:,:) = va(:,:,:)
126            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
127            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
128            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
129            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
130            ztrdu(:,:,:) = ua(:,:,:)
131            ztrdv(:,:,:) = va(:,:,:)
132            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
133            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
134            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
135            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
136            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
137         ELSE
138            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
139         ENDIF
140         !
141      CASE ( 2 )                                       ! mixed ene-ens scheme
142         IF( l_trddyn )   THEN
143            ztrdu(:,:,:) = ua(:,:,:)
144            ztrdv(:,:,:) = va(:,:,:)
145            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
146            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
147            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
148            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
149            ztrdu(:,:,:) = ua(:,:,:)
150            ztrdv(:,:,:) = va(:,:,:)
151            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
152            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
153            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
154            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
155            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
156         ELSE
157            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
158         ENDIF
159         !
160      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
161         IF( l_trddyn )   THEN
162            ztrdu(:,:,:) = ua(:,:,:)
163            ztrdv(:,:,:) = va(:,:,:)
164            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
165            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
166            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
167            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
168            ztrdu(:,:,:) = ua(:,:,:)
169            ztrdv(:,:,:) = va(:,:,:)
170            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
171            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
172            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
173            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
174            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
175         ELSE
176            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
177         ENDIF
178         !
179      END SELECT
180      !
181      !                       ! print sum trends (used for debugging)
182      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
183         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
184      !
185   END SUBROUTINE dyn_vor
186
187
188   SUBROUTINE vor_ene( kt, kvor, pua, pva )
189      !!----------------------------------------------------------------------
190      !!                  ***  ROUTINE vor_ene  ***
191      !!
192      !! ** Purpose :   Compute the now total vorticity trend and add it to
193      !!      the general trend of the momentum equation.
194      !!
195      !! ** Method  :   Trend evaluated using now fields (centered in time)
196      !!      and the Sadourny (1975) flux form formulation : conserves the
197      !!      horizontal kinetic energy.
198      !!      The trend of the vorticity term is given by:
199      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
200      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
201      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
202      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
203      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
204      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
205      !!      Add this trend to the general momentum trend (ua,va):
206      !!          (ua,va) = (ua,va) + ( voru , vorv )
207      !!
208      !! ** Action : - Update (ua,va) with the now vorticity term trend
209      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
210      !!               and planetary vorticity trends) ('key_trddyn')
211      !!
212      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
213      !!----------------------------------------------------------------------
214      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
215      USE wrk_nemo, ONLY:   zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3     ! 2D workspace
216      !
217      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
218      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
219      !                                                           ! =nrvm (relative vorticity or metric)
220      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
221      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
222      !
223      INTEGER  ::   ji, jj, jk   ! dummy loop indices
224      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
225      !!----------------------------------------------------------------------
226
227      IF( wrk_in_use(2, 1,2,3) ) THEN
228         CALL ctl_stop('dyn:vor_ene: requested workspace arrays unavailable')   ;   RETURN
229      ENDIF
230
231      IF( kt == nit000 ) THEN
232         IF(lwp) WRITE(numout,*)
233         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
234         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
235      ENDIF
236
237      zfact2 = 0.5 * 0.5      ! Local constant initialization
238
239!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
240      !                                                ! ===============
241      DO jk = 1, jpkm1                                 ! Horizontal slab
242         !                                             ! ===============
243         !
244         ! Potential vorticity and horizontal fluxes
245         ! -----------------------------------------
246         SELECT CASE( kvor )      ! vorticity considered
247         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
248         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
249         CASE ( 3 )                                                ! metric term
250            DO jj = 1, jpjm1
251               DO ji = 1, fs_jpim1   ! vector opt.
252                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
253                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
254                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
255               END DO
256            END DO
257         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
258         CASE ( 5 )                                                ! total (coriolis + metric)
259            DO jj = 1, jpjm1
260               DO ji = 1, fs_jpim1   ! vector opt.
261                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
262                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
263                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
264                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
265                       &       )
266               END DO
267            END DO
268         END SELECT
269
270         IF( ln_sco ) THEN
271            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
272            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
273            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
274         ELSE
275            zwx(:,:) = e2u(:,:) * un(:,:,jk)
276            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
277         ENDIF
278
279         ! Compute and add the vorticity term trend
280         ! ----------------------------------------
281         DO jj = 2, jpjm1
282            DO ji = fs_2, fs_jpim1   ! vector opt.
283               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
284               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
285               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
286               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
287               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
288               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
289            END DO 
290         END DO 
291         !                                             ! ===============
292      END DO                                           !   End of slab
293      !                                                ! ===============
294      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn:vor_ene: failed to release workspace arrays')
295      !
296   END SUBROUTINE vor_ene
297
298
299   SUBROUTINE vor_mix( kt )
300      !!----------------------------------------------------------------------
301      !!                 ***  ROUTINE vor_mix  ***
302      !!
303      !! ** Purpose :   Compute the now total vorticity trend and add it to
304      !!      the general trend of the momentum equation.
305      !!
306      !! ** Method  :   Trend evaluated using now fields (centered in time)
307      !!      Mixte formulation : conserves the potential enstrophy of a hori-
308      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
309      !!      ticity term and the horizontal kinetic energy for (f x uh), the
310      !!      coriolis term. the now trend of the vorticity term is given by:
311      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
312      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
313      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
314      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
315      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
316      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
317      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
318      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
319      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
320      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
321      !!      Add this now trend to the general momentum trend (ua,va):
322      !!          (ua,va) = (ua,va) + ( voru , vorv )
323      !!
324      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
325      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
326      !!               and planetary vorticity trends) ('key_trddyn')
327      !!
328      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
329      !!----------------------------------------------------------------------
330      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
331      USE wrk_nemo, ONLY:   zwx => wrk_2d_4 , zwy => wrk_2d_5 , zwz => wrk_2d_6 , zww => wrk_2d_7   ! 2D workspace
332      !
333      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
334      !
335      INTEGER  ::   ji, jj, jk   ! dummy loop indices
336      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
337      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
338      !!----------------------------------------------------------------------
339
340      IF( wrk_in_use(2, 4,5,6,7) ) THEN
341         CALL ctl_stop('dyn:vor_mix: requested workspace arrays unavailable')   ;   RETURN
342      ENDIF
343
344      IF( kt == nit000 ) THEN
345         IF(lwp) WRITE(numout,*)
346         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
347         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
348      ENDIF
349
350      zfact1 = 0.5 * 0.25      ! Local constant initialization
351      zfact2 = 0.5 * 0.5
352
353!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
354      !                                                ! ===============
355      DO jk = 1, jpkm1                                 ! Horizontal slab
356         !                                             ! ===============
357         !
358         ! Relative and planetary potential vorticity and horizontal fluxes
359         ! ----------------------------------------------------------------
360         IF( ln_sco ) THEN       
361            IF( ln_dynadv_vec ) THEN
362               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
363            ELSE                       
364               DO jj = 1, jpjm1
365                  DO ji = 1, fs_jpim1   ! vector opt.
366                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
367                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
368                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
369                  END DO
370               END DO
371            ENDIF
372            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
373            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
374            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
375         ELSE
376            IF( ln_dynadv_vec ) THEN
377               zww(:,:) = rotn(:,:,jk)
378            ELSE                       
379               DO jj = 1, jpjm1
380                  DO ji = 1, fs_jpim1   ! vector opt.
381                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
382                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
383                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
384                  END DO
385               END DO
386            ENDIF
387            zwz(:,:) = ff (:,:)
388            zwx(:,:) = e2u(:,:) * un(:,:,jk)
389            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
390         ENDIF
391
392         ! Compute and add the vorticity term trend
393         ! ----------------------------------------
394         DO jj = 2, jpjm1
395            DO ji = fs_2, fs_jpim1   ! vector opt.
396               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
397               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
398               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
399               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
400               ! enstrophy conserving formulation for relative vorticity term
401               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
402               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
403               ! energy conserving formulation for planetary vorticity term
404               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
405               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
406               ! mixed vorticity trend added to the momentum trends
407               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
408               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
409            END DO 
410         END DO 
411         !                                             ! ===============
412      END DO                                           !   End of slab
413      !                                                ! ===============
414      IF( wrk_not_released(2, 4,5,6,7) )   CALL ctl_stop('dyn:vor_mix: failed to release workspace arrays')
415      !
416   END SUBROUTINE vor_mix
417
418
419   SUBROUTINE vor_ens( kt, kvor, pua, pva )
420      !!----------------------------------------------------------------------
421      !!                ***  ROUTINE vor_ens  ***
422      !!
423      !! ** Purpose :   Compute the now total vorticity trend and add it to
424      !!      the general trend of the momentum equation.
425      !!
426      !! ** Method  :   Trend evaluated using now fields (centered in time)
427      !!      and the Sadourny (1975) flux FORM formulation : conserves the
428      !!      potential enstrophy of a horizontally non-divergent flow. the
429      !!      trend of the vorticity term is given by:
430      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
431      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
432      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
433      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
434      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
435      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
436      !!      Add this trend to the general momentum trend (ua,va):
437      !!          (ua,va) = (ua,va) + ( voru , vorv )
438      !!
439      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
440      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
441      !!               and planetary vorticity trends) ('key_trddyn')
442      !!
443      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
444      !!----------------------------------------------------------------------
445      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
446      USE wrk_nemo, ONLY:   zwx => wrk_2d_4, zwy => wrk_2d_5, zwz => wrk_2d_6    ! 2D workspace
447      !
448      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
449      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
450         !                                                        ! =nrvm (relative vorticity or metric)
451      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
452      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
453      !
454      INTEGER  ::   ji, jj, jk           ! dummy loop indices
455      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
456      !!----------------------------------------------------------------------
457     
458      IF( wrk_in_use(2, 4,5,6) ) THEN
459         CALL ctl_stop('dyn:vor_ens: requested workspace arrays unavailable')   ;   RETURN
460      END IF
461
462      IF( kt == nit000 ) THEN
463         IF(lwp) WRITE(numout,*)
464         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
465         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
466      ENDIF
467
468      zfact1 = 0.5 * 0.25      ! Local constant initialization
469
470!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
471      !                                                ! ===============
472      DO jk = 1, jpkm1                                 ! Horizontal slab
473         !                                             ! ===============
474         !
475         ! Potential vorticity and horizontal fluxes
476         ! -----------------------------------------
477         SELECT CASE( kvor )      ! vorticity considered
478         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
479         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
480         CASE ( 3 )                                                ! metric term
481            DO jj = 1, jpjm1
482               DO ji = 1, fs_jpim1   ! vector opt.
483                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
484                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
485                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
486               END DO
487            END DO
488         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
489         CASE ( 5 )                                                ! total (coriolis + metric)
490            DO jj = 1, jpjm1
491               DO ji = 1, fs_jpim1   ! vector opt.
492                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
493                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
494                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
495                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
496                       &       )
497               END DO
498            END DO
499         END SELECT
500         !
501         IF( ln_sco ) THEN
502            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
503               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
504                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
505                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
506                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
507               END DO
508            END DO
509         ELSE
510            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
511               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
512                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
513                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
514               END DO
515            END DO
516         ENDIF
517         !
518         ! Compute and add the vorticity term trend
519         ! ----------------------------------------
520         DO jj = 2, jpjm1
521            DO ji = fs_2, fs_jpim1   ! vector opt.
522               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
523                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
524               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
525                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
526               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
527               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
528            END DO 
529         END DO 
530         !                                             ! ===============
531      END DO                                           !   End of slab
532      !                                                ! ===============
533      IF( wrk_not_released(2, 4,5,6) )   CALL ctl_stop('dyn:vor_ens: failed to release workspace arrays')
534      !
535   END SUBROUTINE vor_ens
536
537
538   SUBROUTINE vor_een( kt, kvor, pua, pva )
539      !!----------------------------------------------------------------------
540      !!                ***  ROUTINE vor_een  ***
541      !!
542      !! ** Purpose :   Compute the now total vorticity trend and add it to
543      !!      the general trend of the momentum equation.
544      !!
545      !! ** Method  :   Trend evaluated using now fields (centered in time)
546      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
547      !!      both the horizontal kinetic energy and the potential enstrophy
548      !!      when horizontal divergence is zero (see the NEMO documentation)
549      !!      Add this trend to the general momentum trend (ua,va).
550      !!
551      !! ** Action : - Update (ua,va) with the now vorticity term trend
552      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
553      !!               and planetary vorticity trends) ('key_trddyn')
554      !!
555      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
556      !!----------------------------------------------------------------------
557      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
558      USE wrk_nemo, ONLY:   zwx  => wrk_2d_1 , zwy  => wrk_2d_2 ,  zwz => wrk_2d_3     ! 2D workspace
559      USE wrk_nemo, ONLY:   ztnw => wrk_2d_4 , ztne => wrk_2d_5 
560      USE wrk_nemo, ONLY:   ztsw => wrk_2d_6 , ztse => wrk_2d_7
561#if defined key_vvl
562      USE wrk_nemo, ONLY:   ze3f => wrk_3d_1                                           ! 3D workspace (lk_vvl=T)
563#endif
564      !
565      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
566      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
567      !                                                           ! =nrvm (relative vorticity or metric)
568      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
569      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
570      !!
571      INTEGER  ::   ji, jj, jk         ! dummy loop indices
572      INTEGER  ::   ierr               ! local integer
573      REAL(wp) ::   zfac12, zua, zva   ! local scalars
574#if ! defined key_vvl
575      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE ::   ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
576#endif
577      !!----------------------------------------------------------------------
578
579      IF( wrk_in_use(2, 1,2,3,4,5,6,7) .OR. wrk_in_use(3, 1) ) THEN
580         CALL ctl_stop('dyn:vor_een: requested workspace arrays unavailable')   ;   RETURN
581      ENDIF
582
583      IF( kt == nit000 ) THEN
584         IF(lwp) WRITE(numout,*)
585         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
586         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
587         IF( .NOT.lk_vvl ) THEN
588            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
589            IF( lk_mpp    )   CALL mpp_sum ( ierr )
590            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
591         ENDIF
592      ENDIF
593
594      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
595         DO jk = 1, jpk
596            DO jj = 1, jpjm1
597               DO ji = 1, jpim1
598                  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)   &
599                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25
600                  IF( ze3f(ji,jj,jk) /= 0._wp )   ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
601               END DO
602            END DO
603         END DO
604         CALL lbc_lnk( ze3f, 'F', 1. )
605      ENDIF
606
607      zfac12 = 1._wp / 12._wp    ! Local constant initialization
608
609     
610!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
611      !                                                ! ===============
612      DO jk = 1, jpkm1                                 ! Horizontal slab
613         !                                             ! ===============
614         
615         ! Potential vorticity and horizontal fluxes
616         ! -----------------------------------------
617         SELECT CASE( kvor )      ! vorticity considered
618         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
619            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
620         CASE ( 2 )                                                ! relative  vorticity
621            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
622         CASE ( 3 )                                                ! metric term
623            DO jj = 1, jpjm1
624               DO ji = 1, fs_jpim1   ! vector opt.
625                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
626                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
627                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
628               END DO
629            END DO
630            CALL lbc_lnk( zwz, 'F', 1. )
631        CASE ( 4 )                                                ! total (relative + planetary vorticity)
632            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
633         CASE ( 5 )                                                ! total (coriolis + metric)
634            DO jj = 1, jpjm1
635               DO ji = 1, fs_jpim1   ! vector opt.
636                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
637                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
638                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
639                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
640                       &       ) * ze3f(ji,jj,jk)
641               END DO
642            END DO
643            CALL lbc_lnk( zwz, 'F', 1. )
644         END SELECT
645
646         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
647         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
648
649         ! Compute and add the vorticity term trend
650         ! ----------------------------------------
651         jj = 2
652         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
653         DO ji = 2, jpi   
654               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
655               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
656               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
657               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
658         END DO
659         DO jj = 3, jpj
660            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
661               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
662               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
663               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
664               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
665            END DO
666         END DO
667         DO jj = 2, jpjm1
668            DO ji = fs_2, fs_jpim1   ! vector opt.
669               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
670                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
671               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
672                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
673               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
674               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
675            END DO 
676         END DO 
677         !                                             ! ===============
678      END DO                                           !   End of slab
679      !                                                ! ===============
680      IF( wrk_not_released(2, 1,2,3,4,5,6,7) .OR.   &
681          wrk_not_released(3, 1)             )   CALL ctl_stop('dyn:vor_een: failed to release workspace arrays')
682      !
683   END SUBROUTINE vor_een
684
685
686   SUBROUTINE dyn_vor_init
687      !!---------------------------------------------------------------------
688      !!                  ***  ROUTINE dyn_vor_init  ***
689      !!
690      !! ** Purpose :   Control the consistency between cpp options for
691      !!              tracer advection schemes
692      !!----------------------------------------------------------------------
693      INTEGER ::   ioptio          ! local integer
694      INTEGER ::   ji, jj, jk      ! dummy loop indices
695      !!
696      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
697      !!----------------------------------------------------------------------
698
699      REWIND ( numnam )               ! Read Namelist namdyn_vor : Vorticity scheme options
700      READ   ( numnam, namdyn_vor )
701
702      IF(lwp) THEN                    ! Namelist print
703         WRITE(numout,*)
704         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
705         WRITE(numout,*) '~~~~~~~~~~~~'
706         WRITE(numout,*) '        Namelist namdyn_vor : oice of the vorticity term scheme'
707         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
708         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
709         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
710         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
711      ENDIF
712
713      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
714      ! at angles with three ocean points and one land point
715      IF( ln_vorlat .AND. (ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix) ) THEN
716          DO jj = 2, jpjm1
717            DO ji = 2, jpim1
718               IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
719                   fmask(ji,jj,jk) = 1._wp
720            END DO
721          END DO
722          !
723          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
724          !
725      ENDIF
726
727      ioptio = 0                     ! Control of vorticity scheme options
728      IF( ln_dynvor_ene )   ioptio = ioptio + 1
729      IF( ln_dynvor_ens )   ioptio = ioptio + 1
730      IF( ln_dynvor_mix )   ioptio = ioptio + 1
731      IF( ln_dynvor_een )   ioptio = ioptio + 1
732      IF( lk_esopa      )   ioptio =          1
733
734      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
735
736      !                              ! Set nvor (type of scheme for vorticity)
737      IF( ln_dynvor_ene )   nvor =  0
738      IF( ln_dynvor_ens )   nvor =  1
739      IF( ln_dynvor_mix )   nvor =  2
740      IF( ln_dynvor_een )   nvor =  3
741      IF( lk_esopa      )   nvor = -1
742     
743      !                              ! Set ncor, nrvm, ntot (type of vorticity)
744      IF(lwp) WRITE(numout,*)
745      ncor = 1
746      IF( ln_dynadv_vec ) THEN     
747         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
748         nrvm = 2
749         ntot = 4
750      ELSE                       
751         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
752         nrvm = 3
753         ntot = 5
754      ENDIF
755     
756      IF(lwp) THEN                   ! Print the choice
757         WRITE(numout,*)
758         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
759         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
760         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
761         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
762         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
763      ENDIF
764      !
765   END SUBROUTINE dyn_vor_init
766
767   !!==============================================================================
768END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.