source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN/dynvor.F90 @ 10293

Last change on this file since 10293 was 10293, checked in by smasson, 23 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 5: reduce communications in oce, see #2133

  • Property svn:keywords set to Id
File size: 53.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   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
23   !!----------------------------------------------------------------------
24
25   !!----------------------------------------------------------------------
26   !!   dyn_vor       : Update the momentum trend with the vorticity trend
27   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
28   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
29   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
30   !!   dyn_vor_init  : set and control of the different vorticity option
31   !!----------------------------------------------------------------------
32   USE oce            ! ocean dynamics and tracers
33   USE dom_oce        ! ocean space and time domain
34   USE dommsk         ! ocean mask
35   USE dynadv         ! momentum advection
36   USE trd_oce        ! trends: ocean variables
37   USE trddyn         ! trend manager: dynamics
38   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
39   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
40   !
41   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl         ! Print control
43   USE in_out_manager ! I/O manager
44   USE lib_mpp        ! MPP library
45   USE timing         ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dyn_vor        ! routine called by step.F90
51   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
52
53   !                                   !!* Namelist namdyn_vor: vorticity term
54   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
55   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
56   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
57   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
58   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
59   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
60   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
61   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
62
63   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
64   !                                       ! associated indices:
65   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
66   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
67   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
68   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
69   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
70   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
71
72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
73   !                               ! associated indices:
74   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
75   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
76   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
77   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
78   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
79
80   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
82   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
84   
85   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
86   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
87   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
88   
89   !! * Substitutions
90#  include "vectopt_loop_substitute.h90"
91   !!----------------------------------------------------------------------
92   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
93   !! $Id$
94   !! Software governed by the CeCILL license (see ./LICENSE)
95   !!----------------------------------------------------------------------
96CONTAINS
97
98   SUBROUTINE dyn_vor( kt )
99      !!----------------------------------------------------------------------
100      !!
101      !! ** Purpose :   compute the lateral ocean tracer physics.
102      !!
103      !! ** Action : - Update (ua,va) with the now vorticity term trend
104      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
105      !!               and planetary vorticity trends) and send them to trd_dyn
106      !!               for futher diagnostics (l_trddyn=T)
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
109      !
110      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
111      !!----------------------------------------------------------------------
112      !
113      IF( ln_timing )   CALL timing_start('dyn_vor')
114      !
115      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
116         !
117         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
118         !
119         ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
120         ztrdv(:,:,:) = va(:,:,:)
121         SELECT CASE( nvor_scheme )
122         CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme
123            IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
124         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme
125            IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
126         CASE( np_ENT )           ;   CALL vor_enT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts)
127            IF( ln_stcor )            CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
128         CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t)
129            IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
130         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme
131            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
132         END SELECT
133         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
134         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
135         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
136         !
137         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
138            ztrdu(:,:,:) = ua(:,:,:)
139            ztrdv(:,:,:) = va(:,:,:)
140            SELECT CASE( nvor_scheme )
141            CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts)
142            CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t)
143            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme
144            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme
145            CASE( np_EEN )           ;   CALL vor_een( kt, nrvm, un , vn , ua, va )  ! energy & enstrophy scheme
146            END SELECT
147            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
148            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
149            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
150         ENDIF
151         !
152         DEALLOCATE( ztrdu, ztrdv )
153         !
154      ELSE              !==  total vorticity trend added to the general trend  ==!
155         !
156         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
157         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
158                             CALL vor_enT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
159            IF( ln_stcor )   CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
160         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
161                             CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
162            IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
163         CASE( np_ENE )                        !* energy conserving scheme
164                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend
165            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
166         CASE( np_ENS )                        !* enstrophy conserving scheme
167                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend
168            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
169         CASE( np_MIX )                        !* mixed ene-ens scheme
170                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
171                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
172            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
173         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
174                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend
175            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
176         END SELECT
177         !
178      ENDIF
179      !
180      !                       ! print sum trends (used for debugging)
181      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
182         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
183      !
184      IF( ln_timing )   CALL timing_stop('dyn_vor')
185      !
186   END SUBROUTINE dyn_vor
187
188
189   SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs )
190      !!----------------------------------------------------------------------
191      !!                  ***  ROUTINE vor_enT  ***
192      !!
193      !! ** Purpose :   Compute the now total vorticity trend and add it to
194      !!      the general trend of the momentum equation.
195      !!
196      !! ** Method  :   Trend evaluated using now fields (centered in time)
197      !!       and t-point evaluation of vorticity (planetary and relative).
198      !!       conserves the horizontal kinetic energy.
199      !!         The general trend of momentum is increased due to the vorticity
200      !!       term which is given by:
201      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
202      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
203      !!       where rvor is the relative vorticity at f-point
204      !!
205      !! ** Action : - Update (ua,va) with the now vorticity term trend
206      !!----------------------------------------------------------------------
207      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
208      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
209      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
210      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
211      !
212      INTEGER  ::   ji, jj, jk           ! dummy loop indices
213      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
214      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
215      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace
216      !!----------------------------------------------------------------------
217      !
218      IF( kt == nit000 ) THEN
219         IF(lwp) WRITE(numout,*)
220         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
221         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
222      ENDIF
223      !
224      !
225      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
226      CASE ( np_RVO )                           !* relative vorticity
227         DO jk = 1, jpkm1                                 ! Horizontal slab
228            DO jj = 1, jpjm1
229               DO ji = 1, jpim1
230                  zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
231                     &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
232               END DO
233            END DO
234            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
235               DO jj = 1, jpjm1
236                  DO ji = 1, jpim1
237                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
238                  END DO
239               END DO
240            ENDIF
241         END DO
242
243         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
244
245      CASE ( np_CRV )                           !* Coriolis + relative vorticity
246         DO jk = 1, jpkm1                                 ! Horizontal slab
247            DO jj = 1, jpjm1
248               DO ji = 1, jpim1                          ! relative vorticity
249                  zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
250                     &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
251               END DO
252            END DO
253            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
254               DO jj = 1, jpjm1
255                  DO ji = 1, jpim1
256                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
257                  END DO
258               END DO
259            ENDIF
260         END DO
261
262         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
263
264      END SELECT
265
266      !                                                ! ===============
267      DO jk = 1, jpkm1                                 ! Horizontal slab
268      !                                                ! ===============
269
270         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
271         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
272            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk)
273         CASE ( np_RVO )                           !* relative vorticity
274            DO jj = 2, jpj
275               DO ji = 2, jpi   ! vector opt.
276                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
277                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
278               END DO
279            END DO
280         CASE ( np_MET )                           !* metric term
281            DO jj = 2, jpj
282               DO ji = 2, jpi
283                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
284                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk)
285               END DO
286            END DO
287         CASE ( np_CRV )                           !* Coriolis + relative vorticity
288            DO jj = 2, jpj
289               DO ji = 2, jpi   ! vector opt.
290                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
291                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
292               END DO
293            END DO
294         CASE ( np_CME )                           !* Coriolis + metric
295            DO jj = 2, jpj
296               DO ji = 2, jpi   ! vector opt.
297                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
298                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
299                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk)
300               END DO
301            END DO
302         CASE DEFAULT                                             ! error
303            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
304         END SELECT
305         !
306         !                                   !==  compute and add the vorticity term trend  =!
307         DO jj = 2, jpjm1
308            DO ji = 2, jpim1   ! vector opt.
309               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    &
310                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
311                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
312                  !
313               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    &
314                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
315                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
316            END DO 
317         END DO 
318         !                                             ! ===============
319      END DO                                           !   End of slab
320      !                                                ! ===============
321   END SUBROUTINE vor_enT
322
323
324   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
325      !!----------------------------------------------------------------------
326      !!                  ***  ROUTINE vor_ene  ***
327      !!
328      !! ** Purpose :   Compute the now total vorticity trend and add it to
329      !!      the general trend of the momentum equation.
330      !!
331      !! ** Method  :   Trend evaluated using now fields (centered in time)
332      !!       and the Sadourny (1975) flux form formulation : conserves the
333      !!       horizontal kinetic energy.
334      !!         The general trend of momentum is increased due to the vorticity
335      !!       term which is given by:
336      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
337      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
338      !!       where rvor is the relative vorticity
339      !!
340      !! ** Action : - Update (ua,va) with the now vorticity term trend
341      !!
342      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
343      !!----------------------------------------------------------------------
344      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
345      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
346      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
347      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
348      !
349      INTEGER  ::   ji, jj, jk           ! dummy loop indices
350      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
351      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
352      !!----------------------------------------------------------------------
353      !
354      IF( kt == nit000 ) THEN
355         IF(lwp) WRITE(numout,*)
356         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
357         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
358      ENDIF
359      !
360      !                                                ! ===============
361      DO jk = 1, jpkm1                                 ! Horizontal slab
362         !                                             ! ===============
363         !
364         SELECT CASE( kvor )                 !==  vorticity considered  ==!
365         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
366            zwz(:,:) = ff_f(:,:) 
367         CASE ( np_RVO )                           !* relative vorticity
368            DO jj = 1, jpjm1
369               DO ji = 1, fs_jpim1   ! vector opt.
370                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
371                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
372               END DO
373            END DO
374         CASE ( np_MET )                           !* metric term
375            DO jj = 1, jpjm1
376               DO ji = 1, fs_jpim1   ! vector opt.
377                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
378                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
379               END DO
380            END DO
381         CASE ( np_CRV )                           !* Coriolis + relative vorticity
382            DO jj = 1, jpjm1
383               DO ji = 1, fs_jpim1   ! vector opt.
384                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
385                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
386               END DO
387            END DO
388         CASE ( np_CME )                           !* Coriolis + metric
389            DO jj = 1, jpjm1
390               DO ji = 1, fs_jpim1   ! vector opt.
391                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
392                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
393               END DO
394            END DO
395         CASE DEFAULT                                             ! error
396            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
397         END SELECT
398         !
399         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
400            DO jj = 1, jpjm1
401               DO ji = 1, fs_jpim1   ! vector opt.
402                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
403               END DO
404            END DO
405         ENDIF
406
407         IF( ln_sco ) THEN
408            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
409            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
410            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
411         ELSE
412            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
413            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
414         ENDIF
415         !                                   !==  compute and add the vorticity term trend  =!
416         DO jj = 2, jpjm1
417            DO ji = fs_2, fs_jpim1   ! vector opt.
418               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
419               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
420               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
421               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
422               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
423               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
424            END DO 
425         END DO 
426         !                                             ! ===============
427      END DO                                           !   End of slab
428      !                                                ! ===============
429   END SUBROUTINE vor_ene
430
431
432   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
433      !!----------------------------------------------------------------------
434      !!                ***  ROUTINE vor_ens  ***
435      !!
436      !! ** Purpose :   Compute the now total vorticity trend and add it to
437      !!      the general trend of the momentum equation.
438      !!
439      !! ** Method  :   Trend evaluated using now fields (centered in time)
440      !!      and the Sadourny (1975) flux FORM formulation : conserves the
441      !!      potential enstrophy of a horizontally non-divergent flow. the
442      !!      trend of the vorticity term is given by:
443      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
444      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
445      !!      Add this trend to the general momentum trend (ua,va):
446      !!          (ua,va) = (ua,va) + ( voru , vorv )
447      !!
448      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
449      !!
450      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
451      !!----------------------------------------------------------------------
452      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
453      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
454      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
455      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
456      !
457      INTEGER  ::   ji, jj, jk   ! dummy loop indices
458      REAL(wp) ::   zuav, zvau   ! local scalars
459      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
460      !!----------------------------------------------------------------------
461      !
462      IF( kt == nit000 ) THEN
463         IF(lwp) WRITE(numout,*)
464         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
465         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
466      ENDIF
467      !                                                ! ===============
468      DO jk = 1, jpkm1                                 ! Horizontal slab
469         !                                             ! ===============
470         !
471         SELECT CASE( kvor )                 !==  vorticity considered  ==!
472         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
473            zwz(:,:) = ff_f(:,:) 
474         CASE ( np_RVO )                           !* relative vorticity
475            DO jj = 1, jpjm1
476               DO ji = 1, fs_jpim1   ! vector opt.
477                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
478                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
479               END DO
480            END DO
481         CASE ( np_MET )                           !* metric term
482            DO jj = 1, jpjm1
483               DO ji = 1, fs_jpim1   ! vector opt.
484                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
485                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
486               END DO
487            END DO
488         CASE ( np_CRV )                           !* Coriolis + relative vorticity
489            DO jj = 1, jpjm1
490               DO ji = 1, fs_jpim1   ! vector opt.
491                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
492                     &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
493               END DO
494            END DO
495         CASE ( np_CME )                           !* Coriolis + metric
496            DO jj = 1, jpjm1
497               DO ji = 1, fs_jpim1   ! vector opt.
498                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
499                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
500               END DO
501            END DO
502         CASE DEFAULT                                             ! error
503            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
504         END SELECT
505         !
506         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
507            DO jj = 1, jpjm1
508               DO ji = 1, fs_jpim1   ! vector opt.
509                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
510               END DO
511            END DO
512         ENDIF
513         !
514         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
515            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
516            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
517            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
518         ELSE
519            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
520            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
521         ENDIF
522         !                                   !==  compute and add the vorticity term trend  =!
523         DO jj = 2, jpjm1
524            DO ji = fs_2, fs_jpim1   ! vector opt.
525               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
526                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
527               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
528                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
529               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
530               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
531            END DO 
532         END DO 
533         !                                             ! ===============
534      END DO                                           !   End of slab
535      !                                                ! ===============
536   END SUBROUTINE vor_ens
537
538
539   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
540      !!----------------------------------------------------------------------
541      !!                ***  ROUTINE vor_een  ***
542      !!
543      !! ** Purpose :   Compute the now total vorticity trend and add it to
544      !!      the general trend of the momentum equation.
545      !!
546      !! ** Method  :   Trend evaluated using now fields (centered in time)
547      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
548      !!      both the horizontal kinetic energy and the potential enstrophy
549      !!      when horizontal divergence is zero (see the NEMO documentation)
550      !!      Add this trend to the general momentum trend (ua,va).
551      !!
552      !! ** Action : - Update (ua,va) with the now vorticity term trend
553      !!
554      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
555      !!----------------------------------------------------------------------
556      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
557      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
558      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
559      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
560      !
561      INTEGER  ::   ji, jj, jk   ! dummy loop indices
562      INTEGER  ::   ierr         ! local integer
563      REAL(wp) ::   zua, zva     ! local scalars
564      REAL(wp) ::   zmsk, ze3f   ! local scalars
565      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
566      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
567      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
568      !!----------------------------------------------------------------------
569      !
570      IF( kt == nit000 ) THEN
571         IF(lwp) WRITE(numout,*)
572         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
573         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
574      ENDIF
575      !
576      !                                                ! ===============
577      DO jk = 1, jpkm1                                 ! Horizontal slab
578         !                                             ! ===============
579         !
580         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
581         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
582            DO jj = 1, jpjm1
583               DO ji = 1, fs_jpim1   ! vector opt.
584                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
585                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
586                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
587                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
588                  ENDIF
589               END DO
590            END DO
591         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
592            DO jj = 1, jpjm1
593               DO ji = 1, fs_jpim1   ! vector opt.
594                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
595                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
596                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
597                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
598                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
599                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
600                  ENDIF
601               END DO
602            END DO
603         END SELECT
604         !
605         SELECT CASE( kvor )                 !==  vorticity considered  ==!
606         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
607            DO jj = 1, jpjm1
608               DO ji = 1, fs_jpim1   ! vector opt.
609                  zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
610               END DO
611            END DO
612         CASE ( np_RVO )                           !* relative vorticity
613            DO jj = 1, jpjm1
614               DO ji = 1, fs_jpim1   ! vector opt.
615                  zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
616                     &            - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
617               END DO
618            END DO
619         CASE ( np_MET )                           !* metric term
620            DO jj = 1, jpjm1
621               DO ji = 1, fs_jpim1   ! vector opt.
622                  zwz(ji,jj,jk) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
623                     &              - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
624               END DO
625            END DO
626         CASE ( np_CRV )                           !* Coriolis + relative vorticity
627            DO jj = 1, jpjm1
628               DO ji = 1, fs_jpim1   ! vector opt.
629                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
630                     &                              - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   &
631                     &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
632               END DO
633            END DO
634         CASE ( np_CME )                           !* Coriolis + metric
635            DO jj = 1, jpjm1
636               DO ji = 1, fs_jpim1   ! vector opt.
637                  zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
638                     &                            - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
639               END DO
640            END DO
641         CASE DEFAULT                                             ! error
642            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
643         END SELECT
644         !
645         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
646            DO jj = 1, jpjm1
647               DO ji = 1, fs_jpim1   ! vector opt.
648                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
649               END DO
650            END DO
651         ENDIF
652      END DO                                           !   End of slab
653         !
654      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
655
656      DO jk = 1, jpkm1                                 ! Horizontal slab
657         !
658         !                                   !==  horizontal fluxes  ==!
659         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
660         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
661
662         !                                   !==  compute and add the vorticity term trend  =!
663         jj = 2
664         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
665         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
666               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
667               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
668               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
669               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
670         END DO
671         DO jj = 3, jpj
672            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
673               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
674               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
675               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
676               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
677            END DO
678         END DO
679         DO jj = 2, jpjm1
680            DO ji = fs_2, fs_jpim1   ! vector opt.
681               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
682                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
683               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
684                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
685               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
686               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
687            END DO 
688         END DO 
689         !                                             ! ===============
690      END DO                                           !   End of slab
691      !                                                ! ===============
692   END SUBROUTINE vor_een
693
694
695
696   SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva )
697      !!----------------------------------------------------------------------
698      !!                ***  ROUTINE vor_eeT  ***
699      !!
700      !! ** Purpose :   Compute the now total vorticity trend and add it to
701      !!      the general trend of the momentum equation.
702      !!
703      !! ** Method  :   Trend evaluated using now fields (centered in time)
704      !!      and the Arakawa and Lamb (1980) vector form formulation using
705      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
706      !!      The change consists in
707      !!      Add this trend to the general momentum trend (ua,va).
708      !!
709      !! ** Action : - Update (ua,va) with the now vorticity term trend
710      !!
711      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
712      !!----------------------------------------------------------------------
713      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
714      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
715      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
716      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
717      !
718      INTEGER  ::   ji, jj, jk     ! dummy loop indices
719      INTEGER  ::   ierr           ! local integer
720      REAL(wp) ::   zua, zva       ! local scalars
721      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
722      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
723      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
724      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
725      !!----------------------------------------------------------------------
726      !
727      IF( kt == nit000 ) THEN
728         IF(lwp) WRITE(numout,*)
729         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
730         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
731      ENDIF
732      !
733      !                                                ! ===============
734      DO jk = 1, jpkm1                                 ! Horizontal slab
735         !                                             ! ===============
736         !
737         !
738         SELECT CASE( kvor )                 !==  vorticity considered  ==!
739         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
740            DO jj = 1, jpjm1
741               DO ji = 1, fs_jpim1   ! vector opt.
742                  zwz(ji,jj,jk) = ff_f(ji,jj)
743               END DO
744            END DO
745         CASE ( np_RVO )                           !* relative vorticity
746            DO jj = 1, jpjm1
747               DO ji = 1, fs_jpim1   ! vector opt.
748                  zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
749                     &             - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
750                     &          * r1_e1e2f(ji,jj)
751               END DO
752            END DO
753         CASE ( np_MET )                           !* metric term
754            DO jj = 1, jpjm1
755               DO ji = 1, fs_jpim1   ! vector opt.
756                  zwz(ji,jj,jk) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
757                     &          - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
758               END DO
759            END DO
760         CASE ( np_CRV )                           !* Coriolis + relative vorticity
761            DO jj = 1, jpjm1
762               DO ji = 1, fs_jpim1   ! vector opt.
763                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
764                     &                              - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
765                     &                         * r1_e1e2f(ji,jj)    )
766               END DO
767            END DO
768         CASE ( np_CME )                           !* Coriolis + metric
769            DO jj = 1, jpjm1
770               DO ji = 1, fs_jpim1   ! vector opt.
771                  zwz(ji,jj,jk) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
772                     &                        - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
773               END DO
774            END DO
775         CASE DEFAULT                                             ! error
776            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
777         END SELECT
778         !
779         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
780            DO jj = 1, jpjm1
781               DO ji = 1, fs_jpim1   ! vector opt.
782                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
783               END DO
784            END DO
785         ENDIF
786      END DO
787      !
788      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
789      !
790      DO jk = 1, jpkm1                                 ! Horizontal slab
791
792      !                                   !==  horizontal fluxes  ==!
793         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
794         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
795
796         !                                   !==  compute and add the vorticity term trend  =!
797         jj = 2
798         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
799         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
800               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
801               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
802               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
803               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
804               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
805         END DO
806         DO jj = 3, jpj
807            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
808               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
809               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
810               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
811               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
812               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
813            END DO
814         END DO
815         DO jj = 2, jpjm1
816            DO ji = fs_2, fs_jpim1   ! vector opt.
817               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
818                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
819               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
820                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
821               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
822               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
823            END DO 
824         END DO 
825         !                                             ! ===============
826      END DO                                           !   End of slab
827      !                                                ! ===============
828   END SUBROUTINE vor_eeT
829
830
831   SUBROUTINE dyn_vor_init
832      !!---------------------------------------------------------------------
833      !!                  ***  ROUTINE dyn_vor_init  ***
834      !!
835      !! ** Purpose :   Control the consistency between cpp options for
836      !!              tracer advection schemes
837      !!----------------------------------------------------------------------
838      INTEGER ::   ji, jj, jk    ! dummy loop indices
839      INTEGER ::   ioptio, ios   ! local integer
840      !!
841      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
842         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
843      !!----------------------------------------------------------------------
844      !
845      IF(lwp) THEN
846         WRITE(numout,*)
847         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
848         WRITE(numout,*) '~~~~~~~~~~~~'
849      ENDIF
850      !
851      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
852      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
853901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
854      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
855      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
856902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
857      IF(lwm) WRITE ( numond, namdyn_vor )
858      !
859      IF(lwp) THEN                    ! Namelist print
860         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
861         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
862         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
863         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
864         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
865         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
866         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
867         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
868         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
869      ENDIF
870
871      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
872
873!!gm  this should be removed when choosing a unique strategy for fmask at the coast
874      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
875      ! at angles with three ocean points and one land point
876      IF(lwp) WRITE(numout,*)
877      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
878      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
879         DO jk = 1, jpk
880            DO jj = 1, jpjm1
881               DO ji = 1, jpim1
882                  IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
883                     & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
884               END DO
885            END DO
886         END DO
887         !
888         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
889         !
890      ENDIF
891!!gm end
892
893      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
894      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
895      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
896      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
897      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
898      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
899      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
900      !
901      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
902      !                     
903      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
904      ncor = np_COR                       ! planetary vorticity
905      SELECT CASE( n_dynadv )
906      CASE( np_LIN_dyn )
907         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
908         nrvm = np_COR        ! planetary vorticity
909         ntot = np_COR        !     -         -
910      CASE( np_VEC_c2  )
911         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
912         nrvm = np_RVO        ! relative vorticity
913         ntot = np_CRV        ! relative + planetary vorticity         
914      CASE( np_FLX_c2 , np_FLX_ubs  )
915         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
916         nrvm = np_MET        ! metric term
917         ntot = np_CME        ! Coriolis + metric term
918         !
919         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
920         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
921            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
922            DO jj = 2, jpjm1
923               DO ji = 2, jpim1
924                  di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
925                  dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
926               END DO
927            END DO
928            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
929            !
930         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
931            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
932            DO jj = 1, jpjm1
933               DO ji = 1, jpim1
934                  di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
935                  dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
936               END DO
937            END DO
938            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
939         END SELECT
940         !
941      END SELECT
942     
943      IF(lwp) THEN                   ! Print the choice
944         WRITE(numout,*)
945         SELECT CASE( nvor_scheme )
946         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
947         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
948         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
949         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
950         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
951         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
952         END SELECT         
953      ENDIF
954      !
955   END SUBROUTINE dyn_vor_init
956
957   !!==============================================================================
958END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.