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_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynvor.F90 @ 12991

Last change on this file since 12991 was 12991, checked in by emanuelaclementi, 4 years ago

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