Context Navigation

thermcell_flux2.F90 @ 972

Last change on this file since 972 was 972, checked in by lmdzadmin, 16 years ago
Version thermique FH/CRio Ajout tests cas physiques non pris en comptes et ajout/enleve prints Nouvelle routine thermcell_flux2.F90 IM
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 14.8 KB

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

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

Download in other formats:

Original Format