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/UKMO/dev_r10037_GPU/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/dev_r10037_GPU/src/OCE/DYN/dynvor.F90 @ 11467

Last change on this file since 11467 was 11467, checked in by andmirek, 5 years ago

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