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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7508

Last change on this file since 7508 was 7508, checked in by mocavero, 7 years ago

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