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

source: branches/2011/dev_LOCEAN_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 2977

Last change on this file since 2977 was 2977, checked in by cetlod, 13 years ago

Add in branch 2011/dev_LOCEAN_2011 changes from 2011/dev_r2787_PISCES_improvment, 2011/dev_r2787_LOCEAN_offline_fldread and 2011/dev_r2787_LOCEAN3_TRA_TRP branches, see ticket #877
  • 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   !!----------------------------------------------------------------------
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:   tsa            ! tsa used as 2 3D workspace
74      !!
75      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
76      !
77      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
78      !!----------------------------------------------------------------------
79      !
80      IF( l_trddyn )   THEN
81         ztrdu => tsa(:,:,:,1) 
82         ztrdv => tsa(:,:,:,2) 
83      END IF
84      !
85      !                                          ! vorticity term
86      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
87      !
88      CASE ( -1 )                                      ! esopa: test all possibility with control print
89         CALL vor_ene( kt, ntot, ua, va )
90         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
91            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
92         CALL vor_ens( kt, ntot, ua, va )
93         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
94            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
95         CALL vor_mix( kt )
96         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
97            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
98         CALL vor_een( kt, ntot, ua, va )
99         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
100            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
101         !
102      CASE ( 0 )                                       ! energy conserving scheme
103         IF( l_trddyn )   THEN
104            ztrdu(:,:,:) = ua(:,:,:)
105            ztrdv(:,:,:) = va(:,:,:)
106            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
107            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
108            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
109            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
110            ztrdu(:,:,:) = ua(:,:,:)
111            ztrdv(:,:,:) = va(:,:,:)
112            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
113            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
114            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
115            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
116            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
117         ELSE
118            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
119         ENDIF
120         !
121      CASE ( 1 )                                       ! enstrophy conserving scheme
122         IF( l_trddyn )   THEN   
123            ztrdu(:,:,:) = ua(:,:,:)
124            ztrdv(:,:,:) = va(:,:,:)
125            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
126            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
127            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
128            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
129            ztrdu(:,:,:) = ua(:,:,:)
130            ztrdv(:,:,:) = va(:,:,:)
131            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
132            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
133            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
134            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
135            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', 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_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', 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_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
154            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
155         ELSE
156            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
157         ENDIF
158         !
159      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
160         IF( l_trddyn )   THEN
161            ztrdu(:,:,:) = ua(:,:,:)
162            ztrdv(:,:,:) = va(:,:,:)
163            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
164            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
165            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
166            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
167            ztrdu(:,:,:) = ua(:,:,:)
168            ztrdv(:,:,:) = va(:,:,:)
169            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
170            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
171            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
172            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
173            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
174         ELSE
175            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
176         ENDIF
177         !
178      END SELECT
179      !
180      !                       ! print sum trends (used for debugging)
181      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
182         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
183      !
184   END SUBROUTINE dyn_vor
185
186
187   SUBROUTINE vor_ene( kt, kvor, pua, pva )
188      !!----------------------------------------------------------------------
189      !!                  ***  ROUTINE vor_ene  ***
190      !!
191      !! ** Purpose :   Compute the now total vorticity trend and add it to
192      !!      the general trend of the momentum equation.
193      !!
194      !! ** Method  :   Trend evaluated using now fields (centered in time)
195      !!      and the Sadourny (1975) flux form formulation : conserves the
196      !!      horizontal kinetic energy.
197      !!      The trend of the vorticity term is given by:
198      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
199      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
200      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
201      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
202      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
203      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
204      !!      Add this trend to the general momentum trend (ua,va):
205      !!          (ua,va) = (ua,va) + ( voru , vorv )
206      !!
207      !! ** Action : - Update (ua,va) with the now vorticity term trend
208      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
209      !!               and planetary vorticity trends) ('key_trddyn')
210      !!
211      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
212      !!----------------------------------------------------------------------
213      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
214      USE wrk_nemo, ONLY:   zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3     ! 2D workspace
215      !
216      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
217      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
218      !                                                           ! =nrvm (relative vorticity or metric)
219      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
220      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
221      !
222      INTEGER  ::   ji, jj, jk   ! dummy loop indices
223      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
224      !!----------------------------------------------------------------------
225
226      IF( wrk_in_use(2, 1,2,3) ) THEN
227         CALL ctl_stop('dyn:vor_ene: requested workspace arrays unavailable')   ;   RETURN
228      ENDIF
229
230      IF( kt == nit000 ) THEN
231         IF(lwp) WRITE(numout,*)
232         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
233         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
234      ENDIF
235
236      zfact2 = 0.5 * 0.5      ! Local constant initialization
237
238!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
239      !                                                ! ===============
240      DO jk = 1, jpkm1                                 ! Horizontal slab
241         !                                             ! ===============
242         !
243         ! Potential vorticity and horizontal fluxes
244         ! -----------------------------------------
245         SELECT CASE( kvor )      ! vorticity considered
246         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
247         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
248         CASE ( 3 )                                                ! metric term
249            DO jj = 1, jpjm1
250               DO ji = 1, fs_jpim1   ! vector opt.
251                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
252                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
253                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
254               END DO
255            END DO
256         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
257         CASE ( 5 )                                                ! total (coriolis + metric)
258            DO jj = 1, jpjm1
259               DO ji = 1, fs_jpim1   ! vector opt.
260                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
261                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
262                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
263                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
264                       &       )
265               END DO
266            END DO
267         END SELECT
268
269         IF( ln_sco ) THEN
270            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
271            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
272            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
273         ELSE
274            zwx(:,:) = e2u(:,:) * un(:,:,jk)
275            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
276         ENDIF
277
278         ! Compute and add the vorticity term trend
279         ! ----------------------------------------
280         DO jj = 2, jpjm1
281            DO ji = fs_2, fs_jpim1   ! vector opt.
282               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
283               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
284               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
285               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
286               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
287               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
288            END DO 
289         END DO 
290         !                                             ! ===============
291      END DO                                           !   End of slab
292      !                                                ! ===============
293      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn:vor_ene: failed to release workspace arrays')
294      !
295   END SUBROUTINE vor_ene
296
297
298   SUBROUTINE vor_mix( kt )
299      !!----------------------------------------------------------------------
300      !!                 ***  ROUTINE vor_mix  ***
301      !!
302      !! ** Purpose :   Compute the now total vorticity trend and add it to
303      !!      the general trend of the momentum equation.
304      !!
305      !! ** Method  :   Trend evaluated using now fields (centered in time)
306      !!      Mixte formulation : conserves the potential enstrophy of a hori-
307      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
308      !!      ticity term and the horizontal kinetic energy for (f x uh), the
309      !!      coriolis term. the now trend of the vorticity term is given by:
310      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
311      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
312      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
313      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
314      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
315      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
316      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
317      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
318      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
319      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
320      !!      Add this now trend to the general momentum trend (ua,va):
321      !!          (ua,va) = (ua,va) + ( voru , vorv )
322      !!
323      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
324      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
325      !!               and planetary vorticity trends) ('key_trddyn')
326      !!
327      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
328      !!----------------------------------------------------------------------
329      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
330      USE wrk_nemo, ONLY:   zwx => wrk_2d_4 , zwy => wrk_2d_5 , zwz => wrk_2d_6 , zww => wrk_2d_7   ! 2D workspace
331      !
332      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
333      !
334      INTEGER  ::   ji, jj, jk   ! dummy loop indices
335      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
336      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
337      !!----------------------------------------------------------------------
338
339      IF( wrk_in_use(2, 4,5,6,7) ) THEN
340         CALL ctl_stop('dyn:vor_mix: requested workspace arrays unavailable')   ;   RETURN
341      ENDIF
342
343      IF( kt == nit000 ) THEN
344         IF(lwp) WRITE(numout,*)
345         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
346         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
347      ENDIF
348
349      zfact1 = 0.5 * 0.25      ! Local constant initialization
350      zfact2 = 0.5 * 0.5
351
352!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
353      !                                                ! ===============
354      DO jk = 1, jpkm1                                 ! Horizontal slab
355         !                                             ! ===============
356         !
357         ! Relative and planetary potential vorticity and horizontal fluxes
358         ! ----------------------------------------------------------------
359         IF( ln_sco ) THEN       
360            IF( ln_dynadv_vec ) THEN
361               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
362            ELSE                       
363               DO jj = 1, jpjm1
364                  DO ji = 1, fs_jpim1   ! vector opt.
365                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
366                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
367                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
368                  END DO
369               END DO
370            ENDIF
371            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
372            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
373            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
374         ELSE
375            IF( ln_dynadv_vec ) THEN
376               zww(:,:) = rotn(:,:,jk)
377            ELSE                       
378               DO jj = 1, jpjm1
379                  DO ji = 1, fs_jpim1   ! vector opt.
380                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
381                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
382                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
383                  END DO
384               END DO
385            ENDIF
386            zwz(:,:) = ff (:,:)
387            zwx(:,:) = e2u(:,:) * un(:,:,jk)
388            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
389         ENDIF
390
391         ! Compute and add the vorticity term trend
392         ! ----------------------------------------
393         DO jj = 2, jpjm1
394            DO ji = fs_2, fs_jpim1   ! vector opt.
395               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
396               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
397               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
398               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
399               ! enstrophy conserving formulation for relative vorticity term
400               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
401               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
402               ! energy conserving formulation for planetary vorticity term
403               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
404               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
405               ! mixed vorticity trend added to the momentum trends
406               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
407               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
408            END DO 
409         END DO 
410         !                                             ! ===============
411      END DO                                           !   End of slab
412      !                                                ! ===============
413      IF( wrk_not_released(2, 4,5,6,7) )   CALL ctl_stop('dyn:vor_mix: failed to release workspace arrays')
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      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
440      !!               and planetary vorticity trends) ('key_trddyn')
441      !!
442      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
443      !!----------------------------------------------------------------------
444      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
445      USE wrk_nemo, ONLY:   zwx => wrk_2d_4, zwy => wrk_2d_5, zwz => wrk_2d_6    ! 2D workspace
446      !
447      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
448      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
449         !                                                        ! =nrvm (relative vorticity or metric)
450      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
451      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
452      !
453      INTEGER  ::   ji, jj, jk           ! dummy loop indices
454      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
455      !!----------------------------------------------------------------------
456     
457      IF( wrk_in_use(2, 4,5,6) ) THEN
458         CALL ctl_stop('dyn:vor_ens: requested workspace arrays unavailable')   ;   RETURN
459      END IF
460
461      IF( kt == nit000 ) THEN
462         IF(lwp) WRITE(numout,*)
463         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
464         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
465      ENDIF
466
467      zfact1 = 0.5 * 0.25      ! Local constant initialization
468
469!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
470      !                                                ! ===============
471      DO jk = 1, jpkm1                                 ! Horizontal slab
472         !                                             ! ===============
473         !
474         ! Potential vorticity and horizontal fluxes
475         ! -----------------------------------------
476         SELECT CASE( kvor )      ! vorticity considered
477         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
478         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
479         CASE ( 3 )                                                ! metric term
480            DO jj = 1, jpjm1
481               DO ji = 1, fs_jpim1   ! vector opt.
482                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
483                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
484                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
485               END DO
486            END DO
487         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
488         CASE ( 5 )                                                ! total (coriolis + metric)
489            DO jj = 1, jpjm1
490               DO ji = 1, fs_jpim1   ! vector opt.
491                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
492                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
493                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
494                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
495                       &       )
496               END DO
497            END DO
498         END SELECT
499         !
500         IF( ln_sco ) THEN
501            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
502               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
503                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
504                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
505                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
506               END DO
507            END DO
508         ELSE
509            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
510               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
511                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
512                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
513               END DO
514            END DO
515         ENDIF
516         !
517         ! Compute and add the vorticity term trend
518         ! ----------------------------------------
519         DO jj = 2, jpjm1
520            DO ji = fs_2, fs_jpim1   ! vector opt.
521               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
522                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
523               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
524                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
525               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
526               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
527            END DO 
528         END DO 
529         !                                             ! ===============
530      END DO                                           !   End of slab
531      !                                                ! ===============
532      IF( wrk_not_released(2, 4,5,6) )   CALL ctl_stop('dyn:vor_ens: failed to release workspace arrays')
533      !
534   END SUBROUTINE vor_ens
535
536
537   SUBROUTINE vor_een( kt, kvor, pua, pva )
538      !!----------------------------------------------------------------------
539      !!                ***  ROUTINE vor_een  ***
540      !!
541      !! ** Purpose :   Compute the now total vorticity trend and add it to
542      !!      the general trend of the momentum equation.
543      !!
544      !! ** Method  :   Trend evaluated using now fields (centered in time)
545      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
546      !!      both the horizontal kinetic energy and the potential enstrophy
547      !!      when horizontal divergence is zero (see the NEMO documentation)
548      !!      Add this trend to the general momentum trend (ua,va).
549      !!
550      !! ** Action : - Update (ua,va) with the now vorticity term trend
551      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
552      !!               and planetary vorticity trends) ('key_trddyn')
553      !!
554      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
555      !!----------------------------------------------------------------------
556      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
557      USE wrk_nemo, ONLY:   zwx  => wrk_2d_1 , zwy  => wrk_2d_2 ,  zwz => wrk_2d_3     ! 2D workspace
558      USE wrk_nemo, ONLY:   ztnw => wrk_2d_4 , ztne => wrk_2d_5 
559      USE wrk_nemo, ONLY:   ztsw => wrk_2d_6 , ztse => wrk_2d_7
560#if defined key_vvl
561      USE wrk_nemo, ONLY:   ze3f => wrk_3d_1                                           ! 3D workspace (lk_vvl=T)
562#endif
563      !
564      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
565      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
566      !                                                           ! =nrvm (relative vorticity or metric)
567      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
568      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
569      !!
570      INTEGER  ::   ji, jj, jk         ! dummy loop indices
571      INTEGER  ::   ierr               ! local integer
572      REAL(wp) ::   zfac12, zua, zva   ! local scalars
573#if ! defined key_vvl
574      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE ::   ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
575#endif
576      !!----------------------------------------------------------------------
577
578      IF( wrk_in_use(2, 1,2,3,4,5,6,7) .OR. wrk_in_use(3, 1) ) THEN
579         CALL ctl_stop('dyn:vor_een: requested workspace arrays unavailable')   ;   RETURN
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      IF( wrk_not_released(2, 1,2,3,4,5,6,7) .OR.   &
680          wrk_not_released(3, 1)             )   CALL ctl_stop('dyn:vor_een: failed to release workspace arrays')
681      !
682   END SUBROUTINE vor_een
683
684
685   SUBROUTINE dyn_vor_init
686      !!---------------------------------------------------------------------
687      !!                  ***  ROUTINE dyn_vor_init  ***
688      !!
689      !! ** Purpose :   Control the consistency between cpp options for
690      !!              tracer advection schemes
691      !!----------------------------------------------------------------------
692      INTEGER ::   ioptio          ! local integer
693      !!
694      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
695      !!----------------------------------------------------------------------
696
697      REWIND ( numnam )               ! Read Namelist namdyn_vor : Vorticity scheme options
698      READ   ( numnam, namdyn_vor )
699
700      IF(lwp) THEN                    ! Namelist print
701         WRITE(numout,*)
702         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
703         WRITE(numout,*) '~~~~~~~~~~~~'
704         WRITE(numout,*) '        Namelist namdyn_vor : oice of the vorticity term scheme'
705         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
706         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
707         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
708         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
709      ENDIF
710
711      ioptio = 0                     ! Control of vorticity scheme options
712      IF( ln_dynvor_ene )   ioptio = ioptio + 1
713      IF( ln_dynvor_ens )   ioptio = ioptio + 1
714      IF( ln_dynvor_mix )   ioptio = ioptio + 1
715      IF( ln_dynvor_een )   ioptio = ioptio + 1
716      IF( lk_esopa      )   ioptio =          1
717
718      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
719
720      !                              ! Set nvor (type of scheme for vorticity)
721      IF( ln_dynvor_ene )   nvor =  0
722      IF( ln_dynvor_ens )   nvor =  1
723      IF( ln_dynvor_mix )   nvor =  2
724      IF( ln_dynvor_een )   nvor =  3
725      IF( lk_esopa      )   nvor = -1
726     
727      !                              ! Set ncor, nrvm, ntot (type of vorticity)
728      IF(lwp) WRITE(numout,*)
729      ncor = 1
730      IF( ln_dynadv_vec ) THEN     
731         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
732         nrvm = 2
733         ntot = 4
734      ELSE                       
735         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
736         nrvm = 3
737         ntot = 5
738      ENDIF
739     
740      IF(lwp) THEN                   ! Print the choice
741         WRITE(numout,*)
742         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
743         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
744         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
745         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
746         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
747      ENDIF
748      !
749   END SUBROUTINE dyn_vor_init
750
751   !!==============================================================================
752END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.