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/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynvor.F90 @ 10371

Last change on this file since 10371 was 10371, checked in by davestorkey, 5 years ago

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