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_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90 @ 13696

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

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

Download in other formats: