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

source: branches/2011/dev_CMCC/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3033

Last change on this file since 3033 was 3033, checked in by vichi, 13 years ago

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