Context Navigation

epflux.F90 @ 2236

Last change on this file since 2236 was 1454, checked in by slebonnois, 10 years ago
SL: corrections in Titan and Venus tools
File size: 14.4 KB

Rev	Line
[1454]	1	subroutine epflux(iip1,jjp1,llm,indefini,latdeg,rbar &
[816]	2	,teta,u3d,v3d,w3d,press &
	3	,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
	4	! ,vpupbar,wpupbar,vpvpbar,wpvpbar,vptetapbar,wptetapbar &
	5	)
	6
	7	IMPLICIT NONE
	8	!=======================================================================
	9	!
	10	! Audrey Crespin Sept 2007
	11	!
	12	! source: Francois Forget Avril 1996
	13	!
	14	!
	15	! MODIF SLebonnois nov 2009
	16	!
	17	! Cette subroutine Calcule Les flux d'Eliassen Palm a partir
	18	! des composantes INTERPOLE EN NIVEAU DE PRESSION
	19	! Toutes les variables sont sur la grille de pression press(1:llm)
	20	! On suit la methode de Peixoto et Oort
	21	!
	22	! On est dans le plan meridien
	23	!
	24	! Calcule la divergence & composantes du flux d'EP
	25	! Calcule la deviation a la circulation meridienne moyenne
	26	! Calcule les TEM
	27	! Calcule les termes de Reynolds
	28	!
	29	!=======================================================================
	30	!-----------------------------------------------------------------------
	31	! declarations:
	32	! -------------
	33
	34	include "planet.h"
	35
	36	! --------
	37	! ARGUMENTS
	38	! ---------
	39
	40	! Inputs:
	41
	42	integer :: iip1,jjp1,llm
[1454]	43	real :: indefini,latdeg(jjp1),rbar(jjp1,llm)
[816]	44	REAL :: teta(iip1,jjp1,llm)
	45	REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
	46	REAL :: w3d(iip1,jjp1,llm)
	47	REAL :: press(llm)
	48
	49	! Outputs:
	50
	51	REAL :: epy2d(jjp1,llm),epz2d(jjp1,llm),divep2d(jjp1,llm)
	52	REAL :: vtem2d(jjp1,llm),wtem2d(jjp1,llm),acc_meridien2d(jjp1,llm)
	53	REAL :: vpupbar(jjp1,llm),wpupbar(jjp1,llm)
	54	REAL :: wpvpbar(jjp1,llm),vpvpbar(jjp1,llm)
	55	REAL :: vptetapbar(jjp1,llm),wptetapbar(jjp1,llm)
	56
	57	! Variables non zonales (vpup signifie : v^prim * u^prim)
	58	! -------------------
	59	REAL :: vpup(iip1,jjp1,llm), wpup(iip1,jjp1,llm)
	60	REAL :: vpvp(iip1,jjp1,llm), wpvp(iip1,jjp1,llm)
	61	REAL :: vptetap(iip1,jjp1,llm), wptetap(iip1,jjp1,llm)
	62
	63	! Moyennes zonales
	64	! -------------
	65	REAL :: ubar(jjp1,llm),vbar(jjp1,llm),tetabar(jjp1,llm)
	66	REAL :: wbar(jjp1,llm)
	67
	68
	69	REAL :: dtetabardp(jjp1,llm),dubardp(jjp1,llm)
	70	REAL :: vz(jjp1,llm),wlat(jjp1,llm),dvzdp(jjp1,llm)
	71	REAL :: dwlatdlat(jjp1,llm)
	72
	73	REAL :: depzdp(jjp1,llm)
	74
	75	! Autres
	76	! ------
[1454]	77	real :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
[816]	78	real :: pi
	79	REAL :: tetas(llm)
	80	logical :: peixoto
	81
	82	REAL :: f(jjp1), ducosdlat(jjp1,llm)
	83	REAL :: depycosdlat
	84	INTEGER :: i,j,l, n, iim, jjm
	85
	86	!-----------------------------------------------------------------------
	87
	88	iim = iip1-1
	89	jjm = jjp1-1
	90	pi = 2.*asin(1.)
	91
[1454]	92	! To avoid rounding effects at poles
	93	rlatu(:)=latdeg(:)*pi/180.00001
	94
[816]	95	! Calcul de moyennes zonales
	96	! ---------------------------
	97	call moyzon(iim,jjp1,llm,indefini,u3d,ubar)
	98	call moyzon(iim,jjp1,llm,indefini,v3d,vbar)
	99	call moyzon(iim,jjp1,llm,indefini,teta,tetabar)
	100	call moyzon(iim,jjp1,llm,indefini,w3d,wbar)
	101
	102
	103	! ------------------------------------------*--
	104	peixoto = .false.
	105
	106	if (peixoto) then
	107
	108	! MODIF speciale Peixoto: on utilise la moyenne
	109	! globale de teta: tetas(l) a la place de tetabar
	110	do l=1,llm
	111	tetas(l) = 0
	112	n= 0
	113	do j=2,jjm
[1454]	114	if (tetabar(j,l).ne.indefini) then
[816]	115	tetas(l) = tetas(l) + tetabar(j,l)
	116	n = n+1
	117	end if
	118	end do
	119	if (n.eq.0) stop 'bug dans elias '
	120	tetas(l) = tetas(l) / float(n)
	121
	122	do j=1,jjp1
	123	tetabar(j,l) = tetas(l)
	124	end do
	125	end do
	126	! write (,) 'tetas(l) ' , tetas
	127	! write(,)
	128
	129	endif
	130	! ------------------------------------------*--
	131
	132	! coriolis
	133	! --------
	134	do j=1,jjp1
	135	f(j) = 2omegasin(rlatu(j))
	136	enddo
	137
	138
	139	! Calcul des termes non zonaux
	140	! ----------------------------
	141	do l=1,llm
	142	do j=1,jjp1
	143	do i=1,iip1
	144	if (( v3d(i,j,l).eq.indefini).or. &
	145	( vbar(j,l).eq.indefini).or. &
	146	( teta(i,j,l).eq.indefini).or. &
	147	(tetabar(j,l).eq.indefini)) then
	148	vptetap(i,j,l)= indefini
	149	else
	150	vptetap(i,j,l)=(v3d(i,j,l)-vbar(j,l)) &
	151	*(teta(i,j,l)-tetabar(j,l))
	152	end if
	153	if (( w3d(i,j,l).eq.indefini).or. &
	154	( wbar(j,l).eq.indefini).or. &
	155	( teta(i,j,l).eq.indefini).or. &
	156	(tetabar(j,l).eq.indefini)) then
	157	wptetap(i,j,l)= indefini
	158	else
	159	wptetap(i,j,l)=(w3d(i,j,l)-wbar(j,l)) &
	160	*(teta(i,j,l)-tetabar(j,l))
	161	end if
	162	if ((v3d(i,j,l).eq.indefini).or. &
	163	( vbar(j,l).eq.indefini).or. &
	164	(u3d(i,j,l).eq.indefini).or. &
	165	( ubar(j,l).eq.indefini)) then
	166	vpup(i,j,l)= indefini
	167	else
	168	vpup(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(u3d(i,j,l)-ubar(j,l))
	169	end if
	170	if ((w3d(i,j,l).eq.indefini).or. &
	171	( wbar(j,l).eq.indefini).or. &
	172	(u3d(i,j,l).eq.indefini).or. &
	173	( ubar(j,l).eq.indefini)) then
	174	wpup(i,j,l)= indefini
	175	else
	176	wpup(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(u3d(i,j,l)-ubar(j,l))
	177	end if
	178	if ((v3d(i,j,l).eq.indefini).or. &
	179	( vbar(j,l).eq.indefini)) then
	180	vpvp(i,j,l)= indefini
	181	else
	182	vpvp(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(v3d(i,j,l)-vbar(j,l))
	183	end if
	184	if ((w3d(i,j,l).eq.indefini).or. &
	185	( wbar(j,l).eq.indefini).or. &
	186	(v3d(i,j,l).eq.indefini).or. &
	187	( vbar(j,l).eq.indefini)) then
	188	wpvp(i,j,l)= indefini
	189	else
	190	wpvp(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(v3d(i,j,l)-vbar(j,l))
	191	end if
	192	end do
	193	end do
	194	end do
	195
	196	! Moyennes zonales des termes non zonaux
	197	! --------------------------------------
	198
	199	! Termes de Reynolds
	200
	201	! flux zonaux de quantite de mouvement (vpupbar et wpupbar)
	202	! flux meridiens de quantite de mouvement (vpvpbar et wpvpbar)
	203	! flux meridien et vertical de chaleur (vptetapbar et wptetapbar)
	204
	205	call moyzon(iim,jjp1,llm,indefini,vptetap,vptetapbar)
	206	call moyzon(iim,jjp1,llm,indefini,wptetap,wptetapbar)
	207	call moyzon(iim,jjp1,llm,indefini,vpup,vpupbar)
	208	call moyzon(iim,jjp1,llm,indefini,wpup,wpupbar)
	209	call moyzon(iim,jjp1,llm,indefini,vpvp,vpvpbar)
	210	call moyzon(iim,jjp1,llm,indefini,wpvp,wpvpbar)
	211
	212	! write(,) 'vptetapbar(2,23) =' , vptetapbar(2,23)
	213	! write(,) 'wptetapbar(2,23) =' , wptetapbar(2,23)
	214	! write(,) 'vpupbar(2,23) =' , vpupbar(2,23)
	215
	216
	217	! Derivees d/dp des moyennes zonales
	218	! --------------------------------------
	219
	220	call dx_dp(jjp1,llm,indefini,press,ubar,dubardp)
	221	call dx_dp(jjp1,llm,indefini,press,tetabar,dtetabardp)
	222	! write (,) 'dtetabardp(l) (K/Pa)',(dtetabardp(6,l),l=1,llm)
	223	! write(,)
	224
	225	! Controle (on triche !) de dtetabardp
	226	do l=1,llm
	227	do j=1,jjp1
	228	if ((dtetabardp(j,l).gt.-1.e-3).and. &
	229	(dtetabardp(j,l).ne.indefini)) then
	230	dtetabardp(j,l) = -1.e-3
	231	! write(,) 'profil presque instable en j = ',j,' l= ',l
	232	end if
	233	end do
	234	end do
	235	! write(,)'dtetabardp(l) (K/Pa) corr ',(dtetabardp(6,l),l=1,llm)
	236	! write(,)
	237
	238
	239	! calculs intermediaires
	240	! ----------------------
	241	do l=1,llm
	242
[1454]	243	if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
[816]	244	ducosdlat(1,l) = ( ubar(2,l)cos(rlatu(2))-ubar(1,l)cos(rlatu(1)) ) &
	245	/( rlatu(2) - rlatu(1))
	246	else
	247	ducosdlat(1,l) = indefini
	248	end if
	249	do j=2,jjm
[1454]	250	if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
[816]	251	ducosdlat(j,l) = ( ubar(j+1,l)cos(rlatu(j+1))-ubar(j-1,l)cos(rlatu(j-1))) &
	252	/( rlatu(j+1) - rlatu(j-1))
	253	else
	254	ducosdlat(j,l) = indefini
	255	end if
	256	enddo
[1454]	257	if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
[816]	258	ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1)) &
	259	-ubar(jjm,l)*cos(rlatu(jjm)) ) &
	260	/( rlatu(jjp1) - rlatu(jjm))
	261	else
	262	ducosdlat(jjp1,l) = indefini
	263	end if
	264
	265	do j=1,jjp1
[1454]	266	if ((vptetapbar(j,l).ne.indefini).and. &
	267	(dtetabardp(j,l).ne.indefini)) then
[816]	268	vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
	269	wlat(j,l) = cos(rlatu(j))*vz(j,l)
	270	else
	271	vz(j,l) = indefini
	272	wlat(j,l) = indefini
	273	end if
	274	enddo
	275
	276	enddo
	277
	278	! Calcul des vecteurs flux Eliassen Palm et divergence
	279	! -----------------------------------------------------
	280	! ref: Read 1986
	281
	282	do l=1,llm
	283	epy2d(1,l) = indefini
	284	epz2d(1,l) = indefini
	285	do j=2,jjm
	286
	287	! Composante y
	288	! ------------
	289
[1454]	290	if ( (rbar(j,l).ne.indefini).and. &
	291	(vpupbar(j,l).ne.indefini).and. &
	292	(dubardp(j,l).ne.indefini).and. &
	293	(vz(j,l).ne.indefini)) then
[816]	294
	295	epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
	296	- dubardp(j,l) * vz(j,l) & ! terme non geo
	297	)
	298
	299	! if ((j.eq.1).and.(l.eq.23))write(,) 'epy2d(1,23) =' , epy2d(1,23)
	300	! if ((j.eq.2).and.(l.eq.23))write(,) 'epy2d(2,23) =' , epy2d(2,23)
	301	! if ((j.eq.3).and.(l.eq.23))write(,) 'epy2d(3,23) =' , epy2d(3,23)
	302	else
	303	epy2d(j,l) = indefini
	304	end if
	305
	306	! Composante z
	307	! ------------
	308
[1454]	309	if ( (rbar(j,l).ne.indefini).and. &
	310	(wpupbar(j,l).ne.indefini).and. &
	311	(ducosdlat(j,l).ne.indefini).and. &
	312	(vz(j,l).ne.indefini)) then
[816]	313
	314	epz2d(j,l) = rbar(j,l) * cos(rlatu(j))* &
	315	( vz(j,l) * ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) & ! terme non geo
	316	- vz(j,l) * f(j) &
	317	+ wpupbar(j,l) &
	318	)
	319	else
	320	epz2d(j,l) = indefini
	321	end if
	322
	323	! if ((j.eq.2).and.(l.eq.23))write(,) 'epz2d(2,23) =' , epz2d(2,23)
	324	end do
	325	epy2d(jjp1,l) = indefini
	326	epz2d(jjp1,l) = indefini
	327	end do
	328
	329	! Moyennes euleriennes transformees (vtem2d et wtem2d)
	330	! -----------------------------------------------
	331
	332	call dx_dp(jjp1,llm,indefini,press,vz,dvzdp)
	333
	334	do l=1,llm
	335
[1454]	336	if ( (wlat(2,l).ne.indefini) .and. &
	337	(wlat(1,l).ne.indefini) ) then
[816]	338	dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
	339	/(rlatu(2)-rlatu(1))
	340	else
	341	dwlatdlat(1,l)=indefini
	342	end if
	343	do j=2,jjm
[1454]	344	if ( (wlat(j+1,l).ne.indefini) .and. &
	345	(wlat(j-1,l).ne.indefini) ) then
[816]	346	dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
	347	/(rlatu(j+1)-rlatu(j-1))
	348	else
	349	dwlatdlat(j,l)=indefini
	350	end if
	351	enddo
[1454]	352	if ( (wlat(jjp1,l).ne.indefini) .and. &
	353	(wlat(jjm, l).ne.indefini) ) then
[816]	354	dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
	355	/(rlatu(jjp1)-rlatu(jjm))
	356	else
	357	dwlatdlat(jjp1,l)=indefini
	358	end if
	359
	360	do j=1,jjp1
[1454]	361	if ( (dvzdp(j,l).ne.indefini) &
	362	.and.(vbar(j,l).ne.indefini)) then
[816]	363	vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
	364	else
	365	vtem2d(j,l) = indefini
	366	endif
	367	enddo
	368
	369	do j=1,jjp1
[1454]	370	if ( (rbar(j,l).ne.indefini).and. &
	371	(dwlatdlat(j,l).ne.indefini).and. &
	372	(wbar(j,l).ne.indefini)) then
[816]	373	wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
	374	else
	375	wtem2d(j,l) = indefini
	376	endif
	377	enddo
	378
	379	enddo
	380
	381	! write(,) 'vtem2d(2,23) =' , vtem2d(2,23)
	382	! write(,) 'wtem2d(2,23) =' , wtem2d(2,23)
	383
	384	! print*,"OK"
	385
	386	! Deviation par rapport a la circulation meridienne moyenne
	387	! --------------------------------------------------------
	388	! U*.Nabla ubar
	389	do l=1,llm
	390	do j=1,jjp1
[1454]	391	if ( (rbar(j,l).ne.indefini) &
	392	.and.(vtem2d(j,l).ne.indefini) &
	393	.and.(wtem2d(j,l).ne.indefini) &
	394	.and.(ducosdlat(j,l).ne.indefini) &
	395	.and.(dubardp(j,l).ne.indefini)) then
[816]	396	acc_meridien2d(j,l) = &
	397	vtem2d(j,l)(ducosdlat(j,l)/(rbar(j,l)cos(rlatu(j))) &
	398	-f(j)) &
	399	+ wtem2d(j,l)*dubardp(j,l)
	400	else
	401	acc_meridien2d(j,l) = indefini
	402	endif
	403	enddo
	404	enddo
	405
	406	! write(,) 'acc_meridien2d(2,23) =' , acc_meridien2d(2,23)
	407
	408	! Divergence du flux
	409	! ------------------
	410	call dx_dp(jjp1,llm,indefini,press,epz2d,depzdp)
	411	! write(,) 'depzdp(2,23)',depzdp(j,l)
	412
	413	do l=1,llm
	414	divep2d(1,l) =0
	415	do j=2,jjm
[1454]	416	if ( (rbar(j,l).ne.indefini) .and. &
	417	(epy2d(j+1,l).ne.indefini) .and. &
	418	(epy2d(j-1,l).ne.indefini) .and. &
	419	(depzdp(j,l).ne.indefini) ) then
[816]	420
	421	depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1)) &
	422	- epy2d(j-1,l)*cos(rlatu(j-1)) ) &
	423	/(rlatu(j+1) - rlatu(j-1))
	424
	425	! DIVERGENCE DU FLUX :
	426
	427	divep2d(j,l) = depycosdlat/(rbar(j,l)*cos(rlatu(j)))+depzdp(j,l)
	428
	429	! if ((j.eq.2).and.(l.eq.23)) then
	430	! write(,) 'depycosdlat(2,23)', depycosdlat
	431	! write(,) 'depzdp(2,23)',depzdp(j,l)
	432	! write(,) 'divergence(2,23)',divep2d(j,l)
	433	! end if
	434
	435	! Wave driving (Du/Dt et non Dm/dt) :
	436	divep2d(j,l) = divep2d(j,l)/(rbar(j,l)*cos(rlatu(j)))
	437	else
	438	divep2d(j,l) = indefini
	439	end if
	440	end do
	441	divep2d(jjp1,l) =0
	442	end do
	443
	444	! write(,) ' fin :divep2d(2,23) =' , divep2d(2,23)
	445
	446	! Preparation pour sortie graphique (flux ``chapeau'' Peixoto et Oort p 391)
	447	! A VOIR...
	448	! do l=1,llm
	449	! do j=1,jjp1
	450	! if (epy2d(j,l).ne.indefini) &
	451	! epy2d(j,l) = 2pia0 (1./g0)cos(rlatu(j))*epy2d(j,l)
	452	! if (epz2d(j,l).ne.indefini) &
	453	! epz2d(j,l) = 2pia0*2 (1./g0)cos(rlatu(j))epz2d(j,l)
	454	! end do
	455	! end do
	456
	457	return
	458	end

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

Download in other formats:

Original Format