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/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code

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