subroutine cooling_hcn_c2h2(ngrid,nlayer,pplay,pt,dtlw)
          
      implicit none

!==================================================================
!     Purpose
!     -------
!     Calculation of cooling rate for C2H2-HCN
!     = f(pplay) * B(lambda,T)
!
!     Inputs
!     ------ 
!     ngrid                 Number of vertical columns
!     nlayer                Number of layers
!     pt
!
!     Outputs
!     -------
!
!     dtlw       ! cooling rate
!
!     Authors
!     ------- 
!     Tanguy Bertrand (2016)
!     FF (2016)
!==================================================================
!     Arguments
      INTEGER ngrid, nlayer
      REAL pplay(ngrid,nlayer)            ! pres. level in GCM mid of layer
      REAL pt(ngrid,nlayer)
      REAL dtlw(ngrid,nlayer)
!     Local variables
      INTEGER l,ig
      REAL lonw
      REAL alpha, alpha_top
      REAL pref, deltap 
      REAL transition
      REAL BB 
!-----------------------------------------------------------------------

      lonw = 14.e-6  ! 14um
      alpha_top=5.e-11  ! 1.e-13  ! cooling constant at top of atmosphere
      !pref = 0.12      ! pressure at mid transition fo alpha_top (Pa)
      pref = 0.02      ! pressure at mid transition fo alpha_top (Pa)
      deltap = 0.1    ! width of transition to alpha_top (Pa)  
      !lonw = 11.e-6  ! 14um
      !alpha_top=1.e-9  ! 1.e-13  ! cooling constant at top of atmosphere
      !pref = 0.06      ! pressure at mid transition fo alpha_top (Pa)
      !deltap = 0.1    ! width of transition to alpha_top (Pa)  

c     transition = 0 if p>pref+deltap/2 and 1 if p< pref-deltap/2
      DO l = 1, nlayer
         DO ig = 1, ngrid
          transition  =  0.5*(1-tanh((pplay(ig,l)-pref)/deltap))
          dtlw(ig,l)=-transition*alpha_top*BB(lonw,pt(ig,l))
         ENDDO
      ENDDO   
      end 

c******************************************************
c   FUNCTION Blackbody  (Planck)
c**********************************************************
      function BB (lw, T)
c Variable declaration
c --------------------
c wavelenght (m), Temperature (K)
      real lw,T
c constant
      real c1,c2
      parameter ( c1=1.19103E-16 )
      parameter (c2=1.43887E-2 )
c function
c---------
      BB= (c1/lw**5)/(-1.+exp(c2/(lw*T)))
      return
      end