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

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

source: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 9176

Last change on this file since 9176 was 9176, checked in by andmirek, 6 years ago

#2001: OMP directives

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