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

source: trunk/NEMO/OPA_SRC/DYN/dynvor.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 16 years ago

get back to the nemo_v2_3 version for trunk

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