Context Navigation

thermcell_flux1.F90 @ 839

Last change on this file since 839 was 814, checked in by Laurent Fairhead, 18 years ago
Rajout de la physique utilisant les thermiques FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 6.4 KB

Rev	Line
[814]	1	SUBROUTINE thermcell_flux1(ngrid,klev,ptimestep,masse, &
	2	& lentr,lmax,alim_star, &
	3	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,alim,entr, &
	4	& detr,zqla,zmax,lev_out,lunout,igout)
	5
	6
	7	!---------------------------------------------------------------------------
	8	!thermcell_flux: deduction des flux
	9	!---------------------------------------------------------------------------
	10
	11	IMPLICIT NONE
	12
	13	INTEGER ig,k,l
	14	INTEGER ngrid,klev
	15
	16	REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
	17	REAL detr_star(ngrid,klev)
	18	REAL zw2(ngrid,klev+1)
	19	REAL zlev(ngrid,klev+1)
	20	REAL masse(ngrid,klev)
	21	REAL ptimestep
	22	REAL rhobarz(ngrid,klev)
	23	REAL f(ngrid)
	24	INTEGER lmax(ngrid)
	25	INTEGER lentr(ngrid)
	26	REAL zqla(ngrid,klev)
	27	REAL zmax(ngrid)
	28
	29
	30	integer lev_out,lunout,igout
	31
	32
	33	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
	34
	35
	36	REAL alim(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
	37	REAL fm(ngrid,klev+1)
	38
	39	REAL f_old,ddd,eee
	40
	41	REAL fracmax
	42	save fracmax
	43
	44	fracmax=0.5
	45
	46	ncorecfm1=0
	47	ncorecfm2=0
	48	ncorecfm3=0
	49	ncorecalpha=0
	50
	51	!initialisation
	52	do ig=1,ngrid
	53	do k=1,klev+1
	54	fm(ig,k)=0.
	55	enddo
	56	enddo
	57
	58	! Calcul de l alimentation
	59	do ig=1,ngrid
	60	do k=1,klev
	61	alim(ig,k)=f(ig)*alim_star(ig,k)
	62	enddo
	63	enddo
	64	!CR:test pour entrainer moins que la masse
	65	! do ig=1,ngrid
	66	! do l=1,lentr(ig)
	67	! if ((alim(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
	68	! alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
	69	! s -0.9*masse(ig,l)/ptimestep
	70	! alim(ig,l)=0.9*masse(ig,l)/ptimestep
	71	! endif
	72	! enddo
	73	! enddo
	74	! calcul du détrainement et de l entrainement
	75	do ig=1,ngrid
	76	do k=1,klev
	77	detr(ig,k)=f(ig)*detr_star(ig,k)
	78	enddo
	79	do k=1,klev
	80	entr(ig,k)=f(ig)*entr_star(ig,k)
	81	enddo
	82	enddo
	83
	84	if (lev_out.ge.10) then
	85	write(lunout,*) 'Dans thermcell_flux 1'
	86	write(lunout,*) 'flux base ',f(igout)
	87	write(lunout,*) 'lmax ',lmax(igout)
	88	endif
	89
	90	!
	91	! Calcul des flux
	92	do ig=1,ngrid
	93	do l=1,lmax(ig)
	94	!calcul du flux
	95	fm(ig,l+1)=fm(ig,l)+alim(ig,l)+entr(ig,l)-detr(ig,l)
	96
	97	if (lev_out.ge.10.and.ig.eq.igout) then
	98	write(lunout,'(i6,i4,3e15.5)') ig,l,entr(igout,l)+alim(igout,l),detr(igout,l) &
	99	& ,fm(igout,l+1)
	100	endif
	101
	102
	103	if (fm(ig,l+1).lt.0.) then
	104	! print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
	105	ncorecfm1=ncorecfm1+1
	106	fm(ig,l+1)=fm(ig,l)
	107	detr(ig,l)=alim(ig,l)+entr(ig,l)
	108	endif
	109	!verification des valeurs du flux
	110	!Test sur fraca croissant
	111	if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
	112	if ( (l.ge.lentr(ig)).and. &
	113	& ((fm(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt. &
	114	& (fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))))) then
	115	f_old=fm(ig,l+1)
	116	fm(ig,l+1)=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
	117	& /(rhobarz(ig,l)*zw2(ig,l))
	118	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	119	endif
	120	endif
	121	!test sur flux de masse croissant
	122	if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lentr(ig))) then
	123	f_old=fm(ig,l+1)
	124	fm(ig,l+1)=fm(ig,l)
	125	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	126	endif
	127	!detr ne peut pas etre superieur a fm
	128	if (detr(ig,l).gt.fm(ig,l)) then
	129	detr(ig,l)=fm(ig,l)
	130	entr(ig,l)=fm(ig,l+1)-alim(ig,l)
	131	endif
	132	!fm ne peut pas etre negatif
	133	if (fm(ig,l+1).lt.0.) then
	134	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
	135	fm(ig,l+1)=0.
	136	! print*,'fm2<0',l+1,lmax(ig)
	137	ncorecfm2=ncorecfm2+1
	138	endif
	139	!la fraction couverte ne peut pas etre superieure a 1
	140	if (zw2(ig,l+1).gt.1.e-10) then
	141	if ((((fm(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt. &
	142	& 1.)) then
	143	f_old=fm(ig,l+1)
	144	fm(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
	145	zw2(ig,l+1)=0.
	146	zqla(ig,l+1)=0.
	147	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	148	lmax(ig)=l+1
	149	zmax(ig)=zlev(ig,lmax(ig))
	150	! print*,'alpha>1',l+1,lmax(ig)
	151	ncorecalpha=ncorecalpha+1
	152	endif
	153	endif
	154	enddo
	155	enddo
	156	!
	157
	158	!-----------------------------------------------------------------------
	159	! On fait en sorte que la quantite totale d'air entraine dans le
	160	! panache ne soit pas trop grande comparee a la masse de la maille
	161	!-----------------------------------------------------------------------
	162
	163	if (1.eq.0) then
	164	do l=1,klev-1
	165	do ig=1,ngrid
	166	eee=entr(ig,l)+alim(ig,l)-masse(ig,l)*fracmax/ptimestep
	167	ddd=detr(ig,l)-eee
	168	if (eee.gt.0.) then
	169	ncorecfm3=ncorecfm3+1
	170	entr(ig,l)=entr(ig,l)-eee
	171	if ( ddd.gt.0.) then
	172	! l'entrainement est trop fort mais l'exces peut etre compense par une
	173	! diminution du detrainement)
	174	detr(ig,l)=ddd
	175	else
	176	! l'entrainement est trop fort mais l'exces doit etre compense en partie
	177	! par un entrainement plus fort dans la couche superieure
	178	if(l.ge.lmax(ig)) stop'probleme dans thermcell_flux'
	179	entr(ig,l+1)=entr(ig,l+1)-ddd
	180	detr(ig,l)=0.
	181	fm(ig,l+1)=fm(ig,l)+entr(ig,l)+alim(ig,l)
	182	detr(ig,l)=0.
	183	endif
	184	endif
	185	enddo
	186	enddo
	187	endif
	188	!
	189	! ddd=detr(ig)-entre
	190	!on s assure que tout s annule bien en zmax
	191	do ig=1,ngrid
	192	fm(ig,lmax(ig)+1)=0.
	193	entr(ig,lmax(ig))=0.
	194	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig)) &
	195	& +alim(ig,lmax(ig))
	196	enddo
	197
	198	!-----------------------------------------------------------------------
	199	! Impression du nombre de bidouilles qui ont ete necessaires
	200	!-----------------------------------------------------------------------
	201
	202	if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecalpha > 0 ) then
	203
	204	print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',ncorecfm2 &
	205	& ,'x fm2',ncorecfm3,'x fm3 et', ncorecalpha,'x alpha'
	206	endif
	207
	208	stop
	209	return
	210	end

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

Download in other formats:

Original Format