source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/thermcell_flux2.F90 @ 1798

Last change on this file since 1798 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

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