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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynvor.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DYN/dynvor.F90 @ 10170

Last change on this file since 10170 was 10170, checked in by smasson, 6 years ago
dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of lbc_lnk, see #2133
Property svn:keywords set to `Id`
File size: 52.4 KB

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

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

Download in other formats: