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/releases/r4.0/r4.0-HEAD/src/OCE/DYN – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynvor.F90 @ 15795

Last change on this file since 15795 was 14596, checked in by clem, 3 years ago

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