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

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

#1260: LDF simplification + bilap iso-neutral for TRA and GYRE
  • Property svn:keywords set to Id
File size: 35.3 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)  ) &
248                     &       / ( e1f(ji,jj) * e2f(ji,jj) )
249               END DO
250            END DO
251         CASE ( 3 )                                         ! metric term
252            DO jj = 1, jpjm1
253               DO ji = 1, fs_jpim1   ! vector opt.
254                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
255                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
256                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
257               END DO
258            END DO
259         CASE ( 4 )                                         ! total ( planetary + relative vorticity)   (no fmask)
260            DO jj = 1, jpjm1
261               DO ji = 1, fs_jpim1   ! vector opt.
262                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
263                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
264                     &                   / ( e1f(ji,jj) * e2f(ji,jj) )
265               END DO
266            END DO
267         CASE ( 5 )                                         ! total (coriolis + metric)
268            DO jj = 1, jpjm1
269               DO ji = 1, fs_jpim1   ! vector opt.
270                  zwz(ji,jj) = ff(ji,jj)                                                                        &
271                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
272                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
273                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
274               END DO
275            END DO
276         CASE DEFAULT                                             ! error
277            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
278         END SELECT
279         !
280         IF( ln_sco ) THEN
281            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
282            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
283            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
284         ELSE
285            zwx(:,:) = e2u(:,:) * un(:,:,jk)
286            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
287         ENDIF
288         !                                   !==  compute and add the vorticity term trend  =!
289         DO jj = 2, jpjm1
290            DO ji = fs_2, fs_jpim1   ! vector opt.
291               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
292               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
293               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
294               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
295               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
296               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
297            END DO 
298         END DO 
299         !                                             ! ===============
300      END DO                                           !   End of slab
301      !                                                ! ===============
302      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
303      !
304      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
305      !
306   END SUBROUTINE vor_ene
307
308
309   SUBROUTINE vor_ens( kt, kvor, pua, pva )
310      !!----------------------------------------------------------------------
311      !!                ***  ROUTINE vor_ens  ***
312      !!
313      !! ** Purpose :   Compute the now total vorticity trend and add it to
314      !!      the general trend of the momentum equation.
315      !!
316      !! ** Method  :   Trend evaluated using now fields (centered in time)
317      !!      and the Sadourny (1975) flux FORM formulation : conserves the
318      !!      potential enstrophy of a horizontally non-divergent flow. the
319      !!      trend of the vorticity term is given by:
320      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
321      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
322      !!      Add this trend to the general momentum trend (ua,va):
323      !!          (ua,va) = (ua,va) + ( voru , vorv )
324      !!
325      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
326      !!
327      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
328      !!----------------------------------------------------------------------
329      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
330      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
331         !                                                        ! =nrvm (relative vorticity or metric)
332      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
333      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
334      !
335      INTEGER  ::   ji, jj, jk   ! dummy loop indices
336      REAL(wp) ::   zuav, zvau   ! local scalars
337      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
338      !!----------------------------------------------------------------------
339      !
340      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
341      !
342      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
343      !
344      IF( kt == nit000 ) THEN
345         IF(lwp) WRITE(numout,*)
346         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
347         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
348      ENDIF
349      !                                                ! ===============
350      DO jk = 1, jpkm1                                 ! Horizontal slab
351         !                                             ! ===============
352         SELECT CASE( kvor )                 !==  vorticity considered  ==!
353         CASE ( 1 )                                         ! planetary vorticity (Coriolis)
354            zwz(:,:) = ff(:,:) 
355         CASE ( 2 )                                         ! relative  vorticity (no fmask)
356            DO jj = 1, jpjm1
357               DO ji = 1, fs_jpim1   ! vector opt.
358                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
359                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
360                     &       / ( e1f(ji,jj) * e2f(ji,jj) )
361               END DO
362            END DO
363         CASE ( 3 )                                         ! metric term
364            DO jj = 1, jpjm1
365               DO ji = 1, fs_jpim1   ! vector opt.
366                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
367                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
368                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
369               END DO
370            END DO
371         CASE ( 4 )                                         ! total ( planetary + relative vorticity)   (no fmask)
372            DO jj = 1, jpjm1
373               DO ji = 1, fs_jpim1   ! vector opt.
374                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
375                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
376                     &                   / ( e1f(ji,jj) * e2f(ji,jj) )
377               END DO
378            END DO
379         CASE ( 5 )                                         ! total (coriolis + metric)
380            DO jj = 1, jpjm1
381               DO ji = 1, fs_jpim1   ! vector opt.
382                  zwz(ji,jj) = ff(ji,jj)                                                                       &
383                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
384                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
385                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
386               END DO
387            END DO
388         CASE DEFAULT                                             ! error
389            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
390         END SELECT
391         !
392         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
393            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
394            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
395            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
396         ELSE
397            zwx(:,:) = e2u(:,:) * un(:,:,jk)
398            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
399         ENDIF
400         !                                   !==  compute and add the vorticity term trend  =!
401         DO jj = 2, jpjm1
402            DO ji = fs_2, fs_jpim1   ! vector opt.
403               zuav = r1_8 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
404                  &                       + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
405               zvau =-r1_8 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
406                  &                       + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
407               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
408               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
409            END DO 
410         END DO 
411         !                                             ! ===============
412      END DO                                           !   End of slab
413      !                                                ! ===============
414      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
415      !
416      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
417      !
418   END SUBROUTINE vor_ens
419
420
421   SUBROUTINE vor_een( kt, kvor, pua, pva )
422      !!----------------------------------------------------------------------
423      !!                ***  ROUTINE vor_een  ***
424      !!
425      !! ** Purpose :   Compute the now total vorticity trend and add it to
426      !!      the general trend of the momentum equation.
427      !!
428      !! ** Method  :   Trend evaluated using now fields (centered in time)
429      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
430      !!      both the horizontal kinetic energy and the potential enstrophy
431      !!      when horizontal divergence is zero (see the NEMO documentation)
432      !!      Add this trend to the general momentum trend (ua,va).
433      !!
434      !! ** Action : - Update (ua,va) with the now vorticity term trend
435      !!
436      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
437      !!----------------------------------------------------------------------
438      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
439      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
440      !                                                           ! =nrvm (relative vorticity or metric)
441      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
442      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
443      !!
444      INTEGER  ::   ji, jj, jk   ! dummy loop indices
445      INTEGER  ::   ierr         ! local integer
446      REAL(wp) ::   zua, zva     ! local scalars
447      !                                                           !  2D workspace
448      REAL(wp), POINTER    , DIMENSION(:,:  )         ::   zwx, zwy, zwz
449      REAL(wp), POINTER    , DIMENSION(:,:  )         ::   ztnw, ztne, ztsw, ztse
450#if defined key_vvl
451      REAL(wp), POINTER    , DIMENSION(:,:,:)         ::   r1_e3f     !  3D workspace (lk_vvl=T)
452#endif
453#if ! defined key_vvl
454      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   ::   r1_e3f     ! lk_vvl=F, r1_e3f=1/e3f saved one for all
455#endif
456      !!----------------------------------------------------------------------
457      !
458      IF( nn_timing == 1 )  CALL timing_start('vor_een')
459      !
460      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
461      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
462#if defined key_vvl
463      CALL wrk_alloc( jpi, jpj, jpk, r1_e3f                 )
464#endif
465      !
466      IF( kt == nit000 ) THEN
467         IF(lwp) WRITE(numout,*)
468         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
469         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
470#if ! defined key_vvl
471         IF( .NOT.ALLOCATED(r1_e3f) ) THEN
472            ALLOCATE( r1_e3f(jpi,jpj,jpk) , STAT=ierr )
473            IF( lk_mpp    )   CALL mpp_sum ( ierr )
474            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
475         ENDIF
476#endif
477      ENDIF
478
479      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
480         DO jk = 1, jpk
481            DO jj = 1, jpjm1
482               DO ji = 1, jpim1
483                  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)   &
484                     &               + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * r1_4
485                  IF( r1_e3f(ji,jj,jk) /= 0._wp )   r1_e3f(ji,jj,jk) = 1._wp / r1_e3f(ji,jj,jk)
486               END DO
487            END DO
488         END DO
489         CALL lbc_lnk( r1_e3f, 'F', 1. )
490      ENDIF
491      !                                                ! ===============
492      DO jk = 1, jpkm1                                 ! Horizontal slab
493         !                                             ! ===============
494         !
495         SELECT CASE( kvor )                 !==  vorticity considered  ==!
496         CASE ( 1 )                                         ! planetary vorticity (Coriolis)
497            zwz(:,:) = ff(:,:) * r1_e3f(:,:,jk)
498         CASE ( 2 )                                         ! relative  vorticity (no fmask)
499            DO jj = 1, jpjm1
500               DO ji = 1, fs_jpim1   ! vector opt.
501                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
502                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
503                     &       / ( e1f(ji,jj) * e2f(ji,jj) ) * r1_e3f(ji,jj,jk)
504               END DO
505            END DO
506            CALL lbc_lnk( zwz, 'F', 1. )
507         CASE ( 3 )                                         ! metric term
508            DO jj = 1, jpjm1
509               DO ji = 1, fs_jpim1   ! vector opt.
510                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
511                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
512                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * r1_e3f(ji,jj,jk)
513               END DO
514            END DO
515            CALL lbc_lnk( zwz, 'F', 1. )
516         CASE ( 4 )                                         ! total ( planetary + relative vorticity)   (no fmask)
517            DO jj = 1, jpjm1
518               DO ji = 1, fs_jpim1   ! vector opt.
519                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
520                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
521                     &                      / ( e1f(ji,jj) * e2f(ji,jj) )    ) * r1_e3f(ji,jj,jk)
522               END DO
523            END DO
524            CALL lbc_lnk( zwz, 'F', 1. )
525         CASE ( 5 )                                         ! total (coriolis + metric)
526            DO jj = 1, jpjm1
527               DO ji = 1, fs_jpim1   ! vector opt.
528                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
529                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
530                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
531                       &        * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )   ) * r1_e3f(ji,jj,jk)
532               END DO
533            END DO
534            CALL lbc_lnk( zwz, 'F', 1. )
535         CASE DEFAULT                                             ! error
536            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
537         END SELECT
538         !
539         !                                   !==  horizontal fluxes  ==!
540         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
541         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
542
543         !                                   !==  compute and add the vorticity term trend  =!
544         jj = 2
545         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
546         DO ji = 2, jpi   
547               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
548               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
549               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
550               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
551         END DO
552         DO jj = 3, jpj
553            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
554               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
555               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
556               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
557               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
558            END DO
559         END DO
560         DO jj = 2, jpjm1
561            DO ji = fs_2, fs_jpim1   ! vector opt.
562               zua = + r1_12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
563                  &                          + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
564               zva = - r1_12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
565                  &                          + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
566               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
567               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
568            END DO 
569         END DO 
570         !                                             ! ===============
571      END DO                                           !   End of slab
572      !                                                ! ===============
573      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
574      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
575#if defined key_vvl
576      CALL wrk_dealloc( jpi, jpj, jpk, r1_e3f                 )
577#endif
578      !
579      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
580      !
581   END SUBROUTINE vor_een
582
583
584   SUBROUTINE dyn_vor_init
585      !!---------------------------------------------------------------------
586      !!                  ***  ROUTINE dyn_vor_init  ***
587      !!
588      !! ** Purpose :   Control the consistency between cpp options for
589      !!              tracer advection schemes
590      !!----------------------------------------------------------------------
591      INTEGER ::   ioptio          ! local integer
592      INTEGER ::   ji, jj, jk      ! dummy loop indices
593      INTEGER ::   ios             ! Local integer output status for namelist read
594      !!
595      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
596      !!----------------------------------------------------------------------
597
598      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
599      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
600901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
601
602      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
603      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
604902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
605      WRITE ( numond, namdyn_vor )
606
607      IF(lwp) THEN                    ! Namelist print
608         WRITE(numout,*)
609         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
610         WRITE(numout,*) '~~~~~~~~~~~~'
611         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
612         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
613         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
614         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
615         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
616      ENDIF
617
618      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
619      ! at angles with three ocean points and one land point
620      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
621         DO jk = 1, jpk
622            DO jj = 2, jpjm1
623               DO ji = 2, jpim1
624                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
625                      fmask(ji,jj,jk) = 1._wp
626               END DO
627            END DO
628         END DO
629          !
630          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
631          !
632      ENDIF
633
634      ioptio = 0                     ! Control of vorticity scheme options
635      IF( ln_dynvor_ene )   ioptio = ioptio + 1
636      IF( ln_dynvor_ens )   ioptio = ioptio + 1
637      IF( ln_dynvor_mix )   ioptio = ioptio + 1
638      IF( ln_dynvor_een )   ioptio = ioptio + 1
639      IF( lk_esopa      )   ioptio =          1
640
641      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
642
643      !                              ! Set nvor (type of scheme for vorticity)
644      IF( ln_dynvor_ene )   nvor =  0
645      IF( ln_dynvor_ens )   nvor =  1
646      IF( ln_dynvor_mix )   nvor =  2
647      IF( ln_dynvor_een )   nvor =  3
648      IF( lk_esopa      )   nvor = -1
649     
650      !                              ! Set ncor, nrvm, ntot (type of vorticity)
651      IF(lwp) WRITE(numout,*)
652      ncor = 1
653      IF( ln_dynadv_vec ) THEN     
654         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
655         nrvm = 2
656         ntot = 4
657      ELSE                       
658         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
659         nrvm = 3
660         ntot = 5
661      ENDIF
662     
663      IF(lwp) THEN                   ! Print the choice
664         WRITE(numout,*)
665         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
666         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
667         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
668         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
669         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
670      ENDIF
671      !
672   END SUBROUTINE dyn_vor_init
673
674   !!==============================================================================
675END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.