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 
      REAL coeftan
!-----------------------------------------------------------------------

      
      !alpha=(/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.e-13,1.e-13,1.e-13,1.e-13,1.e-13,1.e-13,1.e-13/)
!     alpha=(/0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.e-14,8.e-14,5.e-13,7.e-13,7.e-13,7.e-13,7.e-13,7.e-13,7.e-13,7.e-13,7.e-13/)

      lonw = 14.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.70    ! width of transition to alpha_top (Pa)  
      coeftan= 5.      

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(coeftan*(pplay(ig,l)-pref)/deltap))
          dtlw(ig,l)=-transition*alpha_top*BB(lonw,pt(ig,l))
         ENDDO
c        write(*,*) pplay(1,l),transition,dtlw(1,l)
      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