source: LMDZ6/trunk/libf/dyn3dmem/vlspltgen_loc.F90 @ 5300

!
! $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 Write_Field_loc
  USE VAMPIR
  ! ! CRisi: on rajoute variables utiles d'infotrac
  USE infotrac, ONLY : nqtot, tracers, isoCheck
  USE vlspltgen_mod
  USE comconst_mod, ONLY: cpp
  USE logic_mod, ONLY: adv_qsat_liq
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
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

   ! CALL SCOPY(ijp1llm,q,1,zq,1)
   ! CALL SCOPY(ijp1llm,masse,1,zm,1)

   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)
#ifdef _ADV_HALO
  ! CRisi: on ajoute les nombres de fils et tableaux des fils
  ! On suppose qu'on ne peut advecter les fils que par le schéma 10.
      call vlx_loc(zq,pente_max,zm,mu, &
            ij_begin,ij_begin+2*iip1-1,iq)
      call vlx_loc(zq,pente_max,zm,mu, &
            ij_end-2*iip1+1,ij_end,iq)
#else
      call vlx_loc(zq,pente_max,zm,mu, &
            ij_begin,ij_end,iq)
#endif

!$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)
#ifdef _ADV_HALO
      call vlxqs_loc(zq,pente_max,zm,mu, &
            qsat,ij_begin,ij_begin+2*iip1-1,iq)
      call vlxqs_loc(zq,pente_max,zm,mu, &
            qsat,ij_end-2*iip1+1,ij_end,iq)
#else
      call vlxqs_loc(zq,pente_max,zm,mu, &
            qsat,ij_begin,ij_end,iq)
#endif

!$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)
#ifdef _ADV_HALLO
      call vlx_loc(zq,pente_max,zm,mu, &
            ij_begin+2*iip1,ij_end-2*iip1,iq)
#endif
      CASE(14)
#ifdef _ADV_HALLO
      call vlxqs_loc(zq,pente_max,zm,mu, &
            qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
#endif
      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
   ! !write(*,*) 'vlspltgen 321: iq=',iq

    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
  ! !write(*,*) 'vlspltgen 349: iq=',iq
    SELECT CASE(tracers(iq)%iadv)
      CASE(0); CYCLE
      CASE(10,14)
!$OMP BARRIER
#ifdef _ADV_HALLO
      call vlz_loc(zq,pente_max,zm,mw, &
            ij_begin,ij_begin+2*iip1-1,iq)
      call vlz_loc(zq,pente_max,zm,mw, &
            ij_end-2*iip1+1,ij_end,iq)
#else
      call vlz_loc(zq,pente_max,zm,mw, &
            ij_begin,ij_end,iq)
#endif
!$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

#ifdef _ADV_HALLO
      call vlz_loc(zq,pente_max,zm,mw, &
            ij_begin+2*iip1,ij_end-2*iip1,iq)
#endif

!$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
  ! !write(*,*) 'vlspltgen 449: iq=',iq
    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
  ! !write(*,*) 'vlspltgen 477: iq=',iq
    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
  ! !write(*,*) 'vlspltgen_loc 557'
!$OMP BARRIER

  ! !write(*,*) 'vlspltgen_loc 559'
  DO iq=1,nqtot
   ! !write(*,*) 'vlspltgen_loc 561, iq=',iq
!$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
  ! !write(*,*) 'vlspltgen_loc 575'

!$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
  ! !write(*,*) 'vlspltgen_loc 583'
  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'
  RETURN
END SUBROUTINE vlspltgen_loc