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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 5737

Last change on this file since 5737 was 5737, checked in by gm, 9 years ago

#1593: LDF-ADV, step I: Phasing of horizontal scale factors correct 2

  • Property svn:keywords set to Id
File size: 40.6 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   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
54
55   INTEGER ::   nvor = 0   ! type of vorticity trend used
56   INTEGER ::   ncor = 1   ! coriolis
57   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
58   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
59
60   !! * Substitutions
61#  include "domzgr_substitute.h90"
62#  include "vectopt_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE dyn_vor( kt )
71      !!----------------------------------------------------------------------
72      !!
73      !! ** Purpose :   compute the lateral ocean tracer physics.
74      !!
75      !! ** Action : - Update (ua,va) with the now vorticity term trend
76      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
77      !!               and planetary vorticity trends) and send them to trd_dyn
78      !!               for futher diagnostics (l_trddyn=T)
79      !!----------------------------------------------------------------------
80      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
81      !
82      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
83      !!----------------------------------------------------------------------
84      !
85      IF( nn_timing == 1 )  CALL timing_start('dyn_vor')
86      !
87      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
88      !
89      !                                          ! vorticity term
90      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
91      !
92      CASE ( -1 )                                      ! esopa: test all possibility with control print
93         CALL vor_ene( kt, ntot, ua, va )
94         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
95            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
96         CALL vor_ens( kt, ntot, ua, va )
97         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
98            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
99         CALL vor_mix( kt )
100         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
101            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
102         CALL vor_een( kt, ntot, ua, va )
103         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
104            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
105         !
106      CASE ( 0 )                                       ! energy conserving scheme
107         IF( l_trddyn )   THEN
108            ztrdu(:,:,:) = ua(:,:,:)
109            ztrdv(:,:,:) = va(:,:,:)
110            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
111            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
112            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
113            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
114            ztrdu(:,:,:) = ua(:,:,:)
115            ztrdv(:,:,:) = va(:,:,:)
116            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
117            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
118            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
119            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
120         ELSE
121            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
122         ENDIF
123         !
124      CASE ( 1 )                                       ! enstrophy conserving scheme
125         IF( l_trddyn )   THEN   
126            ztrdu(:,:,:) = ua(:,:,:)
127            ztrdv(:,:,:) = va(:,:,:)
128            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
129            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
130            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
131            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
132            ztrdu(:,:,:) = ua(:,:,:)
133            ztrdv(:,:,:) = va(:,:,:)
134            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
135            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
136            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
137            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
138         ELSE
139            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
140         ENDIF
141         !
142      CASE ( 2 )                                       ! mixed ene-ens scheme
143         IF( l_trddyn )   THEN
144            ztrdu(:,:,:) = ua(:,:,:)
145            ztrdv(:,:,:) = va(:,:,:)
146            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
147            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
148            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
149            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
150            ztrdu(:,:,:) = ua(:,:,:)
151            ztrdv(:,:,:) = va(:,:,:)
152            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
153            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
154            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
155            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
156         ELSE
157            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
158         ENDIF
159         !
160      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
161         IF( l_trddyn )   THEN
162            ztrdu(:,:,:) = ua(:,:,:)
163            ztrdv(:,:,:) = va(:,:,:)
164            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
165            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
166            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
167            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
168            ztrdu(:,:,:) = ua(:,:,:)
169            ztrdv(:,:,:) = va(:,:,:)
170            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
171            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
172            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
173            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
174         ELSE
175            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
176         ENDIF
177         !
178      END SELECT
179      !
180      !                       ! print sum trends (used for debugging)
181      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
182         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
183      !
184      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
185      !
186      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor')
187      !
188   END SUBROUTINE dyn_vor
189
190
191   SUBROUTINE vor_ene( kt, kvor, pua, pva )
192      !!----------------------------------------------------------------------
193      !!                  ***  ROUTINE vor_ene  ***
194      !!
195      !! ** Purpose :   Compute the now total vorticity trend and add it to
196      !!      the general trend of the momentum equation.
197      !!
198      !! ** Method  :   Trend evaluated using now fields (centered in time)
199      !!      and the Sadourny (1975) flux form formulation : conserves the
200      !!      horizontal kinetic energy.
201      !!      The trend of the vorticity term is given by:
202      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
203      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
204      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
205      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
206      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
207      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
208      !!      Add this trend to the general momentum trend (ua,va):
209      !!          (ua,va) = (ua,va) + ( voru , vorv )
210      !!
211      !! ** Action : - Update (ua,va) with the now vorticity term trend
212      !!
213      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
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      !                                                ! ===============
239      DO jk = 1, jpkm1                                 ! Horizontal slab
240         !                                             ! ===============
241         !
242         ! Potential vorticity and horizontal fluxes
243         ! -----------------------------------------
244         SELECT CASE( kvor )      ! vorticity considered
245         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
246         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
247         CASE ( 3 )                                                ! metric term
248            DO jj = 1, jpjm1
249               DO ji = 1, fs_jpim1   ! vector opt.
250                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
251                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
252                       &     * 0.5 * r1_e1e2f(ji,jj)
253               END DO
254            END DO
255         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
256         CASE ( 5 )                                                ! total (coriolis + metric)
257            DO jj = 1, jpjm1
258               DO ji = 1, fs_jpim1   ! vector opt.
259                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
260                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
261                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
262                       &       * 0.5 * r1_e1e2f(ji,jj)                                              &
263                       &       )
264               END DO
265            END DO
266         END SELECT
267
268         IF( ln_sco ) THEN
269            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
270            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
271            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
272         ELSE
273            zwx(:,:) = e2u(:,:) * un(:,:,jk)
274            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
275         ENDIF
276
277         ! Compute and add the vorticity term trend
278         ! ----------------------------------------
279         DO jj = 2, jpjm1
280            DO ji = fs_2, fs_jpim1   ! vector opt.
281               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
282               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
283               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
284               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
285               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
286               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
287            END DO 
288         END DO 
289         !                                             ! ===============
290      END DO                                           !   End of slab
291      !                                                ! ===============
292      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
293      !
294      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
295      !
296   END SUBROUTINE vor_ene
297
298
299   SUBROUTINE vor_mix( kt )
300      !!----------------------------------------------------------------------
301      !!                 ***  ROUTINE vor_mix  ***
302      !!
303      !! ** Purpose :   Compute the now total vorticity trend and add it to
304      !!      the general trend of the momentum equation.
305      !!
306      !! ** Method  :   Trend evaluated using now fields (centered in time)
307      !!      Mixte formulation : conserves the potential enstrophy of a hori-
308      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
309      !!      ticity term and the horizontal kinetic energy for (f x uh), the
310      !!      coriolis term. the now trend of the vorticity term is given by:
311      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
312      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
313      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
314      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
315      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
316      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
317      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
318      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
319      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
320      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
321      !!      Add this now trend to the general momentum trend (ua,va):
322      !!          (ua,va) = (ua,va) + ( voru , vorv )
323      !!
324      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
325      !!
326      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
327      !!----------------------------------------------------------------------
328      !
329      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
330      !
331      INTEGER  ::   ji, jj, jk   ! dummy loop indices
332      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
333      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
334      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
335      !!----------------------------------------------------------------------
336      !
337      IF( nn_timing == 1 )  CALL timing_start('vor_mix')
338      !
339      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) 
340      !
341      IF( kt == nit000 ) THEN
342         IF(lwp) WRITE(numout,*)
343         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
344         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
345      ENDIF
346
347      zfact1 = 0.5 * 0.25      ! Local constant initialization
348      zfact2 = 0.5 * 0.5
349
350!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
351      !                                                ! ===============
352      DO jk = 1, jpkm1                                 ! Horizontal slab
353         !                                             ! ===============
354         !
355         ! Relative and planetary potential vorticity and horizontal fluxes
356         ! ----------------------------------------------------------------
357         IF( ln_sco ) THEN       
358            IF( ln_dynadv_vec ) THEN
359               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
360            ELSE                       
361               DO jj = 1, jpjm1
362                  DO ji = 1, fs_jpim1   ! vector opt.
363                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
364                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
365                        &       * 0.5 / ( e1e2f (ji,jj) * fse3f(ji,jj,jk) )
366                  END DO
367               END DO
368            ENDIF
369            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
370            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
371            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
372         ELSE
373            IF( ln_dynadv_vec ) THEN
374               zww(:,:) = rotn(:,:,jk)
375            ELSE                       
376               DO jj = 1, jpjm1
377                  DO ji = 1, fs_jpim1   ! vector opt.
378                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
379                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
380                        &       * 0.5 * r1_e1e2f(ji,jj)
381                  END DO
382               END DO
383            ENDIF
384            zwz(:,:) = ff (:,:)
385            zwx(:,:) = e2u(:,:) * un(:,:,jk)
386            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
387         ENDIF
388
389         ! Compute and add the vorticity term trend
390         ! ----------------------------------------
391         DO jj = 2, jpjm1
392            DO ji = fs_2, fs_jpim1   ! vector opt.
393               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
394               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
395               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
396               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
397               ! enstrophy conserving formulation for relative vorticity term
398               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
399               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
400               ! energy conserving formulation for planetary vorticity term
401               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
402               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
403               ! mixed vorticity trend added to the momentum trends
404               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
405               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
406            END DO 
407         END DO 
408         !                                             ! ===============
409      END DO                                           !   End of slab
410      !                                                ! ===============
411      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) 
412      !
413      IF( nn_timing == 1 )  CALL timing_stop('vor_mix')
414      !
415   END SUBROUTINE vor_mix
416
417
418   SUBROUTINE vor_ens( kt, kvor, pua, pva )
419      !!----------------------------------------------------------------------
420      !!                ***  ROUTINE vor_ens  ***
421      !!
422      !! ** Purpose :   Compute the now total vorticity trend and add it to
423      !!      the general trend of the momentum equation.
424      !!
425      !! ** Method  :   Trend evaluated using now fields (centered in time)
426      !!      and the Sadourny (1975) flux FORM formulation : conserves the
427      !!      potential enstrophy of a horizontally non-divergent flow. the
428      !!      trend of the vorticity term is given by:
429      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
430      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
431      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
432      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
433      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
434      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
435      !!      Add this trend to the general momentum trend (ua,va):
436      !!          (ua,va) = (ua,va) + ( voru , vorv )
437      !!
438      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
439      !!
440      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
441      !!----------------------------------------------------------------------
442      !
443      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
444      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
445         !                                                        ! =nrvm (relative vorticity or metric)
446      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
447      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
448      !
449      INTEGER  ::   ji, jj, jk           ! dummy loop indices
450      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
451      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
452      !!----------------------------------------------------------------------
453      !
454      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
455      !
456      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
457      !
458      IF( kt == nit000 ) THEN
459         IF(lwp) WRITE(numout,*)
460         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
461         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
462      ENDIF
463
464      zfact1 = 0.5 * 0.25      ! Local constant initialization
465
466!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
467      !                                                ! ===============
468      DO jk = 1, jpkm1                                 ! Horizontal slab
469         !                                             ! ===============
470         !
471         ! Potential vorticity and horizontal fluxes
472         ! -----------------------------------------
473         SELECT CASE( kvor )      ! vorticity considered
474         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
475         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
476         CASE ( 3 )                                                ! metric term
477            DO jj = 1, jpjm1
478               DO ji = 1, fs_jpim1   ! vector opt.
479                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
480                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
481                       &     * 0.5 * r1_e1e2f(ji,jj)
482               END DO
483            END DO
484         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
485         CASE ( 5 )                                                ! total (coriolis + metric)
486            DO jj = 1, jpjm1
487               DO ji = 1, fs_jpim1   ! vector opt.
488                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
489                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
490                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
491                       &       * 0.5 * r1_e1e2f(ji,jj)                                                &
492                       &       )
493               END DO
494            END DO
495         END SELECT
496         !
497         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
498            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
499            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
500            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
501         ELSE
502            zwx(:,:) = e2u(:,:) * un(:,:,jk)
503            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
504         ENDIF
505         !                                   !==  compute and add the vorticity term trend  =!
506         DO jj = 2, jpjm1
507            DO ji = fs_2, fs_jpim1   ! vector opt.
508               zuav = zfact1 * r1_e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
509                  &                            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
510               zvau =-zfact1 * r1_e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
511                  &                            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
512               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
513               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
514            END DO 
515         END DO 
516         !                                             ! ===============
517      END DO                                           !   End of slab
518      !                                                ! ===============
519      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
520      !
521      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
522      !
523   END SUBROUTINE vor_ens
524
525
526   SUBROUTINE vor_een( kt, kvor, pua, pva )
527      !!----------------------------------------------------------------------
528      !!                ***  ROUTINE vor_een  ***
529      !!
530      !! ** Purpose :   Compute the now total vorticity trend and add it to
531      !!      the general trend of the momentum equation.
532      !!
533      !! ** Method  :   Trend evaluated using now fields (centered in time)
534      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
535      !!      both the horizontal kinetic energy and the potential enstrophy
536      !!      when horizontal divergence is zero (see the NEMO documentation)
537      !!      Add this trend to the general momentum trend (ua,va).
538      !!
539      !! ** Action : - Update (ua,va) with the now vorticity term trend
540      !!
541      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
542      !!----------------------------------------------------------------------
543      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
544      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
545      !                                                           ! =nrvm (relative vorticity or metric)
546      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
547      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
548      !!
549      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
550      INTEGER  ::   ierr                                          ! local integer
551      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
552      REAL(wp) ::   zmsk, ze3                                     ! local scalars
553      !                                                           !  3D workspace
554      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
555      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
556#if defined key_vvl
557      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
558#else
559      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
560#endif
561      !!----------------------------------------------------------------------
562      !
563      IF( nn_timing == 1 )  CALL timing_start('vor_een')
564      !
565      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
566      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
567#if defined key_vvl
568      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
569#endif
570      !
571      IF( kt == nit000 ) THEN
572         IF(lwp) WRITE(numout,*)
573         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
574         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
575#if ! defined key_vvl
576         IF( .NOT.ALLOCATED(ze3f) ) THEN
577            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
578            IF( lk_mpp    )   CALL mpp_sum ( ierr )
579            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
580         ENDIF
581         ze3f(:,:,:) = 0._wp
582#endif
583      ENDIF
584
585      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
586
587         IF( ln_dynvor_een_old ) THEN ! original formulation
588            DO jk = 1, jpk
589               DO jj = 1, jpjm1
590                  DO ji = 1, jpim1
591                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
592                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
593                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
594                  END DO
595               END DO
596            END DO
597         ELSE ! new formulation from NEMO 3.6
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         ENDIF
610
611         CALL lbc_lnk( ze3f, 'F', 1. )
612      ENDIF
613
614      zfac12 = 1._wp / 12._wp    ! Local constant initialization
615
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 * r1_e1e2f(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 * r1_e1e2f(ji,jj)   ) * ze3f(ji,jj,jk)
645               END DO
646            END DO
647            CALL lbc_lnk( zwz, 'F', 1. )
648         END SELECT
649         !
650         !                                   !==  horizontal fluxes  ==!
651         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
652         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
653
654         !                                   !==  compute and add the vorticity term trend  =!
655         jj = 2
656         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
657         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
658               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
659               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
660               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
661               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
662         END DO
663         DO jj = 3, jpj
664            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
665               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
666               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
667               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
668               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
669            END DO
670         END DO
671         DO jj = 2, jpjm1
672            DO ji = fs_2, fs_jpim1   ! vector opt.
673               zua = + zfac12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
674                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
675               zva = - zfac12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
676                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
677               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
678               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
679            END DO 
680         END DO 
681         !                                             ! ===============
682      END DO                                           !   End of slab
683      !                                                ! ===============
684      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
685      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
686#if defined key_vvl
687      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
688#endif
689      !
690      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
691      !
692   END SUBROUTINE vor_een
693
694
695   SUBROUTINE dyn_vor_init
696      !!---------------------------------------------------------------------
697      !!                  ***  ROUTINE dyn_vor_init  ***
698      !!
699      !! ** Purpose :   Control the consistency between cpp options for
700      !!              tracer advection schemes
701      !!----------------------------------------------------------------------
702      INTEGER ::   ioptio          ! local integer
703      INTEGER ::   ji, jj, jk      ! dummy loop indices
704      INTEGER ::   ios             ! Local integer output status for namelist read
705      !!
706      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
707      !!----------------------------------------------------------------------
708
709      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
710      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
711901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
712
713      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
714      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
715902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
716      IF(lwm) WRITE ( numond, namdyn_vor )
717
718      IF(lwp) THEN                    ! Namelist print
719         WRITE(numout,*)
720         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
721         WRITE(numout,*) '~~~~~~~~~~~~'
722         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
723         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
724         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
725         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
726         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
727         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
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( ln_dynvor_een_old )   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 .or. ln_dynvor_een_old )   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.