source: LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90 @ 1026

Last change on this file since 1026 was 1026, checked in by lmdzadmin, 16 years ago

Modifs thermiques
FH/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.9 KB
Line 
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#include "thermcell.h"
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"
162! et on modifie le flux avec tous les contrôles appliques d'affilee
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!-------------------------------------------------------------------------
204      if (iflag_thermals_optflux==0) then
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!-------------------------------------------------------------------------
230      if (iflag_thermals_optflux==0) then
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)
258               stop'entr negatif'
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)
287               stop'probleme dans thermcell flux'
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)
314               stop'probleme dans thermcell flux'
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)
415                         stop'probleme dans thermcell_flux'
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.