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 branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago
Ediag branche: #927 split TRA/DYN trd computation
Property svn:keywords set to `Id`
File size: 39.5 KB

Rev	Line
[3]	1	MODULE dynvor
	2	!!======================================================================
	3	!! * MODULE dynvor *
	4	!! Ocean dynamics: Update the momentum trend with the relative and
	5	!! planetary vorticity trends
	6	!!======================================================================
[2715]	7	!! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code
	8	!! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code
	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
	16	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[3318]	17	!! 3.5 ! 2012-02 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
[503]	18	!!----------------------------------------------------------------------
[3]	19
	20	!!----------------------------------------------------------------------
[2528]	21	!! dyn_vor : Update the momentum trend with the vorticity trend
	22	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
	23	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
	24	!! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
	25	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
	26	!! dyn_vor_init : set and control of the different vorticity option
[3]	27	!!----------------------------------------------------------------------
[503]	28	USE oce ! ocean dynamics and tracers
	29	USE dom_oce ! ocean space and time domain
[3294]	30	USE dommsk ! ocean mask
[643]	31	USE dynadv ! momentum advection (use ln_dynadv_vec value)
[3318]	32	USE trd_oce ! trends: ocean variables
	33	USE trddyn ! trend manager: dynamics
[503]	34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	35	USE prtctl ! Print control
	36	USE in_out_manager ! I/O manager
[3294]	37	USE lib_mpp ! MPP library
	38	USE wrk_nemo ! Memory Allocation
	39	USE timing ! Timing
[3]	40
[3294]	41
[3]	42	IMPLICIT NONE
	43	PRIVATE
	44
[2528]	45	PUBLIC dyn_vor ! routine called by step.F90
	46	PUBLIC dyn_vor_init ! routine called by opa.F90
[3]	47
[1601]	48	! !!* Namelist namdyn_vor: vorticity term
[32]	49	LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme
	50	LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme
	51	LOGICAL, PUBLIC :: ln_dynvor_mix = .FALSE. !: mixed scheme
[503]	52	LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme
[3]	53
[503]	54	INTEGER :: nvor = 0 ! type of vorticity trend used
[643]	55	INTEGER :: ncor = 1 ! coriolis
	56	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
	57	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
[455]	58
[3]	59	!! * Substitutions
	60	# include "domzgr_substitute.h90"
	61	# include "vectopt_loop_substitute.h90"
	62	!!----------------------------------------------------------------------
[2528]	63	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	64	!! $Id$
[2715]	65	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	66	!!----------------------------------------------------------------------
	67	CONTAINS
	68
[455]	69	SUBROUTINE dyn_vor( kt )
[3]	70	!!----------------------------------------------------------------------
	71	!!
[455]	72	!! ** Purpose : compute the lateral ocean tracer physics.
	73	!!
	74	!! ** Action : - Update (ua,va) with the now vorticity term trend
[503]	75	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[3318]	76	!! and planetary vorticity trends) (l_trddyn=T)
[503]	77	!!----------------------------------------------------------------------
[3294]	78	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	79	!
[3294]	80	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
[455]	81	!!----------------------------------------------------------------------
[2715]	82	!
[3294]	83	IF( nn_timing == 1 ) CALL timing_start('dyn_vor')
	84	!
	85	IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
	86	!
[643]	87	! ! vorticity term
[455]	88	SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend
[643]	89	!
[455]	90	CASE ( -1 ) ! esopa: test all possibility with control print
[643]	91	CALL vor_ene( kt, ntot, ua, va )
[503]	92	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
	93	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	94	CALL vor_ens( kt, ntot, ua, va )
[503]	95	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
	96	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[455]	97	CALL vor_mix( kt )
[503]	98	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
	99	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	100	CALL vor_een( kt, ntot, ua, va )
[503]	101	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
	102	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	103	!
[455]	104	CASE ( 0 ) ! energy conserving scheme
	105	IF( l_trddyn ) THEN
	106	ztrdu(:,:,:) = ua(:,:,:)
	107	ztrdv(:,:,:) = va(:,:,:)
[643]	108	CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	109	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	110	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	111	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	112	ztrdu(:,:,:) = ua(:,:,:)
	113	ztrdv(:,:,:) = va(:,:,:)
[643]	114	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	115	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	116	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	117	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	118	ELSE
[643]	119	CALL vor_ene( kt, ntot, ua, va ) ! total vorticity
[455]	120	ENDIF
[643]	121	!
[455]	122	CASE ( 1 ) ! enstrophy conserving scheme
	123	IF( l_trddyn ) THEN
	124	ztrdu(:,:,:) = ua(:,:,:)
	125	ztrdv(:,:,:) = va(:,:,:)
[643]	126	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	127	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	128	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	129	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	130	ztrdu(:,:,:) = ua(:,:,:)
	131	ztrdv(:,:,:) = va(:,:,:)
[643]	132	CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	133	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	134	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	135	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	136	ELSE
[643]	137	CALL vor_ens( kt, ntot, ua, va ) ! total vorticity
[455]	138	ENDIF
[643]	139	!
[455]	140	CASE ( 2 ) ! mixed ene-ens scheme
	141	IF( l_trddyn ) THEN
	142	ztrdu(:,:,:) = ua(:,:,:)
	143	ztrdv(:,:,:) = va(:,:,:)
[643]	144	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)
[455]	145	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	146	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	147	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	148	ztrdu(:,:,:) = ua(:,:,:)
	149	ztrdv(:,:,:) = va(:,:,:)
[643]	150	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)
[455]	151	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	152	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	153	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	154	ELSE
	155	CALL vor_mix( kt ) ! total vorticity (mix=ens-ene)
	156	ENDIF
[643]	157	!
[455]	158	CASE ( 3 ) ! energy and enstrophy conserving scheme
	159	IF( l_trddyn ) THEN
	160	ztrdu(:,:,:) = ua(:,:,:)
	161	ztrdv(:,:,:) = va(:,:,:)
[643]	162	CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	163	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	164	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	165	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	166	ztrdu(:,:,:) = ua(:,:,:)
	167	ztrdv(:,:,:) = va(:,:,:)
[643]	168	CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	169	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	170	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[3318]	171	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	172	ELSE
[643]	173	CALL vor_een( kt, ntot, ua, va ) ! total vorticity
[455]	174	ENDIF
[643]	175	!
[455]	176	END SELECT
[2715]	177	!
[455]	178	! ! print sum trends (used for debugging)
[2715]	179	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, &
[455]	180	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]	181	!
[3294]	182	IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
	183	!
	184	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor')
	185	!
[455]	186	END SUBROUTINE dyn_vor
	187
	188
[643]	189	SUBROUTINE vor_ene( kt, kvor, pua, pva )
[455]	190	!!----------------------------------------------------------------------
	191	!! * ROUTINE vor_ene *
	192	!!
[3]	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 the Sadourny (1975) flux form formulation : conserves the
	198	!! horizontal kinetic energy.
	199	!! The trend of the vorticity term is given by:
[455]	200	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]	201	!! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ]
	202	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ]
	203	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	204	!! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ]
	205	!! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ]
	206	!! Add this trend to the general momentum trend (ua,va):
	207	!! (ua,va) = (ua,va) + ( voru , vorv )
	208	!!
	209	!! ** Action : - Update (ua,va) with the now vorticity term trend
	210	!!
[503]	211	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	212	!!----------------------------------------------------------------------
[2715]	213	!
[643]	214	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	215	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	216	! ! =nrvm (relative vorticity or metric)
[643]	217	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	218	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	219	!
	220	INTEGER :: ji, jj, jk ! dummy loop indices
	221	REAL(wp) :: zx1, zy1, zfact2, zx2, zy2 ! local scalars
[3294]	222	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
[3]	223	!!----------------------------------------------------------------------
[3294]	224	!
	225	IF( nn_timing == 1 ) CALL timing_start('vor_ene')
	226	!
	227	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	228	!
[52]	229	IF( kt == nit000 ) THEN
	230	IF(lwp) WRITE(numout,*)
[455]	231	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	232	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	233	ENDIF
[3]	234
[1438]	235	zfact2 = 0.5 * 0.5 ! Local constant initialization
[216]	236
[455]	237	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	238	! ! ===============
	239	DO jk = 1, jpkm1 ! Horizontal slab
	240	! ! ===============
[1438]	241	!
[3]	242	! Potential vorticity and horizontal fluxes
	243	! -----------------------------------------
[643]	244	SELECT CASE( kvor ) ! vorticity considered
	245	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	246	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	247	CASE ( 3 ) ! metric term
	248	DO jj = 1, jpjm1
	249	DO ji = 1, fs_jpim1 ! vector opt.
	250	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	251	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	252	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	253	END DO
	254	END DO
	255	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	256	CASE ( 5 ) ! total (coriolis + metric)
	257	DO jj = 1, jpjm1
	258	DO ji = 1, fs_jpim1 ! vector opt.
	259	zwz(ji,jj) = ( ff (ji,jj) &
	260	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	261	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	262	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
	263	& )
	264	END DO
	265	END DO
[455]	266	END SELECT
	267
	268	IF( ln_sco ) THEN
	269	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
[3]	270	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	271	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	272	ELSE
	273	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	274	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	275	ENDIF
	276
	277	! Compute and add the vorticity term trend
	278	! ----------------------------------------
	279	DO jj = 2, jpjm1
	280	DO ji = fs_2, fs_jpim1 ! vector opt.
	281	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	282	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	283	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	284	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
[455]	285	pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	286	pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[3]	287	END DO
	288	END DO
	289	! ! ===============
	290	END DO ! End of slab
	291	! ! ===============
[3294]	292	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	293	!
[3294]	294	IF( nn_timing == 1 ) CALL timing_stop('vor_ene')
	295	!
[455]	296	END SUBROUTINE vor_ene
[216]	297
	298
[455]	299	SUBROUTINE vor_mix( kt )
[3]	300	!!----------------------------------------------------------------------
[455]	301	!! * ROUTINE vor_mix *
[3]	302	!!
	303	!! ** Purpose : Compute the now total vorticity trend and add it to
	304	!! the general trend of the momentum equation.
	305	!!
	306	!! ** Method : Trend evaluated using now fields (centered in time)
	307	!! Mixte formulation : conserves the potential enstrophy of a hori-
	308	!! zontally non-divergent flow for (rotzu x uh), the relative vor-
	309	!! ticity term and the horizontal kinetic energy for (f x uh), the
	310	!! coriolis term. the now trend of the vorticity term is given by:
[455]	311	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]	312	!! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
	313	!! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]
	314	!! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
	315	!! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]
	316	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	317	!! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]
	318	!! +1/e1u mj-1[ f mi(e1v vn) ]
	319	!! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]
	320	!! +1/e2v mi-1[ f mj(e2u un) ]
	321	!! Add this now trend to the general momentum trend (ua,va):
	322	!! (ua,va) = (ua,va) + ( voru , vorv )
	323	!!
	324	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
	325	!!
[503]	326	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	327	!!----------------------------------------------------------------------
[2715]	328	!
[503]	329	INTEGER, INTENT(in) :: kt ! ocean timestep index
[2715]	330	!
[1438]	331	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	332	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars
	333	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -
[3294]	334	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[3]	335	!!----------------------------------------------------------------------
[3294]	336	!
	337	IF( nn_timing == 1 ) CALL timing_start('vor_mix')
	338	!
	339	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )
	340	!
[52]	341	IF( kt == nit000 ) THEN
	342	IF(lwp) WRITE(numout,*)
[455]	343	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
	344	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	345	ENDIF
[3]	346
[1438]	347	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	348	zfact2 = 0.5 * 0.5
	349
[455]	350	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
[3]	351	! ! ===============
	352	DO jk = 1, jpkm1 ! Horizontal slab
	353	! ! ===============
[1438]	354	!
[3]	355	! Relative and planetary potential vorticity and horizontal fluxes
	356	! ----------------------------------------------------------------
[455]	357	IF( ln_sco ) THEN
[643]	358	IF( ln_dynadv_vec ) THEN
	359	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
	360	ELSE
	361	DO jj = 1, jpjm1
	362	DO ji = 1, fs_jpim1 ! vector opt.
	363	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	364	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	365	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
	366	END DO
	367	END DO
	368	ENDIF
[3]	369	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
	370	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	371	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	372	ELSE
[643]	373	IF( ln_dynadv_vec ) THEN
	374	zww(:,:) = rotn(:,:,jk)
	375	ELSE
	376	DO jj = 1, jpjm1
	377	DO ji = 1, fs_jpim1 ! vector opt.
	378	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	379	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	380	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
	381	END DO
	382	END DO
	383	ENDIF
	384	zwz(:,:) = ff (:,:)
[3]	385	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	386	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	387	ENDIF
	388
	389	! Compute and add the vorticity term trend
	390	! ----------------------------------------
	391	DO jj = 2, jpjm1
	392	DO ji = fs_2, fs_jpim1 ! vector opt.
	393	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	394	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	395	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	396	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	397	! enstrophy conserving formulation for relative vorticity term
	398	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
	399	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
	400	! energy conserving formulation for planetary vorticity term
	401	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	402	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[503]	403	! mixed vorticity trend added to the momentum trends
[3]	404	ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
	405	va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
	406	END DO
	407	END DO
	408	! ! ===============
	409	END DO ! End of slab
	410	! ! ===============
[3294]	411	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )
[2715]	412	!
[3294]	413	IF( nn_timing == 1 ) CALL timing_stop('vor_mix')
	414	!
[455]	415	END SUBROUTINE vor_mix
[216]	416
	417
[643]	418	SUBROUTINE vor_ens( kt, kvor, pua, pva )
[3]	419	!!----------------------------------------------------------------------
[455]	420	!! * ROUTINE vor_ens *
[3]	421	!!
	422	!! ** Purpose : Compute the now total vorticity trend and add it to
	423	!! the general trend of the momentum equation.
	424	!!
	425	!! ** Method : Trend evaluated using now fields (centered in time)
	426	!! and the Sadourny (1975) flux FORM formulation : conserves the
	427	!! potential enstrophy of a horizontally non-divergent flow. the
	428	!! trend of the vorticity term is given by:
[455]	429	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
[3]	430	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
	431	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
	432	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	433	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
	434	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
	435	!! Add this trend to the general momentum trend (ua,va):
	436	!! (ua,va) = (ua,va) + ( voru , vorv )
	437	!!
	438	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
	439	!!
[503]	440	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	441	!!----------------------------------------------------------------------
[2715]	442	!
[643]	443	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	444	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
	445	! ! =nrvm (relative vorticity or metric)
	446	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	447	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	448	!
[503]	449	INTEGER :: ji, jj, jk ! dummy loop indices
	450	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
[3294]	451	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[3]	452	!!----------------------------------------------------------------------
[3294]	453	!
	454	IF( nn_timing == 1 ) CALL timing_start('vor_ens')
	455	!
	456	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	457	!
[52]	458	IF( kt == nit000 ) THEN
	459	IF(lwp) WRITE(numout,*)
[455]	460	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	461	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	462	ENDIF
[3]	463
[1438]	464	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	465
[455]	466	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	467	! ! ===============
	468	DO jk = 1, jpkm1 ! Horizontal slab
	469	! ! ===============
[1438]	470	!
[3]	471	! Potential vorticity and horizontal fluxes
	472	! -----------------------------------------
[643]	473	SELECT CASE( kvor ) ! vorticity considered
	474	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	475	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	476	CASE ( 3 ) ! metric term
	477	DO jj = 1, jpjm1
	478	DO ji = 1, fs_jpim1 ! vector opt.
	479	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	480	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	481	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	482	END DO
	483	END DO
	484	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	485	CASE ( 5 ) ! total (coriolis + metric)
	486	DO jj = 1, jpjm1
	487	DO ji = 1, fs_jpim1 ! vector opt.
	488	zwz(ji,jj) = ( ff (ji,jj) &
	489	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	490	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	491	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	492	& )
	493	END DO
	494	END DO
[455]	495	END SELECT
[1438]	496	!
[455]	497	IF( ln_sco ) THEN
[3]	498	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	499	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	500	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
	501	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
	502	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
[3]	503	END DO
	504	END DO
	505	ELSE
	506	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	507	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	508	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
	509	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
[3]	510	END DO
	511	END DO
	512	ENDIF
[1438]	513	!
[3]	514	! Compute and add the vorticity term trend
	515	! ----------------------------------------
	516	DO jj = 2, jpjm1
	517	DO ji = fs_2, fs_jpim1 ! vector opt.
[455]	518	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
[503]	519	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
[455]	520	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
[503]	521	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
[455]	522	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	523	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[3]	524	END DO
	525	END DO
	526	! ! ===============
	527	END DO ! End of slab
	528	! ! ===============
[3294]	529	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	530	!
[3294]	531	IF( nn_timing == 1 ) CALL timing_stop('vor_ens')
	532	!
[455]	533	END SUBROUTINE vor_ens
[216]	534
	535
[643]	536	SUBROUTINE vor_een( kt, kvor, pua, pva )
[108]	537	!!----------------------------------------------------------------------
[455]	538	!! * ROUTINE vor_een *
[108]	539	!!
	540	!! ** Purpose : Compute the now total vorticity trend and add it to
	541	!! the general trend of the momentum equation.
	542	!!
	543	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	544	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	545	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	546	!! when horizontal divergence is zero (see the NEMO documentation)
	547	!! Add this trend to the general momentum trend (ua,va).
[108]	548	!!
	549	!! ** Action : - Update (ua,va) with the now vorticity term trend
	550	!!
[503]	551	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	552	!!----------------------------------------------------------------------
[2715]	553	!
[643]	554	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	555	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	556	! ! =nrvm (relative vorticity or metric)
[643]	557	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	558	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[218]	559	!!
[3294]	560	INTEGER :: ji, jj, jk ! dummy loop indices
	561	INTEGER :: ierr ! local integer
	562	REAL(wp) :: zfac12, zua, zva ! local scalars
	563	! ! 3D workspace
	564	REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz
	565	REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse
	566	#if defined key_vvl
	567	REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T)
	568	#endif
[2715]	569	#if ! defined key_vvl
[3294]	570	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all
[1438]	571	#endif
[108]	572	!!----------------------------------------------------------------------
[3294]	573	!
	574	IF( nn_timing == 1 ) CALL timing_start('vor_een')
	575	!
	576	CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz )
	577	CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
	578	#if defined key_vvl
	579	CALL wrk_alloc( jpi, jpj, jpk, ze3f )
	580	#endif
	581	!
[108]	582	IF( kt == nit000 ) THEN
	583	IF(lwp) WRITE(numout,*)
[455]	584	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	585	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[2715]	586	IF( .NOT.lk_vvl ) THEN
	587	ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
	588	IF( lk_mpp ) CALL mpp_sum ( ierr )
	589	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
	590	ENDIF
[1438]	591	ENDIF
[108]	592
[1438]	593	IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t)
[108]	594	DO jk = 1, jpk
	595	DO jj = 1, jpjm1
	596	DO ji = 1, jpim1
	597	ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	598	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) * 0.25
[2715]	599	IF( ze3f(ji,jj,jk) /= 0._wp ) ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
[108]	600	END DO
	601	END DO
	602	END DO
	603	CALL lbc_lnk( ze3f, 'F', 1. )
	604	ENDIF
	605
[2715]	606	zfac12 = 1._wp / 12._wp ! Local constant initialization
[216]	607
[108]	608
[455]	609	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
[108]	610	! ! ===============
	611	DO jk = 1, jpkm1 ! Horizontal slab
	612	! ! ===============
	613
	614	! Potential vorticity and horizontal fluxes
	615	! -----------------------------------------
[643]	616	SELECT CASE( kvor ) ! vorticity considered
[1438]	617	CASE ( 1 ) ! planetary vorticity (Coriolis)
	618	zwz(:,:) = ff(:,:) * ze3f(:,:,jk)
	619	CASE ( 2 ) ! relative vorticity
	620	zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
[643]	621	CASE ( 3 ) ! metric term
	622	DO jj = 1, jpjm1
	623	DO ji = 1, fs_jpim1 ! vector opt.
	624	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	625	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	626	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
	627	END DO
	628	END DO
[1516]	629	CALL lbc_lnk( zwz, 'F', 1. )
	630	CASE ( 4 ) ! total (relative + planetary vorticity)
[1438]	631	zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
[643]	632	CASE ( 5 ) ! total (coriolis + metric)
	633	DO jj = 1, jpjm1
	634	DO ji = 1, fs_jpim1 ! vector opt.
	635	zwz(ji,jj) = ( ff (ji,jj) &
	636	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	637	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	638	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	639	& ) * ze3f(ji,jj,jk)
	640	END DO
	641	END DO
[1516]	642	CALL lbc_lnk( zwz, 'F', 1. )
[455]	643	END SELECT
	644
[108]	645	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	646	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	647
	648	! Compute and add the vorticity term trend
	649	! ----------------------------------------
[1438]	650	jj = 2
	651	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[108]	652	DO ji = 2, jpi
	653	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	654	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	655	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	656	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	657	END DO
	658	DO jj = 3, jpj
[1694]	659	DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3
[108]	660	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	661	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	662	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	663	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	664	END DO
	665	END DO
	666	DO jj = 2, jpjm1
	667	DO ji = fs_2, fs_jpim1 ! vector opt.
	668	zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	669	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	670	zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	671	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
[455]	672	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
	673	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]	674	END DO
	675	END DO
	676	! ! ===============
	677	END DO ! End of slab
	678	! ! ===============
[3294]	679	CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz )
	680	CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
	681	#if defined key_vvl
	682	CALL wrk_dealloc( jpi, jpj, jpk, ze3f )
	683	#endif
[2715]	684	!
[3294]	685	IF( nn_timing == 1 ) CALL timing_stop('vor_een')
	686	!
[455]	687	END SUBROUTINE vor_een
[216]	688
	689
[2528]	690	SUBROUTINE dyn_vor_init
[3]	691	!!---------------------------------------------------------------------
[2528]	692	!! * ROUTINE dyn_vor_init *
[3]	693	!!
	694	!! ** Purpose : Control the consistency between cpp options for
[1438]	695	!! tracer advection schemes
[3]	696	!!----------------------------------------------------------------------
[2715]	697	INTEGER :: ioptio ! local integer
[3294]	698	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	699	!!
[1601]	700	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
[3]	701	!!----------------------------------------------------------------------
	702
[1601]	703	REWIND ( numnam ) ! Read Namelist namdyn_vor : Vorticity scheme options
	704	READ ( numnam, namdyn_vor )
[3]	705
[503]	706	IF(lwp) THEN ! Namelist print
[3]	707	WRITE(numout,*)
[2528]	708	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	709	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	710	WRITE(numout,*) ' Namelist namdyn_vor : oice of the vorticity term scheme'
[503]	711	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	712	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
	713	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
	714	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
[52]	715	ENDIF
	716
[3294]	717	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	718	! at angles with three ocean points and one land point
	719	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
	720	DO jk = 1, jpk
	721	DO jj = 2, jpjm1
	722	DO ji = 2, jpim1
	723	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
	724	fmask(ji,jj,jk) = 1._wp
	725	END DO
	726	END DO
	727	END DO
	728	!
	729	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	730	!
	731	ENDIF
	732
[503]	733	ioptio = 0 ! Control of vorticity scheme options
	734	IF( ln_dynvor_ene ) ioptio = ioptio + 1
	735	IF( ln_dynvor_ens ) ioptio = ioptio + 1
	736	IF( ln_dynvor_mix ) ioptio = ioptio + 1
	737	IF( ln_dynvor_een ) ioptio = ioptio + 1
	738	IF( lk_esopa ) ioptio = 1
	739
	740	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
	741
[643]	742	! ! Set nvor (type of scheme for vorticity)
[503]	743	IF( ln_dynvor_ene ) nvor = 0
	744	IF( ln_dynvor_ens ) nvor = 1
	745	IF( ln_dynvor_mix ) nvor = 2
	746	IF( ln_dynvor_een ) nvor = 3
	747	IF( lk_esopa ) nvor = -1
	748
[643]	749	! ! Set ncor, nrvm, ntot (type of vorticity)
	750	IF(lwp) WRITE(numout,*)
	751	ncor = 1
	752	IF( ln_dynadv_vec ) THEN
	753	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
	754	nrvm = 2
	755	ntot = 4
	756	ELSE
	757	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
	758	nrvm = 3
	759	ntot = 5
	760	ENDIF
	761
[503]	762	IF(lwp) THEN ! Print the choice
	763	WRITE(numout,*)
[643]	764	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
	765	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
	766	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
	767	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
[503]	768	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
[3]	769	ENDIF
[503]	770	!
[2528]	771	END SUBROUTINE dyn_vor_init
[3]	772
[503]	773	!!==============================================================================
[3]	774	END MODULE dynvor

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

Download in other formats: