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 NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/dev_r9888_GO6_mixing/src/OCE/DYN/dynvor.F90 @ 9900

Last change on this file since 9900 was 9900, checked in by davestorkey, 6 years ago

UKMO dev_r9888_GO6_mixing branch: clear SVN keywords.

File size: 52.3 KB
Line 
1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
17   !!            3.7  ! 2014-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 licence     (./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, zwz, zwt   ! 2D workspace
215      !!----------------------------------------------------------------------
216      !
217      IF( kt == nit000 ) THEN
218         IF(lwp) WRITE(numout,*)
219         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
221      ENDIF
222      !
223      !                                                ! ===============
224      DO jk = 1, jpkm1                                 ! Horizontal slab
225         !                                             ! ===============
226         !
227         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
228         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
229            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk)
230         CASE ( np_RVO )                           !* relative vorticity
231            DO jj = 1, jpjm1
232               DO ji = 1, jpim1
233                  zwz(ji,jj) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
234                     &          - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
235               END DO
236            END DO
237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
238               DO jj = 1, jpjm1
239                  DO ji = 1, jpim1
240                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
241                  END DO
242               END DO
243            ENDIF
244            CALL lbc_lnk( zwz, 'F', 1. )
245            DO jj = 2, jpj
246               DO ji = 2, jpi   ! vector opt.
247                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ) + zwz(ji,jj  )   &
248                     &                  + zwz(ji-1,jj-1) + zwz(ji,jj-1)   ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
249               END DO
250            END DO
251         CASE ( np_MET )                           !* metric term
252            DO jj = 2, jpj
253               DO ji = 2, jpi
254                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
255                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk)
256               END DO
257            END DO
258         CASE ( np_CRV )                           !* Coriolis + relative vorticity
259            DO jj = 1, jpjm1
260               DO ji = 1, jpim1                          ! relative vorticity
261                  zwz(ji,jj) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
262                     &           - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
263               END DO
264            END DO
265            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
266               DO jj = 1, jpjm1
267                  DO ji = 1, jpim1
268                     zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
269                  END DO
270               END DO
271            ENDIF
272            CALL lbc_lnk( zwz, 'F', 1. )
273            DO jj = 2, jpj
274               DO ji = 2, jpi   ! vector opt.
275                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ) + zwz(ji,jj  )    &
276                     &                                 + zwz(ji-1,jj-1) + zwz(ji,jj-1) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
277               END DO
278            END DO
279         CASE ( np_CME )                           !* Coriolis + metric
280            DO jj = 2, jpj
281               DO ji = 2, jpi   ! vector opt.
282                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
283                       &        + ( 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 DEFAULT                                             ! error
288            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
289         END SELECT
290         !
291         !                                   !==  compute and add the vorticity term trend  =!
292         DO jj = 2, jpjm1
293            DO ji = 2, jpim1   ! vector opt.
294               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    &
295                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
296                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
297                  !
298               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    &
299                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
300                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
301            END DO 
302         END DO 
303         !                                             ! ===============
304      END DO                                           !   End of slab
305      !                                                ! ===============
306   END SUBROUTINE vor_enT
307
308
309   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE vor_ene  ***
312      !!
313      !! ** Purpose :   Compute the now total vorticity trend and add it to
314      !!      the general trend of the momentum equation.
315      !!
316      !! ** Method  :   Trend evaluated using now fields (centered in time)
317      !!       and the Sadourny (1975) flux form formulation : conserves the
318      !!       horizontal kinetic energy.
319      !!         The general trend of momentum is increased due to the vorticity
320      !!       term which is given by:
321      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
322      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
323      !!       where rvor is the relative vorticity
324      !!
325      !! ** Action : - Update (ua,va) with the now vorticity term trend
326      !!
327      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
328      !!----------------------------------------------------------------------
329      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
330      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
331      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
332      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
333      !
334      INTEGER  ::   ji, jj, jk           ! dummy loop indices
335      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
336      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
337      !!----------------------------------------------------------------------
338      !
339      IF( kt == nit000 ) THEN
340         IF(lwp) WRITE(numout,*)
341         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
342         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
343      ENDIF
344      !
345      !                                                ! ===============
346      DO jk = 1, jpkm1                                 ! Horizontal slab
347         !                                             ! ===============
348         !
349         SELECT CASE( kvor )                 !==  vorticity considered  ==!
350         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
351            zwz(:,:) = ff_f(:,:) 
352         CASE ( np_RVO )                           !* relative vorticity
353            DO jj = 1, jpjm1
354               DO ji = 1, fs_jpim1   ! vector opt.
355                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
356                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
357               END DO
358            END DO
359         CASE ( np_MET )                           !* metric term
360            DO jj = 1, jpjm1
361               DO ji = 1, fs_jpim1   ! vector opt.
362                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
363                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
364               END DO
365            END DO
366         CASE ( np_CRV )                           !* Coriolis + relative vorticity
367            DO jj = 1, jpjm1
368               DO ji = 1, fs_jpim1   ! vector opt.
369                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
370                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
371               END DO
372            END DO
373         CASE ( np_CME )                           !* Coriolis + metric
374            DO jj = 1, jpjm1
375               DO ji = 1, fs_jpim1   ! vector opt.
376                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
377                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
378               END DO
379            END DO
380         CASE DEFAULT                                             ! error
381            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
382         END SELECT
383         !
384         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
385            DO jj = 1, jpjm1
386               DO ji = 1, fs_jpim1   ! vector opt.
387                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
388               END DO
389            END DO
390         ENDIF
391
392         IF( ln_sco ) THEN
393            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
394            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
395            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
396         ELSE
397            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
398            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
399         ENDIF
400         !                                   !==  compute and add the vorticity term trend  =!
401         DO jj = 2, jpjm1
402            DO ji = fs_2, fs_jpim1   ! vector opt.
403               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
404               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
405               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
406               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
407               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
408               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
409            END DO 
410         END DO 
411         !                                             ! ===============
412      END DO                                           !   End of slab
413      !                                                ! ===============
414   END SUBROUTINE vor_ene
415
416
417   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
418      !!----------------------------------------------------------------------
419      !!                ***  ROUTINE vor_ens  ***
420      !!
421      !! ** Purpose :   Compute the now total vorticity trend and add it to
422      !!      the general trend of the momentum equation.
423      !!
424      !! ** Method  :   Trend evaluated using now fields (centered in time)
425      !!      and the Sadourny (1975) flux FORM formulation : conserves the
426      !!      potential enstrophy of a horizontally non-divergent flow. the
427      !!      trend of the vorticity term is given by:
428      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
429      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
430      !!      Add this trend to the general momentum trend (ua,va):
431      !!          (ua,va) = (ua,va) + ( voru , vorv )
432      !!
433      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
434      !!
435      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
436      !!----------------------------------------------------------------------
437      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
438      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
440      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
441      !
442      INTEGER  ::   ji, jj, jk   ! dummy loop indices
443      REAL(wp) ::   zuav, zvau   ! local scalars
444      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
445      !!----------------------------------------------------------------------
446      !
447      IF( kt == nit000 ) THEN
448         IF(lwp) WRITE(numout,*)
449         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
450         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
451      ENDIF
452      !                                                ! ===============
453      DO jk = 1, jpkm1                                 ! Horizontal slab
454         !                                             ! ===============
455         !
456         SELECT CASE( kvor )                 !==  vorticity considered  ==!
457         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
458            zwz(:,:) = ff_f(:,:) 
459         CASE ( np_RVO )                           !* relative vorticity
460            DO jj = 1, jpjm1
461               DO ji = 1, fs_jpim1   ! vector opt.
462                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
463                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
464               END DO
465            END DO
466         CASE ( np_MET )                           !* metric term
467            DO jj = 1, jpjm1
468               DO ji = 1, fs_jpim1   ! vector opt.
469                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
470                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
471               END DO
472            END DO
473         CASE ( np_CRV )                           !* Coriolis + relative vorticity
474            DO jj = 1, jpjm1
475               DO ji = 1, fs_jpim1   ! vector opt.
476                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
477                     &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
478               END DO
479            END DO
480         CASE ( np_CME )                           !* Coriolis + metric
481            DO jj = 1, jpjm1
482               DO ji = 1, fs_jpim1   ! vector opt.
483                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
484                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
485               END DO
486            END DO
487         CASE DEFAULT                                             ! error
488            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
489         END SELECT
490         !
491         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
492            DO jj = 1, jpjm1
493               DO ji = 1, fs_jpim1   ! vector opt.
494                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
495               END DO
496            END DO
497         ENDIF
498         !
499         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
500            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
501            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
502            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
503         ELSE
504            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
505            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
506         ENDIF
507         !                                   !==  compute and add the vorticity term trend  =!
508         DO jj = 2, jpjm1
509            DO ji = fs_2, fs_jpim1   ! vector opt.
510               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
511                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
512               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
513                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
514               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
515               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
516            END DO 
517         END DO 
518         !                                             ! ===============
519      END DO                                           !   End of slab
520      !                                                ! ===============
521   END SUBROUTINE vor_ens
522
523
524   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
525      !!----------------------------------------------------------------------
526      !!                ***  ROUTINE vor_een  ***
527      !!
528      !! ** Purpose :   Compute the now total vorticity trend and add it to
529      !!      the general trend of the momentum equation.
530      !!
531      !! ** Method  :   Trend evaluated using now fields (centered in time)
532      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
533      !!      both the horizontal kinetic energy and the potential enstrophy
534      !!      when horizontal divergence is zero (see the NEMO documentation)
535      !!      Add this trend to the general momentum trend (ua,va).
536      !!
537      !! ** Action : - Update (ua,va) with the now vorticity term trend
538      !!
539      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
540      !!----------------------------------------------------------------------
541      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
542      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
543      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
544      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
545      !
546      INTEGER  ::   ji, jj, jk   ! dummy loop indices
547      INTEGER  ::   ierr         ! local integer
548      REAL(wp) ::   zua, zva     ! local scalars
549      REAL(wp) ::   zmsk, ze3f   ! local scalars
550      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz , z1_e3f
551      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse
552      !!----------------------------------------------------------------------
553      !
554      IF( kt == nit000 ) THEN
555         IF(lwp) WRITE(numout,*)
556         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
557         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
558      ENDIF
559      !
560      !                                                ! ===============
561      DO jk = 1, jpkm1                                 ! Horizontal slab
562         !                                             ! ===============
563         !
564         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
565         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
566            DO jj = 1, jpjm1
567               DO ji = 1, fs_jpim1   ! vector opt.
568                  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)   &
569                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
570                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
571                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
572                  ENDIF
573               END DO
574            END DO
575         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
576            DO jj = 1, jpjm1
577               DO ji = 1, fs_jpim1   ! vector opt.
578                  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)   &
579                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
580                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
581                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
582                  IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
583                  ELSE                       ;   z1_e3f(ji,jj) = 0._wp
584                  ENDIF
585               END DO
586            END DO
587         END SELECT
588         !
589         SELECT CASE( kvor )                 !==  vorticity considered  ==!
590         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
591            DO jj = 1, jpjm1
592               DO ji = 1, fs_jpim1   ! vector opt.
593                  zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
594               END DO
595            END DO
596         CASE ( np_RVO )                           !* relative vorticity
597            DO jj = 1, jpjm1
598               DO ji = 1, fs_jpim1   ! vector opt.
599                  zwz(ji,jj) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
600                     &         - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
601               END DO
602            END DO
603         CASE ( np_MET )                           !* metric term
604            DO jj = 1, jpjm1
605               DO ji = 1, fs_jpim1   ! vector opt.
606                  zwz(ji,jj) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
607                     &           - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
608               END DO
609            END DO
610         CASE ( np_CRV )                           !* Coriolis + relative vorticity
611            DO jj = 1, jpjm1
612               DO ji = 1, fs_jpim1   ! vector opt.
613                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
614                     &                           - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   &
615                     &                        * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
616               END DO
617            END DO
618         CASE ( np_CME )                           !* Coriolis + metric
619            DO jj = 1, jpjm1
620               DO ji = 1, fs_jpim1   ! vector opt.
621                  zwz(ji,jj) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
622                     &                         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
623               END DO
624            END DO
625         CASE DEFAULT                                             ! error
626            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
627         END SELECT
628         !
629         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
630            DO jj = 1, jpjm1
631               DO ji = 1, fs_jpim1   ! vector opt.
632                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
633               END DO
634            END DO
635         ENDIF
636         !
637         CALL lbc_lnk( zwz, 'F', 1. )
638         !
639         !                                   !==  horizontal fluxes  ==!
640         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
641         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
642
643         !                                   !==  compute and add the vorticity term trend  =!
644         jj = 2
645         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
646         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
647               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
648               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
649               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
650               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
651         END DO
652         DO jj = 3, jpj
653            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
654               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
655               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
656               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
657               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
658            END DO
659         END DO
660         DO jj = 2, jpjm1
661            DO ji = fs_2, fs_jpim1   ! vector opt.
662               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
663                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
664               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
665                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
666               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
667               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
668            END DO 
669         END DO 
670         !                                             ! ===============
671      END DO                                           !   End of slab
672      !                                                ! ===============
673   END SUBROUTINE vor_een
674
675
676
677   SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva )
678      !!----------------------------------------------------------------------
679      !!                ***  ROUTINE vor_eeT  ***
680      !!
681      !! ** Purpose :   Compute the now total vorticity trend and add it to
682      !!      the general trend of the momentum equation.
683      !!
684      !! ** Method  :   Trend evaluated using now fields (centered in time)
685      !!      and the Arakawa and Lamb (1980) vector form formulation using
686      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
687      !!      The change consists in
688      !!      Add this trend to the general momentum trend (ua,va).
689      !!
690      !! ** Action : - Update (ua,va) with the now vorticity term trend
691      !!
692      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
693      !!----------------------------------------------------------------------
694      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
695      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
696      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
697      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
698      !
699      INTEGER  ::   ji, jj, jk     ! dummy loop indices
700      INTEGER  ::   ierr           ! local integer
701      REAL(wp) ::   zua, zva       ! local scalars
702      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
703      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz
704      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse
705      !!----------------------------------------------------------------------
706      !
707      IF( kt == nit000 ) THEN
708         IF(lwp) WRITE(numout,*)
709         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
710         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
711      ENDIF
712      !
713      !                                                ! ===============
714      DO jk = 1, jpkm1                                 ! Horizontal slab
715         !                                             ! ===============
716         !
717         !
718         SELECT CASE( kvor )                 !==  vorticity considered  ==!
719         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
720            DO jj = 1, jpjm1
721               DO ji = 1, fs_jpim1   ! vector opt.
722                  zwz(ji,jj) = ff_f(ji,jj)
723               END DO
724            END DO
725         CASE ( np_RVO )                           !* relative vorticity
726            DO jj = 1, jpjm1
727               DO ji = 1, fs_jpim1   ! vector opt.
728                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
729                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
730                     &       * r1_e1e2f(ji,jj)
731               END DO
732            END DO
733         CASE ( np_MET )                           !* metric term
734            DO jj = 1, jpjm1
735               DO ji = 1, fs_jpim1   ! vector opt.
736                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
737                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
738               END DO
739            END DO
740         CASE ( np_CRV )                           !* Coriolis + relative vorticity
741            DO jj = 1, jpjm1
742               DO ji = 1, fs_jpim1   ! vector opt.
743                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
744                     &                           - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
745                     &                      * r1_e1e2f(ji,jj)    )
746               END DO
747            END DO
748         CASE ( np_CME )                           !* Coriolis + metric
749            DO jj = 1, jpjm1
750               DO ji = 1, fs_jpim1   ! vector opt.
751                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
752                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
753               END DO
754            END DO
755         CASE DEFAULT                                             ! error
756            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
757         END SELECT
758         !
759         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
760            DO jj = 1, jpjm1
761               DO ji = 1, fs_jpim1   ! vector opt.
762                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
763               END DO
764            END DO
765         ENDIF
766         !
767         CALL lbc_lnk( zwz, 'F', 1. )
768         !
769         !                                   !==  horizontal fluxes  ==!
770         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
771         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
772
773         !                                   !==  compute and add the vorticity term trend  =!
774         jj = 2
775         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
776         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
777               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
778               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t
779               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t
780               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
781               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t
782         END DO
783         DO jj = 3, jpj
784            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
785               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
786               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) * z1_e3t
787               ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) * z1_e3t
788               ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
789               ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) * z1_e3t
790            END DO
791         END DO
792         DO jj = 2, jpjm1
793            DO ji = fs_2, fs_jpim1   ! vector opt.
794               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
795                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
796               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
797                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
798               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
799               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
800            END DO 
801         END DO 
802         !                                             ! ===============
803      END DO                                           !   End of slab
804      !                                                ! ===============
805   END SUBROUTINE vor_eeT
806
807
808   SUBROUTINE dyn_vor_init
809      !!---------------------------------------------------------------------
810      !!                  ***  ROUTINE dyn_vor_init  ***
811      !!
812      !! ** Purpose :   Control the consistency between cpp options for
813      !!              tracer advection schemes
814      !!----------------------------------------------------------------------
815      INTEGER ::   ji, jj, jk    ! dummy loop indices
816      INTEGER ::   ioptio, ios   ! local integer
817      !!
818      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
819         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
820      !!----------------------------------------------------------------------
821      !
822      IF(lwp) THEN
823         WRITE(numout,*)
824         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
825         WRITE(numout,*) '~~~~~~~~~~~~'
826      ENDIF
827      !
828      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
829      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
830901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
831      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
832      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
833902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
834      IF(lwm) WRITE ( numond, namdyn_vor )
835      !
836      IF(lwp) THEN                    ! Namelist print
837         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
838         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
839         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
840         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
841         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
842         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
843         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
844         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
845         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
846      ENDIF
847
848      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
849
850!!gm  this should be removed when choosing a unique strategy for fmask at the coast
851      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
852      ! at angles with three ocean points and one land point
853      IF(lwp) WRITE(numout,*)
854      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
855      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
856         DO jk = 1, jpk
857            DO jj = 1, jpjm1
858               DO ji = 1, jpim1
859                  IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
860                     & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
861               END DO
862            END DO
863         END DO
864         !
865         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
866         !
867      ENDIF
868!!gm end
869
870      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
871      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
872      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
873      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
874      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
875      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
876      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
877      !
878      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
879      !                     
880      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
881      ncor = np_COR                       ! planetary vorticity
882      SELECT CASE( n_dynadv )
883      CASE( np_LIN_dyn )
884         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
885         nrvm = np_COR        ! planetary vorticity
886         ntot = np_COR        !     -         -
887      CASE( np_VEC_c2  )
888         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
889         nrvm = np_RVO        ! relative vorticity
890         ntot = np_CRV        ! relative + planetary vorticity         
891      CASE( np_FLX_c2 , np_FLX_ubs  )
892         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
893         nrvm = np_MET        ! metric term
894         ntot = np_CME        ! Coriolis + metric term
895         !
896         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
897         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
898            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
899            DO jj = 2, jpjm1
900               DO ji = 2, jpim1
901                  di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
902                  dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
903               END DO
904            END DO
905            CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
906            !
907         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
908            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
909            DO jj = 1, jpjm1
910               DO ji = 1, jpim1
911                  di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
912                  dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
913               END DO
914            END DO
915            CALL lbc_lnk_multi( di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
916         END SELECT
917         !
918      END SELECT
919     
920      IF(lwp) THEN                   ! Print the choice
921         WRITE(numout,*)
922         SELECT CASE( nvor_scheme )
923         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
924         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
925         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
926         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
927         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
928         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
929         END SELECT         
930      ENDIF
931      !
932   END SUBROUTINE dyn_vor_init
933
934   !!==============================================================================
935END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.