SUBROUTINE LW_venus_ve( 
     S              PPB, PT, PTSURF,
     S              PCOOL, 
     S              PTOPLW,PSOLLW,PSOLLWDN,
     S              ZFLNET)
      
      IMPLICIT none

#include "dimensions.h"
#include "dimphy.h"
#include "raddim.h"
#include "YOMCST.h"
C
C     ------------------------------------------------------------------
C
C     PURPOSE.
C     --------
C
c     This routine loads the longwave matrix of factors Ksi,
c     used 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     This Ksi matrix has been computed by Vincent Eymet
C
c     The NER matrix is then integrated in frequency, and the output
c     are calculated.
c
C     AUTHOR.
C     -------
C        Sebastien Lebonnois
C
C     MODIFICATIONS.
C     --------------
C        ORIGINAL : 27/07/2005
C     ------------------------------------------------------------------
C
C* ARGUMENTS:
C
c inputs

      REAL   PPB(KFLEV+1)  ! inter-couches PRESSURE (bar)
      REAL   PT(KFLEV)     ! Temperature in layer (K)
      REAL   PTSURF        ! Surface temperature
C
c output

      REAL   PCOOL(KFLEV) ! LONGWAVE COOLING (K/VENUSDAY) within each layer
      REAL   PTOPLW       ! LONGWAVE FLUX AT T.O.A. (net, + vers le haut)
      REAL   PSOLLW       ! LONGWAVE FLUX AT SURFACE (net, + vers le haut)
      REAL   PSOLLWDN     ! LONGWAVE FLUX AT SURFACE (down, + vers le bas)
      REAL   ZFLNET(KFLEV+1) ! net thermal flux at ppb levels (+ vers le haut)

C
C* LOCAL VARIABLES:
C
      integer nlve,nnuve
      parameter (nlve=kflev) ! fichiers Vincent: same grid than GCM
      parameter (nnuve=68)   ! fichiers Vincent et Bullock
      real    dureejour
      parameter (dureejour=10.087e6)
      
      integer i,j,p,nl0,nnu0,band
      real   lambda(nnuve)    ! wavelenght in table (mu->m, middle of interval)
      real   ksive(0:nlve+1,0:nlve+1,nnuve) ! ksi factors
      real   bplck(0:nlve+1,nnuve)          ! Planck luminances in table layers
      real   al(nnuve),bl(nnuve)            ! for Planck luminances calculations
      real   psive(0:nlve+1,0:nlve+1,nnuve) ! NER in W/m**2 per wavelength band
      real   psi_1(0:nlve+1,0:nlve+1)       ! NER in W/m**2 (sum on lambda)

      real   ztemp(0:nlve)    ! GCM temperature in table layers
      real   zlnet(nlve+1)    ! net thermal flux (W/m**2)
      real   dzlnet(0:nlve)   ! Radiative budget (W/m**2)
      character*22 nullchar
      real   lambdamin,lambdamax ! in microns
      real   dlambda             ! cm-1  

      real   y(0:nlve,nnuve)    ! intermediaire Planck
      real   pdp(nlve)          ! epaisseur de la couche en pression (Pa)
      real   zdblay(nlve,nnuve)     ! gradient en temperature de planck
      
      logical firstcall
      data    firstcall/.true./

      save   lambda,ksive,al,bl,firstcall
      
c ------------------------
c Loading the files
c ------------------------

      if (firstcall) then

c Matrice Ksi
c------------
      open(13,file='ksi_gccr.txt')
      read(13,*) nl0,nnu0
      if (nl0.ne.nlve) then
         print*,'Probleme de dimension entre ksi.txt et lw'
         print*,'N levels = ',nl0,nlve
         stop
      endif
      if (nnu0.ne.nnuve) then
         print*,'Probleme de dimension entre ksi.txt et lw'
         print*,'N freq = ',nnu0,nnuve
         stop
      endif
      do band=1,nnuve
         read(13,*) lambdamin,lambdamax                ! en microns
         lambda(band)=(lambdamin+lambdamax)/2.*1.e-6   ! en m
         dlambda     =(1./lambdamin-1./lambdamax)*1.e4 ! en cm-1
c        print*,band,lambdamin,dlambda,lambdamax
         do i=0,nlve+1
           read(13,'(52e17.9)') (ksive(i,j,band),j=0,nlve+1)
c ecart-type MC sur les ksi: pas utilise
c          read(13,'(52e17.9)') (psive(i,j,band),j=0,nlve+1) 
c changement de convention (signe) pour ksi,
c et prise en compte de la largeur de bande (en cm-1):
           do j=0,nlve+1
              ksive(i,j,band) = -ksive(i,j,band)*dlambda
           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
      close(13)

      endif  ! firstcall

c --------------------------------------
c Calculation of the Psi matrix
c --------------------------------------

c temperature in the table layers
c -------------------------------
      
      do j=1,nlve
        ztemp(j) = PT(j)
      enddo

      ztemp(0) = PTSURF

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

      do band=1,nnuve
        y(0,band) = exp(bl(band)/ztemp(0))-1.
        bplck(0,band) = al(band)/(y(0,band))
       do j=1,nlve
c Developpement en polynomes ?
c       bplck(j,band) =            xp(1,band)
c    .                  +ztemp(j)*(xp(2,band)
c    .                  +ztemp(j)*(xp(3,band)
c    .                  +ztemp(j)*(xp(4,band)
c    .                  +ztemp(j)*(xp(5,band)
c    .                  +ztemp(j)*(xp(6,band) )))))

c B(T,l) = al/(exp(bl/T)-1)
        y(j,band) = exp(bl(band)/ztemp(j))-1.
        bplck(j,band) = al(band)/(y(j,band))
        zdblay(j,band) = al(band)*bl(band)*exp(bl(band)/ztemp(j))/
     .                  ((ztemp(j)**2)*(y(j,band)**2))
       enddo
       bplck(nlve+1,band) = 0.0
      enddo

c Calculation of Psi
c ------------------

      do band=1,nnuve
       do j=0,nlve+1
        do i=0,nlve+1
         psive(i,j,band)=ksive(i,j,band)*(bplck(i,band)-bplck(j,band))
        enddo
       enddo
      enddo

      do j=0,nlve+1
       do i=0,nlve+1
        psi_1(i,j) = 0.0  ! positif quand nrj de i->j 
       enddo
      enddo
      
      do band=1,nnuve
       do j=0,nlve+1
        do i=0,nlve+1
         psi_1(i,j) = psi_1(i,j)+psive(i,j,band)
        enddo
       enddo
      enddo

c Verif...-----------------------
c     open(11,file="psi.dat")
c     do i=0,nlve+1
c       write(11,'(I3,83E17.9)') i,(psi_1(j,i),j=0,nlve+1)
c     enddo
c     close(11)
c     stop
c -------------------------------

c --------------------------
c Calculation of the fluxes
c --------------------------

c flux aux intercouches:
c zlnet(i+1) est le flux net traversant le plafond de la couche i (+ vers le haut)
      do p=0,nlve ! numero de la couche
        zlnet(p+1) = 0.0
        do j=p+1,nlve+1
         do i=0,p
           zlnet(p+1) = zlnet(p+1)+psi_1(i,j)
         enddo
        enddo
      enddo

c flux net au sol, + vers le haut:
      PSOLLW = zlnet(1)
c flux vers le bas au sol, + vers le bas:
      PSOLLWDN = 0.0
      do i=1,nlve+1
        PSOLLWDN = PSOLLWDN+max(psi_1(i,0),0.0)
      enddo

c dfluxnet = radiative budget (W m-2)
      do p=0,nlve ! numero de la couche
        dzlnet(p) = 0.0
        do j=0,nlve+1
           dzlnet(p) = dzlnet(p)+psi_1(p,j)
        enddo
      enddo

      
c --------------------------------------
c Interpolation in the GCM vertical grid
c --------------------------------------

c Flux net
c --------
      
      do j=1,kflev+1
        ZFLNET(j) =  zlnet(j)
      enddo

      PTOPLW   = ZFLNET(kflev+1)
      
c Heating rates
c -------------

c  cool (K/s) = dfluxnet (W/m2)    ! positif quand nrj sort de la couche
c              *g        (m/s2)
c              /(-dp)  (epaisseur couche, en Pa=kg/m/s2)
c              /cp  (J/kg/K) 
      

      do j=1,kflev
        pdp(j)=(PPB(j)-PPB(j+1))*1.e5
      enddo

c calcul direct OU calcul par schema implicit
      if (1.eq.0) then 
        do j=1,kflev
         PCOOL(j) = dzlnet(j) *RG/RCPD / pdp(j)
        enddo
      else
        call lwi(nlve,nnuve,dzlnet,zdblay,pdp,ksive,PCOOL)
      endif
c     print*,dzlnet
c     print*,pdp
c     print*,PCOOL

      do j=1,kflev
        PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday
      enddo

      firstcall = .false.
      return
      end