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

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

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90 @ 12667

Last change on this file since 12667 was 12667, checked in by gm, 4 years ago

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