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_r13296_HPC-07_mocavero_mpi3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DYN/dynvor.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496
  • Property svn:keywords set to Id
File size: 53.8 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity
18   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!----------------------------------------------------------------------
24
25   !!----------------------------------------------------------------------
26   !!   dyn_vor       : Update the momentum trend with the vorticity trend
27   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
28   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
29   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
30   !!   dyn_vor_init  : set and control of the different vorticity option
31   !!----------------------------------------------------------------------
32   USE oce            ! ocean dynamics and tracers
33   USE dom_oce        ! ocean space and time domain
34   USE dommsk         ! ocean mask
35   USE dynadv         ! momentum advection
36   USE trd_oce        ! trends: ocean variables
37   USE trddyn         ! trend manager: dynamics
38   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
39   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
40   !
41   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl         ! Print control
43   USE in_out_manager ! I/O manager
44   USE lib_mpp        ! MPP library
45   USE timing         ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dyn_vor        ! routine called by step.F90
51   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
52
53   !                                   !!* Namelist namdyn_vor: vorticity term
54   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
55   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
56   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
57   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
58   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
59   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
60   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
61   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
62
63   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
64   !                                       ! associated indices:
65   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
66   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
67   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
68   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
69   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
70   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
71
72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
73   !                               ! associated indices:
74   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
75   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
76   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
77   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
78   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
79
80   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
82   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation
83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       -
84   
85   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
86   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
87   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
88   
89   !! * Substitutions
90#  include "do_loop_substitute.h90"
91#  include "domzgr_substitute.h90"
92
93   !!----------------------------------------------------------------------
94   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
95   !! $Id$
96   !! Software governed by the CeCILL license (see ./LICENSE)
97   !!----------------------------------------------------------------------
98CONTAINS
99
100   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
101      !!----------------------------------------------------------------------
102      !!
103      !! ** Purpose :   compute the lateral ocean tracer physics.
104      !!
105      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
106      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
107      !!               and planetary vorticity trends) and send them to trd_dyn
108      !!               for futher diagnostics (l_trddyn=T)
109      !!----------------------------------------------------------------------
110      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
111      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
112      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
113      !
114      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
115      !!----------------------------------------------------------------------
116      !
117      IF( ln_timing )   CALL timing_start('dyn_vor')
118      !
119      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
120         !
121         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
122         !
123         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
124         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
125         SELECT CASE( nvor_scheme )
126         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
127            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
128         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
129            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
130         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
131            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
132         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
133            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
134         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
135            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
136         END SELECT
137         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
138         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
139         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
140         !
141         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
142            ztrdu(:,:,:) = puu(:,:,:,Krhs)
143            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
144            SELECT CASE( nvor_scheme )
145            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
146            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
147            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
148            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
149            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
150            END SELECT
151            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
152            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
153            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
154         ENDIF
155         !
156         DEALLOCATE( ztrdu, ztrdv )
157         !
158      ELSE              !==  total vorticity trend added to the general trend  ==!
159         !
160         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
161         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
162                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
163            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
164         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
165                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
166            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
167         CASE( np_ENE )                        !* energy conserving scheme
168                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
169            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
170         CASE( np_ENS )                        !* enstrophy conserving scheme
171                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
172            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
173         CASE( np_MIX )                        !* mixed ene-ens scheme
174                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
175                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
176            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
177         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
178                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
179            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
180         END SELECT
181         !
182      ENDIF
183      !
184      !                       ! print sum trends (used for debugging)
185      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
186         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
187      !
188      IF( ln_timing )   CALL timing_stop('dyn_vor')
189      !
190   END SUBROUTINE dyn_vor
191
192
193   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
194      !!----------------------------------------------------------------------
195      !!                  ***  ROUTINE vor_enT  ***
196      !!
197      !! ** Purpose :   Compute the now total vorticity trend and add it to
198      !!      the general trend of the momentum equation.
199      !!
200      !! ** Method  :   Trend evaluated using now fields (centered in time)
201      !!       and t-point evaluation of vorticity (planetary and relative).
202      !!       conserves the horizontal kinetic energy.
203      !!         The general trend of momentum is increased due to the vorticity
204      !!       term which is given by:
205      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
206      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
207      !!       where rvor is the relative vorticity at f-point
208      !!
209      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
210      !!----------------------------------------------------------------------
211      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
212      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
213      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
214      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
215      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
216      !
217      INTEGER  ::   ji, jj, jk           ! dummy loop indices
218      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
219      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx, zwy, zwt   ! 2D workspace
220      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz             ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
221      !!----------------------------------------------------------------------
222      !
223      IF( kt == nit000 ) THEN
224         IF(lwp) WRITE(numout,*)
225         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
226         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
227      ENDIF
228      !
229      !
230      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
231      CASE ( np_RVO )                           !* relative vorticity
232         DO jk = 1, jpkm1                                 ! Horizontal slab
233            DO_2D( 1, 0, 1, 0 )
234               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
235                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
236            END_2D
237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
238               DO_2D( 1, 0, 1, 0 )
239                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
240               END_2D
241            ENDIF
242         END DO
243
244#if defined key_mpi3
245         CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
246#else
247         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
248#endif
249
250      CASE ( np_CRV )                           !* Coriolis + relative vorticity
251         DO jk = 1, jpkm1                                 ! Horizontal slab
252            DO_2D( 1, 0, 1, 0 )                          ! relative vorticity
253               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
254                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
255            END_2D
256            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
257               DO_2D( 1, 0, 1, 0 )
258                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
259               END_2D
260            ENDIF
261         END DO
262
263#if defined key_mpi3
264         CALL lbc_lnk_nc_multi( 'dynvor', zwz, 'F', 1.0_wp )
265#else
266         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
267#endif
268
269      END SELECT
270
271      !                                                ! ===============
272      DO jk = 1, jpkm1                                 ! Horizontal slab
273      !                                                ! ===============
274
275         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
276         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
277            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
278         CASE ( np_RVO )                           !* relative vorticity
279            DO_2D( 0, 1, 0, 1 )
280               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
281                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
282                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
283            END_2D
284         CASE ( np_MET )                           !* metric term
285            DO_2D( 0, 1, 0, 1 )
286               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
287                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
288                  &             * e3t(ji,jj,jk,Kmm)
289            END_2D
290         CASE ( np_CRV )                           !* Coriolis + relative vorticity
291            DO_2D( 0, 1, 0, 1 )
292               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
293                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
294                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
295            END_2D
296         CASE ( np_CME )                           !* Coriolis + metric
297            DO_2D( 0, 1, 0, 1 )
298               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
299                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
300                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
301                    &          * e3t(ji,jj,jk,Kmm)
302            END_2D
303         CASE DEFAULT                                             ! error
304            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
305         END SELECT
306         !
307         !                                   !==  compute and add the vorticity term trend  =!
308         DO_2D( 0, 0, 0, 0 )
309            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
310               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
311               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
312               !
313            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
314               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
315               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
316         END_2D
317         !                                             ! ===============
318      END DO                                           !   End of slab
319      !                                                ! ===============
320   END SUBROUTINE vor_enT
321
322
323   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
324      !!----------------------------------------------------------------------
325      !!                  ***  ROUTINE vor_ene  ***
326      !!
327      !! ** Purpose :   Compute the now total vorticity trend and add it to
328      !!      the general trend of the momentum equation.
329      !!
330      !! ** Method  :   Trend evaluated using now fields (centered in time)
331      !!       and the Sadourny (1975) flux form formulation : conserves the
332      !!       horizontal kinetic energy.
333      !!         The general trend of momentum is increased due to the vorticity
334      !!       term which is given by:
335      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
336      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
337      !!       where rvor is the relative vorticity
338      !!
339      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
340      !!
341      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
342      !!----------------------------------------------------------------------
343      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
344      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
345      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
346      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
347      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
348      !
349      INTEGER  ::   ji, jj, jk           ! dummy loop indices
350      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
351      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
352      !!----------------------------------------------------------------------
353      !
354      IF( kt == nit000 ) THEN
355         IF(lwp) WRITE(numout,*)
356         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
357         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
358      ENDIF
359      !
360      !                                                ! ===============
361      DO jk = 1, jpkm1                                 ! Horizontal slab
362         !                                             ! ===============
363         !
364         SELECT CASE( kvor )                 !==  vorticity considered  ==!
365         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
366            zwz(:,:) = ff_f(:,:) 
367         CASE ( np_RVO )                           !* relative vorticity
368            DO_2D( 1, 0, 1, 0 )
369               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
370                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
371            END_2D
372         CASE ( np_MET )                           !* metric term
373            DO_2D( 1, 0, 1, 0 )
374               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
375                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
376            END_2D
377         CASE ( np_CRV )                           !* Coriolis + relative vorticity
378            DO_2D( 1, 0, 1, 0 )
379               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
380                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
381            END_2D
382         CASE ( np_CME )                           !* Coriolis + metric
383            DO_2D( 1, 0, 1, 0 )
384               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
385                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
386            END_2D
387         CASE DEFAULT                                             ! error
388            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
389         END SELECT
390         !
391         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
392            DO_2D( 1, 0, 1, 0 )
393               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
394            END_2D
395         ENDIF
396
397         IF( ln_sco ) THEN
398            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
399            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
400            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
401         ELSE
402            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
403            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
404         ENDIF
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         CASE ( np_MET )                           !* metric term
469            DO_2D( 1, 0, 1, 0 )
470               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
471                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
472            END_2D
473         CASE ( np_CRV )                           !* Coriolis + relative vorticity
474            DO_2D( 1, 0, 1, 0 )
475               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
476                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
477            END_2D
478         CASE ( np_CME )                           !* Coriolis + metric
479            DO_2D( 1, 0, 1, 0 )
480               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
481                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
482            END_2D
483         CASE DEFAULT                                             ! error
484            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
485         END SELECT
486         !
487         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
488            DO_2D( 1, 0, 1, 0 )
489               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
490            END_2D
491         ENDIF
492         !
493         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
494            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
495            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
496            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
497         ELSE
498            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
499            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
500         ENDIF
501         !                                   !==  compute and add the vorticity term trend  =!
502         DO_2D( 0, 0, 0, 0 )
503            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
504               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
505            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
506               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
507            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
508            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
509         END_2D
510         !                                             ! ===============
511      END DO                                           !   End of slab
512      !                                                ! ===============
513   END SUBROUTINE vor_ens
514
515
516   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
517      !!----------------------------------------------------------------------
518      !!                ***  ROUTINE vor_een  ***
519      !!
520      !! ** Purpose :   Compute the now total vorticity trend and add it to
521      !!      the general trend of the momentum equation.
522      !!
523      !! ** Method  :   Trend evaluated using now fields (centered in time)
524      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
525      !!      both the horizontal kinetic energy and the potential enstrophy
526      !!      when horizontal divergence is zero (see the NEMO documentation)
527      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
528      !!
529      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
530      !!
531      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
532      !!----------------------------------------------------------------------
533      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
534      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
535      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
536      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
537      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
538      !
539      INTEGER  ::   ji, jj, jk   ! dummy loop indices
540      INTEGER  ::   ierr         ! local integer
541      REAL(wp) ::   zua, zva     ! local scalars
542      REAL(wp) ::   zmsk, ze3f   ! local scalars
543      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy , z1_e3f
544      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
545      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
546      !!----------------------------------------------------------------------
547      !
548      IF( kt == nit000 ) THEN
549         IF(lwp) WRITE(numout,*)
550         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
551         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
552      ENDIF
553      !
554      !                                                ! ===============
555      DO jk = 1, jpkm1                                 ! Horizontal slab
556         !                                             ! ===============
557         !
558         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
559         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
560            DO_2D( 1, 0, 1, 0 )
561               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
562                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
563                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
564                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
565               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
566               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
567               ENDIF
568            END_2D
569         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
570            DO_2D( 1, 0, 1, 0 )
571               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
572                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
573                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
574                  &    + 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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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.0_wp )
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( 0, 0, 0, 0 )
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,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
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( 1, 0, 1, 0 )
706               zwz(ji,jj,jk) = ff_f(ji,jj)
707            END_2D
708         CASE ( np_RVO )                           !* relative vorticity
709            DO_2D( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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( 1, 0, 1, 0 )
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.0_wp )
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( 0, 0, 0, 0 )
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( 1, 0, 1, 0, 1, jpk )
829            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
830               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
831         END_3D
832         !
833#if defined key_mpi3
834         CALL lbc_lnk_nc_multi( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
835#else
836         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
837#endif
838         !
839      ENDIF
840!!gm end
841
842      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
843      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
844      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
845      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
846      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
847      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
848      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
849      !
850      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
851      !                     
852      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
853      ncor = np_COR                       ! planetary vorticity
854      SELECT CASE( n_dynadv )
855      CASE( np_LIN_dyn )
856         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
857         nrvm = np_COR        ! planetary vorticity
858         ntot = np_COR        !     -         -
859      CASE( np_VEC_c2  )
860         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
861         nrvm = np_RVO        ! relative vorticity
862         ntot = np_CRV        ! relative + planetary vorticity         
863      CASE( np_FLX_c2 , np_FLX_ubs  )
864         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
865         nrvm = np_MET        ! metric term
866         ntot = np_CME        ! Coriolis + metric term
867         !
868         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
869         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
870            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
871            DO_2D( 0, 0, 0, 0 )
872               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
873               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
874            END_2D
875#if defined key_mpi3
876            CALL lbc_lnk_nc_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
877#else
878            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
879#endif
880            !
881         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
882            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
883            DO_2D( 1, 0, 1, 0 )
884               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
885               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
886            END_2D
887#if defined key_mpi3
888            CALL lbc_lnk_nc_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
889#else
890            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
891#endif
892         END SELECT
893         !
894      END SELECT
895     
896      IF(lwp) THEN                   ! Print the choice
897         WRITE(numout,*)
898         SELECT CASE( nvor_scheme )
899         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
900         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
901         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
902         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
903         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
904         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
905         END SELECT         
906      ENDIF
907      !
908   END SUBROUTINE dyn_vor_init
909
910   !!==============================================================================
911END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.