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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 10473

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

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

Download in other formats: