MODULE climb_hq_mod
!
! Module to solve the vertical diffusion of "q" and "H";
! specific humidity and potential energi.
!
  USE dimphy
  USE compbl_mod_h
#ifdef ISO
  USE infotrac_phy, ONLY: ntraciso=>ntiso ! ajout C Risi pour isos      
#endif

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: climb_hq_down, climb_hq_up, climb_hq_init, climb_hq_finalize

  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
  !$OMP THREADPRIVATE(gamaq,gamah)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
  !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
  !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_Q, Bcoef_Q
  !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_H, Bcoef_H
  !$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
  !$OMP THREADPRIVATE(Kcoefhq)
  REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
  !$OMP THREADPRIVATE(h_old)
  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
  !$OMP THREADPRIVATE(d_h_col_vdf)
  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
  !$OMP THREADPRIVATE(f_h_bnd)
#ifdef ISO
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: gamaxt
  !$OMP THREADPRIVATE(gamaxt)
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: Ccoef_XT, Dcoef_XT
  !$OMP THREADPRIVATE(Ccoef_XT, Dcoef_XT)
  REAL, DIMENSION(:,:), ALLOCATABLE   :: Acoef_XT, Bcoef_XT
  !$OMP THREADPRIVATE(Acoef_XT, Bcoef_XT)
#endif

CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE climb_hq_init
  USE dimphy, ONLY : klon, klev
  IMPLICIT NONE
    INTEGER :: ierr
  
     ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
     Ccoef_Q(:,:) = 0.

     ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
     Dcoef_Q(:,:) = 0.

     ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
     Ccoef_H(:,:) = 0.

     ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
     Dcoef_H(:,:) = 0.

     ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
     Acoef_Q(:)=0. ; Bcoef_Q(:)=0. ; Acoef_H(:)=0. ; Bcoef_H(:)=0. ;

     ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
     Kcoefhq(:,:) = 0.

     ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
     IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
     gamaq(:,:) = 0.

     ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
     IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
     gamah(:,:)=0.

     ALLOCATE(h_old(klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr
     h_old(:,:) = 0.

#ifdef ISO
     ALLOCATE(Ccoef_XT(ntraciso,klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_XT, ierr=', ierr
       
     ALLOCATE(Dcoef_XT(ntraciso,klon,klev), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_XT, ierr=', ierr
       
     ALLOCATE(Acoef_XT(ntraciso,klon), Bcoef_XT(ntraciso,klon), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_XT and Bcoef_XT, ierr=', ierr
       
     ALLOCATE(gamaxt(ntraciso,1:klon,2:klev), STAT=ierr)
     IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaxt, ierr=', ierr
#endif

     ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
     d_h_col_vdf(:) = 0.

     ALLOCATE(f_h_bnd(klon), STAT=ierr)
     IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
     f_h_bnd(:) = 0.

  END SUBROUTINE climb_hq_init

  SUBROUTINE climb_hq_finalize
  IMPLICIT NONE
    INTEGER :: ierr

    !****************************************************************************************
    ! Some deallocations
    !
    !****************************************************************************************
    
    DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)    
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
    DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)    
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
    DEALLOCATE(gamaq, gamah,stat=ierr)
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
    DEALLOCATE(Kcoefhq,stat=ierr)
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
    DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr

  END  SUBROUTINE climb_hq_finalize

  SUBROUTINE climb_hq_down(knon, ni, coefhq, paprs, pplay, &
       delp, temp, q, dtime, &
!!! nrlmd le 02/05/2011
       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
       Kcoef_hq_out, gama_q_out, gama_h_out, &
!!!
       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out &
#ifdef ISO
       ,xt,  &
       Ccoef_XT_out, Dcoef_XT_out, gama_xt_out,  &
       Acoef_XT_out, Bcoef_XT_out & 
#endif               
       )
!$gpum horizontal knon

  USE yomcst_mod_h
#ifdef ISOVERIF
  USE isotopes_mod, ONLY: iso_eau,iso_HDO
 !USE isotopes_verif_mod, ONLY: errmax, errmaxrel
  USE isotopes_verif_mod
#endif

! This routine calculates recursivly the coefficients C and D
! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
! the index of the vertical layer.

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    INTEGER, INTENT(IN)                      :: ni(knon)
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: coefhq
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay 
    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs 
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: temp, delp  ! temperature
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: q
    REAL, INTENT(IN)                         :: dtime
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(IN) :: xt
#endif

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_H_out
    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_Q_out
    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_H_out
    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_Q_out
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon), INTENT(OUT) :: Acoef_XT_out
    REAL, DIMENSION(ntraciso,knon), INTENT(OUT) :: Bcoef_XT_out
#endif

!!! nrlmd le 02/05/2011
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_H_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_Q_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_H_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_Q_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Kcoef_hq_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: gama_q_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: gama_h_out
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(OUT) :: Ccoef_XT_out
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(OUT) :: Dcoef_XT_out
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(OUT) :: gama_xt_out
#endif
!!!

    REAL :: yCcoef_Q(knon,klev)
    REAL :: yDcoef_Q(knon,klev)
    REAL :: yCcoef_H(knon,klev)
    REAL :: yDcoef_H(knon,klev)
    REAL :: yAcoef_Q(knon), yBcoef_Q(knon), yAcoef_H(knon), yBcoef_H(knon)
    REAL :: yKcoefhq(knon,klev)
    REAL :: ygamaq(knon,2:klev)
    REAL :: ygamah(knon,2:klev)
    REAL :: yh_old(knon,klev) ! for diagnostics, h before solving diffusion
    REAL :: yd_h_col_vdf(knon) ! for diagnostics, vertical integral of enthalpy change
    REAL :: yf_h_bnd(knon) ! for diagnostics, enthalpy flux at surface
#ifdef ISO
    REAL :: yCcoef_XT(ntraciso,knon,klev)
    REAL :: yDcoef_XT(ntraciso,knon,klev)
    REAL :: yAcoef_XT(ntraciso,knon), yBcoef_XT(ntraciso,knon)
    REAL :: ygamaxt(ntraciso,knon,2:klev)
#endif
!!!

! Local variables
!****************************************************************************************
! JLD now renamed h_old and declared in module
!    REAL, DIMENSION(knon,klev)               :: local_H
    REAL, DIMENSION(knon)                    :: psref 
    REAL                                     :: delz, pkh
    INTEGER                                  :: k, i, j, ierr

#ifdef ISO
    REAL, DIMENSION(knon,2:klev) :: ygamaxt_tmp
    REAL, DIMENSION(knon,klev)   :: xt_tmp
    REAL, DIMENSION(knon,klev)   :: yCcoef_XT_tmp,yDcoef_XT_tmp
    REAL, DIMENSION(knon)        :: yAcoef_XT_tmp,yBcoef_XT_tmp
    INTEGER :: ixt
#endif

    yd_h_col_vdf(:) = 0.
    yf_h_bnd(:) = 0.

#ifdef ISO
#ifdef ISOVERIF
   IF (iso_eau.GT.0) THEN
     DO k = 1, klev
       DO i = 1, knon
         CALL iso_verif_egalite_choix( &
              xt(iso_eau,i,k),q(i,k), &
              'climb_hq 100',errmax,errmaxrel)
       ENDDO
     ENDDO
   ENDIF ! if (iso_eau.gt.0) then
#endif
#endif

!****************************************************************************************
! 2)
! Definition of the coeficient K 
!
!****************************************************************************************
    yKcoefhq(:,:) = 0.0
    DO k = 2, klev
      DO i = 1, knon
          yKcoefhq(i,k) = &
             coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
             *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
      ENDDO
    ENDDO

!****************************************************************************************
! 3)
! Calculation of gama for "Q" and "H"
!
!****************************************************************************************
!   surface pressure is used as reference
    psref(:) = paprs(:,1) 

!   definition of gama
    IF (iflag_pbl == 1) THEN
      ygamaq(:,:) = 0.0
      ygamah(:,:) = -1.0e-03
      ygamah(:,2) = -2.5e-03
#ifdef ISO
      DO ixt=1,ntraciso
        ygamaxt(:,:,:) = 0.0
      ENDDO ! do ixt=1,ntraciso
#endif
 
! conversion de gama
      DO k = 2, klev
        DO i = 1, knon
          delz = RD * (temp(i,k-1)+temp(i,k)) / & 
                 2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
          pkh  = (psref(i)/paprs(i,k))**RKAPPA
          
! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
          ygamaq(i,k) = ygamaq(i,k) * delz    
! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches 
          ygamah(i,k) = ygamah(i,k) * delz * RCPD * pkh
#ifdef ISO
          DO ixt=1,ntraciso
            ygamaxt(ixt,i,k) = ygamaxt(ixt,i,k) * delz
          ENDDO
#endif
        ENDDO
      ENDDO

    ELSE
      ygamaq(:,:) = 0.0
      ygamah(:,:) = 0.0
#ifdef ISO
      DO ixt=1,ntraciso
        ygamaxt(:,:,:) = 0.0
      ENDDO ! do ixt=1,ntraciso
#endif
    ENDIF
    
#ifdef ISO
#ifdef ISOVERIF
    DO k = 2, klev
      DO i = 1, knon
        CALL iso_verif_egalite_choix( &
             ygamaxt(iso_eau,i,k),ygamaq(i,k), &
             'climb_hq 209',errmax,errmaxrel)
      ENDDO
    ENDDO
#endif
#endif
    

!****************************************************************************************    
! 4)
! Calculate the coefficients C and D for specific humidity, q
!
!****************************************************************************************
    
    CALL calc_coef(knon, yKcoefhq, ygamaq, delp, q, &
         yCcoef_Q, yDcoef_Q, yAcoef_Q, yBcoef_Q)


#ifdef ISO
    DO ixt=1,ntraciso        
      ! compression
      DO k = 2, klev
        DO i = 1, knon
          ygamaxt_tmp(i,k)=ygamaxt(ixt,i,k)
        ENDDO
      ENDDO !do k = 2, klev
      DO k = 1, klev
        DO i = 1, knon
          xt_tmp(i,k)=xt(ixt,i,k)
        ENDDO
      ENDDO !do k = 2, klev
      !appel routine generique
      CALL calc_coef(knon, yKcoefhq, ygamaxt_tmp, delp, xt_tmp, &
           yCcoef_XT_tmp, yDcoef_XT_tmp, yAcoef_XT_tmp, yBcoef_XT_tmp)
      ! decompression
      DO k = 1, klev
        DO i = 1, knon
          yCcoef_XT(ixt,i,k)=yCcoef_XT_tmp(i,k)
          yDcoef_XT(ixt,i,k)=yDcoef_XT_tmp(i,k)
        ENDDO
      ENDDO !do k = 2, klev
      DO i = 1, knon
        yAcoef_XT(ixt,i)=yAcoef_XT_tmp(i)
        yBcoef_XT(ixt,i)=yBcoef_XT_tmp(i)
      ENDDO
    ENDDO ! do ixt=1,ntraciso
#ifdef ISOVERIF
    IF (iso_eau.GT.0) THEN
      DO k = 1, klev
        DO i = 1, knon
          CALL iso_verif_egalite_choix( &
              yCcoef_XT(iso_eau,i,k),yCcoef_Q(i,k), &
                      'climb_hq 234c',errmax,errmaxrel)
          CALL iso_verif_egalite_choix( &
              yDcoef_XT(iso_eau,i,k),yDcoef_Q(i,k), &
                      'climb_hq 234d',errmax,errmaxrel)
        ENDDO !DO i = 1, knon
      ENDDO !do k = 2, klev
      DO i = 1, knon
        CALL iso_verif_egalite_choix( &
            yAcoef_XT(iso_eau,i),yAcoef_Q(i), &
                    'climb_hq 234a',errmax,errmaxrel)
        CALL iso_verif_egalite_choix( &
            yBcoef_XT(iso_eau,i),yBcoef_Q(i), &
                    'climb_hq 234b',errmax,errmaxrel)
      ENDDO !DO i = 1, knon
    ENDIF !if (iso_eau.gt.0) then
#endif
#endif  
!****************************************************************************************
! 5)
! Calculate the coefficients C and D for potentiel enthalpy, H 
!
!****************************************************************************************
    yh_old(:,:) = 0.0

    DO k=1,klev
      DO i = 1, knon
        ! convertie la temperature en entalpie potentielle
          yh_old(i,k) = RCPD * temp(i,k) * &
             (psref(i)/pplay(i,k))**RKAPPA
      ENDDO
    ENDDO

    CALL calc_coef(knon, yKcoefhq, ygamah, delp, yh_old, &
         yCcoef_H, yDcoef_H, yAcoef_H, yBcoef_H)
 
!****************************************************************************************
! 6)
! Return the first layer in output variables
!
!****************************************************************************************
    Acoef_H_out = yAcoef_H
    Bcoef_H_out = yBcoef_H
    Acoef_Q_out = yAcoef_Q
    Bcoef_Q_out = yBcoef_Q
#ifdef ISO
    Acoef_XT_out = yAcoef_XT
    Bcoef_XT_out = yBcoef_XT
#endif

!****************************************************************************************
! 7)
! If Pbl is split, return also the other layers in output variables
!
!****************************************************************************************
!!! jyg le 07/02/2012
!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
  IF (mod(iflag_pbl_split,10) .ge.1) THEN
!!! nrlmd le 02/05/2011
    DO k= 1, klev
      DO i= 1, knon
        Ccoef_H_out(i,k) = yCcoef_H(i,k)
        Dcoef_H_out(i,k) = yDcoef_H(i,k)
        Ccoef_Q_out(i,k) = yCcoef_Q(i,k)
        Dcoef_Q_out(i,k) = yDcoef_Q(i,k)
        Kcoef_hq_out(i,k) = yKcoefhq(i,k)
#ifdef ISO
        DO ixt=1,ntraciso
          Ccoef_XT_out(ixt,i,k) = yCcoef_XT(ixt,i,k)
          Dcoef_XT_out(ixt,i,k) = yDcoef_XT(ixt,i,k)              
        ENDDO   
#endif
        IF (k.eq.1) THEN
          gama_h_out(i,k)  = 0.
          gama_q_out(i,k)  = 0.
#ifdef ISO
          DO ixt=1,ntraciso
            gama_xt_out(ixt,i,k)  = 0. 
          ENDDO   
#endif
        ELSE
          gama_h_out(i,k)  = ygamah(i,k)
          gama_q_out(i,k)  = ygamaq(i,k)
#ifdef ISO
          DO ixt=1,ntraciso
              gama_xt_out(ixt,i,k)  = ygamaxt(ixt,i,k) 
          ENDDO   
#endif
        ENDIF
      ENDDO
    ENDDO
!!!      
  ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!
    DO k=1,klev
      DO j=1,knon
        i=ni(j)
        IF (k==1) THEN
          Acoef_Q(i) = yAcoef_Q(j) 
          Bcoef_Q(i) = yBcoef_Q(j) 
          Acoef_H(i) = yAcoef_H(j) 
          Bcoef_H(i) = yBcoef_H(j)
          d_h_col_vdf(i)= yd_h_col_vdf(j)
          f_h_bnd(i)= yf_h_bnd(j)
#ifdef ISO
          DO ixt=1,ntraciso
            Acoef_XT(ixt,i) = yAcoef_XT(ixt,j)
            Bcoef_XT(ixt,i) = yBcoef_XT(ixt,j)
          ENDDO      
#endif
        ENDIF
        
        IF (k>=2) THEN
          gamaq(i,k)=ygamaq(j,k)
          gamah(i,k)=ygamah(j,k)
#ifdef ISO
          DO ixt=1,ntraciso
            gamaxt(ixt,i,k) = ygamaxt(ixt,j,k)
          ENDDO      
#endif
        ENDIF

        Ccoef_Q(i,k) = yCcoef_Q(j,k)
        Dcoef_Q(i,k) = yDcoef_Q(j,k)
        Ccoef_H(i,k) = yCcoef_H(j,k)
        Dcoef_H(i,k) = yDcoef_H(j,k)
        Kcoefhq(i,k) = yKcoefhq(j,k)
        h_old(i,k) = yh_old(j,k)
#ifdef ISO
        DO ixt=1,ntraciso
          Ccoef_XT(ixt,i,k) = yCcoef_XT(ixt,j,k)
          Dcoef_XT(ixt,i,k) = yDcoef_XT(ixt,j,k)
        ENDDO      
#endif
      ENDDO
    ENDDO

  END SUBROUTINE climb_hq_down
!
!****************************************************************************************
!
  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
!$gpum horizontal knon
!
! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
! where X is H or Q, and k the vertical level k=1,klev
!
  USE yomcst_mod_h

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef, delp
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: X
    REAL, DIMENSION(knon,2:klev), INTENT(IN) :: gama

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef, Bcoef
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef, Dcoef

! Local variables
!****************************************************************************************
    INTEGER                                  :: k, i
    REAL                                     :: buf

!****************************************************************************************
! Niveau au sommet, k=klev
!
!****************************************************************************************
    Ccoef(:,:) = 0.0
    Dcoef(:,:) = 0.0

    DO i = 1, knon
      buf = delp(i,klev) + Kcoef(i,klev)
       
      Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
      Dcoef(i,klev) = Kcoef(i,klev)/buf
    ENDDO


!****************************************************************************************
! Niveau  (klev-1) <= k <= 2
!
!****************************************************************************************

    DO k=(klev-1),2,-1
      DO i = 1, knon
        buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
        Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
        Dcoef(i,k) = Kcoef(i,k)/buf
      ENDDO
    ENDDO

!****************************************************************************************
! Niveau k=1
!
!****************************************************************************************

    DO i = 1, knon
      buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
      Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
      Bcoef(i) = -1. * RG / buf
    ENDDO

  END SUBROUTINE calc_coef
!
!****************************************************************************************
!
  SUBROUTINE climb_hq_up(knon, ni, dtime, t_old, q_old, &
       flx_q1, flx_h1, paprs, pplay, &
!!! nrlmd le 02/05/2011
       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
       Kcoef_hq_in, gama_q_in, gama_h_in, &
!!!
       flux_q, flux_h, d_q, d_t &
#ifdef ISO
       ,xt_old, flx_xt1, &
       Acoef_XT_in,Bcoef_XT_in,Ccoef_XT_in,Dcoef_XT_in,gama_xt_in,  &
       flux_xt, d_xt &
#endif       
       )
!$gpum horizontal knon

! 
! This routine calculates the flux and tendency of the specific humidity q and 
! the potential engergi H. 
! The quantities q and H are calculated according to 
! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients 
! C and D are known from before and k is index of the vertical layer.
!   
  USE yomcst_mod_h
  USE compbl_mod_h
#ifdef ISOVERIF
  USE infotrac_phy, ONLY: nzone
  USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18, ridicule
  USE isotopes_verif_mod
#endif

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    INTEGER, INTENT(IN)                      :: ni(knon)
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: t_old, q_old
    REAL, DIMENSION(knon), INTENT(IN)        :: flx_q1, flx_h1
    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(IN) :: xt_old
    REAL, DIMENSION(ntraciso,knon), INTENT(IN) :: flx_xt1
#endif

!!! nrlmd le 02/05/2011
    REAL, DIMENSION(knon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon), INTENT(IN)      :: Acoef_XT_in,  Bcoef_XT_in
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(IN) :: Ccoef_XT_in,  Dcoef_XT_in
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(IN) :: gama_xt_in
#endif
!!!

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon,klev), INTENT(OUT) :: flux_xt, d_xt
#endif

! Local variables
!****************************************************************************************
    REAL :: yCcoef_Q(knon,klev)
    REAL :: yDcoef_Q(knon,klev)
    REAL :: yCcoef_H(knon,klev)
    REAL :: yDcoef_H(knon,klev)
    REAL :: yAcoef_Q(knon), yBcoef_Q(knon), yAcoef_H(knon), yBcoef_H(knon)
    REAL :: yKcoefhq(knon,klev)
    REAL :: ygamaq(knon,2:klev)
    REAL :: ygamah(knon,2:klev)
    REAL :: yh_old(knon,klev)
    REAL :: yd_h_col_vdf(knon)
    REAL :: yf_h_bnd(knon)
#ifdef ISO
    REAL :: yCcoef_XT(ntraciso,knon,klev)
    REAL :: yDcoef_XT(ntraciso,knon,klev)
    REAL :: yAcoef_XT(ntraciso,knon), yBcoef_XT(ntraciso,knon)
    REAL :: ygamaxt(ntraciso,knon,2:klev)
#endif

    REAL, DIMENSION(knon,klev)               :: h_new, q_new
    REAL, DIMENSION(knon)                    :: psref         
    INTEGER                                  :: k, i, j, ierr

#ifdef ISO    
    REAL, DIMENSION(ntraciso,knon,klev)      :: xt_new
    INTEGER :: ixt
!#ifdef ISOVERIF
!    integer iso_verif_noNaN_nostop
!#endif
#endif
 
!****************************************************************************************
! 1) 
! Definition of some variables
    REAL, DIMENSION(knon,klev)               :: d_h, zairm
!
!****************************************************************************************

#ifdef ISO
#ifdef ISOVERIF
    DO k = 1, klev
      DO i = 1, knon  
        IF (iso_eau.GT.0) THEN
          CALL iso_verif_egalite(xt_old(iso_eau,i,k), &
    &          q_old(i,k),'climb_hq_mod 421')
        ENDIF
#ifdef ISOTRAC
        IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422')
#endif
      ENDDO
    ENDDO
#endif  
#endif

    flux_q(:,:) = 0.0
    flux_h(:,:) = 0.0
    d_q(:,:)    = 0.0
    d_t(:,:)    = 0.0
    d_h(:,:)    = 0.0
!    f_h_bnd(:)= 0.0 !SN plus utile on utilise yf_h_bnd

    psref(1:knon) = paprs(1:knon,1)  

#ifdef ISO
    flux_xt(:,:,:) = 0.0
    d_xt(:,:,:)    = 0.0    
    xt_new(:,:,:)  = 0.0
#endif

    DO k=1,klev
      DO j=1,knon
        i=ni(j)
        IF (k==1) THEN
          yAcoef_Q(j) = Acoef_Q(i) 
          yBcoef_Q(j) = Bcoef_Q(i) 
          yAcoef_H(j) = Acoef_H(i) 
          yBcoef_H(j) = Bcoef_H(i)
          yd_h_col_vdf(j)= d_h_col_vdf(i)
          yf_h_bnd(j)= f_h_bnd(i)
#ifdef ISO
          DO ixt=1, ntraciso
            yAcoef_XT(ixt,j) = Acoef_XT(ixt,i) 
            yBcoef_XT(ixt,j) = Bcoef_XT(ixt,i)
          ENDDO
#endif
        ENDIF
        
        IF (k>=2) THEN
          ygamaq(j,k)=gamaq(i,k)
          ygamah(j,k)=gamah(i,k)
#ifdef ISO
          DO ixt=1, ntraciso
            ygamaxt(ixt,j,k)=gamaxt(ixt,i,k)
          ENDDO
#endif
        ENDIF
        yCcoef_Q(j,k) = Ccoef_Q(i,k)
        yDcoef_Q(j,k) = Dcoef_Q(i,k)
        yCcoef_H(j,k) = Ccoef_H(i,k)
        yDcoef_H(j,k) = Dcoef_H(i,k)
        yKcoefhq(j,k) = Kcoefhq(i,k)
        yh_old(j,k) = h_old(i,k)
#ifdef ISO
          DO ixt=1, ntraciso
            yCcoef_XT(ixt,j,k) = Ccoef_XT(ixt,i,k)
            yDcoef_XT(ixt,j,k) = Dcoef_XT(ixt,i,k)
          ENDDO
#endif
      ENDDO
    ENDDO      
    yf_h_bnd(:)= 0.0
!!! jyg le 07/02/2012
!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
    IF (mod(iflag_pbl_split,10).GE.1) THEN
!!! nrlmd le 02/05/2011
      DO i = 1, knon
        yAcoef_H(i)=Acoef_H_in(i)
        yAcoef_Q(i)=Acoef_Q_in(i)
        yBcoef_H(i)=Bcoef_H_in(i)
        yBcoef_Q(i)=Bcoef_Q_in(i)
#ifdef ISO
        DO ixt=1, ntraciso
          yAcoef_XT(ixt,i)=Acoef_XT_in(ixt,i)
          yBcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i)
        ENDDO    
#endif
      ENDDO
      DO k = 1, klev
        DO i = 1, knon
          yCcoef_H(i,k)=Ccoef_H_in(i,k)
          yCcoef_Q(i,k)=Ccoef_Q_in(i,k)
          yDcoef_H(i,k)=Dcoef_H_in(i,k)
          yDcoef_Q(i,k)=Dcoef_Q_in(i,k)
          yKcoefhq(i,k)=Kcoef_hq_in(i,k)
#ifdef ISO
          DO ixt=1,ntraciso
            yCcoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k)
            yDcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k)
          ENDDO    
#endif
          IF (k.GT.1) THEN
            ygamah(i,k)=gama_h_in(i,k)
            ygamaq(i,k)=gama_q_in(i,k)
#ifdef ISO
            DO ixt=1,ntraciso
              ygamaxt(ixt,i,k)=gama_xt_in(ixt,i,k)
            ENDDO 
#endif
          ENDIF
        ENDDO
      ENDDO
    
!!!      
    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!

!****************************************************************************************
! 2)
! Calculation of Q and H
!
!****************************************************************************************

!- First layer
    q_new(1:knon,1) = yAcoef_Q(1:knon) + yBcoef_Q(1:knon)*flx_q1(1:knon)*dtime
    h_new(1:knon,1) = yAcoef_H(1:knon) + yBcoef_H(1:knon)*flx_h1(1:knon)*dtime
    yf_h_bnd(1:knon) = flx_h1(1:knon)
#ifdef ISO
    DO ixt=1,ntraciso
      xt_new(ixt,1:knon,1) = yAcoef_XT(ixt,1:knon) + yBcoef_XT(ixt,1:knon)*flx_xt1(ixt,1:knon)*dtime
    ENDDO ! do ixt=1,ntraciso
#endif
!- All the other layers 
    DO k = 2, klev
      DO i = 1, knon
        q_new(i,k) = yCcoef_Q(i,k) + yDcoef_Q(i,k)*q_new(i,k-1)
        h_new(i,k) = yCcoef_H(i,k) + yDcoef_H(i,k)*h_new(i,k-1)
#ifdef ISO
        DO ixt=1,ntraciso
          xt_new(ixt,i,k) = yCcoef_XT(ixt,i,k) + yDcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1)
        ENDDO ! do ixt=1,ntraciso
#endif
      ENDDO
    ENDDO

#ifdef ISO
#ifdef ISOVERIF
    DO k = 1, klev
      DO i = 1, knon     
        DO ixt=1,ntraciso
          IF (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').EQ.1) THEN
            WRITE(*,*) 'Acoef_XT(ixt,i)=',yAcoef_XT(ixt,i)
            WRITE(*,*) 'Bcoef_XT(ixt,i)=',yBcoef_XT(ixt,i)
            WRITE(*,*) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i)
            IF (k.GE.2) THEN
              WRITE(*,*) 'Ccoef_XT(ixt,i,k)=',yCcoef_XT(ixt,i,k)
              WRITE(*,*) 'Dcoef_XT(ixt,i,k)=',yDcoef_XT(ixt,i,k)
            ENDIF
            STOP
          ENDIF
        ENDDO !do ixt=1,ntraciso
      ENDDO
    ENDDO        
#endif
#ifdef ISOVERIF
    IF (iso_eau.GT.0) THEN
      CALL iso_verif_egalite_vect2D( &
     &           xt_new,q_new, &
     &           'climb_hq_mod 504',ntraciso,knon,klev)
    ENDIF !if (iso_eau.gt.0) then
    IF ((iso_HDO.GT.0).AND.(iso_O18.GT.0)) THEN  
      DO k=1,klev
        DO i=1,knon 
          IF (q_new(i,k).GT.ridicule) THEN  
            IF (iso_verif_o18_aberrant_nostop( &
     &        xt_new(iso_HDO,i,k)/q_new(i,k), &
     &        xt_new(iso_O18,i,k)/q_new(i,k), &
     &               'climb_hq_mod 690').EQ.1) THEN
              WRITE(*,*) 'i,k,q_new(i,k)=',i,k,q_new(i,k)                       
              STOP
            ENDIF !  if (iso_verif_o18_aberrant_nostop 
          ENDIF !if (q_seri(i,k).gt.errmax) then   
        ENDDO !k=1,klev
      ENDDO  !i=1,klon
    ENDIF !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then   
#endif
#endif   
!****************************************************************************************
! 3)
! Calculation of the flux for Q and H
!
!****************************************************************************************

!- The flux at first layer, k=1
    flux_q(1:knon,1)=flx_q1(1:knon)
    flux_h(1:knon,1)=flx_h1(1:knon)
#ifdef ISO
    DO ixt=1,ntraciso
      flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon)
    ENDDO ! do ixt=1,ntraciso
#endif

!- The flux at all layers above surface
    DO k = 2, klev
      DO i = 1, knon
        flux_q(i,k) = (yKcoefhq(i,k)/RG/dtime) * &
             (q_new(i,k)-q_new(i,k-1)+ygamaq(i,k))

        flux_h(i,k) = (yKcoefhq(i,k)/RG/dtime) * &
             (h_new(i,k)-h_new(i,k-1)+ygamah(i,k)) 
#ifdef ISO
        DO ixt=1,ntraciso
          flux_xt(ixt,i,k) = (yKcoefhq(i,k)/RG/dtime) * &
          (xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+ygamaxt(ixt,i,k))
        ENDDO ! do ixt=1,ntraciso
#endif
      ENDDO
    ENDDO

!****************************************************************************************
! 4)
! Calculation of tendency for Q and H
!
!****************************************************************************************
    yd_h_col_vdf(:) = 0.0
    DO k = 1, klev
      DO i = 1, knon
        d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
        d_q(i,k) = q_new(i,k) - q_old(i,k)
        d_h(i,k) = h_new(i,k) - yh_old(i,k)
!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
    ! layer air mass
        zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
        yd_h_col_vdf(i) = yd_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
#ifdef ISO
        DO ixt=1,ntraciso
          d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k)
        ENDDO ! do ixt=1,ntraciso
#ifdef ISOVERIF
        DO ixt=1,ntraciso
          CALL iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562')
          CALL iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562')
        ENDDO ! do ixt=1,ntraciso
#endif
#ifdef ISOVERIF
        IF (iso_eau.GT.0) THEN
          CALL iso_verif_egalite(d_xt(iso_eau,i,k), &
     &         d_q(i,k),'climb_hq_mod 503') 
          CALL iso_verif_egalite(xt_new(iso_eau,i,k), &
     &         q_new(i,k),'climb_hq_mod 503b')
        ENDIF
#ifdef ISOTRAC
        IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526')
#endif       
#endif        
#endif
      ENDDO
    ENDDO

#ifdef ISO
#ifdef ISOVERIF
!    write(*,*) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1)
    IF (iso_eau.GT.0) THEN
      CALL iso_verif_egalite_vect2D( &
                           d_xt,d_q, &
           'climb_hq_mod 761',ntraciso,knon,klev)
    ENDIF       
#endif
#endif

    DO k=1,klev
      DO j=1,knon
        i=ni(j)
        IF (k==1) THEN
          Acoef_Q(i) = yAcoef_Q(j) 
          Bcoef_Q(i) = yBcoef_Q(j) 
          Acoef_H(i) = yAcoef_H(j) 
          Bcoef_H(i) = yBcoef_H(j)
          d_h_col_vdf(i)= yd_h_col_vdf(j)
          f_h_bnd(i)= yf_h_bnd(j)
#ifdef ISO
          DO ixt=1, ntraciso
            Acoef_XT(ixt,i) = yAcoef_XT(ixt,j)
            Bcoef_XT(ixt,i) = yBcoef_XT(ixt,j)
          ENDDO
#endif
        ENDIF
        
        IF (k>=2) THEN
          gamaq(i,k)=ygamaq(j,k)
          gamah(i,k)=ygamah(j,k)
#ifdef ISO
          DO ixt=1, ntraciso
            gamaxt(ixt,i,k)=ygamaxt(ixt,j,k)
          ENDDO
#endif
        ENDIF

        Ccoef_Q(i,k) = yCcoef_Q(j,k)
        Dcoef_Q(i,k) = yDcoef_Q(j,k)
        Ccoef_H(i,k) = yCcoef_H(j,k)
        Dcoef_H(i,k) = yDcoef_H(j,k)
        Kcoefhq(i,k) = yKcoefhq(j,k)
        h_old(i,k) = yh_old(j,k)
#ifdef ISO
          DO ixt=1, ntraciso
            Ccoef_XT(ixt,i,k) = yCcoef_XT(ixt,j,k)
            Dcoef_XT(ixt,i,k) = yDcoef_XT(ixt,j,k)
          ENDDO
#endif
      ENDDO
    ENDDO

  END SUBROUTINE climb_hq_up
!
!****************************************************************************************
!
END MODULE climb_hq_mod