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

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

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