!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
***S/P LIN_DIFF_VERT1
*

      SUBROUTINE LIN_DIFF_VERT1 (TU, U, KU, GU, ALFA, BETA, S, SK, 4,2
     %                           TAU, F, A, B, C, D, N, NK)
*
#include "impnone.cdk"
      INTEGER N, NK
      REAL TU(N, NK), U(N, NK), KU(N, NK), GU(N, NK)
      REAL ALFA(N), BETA(N), S(N,NK), SK(N,NK), TAU, F
      REAL A(N, NK), B(N, NK), C(N, NK), D(N, NK)
*
*Author
*          S. Laroche (oct 2000)
*
*Object
*          to solve a vertical diffusion equation by finite
*          differences (originally from DIFUVDFj, R. Benoit)
*
*Arguments
*
*          - Output -
* TU       U tendency (D/DT U) due to the vertical diffusion and to
*          term R
*
*          - Input -
* U        variable to diffuse (U,V,T,Q,E)
* KU       diffusion coefficient
* GU       optional countergradient term
* ALFA     inhomogeneous term for the surface flux
* BETA     homogeneous term for the surface flux
* S        sigma coordinates of full levels
* SK       sigma coordinates of diffusion coefficient levels
* TAU      length of timestep
* F        waiting factor for time 'N+1'
* A        work space (N,NK)
* B        work space (N,NK)
* C        work space (N,NK)
* D        work space (N,NK)
* N        number of columns to process
* NK       vertical dimension
*
*Notes
*          D/DT U = D(U)
*          D(U) = D/DS J(U)
*          J(U) = KU*(D/DS U + GU)
*
**
*
      INTEGER I, K
      REAL SC, SCTU, HM, HP, HD, KUM, KUP, SCK1
      EXTERNAL LIN_DIFUVD1, LIN_DIFUVD2
*
      SC   = 1.0
      SCTU = 0
      SCK1 = 1
*
* (1) CONSTRUIRE L'OPERATEUR TRIDIAGONAL DE DIFFUSION N=(A,B,C)
*                ET LE TERME CONTRE-GRADIENT (DANS D)
*
*     K=1
*
      HM=0
      DO I=1,N
          HP=S(I,2)-S(I,1)
          HD=0.5*(S(I,1)+S(I,2))-S(I,1)
          A(I,1)=0
          B(I,1)=-SCK1*KU(I,1)/(HP*HD)
          C(I,1)=-SCK1*B(I,1)
          D(I,1)=SCK1*KU(I,1)*GU(I,1)/HD
      ENDDO
*
*     K=2...NK-1
*
      DO K=2,NK-1,1
         DO I=1,N
            HM=S(I,K)-S(I,K-1)
            HP=S(I,K+1)-S(I,K)
            HD=0.5*(HM+HP)
            A(I,K)=KU(I,K-1)/(HM*HD)
            B(I,K)=-(KU(I,K-1)/HM +KU(I,K)/HP)/HD
            C(I,K)=KU(I,K)/(HP*HD)
            D(I,K)=(KU(I,K)*GU(I,K)-KU(I,K-1)*GU(I,K-1))/HD
         ENDDO
      ENDDO
*
*     K=NK
*
      HP=0
      DO I=1,N
         HM=S(I,NK)-S(I,NK-1)
         HD=SK(I,NK)-0.5*(S(I,NK-1)+S(I,NK))
         A(I,NK)=KU(I,NK-1)/(HM*HD)
         B(I,NK)=-(KU(I,NK-1)/HM + 0)/HD
         C(I,NK)=0
         D(I,NK)=(0-KU(I,NK-1)*GU(I,NK-1))/HD
      ENDDO
*
*
* (2) CALCULER LE COTE DROIT D=TAU*(SC*N(U)+R+D/DS(KU*GU))
*
      CALL LIN_DIFUVD1 (D, SC, A, B, C, U, D, N, NK)
      DO K=1,NK
         DO I=1,N
           D(I,K)  = TAU*D(I,K)
         ENDDO
      ENDDO
*
* (3) CALCULER OPERATEUR DU COTE GAUCHE
*
      DO K=1,NK
         DO I=1,N
            A(I,K)=  -SC*TAU*A(I,K)
            B(I,K)=1.-SC*TAU*B(I,K)
            C(I,K)=  -SC*TAU*C(I,K)
         ENDDO
      ENDDO
*
* (4) AJOUTER TERME DE FLUX DE SURFACE POUR TYPE='U'/'UT'
*
      DO I=1,N
         HD=SK(I,NK)-0.5*(S(I,NK-1)+S(I,NK))
         B(I,NK)=B(I,NK)-F*TAU*BETA(I)/HD
         D(I,NK)=D(I,NK)+(ALFA(I)+BETA(I)*U(I,NK))*TAU/HD
      ENDDO
*
*
* (5) RESOUDRE SYSTEME TRIDIAGONAL [A,B,C] X = D. METTRE X DANS TU.
*
      CALL LIN_DIFUVD2 (TU, A, B, C, D, D, N, NK)
*
* (6) OBTENIR TENDANCE
*
      DO K=1,NK
         DO I=1,N
          TU(I,K)  = TU(I,K)/TAU
         ENDDO
      ENDDO
*
      RETURN
      END