Index: LMDZ6/trunk/libf/phylmdiso/thermcell_flux2.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_flux2.F90	(revision 3939)
+++ LMDZ6/trunk/libf/phylmdiso/thermcell_flux2.F90	(revision 3940)
@@ -1,515 +1,1 @@
-!
-! $Id: thermcell_flux2.F90 3102 2017-12-03 20:27:42Z oboucher $
-!
-      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
-!---------------------------------------------------------------------------
-
-      USE print_control_mod, ONLY: prt_level
-      IMPLICIT NONE
-#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,lout
-      integer lev_out
-      integer lunout1
-
-      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.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,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 = ''
-                    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,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
-
-
-
-         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,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.
-               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,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
-      labort_physic=.false.
-      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
-                         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,klev,lunout1,igout,f,lmax,lalim, &
-!    &    ptimestep,masse,entr,detr,fm,'fin')
-
-
-      return
-      end
+link ../phylmd/thermcell_flux2.F90
