source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_flux2.f90 @ 5300

Last change on this file since 5300 was 5268, checked in by abarral, 5 weeks ago

.f90 <-> .F90 depending on cpp key use

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 KB
Line 
1MODULE lmdz_thermcell_flux2
2!
3! $Id: lmdz_thermcell_flux2.f90 5268 2024-10-23 17:02:39Z abarral $
4!
5CONTAINS
6
7      SUBROUTINE thermcell_flux2(ngrid,nlay,ptimestep,masse, &
8     &       lalim,lmax,alim_star,  &
9     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
10     &       detr,zqla,lev_out,lunout1,igout)
11!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
12
13
14!---------------------------------------------------------------------------
15!thermcell_flux: deduction des flux
16!---------------------------------------------------------------------------
17
18      USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux
19      IMPLICIT NONE
20     
21! arguments
22      INTEGER, intent(in) :: ngrid,nlay
23      REAL, intent(in) :: ptimestep
24      REAL, intent(in), dimension(ngrid,nlay) :: masse
25      INTEGER, intent(in), dimension(ngrid) :: lalim,lmax
26      REAL, intent(in), dimension(ngrid,nlay) :: alim_star,entr_star,detr_star
27      REAL, intent(in), dimension(ngrid) :: f
28      REAL, intent(in), dimension(ngrid,nlay) :: rhobarz
29      REAL, intent(in), dimension(ngrid,nlay+1) :: zw2,zlev
30! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les
31!      lignes faisant apparaitre zqla, zmax ...
32!     REAL, intent(in), dimension(ngrid) :: zmax(ngrid)
33!     enlever aussi zqla
34      REAL, intent(in), dimension(ngrid,nlay) :: zqla  ! not used
35      integer, intent(in) :: lev_out, lunout1
36
37      REAL,intent(out), dimension(ngrid,nlay) :: entr,detr
38      REAL,intent(out), dimension(ngrid,nlay+1) :: fm
39
40! local
41      INTEGER ig,l
42      integer igout,lout
43      REAL zfm
44      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
45      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
46     
47
48      REAL f_old,ddd0,eee0,ddd,eee,zzz
49
50      REAL,SAVE :: fomass_max=0.5
51      REAL,SAVE :: alphamax=0.7
52!$OMP THREADPRIVATE(fomass_max,alphamax)
53
54      logical check_debug,labort_physic
55
56      character (len=20) :: modname='thermcell_flux2'
57      character (len=80) :: abort_message
58
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,nlay
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_physic=.true.
111                    CALL abort_physic (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,nlay
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,nlay
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
153
154! Test provisoire : pour comprendre pourquoi on corrige plein de fois
155! le cas fm6, on commence par regarder une premiere fois avant les
156! autres corrections.
157
158      do l=1,nlay
159         do ig=1,ngrid
160            if (detr(ig,l).gt.fm(ig,l)) then
161               ncorecfm8=ncorecfm8+1
162!              igout=ig
163            endif
164         enddo
165      enddo
166
167!      if (prt_level.ge.10) &
168!    &    call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
169!    &    ptimestep,masse,entr,detr,fm,'2  ')
170
171
172
173!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174! FH Version en cours de test;
175! par rapport a thermcell_flux, on fait une grande boucle sur "l"
176! et on modifie le flux avec tous les contr�les appliques d'affilee
177! pour la meme couche
178! Momentanement, on duplique le calcule du flux pour pouvoir comparer
179! les flux avant et apres modif
180!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
181
182      do l=1,nlay
183
184         do ig=1,ngrid
185            if (l.lt.lmax(ig)) then
186               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
187            elseif(l.eq.lmax(ig)) then
188               fm(ig,l+1)=0.
189               detr(ig,l)=fm(ig,l)+entr(ig,l)
190            else
191               fm(ig,l+1)=0.
192            endif
193         enddo
194
195
196!-------------------------------------------------------------------------
197! Verification de la positivite des flux de masse
198!-------------------------------------------------------------------------
199
200!     do l=1,nlay
201         do ig=1,ngrid
202            if (fm(ig,l+1).lt.0.) then
203!              print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
204                ncorecfm1=ncorecfm1+1
205               fm(ig,l+1)=fm(ig,l)
206               detr(ig,l)=entr(ig,l)
207            endif
208         enddo
209!     enddo
210
211      if (prt_level.ge.10) &
212     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
213     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
214
215!-------------------------------------------------------------------------
216!Test sur fraca croissant
217!-------------------------------------------------------------------------
218      if (iflag_thermals_optflux==0) then
219!     do l=1,nlay
220         do ig=1,ngrid
221          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
222     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
223!  zzz est le flux en l+1 a frac constant
224             zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
225     &                          /(rhobarz(ig,l)*zw2(ig,l))
226             if (fm(ig,l+1).gt.zzz) then
227                detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
228                fm(ig,l+1)=zzz
229                ncorecfm4=ncorecfm4+1
230             endif
231          endif
232        enddo
233!     enddo
234      endif
235
236      if (prt_level.ge.10) &
237     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
238     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
239
240
241!-------------------------------------------------------------------------
242!test sur flux de masse croissant
243!-------------------------------------------------------------------------
244      if (iflag_thermals_optflux==0) then
245!     do l=1,nlay
246         do ig=1,ngrid
247            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
248              f_old=fm(ig,l+1)
249              fm(ig,l+1)=fm(ig,l)
250              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
251               ncorecfm5=ncorecfm5+1
252            endif
253         enddo
254!     enddo
255      endif
256
257      if (prt_level.ge.10) &
258     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
259     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
260
261!fin 1.eq.0
262!-------------------------------------------------------------------------
263!detr ne peut pas etre superieur a fm
264!-------------------------------------------------------------------------
265
266      if(1.eq.1) then
267
268!     do l=1,nlay
269
270
271
272         labort_physic=.false.
273         do ig=1,ngrid
274            if (entr(ig,l)<0.) then
275               labort_physic=.true.
276               igout=ig
277               lout=l
278            endif
279         enddo
280
281         if (labort_physic) then
282            print*,'N1 ig,l,entr',igout,lout,entr(igout,lout)
283            abort_message = 'entr negatif'
284            CALL abort_physic (modname,abort_message,1)
285         endif
286
287         do ig=1,ngrid
288            if (detr(ig,l).gt.fm(ig,l)) then
289               ncorecfm6=ncorecfm6+1
290               detr(ig,l)=fm(ig,l)
291               entr(ig,l)=fm(ig,l+1)
292
293! Dans le cas ou on est au dessus de la couche d'alimentation et que le
294! detrainement est plus fort que le flux de masse, on stope le thermique.
295!test:on commente
296!               if (l.gt.lalim(ig)) then
297!                  lmax(ig)=l
298!                  fm(ig,l+1)=0.
299!                  entr(ig,l)=0.
300!               else
301!                  ncorecfm7=ncorecfm7+1
302!               endif
303            endif
304
305            if(l.gt.lmax(ig)) then
306               detr(ig,l)=0.
307               fm(ig,l+1)=0.
308               entr(ig,l)=0.
309            endif
310         enddo
311
312         labort_physic=.false.
313         do ig=1,ngrid
314            if (entr(ig,l).lt.0.) then
315               labort_physic=.true.
316               igout=ig
317            endif
318         enddo
319         if (labort_physic) then
320            ig=igout
321            print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
322            print*,'entr(ig,l)',entr(ig,l)
323            print*,'fm(ig,l)',fm(ig,l)
324            abort_message = 'probleme dans thermcell flux'
325            CALL abort_physic (modname,abort_message,1)
326         endif
327
328
329!     enddo
330      endif
331
332
333      if (prt_level.ge.10) &
334     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
335     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
336
337!-------------------------------------------------------------------------
338!fm ne peut pas etre negatif
339!-------------------------------------------------------------------------
340
341!     do l=1,nlay
342         do ig=1,ngrid
343            if (fm(ig,l+1).lt.0.) then
344               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
345               fm(ig,l+1)=0.
346               ncorecfm2=ncorecfm2+1
347            endif
348         enddo
349
350         labort_physic=.false.
351         do ig=1,ngrid
352            if (detr(ig,l).lt.0.) then
353               labort_physic=.true.
354               igout=ig
355            endif
356        enddo
357        if (labort_physic) then
358               ig=igout
359               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
360               print*,'detr(ig,l)',detr(ig,l)
361               print*,'fm(ig,l)',fm(ig,l)
362               abort_message = 'probleme dans thermcell flux'
363               CALL abort_physic (modname,abort_message,1)
364        endif
365!    enddo
366
367      if (prt_level.ge.10) &
368     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
369     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
370
371!-----------------------------------------------------------------------
372!la fraction couverte ne peut pas etre superieure a 1           
373!-----------------------------------------------------------------------
374
375
376!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
377! FH Partie a revisiter.
378! Il semble qu'etaient codees ici deux optiques dans le cas
379! F/ (rho *w) > 1
380! soit limiter la hauteur du thermique en considerant que c'est
381! la derniere chouche, soit limiter F a rho w.
382! Dans le second cas, il faut en fait limiter a un peu moins
383! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
384! dans thermcell_main et qu'il semble de toutes facons deraisonable
385! d'avoir des fractions de 1..
386! Ci dessous, et dans l'etat actuel, le premier des  deux if est
387! sans doute inutile.
388!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
389
390!    do l=1,nlay
391        do ig=1,ngrid
392           if (zw2(ig,l+1).gt.1.e-10) then
393           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
394           if ( fm(ig,l+1) .gt. zfm) then
395              f_old=fm(ig,l+1)
396              fm(ig,l+1)=zfm
397              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
398!             lmax(ig)=l+1
399!             zmax(ig)=zlev(ig,lmax(ig))
400!             print*,'alpha>1',l+1,lmax(ig)
401              ncorecalpha=ncorecalpha+1
402           endif
403           endif
404        enddo
405!    enddo
406!
407
408
409      if (prt_level.ge.10) &
410     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
411     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
412
413! Fin de la grande boucle sur les niveaux verticaux
414      enddo
415
416!      if (prt_level.ge.10) &
417!    &    call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
418!    &    ptimestep,masse,entr,detr,fm,'8  ')
419
420
421!-----------------------------------------------------------------------
422! On fait en sorte que la quantite totale d'air entraine dans le
423! panache ne soit pas trop grande comparee a la masse de la maille
424!-----------------------------------------------------------------------
425
426      if (1.eq.1) then
427      labort_physic=.false.
428      do l=1,nlay-1
429         do ig=1,ngrid
430            eee0=entr(ig,l)
431            ddd0=detr(ig,l)
432            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
433            ddd=detr(ig,l)-eee
434            if (eee.gt.0.) then
435                ncorecfm3=ncorecfm3+1
436                entr(ig,l)=entr(ig,l)-eee
437                if ( ddd.gt.0.) then
438!   l'entrainement est trop fort mais l'exces peut etre compense par une
439!   diminution du detrainement)
440                   detr(ig,l)=ddd
441                else
442!   l'entrainement est trop fort mais l'exces doit etre compense en partie
443!   par un entrainement plus fort dans la couche superieure
444                   if(l.eq.lmax(ig)) then
445                      detr(ig,l)=fm(ig,l)+entr(ig,l)
446                   else
447                      if(l.ge.lmax(ig).and.0.eq.1) then
448                         igout=ig
449                         lout=l
450                         labort_physic=.true.
451                      endif
452                      entr(ig,l+1)=entr(ig,l+1)-ddd
453                      detr(ig,l)=0.
454                      fm(ig,l+1)=fm(ig,l)+entr(ig,l)
455                      detr(ig,l)=0.
456                   endif
457                endif
458            endif
459         enddo
460      enddo
461      if (labort_physic) 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_physic (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,nlay,lunout1,igout,f,lmax,lalim, &
510!    &    ptimestep,masse,entr,detr,fm,'fin')
511
512
513 RETURN
514      end
515END MODULE lmdz_thermcell_flux2
Note: See TracBrowser for help on using the repository browser.