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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 4291

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

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

Download in other formats: