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

Last change on this file since 553 was 541, checked in by opalod, 18 years ago

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