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/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7905

Last change on this file since 7905 was 7905, checked in by jcastill, 7 years ago

Series of small bug fixes and stetic changes:

-Fix possible bug in the calculation of Stokes-Coriolis
-Move all the wave control variables to namelist namsbc_wave
-Use one namelist variable instead of two to set Stokes drift velocity coupling
-Cap the values of the Craig and Banner constant as calculated from wave input fields to take into account small values of the friction velocity
-Add new Phillips parametrization for Stokes drift vertical velocity, using the inverse depth scale as in Breivik 2015, instead of the peak wave number as calculated from wave input fields
-Better control of the wave fields that are read from file depending on the wave parameters
File size: 42.8 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   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dyn_vor      : Update the momentum trend with the vorticity trend
23   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
24   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
25   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
26   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
27   !!   dyn_vor_init : set and control of the different vorticity option
28   !!----------------------------------------------------------------------
29   USE oce            ! ocean dynamics and tracers
30   USE dom_oce        ! ocean space and time domain
31   USE dommsk         ! ocean mask
32   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
33   USE trd_oce        ! trends: ocean variables
34   USE trddyn         ! trend manager: dynamics
35   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
36   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
37   !
38   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
39   USE prtctl         ! Print control
40   USE in_out_manager ! I/O manager
41   USE lib_mpp        ! MPP library
42   USE wrk_nemo       ! Memory Allocation
43   USE timing         ! Timing
44
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      !!----------------------------------------------------------------------
88      !
89      IF( nn_timing == 1 )  CALL timing_start('dyn_vor')
90      !
91      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
92      !
93      !                                          ! vorticity term
94      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
95      !
96      CASE ( -1 )                                      ! esopa: test all possibility with control print
97         CALL vor_ene( kt, ntot, un, vn, ua, va )
98         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
99            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
100         CALL vor_ens( kt, ntot, un, vn, ua, va )
101         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
102            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
103         CALL vor_mix( kt )
104         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
105            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
106         CALL vor_een( kt, ntot, un, vn, ua, va )
107         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
108            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
109         !
110      CASE ( 0 )                                       ! energy conserving scheme
111         IF( l_trddyn )   THEN
112            ztrdu(:,:,:) = ua(:,:,:)
113            ztrdv(:,:,:) = va(:,:,:)
114            CALL vor_ene( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
115            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
116            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
117            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
118            ztrdu(:,:,:) = ua(:,:,:)
119            ztrdv(:,:,:) = va(:,:,:)
120            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
121            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
122            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
123            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
124         ELSE
125                             CALL vor_ene( kt, ntot, un, vn, ua, va )    ! total vorticity trend
126            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
127         ENDIF
128         !
129      CASE ( 1 )                                       ! enstrophy conserving scheme
130         IF( l_trddyn )   THEN   
131            ztrdu(:,:,:) = ua(:,:,:)
132            ztrdv(:,:,:) = va(:,:,:)
133            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
134            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
135            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
136            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
137            ztrdu(:,:,:) = ua(:,:,:)
138            ztrdv(:,:,:) = va(:,:,:)
139            CALL vor_ens( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
140            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
141            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
142            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
143         ELSE
144                             CALL vor_ens( kt, ntot, un, vn, ua, va )    ! total vorticity
145            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
146         ENDIF
147         !
148      CASE ( 2 )                                       ! mixed ene-ens scheme
149         IF( l_trddyn )   THEN
150            ztrdu(:,:,:) = ua(:,:,:)
151            ztrdv(:,:,:) = va(:,:,:)
152            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend (ens)
153            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
154            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
155            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
156            ztrdu(:,:,:) = ua(:,:,:)
157            ztrdv(:,:,:) = va(:,:,:)
158            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend (ene)
159            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
160            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
161            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
162         ELSE
163                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
164                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
165            IF( ln_stcor )   CALL vor_ens( kt, nrvm, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
166            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
167         ENDIF
168         !
169      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
170         IF( l_trddyn )   THEN
171            ztrdu(:,:,:) = ua(:,:,:)
172            ztrdv(:,:,:) = va(:,:,:)
173            CALL vor_een( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend
174            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
175            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
176            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
177            ztrdu(:,:,:) = ua(:,:,:)
178            ztrdv(:,:,:) = va(:,:,:)
179            CALL vor_een( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend
180            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
181            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
182            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
183         ELSE
184                             CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity
185            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
186         ENDIF
187         !
188      END SELECT
189      !
190      !                       ! print sum trends (used for debugging)
191      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
192         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
193      !
194      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
195      !
196      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor')
197      !
198   END SUBROUTINE dyn_vor
199
200
201   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
202      !!----------------------------------------------------------------------
203      !!                  ***  ROUTINE vor_ene  ***
204      !!
205      !! ** Purpose :   Compute the now total vorticity trend and add it to
206      !!      the general trend of the momentum equation.
207      !!
208      !! ** Method  :   Trend evaluated using now fields (centered in time)
209      !!      and the Sadourny (1975) flux form formulation : conserves the
210      !!      horizontal kinetic energy.
211      !!      The trend of the vorticity term is given by:
212      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
213      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
214      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
215      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
216      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
217      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
218      !!      Add this trend to the general momentum trend (ua,va):
219      !!          (ua,va) = (ua,va) + ( voru , vorv )
220      !!
221      !! ** Action : - Update (ua,va) with the now vorticity term trend
222      !!
223      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
224      !!----------------------------------------------------------------------
225      !
226      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
227      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
228      !                                                                ! =nrvm (relative vorticity or metric)
229      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend
230      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend
231      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
232      !
233      INTEGER  ::   ji, jj, jk   ! dummy loop indices
234      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
235      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
236      !!----------------------------------------------------------------------
237      !
238      IF( nn_timing == 1 )  CALL timing_start('vor_ene')
239      !
240      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
241      !
242      IF( kt == nit000 ) THEN
243         IF(lwp) WRITE(numout,*)
244         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
245         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
246      ENDIF
247
248      zfact2 = 0.5 * 0.5      ! Local constant initialization
249
250!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
251      !                                                ! ===============
252      DO jk = 1, jpkm1                                 ! Horizontal slab
253         !                                             ! ===============
254         !
255         ! Potential vorticity and horizontal fluxes
256         ! -----------------------------------------
257         SELECT CASE( kvor )      ! vorticity considered
258         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
259         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
260         CASE ( 3 )                                                ! metric term
261            DO jj = 1, jpjm1
262               DO ji = 1, fs_jpim1   ! vector opt.
263                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
264                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
265                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
266               END DO
267            END DO
268         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
269         CASE ( 5 )                                                ! total (coriolis + metric)
270            DO jj = 1, jpjm1
271               DO ji = 1, fs_jpim1   ! vector opt.
272                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
273                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
274                       &           - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
275                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
276                       &       )
277               END DO
278            END DO
279         END SELECT
280
281         IF( ln_sco ) THEN
282            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
283            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk)
284            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk)
285         ELSE
286            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
287            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
288         ENDIF
289
290         ! Compute and add the vorticity term trend
291         ! ----------------------------------------
292         DO jj = 2, jpjm1
293            DO ji = fs_2, fs_jpim1   ! vector opt.
294               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
295               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
296               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
297               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
298               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
299               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
300            END DO 
301         END DO 
302         !                                             ! ===============
303      END DO                                           !   End of slab
304      !                                                ! ===============
305      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
306      !
307      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
308      !
309   END SUBROUTINE vor_ene
310
311
312   SUBROUTINE vor_mix( kt )
313      !!----------------------------------------------------------------------
314      !!                 ***  ROUTINE vor_mix  ***
315      !!
316      !! ** Purpose :   Compute the now total vorticity trend and add it to
317      !!      the general trend of the momentum equation.
318      !!
319      !! ** Method  :   Trend evaluated using now fields (centered in time)
320      !!      Mixte formulation : conserves the potential enstrophy of a hori-
321      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
322      !!      ticity term and the horizontal kinetic energy for (f x uh), the
323      !!      coriolis term. the now trend of the vorticity term is given by:
324      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
325      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
326      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
327      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
328      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
329      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
330      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
331      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
332      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
333      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
334      !!      Add this now trend to the general momentum trend (ua,va):
335      !!          (ua,va) = (ua,va) + ( voru , vorv )
336      !!
337      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
338      !!
339      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
340      !!----------------------------------------------------------------------
341      !
342      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
343      !
344      INTEGER  ::   ji, jj, jk   ! dummy loop indices
345      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
346      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
347      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
348      !!----------------------------------------------------------------------
349      !
350      IF( nn_timing == 1 )  CALL timing_start('vor_mix')
351      !
352      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) 
353      !
354      IF( kt == nit000 ) THEN
355         IF(lwp) WRITE(numout,*)
356         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
357         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
358      ENDIF
359
360      zfact1 = 0.5 * 0.25      ! Local constant initialization
361      zfact2 = 0.5 * 0.5
362
363!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
364      !                                                ! ===============
365      DO jk = 1, jpkm1                                 ! Horizontal slab
366         !                                             ! ===============
367         !
368         ! Relative and planetary potential vorticity and horizontal fluxes
369         ! ----------------------------------------------------------------
370         IF( ln_sco ) THEN       
371            IF( ln_dynadv_vec ) THEN
372               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
373            ELSE                       
374               DO jj = 1, jpjm1
375                  DO ji = 1, fs_jpim1   ! vector opt.
376                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
377                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
378                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
379                  END DO
380               END DO
381            ENDIF
382            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
383            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
384            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
385         ELSE
386            IF( ln_dynadv_vec ) THEN
387               zww(:,:) = rotn(:,:,jk)
388            ELSE                       
389               DO jj = 1, jpjm1
390                  DO ji = 1, fs_jpim1   ! vector opt.
391                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
392                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
393                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
394                  END DO
395               END DO
396            ENDIF
397            zwz(:,:) = ff (:,:)
398            zwx(:,:) = e2u(:,:) * un(:,:,jk)
399            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
400         ENDIF
401
402         ! Compute and add the vorticity term trend
403         ! ----------------------------------------
404         DO jj = 2, jpjm1
405            DO ji = fs_2, fs_jpim1   ! vector opt.
406               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
407               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
408               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
409               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
410               ! enstrophy conserving formulation for relative vorticity term
411               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
412               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
413               ! energy conserving formulation for planetary vorticity term
414               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
415               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
416               ! mixed vorticity trend added to the momentum trends
417               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
418               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
419            END DO 
420         END DO 
421         !                                             ! ===============
422      END DO                                           !   End of slab
423      !                                                ! ===============
424      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) 
425      !
426      IF( nn_timing == 1 )  CALL timing_stop('vor_mix')
427      !
428   END SUBROUTINE vor_mix
429
430
431   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
432      !!----------------------------------------------------------------------
433      !!                ***  ROUTINE vor_ens  ***
434      !!
435      !! ** Purpose :   Compute the now total vorticity trend and add it to
436      !!      the general trend of the momentum equation.
437      !!
438      !! ** Method  :   Trend evaluated using now fields (centered in time)
439      !!      and the Sadourny (1975) flux FORM formulation : conserves the
440      !!      potential enstrophy of a horizontally non-divergent flow. the
441      !!      trend of the vorticity term is given by:
442      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
443      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
444      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
445      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
446      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
447      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
448      !!      Add this trend to the general momentum trend (ua,va):
449      !!          (ua,va) = (ua,va) + ( voru , vorv )
450      !!
451      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
452      !!
453      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
454      !!----------------------------------------------------------------------
455      !
456      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
457      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
458         !                                                             ! =nrvm (relative vorticity or metric)
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(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
462      !
463      INTEGER  ::   ji, jj, jk           ! dummy loop indices
464      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
465      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
466      !!----------------------------------------------------------------------
467      !
468      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
469      !
470      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
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) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
494                       &         - ( pun(ji  ,jj+1,jk) + pun (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                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
504                       &           - ( pun(ji  ,jj+1,jk) + pun (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) * pun(ji,jj,jk)
516                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(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) * pun(ji,jj,jk)
523                  zwy(ji,jj) = e1v(ji,jj) * pvn(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      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
544      !
545      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
546      !
547   END SUBROUTINE vor_ens
548
549
550   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
551      !!----------------------------------------------------------------------
552      !!                ***  ROUTINE vor_een  ***
553      !!
554      !! ** Purpose :   Compute the now total vorticity trend and add it to
555      !!      the general trend of the momentum equation.
556      !!
557      !! ** Method  :   Trend evaluated using now fields (centered in time)
558      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
559      !!      both the horizontal kinetic energy and the potential enstrophy
560      !!      when horizontal divergence is zero (see the NEMO documentation)
561      !!      Add this trend to the general momentum trend (ua,va).
562      !!
563      !! ** Action : - Update (ua,va) with the now vorticity term trend
564      !!
565      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
566      !!----------------------------------------------------------------------
567      !
568      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
569      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
570      !                                                                ! =nrvm (relative vorticity or metric)
571      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend
572      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend
573      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
574      !!
575      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
576      INTEGER  ::   ierr                                          ! local integer
577      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
578      REAL(wp) ::   zmsk, ze3                                     ! local scalars
579      !                                                           !  3D workspace
580      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
581      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
582#if defined key_vvl
583      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
584#else
585      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
586#endif
587      !!----------------------------------------------------------------------
588      !
589      IF( nn_timing == 1 )  CALL timing_start('vor_een')
590      !
591      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
592      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
593#if defined key_vvl
594      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
595#endif
596      !
597      IF( kt == nit000 ) THEN
598         IF(lwp) WRITE(numout,*)
599         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
600         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
601#if ! defined key_vvl
602         IF( .NOT.ALLOCATED(ze3f) ) THEN
603            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
604            IF( lk_mpp    )   CALL mpp_sum ( ierr )
605            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
606         ENDIF
607         ze3f(:,:,:) = 0._wp
608#endif
609      ENDIF
610
611      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
612
613         IF( ln_dynvor_een_old ) THEN ! original formulation
614            DO jk = 1, jpk
615               DO jj = 1, jpjm1
616                  DO ji = 1, jpim1
617                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
618                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
619                     IF( ze3 /= 0._wp ) THEN  ;   ze3f(ji,jj,jk) = 4.0_wp / ze3
620                     ELSE                     ;   ze3f(ji,jj,jk) = 0._wp
621                     ENDIF
622                  END DO
623               END DO
624            END DO
625         ELSE ! new formulation from NEMO 3.6
626            DO jk = 1, jpk
627               DO jj = 1, jpjm1
628                  DO ji = 1, jpim1
629                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
630                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
631                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
632                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
633                     IF( ze3 /= 0._wp ) THEN  ;  ze3f(ji,jj,jk) = zmsk / ze3
634                     ELSE                     ;  ze3f(ji,jj,jk) = 0._wp
635                     ENDIF
636                  END DO
637               END DO
638            END DO
639         ENDIF
640
641         CALL lbc_lnk( ze3f, 'F', 1. )
642      ENDIF
643
644      zfac12 = 1._wp / 12._wp    ! Local constant initialization
645
646     
647!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
648      !                                                ! ===============
649      DO jk = 1, jpkm1                                 ! Horizontal slab
650         !                                             ! ===============
651         
652         ! Potential vorticity and horizontal fluxes
653         ! -----------------------------------------
654         SELECT CASE( kvor )      ! vorticity considered
655         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
656            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
657         CASE ( 2 )                                                ! relative  vorticity
658            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
659         CASE ( 3 )                                                ! metric term
660            DO jj = 1, jpjm1
661               DO ji = 1, fs_jpim1   ! vector opt.
662                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
663                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
664                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
665               END DO
666            END DO
667            CALL lbc_lnk( zwz, 'F', 1. )
668        CASE ( 4 )                                                ! total (relative + planetary vorticity)
669            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
670         CASE ( 5 )                                                ! total (coriolis + metric)
671            DO jj = 1, jpjm1
672               DO ji = 1, fs_jpim1   ! vector opt.
673                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
674                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
675                       &           - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
676                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
677                       &       ) * ze3f(ji,jj,jk)
678               END DO
679            END DO
680            CALL lbc_lnk( zwz, 'F', 1. )
681         END SELECT
682
683         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk)
684         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk)
685
686         ! Compute and add the vorticity term trend
687         ! ----------------------------------------
688         jj = 2
689         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
690         DO ji = 2, jpi   
691               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
692               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
693               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
694               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
695         END DO
696         DO jj = 3, jpj
697            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
698               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
699               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
700               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
701               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
702            END DO
703         END DO
704         DO jj = 2, jpjm1
705            DO ji = fs_2, fs_jpim1   ! vector opt.
706               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
707                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
708               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
709                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
710               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
711               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
712            END DO 
713         END DO 
714         !                                             ! ===============
715      END DO                                           !   End of slab
716      !                                                ! ===============
717      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
718      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
719#if defined key_vvl
720      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
721#endif
722      !
723      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
724      !
725   END SUBROUTINE vor_een
726
727
728   SUBROUTINE dyn_vor_init
729      !!---------------------------------------------------------------------
730      !!                  ***  ROUTINE dyn_vor_init  ***
731      !!
732      !! ** Purpose :   Control the consistency between cpp options for
733      !!              tracer advection schemes
734      !!----------------------------------------------------------------------
735      INTEGER ::   ioptio          ! local integer
736      INTEGER ::   ji, jj, jk      ! dummy loop indices
737      INTEGER ::   ios             ! Local integer output status for namelist read
738      !!
739      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
740      !!----------------------------------------------------------------------
741
742      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
743      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
744901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
745
746      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
747      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
748902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
749      IF(lwm) WRITE ( numond, namdyn_vor )
750
751      IF(lwp) THEN                    ! Namelist print
752         WRITE(numout,*)
753         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
754         WRITE(numout,*) '~~~~~~~~~~~~'
755         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
756         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
757         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
758         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
759         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
760         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
761      ENDIF
762
763      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
764      ! at angles with three ocean points and one land point
765      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
766         DO jk = 1, jpk
767            DO jj = 2, jpjm1
768               DO ji = 2, jpim1
769                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
770                      fmask(ji,jj,jk) = 1._wp
771               END DO
772            END DO
773         END DO
774          !
775          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
776          !
777      ENDIF
778
779      ioptio = 0                     ! Control of vorticity scheme options
780      IF( ln_dynvor_ene )   ioptio = ioptio + 1
781      IF( ln_dynvor_ens )   ioptio = ioptio + 1
782      IF( ln_dynvor_mix )   ioptio = ioptio + 1
783      IF( ln_dynvor_een )   ioptio = ioptio + 1
784      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
785      IF( lk_esopa      )   ioptio =          1
786
787      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
788
789      !                              ! Set nvor (type of scheme for vorticity)
790      IF( ln_dynvor_ene )   nvor =  0
791      IF( ln_dynvor_ens )   nvor =  1
792      IF( ln_dynvor_mix )   nvor =  2
793      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
794      IF( lk_esopa      )   nvor = -1
795     
796      !                              ! Set ncor, nrvm, ntot (type of vorticity)
797      IF(lwp) WRITE(numout,*)
798      ncor = 1
799      IF( ln_dynadv_vec ) THEN     
800         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
801         nrvm = 2
802         ntot = 4
803      ELSE                       
804         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
805         nrvm = 3
806         ntot = 5
807      ENDIF
808     
809      IF(lwp) THEN                   ! Print the choice
810         WRITE(numout,*)
811         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
812         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
813         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
814         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
815         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
816      ENDIF
817      !
818   END SUBROUTINE dyn_vor_init
819
820   !!==============================================================================
821END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.