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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago
ORCA2_LIM_PISCES hybrid version update
Property svn:keywords set to `Id`
File size: 38.4 KB

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

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

Download in other formats: