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/2019/dev_r11943_MERGE_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynvor.F90 @ 12236

Last change on this file since 12236 was 12236, checked in by acc, 4 years ago

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