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/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 5908

Last change on this file since 5908 was 5908, checked in by jchanut, 8 years ago

Move lbclk else unreproductible results with ln_dynvor_een=T and ln_dynvor_msk=T

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