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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7698

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

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