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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 9190

Last change on this file since 9190 was 9190, checked in by gm, 6 years ago

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