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 branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

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