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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3316

Last change on this file since 3316 was 3316, checked in by gm, 12 years ago

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