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=81)    ! fichiers Vincent
      parameter (nnuve=68)   ! fichiers Vincent et Bullock
      real    dureejour
      parameter (dureejour=10.087e6)
      
      integer i,j,p,nl0,nnu0,band,k,l
      real   presve(nlve+1)   ! pressure levels in table (Pa->bar) 
      real   tempve(nlve+1)   ! temperature in table (K) (middle of layer)
      real   altve(nlve+1)    ! altitude in table (km)
      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)
      real   radbudget(kflev) ! Radiative budget on GCM grid
      real   coolrate(nlve)   ! cooling rates (K/s) on table grid
      character*22 nullchar
      real   lambdamin,lambdamax ! in microns
      real   dlambda             ! cm-1  

      real   y(0:nlve,nnuve)    ! intermediaire Planck
      real   pdp(kflev)         ! delta pression (Pa), grille GCM
      real   pdpve(nlve)        ! delta pression (Pa), grille table
      real   zdblay(nlve,nnuve) ! gradient en temperature de planck
      
      real   factflux
      real   facttemp,prT(kflev),prTve(nlve)

      logical firstcall
      data    firstcall/.true./

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

      if (firstcall) then

      print*,"PREMIER APPEL RADIATIF"

c Grilles alt et press
c---------------------
      open(11,file='mesh.txt')
      read(11,*) nl0,nnu0,i
      read(11,*) nullchar
      read(11,'(82(2x,F15.9))') altve
      read(11,'(82(2x,F15.9))') presve
      read(11,'(81(2x,F15.9))') tempve
      tempve(nlve+1)=tempve(nlve)
      close(11)
      if (nl0.ne.nlve) then
         print*,'Probleme de dimension entre mesh.txt et lw'
         print*,'N levels = ',nl0,nlve
         stop
      endif
      if (nnu0.ne.nnuve) then
         print*,'Probleme de dimension entre mesh.txt et lw'
         print*,'N freq = ',nnu0,nnuve
         stop
      endif
      do i=1,nlve+1
        presve(i) = presve(i)*1.e-5     ! convert to bar
      enddo

c Verifs...
c     print*, altve
c     print*, presve
      
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,'(83e17.9)') (ksive(i,j,band),j=0,nlve+1)
c ecart-type MC sur les ksi: pas utilise
c          read(13,'(83e17.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 i=1,kflev
           prT(i) = (PPB(i)+PPB(i+1))/2.
c          prT(i) = 10.**((log10(PPB(i))+log10(PPB(i+1)))/2.)
      enddo

      do j=1,nlve
           prTve(j) = (presve(j)+presve(j+1))/2.
c       prTve(j) = max(10.**((log10(presve(j))+log10(presve(j+1)))/2.)
c     .            ,1.e-5)
      enddo
      
      do j=1,nlve
        nl0 = 2
        do i=1,kflev-1
           if (prT(i).ge.prTve(j)) then
                nl0 = i+1
           endif
        enddo
        
        facttemp = (log10(prTve(j))-log10(prT(nl0-1)))
     .            /(log10(prT(nl0))-log10(prT(nl0-1)))
        ztemp(j) =   facttemp   *PT(nl0)
     .           + (1.-facttemp)*PT(nl0-1)

c       write(100,*) prTve(j),ztemp(j)
      enddo

      ztemp(0) = PTSURF

c     do j=1,kflev
c       write(101,*) prT(j),PT(j)
c     enddo

c     print*,'VERIF TEMP'
c     print*,PTSURF,PT
c     print*,ztemp
c     print*,tempve

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     do p=1,nlve
c       write(102,*) presve(p),zlnet(p),
c    .     (zlnet(p+1)-zlnet(p))/(presve(p)-presve(p+1))
c     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
        nl0 = 2
        do i=1,nlve
           if (presve(i).ge.PPB(j)) then
                nl0 = i+1
           endif
        enddo
        
        factflux = (log10(max(PPB(j),presve(nlve+1)))
     .                          -log10(presve(nl0-1)))
     .            /(log10(presve(nl0))-log10(presve(nl0-1)))
          ZFLNET(j) =  factflux   *zlnet(nl0)
     .             + (1.-factflux)*zlnet(nl0-1)
        
      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) 

c layers thickness on each pressure grid (in Pa)

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

      do j=1,nlve
        pdpve(j)=(presve(j)-presve(j+1))*1.e5
      enddo

c CHOIX CALCUL DIRECT OU IMPLICIT

c Ici, le budget radiatif est en calcul direct. 
c On ne fait rien. Si on veut l'implicit, on autorise le test suivant:

      if (1.eq.0) then 

c Pour calcul par schema implicite, on obtient en sortie de lwi le coolrate.
c Donc on actualise le dzlnet par dzlnet=coolrate*(cp/g)*pdpve

        call lwi(nlve,nnuve,dzlnet,zdblay,pdpve,ksive,coolrate)
	do j=1,nlve
	   dzlnet(j) = coolrate(j) *RCPD/RG *pdpve(j)
	enddo

      endif

c Interpolation on GCM grid of radiative budgets (dzlnet)

c on divise l'energie deposee dans la couche par l'epaisseur
c on moyenne ensuite ces valeurs (creneaux sur grille VE) 
c entre les niveaux de la grille GCM, et on multiplie ensuite par
c l'epaisseur (nouvelle grille) pour avoir l'energie deposee dans les
c couches GCM.

      i=1
      do j=1,kflev
        if (PPB(j+1).ge.presve(i+1)) then
          radbudget(j) = dzlnet(i)/(log10(presve(i+1))-log10(presve(i)))
     .                            *(log10(PPB(j+1))-log10(PPB(j)))
	else
          l=i+1
	  do while ((PPB(j+1).lt.presve(l+1)).and.(l.ne.nlve))
	     l=l+1
	  enddo
	radbudget(j) = dzlnet(i)/(log10(presve(i+1))-log10(presve(i)))*
     .                    (log10(presve(i+1))-log10(PPB(j)))
     .                +dzlnet(l)/(log10(presve(l+1))-log10(presve(l)))*
     .                    (log10(PPB(j+1))-log10(presve(l)))
	  do k=i+2,l
	radbudget(j) = radbudget(j)+dzlnet(k-1)
	  enddo
          i=l
        endif
      enddo

c     do i=1,kflev
c      print*,radbudget(i),prT(i)
c     enddo
c     do i=1,nlve
c      print*,dzlnet(i),prTve(i)
c     enddo
c     stop

c On obtient le coolrate en calculant: PCOOL = radbudget*(g/cp)/pdp

      do j=1,kflev
         PCOOL(j) = radbudget(j) *RG/RCPD / pdp(j)
         PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday
      enddo
      
c     print*,PCOOL

      firstcall = .false.
      return
      end