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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

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