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

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

first Shallow Water Eq. update

  • 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!!an                  &                        * fmask(ji,jj,jk)
464            END_2D
465         CASE ( np_CME )                           !* Coriolis + metric
466            DO_2D_10_10
467               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
468                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
469            END_2D
470         CASE DEFAULT                                             ! error
471            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
472         END SELECT
473         !
474         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
475            DO_2D_10_10
476               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
477            END_2D
478         ENDIF
479         !
480         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
481            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
482            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
483            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
484         ELSE
485            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
486            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
487         ENDIF
488         !                                   !==  compute and add the vorticity term trend  =!
489         DO_2D_00_00
490            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
491               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
492            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
493               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
494            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
495            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
496         END_2D
497         !                                             ! ===============
498      END DO                                           !   End of slab
499      !                                                ! ===============
500   END SUBROUTINE vor_ens
501
502
503   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
504      !!----------------------------------------------------------------------
505      !!                ***  ROUTINE vor_een  ***
506      !!
507      !! ** Purpose :   Compute the now total vorticity trend and add it to
508      !!      the general trend of the momentum equation.
509      !!
510      !! ** Method  :   Trend evaluated using now fields (centered in time)
511      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
512      !!      both the horizontal kinetic energy and the potential enstrophy
513      !!      when horizontal divergence is zero (see the NEMO documentation)
514      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
515      !!
516      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
517      !!
518      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
519      !!----------------------------------------------------------------------
520      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
521      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
522      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
523      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
524      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
525      !
526      INTEGER  ::   ji, jj, jk   ! dummy loop indices
527      INTEGER  ::   ierr         ! local integer
528      REAL(wp) ::   zua, zva     ! local scalars
529      REAL(wp) ::   zmsk, ze3f   ! local scalars
530      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
531      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
532      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
533      !!----------------------------------------------------------------------
534      !
535      IF( kt == nit000 ) THEN
536         IF(lwp) WRITE(numout,*)
537         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
538         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
539      ENDIF
540      !
541      !                                                ! ===============
542      DO jk = 1, jpkm1                                 ! Horizontal slab
543         !                                             ! ===============
544         !
545         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
546         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
547            DO_2D_10_10
548               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)   &
549                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
550               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
551               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
552               ENDIF
553            END_2D
554         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
555            DO_2D_10_10
556               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)   &
557                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
558               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
559                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
560               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
561               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
562               ENDIF
563            END_2D
564         END SELECT
565         !
566         SELECT CASE( kvor )                 !==  vorticity considered  ==!
567         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
568            DO_2D_10_10
569               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
570            END_2D
571         CASE ( np_RVO )                           !* relative vorticity
572            DO_2D_10_10
573               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
574                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
575            END_2D
576         CASE ( np_MET )                           !* metric term
577            DO_2D_10_10
578               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
579                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
580            END_2D
581         CASE ( np_CRV )                           !* Coriolis + relative vorticity
582            DO_2D_10_10
583               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
584                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
585                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
586            END_2D
587         CASE ( np_CME )                           !* Coriolis + metric
588            DO_2D_10_10
589               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
590                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
591            END_2D
592         CASE DEFAULT                                             ! error
593            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
594         END SELECT
595         !
596         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
597            DO_2D_10_10
598               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
599            END_2D
600         ENDIF
601      END DO                                           !   End of slab
602         !
603      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
604
605      DO jk = 1, jpkm1                                 ! Horizontal slab
606         !
607         !                                   !==  horizontal fluxes  ==!
608         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
609         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
610
611         !                                   !==  compute and add the vorticity term trend  =!
612         jj = 2
613         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
614         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
615               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
616               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
617               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
618               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
619         END DO
620         DO jj = 3, jpj
621            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
622               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
623               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
624               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
625               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
626            END DO
627         END DO
628         DO_2D_00_00
629            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
630               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
631            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
632               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
633            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
634            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
635         END_2D
636         !                                             ! ===============
637      END DO                                           !   End of slab
638      !                                                ! ===============
639   END SUBROUTINE vor_een
640
641
642
643   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
644      !!----------------------------------------------------------------------
645      !!                ***  ROUTINE vor_eeT  ***
646      !!
647      !! ** Purpose :   Compute the now total vorticity trend and add it to
648      !!      the general trend of the momentum equation.
649      !!
650      !! ** Method  :   Trend evaluated using now fields (centered in time)
651      !!      and the Arakawa and Lamb (1980) vector form formulation using
652      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
653      !!      The change consists in
654      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
655      !!
656      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
657      !!
658      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
659      !!----------------------------------------------------------------------
660      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
661      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
662      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
663      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
664      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
665      !
666      INTEGER  ::   ji, jj, jk     ! dummy loop indices
667      INTEGER  ::   ierr           ! local integer
668      REAL(wp) ::   zua, zva       ! local scalars
669      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
670      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
671      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
672      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
673      !!----------------------------------------------------------------------
674      !
675      IF( kt == nit000 ) THEN
676         IF(lwp) WRITE(numout,*)
677         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
678         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
679      ENDIF
680      !
681      !                                                ! ===============
682      DO jk = 1, jpkm1                                 ! Horizontal slab
683         !                                             ! ===============
684         !
685         !
686         SELECT CASE( kvor )                 !==  vorticity considered  ==!
687         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
688            DO_2D_10_10
689               zwz(ji,jj,jk) = ff_f(ji,jj)
690            END_2D
691         CASE ( np_RVO )                           !* relative vorticity
692            DO_2D_10_10
693               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
694                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
695                  &          * r1_e1e2f(ji,jj)
696            END_2D
697         CASE ( np_MET )                           !* metric term
698            DO_2D_10_10
699               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
700                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
701            END_2D
702         CASE ( np_CRV )                           !* Coriolis + relative vorticity
703            DO_2D_10_10
704               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
705                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
706                  &                         * r1_e1e2f(ji,jj)    )
707            END_2D
708         CASE ( np_CME )                           !* Coriolis + metric
709            DO_2D_10_10
710               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
711                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
712            END_2D
713         CASE DEFAULT                                             ! error
714            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
715         END SELECT
716         !
717         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
718            DO_2D_10_10
719               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
720            END_2D
721         ENDIF
722      END DO
723      !
724      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
725      !
726      DO jk = 1, jpkm1                                 ! Horizontal slab
727
728      !                                   !==  horizontal fluxes  ==!
729         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
730         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
731
732         !                                   !==  compute and add the vorticity term trend  =!
733         jj = 2
734         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
735         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
736               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
737               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
738               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
739               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
740               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
741         END DO
742         DO jj = 3, jpj
743            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
744               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
745               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
746               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
747               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
748               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
749            END DO
750         END DO
751         DO_2D_00_00
752            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
753               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
754            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
755               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
756            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
757            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
758         END_2D
759         !                                             ! ===============
760      END DO                                           !   End of slab
761      !                                                ! ===============
762   END SUBROUTINE vor_eeT
763
764
765   SUBROUTINE dyn_vor_init
766      !!---------------------------------------------------------------------
767      !!                  ***  ROUTINE dyn_vor_init  ***
768      !!
769      !! ** Purpose :   Control the consistency between cpp options for
770      !!              tracer advection schemes
771      !!----------------------------------------------------------------------
772      INTEGER ::   ji, jj, jk    ! dummy loop indices
773      INTEGER ::   ioptio, ios   ! local integer
774      !!
775      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
776         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
777      !!----------------------------------------------------------------------
778      !
779      IF(lwp) THEN
780         WRITE(numout,*)
781         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
782         WRITE(numout,*) '~~~~~~~~~~~~'
783      ENDIF
784      !
785      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
786901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
787      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
788902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
789      IF(lwm) WRITE ( numond, namdyn_vor )
790      !
791      IF(lwp) THEN                    ! Namelist print
792         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
793         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
794         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
795         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
796         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
797         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
798         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
799         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
800         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
801      ENDIF
802
803      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
804
805!!gm  this should be removed when choosing a unique strategy for fmask at the coast
806      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
807      ! at angles with three ocean points and one land point
808      IF(lwp) WRITE(numout,*)
809      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
810      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
811         DO_3D_10_10( 1, jpk )
812            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
813               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
814         END_3D
815         !
816         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
817         !
818      ENDIF
819!!gm end
820
821      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
822      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
823      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
824      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
825      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
826      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
827      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
828      !
829      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
830      !                     
831      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
832      ncor = np_COR                       ! planetary vorticity
833      SELECT CASE( n_dynadv )
834      CASE( np_LIN_dyn )
835         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
836         nrvm = np_COR        ! planetary vorticity
837         ntot = np_COR        !     -         -
838      CASE( np_VEC_c2  )
839         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
840         nrvm = np_RVO        ! relative vorticity
841         ntot = np_CRV        ! relative + planetary vorticity         
842      CASE( np_FLX_c2 , np_FLX_ubs  )
843         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
844         nrvm = np_MET        ! metric term
845         ntot = np_CME        ! Coriolis + metric term
846         !
847         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
848         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
849            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
850            DO_2D_00_00
851               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
852               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
853            END_2D
854            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
855            !
856         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
857            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
858            DO_2D_10_10
859               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
860               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
861            END_2D
862            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
863         END SELECT
864         !
865      END SELECT
866     
867      IF(lwp) THEN                   ! Print the choice
868         WRITE(numout,*)
869         SELECT CASE( nvor_scheme )
870         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
871         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
872         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
873         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
874         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
875         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
876         END SELECT         
877      ENDIF
878      !
879   END SUBROUTINE dyn_vor_init
880
881   !!==============================================================================
882END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.