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 @ 8300

Last change on this file since 8300 was 8300, checked in by glong, 7 years ago

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