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

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

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