source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 8280

Last change on this file since 8280 was 8280, checked in by timgraham, 3 years ago

331: Merge of MEDUSA stable branch and HadGEM3 coupling branches into GO6 package branch.

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