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/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_vvl/MY_SRC – NEMO

source: branches/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_vvl/MY_SRC/dynvor.F90 @ 4648

Last change on this file since 4648 was 4648, checked in by flavoni, 10 years ago

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