subroutine lwi(nl,nnu,netrad,deltab,deltap,ksi,coolrate)


      implicit none


#include "dimensions.h"
#include "dimphy.h"
#include "comg1d.h"  !pour grads 1D
#include "YOMCST.h"
#include "timerad.h"
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                             -   lwi    -    
C
C     PURPOSE:       Schema semi - implicite 
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


c************************************************************************
c
c        0.    Declarations
c              ------------
c
c-------------------------------------------------------------------------
c        0.1   Arguments
c              ---------
c Inputs
c.......

      integer nl,nnu
      real    ksi(0:nl+1,0:nl+1,nnu) !ksi factors  
      real    netrad(0:nl)     ! radiative budget (W/m**2)
      real    deltab(nl,nnu) !gradient de fonction de Planck
      real    deltap(nl) !epaisseur de la couche en pression (Pa) 
      
c Outputs
c........

      real    coolrate(nl)
      
c-------------------------------------------------------------------------
c        0.2   local arrays
c              ------------
c
      real    di(nl) 
     .      , spec(nl) 
     .      , hi(nl) 
     .      , bi(nl) 

      real    ci(nl) 
     .      , ai(nl) 

      real    deltat

      real    semi,semit,denom

      integer i,jl,band

c************************************************************************
c
c        1.    Initialisations
c              ---------------
c
c-----------------------------------------------------------------------
        deltat = dtimerad
c       print*,'SEMI = ',semi, '(expl:0  semi-implicite:0.5  impl:1)'
c        semi  = 0.5
        semi = 0.5
c       semi = 1.
        semit = semi * deltat

c       print*,'dtimerad,deltat,semit:',dtimerad,deltat,semit

c************************************************************************
c
c        2.    
c              ---------------
c
c-------------------------------------------------------------------------
c        2.1   Calcul des di
c              -------------
c

      do i = 1 , nl-1
        spec(i) = 0.0e0
        do band = 1, nnu
          spec(i) = spec(i) + 
     .      ( ksi(i,nl+1,band)
     .      + ksi(i,i+1,band)
     .      + ksi(i,i-1,band) ) * deltab(i,band)
        enddo
      
        di(i) =  1. + semit * (RG / deltap(i) / RCPD) * spec(i)
c     print*,i,' di(i)=',di(i)
      enddo
         

c couche nl
c ------------
c      , on enleve i,i+1 sinon on a 2 fois le cooling2space
 
        spec(nl) = 0.0e0
        do band = 1, nnu
          spec(nl) = spec(nl) +
     .      ( ksi(nl,nl+1,band)
     .      + ksi(nl,nl-1,band) ) * deltab(nl,band)
        enddo
     
        di(nl) =  1. + semit * (RG / deltap(nl) / RCPD) * spec(nl)

c-------------------------------------------------------------------------
c        2.2   Calcul des hi
c              -------------
c

      do i = 1 , nl-1
      
        spec(i) = 0.0e0
        do band = 1, nnu
          spec(i) = spec(i) + ksi(i,i+1,band) * deltab(i+1,band)
        enddo

        hi(i) =  - semit * (RG / deltap(i) / RCPD) * spec(i)

      enddo
                        
c     print*,'hi(i)',hi(i)

c-------------------------------------------------------------------------
c        2.3   Calcul des bi
c              -------------
c


      do i = 2 , nl
      
        spec(i) = 0.0e0
        do band = 1, nnu
          spec(i) = spec(i) + ksi(i,i-1,band) * deltab(i-1,band)
        enddo
       
        bi(i) =   - semit * (RG / deltap(i) / RCPD) * spec(i)
        
      enddo

c     print*,'bi(i)',bi(i)

c couche 1
c --------
c  tant qu'on a pas un calcul propre de deltab(0) qui tienne compte de 
c    la discontinuite de temperature au sol , on met  b1 = 0

      
        bi(1) = 0 

c-------------------------------------------------------------------------
c        2.4   
c              -------------
c

c couche nl
c ------------

c     ci(nl) = ((RG/RCPD) * netrad(nl) / deltap(nl))
c    .                   / di(nl)
      ci(nl) = netrad(nl) *RG/RCPD / deltap(nl)

      ai(nl) = - bi(nl) / di(nl)

      do i = nl-1 , 1 , -1
        denom = di(i) + hi(i) * ai(i+1)

c       ci(i) = (  (RG/RCPD) * netrad(i) / deltap(i)
c    .               - hi(i) * ci(i+1)  )  / denom
        ci(i) = netrad(i) *RG/RCPD / deltap(i)
 
        ai(i) = -bi(i) / denom
      enddo

c-------------------------------------------------------------------------
c        2.5   
c              -------------
c

      coolrate(1) = ci(1)

      do i = 2 , nl
           coolrate(i) = ci(i) + ai(i) * coolrate(i-1)
c       print*,i,' coolrate(i)=',coolrate(i)
      enddo

c-------------------------------------------------------------------------

      RETURN
      END