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 NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/dynvor.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 3 years ago

Add Mixed Precision support by Oriol Tintó
  • Property svn:keywords set to Id
File size: 59.6 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   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
24   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
25   !!----------------------------------------------------------------------
26
27   !!----------------------------------------------------------------------
28   !!   dyn_vor       : Update the momentum trend with the vorticity trend
29   !!       vor_enT   : energy conserving scheme at T-pt  (ln_dynvor_enT=T)
30   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
31   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
32   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
33   !!       vor_eeT   : energy conserving at T-pt         (ln_dynvor_eeT=T)
34   !!   dyn_vor_init  : set and control of the different vorticity option
35   !!----------------------------------------------------------------------
36   USE oce            ! ocean dynamics and tracers
37   USE dom_oce        ! ocean space and time domain
38   USE dommsk         ! ocean mask
39   USE dynadv         ! momentum advection
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
43   USE sbc_oce,  ONLY : ln_stcor, ln_vortex_force   ! use Stoke-Coriolis force
44   !
45   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
46   USE prtctl         ! Print control
47   USE in_out_manager ! I/O manager
48   USE lib_mpp        ! MPP library
49   USE timing         ! Timing
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   dyn_vor        ! routine called by step.F90
55   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
56
57   !                                   !!* Namelist namdyn_vor: vorticity term
58   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
59   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
60   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
61   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
62   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
63   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
64   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
65   INTEGER, PUBLIC ::   nn_e3f_typ      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
66
67   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
68   !                                       ! associated indices:
69   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
70   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
71   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
72   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
73   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
74   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
75
76   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
77   !                               ! associated indices:
78   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
79   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
80   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
81   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
82   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
83
84   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
85   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
86   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
87   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
88   !
89   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   e3f_0vor   ! e3f used in EEN, ENE and ENS cases (key_qco only)
90
91   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
92   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
93   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
94
95   !! * Substitutions
96#  include "do_loop_substitute.h90"
97#  include "domzgr_substitute.h90"
98#  include "single_precision_substitute.h90"
99
100   !!----------------------------------------------------------------------
101   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
102   !! $Id$
103   !! Software governed by the CeCILL license (see ./LICENSE)
104   !!----------------------------------------------------------------------
105CONTAINS
106
107   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
108      !!----------------------------------------------------------------------
109      !!
110      !! ** Purpose :   compute the lateral ocean tracer physics.
111      !!
112      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
113      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
114      !!               and planetary vorticity trends) and send them to trd_dyn
115      !!               for futher diagnostics (l_trddyn=T)
116      !!----------------------------------------------------------------------
117      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
118      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
119      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
120      !
121      REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
122      !!----------------------------------------------------------------------
123      !
124      IF( ln_timing )   CALL timing_start('dyn_vor')
125      !
126      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
127         !
128         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
129         !
130         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend
131         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
132         SELECT CASE( nvor_scheme )
133         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
134         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
135         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
136         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
137         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
138         END SELECT
139         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
140         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
141         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
142         !
143         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
144            ztrdu(:,:,:) = puu(:,:,:,Krhs)
145            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
146            SELECT CASE( nvor_scheme )
147            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
148            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
149            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
150            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
151            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
152            END SELECT
153            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
154            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
155            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
156         ENDIF
157         !
158         DEALLOCATE( ztrdu, ztrdv )
159         !
160      ELSE              !==  total vorticity trend added to the general trend  ==!
161         !
162         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
163         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
164                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
165            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
166                             CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
167            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
168                             CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
169            ENDIF
170         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
171                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
172            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
173                             CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
174            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
175                             CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
176            ENDIF
177         CASE( np_ENE )                        !* energy conserving scheme
178                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
179            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
180                             CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
181            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
182                             CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
183            ENDIF
184         CASE( np_ENS )                        !* enstrophy conserving scheme
185                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
186
187            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
188                             CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
189            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
190                             CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend and vortex force
191            ENDIF
192         CASE( np_MIX )                        !* mixed ene-ens scheme
193                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
194                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
195            IF( ln_stcor )        CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )        ! add the Stokes-Coriolis trend
196            IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add vortex force
197         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
198                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
199            IF( ln_stcor .AND. .NOT. ln_vortex_force )  THEN
200                             CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
201            ELSE IF( ln_stcor .AND. ln_vortex_force )   THEN
202                             CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend and vortex force
203            ENDIF
204         END SELECT
205         !
206      ENDIF
207      !
208      !                       ! print sum trends (used for debugging)
209      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=CASTWP(puu(:,:,:,Krhs)), clinfo1=' vor  - Ua: ', mask1=umask,               &
210         &                                tab3d_2=CASTWP(pvv(:,:,:,Krhs)), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
211      !
212      IF( ln_timing )   CALL timing_stop('dyn_vor')
213      !
214   END SUBROUTINE dyn_vor
215
216
217   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
218      !!----------------------------------------------------------------------
219      !!                  ***  ROUTINE vor_enT  ***
220      !!
221      !! ** Purpose :   Compute the now total vorticity trend and add it to
222      !!      the general trend of the momentum equation.
223      !!
224      !! ** Method  :   Trend evaluated using now fields (centered in time)
225      !!       and t-point evaluation of vorticity (planetary and relative).
226      !!       conserves the horizontal kinetic energy.
227      !!         The general trend of momentum is increased due to the vorticity
228      !!       term which is given by:
229      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
230      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
231      !!       where rvor is the relative vorticity at f-point
232      !!
233      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
234      !!----------------------------------------------------------------------
235      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
236      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
237      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
238      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
239      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
240      !
241      INTEGER  ::   ji, jj, jk           ! dummy loop indices
242      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
243      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
244      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwz      ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
245      !!----------------------------------------------------------------------
246      !
247      IF( kt == nit000 ) THEN
248         IF(lwp) WRITE(numout,*)
249         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
250         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
251      ENDIF
252      !
253      !
254      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
255      !
256      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
257         ALLOCATE( zwz(jpi,jpj,jpk) )
258         DO jk = 1, jpkm1                                ! Horizontal slab
259            DO_2D( 1, 0, 1, 0 )
260               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
261                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
262            END_2D
263            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity
264               DO_2D( 1, 0, 1, 0 )
265                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
266               END_2D
267            ENDIF
268         END DO
269         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
270         !
271      END SELECT
272
273      !                                                ! ===============
274      DO jk = 1, jpkm1                                 ! Horizontal slab
275         !                                             ! ===============
276         !
277         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
278         !
279         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
280            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
281         CASE ( np_RVO )                           !* relative vorticity
282            DO_2D( 0, 1, 0, 1 )
283               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       &
284                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )   &
285                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
286            END_2D
287         CASE ( np_MET )                           !* metric term
288            DO_2D( 0, 1, 0, 1 )
289               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       &
290                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   &
291                  &       * e3t(ji,jj,jk,Kmm)
292            END_2D
293         CASE ( np_CRV )                           !* Coriolis + relative vorticity
294            DO_2D( 0, 1, 0, 1 )
295               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        &
296                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   &
297                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
298            END_2D
299         CASE ( np_CME )                           !* Coriolis + metric
300            DO_2D( 0, 1, 0, 1 )
301               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               &
302                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      &
303                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   &
304                    &     * e3t(ji,jj,jk,Kmm)
305            END_2D
306         CASE DEFAULT                                             ! error
307            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
308         END SELECT
309         !
310         !                                   !==  compute and add the vorticity term trend  =!
311         DO_2D( 0, 0, 0, 0 )
312            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
313               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
314               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
315               !
316            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
317               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &
318               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )
319         END_2D
320         !                                             ! ===============
321      END DO                                           !   End of slab
322      !                                                ! ===============
323      !
324      SELECT CASE( kvor )        ! deallocate zwz if necessary
325      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz )
326      END SELECT
327      !
328   END SUBROUTINE vor_enT
329
330
331   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
332      !!----------------------------------------------------------------------
333      !!                  ***  ROUTINE vor_ene  ***
334      !!
335      !! ** Purpose :   Compute the now total vorticity trend and add it to
336      !!      the general trend of the momentum equation.
337      !!
338      !! ** Method  :   Trend evaluated using now fields (centered in time)
339      !!       and the Sadourny (1975) flux form formulation : conserves the
340      !!       horizontal kinetic energy.
341      !!         The general trend of momentum is increased due to the vorticity
342      !!       term which is given by:
343      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
344      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
345      !!       where rvor is the relative vorticity
346      !!
347      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
348      !!
349      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
350      !!----------------------------------------------------------------------
351      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
352      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
353      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
354      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
355      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
356      !
357      INTEGER  ::   ji, jj, jk           ! dummy loop indices
358      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars
359      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
360      !!----------------------------------------------------------------------
361      !
362      IF( kt == nit000 ) THEN
363         IF(lwp) WRITE(numout,*)
364         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
365         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
366      ENDIF
367      !
368      !                                                ! ===============
369      DO jk = 1, jpkm1                                 ! Horizontal slab
370         !                                             ! ===============
371         !
372         SELECT CASE( kvor )                 !==  vorticity considered  ==!
373         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
374            zwz(:,:) = ff_f(:,:)
375         CASE ( np_RVO )                           !* relative vorticity
376            DO_2D( 1, 0, 1, 0 )
377               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
378                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
379            END_2D
380            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
381               DO_2D( 1, 0, 1, 0 )
382                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
383               END_2D
384            ENDIF
385         CASE ( np_MET )                           !* metric term
386            DO_2D( 1, 0, 1, 0 )
387               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
388                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
389            END_2D
390         CASE ( np_CRV )                           !* Coriolis + relative vorticity
391            DO_2D( 1, 0, 1, 0 )
392               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
393                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
394            END_2D
395            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
396               DO_2D( 1, 0, 1, 0 )
397                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
398               END_2D
399            ENDIF
400         CASE ( np_CME )                           !* Coriolis + metric
401            DO_2D( 1, 0, 1, 0 )
402               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
403                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
404            END_2D
405         CASE DEFAULT                                             ! error
406            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
407         END SELECT
408         !
409#if defined key_qco   ||   defined key_linssh
410         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
411            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
412         END_2D
413#else
414         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==!
415         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
416            DO_2D( 1, 0, 1, 0 )
417               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
418                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
419                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
420                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
421               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
422               ELSE                       ;   zwz(ji,jj) = 0._wp
423               ENDIF
424            END_2D
425         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
426            DO_2D( 1, 0, 1, 0 )
427               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
428                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
429                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
430                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
431               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
432                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
433               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
434               ELSE                       ;   zwz(ji,jj) = 0._wp
435               ENDIF
436            END_2D
437         END SELECT
438#endif
439         !                                   !==  horizontal fluxes  ==!
440         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
441         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
442         !
443         !                                   !==  compute and add the vorticity term trend  =!
444         DO_2D( 0, 0, 0, 0 )
445            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
446            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
447            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
448            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
449            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
450            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
451         END_2D
452         !                                             ! ===============
453      END DO                                           !   End of slab
454      !                                                ! ===============
455   END SUBROUTINE vor_ene
456
457
458   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
459      !!----------------------------------------------------------------------
460      !!                ***  ROUTINE vor_ens  ***
461      !!
462      !! ** Purpose :   Compute the now total vorticity trend and add it to
463      !!      the general trend of the momentum equation.
464      !!
465      !! ** Method  :   Trend evaluated using now fields (centered in time)
466      !!      and the Sadourny (1975) flux FORM formulation : conserves the
467      !!      potential enstrophy of a horizontally non-divergent flow. the
468      !!      trend of the vorticity term is given by:
469      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
470      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
471      !!      Add this trend to the general momentum trend:
472      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
473      !!
474      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
475      !!
476      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
477      !!----------------------------------------------------------------------
478      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
479      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
480      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
481      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
482      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
483      !
484      INTEGER  ::   ji, jj, jk   ! dummy loop indices
485      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars
486      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
487      !!----------------------------------------------------------------------
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_2D( 1, 0, 1, 0 )
503               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
504                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
505            END_2D
506            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
507               DO_2D( 1, 0, 1, 0 )
508                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
509               END_2D
510            ENDIF
511         CASE ( np_MET )                           !* metric term
512            DO_2D( 1, 0, 1, 0 )
513               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
514                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
515            END_2D
516         CASE ( np_CRV )                           !* Coriolis + relative vorticity
517            DO_2D( 1, 0, 1, 0 )
518               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
519                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
520            END_2D
521            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
522               DO_2D( 1, 0, 1, 0 )
523                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
524               END_2D
525            ENDIF
526         CASE ( np_CME )                           !* Coriolis + metric
527            DO_2D( 1, 0, 1, 0 )
528               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
529                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
530            END_2D
531         CASE DEFAULT                                             ! error
532            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
533         END SELECT
534         !
535         !
536#if defined key_qco   ||   defined key_linssh
537         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
538            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
539         END_2D
540#else
541         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==!
542         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
543            DO_2D( 1, 0, 1, 0 )
544               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
545                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
546                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
547                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
548               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
549               ELSE                       ;   zwz(ji,jj) = 0._wp
550               ENDIF
551            END_2D
552         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
553            DO_2D( 1, 0, 1, 0 )
554               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
555                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
556                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
557                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
558               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
559                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
560               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
561               ELSE                       ;   zwz(ji,jj) = 0._wp
562               ENDIF
563            END_2D
564         END SELECT
565#endif
566         !                                   !==  horizontal fluxes  ==!
567         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
568         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
569         !
570         !                                   !==  compute and add the vorticity term trend  =!
571         DO_2D( 0, 0, 0, 0 )
572            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
573               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
574            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
575               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
576            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
577            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
578         END_2D
579         !                                             ! ===============
580      END DO                                           !   End of slab
581      !                                                ! ===============
582   END SUBROUTINE vor_ens
583
584
585   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
586      !!----------------------------------------------------------------------
587      !!                ***  ROUTINE vor_een  ***
588      !!
589      !! ** Purpose :   Compute the now total vorticity trend and add it to
590      !!      the general trend of the momentum equation.
591      !!
592      !! ** Method  :   Trend evaluated using now fields (centered in time)
593      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
594      !!      both the horizontal kinetic energy and the potential enstrophy
595      !!      when horizontal divergence is zero (see the NEMO documentation)
596      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
597      !!
598      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
599      !!
600      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
601      !!----------------------------------------------------------------------
602      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
603      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
604      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
605      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
606      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
607      !
608      INTEGER  ::   ji, jj, jk   ! dummy loop indices
609      INTEGER  ::   ierr         ! local integer
610      REAL(wp) ::   zua, zva     ! local scalars
611      REAL(wp) ::   zmsk, ze3f   ! local scalars
612      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy , z1_e3f
613      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
614      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
615      !!----------------------------------------------------------------------
616      !
617      IF( kt == nit000 ) THEN
618         IF(lwp) WRITE(numout,*)
619         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
620         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
621      ENDIF
622      !
623      !                                                ! ===============
624      DO jk = 1, jpkm1                                 ! Horizontal slab
625         !                                             ! ===============
626         !
627#if defined key_qco   ||   defined key_linssh
628         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco)
629            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
630         END_2D
631#else
632         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point
633         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
634            DO_2D( 1, 0, 1, 0 )
635               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
636                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
637                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
638                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
639               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
640               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
641               ENDIF
642            END_2D
643         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
644            DO_2D( 1, 0, 1, 0 )
645               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
646                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
647                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
648                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
649               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
650                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
651               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
652               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
653               ENDIF
654            END_2D
655         END SELECT
656#endif
657         !
658         SELECT CASE( kvor )                 !==  vorticity considered  ==!
659         !
660         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
661            DO_2D( 1, 0, 1, 0 )
662               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
663            END_2D
664         CASE ( np_RVO )                           !* relative vorticity
665            DO_2D( 1, 0, 1, 0 )
666               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
667                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
668            END_2D
669            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
670               DO_2D( 1, 0, 1, 0 )
671                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
672               END_2D
673            ENDIF
674         CASE ( np_MET )                           !* metric term
675            DO_2D( 1, 0, 1, 0 )
676               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
677                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
678            END_2D
679         CASE ( np_CRV )                           !* Coriolis + relative vorticity
680            DO_2D( 1, 0, 1, 0 )
681               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
682                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
683                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
684            END_2D
685            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
686               DO_2D( 1, 0, 1, 0 )
687                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
688               END_2D
689            ENDIF
690         CASE ( np_CME )                           !* Coriolis + metric
691            DO_2D( 1, 0, 1, 0 )
692               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
693                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
694            END_2D
695         CASE DEFAULT                                             ! error
696            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
697         END SELECT
698         !                                             ! ===============
699      END DO                                           !   End of slab
700      !                                                ! ===============
701      !
702      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
703      !
704      !                                                ! ===============
705      DO jk = 1, jpkm1                                 ! Horizontal slab
706         !                                             ! ===============
707         !
708         !                                   !==  horizontal fluxes  ==!
709         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
710         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
711         !
712         !                                   !==  compute and add the vorticity term trend  =!
713         DO_2D( 0, 1, 0, 1 )
714            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
715            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
716            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
717            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
718         END_2D
719         !
720         DO_2D( 0, 0, 0, 0 )
721            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
722               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
723            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
724               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
725            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
726            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
727         END_2D
728         !                                             ! ===============
729      END DO                                           !   End of slab
730      !                                                ! ===============
731   END SUBROUTINE vor_een
732
733
734   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
735      !!----------------------------------------------------------------------
736      !!                ***  ROUTINE vor_eeT  ***
737      !!
738      !! ** Purpose :   Compute the now total vorticity trend and add it to
739      !!      the general trend of the momentum equation.
740      !!
741      !! ** Method  :   Trend evaluated using now fields (centered in time)
742      !!      and the Arakawa and Lamb (1980) vector form formulation using
743      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
744      !!      The change consists in
745      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
746      !!
747      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
748      !!
749      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
750      !!----------------------------------------------------------------------
751      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
752      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
753      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
754      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
755      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
756      !
757      INTEGER  ::   ji, jj, jk     ! dummy loop indices
758      INTEGER  ::   ierr           ! local integer
759      REAL(wp) ::   zua, zva       ! local scalars
760      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
761      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy
762      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
763      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
764      !!----------------------------------------------------------------------
765      !
766      IF( kt == nit000 ) THEN
767         IF(lwp) WRITE(numout,*)
768         IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
769         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
770      ENDIF
771      !
772      !                                                ! ===============
773      DO jk = 1, jpkm1                                 ! Horizontal slab
774         !                                             ! ===============
775         !
776         !
777         SELECT CASE( kvor )                 !==  vorticity considered  ==!
778         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
779            DO_2D( 1, 0, 1, 0 )
780               zwz(ji,jj,jk) = ff_f(ji,jj)
781            END_2D
782         CASE ( np_RVO )                           !* relative vorticity
783            DO_2D( 1, 0, 1, 0 )
784               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
785                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
786                  &          * r1_e1e2f(ji,jj)
787            END_2D
788            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
789               DO_2D( 1, 0, 1, 0 )
790                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
791               END_2D
792            ENDIF
793         CASE ( np_MET )                           !* metric term
794            DO_2D( 1, 0, 1, 0 )
795               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
796                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
797            END_2D
798         CASE ( np_CRV )                           !* Coriolis + relative vorticity
799            DO_2D( 1, 0, 1, 0 )
800               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
801                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
802                  &                         * r1_e1e2f(ji,jj)    )
803            END_2D
804            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
805               DO_2D( 1, 0, 1, 0 )
806                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
807               END_2D
808            ENDIF
809         CASE ( np_CME )                           !* Coriolis + metric
810            DO_2D( 1, 0, 1, 0 )
811               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
812                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
813            END_2D
814         CASE DEFAULT                                             ! error
815            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
816         END SELECT
817         !
818         !                                             ! ===============
819      END DO                                           !   End of slab
820      !                                                ! ===============
821      !
822      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
823      !
824      !                                                ! ===============
825      DO jk = 1, jpkm1                                 ! Horizontal slab
826         !                                             ! ===============
827         !
828         !                                   !==  horizontal fluxes  ==!
829         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
830         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
831         !
832         !                                   !==  compute and add the vorticity term trend  =!
833         DO_2D( 0, 1, 0, 1 )
834            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
835            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
836            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
837            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
838            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
839         END_2D
840         !
841         DO_2D( 0, 0, 0, 0 )
842            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
843               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
844            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
845               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
846            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
847            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
848         END_2D
849         !                                             ! ===============
850      END DO                                           !   End of slab
851      !                                                ! ===============
852   END SUBROUTINE vor_eeT
853
854
855   SUBROUTINE dyn_vor_init
856      !!---------------------------------------------------------------------
857      !!                  ***  ROUTINE dyn_vor_init  ***
858      !!
859      !! ** Purpose :   Control the consistency between cpp options for
860      !!              tracer advection schemes
861      !!----------------------------------------------------------------------
862      INTEGER ::   ji, jj, jk    ! dummy loop indices
863      INTEGER ::   ioptio, ios   ! local integer
864      REAL(wp) ::   zmsk    ! local scalars
865      !!
866      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
867         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk
868      !!----------------------------------------------------------------------
869      !
870      IF(lwp) THEN
871         WRITE(numout,*)
872         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
873         WRITE(numout,*) '~~~~~~~~~~~~'
874      ENDIF
875      !
876      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
877901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
878      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
879902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
880      IF(lwm) WRITE ( numond, namdyn_vor )
881      !
882      IF(lwp) THEN                    ! Namelist print
883         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
884         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
885         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
886         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
887         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
888         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
889         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ
890         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
891         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
892      ENDIF
893
894!!gm  this should be removed when choosing a unique strategy for fmask at the coast
895      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
896      ! at angles with three ocean points and one land point
897      IF(lwp) WRITE(numout,*)
898      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
899      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
900         DO_3D( 1, 0, 1, 0, 1, jpk )
901            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
902               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
903         END_3D
904         !
905         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
906         !
907      ENDIF
908!!gm end
909
910      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
911      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
912      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
913      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
914      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
915      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
916      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
917      !
918      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
919      !
920      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
921      ncor = np_COR                       ! planetary vorticity
922      SELECT CASE( n_dynadv )
923      CASE( np_LIN_dyn )
924         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
925         nrvm = np_COR        ! planetary vorticity
926         ntot = np_COR        !     -         -
927      CASE( np_VEC_c2  )
928         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'
929         nrvm = np_RVO        ! relative vorticity
930         ntot = np_CRV        ! relative + planetary vorticity
931      CASE( np_FLX_c2 , np_FLX_ubs  )
932         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
933         nrvm = np_MET        ! metric term
934         ntot = np_CME        ! Coriolis + metric term
935         !
936         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
937         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
938            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
939            DO_2D( 0, 0, 0, 0 )
940               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
941               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
942            END_2D
943            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
944            !
945         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
946            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
947            DO_2D( 1, 0, 1, 0 )
948               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
949               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
950            END_2D
951            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
952         END SELECT
953         !
954      END SELECT
955#if defined key_qco   ||   defined key_linssh
956      SELECT CASE( nvor_scheme )    ! qco or linssh cases : pre-computed a specific e3f_0 for some vorticity schemes
957      CASE( np_ENS , np_ENE , np_EEN , np_MIX )
958         !
959         ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
960         !
961         SELECT CASE( nn_e3f_typ )
962         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4)
963            DO_3D( 0, 0, 0, 0, 1, jpk )
964               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
965                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
966                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
967                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp
968            END_3D
969         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask)
970            DO_3D( 0, 0, 0, 0, 1, jpk )
971               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   &
972                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  )
973               !
974               IF( zmsk /= 0._wp ) THEN
975                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
976                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
977                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
978                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk
979               ENDIF
980            END_3D
981         END SELECT
982         !
983         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
984         !                                 ! insure e3f_0vor /= 0
985         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:)
986         !
987      END SELECT
988      !
989#endif
990      IF(lwp) THEN                   ! Print the choice
991         WRITE(numout,*)
992         SELECT CASE( nvor_scheme )
993         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
994         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
995         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
996                              IF( ln_dynadv_vec )   CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form')
997         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
998         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
999         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
1000         END SELECT
1001      ENDIF
1002      !
1003   END SUBROUTINE dyn_vor_init
1004
1005   !!==============================================================================
1006END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.