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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 7 years ago

ORCA2_LIM_PISCES hybrid version update

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