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 NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90 @ 13005

Last change on this file since 13005 was 13005, checked in by gm, 4 years ago
ADE and more options to AM98 config
Property svn:keywords set to `Id`
File size: 71.2 KB

Rev	Line
[3]	1	MODULE 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
[12667]	23	!! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity
[13005]	24	!! 4.x ! 2020-03 (G. Madec, A. Nasser) alternate direction computation of vorticity tendancy
	25	!! ! for ENS, ENE
[503]	26	!!----------------------------------------------------------------------
[3]	27
	28	!!----------------------------------------------------------------------
[9019]	29	!! dyn_vor : Update the momentum trend with the vorticity trend
	30	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
	31	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
	32	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
	33	!! dyn_vor_init : set and control of the different vorticity option
[3]	34	!!----------------------------------------------------------------------
[503]	35	USE oce ! ocean dynamics and tracers
	36	USE dom_oce ! ocean space and time domain
[3294]	37	USE dommsk ! ocean mask
[9019]	38	USE dynadv ! momentum advection
[4990]	39	USE trd_oce ! trends: ocean variables
	40	USE trddyn ! trend manager: dynamics
[7646]	41	USE sbcwave ! Surface Waves (add Stokes-Coriolis force)
	42	USE sbc_oce , ONLY : ln_stcor ! use Stoke-Coriolis force
[5836]	43	!
[503]	44	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	45	USE prtctl ! Print control
	46	USE in_out_manager ! I/O manager
[3294]	47	USE lib_mpp ! MPP library
	48	USE timing ! Timing
[13005]	49	!!anAD only
	50	USE dynkeg, ONLY : dyn_kegAD
	51	USE dynadv, ONLY : nn_dynkeg
[3]	52
	53	IMPLICIT NONE
	54	PRIVATE
	55
[2528]	56	PUBLIC dyn_vor ! routine called by step.F90
[5836]	57	PUBLIC dyn_vor_init ! routine called by nemogcm.F90
[3]	58
[4147]	59	! !!* Namelist namdyn_vor: vorticity term
[13005]	60	LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS)
	61	LOGICAL, PUBLIC :: ln_dynvor_ens_adVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad)
	62	LOGICAL, PUBLIC :: ln_dynvor_ens_adKE = .FALSE. !: AD enstrophy conserving scheme (ENS_ad)
	63	LOGICAL, PUBLIC :: ln_dynvor_ens_adKEVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad)
	64	LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE)
	65	LOGICAL, PUBLIC :: ln_dynvor_ene_adVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad)
	66	LOGICAL, PUBLIC :: ln_dynvor_ene_adKE = .FALSE. !: f-point AD energy conserving scheme (ENE_ad)
	67	LOGICAL, PUBLIC :: ln_dynvor_ene_adKEVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad)
	68	LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT)
	69	LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET)
	70	LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN)
	71	INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
	72	LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX)
	73	LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
[3]	74
[9528]	75	INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme
	76	! ! associated indices:
	77	INTEGER, PUBLIC, PARAMETER :: np_ENS = 0 ! ENS scheme
[5836]	78	INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme
[9528]	79	INTEGER, PUBLIC, PARAMETER :: np_ENT = 2 ! ENT scheme (t-point vorticity)
	80	INTEGER, PUBLIC, PARAMETER :: np_EET = 3 ! EET scheme (EEN using e3t)
[5836]	81	INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme
[9528]	82	INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme
[13005]	83	!!an
	84	INTEGER, PUBLIC, PARAMETER :: np_ENS_adKE = 11 ! ENS scheme - AD scheme (KE only)
	85	INTEGER, PUBLIC, PARAMETER :: np_ENS_adVO = 12 ! ENS scheme - AD scheme (VOR only)
	86	INTEGER, PUBLIC, PARAMETER :: np_ENS_adKEVO = 13 ! ENS scheme - AD scheme (KE+VOR)
	87	INTEGER, PUBLIC, PARAMETER :: np_ENE_adKE = 21 ! ENE scheme - AD scheme (KE only)
	88	INTEGER, PUBLIC, PARAMETER :: np_ENE_adVO = 22 ! ENE scheme - AD scheme (VOR only)
	89	INTEGER, PUBLIC, PARAMETER :: np_ENE_adKEVO = 23 ! ENE scheme - AD scheme (KE+VOR)
	90	!!an
	91
	92	!!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO
	93	INTEGER, PUBLIC :: ncor, nrvm, ntot ! choice of calculated vorticity
[5836]	94	! ! associated indices:
[9528]	95	INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary)
	96	INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity
	97	INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term
	98	INTEGER, PUBLIC, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity)
	99	INTEGER, PUBLIC, PARAMETER :: np_CME = 5 ! Coriolis + metric term
	100
	101	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation
	102	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - -
	103	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2u)/(2*e1e2f) used in F-point metric term calculation
	104	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1v)/(2*e1e2f) - - - -
[5836]	105
	106	REAL(wp) :: r1_4 = 0.250_wp ! =1/4
	107	REAL(wp) :: r1_8 = 0.125_wp ! =1/8
	108	REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12
	109
[3]	110	!! * Substitutions
[12377]	111	# include "do_loop_substitute.h90"
[3]	112	!!----------------------------------------------------------------------
[9598]	113	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	114	!! $Id$
[10068]	115	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	116	!!----------------------------------------------------------------------
	117	CONTAINS
	118
[13005]	119	SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs)
[3]	120	!!----------------------------------------------------------------------
	121	!!
[455]	122	!! ** Purpose : compute the lateral ocean tracer physics.
	123	!!
[12377]	124	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
[503]	125	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]	126	!! and planetary vorticity trends) and send them to trd_dyn
	127	!! for futher diagnostics (l_trddyn=T)
[503]	128	!!----------------------------------------------------------------------
[13005]	129	INTEGER :: ji, jj, jk ! dummy loop indice
	130	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	131	INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices
	132	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
[2715]	133	!
[9019]	134	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
[13005]	135	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv
[455]	136	!!----------------------------------------------------------------------
[2715]	137	!
[9019]	138	IF( ln_timing ) CALL timing_start('dyn_vor')
[3294]	139	!
[9019]	140	IF( l_trddyn ) THEN !== trend diagnostics case : split the added trend in two parts ==!
	141	!
	142	ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
	143	!
[12377]	144	ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend (including Stokes-Coriolis force)
	145	ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]	146	SELECT CASE( nvor_scheme )
[12377]	147	CASE( np_ENS ) ; CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme
	148	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	149	CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme
	150	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	151	CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts)
	152	IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	153	CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t)
	154	IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	155	CASE( np_EEN ) ; CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme
	156	IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	157	END SELECT
[12377]	158	ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
	159	ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
	160	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
[9019]	161	!
	162	IF( n_dynadv /= np_LIN_dyn ) THEN !* relative vorticity or metric trend (only in non-linear case)
[12377]	163	ztrdu(:,:,:) = puu(:,:,:,Krhs)
	164	ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]	165	SELECT CASE( nvor_scheme )
[12377]	166	CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts)
	167	CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t)
	168	CASE( np_ENE ) ; CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme
	169	CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme
	170	CASE( np_EEN ) ; CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme
[9019]	171	END SELECT
[12377]	172	ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
	173	ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
	174	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
[9019]	175	ENDIF
	176	!
	177	DEALLOCATE( ztrdu, ztrdv )
	178	!
	179	ELSE !== total vorticity trend added to the general trend ==!
	180	!
	181	SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==!
[9528]	182	CASE( np_ENT ) !* energy conserving scheme (T-pts)
[12377]	183	CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	184	IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9528]	185	CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
[12377]	186	CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	187	IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	188	CASE( np_ENE ) !* energy conserving scheme
[13005]	189	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
	190	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	191	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	192
	193	CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction
	194	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	195
	196	!== Alternative direction - VOR only ==!
	197
	198	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	199
	200	ALLOCATE( zuu(jpi,jpj,jpk) )
	201	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	202	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	203	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	204	CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	205	DEALLOCATE( zuu )
	206	ELSE
	207	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	208
	209	ALLOCATE( zvv(jpi,jpj,jpk) )
	210	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	211	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	212	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	213	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	214	DEALLOCATE( zvv )
	215	ENDIF
	216	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	217	CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction
	218	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	219
	220	!== Alternative direction - KEG only ==!
	221
[12377]	222	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
[13005]	223
	224	ALLOCATE( zuu(jpi,jpj,jpk) )
	225	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	226	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	227	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	228	CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	229	DEALLOCATE( zuu )
	230	ELSE
	231
	232	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	233
	234	ALLOCATE( zvv(jpi,jpj,jpk) )
	235	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	236	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	237	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	238	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	239	DEALLOCATE( zvv )
	240	ENDIF
[12377]	241	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[13005]	242
	243	CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction
	244	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	245
	246	!== Alternative direction - KE + VOR ==!
	247
	248	ALLOCATE( zuu(jpi,jpj,jpk) )
	249	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	250	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) !
	251	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	252	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	253	CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	254	CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
	255	DEALLOCATE( zuu )
	256	ELSE
	257
	258	ALLOCATE( zvv(jpi,jpj,jpk) )
	259	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	260	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
	261	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	262	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	263	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	264	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )
	265	DEALLOCATE( zvv )
	266	ENDIF
[9019]	267	CASE( np_ENS ) !* enstrophy conserving scheme
[13005]	268	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
	269	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	270	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	271	CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction
	272	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	273
	274	!== Alternative direction - VOR only ==!
	275
	276	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
	277
	278	ALLOCATE( zuu(jpi,jpj,jpk) )
	279	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	280	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	281	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	282	CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	283	DEALLOCATE( zuu )
	284	ELSE
	285	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )
	286
	287	ALLOCATE( zvv(jpi,jpj,jpk) )
	288	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	289	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	290	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	291	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	292	DEALLOCATE( zvv )
	293	ENDIF
	294	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	295	CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction
	296	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	297
	298	!== Alternative direction - KEG only ==!
	299
[12377]	300	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
[13005]	301
	302	ALLOCATE( zuu(jpi,jpj,jpk) )
	303	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	304	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	305	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	306	CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	307	DEALLOCATE( zuu )
	308	ELSE
	309
	310	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	311
	312	ALLOCATE( zvv(jpi,jpj,jpk) )
	313	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	314	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	315	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	316	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	317	DEALLOCATE( zvv )
	318	ENDIF
	319	CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction
	320	IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components
	321
	322	!== Alternative direction - KE + VOR ==!
	323
	324	ALLOCATE( zuu(jpi,jpj,jpk) )
	325	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	326	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) !
	327	zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)
	328	CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )
	329	CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	330	CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
	331	DEALLOCATE( zuu )
	332	ELSE
	333
	334	ALLOCATE( zvv(jpi,jpj,jpk) )
	335	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend
	336	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )
	337	zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)
	338	CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )
	339	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend
	340	CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )
	341	DEALLOCATE( zvv )
	342	ENDIF
[12377]	343	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	344	CASE( np_MIX ) !* mixed ene-ens scheme
[12377]	345	CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens)
	346	CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! planetary vorticity trend (ene)
	347	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	348	CASE( np_EEN ) !* energy and enstrophy conserving scheme
[12377]	349	CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	350	IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	351	END SELECT
[643]	352	!
[9019]	353	ENDIF
[2715]	354	!
[455]	355	! ! print sum trends (used for debugging)
[12377]	356	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
	357	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]	358	!
[9019]	359	IF( ln_timing ) CALL timing_stop('dyn_vor')
[3294]	360	!
[455]	361	END SUBROUTINE dyn_vor
	362
	363
[12377]	364	SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]	365	!!----------------------------------------------------------------------
	366	!! * ROUTINE vor_enT *
	367	!!
	368	!! ** Purpose : Compute the now total vorticity trend and add it to
	369	!! the general trend of the momentum equation.
	370	!!
	371	!! ** Method : Trend evaluated using now fields (centered in time)
	372	!! and t-point evaluation of vorticity (planetary and relative).
	373	!! conserves the horizontal kinetic energy.
	374	!! The general trend of momentum is increased due to the vorticity
	375	!! term which is given by:
	376	!! voru = 1/bu mj[ ( mi(mj(bfrvor))+btf_t)/e3t mj[vn] ]
	377	!! vorv = 1/bv mi[ ( mi(mj(bfrvor))+btf_t)/e3f mj[un] ]
	378	!! where rvor is the relative vorticity at f-point
	379	!!
[12377]	380	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]	381	!!----------------------------------------------------------------------
	382	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	383	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9528]	384	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
	385	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	386	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
	387	!
	388	INTEGER :: ji, jj, jk ! dummy loop indices
	389	REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
[10425]	390	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace
	391	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz ! 3D workspace
[9528]	392	!!----------------------------------------------------------------------
	393	!
	394	IF( kt == nit000 ) THEN
	395	IF(lwp) WRITE(numout,*)
	396	IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
	397	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	398	ENDIF
	399	!
[10425]	400	!
[12667]	401	SELECT CASE( kvor ) !== relative vorticity considered ==!
	402	CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
	403	DO jk = 1, jpkm1 ! Horizontal slab
[12377]	404	DO_2D_10_10
	405	zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	406	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	407	END_2D
[12667]	408	IF( ln_dynvor_msk ) THEN ! mask relative vorticity
[12377]	409	DO_2D_10_10
	410	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	411	END_2D
[9528]	412	ENDIF
[10425]	413	END DO
	414	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
[12667]	415	!
[10425]	416	END SELECT
	417
	418	! ! ===============
	419	DO jk = 1, jpkm1 ! Horizontal slab
[12667]	420	! ! ===============
	421	!
[10425]	422	SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
[12667]	423	!
[10425]	424	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[12377]	425	zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
[10425]	426	CASE ( np_RVO ) !* relative vorticity
[12377]	427	DO_2D_01_01
	428	zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
	429	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
	430	END_2D
[9528]	431	CASE ( np_MET ) !* metric term
[12377]	432	DO_2D_01_01
	433	zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
	434	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm)
	435	END_2D
[9528]	436	CASE ( np_CRV ) !* Coriolis + relative vorticity
[12377]	437	DO_2D_01_01
	438	zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
	439	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
	440	END_2D
[9528]	441	CASE ( np_CME ) !* Coriolis + metric
[12377]	442	DO_2D_01_01
	443	zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
	444	& + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
	445	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm)
	446	END_2D
[9528]	447	CASE DEFAULT ! error
	448	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
	449	END SELECT
	450	!
	451	! !== compute and add the vorticity term trend =!
[12377]	452	DO_2D_00_00
	453	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
	454	& * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
	455	& + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
	456	!
	457	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
	458	& * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
	459	& + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
	460	END_2D
[9528]	461	! ! ===============
	462	END DO ! End of slab
	463	! ! ===============
	464	END SUBROUTINE vor_enT
	465
	466
[12377]	467	SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[455]	468	!!----------------------------------------------------------------------
	469	!! * ROUTINE vor_ene *
	470	!!
[3]	471	!! ** Purpose : Compute the now total vorticity trend and add it to
	472	!! the general trend of the momentum equation.
	473	!!
	474	!! ** Method : Trend evaluated using now fields (centered in time)
[5836]	475	!! and the Sadourny (1975) flux form formulation : conserves the
	476	!! horizontal kinetic energy.
	477	!! The general trend of momentum is increased due to the vorticity
	478	!! term which is given by:
[12377]	479	!! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	480	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u puu(:,:,:,Kmm)) ]
[5836]	481	!! where rvor is the relative vorticity
[3]	482	!!
[12377]	483	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[3]	484	!!
[503]	485	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	486	!!----------------------------------------------------------------------
[13005]	487	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	488	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
	489	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
	490	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities
	491	REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	492	!
[5836]	493	INTEGER :: ji, jj, jk ! dummy loop indices
	494	REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
[9019]	495	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace
[3]	496	!!----------------------------------------------------------------------
[3294]	497	!
[52]	498	IF( kt == nit000 ) THEN
	499	IF(lwp) WRITE(numout,*)
[455]	500	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	501	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	502	ENDIF
[5836]	503	!
[3]	504	! ! ===============
	505	DO jk = 1, jpkm1 ! Horizontal slab
	506	! ! ===============
[1438]	507	!
[5836]	508	SELECT CASE( kvor ) !== vorticity considered ==!
	509	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[7753]	510	zwz(:,:) = ff_f(:,:)
[5836]	511	CASE ( np_RVO ) !* relative vorticity
[12377]	512	DO_2D_10_10
	513	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	514	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	515	END_2D
[12667]	516	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	517	DO_2D_10_10
	518	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	519	END_2D
	520	ENDIF
[5836]	521	CASE ( np_MET ) !* metric term
[12377]	522	DO_2D_10_10
	523	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	524	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	525	END_2D
[5836]	526	CASE ( np_CRV ) !* Coriolis + relative vorticity
[12377]	527	DO_2D_10_10
	528	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	529	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	530	END_2D
[12667]	531	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity (NOT the Coriolis term)
	532	DO_2D_10_10
	533	zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
	534	END_2D
	535	ENDIF
[5836]	536	CASE ( np_CME ) !* Coriolis + metric
[12377]	537	DO_2D_10_10
	538	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	539	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	540	END_2D
[5836]	541	CASE DEFAULT ! error
	542	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	543	END SELECT
[5836]	544	!
[13005]	545	IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !* NO alternating direction *!
	546	!
	547	! !== horizontal fluxes and potential vorticity ==!
	548	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	549	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
	550	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	551	!
	552	! !== compute and add the vorticity term trend =!
	553	DO_2D_00_00
	554	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	555	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	556	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	557	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
	558	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 )
	559	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 )
	560	END_2D
	561	!
	562	!
	563	ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !* Alternating direction : i-component *!
	564	!
	565	!
	566	! !== horizontal fluxes and potential vorticity ==!
	567	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
	568	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	569	!
	570	! !== compute and add the vorticity term trend =!
	571	DO_2D_00_00
	572	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	573	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	574	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 )
	575	END_2D
	576	!
	577	ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !* Alternating direction : j-component *!
	578	!
	579	!
	580	! !== horizontal fluxes and potential vorticity ==!
	581	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	582	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	583	!
	584	! !== compute and add the vorticity term trend =!
	585	DO_2D_00_00
	586	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	587	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
	588	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 )
	589	END_2D
	590	!
	591	ENDIF
[3]	592	! ! ===============
	593	END DO ! End of slab
	594	! ! ===============
[455]	595	END SUBROUTINE vor_ene
[216]	596
	597
[12377]	598	SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[3]	599	!!----------------------------------------------------------------------
[455]	600	!! * ROUTINE vor_ens *
[3]	601	!!
	602	!! ** Purpose : Compute the now total vorticity trend and add it to
	603	!! the general trend of the momentum equation.
	604	!!
	605	!! ** Method : Trend evaluated using now fields (centered in time)
	606	!! and the Sadourny (1975) flux FORM formulation : conserves the
	607	!! potential enstrophy of a horizontally non-divergent flow. the
	608	!! trend of the vorticity term is given by:
[12377]	609	!! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	610	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
	611	!! Add this trend to the general momentum trend:
	612	!! (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
[3]	613	!!
[12377]	614	!! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
[3]	615	!!
[503]	616	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	617	!!----------------------------------------------------------------------
[13005]	618	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	619	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
	620	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
	621	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities
	622	REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	623	!
[5836]	624	INTEGER :: ji, jj, jk ! dummy loop indices
	625	REAL(wp) :: zuav, zvau ! local scalars
[9019]	626	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! 2D workspace
[3]	627	!!----------------------------------------------------------------------
[3294]	628	!
[52]	629	IF( kt == nit000 ) THEN
	630	IF(lwp) WRITE(numout,*)
[455]	631	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	632	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	633	ENDIF
[3]	634	! ! ===============
	635	DO jk = 1, jpkm1 ! Horizontal slab
	636	! ! ===============
[1438]	637	!
[5836]	638	SELECT CASE( kvor ) !== vorticity considered ==!
	639	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[7646]	640	zwz(:,:) = ff_f(:,:)
[5836]	641	CASE ( np_RVO ) !* relative vorticity
[12377]	642	DO_2D_10_10
	643	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	644	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	645	END_2D
[12667]	646	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	647	DO_2D_10_10
	648	zwz(ji,jj) = ff_f(ji,jj) * fmask(ji,jj,jk)
	649	END_2D
	650	ENDIF
[5836]	651	CASE ( np_MET ) !* metric term
[12377]	652	DO_2D_10_10
	653	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	654	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	655	END_2D
[5836]	656	CASE ( np_CRV ) !* Coriolis + relative vorticity
[12377]	657	DO_2D_10_10
	658	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	659	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	660	END_2D
[12667]	661	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity (NOT the Coriolis term)
	662	DO_2D_10_10
	663	zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
	664	END_2D
	665	ENDIF
[5836]	666	CASE ( np_CME ) !* Coriolis + metric
[12377]	667	DO_2D_10_10
	668	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	669	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	670	END_2D
[5836]	671	CASE DEFAULT ! error
	672	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	673	END SELECT
[1438]	674	!
[5836]	675	!
[13005]	676	!!an wut ? v et u
	677	IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !* NO alternating direction *!
	678	!
	679	! !== horizontal fluxes and potential vorticity ==!
	680	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	681	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
	682	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	683	!
	684	! !== compute and add the vorticity term trend =!
	685	DO_2D_00_00
	686	zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	687	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
	688	zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	689	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
	690	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	691	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	692	END_2D
	693	!
	694	ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !* Alternating direction : i-component *!
	695	!
	696	!
	697	! !== horizontal fluxes and potential vorticity ==!
	698	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
	699	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	700	!
	701	! !== compute and add the vorticity term trend =!
	702	DO_2D_00_00
	703	zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	704	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
	705	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	706	END_2D
	707	!
	708	ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !* Alternating direction : j-component *!
	709	!
	710	!
	711	! !== horizontal fluxes and potential vorticity ==!
	712	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	713	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	714	!
	715	! !== compute and add the vorticity term trend =!
	716	DO_2D_00_00
	717	zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	718	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
	719	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	720	END_2D
	721	!
	722	ENDIF
[3]	723	! ! ===============
	724	END DO ! End of slab
	725	! ! ===============
[455]	726	END SUBROUTINE vor_ens
[216]	727
	728
[12377]	729	SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[108]	730	!!----------------------------------------------------------------------
[455]	731	!! * ROUTINE vor_een *
[108]	732	!!
	733	!! ** Purpose : Compute the now total vorticity trend and add it to
	734	!! the general trend of the momentum equation.
	735	!!
	736	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	737	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	738	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	739	!! when horizontal divergence is zero (see the NEMO documentation)
[12377]	740	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[108]	741	!!
[12377]	742	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[108]	743	!!
[503]	744	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	745	!!----------------------------------------------------------------------
[9019]	746	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	747	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	748	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	749	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	750	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[5836]	751	!
	752	INTEGER :: ji, jj, jk ! dummy loop indices
	753	INTEGER :: ierr ! local integer
	754	REAL(wp) :: zua, zva ! local scalars
[9528]	755	REAL(wp) :: zmsk, ze3f ! local scalars
[10425]	756	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f
	757	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	758	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz
[108]	759	!!----------------------------------------------------------------------
[3294]	760	!
[108]	761	IF( kt == nit000 ) THEN
	762	IF(lwp) WRITE(numout,*)
[455]	763	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	764	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]	765	ENDIF
[5836]	766	!
	767	! ! ===============
	768	DO jk = 1, jpkm1 ! Horizontal slab
	769	! ! ===============
	770	!
	771	SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point
	772	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[12377]	773	DO_2D_10_10
	774	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) &
	775	& + e3t(ji,jj ,jk,Kmm)tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)tmask(ji+1,jj ,jk) )
	776	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
	777	ELSE ; z1_e3f(ji,jj) = 0._wp
	778	ENDIF
	779	END_2D
[5836]	780	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[12377]	781	DO_2D_10_10
	782	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) &
	783	& + e3t(ji,jj ,jk,Kmm)tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)tmask(ji+1,jj ,jk) )
	784	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	785	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	786	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
	787	ELSE ; z1_e3f(ji,jj) = 0._wp
	788	ENDIF
	789	END_2D
[5836]	790	END SELECT
	791	!
	792	SELECT CASE( kvor ) !== vorticity considered ==!
	793	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[12377]	794	DO_2D_10_10
	795	zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
	796	END_2D
[5836]	797	CASE ( np_RVO ) !* relative vorticity
[12377]	798	DO_2D_10_10
	799	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	800	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
	801	END_2D
[5836]	802	CASE ( np_MET ) !* metric term
[12377]	803	DO_2D_10_10
	804	zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	805	& - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	806	END_2D
[5836]	807	CASE ( np_CRV ) !* Coriolis + relative vorticity
[12377]	808	DO_2D_10_10
	809	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	810	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	811	& * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	812	END_2D
[5836]	813	CASE ( np_CME ) !* Coriolis + metric
[12377]	814	DO_2D_10_10
	815	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	816	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	817	END_2D
[5836]	818	CASE DEFAULT ! error
	819	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	820	END SELECT
[5836]	821	!
	822	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[12377]	823	DO_2D_10_10
	824	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	825	END_2D
[5836]	826	ENDIF
[10425]	827	END DO ! End of slab
[5836]	828	!
[10425]	829	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	830
	831	DO jk = 1, jpkm1 ! Horizontal slab
[5907]	832	!
[5836]	833	! !== horizontal fluxes ==!
[12377]	834	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	835	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[108]	836
[5836]	837	! !== compute and add the vorticity term trend =!
[1438]	838	jj = 2
	839	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[5836]	840	DO ji = 2, jpi ! split in 2 parts due to vector opt.
[10425]	841	ztne(ji,jj) = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk)
	842	ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk)
	843	ztse(ji,jj) = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk)
	844	ztsw(ji,jj) = zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk)
[108]	845	END DO
	846	DO jj = 3, jpj
[12377]	847	DO ji = 2, jpi ! vector opt. ok because we start at jj = 3
[10425]	848	ztne(ji,jj) = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk)
	849	ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk)
	850	ztse(ji,jj) = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk)
	851	ztsw(ji,jj) = zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk)
[108]	852	END DO
	853	END DO
[12377]	854	DO_2D_00_00
	855	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	856	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	857	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	858	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	859	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	860	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	861	END_2D
[108]	862	! ! ===============
	863	END DO ! End of slab
	864	! ! ===============
[455]	865	END SUBROUTINE vor_een
[216]	866
	867
[9528]	868
[12377]	869	SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]	870	!!----------------------------------------------------------------------
	871	!! * ROUTINE vor_eeT *
	872	!!
	873	!! ** Purpose : Compute the now total vorticity trend and add it to
	874	!! the general trend of the momentum equation.
	875	!!
	876	!! ** Method : Trend evaluated using now fields (centered in time)
	877	!! and the Arakawa and Lamb (1980) vector form formulation using
	878	!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
	879	!! The change consists in
[12377]	880	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[9528]	881	!!
[12377]	882	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]	883	!!
	884	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	885	!!----------------------------------------------------------------------
	886	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	887	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9528]	888	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	889	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	890	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[9528]	891	!
	892	INTEGER :: ji, jj, jk ! dummy loop indices
	893	INTEGER :: ierr ! local integer
	894	REAL(wp) :: zua, zva ! local scalars
	895	REAL(wp) :: zmsk, z1_e3t ! local scalars
[10425]	896	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy
	897	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	898	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz
[9528]	899	!!----------------------------------------------------------------------
	900	!
	901	IF( kt == nit000 ) THEN
	902	IF(lwp) WRITE(numout,*)
	903	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	904	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	905	ENDIF
	906	!
	907	! ! ===============
	908	DO jk = 1, jpkm1 ! Horizontal slab
	909	! ! ===============
	910	!
	911	!
	912	SELECT CASE( kvor ) !== vorticity considered ==!
	913	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[12377]	914	DO_2D_10_10
	915	zwz(ji,jj,jk) = ff_f(ji,jj)
	916	END_2D
[9528]	917	CASE ( np_RVO ) !* relative vorticity
[12377]	918	DO_2D_10_10
	919	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	920	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	921	& * r1_e1e2f(ji,jj)
	922	END_2D
[9528]	923	CASE ( np_MET ) !* metric term
[12377]	924	DO_2D_10_10
	925	zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	926	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	927	END_2D
[9528]	928	CASE ( np_CRV ) !* Coriolis + relative vorticity
[12377]	929	DO_2D_10_10
	930	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	931	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	932	& * r1_e1e2f(ji,jj) )
	933	END_2D
[9528]	934	CASE ( np_CME ) !* Coriolis + metric
[12377]	935	DO_2D_10_10
	936	zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	937	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	938	END_2D
[9528]	939	CASE DEFAULT ! error
	940	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
	941	END SELECT
	942	!
	943	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[12377]	944	DO_2D_10_10
	945	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	946	END_2D
[9528]	947	ENDIF
[10425]	948	END DO
	949	!
	950	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	951	!
	952	DO jk = 1, jpkm1 ! Horizontal slab
[12667]	953	!
	954	! !== horizontal fluxes ==!
[12377]	955	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	956	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[9528]	957
	958	! !== compute and add the vorticity term trend =!
	959	jj = 2
	960	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
	961	DO ji = 2, jpi ! split in 2 parts due to vector opt.
[12377]	962	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]	963	ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t
	964	ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t
	965	ztse(ji,jj) = ( zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
	966	ztsw(ji,jj) = ( zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) ) * z1_e3t
[9528]	967	END DO
	968	DO jj = 3, jpj
[12377]	969	DO ji = 2, jpi ! vector opt. ok because we start at jj = 3
	970	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]	971	ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t
	972	ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t
	973	ztse(ji,jj) = ( zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
	974	ztsw(ji,jj) = ( zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) ) * z1_e3t
[9528]	975	END DO
	976	END DO
[12377]	977	DO_2D_00_00
	978	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	979	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	980	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	981	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	982	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	983	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	984	END_2D
[9528]	985	! ! ===============
	986	END DO ! End of slab
	987	! ! ===============
	988	END SUBROUTINE vor_eeT
	989
	990
[2528]	991	SUBROUTINE dyn_vor_init
[3]	992	!!---------------------------------------------------------------------
[2528]	993	!! * ROUTINE dyn_vor_init *
[3]	994	!!
	995	!! ** Purpose : Control the consistency between cpp options for
[1438]	996	!! tracer advection schemes
[3]	997	!!----------------------------------------------------------------------
[9528]	998	INTEGER :: ji, jj, jk ! dummy loop indices
	999	INTEGER :: ioptio, ios ! local integer
[2715]	1000	!!
[9528]	1001	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, &
[13005]	1002	& ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, &
	1003	& ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters
	1004	& ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO
[3]	1005	!!----------------------------------------------------------------------
[9528]	1006	!
	1007	IF(lwp) THEN
	1008	WRITE(numout,*)
	1009	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	1010	WRITE(numout,*) '~~~~~~~~~~~~'
	1011	ENDIF
	1012	!
[4147]	1013	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
[11536]	1014	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
[4147]	1015	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
[11536]	1016	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
[4624]	1017	IF(lwm) WRITE ( numond, namdyn_vor )
[9528]	1018	!
[503]	1019	IF(lwp) THEN ! Namelist print
[7646]	1020	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
[13005]	1021	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
	1022	WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	1023	WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT
	1024	WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT
	1025	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
	1026	WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f
	1027	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
	1028	WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk
[52]	1029	ENDIF
	1030
[12667]	1031	!!an IF( ln_dynvor_msk ) CALL ctl_stop( 'dyn_vor_init: masked vorticity is not currently not available')
[9528]	1032
[5836]	1033	!!gm this should be removed when choosing a unique strategy for fmask at the coast
[3294]	1034	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	1035	! at angles with three ocean points and one land point
[5836]	1036	IF(lwp) WRITE(numout,*)
[7646]	1037	IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat
[3294]	1038	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
[12377]	1039	DO_3D_10_10( 1, jpk )
	1040	IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
[13005]	1041	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp
[12377]	1042	END_3D
[9528]	1043	!
[10425]	1044	CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
[9528]	1045	!
[3294]	1046	ENDIF
[5836]	1047	!!gm end
[3294]	1048
[5836]	1049	ioptio = 0 ! type of scheme for vorticity (set nvor_scheme)
[9528]	1050	IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF
[13005]	1051	IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF
	1052	IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF
	1053	IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF
[9528]	1054	IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF
[13005]	1055	IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF
	1056	IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF
	1057	IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF
[9528]	1058	IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF
	1059	IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF
	1060	IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF
	1061	IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF
[5836]	1062	!
[6140]	1063	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[5836]	1064	!
	1065	IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot)
[9019]	1066	ncor = np_COR ! planetary vorticity
	1067	SELECT CASE( n_dynadv )
	1068	CASE( np_LIN_dyn )
[9190]	1069	IF(lwp) WRITE(numout,*) ' ==>>> linear dynamics : total vorticity = Coriolis'
[9019]	1070	nrvm = np_COR ! planetary vorticity
	1071	ntot = np_COR ! - -
	1072	CASE( np_VEC_c2 )
[9190]	1073	IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
[13005]	1074	nrvm = np_RVO ! relative vorticity
[9019]	1075	ntot = np_CRV ! relative + planetary vorticity
	1076	CASE( np_FLX_c2 , np_FLX_ubs )
[9190]	1077	IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
[5836]	1078	nrvm = np_MET ! metric term
	1079	ntot = np_CME ! Coriolis + metric term
[9528]	1080	!
	1081	SELECT CASE( nvor_scheme ) ! pre-computed gradients for the metric term:
	1082	CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2
	1083	ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
[12377]	1084	DO_2D_00_00
	1085	di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp
	1086	dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp
	1087	END_2D
[10425]	1088	CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. ) ! Lateral boundary conditions
[9528]	1089	!
	1090	CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2e1e2f) and dj(e1v)/(2e1e2f)
	1091	ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
[12377]	1092	DO_2D_10_10
	1093	di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	1094	dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	1095	END_2D
[10425]	1096	CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. ) ! Lateral boundary conditions
[9528]	1097	END SELECT
	1098	!
[9019]	1099	END SELECT
[643]	1100
[503]	1101	IF(lwp) THEN ! Print the choice
	1102	WRITE(numout,*)
[9019]	1103	SELECT CASE( nvor_scheme )
[13005]	1104
	1105	CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)'
	1106	CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only'
	1107	CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only'
	1108	CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity'
	1109
	1110	CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)'
	1111	CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only'
	1112	CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only'
	1113	CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity'
	1114
	1115	CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)'
	1116	CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)'
	1117	CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)'
	1118	CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)'
[9019]	1119	END SELECT
[3]	1120	ENDIF
[503]	1121	!
[2528]	1122	END SUBROUTINE dyn_vor_init
[3]	1123
[503]	1124	!!==============================================================================
[3]	1125	END MODULE dynvor

Note: See TracBrowser for help on using the repository browser.

Download in other formats: