source: NEMO/trunk/src/OCE/DYN/dynvor.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 7 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 52.3 KB
RevLine 
[3]1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
[2715]7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
[9528]8   !!            5.0  ! 1991-11  (G. Madec)  vor_ene, vor_mix: Original code
[2715]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
[9019]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
[7646]19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
[9019]20   !!            4.0  ! 2017-07  (G. Madec)  linear dynamics + trends diag. with Stokes-Coriolis
[9528]21   !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
22   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
[503]23   !!----------------------------------------------------------------------
[3]24
25   !!----------------------------------------------------------------------
[9019]26   !!   dyn_vor       : Update the momentum trend with the vorticity trend
27   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
28   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
29   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
30   !!   dyn_vor_init  : set and control of the different vorticity option
[3]31   !!----------------------------------------------------------------------
[503]32   USE oce            ! ocean dynamics and tracers
33   USE dom_oce        ! ocean space and time domain
[3294]34   USE dommsk         ! ocean mask
[9019]35   USE dynadv         ! momentum advection
[4990]36   USE trd_oce        ! trends: ocean variables
37   USE trddyn         ! trend manager: dynamics
[7646]38   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
39   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
[5836]40   !
[503]41   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl         ! Print control
43   USE in_out_manager ! I/O manager
[3294]44   USE lib_mpp        ! MPP library
45   USE timing         ! Timing
[3]46
47   IMPLICIT NONE
48   PRIVATE
49
[2528]50   PUBLIC   dyn_vor        ! routine called by step.F90
[5836]51   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
[3]52
[4147]53   !                                   !!* Namelist namdyn_vor: vorticity term
[9528]54   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
55   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
56   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
57   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
58   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
[5836]59   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
[9528]60   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
[5836]61   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
[3]62
[9528]63   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
64   !                                       ! associated indices:
65   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
[5836]66   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
[9528]67   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
68   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
[5836]69   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
[9528]70   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
[455]71
[5836]72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
73   !                               ! associated indices:
[9528]74   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
75   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
76   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
77   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
78   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
79
80   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
82   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
[5836]84   
85   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
86   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
87   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
88   
[3]89   !! * Substitutions
[12377]90#  include "do_loop_substitute.h90"
[3]91   !!----------------------------------------------------------------------
[9598]92   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]93   !! $Id$
[10068]94   !! Software governed by the CeCILL license (see ./LICENSE)
[3]95   !!----------------------------------------------------------------------
96CONTAINS
97
[12377]98   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
[3]99      !!----------------------------------------------------------------------
100      !!
[455]101      !! ** Purpose :   compute the lateral ocean tracer physics.
102      !!
[12377]103      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
[503]104      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]105      !!               and planetary vorticity trends) and send them to trd_dyn
106      !!               for futher diagnostics (l_trddyn=T)
[503]107      !!----------------------------------------------------------------------
[12377]108      INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index
109      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices
110      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
[2715]111      !
[9019]112      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
[455]113      !!----------------------------------------------------------------------
[2715]114      !
[9019]115      IF( ln_timing )   CALL timing_start('dyn_vor')
[3294]116      !
[9019]117      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
118         !
119         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
120         !
[12377]121         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
122         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]123         SELECT CASE( nvor_scheme )
[12377]124         CASE( np_ENS )           ;   CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
125            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
126         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
127            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
128         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
129            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
130         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
131            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
132         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
133            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9019]134         END SELECT
[12377]135         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
136         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
137         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
[9019]138         !
139         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
[12377]140            ztrdu(:,:,:) = puu(:,:,:,Krhs)
141            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]142            SELECT CASE( nvor_scheme )
[12377]143            CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
144            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
145            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
146            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
147            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
[9019]148            END SELECT
[12377]149            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
150            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
151            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
[9019]152         ENDIF
153         !
154         DEALLOCATE( ztrdu, ztrdv )
155         !
156      ELSE              !==  total vorticity trend added to the general trend  ==!
157         !
158         SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
[9528]159         CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
[12377]160                             CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
161            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9528]162         CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
[12377]163                             CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
164            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9019]165         CASE( np_ENE )                        !* energy conserving scheme
[12377]166                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
167            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9019]168         CASE( np_ENS )                        !* enstrophy conserving scheme
[12377]169                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
170            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend
[9019]171         CASE( np_MIX )                        !* mixed ene-ens scheme
[12377]172                             CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
173                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
174            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9019]175         CASE( np_EEN )                        !* energy and enstrophy conserving scheme
[12377]176                             CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
177            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
[9019]178         END SELECT
[643]179         !
[9019]180      ENDIF
[2715]181      !
[455]182      !                       ! print sum trends (used for debugging)
[12377]183      IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
184         &                                tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]185      !
[9019]186      IF( ln_timing )   CALL timing_stop('dyn_vor')
[3294]187      !
[455]188   END SUBROUTINE dyn_vor
189
190
[12377]191   SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]192      !!----------------------------------------------------------------------
193      !!                  ***  ROUTINE vor_enT  ***
194      !!
195      !! ** Purpose :   Compute the now total vorticity trend and add it to
196      !!      the general trend of the momentum equation.
197      !!
198      !! ** Method  :   Trend evaluated using now fields (centered in time)
199      !!       and t-point evaluation of vorticity (planetary and relative).
200      !!       conserves the horizontal kinetic energy.
201      !!         The general trend of momentum is increased due to the vorticity
202      !!       term which is given by:
203      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
204      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
205      !!       where rvor is the relative vorticity at f-point
206      !!
[12377]207      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]208      !!----------------------------------------------------------------------
209      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
[12377]210      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
[9528]211      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
212      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
213      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
214      !
215      INTEGER  ::   ji, jj, jk           ! dummy loop indices
216      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
[10425]217      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
218      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace
[9528]219      !!----------------------------------------------------------------------
220      !
221      IF( kt == nit000 ) THEN
222         IF(lwp) WRITE(numout,*)
223         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
224         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
225      ENDIF
226      !
[10425]227      !
228      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
229      CASE ( np_RVO )                           !* relative vorticity
230         DO jk = 1, jpkm1                                 ! Horizontal slab
[12377]231            DO_2D_10_10
232               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
233                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
234            END_2D
[9528]235            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
[12377]236               DO_2D_10_10
237                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
238               END_2D
[9528]239            ENDIF
[10425]240         END DO
241
242         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
243
244      CASE ( np_CRV )                           !* Coriolis + relative vorticity
245         DO jk = 1, jpkm1                                 ! Horizontal slab
[12377]246            DO_2D_10_10
247               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
248                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
249            END_2D
[10425]250            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity
[12377]251               DO_2D_10_10
252                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
253               END_2D
[10425]254            ENDIF
255         END DO
256
257         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
258
259      END SELECT
260
261      !                                                ! ===============
262      DO jk = 1, jpkm1                                 ! Horizontal slab
263      !                                                ! ===============
264
265         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
266         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[12377]267            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
[10425]268         CASE ( np_RVO )                           !* relative vorticity
[12377]269            DO_2D_01_01
270               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
271                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
272            END_2D
[9528]273         CASE ( np_MET )                           !* metric term
[12377]274            DO_2D_01_01
275               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
276                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,Kmm)
277            END_2D
[9528]278         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[12377]279            DO_2D_01_01
280               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
281                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
282            END_2D
[9528]283         CASE ( np_CME )                           !* Coriolis + metric
[12377]284            DO_2D_01_01
285               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
286                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
287                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,Kmm)
288            END_2D
[9528]289         CASE DEFAULT                                             ! error
290            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
291         END SELECT
292         !
293         !                                   !==  compute and add the vorticity term trend  =!
[12377]294         DO_2D_00_00
295            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
296               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
297               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
298               !
299            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
300               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
301               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
302         END_2D
[9528]303         !                                             ! ===============
304      END DO                                           !   End of slab
305      !                                                ! ===============
306   END SUBROUTINE vor_enT
307
308
[12377]309   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[455]310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE vor_ene  ***
312      !!
[3]313      !! ** Purpose :   Compute the now total vorticity trend and add it to
314      !!      the general trend of the momentum equation.
315      !!
316      !! ** Method  :   Trend evaluated using now fields (centered in time)
[5836]317      !!       and the Sadourny (1975) flux form formulation : conserves the
318      !!       horizontal kinetic energy.
319      !!         The general trend of momentum is increased due to the vorticity
320      !!       term which is given by:
[12377]321      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v pvv(:,:,:,Kmm)) ]
322      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
[5836]323      !!       where rvor is the relative vorticity
[3]324      !!
[12377]325      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[3]326      !!
[503]327      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]328      !!----------------------------------------------------------------------
[9019]329      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
[12377]330      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
[9019]331      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
[12377]332      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
333      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
[2715]334      !
[5836]335      INTEGER  ::   ji, jj, jk           ! dummy loop indices
336      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
[9019]337      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
[3]338      !!----------------------------------------------------------------------
[3294]339      !
[52]340      IF( kt == nit000 ) THEN
341         IF(lwp) WRITE(numout,*)
[455]342         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
343         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]344      ENDIF
[5836]345      !
[3]346      !                                                ! ===============
347      DO jk = 1, jpkm1                                 ! Horizontal slab
348         !                                             ! ===============
[1438]349         !
[5836]350         SELECT CASE( kvor )                 !==  vorticity considered  ==!
351         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[7753]352            zwz(:,:) = ff_f(:,:) 
[5836]353         CASE ( np_RVO )                           !* relative vorticity
[12377]354            DO_2D_10_10
355               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
356                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
357            END_2D
[5836]358         CASE ( np_MET )                           !* metric term
[12377]359            DO_2D_10_10
360               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
361                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
362            END_2D
[5836]363         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[12377]364            DO_2D_10_10
365               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
366                  &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
367            END_2D
[5836]368         CASE ( np_CME )                           !* Coriolis + metric
[12377]369            DO_2D_10_10
370               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
371                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
372            END_2D
[5836]373         CASE DEFAULT                                             ! error
374            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
[455]375         END SELECT
[5836]376         !
377         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
[12377]378            DO_2D_10_10
379               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
380            END_2D
[5836]381         ENDIF
[455]382
383         IF( ln_sco ) THEN
[12377]384            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
385            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
386            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[3]387         ELSE
[12377]388            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
389            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
[3]390         ENDIF
[5836]391         !                                   !==  compute and add the vorticity term trend  =!
[12377]392         DO_2D_00_00
393            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
394            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
395            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
396            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
397            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
398            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
399         END_2D
[3]400         !                                             ! ===============
401      END DO                                           !   End of slab
402      !                                                ! ===============
[455]403   END SUBROUTINE vor_ene
[216]404
405
[12377]406   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[3]407      !!----------------------------------------------------------------------
[455]408      !!                ***  ROUTINE vor_ens  ***
[3]409      !!
410      !! ** Purpose :   Compute the now total vorticity trend and add it to
411      !!      the general trend of the momentum equation.
412      !!
413      !! ** Method  :   Trend evaluated using now fields (centered in time)
414      !!      and the Sadourny (1975) flux FORM formulation : conserves the
415      !!      potential enstrophy of a horizontally non-divergent flow. the
416      !!      trend of the vorticity term is given by:
[12377]417      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
418      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
419      !!      Add this trend to the general momentum trend:
420      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
[3]421      !!
[12377]422      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
[3]423      !!
[503]424      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]425      !!----------------------------------------------------------------------
[9019]426      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
[12377]427      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
[9019]428      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
[12377]429      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
430      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
[2715]431      !
[5836]432      INTEGER  ::   ji, jj, jk   ! dummy loop indices
433      REAL(wp) ::   zuav, zvau   ! local scalars
[9019]434      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
[3]435      !!----------------------------------------------------------------------
[3294]436      !
[52]437      IF( kt == nit000 ) THEN
438         IF(lwp) WRITE(numout,*)
[455]439         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
440         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]441      ENDIF
[3]442      !                                                ! ===============
443      DO jk = 1, jpkm1                                 ! Horizontal slab
444         !                                             ! ===============
[1438]445         !
[5836]446         SELECT CASE( kvor )                 !==  vorticity considered  ==!
447         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[7646]448            zwz(:,:) = ff_f(:,:) 
[5836]449         CASE ( np_RVO )                           !* relative vorticity
[12377]450            DO_2D_10_10
451               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
452                  &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
453            END_2D
[5836]454         CASE ( np_MET )                           !* metric term
[12377]455            DO_2D_10_10
456               zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
457                  &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
458            END_2D
[5836]459         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[12377]460            DO_2D_10_10
461               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
462                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
463            END_2D
[5836]464         CASE ( np_CME )                           !* Coriolis + metric
[12377]465            DO_2D_10_10
466               zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
467                  &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
468            END_2D
[5836]469         CASE DEFAULT                                             ! error
470            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
[455]471         END SELECT
[1438]472         !
[5836]473         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
[12377]474            DO_2D_10_10
475               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
476            END_2D
[5836]477         ENDIF
478         !
479         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
[12377]480            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
481            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
482            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[3]483         ELSE
[12377]484            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
485            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
[3]486         ENDIF
[5836]487         !                                   !==  compute and add the vorticity term trend  =!
[12377]488         DO_2D_00_00
489            zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  &
490               &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  )
491            zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
492               &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
493            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
494            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
495         END_2D
[3]496         !                                             ! ===============
497      END DO                                           !   End of slab
498      !                                                ! ===============
[455]499   END SUBROUTINE vor_ens
[216]500
501
[12377]502   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[108]503      !!----------------------------------------------------------------------
[455]504      !!                ***  ROUTINE vor_een  ***
[108]505      !!
506      !! ** Purpose :   Compute the now total vorticity trend and add it to
507      !!      the general trend of the momentum equation.
508      !!
509      !! ** Method  :   Trend evaluated using now fields (centered in time)
[1438]510      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]511      !!      both the horizontal kinetic energy and the potential enstrophy
[1438]512      !!      when horizontal divergence is zero (see the NEMO documentation)
[12377]513      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[108]514      !!
[12377]515      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[108]516      !!
[503]517      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
518      !!----------------------------------------------------------------------
[9019]519      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
[12377]520      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
[9019]521      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
[12377]522      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
523      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
[5836]524      !
525      INTEGER  ::   ji, jj, jk   ! dummy loop indices
526      INTEGER  ::   ierr         ! local integer
527      REAL(wp) ::   zua, zva     ! local scalars
[9528]528      REAL(wp) ::   zmsk, ze3f   ! local scalars
[10425]529      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f
530      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
531      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
[108]532      !!----------------------------------------------------------------------
[3294]533      !
[108]534      IF( kt == nit000 ) THEN
535         IF(lwp) WRITE(numout,*)
[455]536         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
537         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]538      ENDIF
[5836]539      !
540      !                                                ! ===============
541      DO jk = 1, jpkm1                                 ! Horizontal slab
542         !                                             ! ===============
543         !
544         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
545         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
[12377]546            DO_2D_10_10
547               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
548                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
549               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
550               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
551               ENDIF
552            END_2D
[5836]553         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
[12377]554            DO_2D_10_10
555               ze3f = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
556                  &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
557               zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
558                  &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
559               IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3f
560               ELSE                       ;   z1_e3f(ji,jj) = 0._wp
561               ENDIF
562            END_2D
[5836]563         END SELECT
564         !
565         SELECT CASE( kvor )                 !==  vorticity considered  ==!
566         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[12377]567            DO_2D_10_10
568               zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
569            END_2D
[5836]570         CASE ( np_RVO )                           !* relative vorticity
[12377]571            DO_2D_10_10
572               zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
573                  &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
574            END_2D
[5836]575         CASE ( np_MET )                           !* metric term
[12377]576            DO_2D_10_10
577               zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
578                  &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
579            END_2D
[5836]580         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[12377]581            DO_2D_10_10
582               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
583                  &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
584                  &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
585            END_2D
[5836]586         CASE ( np_CME )                           !* Coriolis + metric
[12377]587            DO_2D_10_10
588               zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
589                  &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
590            END_2D
[5836]591         CASE DEFAULT                                             ! error
592            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
[455]593         END SELECT
[5836]594         !
595         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
[12377]596            DO_2D_10_10
597               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
598            END_2D
[5836]599         ENDIF
[10425]600      END DO                                           !   End of slab
[5836]601         !
[10425]602      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
603
604      DO jk = 1, jpkm1                                 ! Horizontal slab
[5907]605         !
[5836]606         !                                   !==  horizontal fluxes  ==!
[12377]607         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
608         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[108]609
[5836]610         !                                   !==  compute and add the vorticity term trend  =!
[1438]611         jj = 2
612         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
[5836]613         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
[10425]614               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
615               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
616               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
617               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
[108]618         END DO
619         DO jj = 3, jpj
[12377]620            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
[10425]621               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
622               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
623               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
624               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
[108]625            END DO
626         END DO
[12377]627         DO_2D_00_00
628            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
629               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
630            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
631               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
632            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
633            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
634         END_2D
[108]635         !                                             ! ===============
636      END DO                                           !   End of slab
637      !                                                ! ===============
[455]638   END SUBROUTINE vor_een
[216]639
640
[9528]641
[12377]642   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]643      !!----------------------------------------------------------------------
644      !!                ***  ROUTINE vor_eeT  ***
645      !!
646      !! ** Purpose :   Compute the now total vorticity trend and add it to
647      !!      the general trend of the momentum equation.
648      !!
649      !! ** Method  :   Trend evaluated using now fields (centered in time)
650      !!      and the Arakawa and Lamb (1980) vector form formulation using
651      !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
652      !!      The change consists in
[12377]653      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[9528]654      !!
[12377]655      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]656      !!
657      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
658      !!----------------------------------------------------------------------
659      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
[12377]660      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
[9528]661      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
[12377]662      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
663      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
[9528]664      !
665      INTEGER  ::   ji, jj, jk     ! dummy loop indices
666      INTEGER  ::   ierr           ! local integer
667      REAL(wp) ::   zua, zva       ! local scalars
668      REAL(wp) ::   zmsk, z1_e3t   ! local scalars
[10425]669      REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy 
670      REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse
671      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz
[9528]672      !!----------------------------------------------------------------------
673      !
674      IF( kt == nit000 ) THEN
675         IF(lwp) WRITE(numout,*)
676         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
677         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
678      ENDIF
679      !
680      !                                                ! ===============
681      DO jk = 1, jpkm1                                 ! Horizontal slab
682         !                                             ! ===============
683         !
684         !
685         SELECT CASE( kvor )                 !==  vorticity considered  ==!
686         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
[12377]687            DO_2D_10_10
688               zwz(ji,jj,jk) = ff_f(ji,jj)
689            END_2D
[9528]690         CASE ( np_RVO )                           !* relative vorticity
[12377]691            DO_2D_10_10
692               zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
693                  &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
694                  &          * r1_e1e2f(ji,jj)
695            END_2D
[9528]696         CASE ( np_MET )                           !* metric term
[12377]697            DO_2D_10_10
698               zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
699                  &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
700            END_2D
[9528]701         CASE ( np_CRV )                           !* Coriolis + relative vorticity
[12377]702            DO_2D_10_10
703               zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
704                  &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
705                  &                         * r1_e1e2f(ji,jj)    )
706            END_2D
[9528]707         CASE ( np_CME )                           !* Coriolis + metric
[12377]708            DO_2D_10_10
709               zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
710                  &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
711            END_2D
[9528]712         CASE DEFAULT                                             ! error
713            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
714         END SELECT
715         !
716         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
[12377]717            DO_2D_10_10
718               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
719            END_2D
[9528]720         ENDIF
[10425]721      END DO
722      !
723      CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
724      !
725      DO jk = 1, jpkm1                                 ! Horizontal slab
726
727      !                                   !==  horizontal fluxes  ==!
[12377]728         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
729         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[9528]730
731         !                                   !==  compute and add the vorticity term trend  =!
732         jj = 2
733         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
734         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
[12377]735               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]736               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
737               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
738               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
739               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
[9528]740         END DO
741         DO jj = 3, jpj
[12377]742            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
743               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]744               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
745               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
746               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
747               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
[9528]748            END DO
749         END DO
[12377]750         DO_2D_00_00
751            zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
752               &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
753            zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
754               &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
755            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
756            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
757         END_2D
[9528]758         !                                             ! ===============
759      END DO                                           !   End of slab
760      !                                                ! ===============
761   END SUBROUTINE vor_eeT
762
763
[2528]764   SUBROUTINE dyn_vor_init
[3]765      !!---------------------------------------------------------------------
[2528]766      !!                  ***  ROUTINE dyn_vor_init  ***
[3]767      !!
768      !! ** Purpose :   Control the consistency between cpp options for
[1438]769      !!              tracer advection schemes
[3]770      !!----------------------------------------------------------------------
[9528]771      INTEGER ::   ji, jj, jk    ! dummy loop indices
772      INTEGER ::   ioptio, ios   ! local integer
[2715]773      !!
[9528]774      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
775         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
[3]776      !!----------------------------------------------------------------------
[9528]777      !
778      IF(lwp) THEN
779         WRITE(numout,*)
780         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
781         WRITE(numout,*) '~~~~~~~~~~~~'
782      ENDIF
783      !
[4147]784      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
[11536]785901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
[4147]786      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
[11536]787902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
[4624]788      IF(lwm) WRITE ( numond, namdyn_vor )
[9528]789      !
[503]790      IF(lwp) THEN                    ! Namelist print
[7646]791         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
792         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
[9528]793         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
794         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
795         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
[7646]796         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
797         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
[9528]798         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
[7646]799         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
[52]800      ENDIF
801
[9528]802      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
803
[5836]804!!gm  this should be removed when choosing a unique strategy for fmask at the coast
[3294]805      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
806      ! at angles with three ocean points and one land point
[5836]807      IF(lwp) WRITE(numout,*)
[7646]808      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
[3294]809      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
[12377]810         DO_3D_10_10( 1, jpk )
811            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
812               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
813         END_3D
[9528]814         !
[10425]815         CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
[9528]816         !
[3294]817      ENDIF
[5836]818!!gm end
[3294]819
[5836]820      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
[9528]821      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
822      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
823      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
824      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
825      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
826      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
[5836]827      !
[6140]828      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[5836]829      !                     
830      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
[9019]831      ncor = np_COR                       ! planetary vorticity
832      SELECT CASE( n_dynadv )
833      CASE( np_LIN_dyn )
[9190]834         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
[9019]835         nrvm = np_COR        ! planetary vorticity
836         ntot = np_COR        !     -         -
837      CASE( np_VEC_c2  )
[9190]838         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
[5836]839         nrvm = np_RVO        ! relative vorticity
[9019]840         ntot = np_CRV        ! relative + planetary vorticity         
841      CASE( np_FLX_c2 , np_FLX_ubs  )
[9190]842         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
[5836]843         nrvm = np_MET        ! metric term
844         ntot = np_CME        ! Coriolis + metric term
[9528]845         !
846         SELECT CASE( nvor_scheme )    ! pre-computed gradients for the metric term:
847         CASE( np_ENT )                      !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
848            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
[12377]849            DO_2D_00_00
850               di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
851               dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
852            END_2D
[10425]853            CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
[9528]854            !
855         CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
856            ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
[12377]857            DO_2D_10_10
858               di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
859               dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
860            END_2D
[10425]861            CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. )   ! Lateral boundary conditions
[9528]862         END SELECT
863         !
[9019]864      END SELECT
[643]865     
[503]866      IF(lwp) THEN                   ! Print the choice
867         WRITE(numout,*)
[9019]868         SELECT CASE( nvor_scheme )
[9528]869         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)'
870         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
871         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)'
872         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)'
873         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)'
874         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
[9019]875         END SELECT         
[3]876      ENDIF
[503]877      !
[2528]878   END SUBROUTINE dyn_vor_init
[3]879
[503]880   !!==============================================================================
[3]881END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.