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

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

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SWE/dynvor.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

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