!=================================================================
!  地表面温度計算プログラム: 拡散方程式, Crank-Nicolson 法
!
!  2000/09/21  小高正嗣
!=================================================================

      subroutine CLTSFC
     I   ( FHEAT , FXRADD, FXSOL,
     O     TEMPG , TSFC  ,
     C     SIGMA , ALBEDO, CPGRND, 
     C     CDTURBG,DTIME , DENSG ,
     G     GDZ   , GDZM  , 
     D     NX    , NXB   , NZ   , NZB   ,
     D     NZG   , NZGB                   )

!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      

      dimension FHEAT ( 2 )
      dimension FXRADD( -2:NZB )
      dimension TSFC ( -2:NXB )
      dimension TEMPG( NZB )
      dimension GDZ(0:NZGB), GDZM(0:NZG)

      dimension
     \     ELMG_FOR(NZG,NZG),   ! 係数行列1
     \     ELMG_BCK(NZG,NZG)    ! 係数行列2

      dimension                 ! LAPACK 用作業配列
     \     DGL(NZG-1), 
     \     DG(NZG), 
     \     DGU(NZG-1), 
     \     DGU2(NZG-2),
     \     BG(NZG),
     \     IPIVG(NZG)

!---------------------------------------------------------------------
! 係数行列の設定

      call SETMTRXG
     I   ( CDTURBG, 
     O     ELMG_FOR, ELMG_BCK, DGL, DG, DGU,
     G     GDZ , GDZM,
     D     DTIME, NZG, NZGB      )

!---------------------------------------------------------------------
! 係数行列 A の分解

      call DGTTRF
     \   ( NZG, DGL, DG, DGU, DGU2, IPIVG, INFO)

!-----------------------------------------------------------------      
! 境界の熱フラックス

      HEAT0 = FXSOL * ( 1.0D0 - ALBEDO )
     \        - SIGMA * TEMPG( NZG )**4.0D0 
     \        + FXRADD( 0 )
     \        - FHEAT( 1 )

!---------------------------------------------------------------------
! 境界値の寄与を先に計算

      TEMPG(NZG) = TEMPG(NZG) 
     \     + DTIME / DENSG / CPGRND / GDZ(NZG) * HEAT0

!---------------------------------------------------------------------
! 拡散項の計算

      call CLDIFG
     V   ( TEMPG, 
     I     ELMG_BCK,
     W     DGL, DG, DGU, DGU2, BG, IPIVG,
     D     NZG )

!---------------------------------------------------------------------
      do IX = -2, NXB
         TSFC( IX ) = TEMPG( NZG )
      end do
!---------------------------------------------------------------------

      return

      end

