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_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

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