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/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7422

Last change on this file since 7422 was 7422, checked in by gm, 7 years ago

#1805 dev_INGV_UKMO_2016: improve Stokes drift (including dynspg_ts , Stokes-Coriolis force , and GLS surface roughness

  • Property svn:keywords set to Id
File size: 37.5 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(:,:) 
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(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(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         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
292         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
293         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
294
295         !                                   !==  compute and add the vorticity term trend  =!
296         DO jj = 2, jpjm1
297            DO ji = fs_2, fs_jpim1   ! vector opt.
298               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
299               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
300               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
301               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
302               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
303               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
304            END DO 
305         END DO 
306         !                                             ! ===============
307      END DO                                           !   End of slab
308      !                                                ! ===============
309      CALL wrk_dealloc( jpi,jpj,   zwx, zwy, zwz ) 
310      !
311      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
312      !
313   END SUBROUTINE vor_ene
314
315
316   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
317      !!----------------------------------------------------------------------
318      !!                ***  ROUTINE vor_ens  ***
319      !!
320      !! ** Purpose :   Compute the now total vorticity trend and add it to
321      !!      the general trend of the momentum equation.
322      !!
323      !! ** Method  :   Trend evaluated using now fields (centered in time)
324      !!      and the Sadourny (1975) flux FORM formulation : conserves the
325      !!      potential enstrophy of a horizontally non-divergent flow. the
326      !!      trend of the vorticity term is given by:
327      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
328      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
329      !!      Add this trend to the general momentum trend (ua,va):
330      !!          (ua,va) = (ua,va) + ( voru , vorv )
331      !!
332      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
333      !!
334      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
335      !!----------------------------------------------------------------------
336      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
337      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
338         !                                                             ! =nrvm (relative vorticity or metric)
339      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
340      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
341      !
342      INTEGER  ::   ji, jj, jk   ! dummy loop indices
343      REAL(wp) ::   zuav, zvau   ! local scalars
344      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
345      !!----------------------------------------------------------------------
346      !
347      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
348      !
349      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zwz ) 
350      !
351      IF( kt == nit000 ) THEN
352         IF(lwp) WRITE(numout,*)
353         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
354         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
355      ENDIF
356      !                                                ! ===============
357      DO jk = 1, jpkm1                                 ! Horizontal slab
358         !                                             ! ===============
359         !
360         SELECT CASE( kvor )                 !==  vorticity considered  ==!
361         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
362            zwz(:,:) = ff(:,:) 
363         CASE ( np_RVO )                           !* relative vorticity
364            DO jj = 1, jpjm1
365               DO ji = 1, fs_jpim1   ! vector opt.
366                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
367                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
368               END DO
369            END DO
370         CASE ( np_MET )                           !* metric term
371            DO jj = 1, jpjm1
372               DO ji = 1, fs_jpim1   ! vector opt.
373                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
374                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
375                       &     * 0.5 * r1_e1e2f(ji,jj)
376               END DO
377            END DO
378         CASE ( np_CRV )                           !* Coriolis + relative vorticity
379            DO jj = 1, jpjm1
380               DO ji = 1, fs_jpim1   ! vector opt.
381                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
382                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
383                     &                   * r1_e1e2f(ji,jj)
384               END DO
385            END DO
386         CASE ( np_CME )                           !* Coriolis + metric
387            DO jj = 1, jpjm1
388               DO ji = 1, fs_jpim1   ! vector opt.
389                  zwz(ji,jj) = ff(ji,jj)                                                                       &
390                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
391                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
392                       &     * 0.5 * r1_e1e2f(ji,jj)
393               END DO
394            END DO
395         CASE DEFAULT                                             ! error
396            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
397         END SELECT
398         !
399         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
400            DO jj = 1, jpjm1
401               DO ji = 1, fs_jpim1   ! vector opt.
402                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
403               END DO
404            END DO
405         ENDIF
406         !
407         !                                 !==  horizontal fluxes  ==!
408         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
409         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
410         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
411         !
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, pun, pvn, pua, pva )
434      !!----------------------------------------------------------------------
435      !!                ***  ROUTINE vor_een  ***
436      !!
437      !! ** Purpose :   Compute the now total vorticity trend and add it to
438      !!      the general trend of the momentum equation.
439      !!
440      !! ** Method  :   Trend evaluated using now fields (centered in time)
441      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
442      !!      both the horizontal kinetic energy and the potential enstrophy
443      !!      when horizontal divergence is zero (see the NEMO documentation)
444      !!      Add this trend to the general momentum trend (ua,va).
445      !!
446      !! ** Action : - Update (ua,va) with the now vorticity term trend
447      !!
448      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
449      !!----------------------------------------------------------------------
450      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index
451      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;
452         !                                                             ! =nrvm (relative vorticity or metric)
453      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities
454      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend
455      !
456      INTEGER  ::   ji, jj, jk   ! dummy loop indices
457      INTEGER  ::   ierr         ! local integer
458      REAL(wp) ::   zua, zva     ! local scalars
459      REAL(wp) ::   zmsk, ze3    ! local scalars
460      !
461      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f
462      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse
463      !!----------------------------------------------------------------------
464      !
465      IF( nn_timing == 1 )  CALL timing_start('vor_een')
466      !
467      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
468      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
469      !
470      IF( kt == nit000 ) THEN
471         IF(lwp) WRITE(numout,*)
472         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
473         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
474      ENDIF
475      !
476      !                                                ! ===============
477      DO jk = 1, jpkm1                                 ! Horizontal slab
478         !                                             ! ===============
479         !
480         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
481         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
482            DO jj = 1, jpjm1
483               DO ji = 1, fs_jpim1   ! vector opt.
484                  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)   &
485                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
486                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3
487                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp
488                  ENDIF
489               END DO
490            END DO
491         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
492            DO jj = 1, jpjm1
493               DO ji = 1, fs_jpim1   ! vector opt.
494                  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)   &
495                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
496                  zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
497                     &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
498                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3
499                  ELSE                      ;   z1_e3f(ji,jj) = 0._wp
500                  ENDIF
501               END DO
502            END DO
503         END SELECT
504         !
505         SELECT CASE( kvor )                 !==  vorticity considered  ==!
506         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
507            DO jj = 1, jpjm1
508               DO ji = 1, fs_jpim1   ! vector opt.
509                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
510               END DO
511            END DO
512         CASE ( np_RVO )                           !* relative vorticity
513            DO jj = 1, jpjm1
514               DO ji = 1, fs_jpim1   ! vector opt.
515                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
516                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
517                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
518               END DO
519            END DO
520         CASE ( np_MET )                           !* metric term
521            DO jj = 1, jpjm1
522               DO ji = 1, fs_jpim1   ! vector opt.
523                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
524                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
525                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
526               END DO
527            END DO
528         CASE ( np_CRV )                           !* Coriolis + relative vorticity
529            DO jj = 1, jpjm1
530               DO ji = 1, fs_jpim1   ! vector opt.
531                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
532                     &                         - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
533                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
534               END DO
535            END DO
536         CASE ( np_CME )                           !* Coriolis + metric
537            DO jj = 1, jpjm1
538               DO ji = 1, fs_jpim1   ! vector opt.
539                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
540                       &        + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
541                       &            - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
542                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
543               END DO
544            END DO
545         CASE DEFAULT                                             ! error
546            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
547         END SELECT
548         !
549         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
550            DO jj = 1, jpjm1
551               DO ji = 1, fs_jpim1   ! vector opt.
552                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
553               END DO
554            END DO
555         ENDIF
556         !
557         CALL lbc_lnk( zwz, 'F', 1. )
558         !
559         !                                   !==  horizontal fluxes  ==!
560         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
561         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
562
563         !                                   !==  compute and add the vorticity term trend  =!
564         jj = 2
565         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
566         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
567               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
568               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
569               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
570               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
571         END DO
572         DO jj = 3, jpj
573            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
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         END DO
580         DO jj = 2, jpjm1
581            DO ji = fs_2, fs_jpim1   ! vector opt.
582               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
583                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
584               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
585                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
586               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
587               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
588            END DO 
589         END DO 
590         !                                             ! ===============
591      END DO                                           !   End of slab
592      !                                                ! ===============
593      !
594      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
595      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
596      !
597      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
598      !
599   END SUBROUTINE vor_een
600
601
602   SUBROUTINE dyn_vor_init
603      !!---------------------------------------------------------------------
604      !!                  ***  ROUTINE dyn_vor_init  ***
605      !!
606      !! ** Purpose :   Control the consistency between cpp options for
607      !!              tracer advection schemes
608      !!----------------------------------------------------------------------
609      INTEGER ::   ioptio          ! local integer
610      INTEGER ::   ji, jj, jk      ! dummy loop indices
611      INTEGER ::   ios             ! Local integer output status for namelist read
612      !!
613      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
614      !!----------------------------------------------------------------------
615
616      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
617      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
618901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
619
620      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
621      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
622902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
623      IF(lwm) WRITE ( numond, namdyn_vor )
624
625      IF(lwp) THEN                    ! Namelist print
626         WRITE(numout,*)
627         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
628         WRITE(numout,*) '~~~~~~~~~~~~'
629         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
630         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
631         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
632         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
633         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
634         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
635         WRITE(numout,*) '           masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
636      ENDIF
637
638!!gm  this should be removed when choosing a unique strategy for fmask at the coast
639      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
640      ! at angles with three ocean points and one land point
641      IF(lwp) WRITE(numout,*)
642      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
643      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
644         DO jk = 1, jpk
645            DO jj = 2, jpjm1
646               DO ji = 2, jpim1
647                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
648                      fmask(ji,jj,jk) = 1._wp
649               END DO
650            END DO
651         END DO
652          !
653          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
654          !
655      ENDIF
656!!gm end
657
658      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
659      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF
660      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF
661      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF
662      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF
663      !
664      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
665      !                     
666      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
667      ncor = np_COR
668      IF( ln_dynadv_vec ) THEN     
669         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
670         nrvm = np_RVO        ! relative vorticity
671         ntot = np_CRV        ! relative + planetary vorticity
672      ELSE                       
673         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
674         nrvm = np_MET        ! metric term
675         ntot = np_CME        ! Coriolis + metric term
676      ENDIF
677     
678      IF(lwp) THEN                   ! Print the choice
679         WRITE(numout,*)
680         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
681         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
682         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
683         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
684      ENDIF
685      !
686   END SUBROUTINE dyn_vor_init
687
688   !!==============================================================================
689END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.