source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltgen_loc.F90 @ 5423

! $Header$

SUBROUTINE vlspltgen_loc(q, pente_max, masse, w, pbaru, pbarv, &
        pdt, p, pk, teta)


  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron

  !    ********************************************************************
  !      Schema  d'advection " pseudo amont " .
  !  + test sur humidite specifique: Q advecte< Qsat aval
  !               (F. Codron, 10/99)
  !    ********************************************************************
  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....

  ! pente_max facteur de limitation des pentes: 2 en general
  !                                            0 pour un schema amont
  ! pbaru,pbarv,w flux de masse en u ,v ,w
  ! pdt pas de temps

  ! teta temperature potentielle, p pression aux interfaces,
  ! pk exner au milieu des couches necessaire pour calculer Qsat
  !   --------------------------------------------------------------------
  USE parallel_lmdz
  USE mod_hallo
  USE lmdz_vampir
  ! CRisi: on rajoute variables utiles d'infotrac
  USE lmdz_infotrac, ONLY: nqtot, tracers, isoCheck
  USE vlspltgen_mod
  USE comconst_mod, ONLY: cpp
  USE logic_mod, ONLY: adv_qsat_liq


USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet
  IMPLICIT NONE

  !




  !   Arguments:
  !   ----------
  REAL :: masse(ijb_u:ije_u, llm), pente_max
  REAL :: pbaru(ijb_u:ije_u, llm), pbarv(ijb_v:ije_v, llm)
  REAL :: q(ijb_u:ije_u, llm, nqtot)
  REAL :: w(ijb_u:ije_u, llm), pdt
  REAL :: p(ijb_u:ije_u, llmp1), teta(ijb_u:ije_u, llm)
  REAL :: pk(ijb_u:ije_u, llm)

  !  Local
  !   ---------

  INTEGER :: ij, l

  REAL :: zzpbar, zzw

  REAL :: qmin, qmax
  DATA qmin, qmax/0., 1.e33/

  !--pour rapport de melange saturant--

  REAL :: rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
  REAL :: ptarg, pdelarg, foeew, zdelta
  REAL :: tempe(ijb_u:ije_u)
  INTEGER :: ijb, ije, iq, iq2, ifils
  LOGICAL, SAVE :: firstcall = .TRUE.
  !$OMP THREADPRIVATE(firstcall)
  type(request), SAVE :: MyRequest1
  !$OMP THREADPRIVATE(MyRequest1)
  type(request), SAVE :: MyRequest2
  !$OMP THREADPRIVATE(MyRequest2)
  !    fonction psat(T)

  FOEEW (PTARG, PDELARG) = EXP (&
          (R3LES * (1. - PDELARG) + R3IES * PDELARG) * (PTARG - RTT) &
                  / (PTARG - (R4LES * (1. - PDELARG) + R4IES * PDELARG)))

  r2es = 380.11733
  r3les = 17.269
  r3ies = 21.875
  r4les = 35.86
  r4ies = 7.66
  retv = 0.6077667
  rtt = 273.16

  ! Allocate variables depending on dynamic variable nqtot

  IF (firstcall) THEN
    firstcall = .FALSE.
  END IF
  !-- Calcul de Qsat en chaque point
  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
  !   pour eviter une exponentielle.

  CALL SetTag(MyRequest1, 100)
  CALL SetTag(MyRequest2, 101)

  ijb = ij_begin - iip1
  ije = ij_end + iip1
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud) ije = ij_end

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      tempe(ij) = teta(ij, l) * pk(ij, l) / cpp
    ENDDO
    DO ij = ijb, ije
      IF (adv_qsat_liq) THEN
        zdelta = 0.
      ELSE
        zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
      ENDIF
      play = 0.5 * (p(ij, l) + p(ij, l + 1))
      qsat(ij, l) = MIN(0.5, r2es * FOEEW(tempe(ij), zdelta) / play)
      qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  ! PRINT*,'Debut vlsplt version debug sans vlyqs'

  zzpbar = 0.5 * pdt
  zzw = pdt

  ijb = ij_begin
  ije = ij_end
  IF (pole_nord) ijb = ijb + iip1
  IF (pole_sud)  ije = ije - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      mu(ij, l) = pbaru(ij, l) * zzpbar
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ijb = ij_begin - iip1
  ije = ij_end
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud)  ije = ij_end - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      mv(ij, l) = pbarv(ij, l) * zzpbar
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ijb = ij_begin
  ije = ij_end

  DO iq = 1, nqtot
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        mw(ij, l, iq) = w(ij, l) * zzw
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT
  ENDDO

  DO iq = 1, nqtot
    !$OMP MASTER
    DO ij = ijb, ije
      mw(ij, llm + 1, iq) = 0.
    ENDDO
    !$OMP END MASTER
  ENDDO

  ijb = ij_begin
  ije = ij_end

  DO iq = 1, nqtot
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      zq(ijb:ije, l, iq) = q(ijb:ije, l, iq)
      zm(ijb:ije, l, iq) = masse(ijb:ije, l)
    ENDDO
    !$OMP END DO NOWAIT
  ENDDO

  ! verif temporaire
  ijb = ij_begin
  ije = ij_end
  CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 191')

  !$OMP BARRIER
  DO iq = 1, nqtot
    ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
    IF(tracers(iq)%parent /= 'air') CYCLE
    !WRITE(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10)
      CALL vlx_loc(zq, pente_max, zm, mu, &
              ij_begin, ij_end, iq)

      !$OMP MASTER
      CALL VTb(VTHallo)
      !$OMP END MASTER
      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
      ! CRisi
      DO ifils = 1, tracers(iq)%nqDescen
        iq2 = tracers(iq)%iqDescen(ifils)
        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
      enddo

      !$OMP MASTER
      CALL VTe(VTHallo)
      !$OMP END MASTER
    CASE(14)
      CALL vlxqs_loc(zq, pente_max, zm, mu, &
              qsat, ij_begin, ij_end, iq)

      !$OMP MASTER
      CALL VTb(VTHallo)
      !$OMP END MASTER

      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest1)
      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest1)
      DO ifils = 1, tracers(iq)%nqDescen
        iq2 = tracers(iq)%iqDescen(ifils)
        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest1)
        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest1)
      enddo

      !$OMP MASTER
      CALL VTe(VTHallo)
      !$OMP END MASTER
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_loc", "schema non parallelise", 1)
    END SELECT

  enddo !DO iq=1,nqtot


  !$OMP BARRIER
  !$OMP MASTER
  CALL VTb(VTHallo)
  !$OMP END MASTER

  CALL SendRequest(MyRequest1)

  !$OMP MASTER
  CALL VTe(VTHallo)
  !$OMP END MASTER
  !$OMP BARRIER

  ! verif temporaire
  ijb = ij_begin - 2 * iip1
  ije = ij_end + 2 * iip1
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud)  ije = ij_end
  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 280')

  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE
    !WRITE(*,*) 'vlspltgen 279: iq=',iq

    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10)
    CASE(14)
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo
  !$OMP BARRIER
  !$OMP MASTER
  CALL VTb(VTHallo)
  !$OMP END MASTER

  ! CALL WaitRecvRequest(MyRequest1)
  ! CALL WaitSendRequest(MyRequest1)
  !$OMP BARRIER
  CALL WaitRequest(MyRequest1)


  !$OMP MASTER
  CALL VTe(VTHallo)
  !$OMP END MASTER
  !$OMP BARRIER

  IF(isoCheck) THEN
    CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 326')
    ijb = ij_begin - 2 * iip1
    ije = ij_end + 2 * iip1
    IF (pole_nord) ijb = ij_begin
    IF (pole_sud)  ije = ij_end
    CALL check_isotopes(zq, ijb, ije, 'vlspltgen_loc 336')
  END IF

  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE

    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo

  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 357')

  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE
    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10, 14)
      !$OMP BARRIER
      CALL vlz_loc(zq, pente_max, zm, mw, &
              ij_begin, ij_end, iq)
      !$OMP BARRIER

      !$OMP MASTER
      CALL VTb(VTHallo)
      !$OMP END MASTER

      CALL Register_Hallo_u(zq(:, :, iq), llm, 2, 2, 2, 2, MyRequest2)
      CALL Register_Hallo_u(zm(:, :, iq), llm, 1, 1, 1, 1, MyRequest2)
      ! CRisi
      DO ifils = 1, tracers(iq)%nqDescen
        iq2 = tracers(iq)%iqDescen(ifils)
        CALL Register_Hallo_u(zq(:, :, iq2), llm, 2, 2, 2, 2, MyRequest2)
        CALL Register_Hallo_u(zm(:, :, iq2), llm, 1, 1, 1, 1, MyRequest2)
      enddo
      !$OMP MASTER
      CALL VTe(VTHallo)
      !$OMP END MASTER
      !$OMP BARRIER
    CASE DEFAULT

      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo
  !$OMP BARRIER

  !$OMP MASTER
  CALL VTb(VTHallo)
  !$OMP END MASTER

  CALL SendRequest(MyRequest2)

  !$OMP MASTER
  CALL VTe(VTHallo)
  !$OMP END MASTER

  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 429')

  !$OMP BARRIER
  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE
    !WRITE(*,*) 'vlspltgen 409: iq=',iq

    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10, 14)
      !$OMP BARRIER
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo
  !WRITE(*,*) 'vlspltgen_loc 476'

  !$OMP BARRIER
  !WRITE(*,*) 'vlspltgen_loc 477'
  !$OMP MASTER
  CALL VTb(VTHallo)
  !$OMP END MASTER

  ! CALL WaitRecvRequest(MyRequest2)
  ! CALL WaitSendRequest(MyRequest2)
  !$OMP BARRIER
  CALL WaitRequest(MyRequest2)

  !$OMP MASTER
  CALL VTe(VTHallo)
  !$OMP END MASTER
  !$OMP BARRIER


  !WRITE(*,*) 'vlspltgen_loc 494'
  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 461')

  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE
    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10); CALL   vly_loc(zq, pente_max, zm, mv, iq)
    CASE(14); CALL vlyqs_loc(zq, pente_max, zm, mv, qsat, iq)
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo !do iq=1,nqtot

  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 493')

  DO iq = 1, nqtot
    IF(tracers(iq)%parent /= 'air') CYCLE
    SELECT CASE(tracers(iq)%iadv)
    CASE(0); CYCLE
    CASE(10); CALL   vlx_loc(zq, pente_max, zm, mu, &
            ij_begin, ij_end, iq)
    CASE(14); CALL vlxqs_loc(zq, pente_max, zm, mu, &
            qsat, ij_begin, ij_end, iq)
    CASE DEFAULT
      CALL abort_gcm("vlspltgen_p", "schema non parallelise", 1)
    END SELECT

  enddo !do iq=1,nqtot

  !WRITE(*,*) 'vlspltgen 550: apres derniere serie de CALL vlx'
  CALL check_isotopes(zq, ij_begin, ij_end, 'vlspltgen_loc 521')

  ijb = ij_begin
  ije = ij_end
  !$OMP BARRIER

  DO iq = 1, nqtot
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        ! PRINT *,'zq-->',ij,l,iq,zq(ij,l,iq)
        ! PRINT *,'q-->',ij,l,iq,q(ij,l,iq)
        q(ij, l, iq) = zq(ij, l, iq)
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije - iip1 + 1, iip1
        q(ij + iim, l, iq) = q(ij, l, iq)
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT
  ENDDO !DO iq=1,nqtot

  CALL check_isotopes(q, ij_begin, ij_end, 'vlspltgen_loc 557')

  !$OMP BARRIER

  !c$OMP MASTER
  ! CALL WaitSendRequest(MyRequest1)
  ! CALL WaitSendRequest(MyRequest2)
  !c$OMP END MASTER
  !c$OMP BARRIER

  !WRITE(*,*) 'vlspltgen 597: sortie'

END SUBROUTINE vlspltgen_loc