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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 40.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
[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
[643]	29	USE dynadv ! momentum advection (use ln_dynadv_vec value)
[503]	30	USE trdmod ! ocean dynamics trends
	31	USE trdmod_oce ! ocean variables trends
	32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	33	USE prtctl ! Print control
	34	USE in_out_manager ! I/O manager
[2715]	35	USE lib_mpp
[3]	36
	37	IMPLICIT NONE
	38	PRIVATE
	39
[2528]	40	PUBLIC dyn_vor ! routine called by step.F90
	41	PUBLIC dyn_vor_init ! routine called by opa.F90
[3]	42
[1601]	43	! !!* Namelist namdyn_vor: vorticity term
[32]	44	LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme
	45	LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme
	46	LOGICAL, PUBLIC :: ln_dynvor_mix = .FALSE. !: mixed scheme
[503]	47	LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme
[3]	48
[503]	49	INTEGER :: nvor = 0 ! type of vorticity trend used
[643]	50	INTEGER :: ncor = 1 ! coriolis
	51	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
	52	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
[455]	53
[3211]	54	!! * Control permutation of array indices
	55	# include "oce_ftrans.h90"
	56	# include "dom_oce_ftrans.h90"
	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	!!----------------------------------------------------------------------
[2715]	77	USE oce, ONLY: ztrdu => ta , ztrdv => sa ! (ta,sa) used as 3D workspace
[3211]	78	!! DCSE_NEMO: need additional directives for renamed module variables
	79	!FTRANS ztrdu :I :I :z
	80	!FTRANS ztrdv :I :I :z
	81
[2715]	82	!
[503]	83	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[455]	84	!!----------------------------------------------------------------------
[2715]	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	!
[455]	185	END SUBROUTINE dyn_vor
	186
	187
[643]	188	SUBROUTINE vor_ene( kt, kvor, pua, pva )
[455]	189	!!----------------------------------------------------------------------
	190	!! * ROUTINE vor_ene *
	191	!!
[3]	192	!! ** Purpose : Compute the now total vorticity trend and add it to
	193	!! the general trend of the momentum equation.
	194	!!
	195	!! ** Method : Trend evaluated using now fields (centered in time)
	196	!! and the Sadourny (1975) flux form formulation : conserves the
	197	!! horizontal kinetic energy.
	198	!! The trend of the vorticity term is given by:
[455]	199	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]	200	!! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ]
	201	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ]
	202	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	203	!! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ]
	204	!! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ]
	205	!! Add this trend to the general momentum trend (ua,va):
	206	!! (ua,va) = (ua,va) + ( voru , vorv )
	207	!!
	208	!! ** Action : - Update (ua,va) with the now vorticity term trend
[503]	209	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[3]	210	!! and planetary vorticity trends) ('key_trddyn')
	211	!!
[503]	212	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	213	!!----------------------------------------------------------------------
[2715]	214	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	215	USE wrk_nemo, ONLY: zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3 ! 2D workspace
	216	!
[643]	217	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	218	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	219	! ! =nrvm (relative vorticity or metric)
[3211]	220	!FTRANS pua :I :I :z
	221	!FTRANS pva :I :I :z
	222	!! DCSE_NEMO: work around a deficiency in ftrans
	223	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	224	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
	225	REAL(wp), INTENT(inout) :: pua(jpi,jpj,jpk) ! total u-trend
	226	REAL(wp), INTENT(inout) :: pva(jpi,jpj,jpk) ! total v-trend
[2715]	227	!
	228	INTEGER :: ji, jj, jk ! dummy loop indices
	229	REAL(wp) :: zx1, zy1, zfact2, zx2, zy2 ! local scalars
[3]	230	!!----------------------------------------------------------------------
	231
[2715]	232	IF( wrk_in_use(2, 1,2,3) ) THEN
	233	CALL ctl_stop('dyn:vor_ene: requested workspace arrays unavailable') ; RETURN
	234	ENDIF
	235
[52]	236	IF( kt == nit000 ) THEN
	237	IF(lwp) WRITE(numout,*)
[455]	238	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	239	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	240	ENDIF
[3]	241
[1438]	242	zfact2 = 0.5 * 0.5 ! Local constant initialization
[216]	243
[455]	244	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	245	! ! ===============
	246	DO jk = 1, jpkm1 ! Horizontal slab
	247	! ! ===============
[1438]	248	!
[3]	249	! Potential vorticity and horizontal fluxes
	250	! -----------------------------------------
[643]	251	SELECT CASE( kvor ) ! vorticity considered
	252	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	253	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	254	CASE ( 3 ) ! metric term
	255	DO jj = 1, jpjm1
	256	DO ji = 1, fs_jpim1 ! vector opt.
	257	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	258	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	259	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	260	END DO
	261	END DO
	262	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	263	CASE ( 5 ) ! total (coriolis + metric)
	264	DO jj = 1, jpjm1
	265	DO ji = 1, fs_jpim1 ! vector opt.
	266	zwz(ji,jj) = ( ff (ji,jj) &
	267	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	268	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	269	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
	270	& )
	271	END DO
	272	END DO
[455]	273	END SELECT
	274
	275	IF( ln_sco ) THEN
	276	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
[3]	277	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	278	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	279	ELSE
	280	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	281	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	282	ENDIF
	283
	284	! Compute and add the vorticity term trend
	285	! ----------------------------------------
	286	DO jj = 2, jpjm1
	287	DO ji = fs_2, fs_jpim1 ! vector opt.
	288	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	289	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	290	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	291	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
[455]	292	pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	293	pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[3]	294	END DO
	295	END DO
	296	! ! ===============
	297	END DO ! End of slab
	298	! ! ===============
[2715]	299	IF( wrk_not_released(2, 1,2,3) ) CALL ctl_stop('dyn:vor_ene: failed to release workspace arrays')
	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	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	336	USE wrk_nemo, ONLY: zwx => wrk_2d_4 , zwy => wrk_2d_5 , zwz => wrk_2d_6 , zww => wrk_2d_7 ! 2D workspace
	337	!
[503]	338	INTEGER, INTENT(in) :: kt ! ocean timestep index
[2715]	339	!
[1438]	340	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	341	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars
	342	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -
[3]	343	!!----------------------------------------------------------------------
	344
[2715]	345	IF( wrk_in_use(2, 4,5,6,7) ) THEN
	346	CALL ctl_stop('dyn:vor_mix: requested workspace arrays unavailable') ; RETURN
	347	ENDIF
	348
[52]	349	IF( kt == nit000 ) THEN
	350	IF(lwp) WRITE(numout,*)
[455]	351	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
	352	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	353	ENDIF
[3]	354
[1438]	355	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	356	zfact2 = 0.5 * 0.5
	357
[455]	358	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
[3]	359	! ! ===============
	360	DO jk = 1, jpkm1 ! Horizontal slab
	361	! ! ===============
[1438]	362	!
[3]	363	! Relative and planetary potential vorticity and horizontal fluxes
	364	! ----------------------------------------------------------------
[455]	365	IF( ln_sco ) THEN
[643]	366	IF( ln_dynadv_vec ) THEN
	367	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
	368	ELSE
	369	DO jj = 1, jpjm1
	370	DO ji = 1, fs_jpim1 ! vector opt.
	371	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	372	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	373	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
	374	END DO
	375	END DO
	376	ENDIF
[3]	377	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
	378	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	379	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	380	ELSE
[643]	381	IF( ln_dynadv_vec ) THEN
	382	zww(:,:) = rotn(:,:,jk)
	383	ELSE
	384	DO jj = 1, jpjm1
	385	DO ji = 1, fs_jpim1 ! vector opt.
	386	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	387	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	388	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
	389	END DO
	390	END DO
	391	ENDIF
	392	zwz(:,:) = ff (:,:)
[3]	393	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	394	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	395	ENDIF
	396
	397	! Compute and add the vorticity term trend
	398	! ----------------------------------------
	399	DO jj = 2, jpjm1
	400	DO ji = fs_2, fs_jpim1 ! vector opt.
	401	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	402	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	403	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	404	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	405	! enstrophy conserving formulation for relative vorticity term
	406	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
	407	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
	408	! energy conserving formulation for planetary vorticity term
	409	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	410	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[503]	411	! mixed vorticity trend added to the momentum trends
[3]	412	ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
	413	va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
	414	END DO
	415	END DO
	416	! ! ===============
	417	END DO ! End of slab
	418	! ! ===============
[2715]	419	IF( wrk_not_released(2, 4,5,6,7) ) CALL ctl_stop('dyn:vor_mix: failed to release workspace arrays')
	420	!
[455]	421	END SUBROUTINE vor_mix
[216]	422
	423
[643]	424	SUBROUTINE vor_ens( kt, kvor, pua, pva )
[3]	425	!!----------------------------------------------------------------------
[455]	426	!! * ROUTINE vor_ens *
[3]	427	!!
	428	!! ** Purpose : Compute the now total vorticity trend and add it to
	429	!! the general trend of the momentum equation.
	430	!!
	431	!! ** Method : Trend evaluated using now fields (centered in time)
	432	!! and the Sadourny (1975) flux FORM formulation : conserves the
	433	!! potential enstrophy of a horizontally non-divergent flow. the
	434	!! trend of the vorticity term is given by:
[455]	435	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
[3]	436	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
	437	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
	438	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	439	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
	440	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
	441	!! Add this trend to the general momentum trend (ua,va):
	442	!! (ua,va) = (ua,va) + ( voru , vorv )
	443	!!
	444	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
[503]	445	!! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
[3]	446	!! and planetary vorticity trends) ('key_trddyn')
	447	!!
[503]	448	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	449	!!----------------------------------------------------------------------
[2715]	450	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	451	USE wrk_nemo, ONLY: zwx => wrk_2d_4, zwy => wrk_2d_5, zwz => wrk_2d_6 ! 2D workspace
	452	!
[643]	453	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	454	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
	455	! ! =nrvm (relative vorticity or metric)
[3211]	456	!FTRANS pua :I :I :z
	457	!FTRANS pva :I :I :z
	458	!! DCSE_NEMO: work around a deficiency in ftrans
	459	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	460	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
	461	REAL(wp), INTENT(inout) :: pua(jpi,jpj,jpk) ! total u-trend
	462	REAL(wp), INTENT(inout) :: pva(jpi,jpj,jpk) ! total v-trend
[2715]	463	!
[503]	464	INTEGER :: ji, jj, jk ! dummy loop indices
	465	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
[3]	466	!!----------------------------------------------------------------------
	467
[2715]	468	IF( wrk_in_use(2, 4,5,6) ) THEN
	469	CALL ctl_stop('dyn:vor_ens: requested workspace arrays unavailable') ; RETURN
	470	END IF
	471
[52]	472	IF( kt == nit000 ) THEN
	473	IF(lwp) WRITE(numout,*)
[455]	474	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	475	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	476	ENDIF
[3]	477
[1438]	478	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	479
[455]	480	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	481	! ! ===============
	482	DO jk = 1, jpkm1 ! Horizontal slab
	483	! ! ===============
[1438]	484	!
[3]	485	! Potential vorticity and horizontal fluxes
	486	! -----------------------------------------
[643]	487	SELECT CASE( kvor ) ! vorticity considered
	488	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	489	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	490	CASE ( 3 ) ! metric term
	491	DO jj = 1, jpjm1
	492	DO ji = 1, fs_jpim1 ! vector opt.
	493	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	494	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	495	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	496	END DO
	497	END DO
	498	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	499	CASE ( 5 ) ! total (coriolis + metric)
	500	DO jj = 1, jpjm1
	501	DO ji = 1, fs_jpim1 ! vector opt.
	502	zwz(ji,jj) = ( ff (ji,jj) &
	503	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	504	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	505	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	506	& )
	507	END DO
	508	END DO
[455]	509	END SELECT
[1438]	510	!
[455]	511	IF( ln_sco ) THEN
[3]	512	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	513	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	514	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
	515	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
	516	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
[3]	517	END DO
	518	END DO
	519	ELSE
	520	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	521	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	522	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
	523	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
[3]	524	END DO
	525	END DO
	526	ENDIF
[1438]	527	!
[3]	528	! Compute and add the vorticity term trend
	529	! ----------------------------------------
	530	DO jj = 2, jpjm1
	531	DO ji = fs_2, fs_jpim1 ! vector opt.
[455]	532	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
[503]	533	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
[455]	534	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
[503]	535	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
[455]	536	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	537	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[3]	538	END DO
	539	END DO
	540	! ! ===============
	541	END DO ! End of slab
	542	! ! ===============
[2715]	543	IF( wrk_not_released(2, 4,5,6) ) CALL ctl_stop('dyn:vor_ens: failed to release workspace arrays')
	544	!
[455]	545	END SUBROUTINE vor_ens
[216]	546
	547
[643]	548	SUBROUTINE vor_een( kt, kvor, pua, pva )
[108]	549	!!----------------------------------------------------------------------
[455]	550	!! * ROUTINE vor_een *
[108]	551	!!
	552	!! ** Purpose : Compute the now total vorticity trend and add it to
	553	!! the general trend of the momentum equation.
	554	!!
	555	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	556	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	557	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	558	!! when horizontal divergence is zero (see the NEMO documentation)
	559	!! Add this trend to the general momentum trend (ua,va).
[108]	560	!!
	561	!! ** Action : - Update (ua,va) with the now vorticity term trend
[503]	562	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[108]	563	!! and planetary vorticity trends) ('key_trddyn')
	564	!!
[503]	565	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	566	!!----------------------------------------------------------------------
[2715]	567	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	568	USE wrk_nemo, ONLY: zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3 ! 2D workspace
	569	USE wrk_nemo, ONLY: ztnw => wrk_2d_4 , ztne => wrk_2d_5
	570	USE wrk_nemo, ONLY: ztsw => wrk_2d_6 , ztse => wrk_2d_7
	571	#if defined key_vvl
[3211]	572	!FTRANS ze3f :I :I :z
[2715]	573	USE wrk_nemo, ONLY: ze3f => wrk_3d_1 ! 3D workspace (lk_vvl=T)
	574	#endif
	575	!
[643]	576	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	577	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	578	! ! =nrvm (relative vorticity or metric)
[3211]	579	!FTRANS pua :I :I :z
	580	!FTRANS pva :I :I :z
	581	!! DCSE_NEMO: work around a deficiency in ftrans
	582	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	583	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
	584	REAL(wp), INTENT(inout) :: pua(jpi,jpj,jpk) ! total u-trend
	585	REAL(wp), INTENT(inout) :: pva(jpi,jpj,jpk) ! total v-trend
[218]	586	!!
[1438]	587	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	588	INTEGER :: ierr ! local integer
	589	REAL(wp) :: zfac12, zua, zva ! local scalars
	590	#if ! defined key_vvl
[3211]	591	!FTRANS ze3f :I :I :z
[2715]	592	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all
[1438]	593	#endif
[108]	594	!!----------------------------------------------------------------------
	595
[2715]	596	IF( wrk_in_use(2, 1,2,3,4,5,6,7) .OR. wrk_in_use(3, 1) ) THEN
	597	CALL ctl_stop('dyn:vor_een: requested workspace arrays unavailable') ; RETURN
	598	ENDIF
	599
[108]	600	IF( kt == nit000 ) THEN
	601	IF(lwp) WRITE(numout,*)
[455]	602	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	603	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[2715]	604	IF( .NOT.lk_vvl ) THEN
	605	ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
	606	IF( lk_mpp ) CALL mpp_sum ( ierr )
	607	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
	608	ENDIF
[1438]	609	ENDIF
[108]	610
[1438]	611	IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t)
[108]	612	DO jk = 1, jpk
	613	DO jj = 1, jpjm1
	614	DO ji = 1, jpim1
	615	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) &
	616	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) ) * 0.25
[2715]	617	IF( ze3f(ji,jj,jk) /= 0._wp ) ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
[108]	618	END DO
	619	END DO
	620	END DO
	621	CALL lbc_lnk( ze3f, 'F', 1. )
	622	ENDIF
	623
[2715]	624	zfac12 = 1._wp / 12._wp ! Local constant initialization
[216]	625
[455]	626	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
[108]	627	! ! ===============
	628	DO jk = 1, jpkm1 ! Horizontal slab
	629	! ! ===============
	630
	631	! Potential vorticity and horizontal fluxes
	632	! -----------------------------------------
[643]	633	SELECT CASE( kvor ) ! vorticity considered
[1438]	634	CASE ( 1 ) ! planetary vorticity (Coriolis)
	635	zwz(:,:) = ff(:,:) * ze3f(:,:,jk)
	636	CASE ( 2 ) ! relative vorticity
	637	zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
[643]	638	CASE ( 3 ) ! metric term
	639	DO jj = 1, jpjm1
	640	DO ji = 1, fs_jpim1 ! vector opt.
	641	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	642	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	643	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
	644	END DO
	645	END DO
[1516]	646	CALL lbc_lnk( zwz, 'F', 1. )
	647	CASE ( 4 ) ! total (relative + planetary vorticity)
[1438]	648	zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
[643]	649	CASE ( 5 ) ! total (coriolis + metric)
	650	DO jj = 1, jpjm1
	651	DO ji = 1, fs_jpim1 ! vector opt.
	652	zwz(ji,jj) = ( ff (ji,jj) &
	653	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	654	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	655	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	656	& ) * ze3f(ji,jj,jk)
	657	END DO
	658	END DO
[1516]	659	CALL lbc_lnk( zwz, 'F', 1. )
[455]	660	END SELECT
	661
[108]	662	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	663	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	664
	665	! Compute and add the vorticity term trend
	666	! ----------------------------------------
[1438]	667	jj = 2
	668	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[108]	669	DO ji = 2, jpi
	670	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	671	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	672	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	673	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	674	END DO
	675	DO jj = 3, jpj
[1694]	676	DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3
[108]	677	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	678	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	679	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	680	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	681	END DO
	682	END DO
	683	DO jj = 2, jpjm1
	684	DO ji = fs_2, fs_jpim1 ! vector opt.
	685	zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	686	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	687	zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	688	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
[455]	689	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
	690	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]	691	END DO
	692	END DO
	693	! ! ===============
	694	END DO ! End of slab
	695	! ! ===============
[2715]	696	IF( wrk_not_released(2, 1,2,3,4,5,6,7) .OR. &
	697	wrk_not_released(3, 1) ) CALL ctl_stop('dyn:vor_een: failed to release workspace arrays')
	698	!
[455]	699	END SUBROUTINE vor_een
[216]	700
	701
[2528]	702	SUBROUTINE dyn_vor_init
[3]	703	!!---------------------------------------------------------------------
[2528]	704	!! * ROUTINE dyn_vor_init *
[3]	705	!!
	706	!! ** Purpose : Control the consistency between cpp options for
[1438]	707	!! tracer advection schemes
[3]	708	!!----------------------------------------------------------------------
[2715]	709	INTEGER :: ioptio ! local integer
	710	!!
[1601]	711	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
[3]	712	!!----------------------------------------------------------------------
	713
[1601]	714	REWIND ( numnam ) ! Read Namelist namdyn_vor : Vorticity scheme options
	715	READ ( numnam, namdyn_vor )
[3]	716
[503]	717	IF(lwp) THEN ! Namelist print
[3]	718	WRITE(numout,*)
[2528]	719	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	720	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	721	WRITE(numout,*) ' Namelist namdyn_vor : oice of the vorticity term scheme'
[503]	722	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	723	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
	724	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
	725	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
[52]	726	ENDIF
	727
[503]	728	ioptio = 0 ! Control of vorticity scheme options
	729	IF( ln_dynvor_ene ) ioptio = ioptio + 1
	730	IF( ln_dynvor_ens ) ioptio = ioptio + 1
	731	IF( ln_dynvor_mix ) ioptio = ioptio + 1
	732	IF( ln_dynvor_een ) ioptio = ioptio + 1
	733	IF( lk_esopa ) ioptio = 1
	734
	735	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
	736
[643]	737	! ! Set nvor (type of scheme for vorticity)
[503]	738	IF( ln_dynvor_ene ) nvor = 0
	739	IF( ln_dynvor_ens ) nvor = 1
	740	IF( ln_dynvor_mix ) nvor = 2
	741	IF( ln_dynvor_een ) nvor = 3
	742	IF( lk_esopa ) nvor = -1
	743
[643]	744	! ! Set ncor, nrvm, ntot (type of vorticity)
	745	IF(lwp) WRITE(numout,*)
	746	ncor = 1
	747	IF( ln_dynadv_vec ) THEN
	748	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
	749	nrvm = 2
	750	ntot = 4
	751	ELSE
	752	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
	753	nrvm = 3
	754	ntot = 5
	755	ENDIF
	756
[503]	757	IF(lwp) THEN ! Print the choice
	758	WRITE(numout,*)
[643]	759	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
	760	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
	761	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
	762	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
[503]	763	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
[3]	764	ENDIF
[503]	765	!
[2528]	766	END SUBROUTINE dyn_vor_init
[3]	767
[503]	768	!!==============================================================================
[3]	769	END MODULE dynvor

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

Download in other formats: