source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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