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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 4460

Last change on this file since 4460 was 4450, checked in by trackstand2, 10 years ago

Add mbkmax in dynvor
  • Property svn:keywords set to Id
File size: 45.3 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec) vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   dyn_vor      : Update the momentum trend with the vorticity trend
21   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
22   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
23   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
24   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
25   !!   dyn_vor_init : set and control of the different vorticity option
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and tracers
28   USE dom_oce        ! ocean space and time domain
29   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
30   USE trdmod         ! ocean dynamics trends
31   USE trdmod_oce     ! ocean variables trends
32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33   USE prtctl         ! Print control
34   USE in_out_manager ! I/O manager
35   USE lib_mpp
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   dyn_vor        ! routine called by step.F90
41   PUBLIC   dyn_vor_init   ! routine called by opa.F90
42
43   !                                             !!* Namelist namdyn_vor: vorticity term
44   LOGICAL, PUBLIC ::   ln_dynvor_ene = .FALSE.   !: energy conserving scheme
45   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
46   LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
47   LOGICAL, PUBLIC ::   ln_dynvor_een = .FALSE.   !: energy and enstrophy conserving scheme
48
49   INTEGER ::   nvor = 0   ! type of vorticity trend used
50   INTEGER ::   ncor = 1   ! coriolis
51   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
52   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
53
54   !! * Control permutation of array indices
55#  include "oce_ftrans.h90"
56#  include "dom_oce_ftrans.h90"
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE dyn_vor( kt )
69      !!----------------------------------------------------------------------
70      !!
71      !! ** Purpose :   compute the lateral ocean tracer physics.
72      !!
73      !! ** Action : - Update (ua,va) with the now vorticity term trend
74      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
75      !!               and planetary vorticity trends) ('key_trddyn')
76      !!----------------------------------------------------------------------
77      USE oce, ONLY:   ztrdu => ta , ztrdv => sa   ! (ta,sa) used as 3D workspace
78      !! DCSE_NEMO: need additional directives for renamed module variables
79!FTRANS ztrdu :I :I :z
80!FTRANS ztrdv :I :I :z
81
82      !
83      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
84      !!----------------------------------------------------------------------
85      !
86      !                                          ! vorticity term
87      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
88      !
89      CASE ( -1 )                                      ! esopa: test all possibility with control print
90         CALL vor_ene( kt, ntot, ua, va )
91         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
92            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
93         CALL vor_ens( kt, ntot, ua, va )
94         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
95            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
96         CALL vor_mix( kt )
97         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
98            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
99         CALL vor_een( kt, ntot, ua, va )
100         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
101            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
102         !
103      CASE ( 0 )                                       ! energy conserving scheme
104         IF( l_trddyn )   THEN
105            ztrdu(:,:,:) = ua(:,:,:)
106            ztrdv(:,:,:) = va(:,:,:)
107            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
108            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
109            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
110            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
111            ztrdu(:,:,:) = ua(:,:,:)
112            ztrdv(:,:,:) = va(:,:,:)
113            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
114            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
115            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
116            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
117            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
118         ELSE
119            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
120         ENDIF
121         !
122      CASE ( 1 )                                       ! enstrophy conserving scheme
123         IF( l_trddyn )   THEN   
124            ztrdu(:,:,:) = ua(:,:,:)
125            ztrdv(:,:,:) = va(:,:,:)
126            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
127            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
128            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
129            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
130            ztrdu(:,:,:) = ua(:,:,:)
131            ztrdv(:,:,:) = va(:,:,:)
132            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
133            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
134            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
135            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
136            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
137         ELSE
138            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
139         ENDIF
140         !
141      CASE ( 2 )                                       ! mixed ene-ens scheme
142         IF( l_trddyn )   THEN
143            ztrdu(:,:,:) = ua(:,:,:)
144            ztrdv(:,:,:) = va(:,:,:)
145            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
146            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
147            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
148            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
149            ztrdu(:,:,:) = ua(:,:,:)
150            ztrdv(:,:,:) = va(:,:,:)
151            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
152            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
153            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
154            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
155            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
156         ELSE
157            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
158         ENDIF
159         !
160      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
161         IF( l_trddyn )   THEN
162            ztrdu(:,:,:) = ua(:,:,:)
163            ztrdv(:,:,:) = va(:,:,:)
164            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
165            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
166            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
167            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
168            ztrdu(:,:,:) = ua(:,:,:)
169            ztrdv(:,:,:) = va(:,:,:)
170            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
171            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
172            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
173            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
174            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
175         ELSE
176            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
177         ENDIF
178         !
179      END SELECT
180      !
181      !                       ! print sum trends (used for debugging)
182      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
183         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
184      !
185   END SUBROUTINE dyn_vor
186
187
188   SUBROUTINE vor_ene( kt, kvor, pua, pva )
189      !!----------------------------------------------------------------------
190      !!                  ***  ROUTINE vor_ene  ***
191      !!
192      !! ** Purpose :   Compute the now total vorticity trend and add it to
193      !!      the general trend of the momentum equation.
194      !!
195      !! ** Method  :   Trend evaluated using now fields (centered in time)
196      !!      and the Sadourny (1975) flux form formulation : conserves the
197      !!      horizontal kinetic energy.
198      !!      The trend of the vorticity term is given by:
199      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
200      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
201      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
202      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
203      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
204      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
205      !!      Add this trend to the general momentum trend (ua,va):
206      !!          (ua,va) = (ua,va) + ( voru , vorv )
207      !!
208      !! ** Action : - Update (ua,va) with the now vorticity term trend
209      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
210      !!               and planetary vorticity trends) ('key_trddyn')
211      !!
212      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
213      !!----------------------------------------------------------------------
214      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
215      USE wrk_nemo, ONLY:   zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3     ! 2D workspace
216      !
217      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
218      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
219      !                                                           ! =nrvm (relative vorticity or metric)
220!FTRANS pua :I :I :z
221!FTRANS pva :I :I :z
222!! DCSE_NEMO: work around a deficiency in ftrans
223!     REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
224!     REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
225      REAL(wp), INTENT(inout) ::   pua(jpi,jpj,jpkorig)    ! total u-trend
226      REAL(wp), INTENT(inout) ::   pva(jpi,jpj,jpkorig)    ! total v-trend
227      !
228      INTEGER  ::   ji, jj, jk   ! dummy loop indices
229      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
230      !!----------------------------------------------------------------------
231
232      IF( wrk_in_use(2, 1,2,3) ) THEN
233         CALL ctl_stop('dyn:vor_ene: requested workspace arrays unavailable')   ;   RETURN
234      ENDIF
235
236      IF( kt == nit000 ) THEN
237         IF(lwp) WRITE(numout,*)
238         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
239         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
240      ENDIF
241
242      zfact2 = 0.5 * 0.5      ! Local constant initialization
243
244!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
245      !                                                ! ===============
246      DO jk = 1, jpkm1                                 ! Horizontal slab
247         !                                             ! ===============
248         !
249         ! Potential vorticity and horizontal fluxes
250         ! -----------------------------------------
251         SELECT CASE( kvor )      ! vorticity considered
252         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
253         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
254         CASE ( 3 )                                                ! metric term
255            DO jj = 1, jpjm1
256               DO ji = 1, fs_jpim1   ! vector opt.
257                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
258                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
259                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
260               END DO
261            END DO
262         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
263         CASE ( 5 )                                                ! total (coriolis + metric)
264            DO jj = 1, jpjm1
265               DO ji = 1, fs_jpim1   ! vector opt.
266                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
267                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
268                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
269                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
270                       &       )
271               END DO
272            END DO
273         END SELECT
274
275         IF( ln_sco ) THEN
276            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
277            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
278            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
279         ELSE
280            zwx(:,:) = e2u(:,:) * un(:,:,jk)
281            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
282         ENDIF
283
284         ! Compute and add the vorticity term trend
285         ! ----------------------------------------
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1   ! vector opt.
288               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
289               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
290               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
291               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
292               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
293               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
294            END DO 
295         END DO 
296         !                                             ! ===============
297      END DO                                           !   End of slab
298      !                                                ! ===============
299      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn:vor_ene: failed to release workspace arrays')
300      !
301   END SUBROUTINE vor_ene
302
303
304   SUBROUTINE vor_mix( kt )
305      !!----------------------------------------------------------------------
306      !!                 ***  ROUTINE vor_mix  ***
307      !!
308      !! ** Purpose :   Compute the now total vorticity trend and add it to
309      !!      the general trend of the momentum equation.
310      !!
311      !! ** Method  :   Trend evaluated using now fields (centered in time)
312      !!      Mixte formulation : conserves the potential enstrophy of a hori-
313      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
314      !!      ticity term and the horizontal kinetic energy for (f x uh), the
315      !!      coriolis term. the now trend of the vorticity term is given by:
316      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
317      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
318      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
319      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
320      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
321      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
322      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
323      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
324      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
325      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
326      !!      Add this now trend to the general momentum trend (ua,va):
327      !!          (ua,va) = (ua,va) + ( voru , vorv )
328      !!
329      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
330      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
331      !!               and planetary vorticity trends) ('key_trddyn')
332      !!
333      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
334      !!----------------------------------------------------------------------
335      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
336      USE wrk_nemo, ONLY:   zwx => wrk_2d_4 , zwy => wrk_2d_5 , zwz => wrk_2d_6 , zww => wrk_2d_7   ! 2D workspace
337      !
338      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
339      !
340      INTEGER  ::   ji, jj, jk   ! dummy loop indices
341      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
342      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
343      !!----------------------------------------------------------------------
344
345      IF( wrk_in_use(2, 4,5,6,7) ) THEN
346         CALL ctl_stop('dyn:vor_mix: requested workspace arrays unavailable')   ;   RETURN
347      ENDIF
348
349      IF( kt == nit000 ) THEN
350         IF(lwp) WRITE(numout,*)
351         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
352         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
353      ENDIF
354
355      zfact1 = 0.5 * 0.25      ! Local constant initialization
356      zfact2 = 0.5 * 0.5
357
358!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
359      !                                                ! ===============
360      DO jk = 1, jpkm1                                 ! Horizontal slab
361         !                                             ! ===============
362         !
363         ! Relative and planetary potential vorticity and horizontal fluxes
364         ! ----------------------------------------------------------------
365         IF( ln_sco ) THEN       
366            IF( ln_dynadv_vec ) THEN
367               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
368            ELSE                       
369               DO jj = 1, jpjm1
370                  DO ji = 1, fs_jpim1   ! vector opt.
371                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
372                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
373                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
374                  END DO
375               END DO
376            ENDIF
377            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
378            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
379            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
380         ELSE
381            IF( ln_dynadv_vec ) THEN
382               zww(:,:) = rotn(:,:,jk)
383            ELSE                       
384               DO jj = 1, jpjm1
385                  DO ji = 1, fs_jpim1   ! vector opt.
386                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
387                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
388                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
389                  END DO
390               END DO
391            ENDIF
392            zwz(:,:) = ff (:,:)
393            zwx(:,:) = e2u(:,:) * un(:,:,jk)
394            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
395         ENDIF
396
397         ! Compute and add the vorticity term trend
398         ! ----------------------------------------
399         DO jj = 2, jpjm1
400            DO ji = fs_2, fs_jpim1   ! vector opt.
401               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
402               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
403               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
404               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
405               ! enstrophy conserving formulation for relative vorticity term
406               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
407               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
408               ! energy conserving formulation for planetary vorticity term
409               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
410               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
411               ! mixed vorticity trend added to the momentum trends
412               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
413               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
414            END DO 
415         END DO 
416         !                                             ! ===============
417      END DO                                           !   End of slab
418      !                                                ! ===============
419      IF( wrk_not_released(2, 4,5,6,7) )   CALL ctl_stop('dyn:vor_mix: failed to release workspace arrays')
420      !
421   END SUBROUTINE vor_mix
422
423
424   SUBROUTINE vor_ens( kt, kvor, pua, pva )
425      !!----------------------------------------------------------------------
426      !!                ***  ROUTINE vor_ens  ***
427      !!
428      !! ** Purpose :   Compute the now total vorticity trend and add it to
429      !!      the general trend of the momentum equation.
430      !!
431      !! ** Method  :   Trend evaluated using now fields (centered in time)
432      !!      and the Sadourny (1975) flux FORM formulation : conserves the
433      !!      potential enstrophy of a horizontally non-divergent flow. the
434      !!      trend of the vorticity term is given by:
435      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
436      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
437      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
438      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
439      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
440      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
441      !!      Add this trend to the general momentum trend (ua,va):
442      !!          (ua,va) = (ua,va) + ( voru , vorv )
443      !!
444      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
445      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
446      !!               and planetary vorticity trends) ('key_trddyn')
447      !!
448      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
449      !!----------------------------------------------------------------------
450      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
451      USE wrk_nemo, ONLY:   zwx => wrk_2d_4, zwy => wrk_2d_5, zwz => wrk_2d_6    ! 2D workspace
452      !
453      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
454      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
455         !                                                        ! =nrvm (relative vorticity or metric)
456!FTRANS pua :I :I :z
457!FTRANS pva :I :I :z
458!! DCSE_NEMO: work around a deficiency in ftrans
459!     REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
460!     REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
461      REAL(wp), INTENT(inout) ::   pua(jpi,jpj,jpkorig)    ! total u-trend
462      REAL(wp), INTENT(inout) ::   pva(jpi,jpj,jpkorig)    ! total v-trend
463      !
464      INTEGER  ::   ji, jj, jk           ! dummy loop indices
465      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
466      !!----------------------------------------------------------------------
467     
468      IF( wrk_in_use(2, 4,5,6) ) THEN
469         CALL ctl_stop('dyn:vor_ens: requested workspace arrays unavailable')   ;   RETURN
470      END IF
471
472      IF( kt == nit000 ) THEN
473         IF(lwp) WRITE(numout,*)
474         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
475         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
476      ENDIF
477
478      zfact1 = 0.5 * 0.25      ! Local constant initialization
479
480!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
481      !                                                ! ===============
482      DO jk = 1, jpkm1                                 ! Horizontal slab
483         !                                             ! ===============
484         !
485         ! Potential vorticity and horizontal fluxes
486         ! -----------------------------------------
487         SELECT CASE( kvor )      ! vorticity considered
488         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
489         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
490         CASE ( 3 )                                                ! metric term
491            DO jj = 1, jpjm1
492               DO ji = 1, fs_jpim1   ! vector opt.
493                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
494                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
495                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
496               END DO
497            END DO
498         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
499         CASE ( 5 )                                                ! total (coriolis + metric)
500            DO jj = 1, jpjm1
501               DO ji = 1, fs_jpim1   ! vector opt.
502                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
503                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
504                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
505                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
506                       &       )
507               END DO
508            END DO
509         END SELECT
510         !
511         IF( ln_sco ) THEN
512            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
513               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
514                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
515                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
516                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
517               END DO
518            END DO
519         ELSE
520            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
521               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
522                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
523                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
524               END DO
525            END DO
526         ENDIF
527         !
528         ! Compute and add the vorticity term trend
529         ! ----------------------------------------
530         DO jj = 2, jpjm1
531            DO ji = fs_2, fs_jpim1   ! vector opt.
532               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
533                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
534               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
535                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
536               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
537               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
538            END DO 
539         END DO 
540         !                                             ! ===============
541      END DO                                           !   End of slab
542      !                                                ! ===============
543      IF( wrk_not_released(2, 4,5,6) )   CALL ctl_stop('dyn:vor_ens: failed to release workspace arrays')
544      !
545   END SUBROUTINE vor_ens
546
547
548   SUBROUTINE vor_een( kt, kvor, pua, pva )
549      !!----------------------------------------------------------------------
550      !!                ***  ROUTINE vor_een  ***
551      !!
552      !! ** Purpose :   Compute the now total vorticity trend and add it to
553      !!      the general trend of the momentum equation.
554      !!
555      !! ** Method  :   Trend evaluated using now fields (centered in time)
556      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
557      !!      both the horizontal kinetic energy and the potential enstrophy
558      !!      when horizontal divergence is zero (see the NEMO documentation)
559      !!      Add this trend to the general momentum trend (ua,va).
560      !!
561      !! ** Action : - Update (ua,va) with the now vorticity term trend
562      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
563      !!               and planetary vorticity trends) ('key_trddyn')
564      !!
565      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
566      !!----------------------------------------------------------------------
567      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
568#if defined key_z_first
569!FTRANS zwx :I :I :z
570!FTRANS zwy :I :I :z
571!FTRANS zwz :I :I :z
572!FTRANS ztnw :I :I :z
573!FTRANS ztne :I :I :z
574!FTRANS ztsw :I :I :z
575!FTRANS ztse :I :I :z
576      USE wrk_nemo, ONLY:   zwx  => wrk_3d_8 , zwy  => wrk_3d_2 ,  zwz => wrk_3d_3     ! 2D workspace
577      USE wrk_nemo, ONLY:   ztnw => wrk_3d_4 , ztne => wrk_3d_5 
578      USE wrk_nemo, ONLY:   ztsw => wrk_3d_6 , ztse => wrk_3d_7
579#else
580      USE wrk_nemo, ONLY:   zwx  => wrk_2d_1 , zwy  => wrk_2d_2 ,  zwz => wrk_2d_3     ! 2D workspace
581      USE wrk_nemo, ONLY:   ztnw => wrk_2d_4 , ztne => wrk_2d_5 
582      USE wrk_nemo, ONLY:   ztsw => wrk_2d_6 , ztse => wrk_2d_7
583#endif
584#if defined key_vvl
585!FTRANS ze3f :I :I :z
586      USE wrk_nemo, ONLY:   ze3f => wrk_3d_1                                           ! 3D workspace (lk_vvl=T)
587#endif
588      !
589      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
590      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
591      !                                                           ! =nrvm (relative vorticity or metric)
592!FTRANS pua :I :I :z
593!FTRANS pva :I :I :z
594!! DCSE_NEMO: work around a deficiency in ftrans
595      REAL(wp), INTENT(inout) ::   pua(jpi,jpj,jpkorig)   ! total u-trend
596      REAL(wp), INTENT(inout) ::   pva(jpi,jpj,jpkorig)   ! total v-trend
597      !!
598      INTEGER  ::   ji, jj, jk         ! dummy loop indices
599      INTEGER  ::   ierr               ! local integer
600      REAL(wp) ::   zfac12, zua, zva   ! local scalars
601#if ! defined key_vvl
602!FTRANS ze3f :I :I :z
603      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE ::   ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
604#endif
605      !!----------------------------------------------------------------------
606
607#if defined key_z_first
608      IF( wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN
609#else
610      IF( wrk_in_use(2, 1,2,3,4,5,6,7) .OR. wrk_in_use(3, 1) ) THEN
611#endif
612         CALL ctl_stop('dyn:vor_een: requested workspace arrays unavailable')   ;   RETURN
613      ENDIF
614
615      IF( kt == nit000 ) THEN
616         IF(lwp) WRITE(numout,*)
617         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
618         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
619         IF( .NOT.lk_vvl ) THEN
620            ALLOCATE( ze3f(jpi,jpj,jpkorig) , STAT=ierr )
621            IF( lk_mpp    )   CALL mpp_sum ( ierr )
622            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
623         ENDIF
624      ENDIF
625
626      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
627#if defined key_z_first
628         DO jj = 1, jpjm1
629            DO ji = 1, jpim1
630               DO jk = 1, mbkmax(ji,jj) ! jpk
631#else
632         DO jk = 1, jpk
633            DO jj = 1, jpjm1
634               DO ji = 1, jpim1
635#endif
636                  ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
637                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25
638                  IF( ze3f(ji,jj,jk) /= 0._wp )   ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
639               END DO
640            END DO
641         END DO
642         CALL lbc_lnk( ze3f, 'F', 1. )
643      ENDIF
644
645      zfac12 = 1._wp / 12._wp    ! Local constant initialization
646
647#if defined key_z_first
648         
649      ! Potential vorticity and horizontal fluxes
650      ! -----------------------------------------
651      SELECT CASE( kvor )      ! vorticity considered
652      CASE ( 1 )                                                ! planetary vorticity (Coriolis)
653         DO jj = 1, jpj, 1
654            DO ji = 1, jpi, 1
655               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
656                  zwz(ji,jj,jk) = ff(ji,jj)      * ze3f(ji,jj,jk)
657               END DO
658            END DO
659         END DO
660      CASE ( 2 )                                                ! relative  vorticity
661         DO jj = 1, jpj, 1
662            DO ji = 1, jpi, 1
663               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
664                  zwz(ji,jj,jk) = rotn(ji,jj,jk) * ze3f(ji,jj,jk)
665               END DO
666            END DO
667         END DO
668      CASE ( 3 )                                                ! metric term
669         DO jj = 1, jpjm1
670            DO ji = 1, jpim1 
671               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
672                  zwz(ji,jj,jk) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) ) &
673                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
674                       &     * 0.5_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
675               END DO
676            END DO
677         END DO
678         CALL lbc_lnk( zwz, 'F', 1. )
679      CASE ( 4 )                                                ! total (relative + planetary vorticity)
680         DO jj = 1, jpj, 1
681            DO ji = 1, jpi, 1
682               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
683                  zwz(ji,jj,jk) = ( rotn(ji,jj,jk) + ff(ji,jj) ) * ze3f(ji,jj,jk)
684               END DO
685            END DO
686         END DO
687      CASE ( 5 )                                                ! total (coriolis + metric)
688         DO jj = 1, jpjm1
689            DO ji = 1, fs_jpim1   ! vector opt.
690               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
691                  zwz(ji,jj,jk) = ( ff (ji,jj)                                                                 &
692                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
693                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
694                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
695                       &       ) * ze3f(ji,jj,jk)
696               END DO
697            END DO
698         END DO
699         CALL lbc_lnk( zwz, 'F', 1. )
700      END SELECT
701
702      DO jj = 1, jpj, 1
703         DO ji = 1, jpi, 1
704            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
705               zwx(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
706               zwy(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
707            END DO
708         END DO
709      END DO
710
711      ! Compute and add the vorticity term trend
712      ! ----------------------------------------
713      DO jj = 1,jpj,1
714         DO jk = 1,mbkmax(1,jj) ! jpk,1
715            ztne(1,jj,jk) = 0.0_wp ; ztnw(1,jj,jk) = 0.0_wp 
716            ztse(1,jj,jk) = 0.0_wp ; ztsw(1,jj,jk) = 0.0_wp
717         END DO
718      END DO
719
720      DO jj = 2, jpj
721         DO ji = 2, jpi
722            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
723               ztne(ji,jj,jk) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
724               ztnw(ji,jj,jk) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
725               ztse(ji,jj,jk) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
726               ztsw(ji,jj,jk) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
727            END DO
728         END DO
729      END DO
730      DO jj = 2, jpjm1
731         DO ji = fs_2, jpim1
732            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
733               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ,jk) * zwy(ji  ,jj  ,jk) + &
734                    ztnw(ji+1,jj,jk) * zwy(ji+1,jj  ,jk)   &
735                    &     + ztse(ji,jj  ,jk) * zwy(ji  ,jj-1,jk) + ztsw(ji+1,jj,jk) * zwy(ji+1,jj-1,jk) )
736               
737               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1,jk) * zwx(ji-1,jj+1,jk) + &
738                    ztse(ji,jj+1,jk) * zwx(ji  ,jj+1,jk)   &
739                    &     + ztnw(ji,jj  ,jk) * zwx(ji-1,jj  ,jk) + ztne(ji,jj  ,jk) * zwx(ji  ,jj  ,jk) )
740
741               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
742               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
743            END DO
744         END DO
745      END DO
746
747#else
748!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
749      !                                                ! ===============
750      DO jk = 1, jpkm1                                 ! Horizontal slab
751         !                                             ! ===============
752         
753         ! Potential vorticity and horizontal fluxes
754         ! -----------------------------------------
755         SELECT CASE( kvor )      ! vorticity considered
756         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
757            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
758         CASE ( 2 )                                                ! relative  vorticity
759            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
760         CASE ( 3 )                                                ! metric term
761            DO jj = 1, jpjm1
762               DO ji = 1, fs_jpim1   ! vector opt.
763                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
764                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
765                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
766               END DO
767            END DO
768            CALL lbc_lnk( zwz, 'F', 1. )
769        CASE ( 4 )                                                ! total (relative + planetary vorticity)
770            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
771         CASE ( 5 )                                                ! total (coriolis + metric)
772            DO jj = 1, jpjm1
773               DO ji = 1, fs_jpim1   ! vector opt.
774                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
775                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
776                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
777                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
778                       &       ) * ze3f(ji,jj,jk)
779               END DO
780            END DO
781            CALL lbc_lnk( zwz, 'F', 1. )
782         END SELECT
783
784         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
785         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
786
787         ! Compute and add the vorticity term trend
788         ! ----------------------------------------
789         jj = 2
790         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
791         DO ji = 2, jpi   
792               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
793               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
794               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
795               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
796         END DO
797         DO jj = 3, jpj
798            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
799               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
800               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
801               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
802               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
803            END DO
804         END DO
805         DO jj = 2, jpjm1
806            DO ji = fs_2, fs_jpim1   ! vector opt.
807               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
808                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
809               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
810                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
811               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
812               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
813            END DO 
814         END DO 
815         !                                             ! ===============
816      END DO                                           !   End of slab
817      !                                                ! ===============
818#endif
819
820#if defined key_z_first
821      IF( wrk_not_released(3, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:vor_een: failed to release workspace arrays')
822#else
823      IF( wrk_not_released(2, 1,2,3,4,5,6,7) .OR.   &
824          wrk_not_released(3, 1)             )   CALL ctl_stop('dyn:vor_een: failed to release workspace arrays')
825#endif
826      !
827   END SUBROUTINE vor_een
828
829
830   SUBROUTINE dyn_vor_init
831      !!---------------------------------------------------------------------
832      !!                  ***  ROUTINE dyn_vor_init  ***
833      !!
834      !! ** Purpose :   Control the consistency between cpp options for
835      !!              tracer advection schemes
836      !!----------------------------------------------------------------------
837      INTEGER ::   ioptio          ! local integer
838      !!
839      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
840      !!----------------------------------------------------------------------
841
842      REWIND ( numnam )               ! Read Namelist namdyn_vor : Vorticity scheme options
843      READ   ( numnam, namdyn_vor )
844
845      IF(lwp) THEN                    ! Namelist print
846         WRITE(numout,*)
847         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
848         WRITE(numout,*) '~~~~~~~~~~~~'
849         WRITE(numout,*) '        Namelist namdyn_vor : oice of the vorticity term scheme'
850         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
851         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
852         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
853         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
854      ENDIF
855
856      ioptio = 0                     ! Control of vorticity scheme options
857      IF( ln_dynvor_ene )   ioptio = ioptio + 1
858      IF( ln_dynvor_ens )   ioptio = ioptio + 1
859      IF( ln_dynvor_mix )   ioptio = ioptio + 1
860      IF( ln_dynvor_een )   ioptio = ioptio + 1
861      IF( lk_esopa      )   ioptio =          1
862
863      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
864
865      !                              ! Set nvor (type of scheme for vorticity)
866      IF( ln_dynvor_ene )   nvor =  0
867      IF( ln_dynvor_ens )   nvor =  1
868      IF( ln_dynvor_mix )   nvor =  2
869      IF( ln_dynvor_een )   nvor =  3
870      IF( lk_esopa      )   nvor = -1
871     
872      !                              ! Set ncor, nrvm, ntot (type of vorticity)
873      IF(lwp) WRITE(numout,*)
874      ncor = 1
875      IF( ln_dynadv_vec ) THEN     
876         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
877         nrvm = 2
878         ntot = 4
879      ELSE                       
880         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
881         nrvm = 3
882         ntot = 5
883      ENDIF
884     
885      IF(lwp) THEN                   ! Print the choice
886         WRITE(numout,*)
887         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
888         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
889         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
890         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
891         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
892      ENDIF
893      !
894   END SUBROUTINE dyn_vor_init
895
896   !!==============================================================================
897END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.