Context Navigation

thermcell_flux.F90 @ 985

Last change on this file since 985 was 938, checked in by lmdzadmin, 17 years ago
Enleve prints par defaut IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.1 KB

Rev	Line
[878]	1	!
	2	! $Header$
	3	!
	4
	5
	6	SUBROUTINE thermcell_flux(ngrid,klev,ptimestep,masse, &
	7	& lalim,lmax,alim_star, &
	8	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
[938]	9	& detr,zqla,zmax,lev_out,lunout1,igout)
[878]	10
	11
	12	!---------------------------------------------------------------------------
	13	!thermcell_flux: deduction des flux
	14	!---------------------------------------------------------------------------
	15
	16	IMPLICIT NONE
[938]	17	#include "iniprint.h"
[878]	18
	19	INTEGER ig,l
	20	INTEGER ngrid,klev
	21
	22	REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
	23	REAL detr_star(ngrid,klev)
	24	REAL zw2(ngrid,klev+1)
	25	REAL zlev(ngrid,klev+1)
	26	REAL masse(ngrid,klev)
	27	REAL ptimestep
	28	REAL rhobarz(ngrid,klev)
	29	REAL f(ngrid)
	30	INTEGER lmax(ngrid)
	31	INTEGER lalim(ngrid)
	32	REAL zqla(ngrid,klev)
	33	REAL zmax(ngrid)
	34
	35	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
	36	integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
	37
	38
	39	REAL entr(ngrid,klev),detr(ngrid,klev)
	40	REAL fm(ngrid,klev+1)
	41	REAL zfm
	42
	43	integer igout
	44	integer lev_out
[938]	45	integer lunout1
[878]	46
	47	REAL f_old,ddd0,eee0,ddd,eee,zzz
	48
	49	REAL fomass_max,alphamax
	50	save fomass_max,alphamax
	51
	52	fomass_max=0.5
	53	alphamax=0.7
	54
	55	ncorecfm1=0
	56	ncorecfm2=0
	57	ncorecfm3=0
	58	ncorecfm4=0
	59	ncorecfm5=0
	60	ncorecfm6=0
	61	ncorecfm7=0
	62	ncorecfm8=0
	63	ncorecalpha=0
	64
	65	!initialisation
	66	fm(:,:)=0.
	67
[938]	68	if (prt_level.ge.10) then
[878]	69	write(lunout,*) 'Dans thermcell_flux 0'
	70	write(lunout,*) 'flux base ',f(igout)
	71	write(lunout,*) 'lmax ',lmax(igout)
	72	write(lunout,*) 'lalim ',lalim(igout)
	73	write(lunout,*) 'ig= ',igout
	74	write(lunout,) ' l E A* D* '
	75	write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
	76	& ,l=1,lmax(igout))
	77	endif
	78
	79
	80	!-------------------------------------------------------------------------
	81	! Verification de la nullite des entrainement et detrainement au dessus
	82	! de lmax(ig)
	83	!-------------------------------------------------------------------------
[938]	84	if ( prt_level > 1 ) THEN
[878]	85	do l=1,klev
	86	do ig=1,ngrid
	87	if (l.le.lmax(ig)) then
	88	if (entr_star(ig,l).gt.1.) then
	89	print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
	90	print*,'entr_star(ig,l)',entr_star(ig,l)
	91	print*,'alim_star(ig,l)',alim_star(ig,l)
	92	print*,'detr_star(ig,l)',detr_star(ig,l)
	93	! stop
	94	endif
	95	else
	96	if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
	97	print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
	98	print*,'entr_star(ig,l)',entr_star(ig,l)
	99	print*,'alim_star(ig,l)',alim_star(ig,l)
	100	print*,'detr_star(ig,l)',detr_star(ig,l)
	101	stop
	102	endif
	103	endif
	104	enddo
	105	enddo
[938]	106	endif !( prt_level > 1 ) THEN
[878]	107	!-------------------------------------------------------------------------
	108	! Multiplication par le flux de masse issu de la femreture
	109	!-------------------------------------------------------------------------
	110
	111	do l=1,klev
	112	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
	113	detr(:,l)=f(:)*detr_star(:,l)
	114	enddo
	115
[938]	116	if (prt_level.ge.10) then
[878]	117	write(lunout,*) 'Dans thermcell_flux 1'
	118	write(lunout,*) 'flux base ',f(igout)
	119	write(lunout,*) 'lmax ',lmax(igout)
	120	write(lunout,*) 'lalim ',lalim(igout)
	121	write(lunout,*) 'ig= ',igout
	122	write(lunout,*) ' l E D W2'
	123	write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
	124	& ,zw2(igout,l+1),l=1,lmax(igout))
	125	endif
	126
	127	fm(:,1)=0.
	128	do l=1,klev
	129	do ig=1,ngrid
	130	if (l.lt.lmax(ig)) then
	131	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	132	elseif(l.eq.lmax(ig)) then
	133	fm(ig,l+1)=0.
	134	detr(ig,l)=fm(ig,l)+entr(ig,l)
	135	else
	136	fm(ig,l+1)=0.
	137	endif
	138	enddo
	139	enddo
	140
	141
	142
	143	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
	144	! le cas fm6, on commence par regarder une premiere fois avant les
	145	! autres corrections.
	146
	147	do l=1,klev
	148	do ig=1,ngrid
	149	if (detr(ig,l).gt.fm(ig,l)) then
	150	ncorecfm8=ncorecfm8+1
	151	! igout=ig
	152	endif
	153	enddo
	154	enddo
	155
[938]	156	if (prt_level.ge.10) &
[878]	157	& call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
	158	& ptimestep,masse,entr,detr,fm,'2 ')
	159
	160
	161
	162	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	163	! FH Version en cours de test;
	164	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
	165	! et on modifie le flux avec tous les contrôles appliques d'affilee
	166	! pour la meme couche
	167	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
	168	! les flux avant et apres modif
	169	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	170
	171	do l=1,klev
	172
	173	do ig=1,ngrid
	174	if (l.lt.lmax(ig)) then
	175	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	176	elseif(l.eq.lmax(ig)) then
	177	fm(ig,l+1)=0.
	178	detr(ig,l)=fm(ig,l)+entr(ig,l)
	179	else
	180	fm(ig,l+1)=0.
	181	endif
	182	enddo
	183
	184
	185	!-------------------------------------------------------------------------
	186	! Verification de la positivite des flux de masse
	187	!-------------------------------------------------------------------------
	188
	189	! do l=1,klev
	190	do ig=1,ngrid
	191	if (fm(ig,l+1).lt.0.) then
	192	! print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
	193	ncorecfm1=ncorecfm1+1
	194	fm(ig,l+1)=fm(ig,l)
	195	detr(ig,l)=entr(ig,l)
	196	endif
	197	enddo
	198	! enddo
	199
[938]	200	if (prt_level.ge.10) &
[878]	201	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	202	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	203
	204	!-------------------------------------------------------------------------
	205	!Test sur fraca croissant
	206	!-------------------------------------------------------------------------
	207
	208
[883]	209	if (1.eq.1) then
[878]	210	! do l=1,klev
	211	do ig=1,ngrid
	212	if (l.ge.lalim(ig).and.l.le.lmax(ig) &
	213	& .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
	214	! zzz est le flux en l+1 a frac constant
	215	zzz=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
	216	& /(rhobarz(ig,l)*zw2(ig,l))
	217	if (fm(ig,l+1).gt.zzz) then
	218	detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
	219	fm(ig,l+1)=zzz
	220	ncorecfm4=ncorecfm4+1
	221	endif
	222	endif
	223	enddo
	224	! enddo
	225
[938]	226	if (prt_level.ge.10) &
[878]	227	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	228	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	229	else
[938]	230	if (l.eq.1) then
[878]	231	print*,'Test sur les fractions croissantes inhibe dans thermcell_flux2'
[938]	232	endif
[878]	233	endif
	234
	235
	236	!-------------------------------------------------------------------------
	237	!test sur flux de masse croissant
	238	!-------------------------------------------------------------------------
	239
	240	! do l=1,klev
	241	do ig=1,ngrid
	242	if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
	243	f_old=fm(ig,l+1)
	244	fm(ig,l+1)=fm(ig,l)
	245	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	246	ncorecfm5=ncorecfm5+1
	247	endif
	248	enddo
	249	! enddo
	250
[938]	251	if (prt_level.ge.10) &
[878]	252	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	253	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	254
	255	!-------------------------------------------------------------------------
	256	!detr ne peut pas etre superieur a fm
	257	!-------------------------------------------------------------------------
	258
	259	if(1.eq.1) then
	260
	261	! do l=1,klev
	262	do ig=1,ngrid
	263	if (entr(ig,l)<0.) then
	264	print*,'N1 ig,l,entr',ig,l,entr(ig,l)
	265	stop'entr negatif'
	266	endif
	267	if (detr(ig,l).gt.fm(ig,l)) then
	268	ncorecfm6=ncorecfm6+1
	269	detr(ig,l)=fm(ig,l)
	270	! entr(ig,l)=fm(ig,l+1)
	271
	272	! Dans le cas ou on est au dessus de la couche d'alimentation et que le
	273	! detrainement est plus fort que le flux de masse, on stope le thermique.
	274	if (l.gt.lalim(ig)) then
	275	lmax(ig)=l
	276	fm(ig,l+1)=0.
	277	entr(ig,l)=0.
	278	else
	279	ncorecfm7=ncorecfm7+1
	280	endif
	281	endif
	282
	283	if(l.gt.lmax(ig)) then
	284	detr(ig,l)=0.
	285	fm(ig,l+1)=0.
	286	entr(ig,l)=0.
	287	endif
	288
	289	if (entr(ig,l).lt.0.) then
	290	print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
	291	print*,'entr(ig,l)',entr(ig,l)
	292	print*,'fm(ig,l)',fm(ig,l)
	293	stop'probleme dans thermcell flux'
	294	endif
	295	enddo
	296	! enddo
	297	endif
	298
	299
[938]	300	if (prt_level.ge.10) &
[878]	301	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	302	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	303
	304	!-------------------------------------------------------------------------
	305	!fm ne peut pas etre negatif
	306	!-------------------------------------------------------------------------
	307
	308	! do l=1,klev
	309	do ig=1,ngrid
	310	if (fm(ig,l+1).lt.0.) then
	311	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
	312	fm(ig,l+1)=0.
	313	! print*,'fm2<0',l+1,lmax(ig)
	314	ncorecfm2=ncorecfm2+1
	315	endif
	316	if (detr(ig,l).lt.0.) then
	317	print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
	318	print*,'detr(ig,l)',detr(ig,l)
	319	print*,'fm(ig,l)',fm(ig,l)
	320	stop'probleme dans thermcell flux'
	321	endif
	322	enddo
	323	! enddo
	324
[938]	325	if (prt_level.ge.10) &
[878]	326	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	327	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	328
	329	!-----------------------------------------------------------------------
	330	!la fraction couverte ne peut pas etre superieure a 1
	331	!-----------------------------------------------------------------------
	332
	333
	334	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	335	! FH Partie a revisiter.
	336	! Il semble qu'etaient codees ici deux optiques dans le cas
	337	! F/ (rho *w) > 1
	338	! soit limiter la hauteur du thermique en considerant que c'est
	339	! la derniere chouche, soit limiter F a rho w.
	340	! Dans le second cas, il faut en fait limiter a un peu moins
	341	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
	342	! dans thermcell_main et qu'il semble de toutes facons deraisonable
	343	! d'avoir des fractions de 1..
	344	! Ci dessous, et dans l'etat actuel, le premier des deux if est
	345	! sans doute inutile.
	346	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	347
	348	! do l=1,klev
	349	do ig=1,ngrid
	350	if (zw2(ig,l+1).gt.1.e-10) then
	351	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
	352	if ( fm(ig,l+1) .gt. zfm) then
	353	f_old=fm(ig,l+1)
	354	fm(ig,l+1)=zfm
	355	! zw2(ig,l+1)=0.
	356	! zqla(ig,l+1)=0.
	357	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	358	! lmax(ig)=l+1
	359	! zmax(ig)=zlev(ig,lmax(ig))
	360	! print*,'alpha>1',l+1,lmax(ig)
	361	ncorecalpha=ncorecalpha+1
	362	endif
	363	endif
	364	enddo
	365	! enddo
	366	!
	367
	368
[938]	369	if (prt_level.ge.10) &
[878]	370	& write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	371	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	372
	373	! Fin de la grande boucle sur les niveaux verticaux
	374	enddo
	375
[938]	376	if (prt_level.ge.10) &
[878]	377	& call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
	378	& ptimestep,masse,entr,detr,fm,'8 ')
	379
	380
	381	!-----------------------------------------------------------------------
	382	! On fait en sorte que la quantite totale d'air entraine dans le
	383	! panache ne soit pas trop grande comparee a la masse de la maille
	384	!-----------------------------------------------------------------------
	385
	386	if (1.eq.1) then
	387	do l=1,klev-1
	388	do ig=1,ngrid
	389	eee0=entr(ig,l)
	390	ddd0=detr(ig,l)
	391	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
	392	ddd=detr(ig,l)-eee
	393	if (eee.gt.0.) then
	394	ncorecfm3=ncorecfm3+1
	395	entr(ig,l)=entr(ig,l)-eee
	396	if ( ddd.gt.0.) then
	397	! l'entrainement est trop fort mais l'exces peut etre compense par une
	398	! diminution du detrainement)
	399	detr(ig,l)=ddd
	400	else
	401	! l'entrainement est trop fort mais l'exces doit etre compense en partie
	402	! par un entrainement plus fort dans la couche superieure
	403	if(l.eq.lmax(ig)) then
	404	detr(ig,l)=fm(ig,l)+entr(ig,l)
	405	else
	406	if(l.ge.lmax(ig).and.0.eq.1) then
	407	print*,'ig,l',ig,l
	408	print*,'eee0',eee0
	409	print*,'ddd0',ddd0
	410	print*,'eee',eee
	411	print*,'ddd',ddd
	412	print*,'entr',entr(ig,l)
	413	print*,'detr',detr(ig,l)
	414	print*,'masse',masse(ig,l)
	415	print*,'fomass_max',fomass_max
	416	print,'masse(ig,l)fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
	417	print*,'ptimestep',ptimestep
	418	print*,'lmax(ig)',lmax(ig)
	419	print*,'fm(ig,l+1)',fm(ig,l+1)
	420	print*,'fm(ig,l)',fm(ig,l)
	421	stop'probleme dans thermcell_flux'
	422	endif
	423	entr(ig,l+1)=entr(ig,l+1)-ddd
	424	detr(ig,l)=0.
	425	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
	426	detr(ig,l)=0.
	427	endif
	428	endif
	429	endif
	430	enddo
	431	enddo
	432	endif
	433	!
	434	! ddd=detr(ig)-entre
	435	!on s assure que tout s annule bien en zmax
	436	do ig=1,ngrid
	437	fm(ig,lmax(ig)+1)=0.
	438	entr(ig,lmax(ig))=0.
	439	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
	440	enddo
	441
	442	!-----------------------------------------------------------------------
	443	! Impression du nombre de bidouilles qui ont ete necessaires
	444	!-----------------------------------------------------------------------
	445
	446	if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
[938]	447	if (prt_level.ge.10) then
[878]	448	print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
	449	& ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
	450	& ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
	451	& ncorecfm6,'x fm6', &
	452	& ncorecfm7,'x fm7', &
	453	& ncorecfm8,'x fm8', &
	454	& ncorecalpha,'x alpha'
[938]	455	endif
[878]	456	endif
	457
[938]	458	if (prt_level.ge.10) &
[878]	459	& call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
	460	& ptimestep,masse,entr,detr,fm,'fin')
	461
	462
	463	return
	464	end
	465
	466	subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
	467	& ptimestep,masse,entr,detr,fm,descr)
	468
	469	implicit none
	470
	471	integer ngrid,klev,lunout,igout,l,lm
	472
	473	integer lmax(klev),lalim(klev)
	474	real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
	475	real fm(ngrid,klev+1),f(ngrid)
	476
	477	character*3 descr
	478
	479	lm=lmax(igout)+5
	480	if(lm.gt.klev) lm=klev
	481
	482	print*,'Impression jusque lm=',lm
	483
	484	write(lunout,*) 'Dans thermcell_flux '//descr
	485	write(lunout,*) 'flux base ',f(igout)
	486	write(lunout,*) 'lmax ',lmax(igout)
	487	write(lunout,*) 'lalim ',lalim(igout)
	488	write(lunout,*) 'ig= ',igout
	489	write(lunout,'(a3,4a14)') 'l','M','E','D','F'
	490	write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, &
	491	& entr(igout,l),detr(igout,l) &
	492	& ,fm(igout,l+1),l=1,lm)
	493
	494
	495	do l=lmax(igout)+1,klev
	496	if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then
	497	print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout)
	498	print*,'entr(igout,l)',entr(igout,l)
	499	print*,'detr(igout,l)',detr(igout,l)
	500	print*,'fm(igout,l)',fm(igout,l)
	501	stop
	502	endif
	503	enddo
	504
	505	return
	506	end
	507

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

Download in other formats:

Original Format