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

Context Navigation

source: branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 3084

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

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

Download in other formats: