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

Context Navigation

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

Last change on this file since 14072 was 14072, checked in by laurent, 4 years ago
Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
Property svn:keywords set to `Id`
File size: 59.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
[14053]	242	REAL(wp), DIMENSION(jpi,jpj) :: 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	!
	246	IF( kt == nit000 ) THEN
	247	IF(lwp) WRITE(numout,*)
	248	IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
	249	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	250	ENDIF
	251	!
[10425]	252	!
[14053]	253	SELECT CASE( kvor ) !== relative vorticity considered ==!
	254	!
	255	CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
	256	ALLOCATE( zwz(jpi,jpj,jpk) )
	257	DO jk = 1, jpkm1 ! Horizontal slab
[13295]	258	DO_2D( 1, 0, 1, 0 )
[12377]	259	zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	260	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	261	END_2D
[14072]	262	IF( ln_dynvor_msk ) THEN ! mask relative vorticity
[13295]	263	DO_2D( 1, 0, 1, 0 )
[12377]	264	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	265	END_2D
[9528]	266	ENDIF
[10425]	267	END DO
[13226]	268	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[14053]	269	!
[10425]	270	END SELECT
	271
	272	! ! ===============
	273	DO jk = 1, jpkm1 ! Horizontal slab
[14053]	274	! ! ===============
	275	!
[10425]	276	SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
[14053]	277	!
[10425]	278	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[12377]	279	zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
[10425]	280	CASE ( np_RVO ) !* relative vorticity
[13295]	281	DO_2D( 0, 1, 0, 1 )
[14053]	282	zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
	283	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
	284	& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
[12377]	285	END_2D
[9528]	286	CASE ( np_MET ) !* metric term
[13295]	287	DO_2D( 0, 1, 0, 1 )
[14053]	288	zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
	289	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
	290	& * e3t(ji,jj,jk,Kmm)
[12377]	291	END_2D
[9528]	292	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	293	DO_2D( 0, 1, 0, 1 )
[14053]	294	zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
	295	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) &
	296	& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
[12377]	297	END_2D
[9528]	298	CASE ( np_CME ) !* Coriolis + metric
[13295]	299	DO_2D( 0, 1, 0, 1 )
[14053]	300	zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
	301	& + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
	302	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
	303	& * e3t(ji,jj,jk,Kmm)
[12377]	304	END_2D
[9528]	305	CASE DEFAULT ! error
[14053]	306	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
[9528]	307	END SELECT
	308	!
	309	! !== compute and add the vorticity term trend =!
[13295]	310	DO_2D( 0, 0, 0, 0 )
[12377]	311	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
	312	& * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
	313	& + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
	314	!
	315	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
[14072]	316	& * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
	317	& + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
[12377]	318	END_2D
[9528]	319	! ! ===============
	320	END DO ! End of slab
	321	! ! ===============
[14053]	322	!
	323	SELECT CASE( kvor ) ! deallocate zwz if necessary
	324	CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
	325	END SELECT
	326	!
[9528]	327	END SUBROUTINE vor_enT
	328
	329
[12377]	330	SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[455]	331	!!----------------------------------------------------------------------
	332	!! * ROUTINE vor_ene *
	333	!!
[14072]	334	!! ** Purpose : Compute the now total vorticity trend and add it to
[3]	335	!! the general trend of the momentum equation.
	336	!!
[14072]	337	!! ** Method : Trend evaluated using now fields (centered in time)
[5836]	338	!! and the Sadourny (1975) flux form formulation : conserves the
	339	!! horizontal kinetic energy.
[14072]	340	!! The general trend of momentum is increased due to the vorticity
[5836]	341	!! term which is given by:
[12377]	342	!! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	343	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u puu(:,:,:,Kmm)) ]
[5836]	344	!! where rvor is the relative vorticity
[3]	345	!!
[12377]	346	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[3]	347	!!
[503]	348	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	349	!!----------------------------------------------------------------------
[9019]	350	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	351	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	352	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	353	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	354	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	355	!
[5836]	356	INTEGER :: ji, jj, jk ! dummy loop indices
[14053]	357	REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars
[9019]	358	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace
[3]	359	!!----------------------------------------------------------------------
[3294]	360	!
[52]	361	IF( kt == nit000 ) THEN
	362	IF(lwp) WRITE(numout,*)
[455]	363	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	364	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	365	ENDIF
[5836]	366	!
[3]	367	! ! ===============
	368	DO jk = 1, jpkm1 ! Horizontal slab
	369	! ! ===============
[1438]	370	!
[5836]	371	SELECT CASE( kvor ) !== vorticity considered ==!
	372	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[14072]	373	zwz(:,:) = ff_f(:,:)
[5836]	374	CASE ( np_RVO ) !* relative vorticity
[13295]	375	DO_2D( 1, 0, 1, 0 )
[12377]	376	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	377	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	378	END_2D
[14053]	379	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	380	DO_2D( 1, 0, 1, 0 )
	381	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	382	END_2D
	383	ENDIF
[5836]	384	CASE ( np_MET ) !* metric term
[13295]	385	DO_2D( 1, 0, 1, 0 )
[12377]	386	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	387	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	388	END_2D
[5836]	389	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	390	DO_2D( 1, 0, 1, 0 )
[12377]	391	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	392	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	393	END_2D
[14053]	394	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity (NOT the Coriolis term)
	395	DO_2D( 1, 0, 1, 0 )
	396	zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
	397	END_2D
	398	ENDIF
[5836]	399	CASE ( np_CME ) !* Coriolis + metric
[13295]	400	DO_2D( 1, 0, 1, 0 )
[12377]	401	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	402	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	403	END_2D
[5836]	404	CASE DEFAULT ! error
	405	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	406	END SELECT
[5836]	407	!
[14053]	408	#if defined key_qco
	409	DO_2D( 1, 0, 1, 0 ) !== potential vorticity ==! (key_qco)
	410	zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
	411	END_2D
	412	#else
	413	SELECT CASE( nn_e3f_typ ) !== potential vorticity ==!
	414	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	415	DO_2D( 1, 0, 1, 0 )
[14053]	416	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	417	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	418	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	419	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
	420	IF( ze3f /= 0._wp ) THEN ; zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
	421	ELSE ; zwz(ji,jj) = 0._wp
	422	ENDIF
[12377]	423	END_2D
[14053]	424	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
	425	DO_2D( 1, 0, 1, 0 )
	426	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	427	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	428	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	429	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
	430	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	431	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	432	IF( ze3f /= 0._wp ) THEN ; zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
	433	ELSE ; zwz(ji,jj) = 0._wp
	434	ENDIF
	435	END_2D
	436	END SELECT
	437	#endif
	438	! !== horizontal fluxes ==!
	439	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	440	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
	441	!
[5836]	442	! !== compute and add the vorticity term trend =!
[13295]	443	DO_2D( 0, 0, 0, 0 )
[12377]	444	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	445	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	446	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	447	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
	448	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]	449	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]	450	END_2D
[3]	451	! ! ===============
	452	END DO ! End of slab
	453	! ! ===============
[455]	454	END SUBROUTINE vor_ene
[216]	455
	456
[12377]	457	SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[3]	458	!!----------------------------------------------------------------------
[455]	459	!! * ROUTINE vor_ens *
[3]	460	!!
	461	!! ** Purpose : Compute the now total vorticity trend and add it to
	462	!! the general trend of the momentum equation.
	463	!!
	464	!! ** Method : Trend evaluated using now fields (centered in time)
	465	!! and the Sadourny (1975) flux FORM formulation : conserves the
	466	!! potential enstrophy of a horizontally non-divergent flow. the
	467	!! trend of the vorticity term is given by:
[12377]	468	!! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	469	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
	470	!! Add this trend to the general momentum trend:
	471	!! (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
[3]	472	!!
[12377]	473	!! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
[3]	474	!!
[503]	475	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	476	!!----------------------------------------------------------------------
[9019]	477	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	478	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	479	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	480	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	481	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	482	!
[5836]	483	INTEGER :: ji, jj, jk ! dummy loop indices
[14053]	484	REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars
[9019]	485	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! 2D workspace
[3]	486	!!----------------------------------------------------------------------
[3294]	487	!
[52]	488	IF( kt == nit000 ) THEN
	489	IF(lwp) WRITE(numout,*)
[455]	490	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	491	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	492	ENDIF
[3]	493	! ! ===============
	494	DO jk = 1, jpkm1 ! Horizontal slab
	495	! ! ===============
[1438]	496	!
[5836]	497	SELECT CASE( kvor ) !== vorticity considered ==!
	498	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[14072]	499	zwz(:,:) = ff_f(:,:)
[5836]	500	CASE ( np_RVO ) !* relative vorticity
[13295]	501	DO_2D( 1, 0, 1, 0 )
[12377]	502	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	503	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	504	END_2D
[14053]	505	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	506	DO_2D( 1, 0, 1, 0 )
	507	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	508	END_2D
	509	ENDIF
[5836]	510	CASE ( np_MET ) !* metric term
[13295]	511	DO_2D( 1, 0, 1, 0 )
[12377]	512	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	513	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	514	END_2D
[5836]	515	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	516	DO_2D( 1, 0, 1, 0 )
[12377]	517	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	518	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	519	END_2D
[14053]	520	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity (NOT the Coriolis term)
	521	DO_2D( 1, 0, 1, 0 )
	522	zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
	523	END_2D
	524	ENDIF
[5836]	525	CASE ( np_CME ) !* Coriolis + metric
[13295]	526	DO_2D( 1, 0, 1, 0 )
[12377]	527	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	528	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	529	END_2D
[5836]	530	CASE DEFAULT ! error
	531	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	532	END SELECT
[1438]	533	!
[14053]	534	!
	535	#if defined key_qco
	536	DO_2D( 1, 0, 1, 0 ) !== potential vorticity ==! (key_qco)
	537	zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
	538	END_2D
	539	#else
	540	SELECT CASE( nn_e3f_typ ) !== potential vorticity ==!
	541	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	542	DO_2D( 1, 0, 1, 0 )
[14053]	543	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	544	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	545	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	546	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
	547	IF( ze3f /= 0._wp ) THEN ; zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
	548	ELSE ; zwz(ji,jj) = 0._wp
	549	ENDIF
[12377]	550	END_2D
[14053]	551	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
	552	DO_2D( 1, 0, 1, 0 )
	553	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	554	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	555	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	556	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
	557	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	558	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	559	IF( ze3f /= 0._wp ) THEN ; zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
	560	ELSE ; zwz(ji,jj) = 0._wp
	561	ENDIF
	562	END_2D
	563	END SELECT
	564	#endif
	565	! !== horizontal fluxes ==!
	566	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	567	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[5836]	568	!
	569	! !== compute and add the vorticity term trend =!
[13295]	570	DO_2D( 0, 0, 0, 0 )
[12377]	571	zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	572	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
	573	zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	574	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
	575	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	576	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	577	END_2D
[3]	578	! ! ===============
	579	END DO ! End of slab
	580	! ! ===============
[455]	581	END SUBROUTINE vor_ens
[216]	582
	583
[12377]	584	SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[108]	585	!!----------------------------------------------------------------------
[455]	586	!! * ROUTINE vor_een *
[108]	587	!!
[14072]	588	!! ** Purpose : Compute the now total vorticity trend and add it to
[108]	589	!! the general trend of the momentum equation.
	590	!!
[14072]	591	!! ** Method : Trend evaluated using now fields (centered in time)
	592	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	593	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	594	!! when horizontal divergence is zero (see the NEMO documentation)
[12377]	595	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[108]	596	!!
[12377]	597	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[108]	598	!!
[503]	599	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	600	!!----------------------------------------------------------------------
[9019]	601	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	602	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	603	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	604	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	605	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[5836]	606	!
	607	INTEGER :: ji, jj, jk ! dummy loop indices
	608	INTEGER :: ierr ! local integer
	609	REAL(wp) :: zua, zva ! local scalars
[9528]	610	REAL(wp) :: zmsk, ze3f ! local scalars
[13546]	611	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f
	612	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	613	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
[108]	614	!!----------------------------------------------------------------------
[3294]	615	!
[108]	616	IF( kt == nit000 ) THEN
	617	IF(lwp) WRITE(numout,*)
[455]	618	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	619	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]	620	ENDIF
[5836]	621	!
	622	! ! ===============
	623	DO jk = 1, jpkm1 ! Horizontal slab
	624	! ! ===============
	625	!
[14053]	626	#if defined key_qco
	627	DO_2D( 1, 0, 1, 0 ) ! == reciprocal of e3 at F-point (key_qco)
	628	z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
	629	END_2D
	630	#else
	631	SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point
[5836]	632	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	633	DO_2D( 1, 0, 1, 0 )
[13237]	634	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	635	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	636	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	637	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
[12377]	638	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
	639	ELSE ; z1_e3f(ji,jj) = 0._wp
	640	ENDIF
	641	END_2D
[5836]	642	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[13295]	643	DO_2D( 1, 0, 1, 0 )
[13237]	644	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	645	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	646	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	647	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
[12377]	648	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	649	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	650	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
	651	ELSE ; z1_e3f(ji,jj) = 0._wp
	652	ENDIF
	653	END_2D
[5836]	654	END SELECT
[14053]	655	#endif
[5836]	656	!
	657	SELECT CASE( kvor ) !== vorticity considered ==!
[14053]	658	!
[5836]	659	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[13295]	660	DO_2D( 1, 0, 1, 0 )
[12377]	661	zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
	662	END_2D
[5836]	663	CASE ( np_RVO ) !* relative vorticity
[13295]	664	DO_2D( 1, 0, 1, 0 )
[12377]	665	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	666	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
	667	END_2D
[14053]	668	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	669	DO_2D( 1, 0, 1, 0 )
	670	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	671	END_2D
	672	ENDIF
[5836]	673	CASE ( np_MET ) !* metric term
[13295]	674	DO_2D( 1, 0, 1, 0 )
[12377]	675	zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	676	& - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	677	END_2D
[5836]	678	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	679	DO_2D( 1, 0, 1, 0 )
[12377]	680	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	681	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	682	& * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	683	END_2D
[14053]	684	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	685	DO_2D( 1, 0, 1, 0 )
[14072]	686	zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
[14053]	687	END_2D
	688	ENDIF
[5836]	689	CASE ( np_CME ) !* Coriolis + metric
[13295]	690	DO_2D( 1, 0, 1, 0 )
[12377]	691	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	692	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	693	END_2D
[5836]	694	CASE DEFAULT ! error
	695	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	696	END SELECT
[14053]	697	! ! ===============
[10425]	698	END DO ! End of slab
[14053]	699	! ! ===============
	700	!
[13226]	701	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[14053]	702	!
	703	! ! ===============
[10425]	704	DO jk = 1, jpkm1 ! Horizontal slab
[14053]	705	! ! ===============
[5907]	706	!
[5836]	707	! !== horizontal fluxes ==!
[12377]	708	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	709	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[14053]	710	!
[5836]	711	! !== compute and add the vorticity term trend =!
[14053]	712	DO_2D( 0, 1, 0, 1 )
	713	ztne(ji,jj) = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk)
	714	ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk)
	715	ztse(ji,jj) = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk)
	716	ztsw(ji,jj) = zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk)
	717	END_2D
	718	!
[13295]	719	DO_2D( 0, 0, 0, 0 )
[12377]	720	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	721	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	722	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	723	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	724	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	725	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	726	END_2D
[108]	727	! ! ===============
	728	END DO ! End of slab
	729	! ! ===============
[455]	730	END SUBROUTINE vor_een
[216]	731
	732
[12377]	733	SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]	734	!!----------------------------------------------------------------------
	735	!! * ROUTINE vor_eeT *
	736	!!
[14072]	737	!! ** Purpose : Compute the now total vorticity trend and add it to
[9528]	738	!! the general trend of the momentum equation.
	739	!!
[14072]	740	!! ** Method : Trend evaluated using now fields (centered in time)
	741	!! and the Arakawa and Lamb (1980) vector form formulation using
[9528]	742	!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
[14072]	743	!! The change consists in
[12377]	744	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[9528]	745	!!
[12377]	746	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]	747	!!
	748	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	749	!!----------------------------------------------------------------------
[14053]	750	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	751	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[14053]	752	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
	753	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	754	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[9528]	755	!
	756	INTEGER :: ji, jj, jk ! dummy loop indices
	757	INTEGER :: ierr ! local integer
	758	REAL(wp) :: zua, zva ! local scalars
	759	REAL(wp) :: zmsk, z1_e3t ! local scalars
[14072]	760	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy
[13546]	761	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	762	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
[9528]	763	!!----------------------------------------------------------------------
	764	!
	765	IF( kt == nit000 ) THEN
	766	IF(lwp) WRITE(numout,*)
[14053]	767	IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
[9528]	768	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	769	ENDIF
	770	!
	771	! ! ===============
	772	DO jk = 1, jpkm1 ! Horizontal slab
	773	! ! ===============
	774	!
	775	!
	776	SELECT CASE( kvor ) !== vorticity considered ==!
	777	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[13295]	778	DO_2D( 1, 0, 1, 0 )
[12377]	779	zwz(ji,jj,jk) = ff_f(ji,jj)
	780	END_2D
[9528]	781	CASE ( np_RVO ) !* relative vorticity
[13295]	782	DO_2D( 1, 0, 1, 0 )
[12377]	783	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	784	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	785	& * r1_e1e2f(ji,jj)
	786	END_2D
[14053]	787	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	788	DO_2D( 1, 0, 1, 0 )
	789	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	790	END_2D
	791	ENDIF
[9528]	792	CASE ( np_MET ) !* metric term
[13295]	793	DO_2D( 1, 0, 1, 0 )
[12377]	794	zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	795	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	796	END_2D
[9528]	797	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	798	DO_2D( 1, 0, 1, 0 )
[12377]	799	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	800	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	801	& * r1_e1e2f(ji,jj) )
	802	END_2D
[14053]	803	IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
	804	DO_2D( 1, 0, 1, 0 )
[14072]	805	zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
[14053]	806	END_2D
	807	ENDIF
[9528]	808	CASE ( np_CME ) !* Coriolis + metric
[13295]	809	DO_2D( 1, 0, 1, 0 )
[12377]	810	zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	811	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	812	END_2D
[9528]	813	CASE DEFAULT ! error
	814	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
	815	END SELECT
	816	!
[14053]	817	! ! ===============
	818	END DO ! End of slab
	819	! ! ===============
[10425]	820	!
[13226]	821	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[10425]	822	!
[14053]	823	! ! ===============
[10425]	824	DO jk = 1, jpkm1 ! Horizontal slab
[14053]	825	! ! ===============
	826	!
	827	! !== horizontal fluxes ==!
[12377]	828	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	829	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[14053]	830	!
[9528]	831	! !== compute and add the vorticity term trend =!
[14053]	832	DO_2D( 0, 1, 0, 1 )
	833	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
	834	ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t
	835	ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t
	836	ztse(ji,jj) = ( zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
	837	ztsw(ji,jj) = ( zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) ) * z1_e3t
	838	END_2D
	839	!
[13295]	840	DO_2D( 0, 0, 0, 0 )
[12377]	841	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	842	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	843	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	844	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	845	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	846	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	847	END_2D
[9528]	848	! ! ===============
	849	END DO ! End of slab
	850	! ! ===============
	851	END SUBROUTINE vor_eeT
	852
	853
[2528]	854	SUBROUTINE dyn_vor_init
[3]	855	!!---------------------------------------------------------------------
[2528]	856	!! * ROUTINE dyn_vor_init *
[3]	857	!!
	858	!! ** Purpose : Control the consistency between cpp options for
[1438]	859	!! tracer advection schemes
[3]	860	!!----------------------------------------------------------------------
[9528]	861	INTEGER :: ji, jj, jk ! dummy loop indices
	862	INTEGER :: ioptio, ios ! local integer
[14053]	863	REAL(wp) :: zmsk ! local scalars
[2715]	864	!!
[9528]	865	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, &
[14053]	866	& ln_dynvor_een, nn_e3f_typ , ln_dynvor_mix, ln_dynvor_msk
[3]	867	!!----------------------------------------------------------------------
[9528]	868	!
	869	IF(lwp) THEN
	870	WRITE(numout,*)
	871	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	872	WRITE(numout,*) '~~~~~~~~~~~~'
	873	ENDIF
	874	!
[4147]	875	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
[11536]	876	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
[4147]	877	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
[11536]	878	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
[4624]	879	IF(lwm) WRITE ( numond, namdyn_vor )
[9528]	880	!
[503]	881	IF(lwp) THEN ! Namelist print
[7646]	882	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
	883	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
[9528]	884	WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	885	WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT
	886	WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT
[7646]	887	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
[14053]	888	WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_e3f_typ = ', nn_e3f_typ
[9528]	889	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
[7646]	890	WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk
[52]	891	ENDIF
	892
[5836]	893	!!gm this should be removed when choosing a unique strategy for fmask at the coast
[3294]	894	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	895	! at angles with three ocean points and one land point
[5836]	896	IF(lwp) WRITE(numout,*)
[7646]	897	IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat
[3294]	898	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
[13295]	899	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	900	IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
[12793]	901	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp
[12377]	902	END_3D
[9528]	903	!
[10425]	904	CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
[9528]	905	!
[3294]	906	ENDIF
[5836]	907	!!gm end
[3294]	908
[5836]	909	ioptio = 0 ! type of scheme for vorticity (set nvor_scheme)
[9528]	910	IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF
	911	IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF
	912	IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF
	913	IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF
	914	IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF
	915	IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF
[5836]	916	!
[6140]	917	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[14072]	918	!
[5836]	919	IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot)
[9019]	920	ncor = np_COR ! planetary vorticity
	921	SELECT CASE( n_dynadv )
	922	CASE( np_LIN_dyn )
[9190]	923	IF(lwp) WRITE(numout,*) ' ==>>> linear dynamics : total vorticity = Coriolis'
[9019]	924	nrvm = np_COR ! planetary vorticity
	925	ntot = np_COR ! - -
	926	CASE( np_VEC_c2 )
[14072]	927	IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
[5836]	928	nrvm = np_RVO ! relative vorticity
[14072]	929	ntot = np_CRV ! relative + planetary vorticity
[9019]	930	CASE( np_FLX_c2 , np_FLX_ubs )
[9190]	931	IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
[5836]	932	nrvm = np_MET ! metric term
	933	ntot = np_CME ! Coriolis + metric term
[9528]	934	!
	935	SELECT CASE( nvor_scheme ) ! pre-computed gradients for the metric term:
	936	CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2
	937	ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
[13295]	938	DO_2D( 0, 0, 0, 0 )
[12377]	939	di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp
	940	dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp
	941	END_2D
[13226]	942	CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp ) ! Lateral boundary conditions
[9528]	943	!
	944	CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2e1e2f) and dj(e1v)/(2e1e2f)
	945	ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
[13295]	946	DO_2D( 1, 0, 1, 0 )
[12377]	947	di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	948	dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	949	END_2D
[13226]	950	CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp ) ! Lateral boundary conditions
[9528]	951	END SELECT
	952	!
[9019]	953	END SELECT
[14053]	954	#if defined key_qco
	955	SELECT CASE( nvor_scheme ) ! qco case: pre-computed a specific e3f_0 for some vorticity schemes
	956	CASE( np_ENS , np_ENE , np_EEN , np_MIX )
	957	!
	958	ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
	959	!
	960	SELECT CASE( nn_e3f_typ )
	961	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
	962	DO_3D( 0, 0, 0, 0, 1, jpk )
	963	e3f_0vor(ji,jj,jk) = ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) &
	964	& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &
	965	& + e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
	966	& + e3t_0(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) 0.25_wp
	967	END_3D
	968	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
	969	DO_3D( 0, 0, 0, 0, 1, jpk )
	970	zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk) &
	971	& + tmask(ji,jj ,jk) +tmask(ji+1,jj ,jk) )
	972	!
[14072]	973	IF( zmsk /= 0._wp ) THEN
[14053]	974	e3f_0vor(ji,jj,jk) = ( e3t_0(ji ,jj+1,jk)*tmask(ji ,jj+1,jk) &
	975	& + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &
	976	& + e3t_0(ji ,jj ,jk)*tmask(ji ,jj ,jk) &
	977	& + e3t_0(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) / zmsk
	978	ENDIF
	979	END_3D
	980	END SELECT
	981	!
	982	CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
	983	! ! insure e3f_0vor /= 0
	984	WHERE( e3f_0vor(:,:,:) == 0._wp ) e3f_0vor(:,:,:) = e3f_0(:,:,:)
	985	!
	986	END SELECT
	987	!
	988	#endif
[503]	989	IF(lwp) THEN ! Print the choice
	990	WRITE(numout,*)
[9019]	991	SELECT CASE( nvor_scheme )
[9528]	992	CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)'
	993	CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)'
	994	CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)'
[14053]	995	IF( ln_dynadv_vec ) CALL ctl_warn('dyn_vor_init: ENT scheme may not work in vector form')
[9528]	996	CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)'
	997	CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)'
	998	CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)'
[14072]	999	END SELECT
[3]	1000	ENDIF
[503]	1001	!
[2528]	1002	END SUBROUTINE dyn_vor_init
[3]	1003
[503]	1004	!!==============================================================================
[3]	1005	END MODULE dynvor

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

Download in other formats: