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

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 5239

Last change on this file since 5239 was 5239, checked in by davestorkey, 9 years ago

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