MODULE lmdz_wake_dadv
PUBLIC wake_dadv
CONTAINS

  SUBROUTINE wake_dadv(klon, klev, dtime, ph, ppi, wk_adv, kupper,  &
                         deltomg, dp_deltomg, sigmaw, dsigspread,  &
                         thw, thx, qw, qx, &
                         d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv)

  USE lmdz_wake_ini , ONLY : flag_dadv_implicit

IMPLICIT NONE

  INTEGER, INTENT(IN)                                     :: klon, klev
  REAL,                               INTENT(IN)          :: dtime
  REAL, DIMENSION (klon, klev+1),     INTENT(IN)          :: ph
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: ppi
  LOGICAL, DIMENSION (klon),          INTENT(IN)          :: wk_adv
  INTEGER, DIMENSION (klon),          INTENT(IN)          :: kupper
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: deltomg
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: dp_deltomg
  REAL, DIMENSION (klon),             INTENT(IN)          :: sigmaw
  REAL, DIMENSION (klon),             INTENT(IN)          :: dsigspread
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: thw           ! component # 1
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: thx           ! component # 2
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: qw            ! component # 1
  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: qx            ! component # 2

  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_deltat_dadv
  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_deltaq_dadv
  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_tb_dadv
  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_qb_dadv

! Internal variables
  INTEGER                               :: i, k
  REAL, DIMENSION (klon, klev)          :: entr_s    ! entrainment into wakes due to spread
  REAL, DIMENSION (klon, klev)          :: thb, qb
  REAL, DIMENSION (klon, klev)          :: delta_th, delta_q

! Tests

! Arrays used in the implicit scheme 
  REAL, DIMENSION (klon)                :: rr11, rr12, rr21, rr22 

  REAL, DIMENSION (klon, klev)          :: aa11, aa12, aa21, aa22 
  REAL, DIMENSION (klon, klev)          :: bb11, bb12, bb21, bb22 
  REAL, DIMENSION (klon, klev)          :: cc11, cc12, cc21, cc22 

  REAL, DIMENSION (klon, klev)          :: alpha11, alpha12, alpha21, alpha22 
  REAL, DIMENSION (klon, klev)          :: beta11, beta12, beta21, beta22 
  REAL, DIMENSION (klon, klev)          :: gamma11, gamma12, gamma21, gamma22 
  REAL, DIMENSION (klon, klev)          :: ai11, ai12, ai21, ai22             ! inverse of alpha

  REAL, DIMENSION (klon, klev)          :: xt1, xt2, xq1, xq2
  REAL, DIMENSION (klon, klev)          :: yt1, yt2, yq1, yq2
  REAL, DIMENSION (klon, klev)          :: zt1, zt2, zq1, zq2
  REAL, DIMENSION (klon, klev)          :: th1, th2, q1, q2

  REAL                                  :: coef, det

  REAL, DIMENSION (klon,klev)           :: xt1inv, xt2inv, xq1inv, xq2inv

! Arrays used in the explicit scheme (vertical gradients)
  REAL, DIMENSION (klon, klev)          :: d_thx, d_qx
  REAL, DIMENSION (klon, klev)          :: d_thw, d_qw
  REAL, DIMENSION (klon, klev)          :: d_dth, d_dq

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! print *,'ZZZwake_dadv_IN wk_adv(1) ', wk_adv(1)
! print *,'ZZZwake_dadv_IN kupper(1) ', kupper(1)
! print *,'ZZZwake_dadv_IN k, thw(1,k), thx(1,k) ', (k, thw(1,k), thx(1,k), k = 1,3)
! print *,'ZZZwake_dadv_IN k, deltomg(1,k) ', (k, deltomg(1,k), k = 1,3)
! print *,'ZZZwake_dadv_IN k, dp_deltomg(1,k) ', (k, dp_deltomg(1,k), k = 1,3)
! print *,'ZZZwake_dadv_IN sigmaw(1) ', sigmaw(1)
! print *,'ZZZwake_dadv_IN dsigspread(1) ', dsigspread(1)

    entr_s(:,:) = 0.
    delta_th(:,:) = 0.
   
    d_deltat_dadv(:,:) = 0.
    d_deltaq_dadv(:,:) = 0.
    d_tb_dadv(:,:) = 0.
    d_qb_dadv(:,:) = 0.


    rr11(:) = sigmaw(:)
    rr12(:) = 1.-sigmaw(:)
    rr21(:) = 1.
    rr22(:) = -1.

    DO k = 1, klev
      DO i = 1,klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
         thb(i,k)      = rr11(i)*thw(i,k)+rr12(i)*thx(i,k)
         delta_th(i,k) = rr21(i)*thw(i,k)+rr22(i)*thx(i,k)
       
         qb(i,k)      = rr11(i)*qw(i,k)+rr12(i)*qx(i,k)
         delta_q(i,k) = rr21(i)*qw(i,k)+rr22(i)*qx(i,k)
        ENDIF
      ENDDO
    ENDDO

    DO i = 1, klon
        entr_s(i,klev) = 0.
    ENDDO

    DO k = 1, klev-1
      DO i = 1,klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
!!        entr_s(i,k) = dsigspread(i) - sigmaw(i)*(1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) / &
!!                     (ph(i,k)-ph(i,k+1))   

          entr_s(i,k) = dsigspread(i) + sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k)
!!  print *,'dadv, k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1)) ', &
!!                 k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1))

        ENDIF
      ENDDO
    ENDDO

! -------------------------------------------------------------------------------------------
!   Depending on flag_dadv_implicit, use implicit upstream scheme or explicit upstream scheme
! -------------------------------------------------------------------------------------------

  IF (flag_dadv_implicit) THEN

!   Implicit scheme : solve for d_deltat_dadv and d_tb_dadv 
!                     (and similarly for d_deltaq_dadv and d_qb_dadv).
!                     The system to be inverted is block-tridiagonal with 2x2 blocks.
! -----------------------------------------------------------------------------------------

!        Matrix indexing:                   Theta_w     Theta_x
!
!                                           /               
!                           Theta_b        |  A11        A12 |
!                                          |                 |
!                           delta_theta    |  A21        A22 |
!                                                           /
!       Tridiagonal matrix
!         /                                                          
!         |   aa(1)   bb(1)  0                                       |
!         |   cc(2)   aa(2)  bb(2)  0                                |
!         |   0       cc(3)  aa(3)  bb(3)                            |
!         |                                                          |
!                     .      .      .       .                         
!                                                                     
!                            .      .       .       .                 
!         |                                                          |
!         |                         cc(n-2) aa(n-2) bb(n-2) 0        |    
!         |                         0       cc(n-1) aa(n-1) bb(n-1)  |             
!                                           0       cc(n)   aa(n)    /             
! -----------------------------------------------------------------------------------------

!! Building the tridiagonal matrix
    DO i = 1,klon
      IF (wk_adv(i)) THEN 
        k = kupper(i)
        coef = dtime/(ph(i,k)-ph(i,k+1))
        aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
        aa12(i,k) = rr12(i)
        aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
                                     (1.-sigmaw(i))*deltomg(i,k) )
        aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) - &
                                     deltomg(i,k) )
   
        cc11(i,k) = 0.
        cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
        cc21(i,k) = 0.
        cc22(i,k) = coef*sigmaw(i)*deltomg(i,k)
      ENDIF  ! (wk_adv(i))
    ENDDO
    DO k = 2, klev-1
      DO i = 1,klon
        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
          coef = dtime/(ph(i,k)-ph(i,k+1))
          aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
          aa12(i,k) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1)
          aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
                                     (1.-sigmaw(i))*deltomg(i,k) )
          aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
                                     (1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) - &
                                     sigmaw(i)*deltomg(i,k) )
   
          bb11(i,k) =  -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1)
          bb12(i,k) =  0.
          bb21(i,k) =  -coef*(1.-sigmaw(i))*deltomg(i,k+1)
          bb22(i,k) =  0.
   
          cc11(i,k) = 0.
          cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
          cc21(i,k) = 0.
          cc22(i,k) = coef*sigmaw(i)*deltomg(i,k)
        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
      ENDDO
    ENDDO
    DO i = 1,klon
      IF (wk_adv(i)) THEN 
        coef = dtime/(ph(i,1)-ph(i,2))
        aa11(i,1) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,1)
        aa12(i,1) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2)
        aa21(i,1) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + &
                                   (1.-sigmaw(i))*deltomg(i,1) )
        aa22(i,1) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + &
                                   (1.-sigmaw(i))*(deltomg(i,2)-deltomg(i,1)) - &
                                   sigmaw(i)*deltomg(i,1) )
   
        bb11(i,1) =  -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2)
        bb12(i,1) =  0.
        bb21(i,1) =  -coef*(1.-sigmaw(i))*deltomg(i,2)
        bb22(i,1) =  0.
      ENDIF  ! (wk_adv(i))
    ENDDO

!!  printing the tridiagonal matrix
!!!  First row
!!   k = 1
!!   print 1789, k, aa11(1,1), aa12(1,1), bb11(1,1), bb12(1,1)
!!   print 1789, k, aa21(1,1), aa22(1,1), bb21(1,1), bb22(1,1)
!!1789 FORMAT(1X, I3, 3(4X, 2E13.5))
!!        coef = dtime/(ph(1,k)-ph(1,k+1))
!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2) ', &
!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2)
!!
!!!  Rows 2 to klev-1
!!   DO k = 2, klev-1
!!     print 1789, k, cc11(1,k), cc12(1,k), aa11(1,k), aa12(1,k), bb11(1,k), bb12(1,k)
!!     print 1789, k, cc21(1,k), cc22(1,k), aa21(1,k), aa22(1,k), bb21(1,k), bb22(1,k)
!!        coef = dtime/(ph(1,k)-ph(1,k+1))
!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1) ', &
!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1)
!!   ENDDO
!!
!!!  Row klev
!!     print 1789, klev, cc11(1,klev), cc12(1,klev), aa11(1,klev), aa12(1,klev)
!!     print 1789, klev, cc21(1,klev), cc22(1,klev), aa21(1,klev), aa22(1,klev)
!!        coef = dtime/(ph(1,klev)-ph(1,klev+1))
!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev) ', &
!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev)


!! Downward loop 

   xt1(:,:) = thb(:,:)     
   xt2(:,:) = delta_th(:,:)
   xq1(:,:) = qb(:,:)      
   xq2(:,:) = delta_q(:,:) 

    DO i = 1,klon
      IF (wk_adv(i)) THEN 
        k = kupper(i)
        alpha11(:,k)=aa11(:,k)
        alpha12(:,k)=aa12(:,k)
        alpha21(:,k)=aa21(:,k)
        alpha22(:,k)=aa22(:,k)
        beta11(:,k)=0.
        beta12(:,k)=0.
        beta21(:,k)=0.
        beta22(:,k)=0.
        yt1(i,k) = xt1(i,k)
        yt2(i,k) = xt2(i,k)
        yq1(i,k) = xq1(i,k)
        yq2(i,k) = xq2(i,k)
      ENDIF  ! (wk_adv(i))
    ENDDO
    DO i = 1,klon
      IF (wk_adv(i)) THEN 
        k = kupper(i)
        det=alpha11(i,k)*alpha22(i,k) - alpha12(i,k)*alpha21(i,k)
        ai11(i,k)= alpha22(i,k)/det
        ai12(i,k)=-alpha12(i,k)/det
        ai21(i,k)=-alpha21(i,k)/det
        ai22(i,k)= alpha11(i,k)/det 
        zt1(i,k) = ai11(i,k)*yt1(i,k) + ai12(i,k)*yt2(i,k)
        zt2(i,k) = ai21(i,k)*yt1(i,k) + ai22(i,k)*yt2(i,k)
        zq1(i,k) = ai11(i,k)*yq1(i,k) + ai12(i,k)*yq2(i,k)
        zq2(i,k) = ai21(i,k)*yq1(i,k) + ai22(i,k)*yq2(i,k)
      ENDIF  ! (wk_adv(i))
    ENDDO

    DO k = klev, 2, -1
      DO i = 1,klon
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          gamma11(i,k) = ai11(i,k)*cc11(i,k) + ai12(i,k)*cc21(i,k)
          gamma12(i,k) = ai11(i,k)*cc12(i,k) + ai12(i,k)*cc22(i,k)
          gamma21(i,k) = ai21(i,k)*cc11(i,k) + ai22(i,k)*cc21(i,k)
          gamma22(i,k) = ai21(i,k)*cc12(i,k) + ai22(i,k)*cc22(i,k)
   
          alpha11(i,k-1) = aa11(i,k-1) - ( bb11(i,k-1)*gamma11(i,k)+bb12(i,k-1)*gamma21(i,k) )
          alpha12(i,k-1) = aa12(i,k-1) - ( bb11(i,k-1)*gamma12(i,k)+bb12(i,k-1)*gamma22(i,k) )
          alpha21(i,k-1) = aa21(i,k-1) - ( bb21(i,k-1)*gamma11(i,k)+bb22(i,k-1)*gamma21(i,k) )
          alpha22(i,k-1) = aa22(i,k-1) - ( bb21(i,k-1)*gamma12(i,k)+bb22(i,k-1)*gamma22(i,k) )
   
          beta11(i,k-1) = bb11(i,k-1)*ai11(i,k)+bb12(i,k-1)*ai21(i,k)
          beta12(i,k-1) = bb11(i,k-1)*ai12(i,k)+bb12(i,k-1)*ai22(i,k)
          beta21(i,k-1) = bb21(i,k-1)*ai11(i,k)+bb22(i,k-1)*ai21(i,k)
          beta22(i,k-1) = bb21(i,k-1)*ai12(i,k)+bb22(i,k-1)*ai22(i,k)
   
          yt1(i,k-1) = xt1(i,k-1) - ( beta11(i,k-1)*yt1(i,k) +beta12(i,k-1)*yt2(i,k) )
          yt2(i,k-1) = xt2(i,k-1) - ( beta21(i,k-1)*yt1(i,k) +beta22(i,k-1)*yt2(i,k) )
          yq1(i,k-1) = xq1(i,k-1) - ( beta11(i,k-1)*yq1(i,k) +beta12(i,k-1)*yq2(i,k) )
          yq2(i,k-1) = xq2(i,k-1) - ( beta21(i,k-1)*yq1(i,k) +beta22(i,k-1)*yq2(i,k) )
   
          det=alpha11(i,k-1)*alpha22(i,k-1) - alpha12(i,k-1)*alpha21(i,k-1)
          ai11(i,k-1)= alpha22(i,k-1)/det
          ai12(i,k-1)=-alpha12(i,k-1)/det
          ai21(i,k-1)=-alpha21(i,k-1)/det
          ai22(i,k-1)= alpha11(i,k-1)/det 
   
          zt1(i,k-1) = ai11(i,k-1)*yt1(i,k-1)+ai12(i,k-1)*yt2(i,k-1)
          zt2(i,k-1) = ai21(i,k-1)*yt1(i,k-1)+ai22(i,k-1)*yt2(i,k-1)
          zq1(i,k-1) = ai11(i,k-1)*yq1(i,k-1)+ai12(i,k-1)*yq2(i,k-1)
          zq2(i,k-1) = ai21(i,k-1)*yq1(i,k-1)+ai22(i,k-1)*yq2(i,k-1)
        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
      ENDDO
    ENDDO
        
!! Upward loop 

    DO i = 1,klon
      IF (wk_adv(i)) THEN 
       th1(i,1) = zt1(i,1)
       th2(i,1) = zt2(i,1)
       q1(i,1)  = zq1(i,1)
       q2(i,1)  = zq2(i,1)
     
       d_tb_dadv(i,1) =     ( rr11(i)*(th1(i,1)-thw(i,1))+rr12(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1)
       d_deltat_dadv(i,1) = ( rr21(i)*(th1(i,1)-thw(i,1))+rr22(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1)
       d_qb_dadv(i,1) =       rr11(i)*(q1(i,1) -qw(i,1)) +rr12(i)*(q2(i,1)-qx(i,1))
       d_deltaq_dadv(i,1) =   rr21(i)*(q1(i,1) -qw(i,1)) +rr22(i)*(q2(i,1)-qx(i,1))
      ENDIF  ! (wk_adv(i))
    ENDDO

    DO k = 2, klev
      DO i = 1,klon
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          th1(i,k) = zt1(i,k) - ( gamma11(i,k)*th1(i,k-1)+gamma12(i,k)*th2(i,k-1) )
          th2(i,k) = zt2(i,k) - ( gamma21(i,k)*th1(i,k-1)+gamma22(i,k)*th2(i,k-1) )
          q1(i,k)  = zq1(i,k) - ( gamma11(i,k)*q1(i,k-1) +gamma12(i,k)*q2(i,k-1) )
          q2(i,k)  = zq2(i,k) - ( gamma21(i,k)*q1(i,k-1) +gamma22(i,k)*q2(i,k-1) )
   
          d_tb_dadv(i,k) =     ( rr11(i)*(th1(i,k)-thw(i,k))+rr12(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k)
          d_deltat_dadv(i,k) = ( rr21(i)*(th1(i,k)-thw(i,k))+rr22(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k)
          d_qb_dadv(i,k) =       rr11(i)*(q1(i,k)-qw(i,k))  +rr12(i)*(q2(i,k)-qx(i,k))
          d_deltaq_dadv(i,k) =   rr21(i)*(q1(i,k)-qw(i,k))  +rr22(i)*(q2(i,k)-qx(i,k))
        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
      ENDDO
    ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!       Verification de l'inversion                !!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!!    DO i = 1,klon
!!        xt1inv(i,1) = aa11(i,1)*th1(i,1) + aa12(i,1)*th2(i,1) + bb11(i,1)*th1(i,2) + bb12(i,1)*th2(i,2)
!!        xt2inv(i,1) = aa21(i,1)*th1(i,1) + aa22(i,1)*th2(i,1) + bb21(i,1)*th1(i,2) + bb22(i,1)*th2(i,2)
!!        xq1inv(i,1)  = aa11(i,1)*q1(i,1)  + aa12(i,1)*q2(i,1)  + bb11(i,1)*q1(i,2)  + bb12(i,1)*q2(i,2) 
!!        xq2inv(i,1)  = aa21(i,1)*q1(i,1)  + aa22(i,1)*q2(i,1)  + bb21(i,1)*q1(i,2)  + bb22(i,1)*q2(i,2) 
!!    ENDDO
!!   
!!      DO k = 2, klev-1
!!        DO i = 1,klon
!!        xt1inv(i,k) = aa11(i,k)*th1(i,k) + aa12(i,k)*th2(i,k) + bb11(i,k)*th1(i,k+1) + bb12(i,k)*th2(i,k+1) &
!!                                                              + cc11(i,k)*th1(i,k-1) + cc12(i,k)*th2(i,k-1)
!!        xt2inv(i,k) = aa21(i,k)*th1(i,k) + aa22(i,k)*th2(i,k) + bb21(i,k)*th1(i,k+1) + bb22(i,k)*th2(i,k+1) &
!!                                                              + cc21(i,k)*th1(i,k-1) + cc22(i,k)*th2(i,k-1)
!!        xq1inv(i,k)  = aa11(i,k)*q1(i,k)  + aa12(i,k)*q2(i,k)  + bb11(i,k)*q1(i,k+1)  + bb12(i,k)*q2(i,k+1)  &
!!                                                              + cc11(i,k)*q1(i,k-1)  + cc12(i,k)*q2(i,k-1)
!!        xq2inv(i,k)  = aa21(i,k)*q1(i,k)  + aa22(i,k)*q2(i,k)  + bb21(i,k)*q1(i,k+1)  + bb22(i,k)*q2(i,k+1)  &
!!                                                              + cc21(i,k)*q1(i,k-1)  + cc22(i,k)*q2(i,k-1)
!!        ENDDO
!!      ENDDO
!!   
!!    DO i = 1,klon
!!        xt1inv(i,klev) = aa11(i,klev)*th1(i,klev) + aa12(i,klev)*th2(i,klev) + cc11(i,klev)*th1(i,klev-1) + cc12(i,klev)*th2(i,klev-1)
!!        xt2inv(i,klev) = aa21(i,klev)*th1(i,klev) + aa22(i,klev)*th2(i,klev) + cc21(i,klev)*th1(i,klev-1) + cc22(i,klev)*th2(i,klev-1)
!!        xq1inv(i,klev)  = aa11(i,klev)*q1(i,klev)  + aa12(i,klev)*q2(i,klev)  + cc11(i,klev)*q1(i,klev-1)  + cc12(i,klev)*q2(i,klev-1) 
!!        xq2inv(i,klev)  = aa21(i,klev)*q1(i,klev)  + aa22(i,klev)*q2(i,klev)  + cc21(i,klev)*q1(i,klev-1)  + cc22(i,klev)*q2(i,klev-1) 
!!    ENDDO
!!   
!!    DO k = 1, 20
!!      IF (abs(xt1inv(1,k)-xt1(1,k)) .GT. 1.e-15*xt1(1,k) ) THEN
!!        print *,'wake_dadv, k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k) ', &
!!                            k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k)
!!      ENDIF
!!      IF (abs(xt2inv(1,k)-xt2(1,k)) .GT. 1.e-15*xt2(1,k) ) THEN
!!        print *,'wake_dadv, k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k) ', &
!!                            k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k)
!!      ENDIF
!!      IF (abs(xq1inv(1,k)-xq1(1,k)) .GT. 1.e-15*xq1(1,k) ) THEN
!!        print *,'wake_dadv, k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k) ', &
!!                            k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k)
!!      ENDIF
!!      IF (abs(xq2inv(1,k)-xq2(1,k)) .GT. 1.e-15*xq2(1,k) ) THEN
!!        print *,'wake_dadv, k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k) ', &
!!                          k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k)
!!      ENDIF
!!    ENDDO

  ELSE  ! (flag_dadv_implicit)

!   Explicit scheme : compute directly d_deltat_dadv and d_tb_dadv 
!                     (and similarly for d_deltaq_dadv and d_qb_dadv).
! -----------------------------------------------------------------------------------------

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        d_thx(i, 1) = 0.
        d_thw(i, 1) = 0.
        d_dth(i, 1) = 0.
        d_qx(i, 1) = 0.
        d_qw(i, 1) = 0.
        d_dq(i, 1) = 0.
      END IF
    END DO

    DO k = 2, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
          d_thx(i, k) = thx(i, k-1) - thx(i, k)
          d_thw(i, k) = thw(i, k-1) - thw(i, k)
          d_dth(i, k) = delta_th(i, k-1) - delta_th(i, k)
          d_qx(i, k) = qx(i, k-1) - qx(i, k)
          d_qw(i, k) = qw(i, k-1) - qw(i, k)
          d_dq(i, k) = delta_q(i, k-1) - delta_q(i, k)
        END IF ! (wk_adv(i) .AND. k<=kupper(i)+1)
      END DO
    END DO

    DO k = 1, klev-1
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
          d_deltat_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
            (rr22(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k) - &
             rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1) )*ppi(i, k) - &
             dtime*entr_s(i,k)*delta_th(i,k)/sigmaw(i)*ppi(i, k)
!
          d_deltaq_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
            (rr22(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- &
             rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1) ) - &
             dtime*entr_s(i,k)*delta_q(i,k)/sigmaw(i)

          ! and increment large scale tendencies
          d_tb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k)- &
                                  rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1))/ &
                                 (ph(i,k)-ph(i,k+1)) &
                                 -sigmaw(i)*(1.-sigmaw(i))*delta_th(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)

          d_qb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- &
                                  rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1))/ &
                                 (ph(i,k)-ph(i,k+1)) &
                                 -sigmaw(i)*(1.-sigmaw(i))*delta_q(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
                                 (ph(i,k)-ph(i,k+1)) )
        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
          d_tb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_thx(i, k)/(ph(i, k)-ph(i, k+1)))*ppi(i, k)

          d_qb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_qx(i, k)/(ph(i, k)-ph(i, k+1)))
        END IF ! (wk_adv(i) .AND. k<=kupper(i)-1)
      END DO
    END DO

  ENDIF! (flag_dadv_implicit)

!print *,'ZZZwake_dadv k, d_deltat_dadv(1,k) ', (k, d_deltat_dadv(1,k), k = 1,3)

    END SUBROUTINE wake_dadv
END MODULE lmdz_wake_dadv
