SUBROUTINE load_psi( 
     S              psurf, ztop, ksive,
     S              temp, psimap, deltapsimap)
      
      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     (interpolated for a given column)
c     to build the Net Exchange Rates matrix Psi.
c     Psi(i,j,nu) = Ksi(i,j,nu) * ( B(i,nu)-B(j,nu) )
c     
c     The Ksi matrixes have been computed by Vincent Eymet
C
c     The NER matrix is then integrated in frequency.
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)
      real   ksive(0:kflev+1,0:kflev+1,nnuve,nbmat)  ! ksi matrixes in Vincent's file
      real   temp(klon,0:kflev+1)  ! Temperature in layer (K)
c outputs
      real   psimap(0:kflev+1,0:kflev+1,klon)
      real   deltapsimap(0:kflev+1,0:kflev+1,klon)

c local variables
      integer i,j,ig,band,nlve
      integer mat,m,mat0
      character*100 file
      real   bplck(0:kflev+1,nnuve)    ! Planck luminances in table layers
      real   y(0:kflev,nnuve)          ! intermediaire Planck
      real   zdblay(0:kflev+1,nnuve)   ! gradient en temperature de planck
      real   ksi
      real   factp,factz

      nlve = kflev  ! (doit correspondre, pour bplck, y et zdblay)

c -----------------------
c Main loop on grid point
c -----------------------
      do ig=1,klon

c Planck function
c ---------------

       do band=1,nnuve
         do j=0,nlve
c B(T,l) = al/(exp(bl/T)-1)
           y(j,band) = exp(bl(band)/temp(ig,j))-1.
           bplck(j,band) = al(band)/(y(j,band))
           zdblay(j,band) = al(band)*bl(band)*exp(bl(band)/temp(ig,j))/
     .                 ((temp(ig,j)*temp(ig,j))*(y(j,band)*y(j,band)))
         enddo
         bplck(nlve+1,band) = 0.0
         zdblay(nlve+1,band)= 0.0
       enddo

c interpolating ksi 
c  and
c computing psi and deltapsi
c ---------------------------

c init
       do j=0,nlve+1
        do i=0,nlve+1
         psimap(i,j,ig) = 0.0  ! positif quand nrj de i->j 
         deltapsimap(i,j,ig) = 0.0
        enddo
       enddo

c finding the right matrixes
       mat0 = 0
       do mat=1,nbmat-nbztopve
         if (  (psurfve(mat).ge.psurf(ig))
     .    .and.(psurfve(mat+nbztopve).lt.psurf(ig))
     .    .and.(ztopve(mat).lt.ztop(ig))
     .    .and.(ztopve(mat+1).ge.ztop(ig)) ) then
              mat0  = mat
c             print*,'ig=',ig,'  mat0=',mat
              factp = (psurf(ig)            -psurfve(mat))
     .               /(psurfve(mat+nbztopve)-psurfve(mat))
              factz = (ztop(ig)     -ztopve(mat))
     .               /(ztopve(mat+1)-ztopve(mat))
         endif
       enddo
       if (mat0.eq.0) then
         write(*,*) 'This is subroutine load_psi'
         print*,'Probleme pour interpolation au point ig=',ig
         print*,'psurf = ',psurf(ig),' ztop = ',ztop(ig)
         stop
       endif

c------------------
c !!TEST!! Matrice unique fixee: psurf = 90 bar, ztop = 70
c      mat0 = 24
c      print*,'MATRICE UNIQUE: ',mat0,' / ps=',psurfve(mat0),
c    .                                ' / ztop=',ztopve(mat0)
c------------------

c interpolation of ksi and computation of psi,deltapsi
       do band=1,nnuve
        do j=0,nlve+1
         do i=0,nlve+1
          ksi = ksive(i,j,band,mat0)*(1-factz)*(1-factp)
     .         +ksive(i,j,band,mat0+1)*factz  *(1-factp)
     .         +ksive(i,j,band,mat0+nbztopve)*(1-factz)*factp
     .         +ksive(i,j,band,mat0+nbztopve+1)*factz  *factp
         psimap(i,j,ig) = psimap(i,j,ig) + 
     .          ksi*(bplck(i,band)-bplck(j,band))
         deltapsimap(i,j,ig) = deltapsimap(i,j,ig) + 
     .          ksi*zdblay(i,band)
         enddo
        enddo
       enddo

      enddo !ig
c -----------------------
c End loop on grid point 
c -----------------------

c     print*,"LOAD_PSI OK"

      return
      end