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
RevLine 
[3]1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
[2715]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
[4990]17   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
[5836]18   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory
[7646]19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
[503]20   !!----------------------------------------------------------------------
[3]21
22   !!----------------------------------------------------------------------
[2528]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
[3]28   !!----------------------------------------------------------------------
[503]29   USE oce            ! ocean dynamics and tracers
30   USE dom_oce        ! ocean space and time domain
[3294]31   USE dommsk         ! ocean mask
[643]32   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
[4990]33   USE trd_oce        ! trends: ocean variables
34   USE trddyn         ! trend manager: dynamics
[7646]35   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
36   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
[5836]37   !
[503]38   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
39   USE prtctl         ! Print control
40   USE in_out_manager ! I/O manager
[3294]41   USE lib_mpp        ! MPP library
42   USE wrk_nemo       ! Memory Allocation
43   USE timing         ! Timing
[3]44
[3294]45
[3]46   IMPLICIT NONE
47   PRIVATE
48
[2528]49   PUBLIC   dyn_vor        ! routine called by step.F90
[5836]50   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
[3]51
[4147]52   !                                   !!* Namelist namdyn_vor: vorticity term
[5836]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)
[3]59
[5836]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
[455]66
[5836]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   
[3]79   !! * Substitutions
80#  include "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
[7646]82   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
[1152]83   !! $Id$
[2715]84   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]85   !!----------------------------------------------------------------------
86CONTAINS
87
[455]88   SUBROUTINE dyn_vor( kt )
[3]89      !!----------------------------------------------------------------------
90      !!
[455]91      !! ** Purpose :   compute the lateral ocean tracer physics.
92      !!
93      !! ** Action : - Update (ua,va) with the now vorticity term trend
[503]94      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]95      !!               and planetary vorticity trends) and send them to trd_dyn
96      !!               for futher diagnostics (l_trddyn=T)
[503]97      !!----------------------------------------------------------------------
[3294]98      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
[7698]99      INTEGER ::   jk, jj, ji
[2715]100      !
[3294]101      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
[455]102      !!----------------------------------------------------------------------
[2715]103      !
[3294]104      IF( nn_timing == 1 )  CALL timing_start('dyn_vor')
105      !
106      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
107      !
[5836]108      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==!
[643]109      !
[5836]110      CASE ( np_ENE )                                 !* energy conserving scheme
111         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
[7698]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
[7646]121            CALL vor_ene( kt, nrvm, un , vn , ua, va )                    ! relative vorticity or metric trend
[7698]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
[4990]131            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[7698]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
[7646]141            CALL vor_ene( kt, ncor, un , vn , ua, va )                    ! planetary vorticity trend
[7698]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
[4990]151            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[7646]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
[455]155         ENDIF
[643]156         !
[5836]157      CASE ( np_ENS )                                 !* enstrophy conserving scheme
158         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two   
[7698]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
[7646]168            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
[7698]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
[4990]178            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[7698]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
[7646]188            CALL vor_ens( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
[7698]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
[4990]198            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[7646]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
[455]202         ENDIF
[643]203         !
[5836]204      CASE ( np_MIX )                                 !* mixed ene-ens scheme
205         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
[7698]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
[7646]215            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend (ens)
[7698]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
[4990]225            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[7698]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
[7646]235            CALL vor_ene( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend (ene)
[7698]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
[4990]245            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[7646]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
[5836]250        ENDIF
[643]251         !
[5836]252      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme
253         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
[7698]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
[7646]263            CALL vor_een( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend
[7698]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
[4990]273            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[7698]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
[7646]283            CALL vor_een( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend
[7698]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
[4990]293            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[7646]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
[455]297         ENDIF
[643]298         !
[455]299      END SELECT
[2715]300      !
[455]301      !                       ! print sum trends (used for debugging)
[2715]302      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
[455]303         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]304      !
[3294]305      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
306      !
307      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor')
308      !
[455]309   END SUBROUTINE dyn_vor
310
311
[7646]312   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
[455]313      !!----------------------------------------------------------------------
314      !!                  ***  ROUTINE vor_ene  ***
315      !!
[3]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)
[5836]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
[3]327      !!
328      !! ** Action : - Update (ua,va) with the now vorticity term trend
329      !!
[503]330      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]331      !!----------------------------------------------------------------------
[7646]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
[2715]337      !
[5836]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
[3]341      !!----------------------------------------------------------------------
[3294]342      !
343      IF( nn_timing == 1 )  CALL timing_start('vor_ene')
344      !
[7646]345      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
[3294]346      !
[52]347      IF( kt == nit000 ) THEN
348         IF(lwp) WRITE(numout,*)
[455]349         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
350         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]351      ENDIF
[5836]352      !
[3]353      !                                                ! ===============
354      DO jk = 1, jpkm1                                 ! Horizontal slab
355         !                                             ! ===============
[1438]356         !
[5836]357         SELECT CASE( kvor )                 !==  vorticity considered  ==!
358         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[7698]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
[5836]365         CASE ( np_RVO )                           !* relative vorticity
[7698]366!$OMP PARALLEL DO schedule(static) private(jj,ji)
[643]367            DO jj = 1, jpjm1
368               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)
[5836]371               END DO
372            END DO
373         CASE ( np_MET )                           !* metric term
[7698]374!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]375            DO jj = 1, jpjm1
376               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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) )   )   &
[5836]379                       &     * 0.5 * r1_e1e2f(ji,jj)
[643]380               END DO
381            END DO
[5836]382         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[7698]383!$OMP PARALLEL DO schedule(static) private(jj,ji)
[643]384            DO jj = 1, jpjm1
385               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)  ) &
[5836]388                     &                   * r1_e1e2f(ji,jj)
[643]389               END DO
390            END DO
[5836]391         CASE ( np_CME )                           !* Coriolis + metric
[7698]392!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]393            DO jj = 1, jpjm1
394               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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) )   )   &
[5836]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'  )
[455]403         END SELECT
[5836]404         !
405         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
[7698]406!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]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
[455]413
414         IF( ln_sco ) THEN
[7698]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
[3]423         ELSE
[7698]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
[3]431         ENDIF
[5836]432         !                                   !==  compute and add the vorticity term trend  =!
[7698]433!$OMP PARALLEL DO schedule(static) private(jj, ji, zy1, zy2, zx1, zx2)
[3]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)
[5836]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 ) 
[3]442            END DO 
443         END DO 
444         !                                             ! ===============
445      END DO                                           !   End of slab
446      !                                                ! ===============
[3294]447      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
[2715]448      !
[3294]449      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
450      !
[455]451   END SUBROUTINE vor_ene
[216]452
453
[7646]454   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
[3]455      !!----------------------------------------------------------------------
[455]456      !!                ***  ROUTINE vor_ens  ***
[3]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:
[5836]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) ]
[3]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      !!
[503]472      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]473      !!----------------------------------------------------------------------
[7646]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
[2715]479      !
[5836]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
[3]483      !!----------------------------------------------------------------------
[3294]484      !
485      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
486      !
[7646]487      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
[3294]488      !
[52]489      IF( kt == nit000 ) THEN
490         IF(lwp) WRITE(numout,*)
[455]491         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
492         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]493      ENDIF
[3]494      !                                                ! ===============
495      DO jk = 1, jpkm1                                 ! Horizontal slab
496         !                                             ! ===============
[1438]497         !
[5836]498         SELECT CASE( kvor )                 !==  vorticity considered  ==!
499         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[7646]500            zwz(:,:) = ff_f(:,:) 
[5836]501         CASE ( np_RVO )                           !* relative vorticity
[643]502            DO jj = 1, jpjm1
503               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)
[5836]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.
[7646]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) )   )   &
[5836]513                       &     * 0.5 * r1_e1e2f(ji,jj)
[643]514               END DO
515            END DO
[5836]516         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[643]517            DO jj = 1, jpjm1
518               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)  ) &
[5836]521                     &                   * r1_e1e2f(ji,jj)
[643]522               END DO
523            END DO
[5836]524         CASE ( np_CME )                           !* Coriolis + metric
525            DO jj = 1, jpjm1
526               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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) )   )   &
[5836]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'  )
[455]535         END SELECT
[1438]536         !
[5836]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)
[3]541               END DO
542            END DO
[5836]543         ENDIF
544         !
545         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
[6140]546            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
[7646]547            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
548            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
[3]549         ELSE
[7646]550            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
551            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
[3]552         ENDIF
[5836]553         !                                   !==  compute and add the vorticity term trend  =!
[3]554         DO jj = 2, jpjm1
555            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]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)  )
[455]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) )
[3]562            END DO 
563         END DO 
564         !                                             ! ===============
565      END DO                                           !   End of slab
566      !                                                ! ===============
[3294]567      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
[2715]568      !
[3294]569      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
570      !
[455]571   END SUBROUTINE vor_ens
[216]572
573
[7646]574   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
[108]575      !!----------------------------------------------------------------------
[455]576      !!                ***  ROUTINE vor_een  ***
[108]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)
[1438]582      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]583      !!      both the horizontal kinetic energy and the potential enstrophy
[1438]584      !!      when horizontal divergence is zero (see the NEMO documentation)
585      !!      Add this trend to the general momentum trend (ua,va).
[108]586      !!
587      !! ** Action : - Update (ua,va) with the now vorticity term trend
588      !!
[503]589      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
590      !!----------------------------------------------------------------------
[7646]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
[5836]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
[108]604      !!----------------------------------------------------------------------
[3294]605      !
606      IF( nn_timing == 1 )  CALL timing_start('vor_een')
607      !
[5836]608      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
609      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
[3294]610      !
[108]611      IF( kt == nit000 ) THEN
612         IF(lwp) WRITE(numout,*)
[455]613         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
614         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]615      ENDIF
[5836]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)
[7698]623!$OMP PARALLEL DO schedule(static) private(jj,ji,ze3)
[5836]624            DO jj = 1, jpjm1
625               DO ji = 1, fs_jpim1   ! vector opt.
[6140]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
[5836]630                  ENDIF
[108]631               END DO
632            END DO
[5836]633         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
[7698]634!$OMP PARALLEL DO schedule(static) private(jj,ji,ze3,zmsk)
[5836]635            DO jj = 1, jpjm1
636               DO ji = 1, fs_jpim1   ! vector opt.
[6140]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)  )
[5836]641                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3
[6140]642                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp
[5836]643                  ENDIF
[5029]644               END DO
645            END DO
[5836]646         END SELECT
647         !
648         SELECT CASE( kvor )                 !==  vorticity considered  ==!
649         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[7698]650!$OMP PARALLEL DO schedule(static) private(jj,ji)
[643]651            DO jj = 1, jpjm1
652               DO ji = 1, fs_jpim1   ! vector opt.
[7646]653                  zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
[5836]654               END DO
655            END DO
656         CASE ( np_RVO )                           !* relative vorticity
[7698]657!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]658            DO jj = 1, jpjm1
659               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)  ) &
[5836]662                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
663               END DO
664            END DO
665         CASE ( np_MET )                           !* metric term
[7698]666!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]667            DO jj = 1, jpjm1
668               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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) )   )   &
[5836]671                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
[643]672               END DO
673            END DO
[5836]674         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[7698]675!$OMP PARALLEL DO schedule(static) private(jj,ji)
[643]676            DO jj = 1, jpjm1
677               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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)  ) &
[5836]680                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
[643]681               END DO
682            END DO
[5836]683         CASE ( np_CME )                           !* Coriolis + metric
[7698]684!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]685            DO jj = 1, jpjm1
686               DO ji = 1, fs_jpim1   ! vector opt.
[7646]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) )   )   &
[5836]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'  )
[455]695         END SELECT
[5836]696         !
697         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
[7698]698!$OMP PARALLEL DO schedule(static) private(jj,ji)
[5836]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         !
[5907]706         CALL lbc_lnk( zwz, 'F', 1. )
707         !
[5836]708         !                                   !==  horizontal fluxes  ==!
[7698]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
[108]716
[5836]717         !                                   !==  compute and add the vorticity term trend  =!
[1438]718         jj = 2
719         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
[7698]720
[5836]721         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
[108]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
[7698]727!$OMP PARALLEL
728!$OMP DO schedule(static) private(jj,ji)
[108]729         DO jj = 3, jpj
[1694]730            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
[108]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
[7698]737!$OMP DO schedule(static) private(jj,ji,zua,zva)
[108]738         DO jj = 2, jpjm1
739            DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]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  ) )
[455]744               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
745               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]746            END DO 
747         END DO 
[7698]748!$OMP END PARALLEL
[108]749         !                                             ! ===============
750      END DO                                           !   End of slab
751      !                                                ! ===============
[2715]752      !
[5836]753      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
754      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
755      !
[3294]756      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
757      !
[455]758   END SUBROUTINE vor_een
[216]759
760
[2528]761   SUBROUTINE dyn_vor_init
[3]762      !!---------------------------------------------------------------------
[2528]763      !!                  ***  ROUTINE dyn_vor_init  ***
[3]764      !!
765      !! ** Purpose :   Control the consistency between cpp options for
[1438]766      !!              tracer advection schemes
[3]767      !!----------------------------------------------------------------------
[2715]768      INTEGER ::   ioptio          ! local integer
[3294]769      INTEGER ::   ji, jj, jk      ! dummy loop indices
[4147]770      INTEGER ::   ios             ! Local integer output status for namelist read
[2715]771      !!
[5836]772      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
[3]773      !!----------------------------------------------------------------------
774
[4147]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 )
[3]778
[4147]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 )
[4624]782      IF(lwm) WRITE ( numond, namdyn_vor )
[4147]783
[503]784      IF(lwp) THEN                    ! Namelist print
[3]785         WRITE(numout,*)
[2528]786         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
787         WRITE(numout,*) '~~~~~~~~~~~~'
[7646]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
[52]795      ENDIF
796
[5836]797!!gm  this should be removed when choosing a unique strategy for fmask at the coast
[3294]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
[5836]800      IF(lwp) WRITE(numout,*)
[7646]801      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
[3294]802      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
[7698]803!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
[3294]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
[5836]816!!gm end
[3294]817
[5836]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      !
[6140]824      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[5836]825      !                     
826      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
827      ncor = np_COR
[643]828      IF( ln_dynadv_vec ) THEN     
[7646]829         IF(lwp) WRITE(numout,*) '      ===>>   Vector form advection : vorticity = Coriolis + relative vorticity'
[5836]830         nrvm = np_RVO        ! relative vorticity
831         ntot = np_CRV        ! relative + planetary vorticity
[643]832      ELSE                       
[7646]833         IF(lwp) WRITE(numout,*) '      ===>>   Flux form advection   : vorticity = Coriolis + metric term'
[5836]834         nrvm = np_MET        ! metric term
835         ntot = np_CME        ! Coriolis + metric term
[643]836      ENDIF
837     
[503]838      IF(lwp) THEN                   ! Print the choice
839         WRITE(numout,*)
[7646]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'
[3]844      ENDIF
[503]845      !
[2528]846   END SUBROUTINE dyn_vor_init
[3]847
[503]848   !!==============================================================================
[3]849END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.