New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynvor.F90 in branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7649

Last change on this file since 7649 was 7649, checked in by glong, 8 years ago

dyn_vrt_dia subroutine added and calls added for most processes - spg, zdf, and bfr not yet done
File size: 42.2 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec) vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
17   !!            3.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 divcur         ! For dyn_vrt_dia
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  :: id_dia_vrt_vor_int = 1                          ! TODO remove once flags set properly
564      INTEGER  :: id_dia_vrt_vor_mn  = 1                          ! TODO remove once flags set properly
565      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
566      INTEGER  ::   ierr                                          ! local integer
567      REAL(wp) ::   zfac12                                        ! local scalars
568      REAL(wp) ::   zmsk, ze3                                     ! local scalars
569      !                                                           !  3D workspace
570      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
571      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
572      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: zua, zva
573#if defined key_vvl
574      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
575#else
576      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
577#endif
578      !!----------------------------------------------------------------------
579      !
580      IF( nn_timing == 1 )  CALL timing_start('vor_een')
581      !
582      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
583      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
584      CALL wrk_alloc( jpi, jpj, jpk, zua, zva               ) 
585#if defined key_vvl
586      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
587#endif
588      !
589      IF( kt == nit000 ) THEN
590         IF(lwp) WRITE(numout,*)
591         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
592         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
593#if ! defined key_vvl
594         IF( .NOT.ALLOCATED(ze3f) ) THEN
595            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
596            IF( lk_mpp    )   CALL mpp_sum ( ierr )
597            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
598         ENDIF
599         ze3f(:,:,:) = 0._wp
600#endif
601      ENDIF
602
603      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
604
605         IF( ln_dynvor_een_old ) THEN ! original formulation
606            DO jk = 1, jpk
607               DO jj = 1, jpjm1
608                  DO ji = 1, jpim1
609                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
610                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
611                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
612                  END DO
613               END DO
614            END DO
615         ELSE ! new formulation from NEMO 3.6
616            DO jk = 1, jpk
617               DO jj = 1, jpjm1
618                  DO ji = 1, jpim1
619                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
620                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
621                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
622                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
623                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
624                  END DO
625               END DO
626            END DO
627         ENDIF
628
629         CALL lbc_lnk( ze3f, 'F', 1. )
630      ENDIF
631
632      zfac12 = 1._wp / 12._wp    ! Local constant initialization
633
634     
635!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
636      !                                                ! ===============
637      DO jk = 1, jpkm1                                 ! Horizontal slab
638         !                                             ! ===============
639         
640         ! Potential vorticity and horizontal fluxes
641         ! -----------------------------------------
642         SELECT CASE( kvor )      ! vorticity considered
643         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
644            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
645         CASE ( 2 )                                                ! relative  vorticity
646            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
647         CASE ( 3 )                                                ! metric term
648            DO jj = 1, jpjm1
649               DO ji = 1, fs_jpim1   ! vector opt.
650                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
651                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
652                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
653               END DO
654            END DO
655            CALL lbc_lnk( zwz, 'F', 1. )
656        CASE ( 4 )                                                ! total (relative + planetary vorticity)
657            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
658         CASE ( 5 )                                                ! total (coriolis + metric)
659            DO jj = 1, jpjm1
660               DO ji = 1, fs_jpim1   ! vector opt.
661                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
662                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
663                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
664                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
665                       &       ) * ze3f(ji,jj,jk)
666               END DO
667            END DO
668            CALL lbc_lnk( zwz, 'F', 1. )
669         END SELECT
670
671         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
672         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
673
674         ! Compute and add the vorticity term trend
675         ! ----------------------------------------
676         jj = 2
677         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
678         DO ji = 2, jpi   
679               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
680               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
681               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
682               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
683         END DO
684         DO jj = 3, jpj
685            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
686               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
687               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
688               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
689               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
690            END DO
691         END DO
692         DO jj = 2, jpjm1
693            DO ji = fs_2, fs_jpim1   ! vector opt.
694               zua(ji,jj,jk) = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  )     &
695                  &                                     + ztnw(ji+1,jj) * zwy(ji+1,jj  )     &
696                  &                                     + ztse(ji,jj  ) * zwy(ji  ,jj-1)     &
697                  &                                     + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
698               zva(ji,jj,jk) = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1)     &
699                  &                                     + ztse(ji,jj+1) * zwx(ji  ,jj+1)     &
700                  &                                     + ztnw(ji,jj  ) * zwx(ji-1,jj  )     &
701                  &                                     + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
702               pua(ji,jj,jk) = pua(ji,jj,jk) + zua(ji,jj,jk)
703               pva(ji,jj,jk) = pva(ji,jj,jk) + zva(ji,jj,jk)
704            END DO 
705         END DO 
706         !                                             ! ===============
707      END DO                                           !   End of slab
708      !                                                ! ===============
709      IF ( ( id_dia_vrt_vor_int == 1 ) .or. ( id_dia_vrt_vor_mn == 1 ) ) THEN
710          CALL dyn_vrt_dia(zua, zva, id_dia_vrt_vor_int, id_dia_vrt_vor_mn)
711      END IF
712      !
713      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
714      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
715      CALL wrk_dealloc( jpi, jpj, jpk, zua, zva               ) 
716#if defined key_vvl
717      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
718#endif
719      !
720      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
721      !
722   END SUBROUTINE vor_een
723
724
725   SUBROUTINE dyn_vor_init
726      !!---------------------------------------------------------------------
727      !!                  ***  ROUTINE dyn_vor_init  ***
728      !!
729      !! ** Purpose :   Control the consistency between cpp options for
730      !!              tracer advection schemes
731      !!----------------------------------------------------------------------
732      INTEGER ::   ioptio          ! local integer
733      INTEGER ::   ji, jj, jk      ! dummy loop indices
734      INTEGER ::   ios             ! Local integer output status for namelist read
735      !!
736      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
737      !!----------------------------------------------------------------------
738
739      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
740      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
741901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
742
743      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
744      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
745902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
746      IF(lwm) WRITE ( numond, namdyn_vor )
747
748      IF(lwp) THEN                    ! Namelist print
749         WRITE(numout,*)
750         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
751         WRITE(numout,*) '~~~~~~~~~~~~'
752         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
753         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
754         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
755         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
756         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
757         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
758      ENDIF
759
760      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
761      ! at angles with three ocean points and one land point
762      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
763         DO jk = 1, jpk
764            DO jj = 2, jpjm1
765               DO ji = 2, jpim1
766                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
767                      fmask(ji,jj,jk) = 1._wp
768               END DO
769            END DO
770         END DO
771          !
772          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
773          !
774      ENDIF
775
776      ioptio = 0                     ! Control of vorticity scheme options
777      IF( ln_dynvor_ene )   ioptio = ioptio + 1
778      IF( ln_dynvor_ens )   ioptio = ioptio + 1
779      IF( ln_dynvor_mix )   ioptio = ioptio + 1
780      IF( ln_dynvor_een )   ioptio = ioptio + 1
781      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
782      IF( lk_esopa      )   ioptio =          1
783
784      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
785
786      !                              ! Set nvor (type of scheme for vorticity)
787      IF( ln_dynvor_ene )   nvor =  0
788      IF( ln_dynvor_ens )   nvor =  1
789      IF( ln_dynvor_mix )   nvor =  2
790      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
791      IF( lk_esopa      )   nvor = -1
792     
793      !                              ! Set ncor, nrvm, ntot (type of vorticity)
794      IF(lwp) WRITE(numout,*)
795      ncor = 1
796      IF( ln_dynadv_vec ) THEN     
797         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
798         nrvm = 2
799         ntot = 4
800      ELSE                       
801         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
802         nrvm = 3
803         ntot = 5
804      ENDIF
805     
806      IF(lwp) THEN                   ! Print the choice
807         WRITE(numout,*)
808         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
809         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
810         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
811         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
812         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
813      ENDIF
814      !
815   END SUBROUTINE dyn_vor_init
816
817   !!==============================================================================
818END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.