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 @ 2676

Last change on this file since 2676 was 2636, checked in by gm, 13 years ago

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