! $Id: qminimum_loc.f90 5113 2024-07-24 11:17:08Z abarral $

SUBROUTINE qminimum_loc( q,nqtot,deltap )
  USE parallel_lmdz
  USE infotrac, ONLY: niso, ntiso, iqIsoPha, tracers, &
        isoCheck, min_qParent
  USE strings_mod, ONLY: strIdx
  USE readTracFiles_mod, ONLY: addPhase
  IMPLICIT none
  !
  !  -- Objet : Traiter les valeurs trop petites (meme negatives)
  !         pour l'eau vapeur et l'eau liquide
  !
  include "dimensions.h"
  include "paramet.h"
  include "iniprint.h"
  !
  INTEGER :: nqtot ! CRisi: on remplace nq par nqtot
  REAL :: q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
  !
  LOGICAL, SAVE :: first=.TRUE.
  INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
!$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
  !
  !  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
  !        parametres seuil_vap, seuil_liq soient pareilles a celles
  !        qui  sont utilisees dans la routine    ADDFI       )
  ! .................................................................
  !
  !DC iq_val and iq_liq are usable for q only, NOT for q_follow
  !   and zx_defau_diag (crash if iq_val/liq==3) => vapor/liquid
  !   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)
  !
  REAL :: SSUM
  EXTERNAL SSUM
  !
  INTEGER :: imprim
  SAVE imprim
  DATA imprim /0/
!$OMP THREADPRIVATE(imprim)
  INTEGER :: ijb,ije
  INTEGER :: Index_pump(ij_end-ij_begin+1)
  INTEGER :: nb_pump
  INTEGER :: ixt
  INTEGER :: iso_verif_noNaN_nostop

!$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
  !
  ! Quand l'eau liquide est trop petite (ou negative), on prend
  ! l'eau vapeur de la meme couche et la convertit en eau liquide
  ! (sans changer la temperature !)
  !

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

  ijb=ij_begin
  ije=ij_end

  DO k = 1, llm
!$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
!$OMP END DO NOWAIT
  ENDDO

  !write(lunout,*) 'qminimum 57'
  DO k = 1, llm
!$OMP DO SCHEDULE(STATIC)
    DO i = ijb, ije
      if (seuil_liq - q(i,k,iq_liq) > 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
!$OMP END DO NOWAIT
  END DO

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

      if ( seuil_vap - q(i,k,iq_vap) > 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
!$OMP END DO NOWAIT
  ENDDO

  !
  ! Quand il s'agit de la premiere couche au-dessus du sol, on
  ! doit imprimer un message d'avertissement (saturation possible).
  !
  !write(lunout,*) 'qminimum 106'
  nb_pump=0
!$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
!$OMP END DO NOWAIT
   ! pompe = SSUM(ije-ijb+1,zx_pump(ijb),1)

  IF (imprim<=100 .AND. nb_pump > 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?
!$OMP DO SCHEDULE(STATIC)
  DO i = ijb, ije
    if (zx_pump(i)>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
!$OMP END DO NOWAIT

  ! 2) transfert de vap vers les couches plus hautes
  !write(lunout,*) 'qminimum 158'
  do k=2,llm
!$OMP DO SCHEDULE(STATIC)
    DO i = ijb, ije
      if (zx_defau_diag(i,k,1)>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)<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')==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')==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
!$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
!$OMP DO SCHEDULE(STATIC)
    DO i = ijb, ije
      if (zx_defau_diag(i,k,2)>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
!$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'
!$OMP BARRIER

  !

END SUBROUTINE qminimum_loc