source: LMDZ6/trunk/libf/dyn3dmem/qminimum_loc.F @ 5206

!
!     $Id: qminimum_loc.F 5199 2024-09-18 21:22:37Z fairhead $
!
      SUBROUTINE qminimum_loc( q,nqtot,deltap )
      USE parallel_lmdz
      USE infotrac, ONLY: niso, ntiso, iqIsoPha, tracers, addPhase,
     &                    isoCheck, min_qParent
      USE strings_mod, ONLY: strIdx
      IMPLICIT none
c
c  -- Objet : Traiter les valeurs trop petites (meme negatives)
c             pour l'eau vapeur et l'eau liquide
c
      include "dimensions.h"
      include "paramet.h"
      include "iniprint.h"
c
      INTEGER nqtot ! CRisi: on remplace nq par nqtot
      REAL q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
c
      LOGICAL, SAVE :: first=.TRUE.
      INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
c$OMP THREADPRIVATE(iq_vap, iq_liq, first)
      REAL, PARAMETER :: seuil_vap = 1.0e-10 ! seuil pour l'eau vapeur
      REAL, PARAMETER :: seuil_liq = 1.0e-11 ! seuil pour l'eau liquide
c
c  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
c            parametres seuil_vap, seuil_liq soient pareilles a celles 
c            qui  sont utilisees dans la routine    ADDFI       )
c     .................................................................
c
cDC iq_val and iq_liq are usable for q only, NOT for q_follow
c   and zx_defau_diag (crash if iq_val/liq==3) => vapor/liquid
c   water at hardcoded indices 1/2 in these variables
      INTEGER i, k, iq
      REAL zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe

      real zx_defau_diag(ijb_u:ije_u,llm,2) 
      real q_follow(ijb_u:ije_u,llm,2) 
c
      REAL SSUM
      EXTERNAL SSUM
c
      INTEGER imprim
      SAVE imprim
      DATA imprim /0/
c$OMP THREADPRIVATE(imprim)
      INTEGER ijb,ije
      INTEGER Index_pump(ij_end-ij_begin+1)
      INTEGER nb_pump
      INTEGER ixt
      INTEGER iso_verif_noNaN_nostop

c$OMP BARRIER

      !write(lunout,*) 'qminimum 52: entree'
      IF(first) THEN
         iq_vap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
         iq_liq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
         first = .FALSE.
      END IF
c
c Quand l'eau liquide est trop petite (ou negative), on prend
c l'eau vapeur de la meme couche et la convertit en eau liquide
c (sans changer la temperature !)
c

      call check_isotopes(q,ij_begin,ij_end,'qminimum 52')   

      ijb=ij_begin
      ije=ij_end

      DO k = 1, llm
c$OMP DO SCHEDULE(STATIC)
        DO i = ijb, ije
          zx_defau_diag(i,k,1)=0.0
          zx_defau_diag(i,k,2)=0.0
          q_follow(i,k,1)=q(i,k,iq_vap)
          q_follow(i,k,2)=q(i,k,iq_liq)
        ENDDO
c$OMP END DO NOWAIT
      ENDDO

      !write(lunout,*) 'qminimum 57'
      DO k = 1, llm
c$OMP DO SCHEDULE(STATIC)       
        DO i = ijb, ije
          if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then

            if (niso > 0) zx_defau_diag(i,k,2)=AMAX1
     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )

            q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
            q(i,k,iq_liq) = seuil_liq
          endif
        END DO
c$OMP END DO NOWAIT
      END DO

c
c Quand l'eau vapeur est trop faible (ou negative), on complete
c le defaut en prennant de l'eau vapeur de la couche au-dessous.
c
      !write(lunout,*) 'qminimum 81'
      DO k = llm, 2, -1
ccc      zx_abc = dpres(k) / dpres(k-1)
c$OMP DO SCHEDULE(STATIC)
        DO i = ijb, ije

          if ( seuil_vap - q(i,k,iq_vap) .gt. 0.d0 ) then

            if (niso > 0) zx_defau_diag(i,k,1)
     &           = AMAX1( seuil_vap - q(i,k,iq_vap), 0.0 )

            q(i,k-1,iq_vap) = q(i,k-1,iq_vap) - (seuil_vap
     &           -q(i,k,iq_vap)) * deltap(i,k)/deltap(i,k-1)
            q(i,k,iq_vap)   =  seuil_vap  

          endif
        ENDDO
c$OMP END DO NOWAIT
      ENDDO

c
c Quand il s'agit de la premiere couche au-dessus du sol, on
c doit imprimer un message d'avertissement (saturation possible).
c
      !write(lunout,*) 'qminimum 106'
      nb_pump=0
c$OMP DO SCHEDULE(STATIC)
      DO i = ijb, ije
         zx_pump(i) = AMAX1( 0.0, seuil_vap - q(i,1,iq_vap) )
         q(i,1,iq_vap)  = AMAX1( q(i,1,iq_vap), seuil_vap )
         IF (zx_pump(i) > 0.0) THEN
            nb_pump = nb_pump+1
            Index_pump(nb_pump)=i
         ENDIF
      ENDDO
c$OMP END DO NOWAIT 
!      pompe = SSUM(ije-ijb+1,zx_pump(ijb),1)

      IF (imprim.LE.100 .AND. nb_pump .GT. 0 ) THEN
         PRINT *, 'ATT!:on pompe de l eau au sol'
         DO i = 1, nb_pump
               imprim = imprim + 1
               PRINT*,'  en ',index_pump(i),zx_pump(index_pump(i))
         ENDDO
      ENDIF

      !write(lunout,*) 'qminimum 128'
      if (niso > 0) then
              !write(lunout,*) 'qminimum 140'
      ! CRisi: traiter de même les traceurs d'eau
      ! Mais il faut les prendre à l'envers pour essayer de conserver la
      ! masse.
      ! 1) pompage dans le sol  
      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
      ! rien ici et on croise les doigts pour que ça ne soit pas trop
      ! génant
      ! en fait, si, c'est genant quand les isotopes doivent eux même transporter des
      ! traceurs -> apporter aussi un peu d'isotopes... Combien?
      ! Essayer tnat/2 = -500 permil? C'est déjà mieux que -1000
      ! permil...
      ! pb: que faire pour les traceurs?
c$OMP DO SCHEDULE(STATIC)      
      DO i = ijb, ije
        if (zx_pump(i).gt.0.0) then
          q_follow(i,1,1)=q_follow(i,1,1)+zx_pump(i)
        endif !if (zx_pump(i).gt.0.0) then
      enddo !DO i = ijb, ije  
c$OMP END DO NOWAIT

      ! 2) transfert de vap vers les couches plus hautes
      !write(lunout,*) 'qminimum 158'
      do k=2,llm
c$OMP DO SCHEDULE(STATIC)      
        DO i = ijb, ije
          if (zx_defau_diag(i,k,1).gt.0.0) then             
              ! on ajoute la vapeur en k     
!              write(lunout,*) 'i,k,q_follow(i,k-1,ivap)=',
!     :                 i,k,q_follow(i,k-1,1)          
              if (q_follow(i,k-1,1).lt.min_qParent) then
                write(lunout,*) 'tmp qmin: on stoppe'
                write(lunout,*) 'zx_pump(i)=',zx_pump(i)
                write(lunout,*) 'q_follow(i,:,ivap)=',
     :                   q_follow(i,:,1)
                write(lunout,*) 'k=',k
                call abort_gcm("qminimum","not enough vapor",1)
              endif  
            do ixt=1,ntiso
!                write(lunout,*) 'qmin 168: ixt=',ixt
!                write(lunout,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
!     :             q(i,k,iqIsoPha(ixt,iq_vap))
!                write(lunout,*) 'zx_defau_diag(i,k,ivap)=',
!     :                  zx_defau_diag(i,k,1)
!                write(lunout,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
!     :                   q(i,k-1,iqIsoPha(ixt,iq_vap))     

               q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
     :           +zx_defau_diag(i,k,1)
     :           *q(i,k-1,iqIsoPha(ixt,iq_vap))/q_follow(i,k-1,1)
                
              if (isoCheck) then
                if(iso_verif_noNaN_nostop(q(i,k,iqIsoPha(ixt,iq_vap)),
     :                   'qminimum 155').eq.1) then
                   write(*,*) 'i,k,ixt=',i,k,ixt
                   write(*,*) 'q_follow(i,k-1,ivap)=',
     :                   q_follow(i,k-1,1)
                   write(*,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
     :                   q(i,k,iqIsoPha(ixt,iq_vap))
                   write(*,*) 'zx_defau_diag(i,k,ivap)=',
     :                   zx_defau_diag(i,k,1)
                   write(*,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
     :                   q(i,k-1,iqIsoPha(ixt,iq_vap))
                CALL abort_gcm("qminimum_loc","stopped",1)
                endif
              endif

              ! et on la retranche en k-1
               q(i,k-1,iqIsoPha(ixt,iq_vap)) =
     :            q(i,k-1,iqIsoPha(ixt,iq_vap))
     :              -zx_defau_diag(i,k,1)
     :              *deltap(i,k)/deltap(i,k-1)
     :              *q(i,k-1,iqIsoPha(ixt,iq_vap))
     :              /q_follow(i,k-1,1)

               if (isoCheck) then
                if (iso_verif_noNaN_nostop(
     :              q(i,k-1,iqIsoPha(ixt,iq_vap)),
     :                   'qminimum 175').eq.1) then
                   write(*,*) 'k,i,ixt=',k,i,ixt
                   write(*,*) 'q_follow(i,k-1,ivap)=',
     :                   q_follow(i,k-1,1)
                   write(*,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
     :                   q(i,k,iqIsoPha(ixt,iq_vap))
                   write(*,*) 'zx_defau_diag(i,k,ivap)=',
     :                   zx_defau_diag(i,k,1)
                   write(*,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
     :                   q(i,k-1,iqIsoPha(ixt,iq_vap))
                   CALL abort_gcm("qminimum_loc","stopped",1)
                endif
              endif 

              enddo !do ixt=1,niso
              q_follow(i,k,1)=   q_follow(i,k,1)
     :               +zx_defau_diag(i,k,1)
              q_follow(i,k-1,1)=   q_follow(i,k-1,1)
     :               -zx_defau_diag(i,k,1)
     :              *deltap(i,k)/deltap(i,k-1)
          endif !if (zx_defau_diag(i,k,1).gt.0.0) then
        enddo !DO i = 1, ip1jmp1        
c$OMP END DO NOWAIT
        enddo !do k=2,llm

        call check_isotopes(q,ijb,ije,'qminimum 168')
        
      
        ! 3) transfert d'eau de la vapeur au liquide
        !write(*,*) 'qminimum 164'
        do k=1,llm
c$OMP DO SCHEDULE(STATIC)
        DO i = ijb, ije
          if (zx_defau_diag(i,k,2).gt.0.0) then

              ! on ajoute eau liquide en k en k              
              do ixt=1,ntiso
               q(i,k,iqIsoPha(ixt,iq_liq))=q(i,k,iqIsoPha(ixt,iq_liq))
     :              +zx_defau_diag(i,k,2)
     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,1)
              ! et on la retranche à la vapeur en k
               q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
     :              -zx_defau_diag(i,k,2)
     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,1)   
              enddo !do ixt=1,niso
              q_follow(i,k,2)=   q_follow(i,k,2)
     :               +zx_defau_diag(i,k,2)
              q_follow(i,k,1)=   q_follow(i,k,1)
     :               -zx_defau_diag(i,k,2)
          endif !if (zx_defau_diag(i,k,1).gt.0.0) then
        enddo !DO i = ijb, ije
c$OMP END DO NOWAIT       
       enddo !do k=2,llm  

       call check_isotopes(q,ijb,ije,'qminimum 197')

      endif !if (niso > 0) then
      !write(*,*) 'qminimum 188'
c$OMP BARRIER

c
      RETURN
      END