source: lmdz_wrf/trunk/WRFV3/lmdz/thermcell_flux2.F90 @ 1544

Last change on this file since 1544 was 142, checked in by lfita, 10 years ago

Incorporing all the work done in LMDZ_WRFmeas about:

1.- Chekcing for NaNs?
2.- Instabilities in thermcell (negative SQRT & robustness IF)

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