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

Last change on this file since 109 was 108, checked in by opalod, 20 years ago

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