SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,
     $             pp,ps_rad,ptsurf,pt,
     $             pfluxir,pdtlw,
     $             lwrite)

      IMPLICIT NONE
c=======================================================================
c
c   calcul de l'evolution de la temperature sous l'effet du rayonnement
c   infra-rouge.
c   Pour simplifier, les transmissions sont precalculees et ne
c   dependent que de l'altitude.
c
c   arguments:
c   ----------
c
c   entree:
c   -------
c      ngrid             nombres de points de la grille horizontale
c      nlayer              nombre de couches
c      ptsurf(ngrid)     temperature de la surface
c      pt(ngrid,nlayer)    temperature des couches
c      pp(ngrid,nlayer+1)  pression entre les couches
c      lwrite            variable logique pour sorties
c
c   sortie:
c   -------
c      pdtlw(ngrid,nlayer) taux de refroidissement
c      pfluxir(ngrid)    flux infrarouge sur le sol
c
c=======================================================================

c   declarations:
c   -------------

#include "comcstfi.h"

c   arguments:
c   ----------

      INTEGER ngrid,nlayer
      REAL coefir,emissiv(ngrid),ps_rad
      REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
      REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)
      LOGICAL lwrite

c   variables locales:
c   ------------------

      INTEGER nlevel,ilev,ig,i,il
      REAL zplanck(ngrid,nlayer+1),zcoef
      REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1)
      REAL zflux(ngrid,nlayer+1)
      REAL zlwtr1(ngrid),zlwtr2(ngrid)
      REAL zup(ngrid,nlayer+1),zdup(ngrid)
      REAL stephan

      LOGICAL lstrong
      SAVE lstrong,stephan
      DATA lstrong/.true./
!$OMP THREADPRIVATE(lstrong,stephan)

c-----------------------------------------------------------------------
c   initialisations:
c   ----------------

      nlevel=nlayer+1
      stephan=5.67e-08

c-----------------------------------------------------------------------
c   2. calcul des quantites d'absorbants:
c   -------------------------------------

c   absorption forte
      IF(lstrong) THEN
         DO ilev=1,nlevel
            DO ig=1,ngrid
               zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
            ENDDO
         ENDDO
         IF(lwrite) THEN
            DO ilev=1,nlayer
             PRINT*,' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
            ENDDO
         ENDIF
         zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))

c   absorption faible
      ELSE
         DO ilev=1,nlevel
            DO ig=1,ngrid
               zup(ig,ilev)=pp(ig,ilev)
            ENDDO
         ENDDO
         zcoef=-log(coefir)/ps_rad
      ENDIF


c-----------------------------------------------------------------------
c   2. calcul de la fonction de corps noir:
c   ---------------------------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
            zplanck(ig,ilev)=stephan*
     $      zplanck(ig,ilev)*zplanck(ig,ilev)
		 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   4. flux descendants:
c   --------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zfluxdn(ig,ilev)=0.
		 ENDDO
         DO ig=1,ngrid
            zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
         ENDDO
         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)

         DO il=nlayer,ilev,-1
            zlwtr2(:)=zlwtr1(:)
            DO ig=1,ngrid
               zdup(ig)=zup(ig,ilev)-zup(ig,il)
            ENDDO
            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
            DO ig=1,ngrid
               zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+
     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
			ENDDO
		 ENDDO
      ENDDO

      DO ig=1,ngrid
         zfluxdn(ig,nlevel)=0.
         pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
      ENDDO

      DO ig=1,ngrid
         zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
         zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)
     $   +(1.-emissiv(ig))*zfluxdn(ig,1)
      ENDDO

c-----------------------------------------------------------------------
c   3. flux montants:
c   ------------------

      DO ilev=1,nlayer
         DO ig=1,ngrid
            zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
         ENDDO
         CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
         DO ig=1,ngrid
            zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
		 ENDDO
         DO il=1,ilev
            zlwtr2(:)=zlwtr1(:)
            DO ig=1,ngrid
               zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
            ENDDO
            CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
            DO ig=1,ngrid
               zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+
     $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
			ENDDO
         ENDDO

      ENDDO

c-----------------------------------------------------------------------
c   5. calcul des flux nets:
c   ------------------------

      DO ilev=1,nlevel
         DO ig=1,ngrid
            zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
		 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   6. Calcul des taux de refroidissement:
c   --------------------------------------
   
      DO ilev=1,nlayer
         DO ig=1,ngrid
            pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*
     $           g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
		 ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   10. sorties eventuelles:
c   ------------------------

      IF (lwrite) THEN

         PRINT*
         PRINT*,'Diagnostique rayonnement thermique'
         PRINT*
         PRINT*,'temperature     ',
     $   'flux montant    flux desc.     taux de refroid.'
         i=ngrid/2+1
         WRITE(6,9000) ptsurf(i)
         DO ilev=1,nlayer
            WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),
     $      zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
		 ENDDO
         WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)

      ENDIF

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

      RETURN
9000  FORMAT(4e18.4)
      END