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/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 5974

Last change on this file since 5974 was 5974, checked in by timgraham, 8 years ago

Upgrade to head of trunk (r5936)

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