source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

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