MODULE lmdz_thermcell_flux2 ! ! $Id: lmdz_thermcell_flux2.f90 5512 2025-01-28 18:07:51Z ymeurdesoif $ ! CONTAINS SUBROUTINE thermcell_flux2(ngrid,nlay,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 !--------------------------------------------------------------------------- USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux, thermals_fomass_max, thermals_alphamax IMPLICIT NONE ! arguments INTEGER, intent(in) :: ngrid,nlay REAL, intent(in) :: ptimestep REAL, intent(in), dimension(ngrid,nlay) :: masse INTEGER, intent(in), dimension(ngrid) :: lalim,lmax REAL, intent(in), dimension(ngrid,nlay) :: alim_star,entr_star,detr_star REAL, intent(in), dimension(ngrid) :: f REAL, intent(in), dimension(ngrid,nlay) :: rhobarz REAL, intent(in), dimension(ngrid,nlay+1) :: zw2,zlev ! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les ! lignes faisant apparaitre zqla, zmax ... ! REAL, intent(in), dimension(ngrid) :: zmax(ngrid) ! enlever aussi zqla REAL, intent(in), dimension(ngrid,nlay) :: zqla ! not used integer, intent(in) :: lev_out, lunout1 REAL,intent(out), dimension(ngrid,nlay) :: entr,detr REAL,intent(out), dimension(ngrid,nlay+1) :: fm ! local INTEGER ig,l integer igout,lout REAL zfm integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 REAL f_old,ddd0,eee0,ddd,eee,zzz logical check_debug,labort_physic character (len=20), PARAMETER :: modname='thermcell_flux2' character (len=80) :: abort_message 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) ! Active uniquement si check_debug=.true. ou prt_level>=10 !------------------------------------------------------------------------- check_debug=.false..or.prt_level>=10 if (check_debug) then do l=1,nlay 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 = '' labort_physic=.true. CALL abort_physic (modname,abort_message,1) endif endif enddo enddo endif !------------------------------------------------------------------------- ! Multiplication par le flux de masse issu de la femreture !------------------------------------------------------------------------- do l=1,nlay 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,nlay 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,nlay 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,nlay,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,nlay 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,nlay 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,nlay 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,nlay 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,nlay labort_physic=.false. do ig=1,ngrid if (entr(ig,l)<0.) then labort_physic=.true. igout=ig lout=l endif enddo if (labort_physic) then print*,'N1 ig,l,entr',igout,lout,entr(igout,lout) abort_message = 'entr negatif' CALL abort_physic (modname,abort_message,1) endif do ig=1,ngrid 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 enddo labort_physic=.false. do ig=1,ngrid if (entr(ig,l).lt.0.) then labort_physic=.true. igout=ig endif enddo if (labort_physic) then ig=igout 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_physic (modname,abort_message,1) endif ! 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,nlay 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. ncorecfm2=ncorecfm2+1 endif enddo labort_physic=.false. do ig=1,ngrid if (detr(ig,l).lt.0.) then labort_physic=.true. igout=ig endif enddo if (labort_physic) then ig=igout 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_physic (modname,abort_message,1) endif ! 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,nlay do ig=1,ngrid if (zw2(ig,l+1).gt.1.e-10) then zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*thermals_alphamax if ( fm(ig,l+1) .gt. zfm) then f_old=fm(ig,l+1) fm(ig,l+1)=zfm 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,nlay,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 labort_physic=.false. do l=1,nlay-1 do ig=1,ngrid eee0=entr(ig,l) ddd0=detr(ig,l) eee=entr(ig,l)-masse(ig,l)*thermals_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 igout=ig lout=l labort_physic=.true. 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 if (labort_physic) then ig=igout l=lout 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*,'thermal_fomass_max',thermals_fomass_max print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*thermals_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_physic (modname,abort_message,1) endif 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,nlay,lunout1,igout,f,lmax,lalim, & ! & ptimestep,masse,entr,detr,fm,'fin') RETURN END SUBROUTINE thermcell_flux2 END MODULE lmdz_thermcell_flux2