Context Navigation

dynvor.F90

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

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

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

Context Navigation

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

Download in other formats: