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

source: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7606

Last change on this file since 7606 was 7606, checked in by jcastill, 7 years ago

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