MODULE lmdz_thermcell_flux2

  ! $Id: lmdz_thermcell_flux2.F90 5119 2024-07-24 16:46:45Z 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
    USE lmdz_abort_physic, ONLY: abort_physic
    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