source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/qminimum_loc.f90 @ 5153

! $Id: qminimum_loc.f90 5134 2024-07-26 15:56:37Z abarral $

SUBROUTINE qminimum_loc(q, nqtot, deltap)
  USE parallel_lmdz
  USE infotrac, ONLY: niso, ntiso, iqIsoPha, tracers, &
          isoCheck, min_qParent
  USE lmdz_strings, ONLY: strIdx
  USE lmdz_readTracFiles, ONLY: addPhase
  USE lmdz_iniprint, ONLY: lunout, prt_level


  IMPLICIT NONE
  !
  !  -- Objet : Traiter les valeurs trop petites (meme negatives)
  !         pour l'eau vapeur et l'eau liquide
  !
  INCLUDE "dimensions.h"
  INCLUDE "paramet.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)
  !
  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