New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynvor.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

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

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

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