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

Last change on this file since 12377 was 12377, checked in by acc, 7 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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