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/2021/dev_r14318_RK3_stage1/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90 @ 14618

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

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

Download in other formats: