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/trunk/src/SWE – NEMO

source: NEMO/trunk/src/SWE/dynvor.F90 @ 13237

Last change on this file since 13237 was 12983, checked in by techene, 4 years ago

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