! ! $Id: thermcell_flux2.F90 1299 2010-01-20 14:27:21Z abarral $ ! SUBROUTINE thermcell_flux2(ngrid,klev,ptimestep,masse, & & lalim,lmax,alim_star, & & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, & & detr,zqla,lev_out,lunout1,igout) !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout) !--------------------------------------------------------------------------- !thermcell_flux: deduction des flux !--------------------------------------------------------------------------- IMPLICIT NONE #include "iniprint.h" #include "thermcell.h" INTEGER ig,l INTEGER ngrid,klev REAL alim_star(ngrid,klev),entr_star(ngrid,klev) REAL detr_star(ngrid,klev) REAL zw2(ngrid,klev+1) REAL zlev(ngrid,klev+1) REAL masse(ngrid,klev) REAL ptimestep REAL rhobarz(ngrid,klev) REAL f(ngrid) INTEGER lmax(ngrid) INTEGER lalim(ngrid) REAL zqla(ngrid,klev) REAL zmax(ngrid) integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 REAL entr(ngrid,klev),detr(ngrid,klev) REAL fm(ngrid,klev+1) REAL zfm integer igout integer lev_out integer lunout1 REAL f_old,ddd0,eee0,ddd,eee,zzz REAL fomass_max,alphamax save fomass_max,alphamax character (len=20) :: modname='thermcell_flux2' character (len=80) :: abort_message fomass_max=0.5 alphamax=0.7 ncorecfm1=0 ncorecfm2=0 ncorecfm3=0 ncorecfm4=0 ncorecfm5=0 ncorecfm6=0 ncorecfm7=0 ncorecfm8=0 ncorecalpha=0 !initialisation fm(:,:)=0. if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_flux 0' write(lunout1,*) 'flux base ',f(igout) write(lunout1,*) 'lmax ',lmax(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) 'ig= ',igout write(lunout1,*) ' l E* A* D* ' write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) & & ,l=1,lmax(igout)) endif !------------------------------------------------------------------------- ! Verification de la nullite des entrainement et detrainement au dessus ! de lmax(ig) !------------------------------------------------------------------------- do l=1,klev do ig=1,ngrid if (l.le.lmax(ig)) then if (entr_star(ig,l).gt.1.) then print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig) print*,'entr_star(ig,l)',entr_star(ig,l) print*,'alim_star(ig,l)',alim_star(ig,l) print*,'detr_star(ig,l)',detr_star(ig,l) endif else if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig) print*,'entr_star(ig,l)',entr_star(ig,l) print*,'alim_star(ig,l)',alim_star(ig,l) print*,'detr_star(ig,l)',detr_star(ig,l) abort_message = '' CALL abort_gcm (modname,abort_message,1) endif endif enddo enddo !------------------------------------------------------------------------- ! Multiplication par le flux de masse issu de la femreture !------------------------------------------------------------------------- do l=1,klev entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l)) detr(:,l)=f(:)*detr_star(:,l) enddo if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_flux 1' write(lunout1,*) 'flux base ',f(igout) write(lunout1,*) 'lmax ',lmax(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) 'ig= ',igout write(lunout1,*) ' l E D W2' write(lunout1,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) & & ,zw2(igout,l+1),l=1,lmax(igout)) endif fm(:,1)=0. do l=1,klev do ig=1,ngrid if (l.lt.lmax(ig)) then fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) elseif(l.eq.lmax(ig)) then fm(ig,l+1)=0. detr(ig,l)=fm(ig,l)+entr(ig,l) else fm(ig,l+1)=0. endif enddo enddo ! Test provisoire : pour comprendre pourquoi on corrige plein de fois ! le cas fm6, on commence par regarder une premiere fois avant les ! autres corrections. do l=1,klev do ig=1,ngrid if (detr(ig,l).gt.fm(ig,l)) then ncorecfm8=ncorecfm8+1 ! igout=ig endif enddo enddo ! if (prt_level.ge.10) & ! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, & ! & ptimestep,masse,entr,detr,fm,'2 ') !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FH Version en cours de test; ! par rapport a thermcell_flux, on fait une grande boucle sur "l" ! et on modifie le flux avec tous les contr�les appliques d'affilee ! pour la meme couche ! Momentanement, on duplique le calcule du flux pour pouvoir comparer ! les flux avant et apres modif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do l=1,klev do ig=1,ngrid if (l.lt.lmax(ig)) then fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) elseif(l.eq.lmax(ig)) then fm(ig,l+1)=0. detr(ig,l)=fm(ig,l)+entr(ig,l) else fm(ig,l+1)=0. endif enddo !------------------------------------------------------------------------- ! Verification de la positivite des flux de masse !------------------------------------------------------------------------- ! do l=1,klev do ig=1,ngrid if (fm(ig,l+1).lt.0.) then ! print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1) ncorecfm1=ncorecfm1+1 fm(ig,l+1)=fm(ig,l) detr(ig,l)=entr(ig,l) endif enddo ! enddo if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) !------------------------------------------------------------------------- !Test sur fraca croissant !------------------------------------------------------------------------- if (iflag_thermals_optflux==0) then ! do l=1,klev do ig=1,ngrid if (l.ge.lalim(ig).and.l.le.lmax(ig) & & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then ! zzz est le flux en l+1 a frac constant zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1) & & /(rhobarz(ig,l)*zw2(ig,l)) if (fm(ig,l+1).gt.zzz) then detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz fm(ig,l+1)=zzz ncorecfm4=ncorecfm4+1 endif endif enddo ! enddo endif if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) !------------------------------------------------------------------------- !test sur flux de masse croissant !------------------------------------------------------------------------- if (iflag_thermals_optflux==0) then ! do l=1,klev do ig=1,ngrid if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then f_old=fm(ig,l+1) fm(ig,l+1)=fm(ig,l) detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) ncorecfm5=ncorecfm5+1 endif enddo ! enddo endif if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) !fin 1.eq.0 !------------------------------------------------------------------------- !detr ne peut pas etre superieur a fm !------------------------------------------------------------------------- if(1.eq.1) then ! do l=1,klev do ig=1,ngrid if (entr(ig,l)<0.) then print*,'N1 ig,l,entr',ig,l,entr(ig,l) abort_message = 'entr negatif' CALL abort_gcm (modname,abort_message,1) endif if (detr(ig,l).gt.fm(ig,l)) then ncorecfm6=ncorecfm6+1 detr(ig,l)=fm(ig,l) entr(ig,l)=fm(ig,l+1) ! Dans le cas ou on est au dessus de la couche d'alimentation et que le ! detrainement est plus fort que le flux de masse, on stope le thermique. !test:on commente ! if (l.gt.lalim(ig)) then ! lmax(ig)=l ! fm(ig,l+1)=0. ! entr(ig,l)=0. ! else ! ncorecfm7=ncorecfm7+1 ! endif endif if(l.gt.lmax(ig)) then detr(ig,l)=0. fm(ig,l+1)=0. entr(ig,l)=0. endif if (entr(ig,l).lt.0.) then print*,'ig,l,lmax(ig)',ig,l,lmax(ig) print*,'entr(ig,l)',entr(ig,l) print*,'fm(ig,l)',fm(ig,l) abort_message = 'probleme dans thermcell flux' CALL abort_gcm (modname,abort_message,1) endif enddo ! enddo endif if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) !------------------------------------------------------------------------- !fm ne peut pas etre negatif !------------------------------------------------------------------------- ! do l=1,klev do ig=1,ngrid if (fm(ig,l+1).lt.0.) then detr(ig,l)=detr(ig,l)+fm(ig,l+1) fm(ig,l+1)=0. ! print*,'fm2<0',l+1,lmax(ig) ncorecfm2=ncorecfm2+1 endif if (detr(ig,l).lt.0.) then print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig) print*,'detr(ig,l)',detr(ig,l) print*,'fm(ig,l)',fm(ig,l) abort_message = 'probleme dans thermcell flux' CALL abort_gcm (modname,abort_message,1) endif enddo ! enddo if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) !----------------------------------------------------------------------- !la fraction couverte ne peut pas etre superieure a 1 !----------------------------------------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FH Partie a revisiter. ! Il semble qu'etaient codees ici deux optiques dans le cas ! F/ (rho *w) > 1 ! soit limiter la hauteur du thermique en considerant que c'est ! la derniere chouche, soit limiter F a rho w. ! Dans le second cas, il faut en fait limiter a un peu moins ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin ! dans thermcell_main et qu'il semble de toutes facons deraisonable ! d'avoir des fractions de 1.. ! Ci dessous, et dans l'etat actuel, le premier des deux if est ! sans doute inutile. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! do l=1,klev do ig=1,ngrid if (zw2(ig,l+1).gt.1.e-10) then zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax if ( fm(ig,l+1) .gt. zfm) then f_old=fm(ig,l+1) fm(ig,l+1)=zfm ! zw2(ig,l+1)=0. ! zqla(ig,l+1)=0. detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) ! lmax(ig)=l+1 ! zmax(ig)=zlev(ig,lmax(ig)) ! print*,'alpha>1',l+1,lmax(ig) ncorecalpha=ncorecalpha+1 endif endif enddo ! enddo ! if (prt_level.ge.10) & & write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & & entr(igout,l),detr(igout,l),fm(igout,l+1) ! Fin de la grande boucle sur les niveaux verticaux enddo ! if (prt_level.ge.10) & ! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, & ! & ptimestep,masse,entr,detr,fm,'8 ') !----------------------------------------------------------------------- ! On fait en sorte que la quantite totale d'air entraine dans le ! panache ne soit pas trop grande comparee a la masse de la maille !----------------------------------------------------------------------- if (1.eq.1) then do l=1,klev-1 do ig=1,ngrid eee0=entr(ig,l) ddd0=detr(ig,l) eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep ddd=detr(ig,l)-eee if (eee.gt.0.) then ncorecfm3=ncorecfm3+1 entr(ig,l)=entr(ig,l)-eee if ( ddd.gt.0.) then ! l'entrainement est trop fort mais l'exces peut etre compense par une ! diminution du detrainement) detr(ig,l)=ddd else ! l'entrainement est trop fort mais l'exces doit etre compense en partie ! par un entrainement plus fort dans la couche superieure if(l.eq.lmax(ig)) then detr(ig,l)=fm(ig,l)+entr(ig,l) else if(l.ge.lmax(ig).and.0.eq.1) then print*,'ig,l',ig,l print*,'eee0',eee0 print*,'ddd0',ddd0 print*,'eee',eee print*,'ddd',ddd print*,'entr',entr(ig,l) print*,'detr',detr(ig,l) print*,'masse',masse(ig,l) print*,'fomass_max',fomass_max print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep print*,'ptimestep',ptimestep print*,'lmax(ig)',lmax(ig) print*,'fm(ig,l+1)',fm(ig,l+1) print*,'fm(ig,l)',fm(ig,l) abort_message = 'probleme dans thermcell_flux' CALL abort_gcm (modname,abort_message,1) endif entr(ig,l+1)=entr(ig,l+1)-ddd detr(ig,l)=0. fm(ig,l+1)=fm(ig,l)+entr(ig,l) detr(ig,l)=0. endif endif endif enddo enddo endif ! ! ddd=detr(ig)-entre !on s assure que tout s annule bien en zmax do ig=1,ngrid fm(ig,lmax(ig)+1)=0. entr(ig,lmax(ig))=0. detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig)) enddo !----------------------------------------------------------------------- ! Impression du nombre de bidouilles qui ont ete necessaires !----------------------------------------------------------------------- !IM 090508 beg ! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then ! ! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',& ! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', & ! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', & ! & ncorecfm6,'x fm6', & ! & ncorecfm7,'x fm7', & ! & ncorecfm8,'x fm8', & ! & ncorecalpha,'x alpha' ! endif !IM 090508 end ! if (prt_level.ge.10) & ! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, & ! & ptimestep,masse,entr,detr,fm,'fin') return end