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/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 8197

Last change on this file since 8197 was 8197, checked in by glong, 7 years ago
Changed id's to be chars e.g. hpg to more easily identify output (and updated field_def.xml accordingly). Also rearranged scaling factors in dyn_vrt_dia subroutines in divcur.F90
File size: 42.1 KB

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

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

Download in other formats: