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

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago

GYRE hybrid parallelization
  • Property svn:keywords set to Id
File size: 37.8 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            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
325            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
326            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
327         ELSE
328!$OMP PARALLEL WORKSHARE
329            zwx(:,:) = e2u(:,:) * un(:,:,jk)
330            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
331!$OMP END PARALLEL WORKSHARE
332         ENDIF
333         !                                   !==  compute and add the vorticity term trend  =!
334!$OMP PARALLEL DO private(jj, ji, zy1, zy2, zx1, zx2)
335         DO jj = 2, jpjm1
336            DO ji = fs_2, fs_jpim1   ! vector opt.
337               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
338               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
339               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
340               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
341               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
342               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
343            END DO 
344         END DO 
345         !                                             ! ===============
346      END DO                                           !   End of slab
347      !                                                ! ===============
348      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
349      !
350      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
351      !
352   END SUBROUTINE vor_ene
353
354
355   SUBROUTINE vor_ens( kt, kvor, pua, pva )
356      !!----------------------------------------------------------------------
357      !!                ***  ROUTINE vor_ens  ***
358      !!
359      !! ** Purpose :   Compute the now total vorticity trend and add it to
360      !!      the general trend of the momentum equation.
361      !!
362      !! ** Method  :   Trend evaluated using now fields (centered in time)
363      !!      and the Sadourny (1975) flux FORM formulation : conserves the
364      !!      potential enstrophy of a horizontally non-divergent flow. the
365      !!      trend of the vorticity term is given by:
366      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
367      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
368      !!      Add this trend to the general momentum trend (ua,va):
369      !!          (ua,va) = (ua,va) + ( voru , vorv )
370      !!
371      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
372      !!
373      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
374      !!----------------------------------------------------------------------
375      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
376      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
377         !                                                        ! =nrvm (relative vorticity or metric)
378      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
379      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
380      !
381      INTEGER  ::   ji, jj, jk   ! dummy loop indices
382      REAL(wp) ::   zuav, zvau   ! local scalars
383      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
384      !!----------------------------------------------------------------------
385      !
386      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
387      !
388      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
389      !
390      IF( kt == nit000 ) THEN
391         IF(lwp) WRITE(numout,*)
392         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
393         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
394      ENDIF
395      !                                                ! ===============
396      DO jk = 1, jpkm1                                 ! Horizontal slab
397         !                                             ! ===============
398         !
399         SELECT CASE( kvor )                 !==  vorticity considered  ==!
400         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
401            zwz(:,:) = ff(:,:) 
402         CASE ( np_RVO )                           !* relative vorticity
403            DO jj = 1, jpjm1
404               DO ji = 1, fs_jpim1   ! vector opt.
405                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
406                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
407               END DO
408            END DO
409         CASE ( np_MET )                           !* metric term
410            DO jj = 1, jpjm1
411               DO ji = 1, fs_jpim1   ! vector opt.
412                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
413                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
414                       &     * 0.5 * r1_e1e2f(ji,jj)
415               END DO
416            END DO
417         CASE ( np_CRV )                           !* Coriolis + relative vorticity
418            DO jj = 1, jpjm1
419               DO ji = 1, fs_jpim1   ! vector opt.
420                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
421                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
422                     &                   * r1_e1e2f(ji,jj)
423               END DO
424            END DO
425         CASE ( np_CME )                           !* Coriolis + metric
426            DO jj = 1, jpjm1
427               DO ji = 1, fs_jpim1   ! vector opt.
428                  zwz(ji,jj) = ff(ji,jj)                                                                       &
429                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
430                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
431                       &     * 0.5 * r1_e1e2f(ji,jj)
432               END DO
433            END DO
434         CASE DEFAULT                                             ! error
435            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
436         END SELECT
437         !
438         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
439            DO jj = 1, jpjm1
440               DO ji = 1, fs_jpim1   ! vector opt.
441                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
442               END DO
443            END DO
444         ENDIF
445         !
446         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
447            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
448            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
449            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
450         ELSE
451            zwx(:,:) = e2u(:,:) * un(:,:,jk)
452            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
453         ENDIF
454         !                                   !==  compute and add the vorticity term trend  =!
455         DO jj = 2, jpjm1
456            DO ji = fs_2, fs_jpim1   ! vector opt.
457               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
458                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
459               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
460                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
461               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
462               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
463            END DO 
464         END DO 
465         !                                             ! ===============
466      END DO                                           !   End of slab
467      !                                                ! ===============
468      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
469      !
470      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
471      !
472   END SUBROUTINE vor_ens
473
474
475   SUBROUTINE vor_een( kt, kvor, pua, pva )
476      !!----------------------------------------------------------------------
477      !!                ***  ROUTINE vor_een  ***
478      !!
479      !! ** Purpose :   Compute the now total vorticity trend and add it to
480      !!      the general trend of the momentum equation.
481      !!
482      !! ** Method  :   Trend evaluated using now fields (centered in time)
483      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
484      !!      both the horizontal kinetic energy and the potential enstrophy
485      !!      when horizontal divergence is zero (see the NEMO documentation)
486      !!      Add this trend to the general momentum trend (ua,va).
487      !!
488      !! ** Action : - Update (ua,va) with the now vorticity term trend
489      !!
490      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
491      !!----------------------------------------------------------------------
492      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
493      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric)
494      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
495      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
496      !
497      INTEGER  ::   ji, jj, jk   ! dummy loop indices
498      INTEGER  ::   ierr         ! local integer
499      REAL(wp) ::   zua, zva     ! local scalars
500      REAL(wp) ::   zmsk, ze3    ! local scalars
501      !
502      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f
503      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse
504      !!----------------------------------------------------------------------
505      !
506      IF( nn_timing == 1 )  CALL timing_start('vor_een')
507      !
508      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
509      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
510      !
511      IF( kt == nit000 ) THEN
512         IF(lwp) WRITE(numout,*)
513         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
514         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
515      ENDIF
516      !
517      !                                                ! ===============
518      DO jk = 1, jpkm1                                 ! Horizontal slab
519         !                                             ! ===============
520         !
521         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
522         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
523            DO jj = 1, jpjm1
524               DO ji = 1, fs_jpim1   ! vector opt.
525                  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)   &
526                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
527                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3
528                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp
529                  ENDIF
530               END DO
531            END DO
532         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
533            DO jj = 1, jpjm1
534               DO ji = 1, fs_jpim1   ! vector opt.
535                  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)   &
536                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
537                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
538                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
539                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3
540                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp
541                  ENDIF
542               END DO
543            END DO
544         END SELECT
545         !
546         SELECT CASE( kvor )                 !==  vorticity considered  ==!
547         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
548            DO jj = 1, jpjm1
549               DO ji = 1, fs_jpim1   ! vector opt.
550                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
551               END DO
552            END DO
553         CASE ( np_RVO )                           !* relative vorticity
554            DO jj = 1, jpjm1
555               DO ji = 1, fs_jpim1   ! vector opt.
556                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
557                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
558                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
559               END DO
560            END DO
561         CASE ( np_MET )                           !* metric term
562            DO jj = 1, jpjm1
563               DO ji = 1, fs_jpim1   ! vector opt.
564                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
565                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
566                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
567               END DO
568            END DO
569         CASE ( np_CRV )                           !* Coriolis + relative vorticity
570            DO jj = 1, jpjm1
571               DO ji = 1, fs_jpim1   ! vector opt.
572                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
573                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
574                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
575               END DO
576            END DO
577         CASE ( np_CME )                           !* Coriolis + metric
578            DO jj = 1, jpjm1
579               DO ji = 1, fs_jpim1   ! vector opt.
580                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
581                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
582                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
583                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
584               END DO
585            END DO
586         CASE DEFAULT                                             ! error
587            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
588         END SELECT
589         !
590         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
591            DO jj = 1, jpjm1
592               DO ji = 1, fs_jpim1   ! vector opt.
593                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
594               END DO
595            END DO
596         ENDIF
597         !
598         CALL lbc_lnk( zwz, 'F', 1. )
599         !
600         !                                   !==  horizontal fluxes  ==!
601         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
602         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
603
604         !                                   !==  compute and add the vorticity term trend  =!
605         jj = 2
606         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
607         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
608               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
609               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
610               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
611               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
612         END DO
613         DO jj = 3, jpj
614            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
615               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
616               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
617               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
618               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
619            END DO
620         END DO
621         DO jj = 2, jpjm1
622            DO ji = fs_2, fs_jpim1   ! vector opt.
623               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
624                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
625               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
626                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
627               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
628               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
629            END DO 
630         END DO 
631         !                                             ! ===============
632      END DO                                           !   End of slab
633      !                                                ! ===============
634      !
635      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
636      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
637      !
638      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
639      !
640   END SUBROUTINE vor_een
641
642
643   SUBROUTINE dyn_vor_init
644      !!---------------------------------------------------------------------
645      !!                  ***  ROUTINE dyn_vor_init  ***
646      !!
647      !! ** Purpose :   Control the consistency between cpp options for
648      !!              tracer advection schemes
649      !!----------------------------------------------------------------------
650      INTEGER ::   ioptio          ! local integer
651      INTEGER ::   ji, jj, jk      ! dummy loop indices
652      INTEGER ::   ios             ! Local integer output status for namelist read
653      !!
654      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
655      !!----------------------------------------------------------------------
656
657      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
658      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
659901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
660
661      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
662      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
663902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
664      IF(lwm) WRITE ( numond, namdyn_vor )
665
666      IF(lwp) THEN                    ! Namelist print
667         WRITE(numout,*)
668         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
669         WRITE(numout,*) '~~~~~~~~~~~~'
670         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
671         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
672         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
673         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
674         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
675         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
676         WRITE(numout,*) '           masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
677      ENDIF
678
679!!gm  this should be removed when choosing a unique strategy for fmask at the coast
680      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
681      ! at angles with three ocean points and one land point
682      IF(lwp) WRITE(numout,*)
683      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
684      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
685         DO jk = 1, jpk
686            DO jj = 2, jpjm1
687               DO ji = 2, jpim1
688                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
689                      fmask(ji,jj,jk) = 1._wp
690               END DO
691            END DO
692         END DO
693          !
694          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
695          !
696      ENDIF
697!!gm end
698
699      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
700      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF
701      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF
702      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF
703      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF
704      !
705      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
706      !                     
707      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
708      ncor = np_COR
709      IF( ln_dynadv_vec ) THEN     
710         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
711         nrvm = np_RVO        ! relative vorticity
712         ntot = np_CRV        ! relative + planetary vorticity
713      ELSE                       
714         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
715         nrvm = np_MET        ! metric term
716         ntot = np_CME        ! Coriolis + metric term
717      ENDIF
718     
719      IF(lwp) THEN                   ! Print the choice
720         WRITE(numout,*)
721         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
722         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
723         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
724         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
725      ENDIF
726      !
727   END SUBROUTINE dyn_vor_init
728
729   !!==============================================================================
730END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.