Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

dynvor.F90

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: