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_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90 @ 13644

Last change on this file since 13644 was 13644, checked in by techene, 4 years ago

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