MODULE lmdz_thermcell_flux2

! $Id: lmdz_thermcell_flux2.F90 5103 2024-07-23 13:29:36Z abarral $

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
      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

      REAL,SAVE :: fomass_max=0.5
      REAL,SAVE :: alphamax=0.7
!$OMP THREADPRIVATE(fomass_max,alphamax)

      logical check_debug,labort_physic

      character (len=20) :: 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>=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<=lmax(ig)) then
               if (entr_star(ig,l)>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))>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>=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<lmax(ig)) then
               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
            elseif(l==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)>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<lmax(ig)) then
               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
            elseif(l==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)<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>=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>=lalim(ig).and.l<=lmax(ig) &
      .and.(zw2(ig,l+1)>1.e-10).and.(zw2(ig,l)>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)>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>=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)>fm(ig,l)).and.(l>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>=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==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)>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>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)<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>=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)<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)<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>=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)>1.e-10) then
           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
           if ( fm(ig,l+1) > 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>=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==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)*fomass_max/ptimestep
            ddd=detr(ig,l)-eee
            if (eee>0.) then
                ncorecfm3=ncorecfm3+1
                entr(ig,l)=entr(ig,l)-eee
                if ( ddd>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==lmax(ig)) then
                      detr(ig,l)=fm(ig,l)+entr(ig,l)
                   else
                      if(l>=lmax(ig).and.0==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*,'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_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
END MODULE lmdz_thermcell_flux2
