Context Navigation

thermcell_flux1.F90 @ 832

Last change on this file since 832 was 814, checked in by Laurent Fairhead, 17 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

Line
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