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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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