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/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

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