New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynvor.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

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