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 @ 14007

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

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