source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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