Context Navigation

thermcell_flux2.F90 @ 1357

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

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

Download in other formats:

Original Format