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

source: NEMO/trunk/src/OCE/DYN/dynvor.F90 @ 13546

Last change on this file since 13546 was 13546, checked in by smasson, 4 years ago

trunk: supress the use of undefined values, see #2535

  • Property svn:keywords set to Id
File size: 53.1 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity
18   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory
19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!----------------------------------------------------------------------
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         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
245
246      CASE ( np_CRV )                           !* Coriolis + relative vorticity
247         DO jk = 1, jpkm1                                 ! Horizontal slab
248            DO_2D( 1, 0, 1, 0 )                          ! relative vorticity
249               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
250                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
251            END_2D
252            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
253               DO_2D( 1, 0, 1, 0 )
254                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
255               END_2D
256            ENDIF
257         END DO
258
259         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
260
261      END SELECT
262
263      !                                                ! ===============
264      DO jk = 1, jpkm1                                 ! Horizontal slab
265      !                                                ! ===============
266
267         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
268         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
269            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
270         CASE ( np_RVO )                           !* relative vorticity
271            DO_2D( 0, 1, 0, 1 )
272               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
273                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
274                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
275            END_2D
276         CASE ( np_MET )                           !* metric term
277            DO_2D( 0, 1, 0, 1 )
278               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
279                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
280                  &             * e3t(ji,jj,jk,Kmm)
281            END_2D
282         CASE ( np_CRV )                           !* Coriolis + relative vorticity
283            DO_2D( 0, 1, 0, 1 )
284               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
285                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
286                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
287            END_2D
288         CASE ( np_CME )                           !* Coriolis + metric
289            DO_2D( 0, 1, 0, 1 )
290               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
291                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
292                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
293                    &          * e3t(ji,jj,jk,Kmm)
294            END_2D
295         CASE DEFAULT                                             ! error
296            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
297         END SELECT
298         !
299         !                                   !==  compute and add the vorticity term trend  =!
300         DO_2D( 0, 0, 0, 0 )
301            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
302               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
303               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
304               !
305            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
306               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
307               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
308         END_2D
309         !                                             ! ===============
310      END DO                                           !   End of slab
311      !                                                ! ===============
312   END SUBROUTINE vor_enT
313
314
315   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
316      !!----------------------------------------------------------------------
317      !!                  ***  ROUTINE vor_ene  ***
318      !!
319      !! ** Purpose :   Compute the now total vorticity trend and add it to
320      !!      the general trend of the momentum equation.
321      !!
322      !! ** Method  :   Trend evaluated using now fields (centered in time)
323      !!       and the Sadourny (1975) flux form formulation : conserves the
324      !!       horizontal kinetic energy.
325      !!         The general trend of momentum is increased due to the vorticity
326      !!       term which is given by:
327      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
328      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
329      !!       where rvor is the relative vorticity
330      !!
331      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
332      !!
333      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
334      !!----------------------------------------------------------------------
335      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
336      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
337      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
338      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
339      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
340      !
341      INTEGER  ::   ji, jj, jk           ! dummy loop indices
342      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
343      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
344      !!----------------------------------------------------------------------
345      !
346      IF( kt == nit000 ) THEN
347         IF(lwp) WRITE(numout,*)
348         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
349         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
350      ENDIF
351      !
352      !                                                ! ===============
353      DO jk = 1, jpkm1                                 ! Horizontal slab
354         !                                             ! ===============
355         !
356         SELECT CASE( kvor )                 !==  vorticity considered  ==!
357         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
358            zwz(:,:) = ff_f(:,:) 
359         CASE ( np_RVO )                           !* relative vorticity
360            DO_2D( 1, 0, 1, 0 )
361               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
362                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
363            END_2D
364         CASE ( np_MET )                           !* metric term
365            DO_2D( 1, 0, 1, 0 )
366               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
367                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
368            END_2D
369         CASE ( np_CRV )                           !* Coriolis + relative vorticity
370            DO_2D( 1, 0, 1, 0 )
371               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
372                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
373            END_2D
374         CASE ( np_CME )                           !* Coriolis + metric
375            DO_2D( 1, 0, 1, 0 )
376               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
377                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
378            END_2D
379         CASE DEFAULT                                             ! error
380            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
381         END SELECT
382         !
383         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
384            DO_2D( 1, 0, 1, 0 )
385               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
386            END_2D
387         ENDIF
388
389         IF( ln_sco ) THEN
390            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
391            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
392            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
393         ELSE
394            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
395            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
396         ENDIF
397         !                                   !==  compute and add the vorticity term trend  =!
398         DO_2D( 0, 0, 0, 0 )
399            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
400            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
401            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
402            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
403            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 )
404            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 ) 
405         END_2D
406         !                                             ! ===============
407      END DO                                           !   End of slab
408      !                                                ! ===============
409   END SUBROUTINE vor_ene
410
411
412   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
413      !!----------------------------------------------------------------------
414      !!                ***  ROUTINE vor_ens  ***
415      !!
416      !! ** Purpose :   Compute the now total vorticity trend and add it to
417      !!      the general trend of the momentum equation.
418      !!
419      !! ** Method  :   Trend evaluated using now fields (centered in time)
420      !!      and the Sadourny (1975) flux FORM formulation : conserves the
421      !!      potential enstrophy of a horizontally non-divergent flow. the
422      !!      trend of the vorticity term is given by:
423      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
424      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
425      !!      Add this trend to the general momentum trend:
426      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
427      !!
428      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
429      !!
430      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
431      !!----------------------------------------------------------------------
432      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
433      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
434      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
435      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
436      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
437      !
438      INTEGER  ::   ji, jj, jk   ! dummy loop indices
439      REAL(wp) ::   zuav, zvau   ! local scalars
440      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
441      !!----------------------------------------------------------------------
442      !
443      IF( kt == nit000 ) THEN
444         IF(lwp) WRITE(numout,*)
445         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
446         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
447      ENDIF
448      !                                                ! ===============
449      DO jk = 1, jpkm1                                 ! Horizontal slab
450         !                                             ! ===============
451         !
452         SELECT CASE( kvor )                 !==  vorticity considered  ==!
453         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
454            zwz(:,:) = ff_f(:,:) 
455         CASE ( np_RVO )                           !* relative vorticity
456            DO_2D( 1, 0, 1, 0 )
457               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
458                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
459            END_2D
460         CASE ( np_MET )                           !* metric term
461            DO_2D( 1, 0, 1, 0 )
462               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
463                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
464            END_2D
465         CASE ( np_CRV )                           !* Coriolis + relative vorticity
466            DO_2D( 1, 0, 1, 0 )
467               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
468                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
469            END_2D
470         CASE ( np_CME )                           !* Coriolis + metric
471            DO_2D( 1, 0, 1, 0 )
472               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
473                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
474            END_2D
475         CASE DEFAULT                                             ! error
476            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
477         END SELECT
478         !
479         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
480            DO_2D( 1, 0, 1, 0 )
481               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
482            END_2D
483         ENDIF
484         !
485         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
486            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
487            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
488            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
489         ELSE
490            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
491            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
492         ENDIF
493         !                                   !==  compute and add the vorticity term trend  =!
494         DO_2D( 0, 0, 0, 0 )
495            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
496               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
497            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
498               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
499            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
500            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
501         END_2D
502         !                                             ! ===============
503      END DO                                           !   End of slab
504      !                                                ! ===============
505   END SUBROUTINE vor_ens
506
507
508   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
509      !!----------------------------------------------------------------------
510      !!                ***  ROUTINE vor_een  ***
511      !!
512      !! ** Purpose :   Compute the now total vorticity trend and add it to
513      !!      the general trend of the momentum equation.
514      !!
515      !! ** Method  :   Trend evaluated using now fields (centered in time)
516      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
517      !!      both the horizontal kinetic energy and the potential enstrophy
518      !!      when horizontal divergence is zero (see the NEMO documentation)
519      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
520      !!
521      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
522      !!
523      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
524      !!----------------------------------------------------------------------
525      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
526      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
527      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
528      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
529      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
530      !
531      INTEGER  ::   ji, jj, jk   ! dummy loop indices
532      INTEGER  ::   ierr         ! local integer
533      REAL(wp) ::   zua, zva     ! local scalars
534      REAL(wp) ::   zmsk, ze3f   ! local scalars
535      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy , z1_e3f
536      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
537      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
538      !!----------------------------------------------------------------------
539      !
540      IF( kt == nit000 ) THEN
541         IF(lwp) WRITE(numout,*)
542         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
543         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
544      ENDIF
545      !
546      !                                                ! ===============
547      DO jk = 1, jpkm1                                 ! Horizontal slab
548         !                                             ! ===============
549         !
550         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
551         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
552            DO_2D( 1, 0, 1, 0 )
553               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
554                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
555                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
556                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
557               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
558               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
559               ENDIF
560            END_2D
561         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
562            DO_2D( 1, 0, 1, 0 )
563               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
564                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
565                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
566                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
567               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
568                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
569               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
570               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
571               ENDIF
572            END_2D
573         END SELECT
574         !
575         SELECT CASE( kvor )                 !==  vorticity considered  ==!
576         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
577            DO_2D( 1, 0, 1, 0 )
578               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
579            END_2D
580         CASE ( np_RVO )                           !* relative vorticity
581            DO_2D( 1, 0, 1, 0 )
582               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
583                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
584            END_2D
585         CASE ( np_MET )                           !* metric term
586            DO_2D( 1, 0, 1, 0 )
587               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
588                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
589            END_2D
590         CASE ( np_CRV )                           !* Coriolis + relative vorticity
591            DO_2D( 1, 0, 1, 0 )
592               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
593                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
594                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
595            END_2D
596         CASE ( np_CME )                           !* Coriolis + metric
597            DO_2D( 1, 0, 1, 0 )
598               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
599                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
600            END_2D
601         CASE DEFAULT                                             ! error
602            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
603         END SELECT
604         !
605         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
606            DO_2D( 1, 0, 1, 0 )
607               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
608            END_2D
609         ENDIF
610      END DO                                           !   End of slab
611         !
612      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
613
614      DO jk = 1, jpkm1                                 ! Horizontal slab
615         !
616         !                                   !==  horizontal fluxes  ==!
617         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
618         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
619
620         !                                   !==  compute and add the vorticity term trend  =!
621         jj = 2
622         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
623         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
624               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
625               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
626               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
627               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
628         END DO
629         DO jj = 3, jpj
630            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
631               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
632               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
633               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
634               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
635            END DO
636         END DO
637         DO_2D( 0, 0, 0, 0 )
638            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
639               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
640            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
641               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
642            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
643            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
644         END_2D
645         !                                             ! ===============
646      END DO                                           !   End of slab
647      !                                                ! ===============
648   END SUBROUTINE vor_een
649
650
651
652   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
653      !!----------------------------------------------------------------------
654      !!                ***  ROUTINE vor_eeT  ***
655      !!
656      !! ** Purpose :   Compute the now total vorticity trend and add it to
657      !!      the general trend of the momentum equation.
658      !!
659      !! ** Method  :   Trend evaluated using now fields (centered in time)
660      !!      and the Arakawa and Lamb (1980) vector form formulation using
661      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
662      !!      The change consists in
663      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
664      !!
665      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
666      !!
667      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
668      !!----------------------------------------------------------------------
669      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
670      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
671      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
672      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
673      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
674      !
675      INTEGER  ::   ji, jj, jk     ! dummy loop indices
676      INTEGER  ::   ierr           ! local integer
677      REAL(wp) ::   zua, zva       ! local scalars
678      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
679      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy 
680      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse
681      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined
682      !!----------------------------------------------------------------------
683      !
684      IF( kt == nit000 ) THEN
685         IF(lwp) WRITE(numout,*)
686         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
687         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
688      ENDIF
689      !
690      !                                                ! ===============
691      DO jk = 1, jpkm1                                 ! Horizontal slab
692         !                                             ! ===============
693         !
694         !
695         SELECT CASE( kvor )                 !==  vorticity considered  ==!
696         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
697            DO_2D( 1, 0, 1, 0 )
698               zwz(ji,jj,jk) = ff_f(ji,jj)
699            END_2D
700         CASE ( np_RVO )                           !* relative vorticity
701            DO_2D( 1, 0, 1, 0 )
702               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
703                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
704                  &          * r1_e1e2f(ji,jj)
705            END_2D
706         CASE ( np_MET )                           !* metric term
707            DO_2D( 1, 0, 1, 0 )
708               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
709                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
710            END_2D
711         CASE ( np_CRV )                           !* Coriolis + relative vorticity
712            DO_2D( 1, 0, 1, 0 )
713               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
714                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
715                  &                         * r1_e1e2f(ji,jj)    )
716            END_2D
717         CASE ( np_CME )                           !* Coriolis + metric
718            DO_2D( 1, 0, 1, 0 )
719               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
720                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
721            END_2D
722         CASE DEFAULT                                             ! error
723            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
724         END SELECT
725         !
726         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
727            DO_2D( 1, 0, 1, 0 )
728               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
729            END_2D
730         ENDIF
731      END DO
732      !
733      CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
734      !
735      DO jk = 1, jpkm1                                 ! Horizontal slab
736
737      !                                   !==  horizontal fluxes  ==!
738         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
739         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
740
741         !                                   !==  compute and add the vorticity term trend  =!
742         jj = 2
743         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
744         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
745               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
746               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
747               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
748               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
749               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
750         END DO
751         DO jj = 3, jpj
752            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
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         END DO
760         DO_2D( 0, 0, 0, 0 )
761            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
762               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
763            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
764               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
765            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
766            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
767         END_2D
768         !                                             ! ===============
769      END DO                                           !   End of slab
770      !                                                ! ===============
771   END SUBROUTINE vor_eeT
772
773
774   SUBROUTINE dyn_vor_init
775      !!---------------------------------------------------------------------
776      !!                  ***  ROUTINE dyn_vor_init  ***
777      !!
778      !! ** Purpose :   Control the consistency between cpp options for
779      !!              tracer advection schemes
780      !!----------------------------------------------------------------------
781      INTEGER ::   ji, jj, jk    ! dummy loop indices
782      INTEGER ::   ioptio, ios   ! local integer
783      !!
784      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
785         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
786      !!----------------------------------------------------------------------
787      !
788      IF(lwp) THEN
789         WRITE(numout,*)
790         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
791         WRITE(numout,*) '~~~~~~~~~~~~'
792      ENDIF
793      !
794      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
795901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
796      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
797902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
798      IF(lwm) WRITE ( numond, namdyn_vor )
799      !
800      IF(lwp) THEN                    ! Namelist print
801         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
802         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
803         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
804         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
805         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
806         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
807         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
808         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
809         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
810      ENDIF
811
812      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
813
814!!gm  this should be removed when choosing a unique strategy for fmask at the coast
815      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
816      ! at angles with three ocean points and one land point
817      IF(lwp) WRITE(numout,*)
818      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
819      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
820         DO_3D( 1, 0, 1, 0, 1, jpk )
821            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
822               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
823         END_3D
824         !
825         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
826         !
827      ENDIF
828!!gm end
829
830      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
831      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
832      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
833      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
834      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
835      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
836      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
837      !
838      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
839      !                     
840      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
841      ncor = np_COR                       ! planetary vorticity
842      SELECT CASE( n_dynadv )
843      CASE( np_LIN_dyn )
844         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
845         nrvm = np_COR        ! planetary vorticity
846         ntot = np_COR        !     -         -
847      CASE( np_VEC_c2  )
848         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
849         nrvm = np_RVO        ! relative vorticity
850         ntot = np_CRV        ! relative + planetary vorticity         
851      CASE( np_FLX_c2 , np_FLX_ubs  )
852         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
853         nrvm = np_MET        ! metric term
854         ntot = np_CME        ! Coriolis + metric term
855         !
856         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
857         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
858            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
859            DO_2D( 0, 0, 0, 0 )
860               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
861               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
862            END_2D
863            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp )   ! Lateral boundary conditions
864            !
865         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
866            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
867            DO_2D( 1, 0, 1, 0 )
868               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
869               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
870            END_2D
871            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp )   ! Lateral boundary conditions
872         END SELECT
873         !
874      END SELECT
875     
876      IF(lwp) THEN                   ! Print the choice
877         WRITE(numout,*)
878         SELECT CASE( nvor_scheme )
879         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
880         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
881         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
882         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
883         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
884         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
885         END SELECT         
886      ENDIF
887      !
888   END SUBROUTINE dyn_vor_init
889
890   !!==============================================================================
891END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.