SUBROUTINE load_ksi(ksive)
      
      use dimphy
      IMPLICIT none

#include "dimensions.h"
#include "YOMCST.h"
#include "comcstVE.h"
C
C     ------------------------------------------------------------------
C
C     PURPOSE.
C     --------
C
c     This routine loads the longwave matrix of factors Ksi
c     
c     The Ksi matrixes have been computed by Vincent Eymet
C
C     AUTHOR.
C     -------
C        Sebastien Lebonnois
C
C     MODIFICATIONS.
C     --------------
C        version multimatrice (topographie, sommet nuages): 20/12/2006
C     ------------------------------------------------------------------
C
C* ARGUMENTS:
C
c inputs
      real   psurf(klon)           ! Surface pressure
      real   ztop(klon)            ! Altitude of the top of cloud deck (km)
c outputs
      real   ksive(0:klev+1,0:klev+1,nnuve,nbmat)  ! ksi matrixes in Vincent's file

c local variables
      integer i,j,ig,band,pve,nlve
      integer mat,Nb,m,Nmat,nl_init,mat0
      parameter(nl_init=8)
      character*9 tmp1
      character*100 file
      real   lambda(nnuve)            ! wavelenght in table (mu->m, middle of interval)
      real   lambdamin(nnuve),lambdamax(nnuve) ! in microns
      real   dlambda                  ! cm-1  

      nlve = klev

c ------------------------
c Loading the ksi file
c ------------------------

      file = "ksi_global.txt"
      open(10,file=file)
      
      do i=1,nl_init-1
         read(10,*)
      enddo
      read(10,*) (tmp1,i=1,4),Nmat

      if (nbmat.ne.Nmat) then
         write(*,*) 'This is subroutine load_ksi'
         print*,'Probleme de dimension entre ksi.txt et le param nbmat'
         print*,'Nb matrices = ',nbmat,Nmat
         stop
      endif

      do mat=1,nbmat
        read(10,*)
        read(10,*) 
        read(10,*) (tmp1,j=1,2),pve
        psurfve(mat) = pve*1.e5  ! pve en bar, psurfve en Pa
        read(10,*) (tmp1,j=1,7),ztopve(mat)
        ztopve(mat) = ztopve(mat)*1.e-3 ! passage en km
        read(10,*)
        read(10,*) m,Nb
        if (m.ne.nlve) then
         write(*,*) 'This is subroutine load_ksi'
         print*,'Probleme de dimension entre ksi.txt et le param nlve'
         print*,'N levels = ',m,nlve
         stop
        endif
        if (Nb.ne.nnuve) then
         write(*,*) 'This is subroutine load_ksi'
         print*,'Probleme de dimension entre ksi.txt et le param nnuve'
         print*,'N freq = ',Nb,nnuve
         stop
        endif
c     Now reading ksi matrix index "mat"
        do band=1,Nb
         read(10,*) lambdamin(band),lambdamax(band)
         do i=0,m+1
            read(10,'(100e17.9)') (ksive(i,j,band,mat),j=0,m+1) ! sr/µm/cm¯¹
         enddo                  ! i
        enddo                     ! band
c       print*,"Matrice ",mat," lue"
c       print*,"   psurf=",psurfve(mat)," bars, Ztop=",ztopve(mat)," km"
      enddo  ! mat
      
      close(10)

c longueur d'onde centrale et largeur de chaque bande
      do band=1,nnuve
         lambda(band)=(lambdamin(band)+lambdamax(band))/2.*1.e-6   ! en m
         dlambda     =(1./lambdamin(band)-1./lambdamax(band))*1.e4 ! en cm-1
c        print*,band,lambdamin(band),dlambda,lambdamax(band)

c changement de convention (signe) pour ksi,
c et prise en compte de la largeur de bande (en cm-1):
         do mat=1,nbmat
         do i=0,nlve+1
           do j=0,nlve+1
              ksive(i,j,band,mat) = -ksive(i,j,band,mat)*dlambda
           enddo
         enddo
         enddo
c calcul des coeff al et bl pour luminance Planck 
         al(band) = 2.*RHPLA*RCLUM*RCLUM/(lambda(band))**5.
c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul
c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6
     .                * 1.e-6
         bl(band) = RHPLA*RCLUM/(RKBOL*lambda(band))
      enddo
      
      print*,"LOAD_KSI OK"

      return
      end