subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,eps23,eps33,eps76)

      use comgeomfi_h
      use callkeys_mod, only: strobel
      use tracer_h, only: igcm_n2,mmol
      use comcstfi_mod, only: pi, g, cpp,r
      implicit none

!==================================================================
!     Purpose
!     -------
!     Calculation of epsilon - non LTE efficiency factor
!     For CH4 heating and cooling rate calculation
!     Cf Strobel et al 1996
!
!     Inputs
!     ------
!     ngrid                 Number of vertical columns
!     nlayer                Number of layers
!     pplay(ngrid,nlayer)   Pressure layers
!     pt
!
!     Outputs
!     -------
!
!     eps_sw       ! efficiency factor for heating rate
!     eps_lw       ! efficiency factor for cooling rate
!
!     Authors
!     -------
!     Tanguy Bertrand (2015,2022)
!
!==================================================================

!-----------------------------------------------------------------------

!     Arguments

      INTEGER ngrid, nlayer, nq
      REAL pplay(ngrid,nlayer)
      REAL pplev(ngrid,nlayer)
      REAL pt(ngrid,nlayer)
      REAL vmrch4(ngrid,nlayer)

!-----------------------------------------------------------------------
!     Local variables

      INTEGER l,ig

      REAL cbol         ! Boltzman constant
      REAL avog         ! avogadro constant
      REAL beta23
      REAL beta33
      REAL nd(ngrid,nlayer)      ! density number  cm-3
      REAL na(ngrid,nlayer)      ! column density number  cm-2
      REAL rho(ngrid,nlayer)     ! density of atmosphere
      REAL kn3,k3,k2n4,k3n4,k4
      REAL An3,A3,A2n4,A3n4,A4
      REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4
      REAL phin3(ngrid,nlayer)
      REAL phi3(ngrid,nlayer)
      REAL phi2n4(ngrid,nlayer)
      REAL phi3n4(ngrid,nlayer)
      REAL phi4(ngrid,nlayer)
      REAL phi4T(ngrid,nlayer)

!     DATA from Strobel et al 1996
      REAL, DIMENSION(22) :: gam2       ! band escape fonction 2.3
      REAL, DIMENSION(22) :: Naval2     ! column density number 2.3
      REAL, DIMENSION(22) :: gam3       ! band escape fonction 3.3
      REAL, DIMENSION(22) :: Naval3    ! column density number 3.3
      REAL gammafi23(nlayer),gammafi33(nlayer)       ! Band espace fonction
      REAL gammafi23mod(nlayer),gammafi33mod(nlayer)       ! Band espace fonction
      REAL vmr
      REAL k4zT  ! function to calculate k4 depending on T and vmr
      REAL log_square ! log function
      REAL log_affine ! log function
      REAL ratio_affine ! linear approx function to calculate ratio for gamma calculation
      REAL ratio_square !
      REAL ratio_sq,ratio_af

!     Output
      REAL,intent(out) :: eps23(ngrid,nlayer)
      REAL,intent(out) :: eps33(ngrid,nlayer)
      REAL,intent(out) :: eps76(ngrid,nlayer)

!     Option simplified coefficient
      REAL klw, ksw,  eps_sw1Pa

!-----------------------------------------------------------------------
      eps23(:,:)=1.
      eps33(:,:)=1.
      eps76(:,:)=1.

      IF (strobel) then

!       Constantes
        cbol=1.3806488e-23 !Boltzman constant
        avog=6.022141e23

        ! band 3
        k3=2.8e-11    !cm-3
        A3=25.2
        dphi3=k3/A3
        ! band n3
        kn3=2.8e-11    !cm-3
        An3=27.3
        dphin3=kn3/An3
        ! band 3n4
        k3n4=1.16e-11  !cm-3
        A3n4=2.12
        dphi3n4=k3n4/(3.*A3n4)
        ! band 4
        !! Here we use function k4zT instead of constant k4
        k4=4.5e-16        !cm-3
        A4=2.12
        dphi4=k4/A4
        ! band 2n4
        k2n4=1.16e-11
        A2n4=2.12
        dphi2n4=k2n4/(2.*A2n4)

!       2.3 band :
        beta23=0.09
!       3.3 band :
        beta33=0.13

!       Escape Function Gamma (Data Strobel et al) :
        ! 100 K 2.3
        Naval2=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
                7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
                9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
                1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
        gam2=(/1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,9.90000000e-01, &
              9.80000000e-01,9.50000000e-0,9.00000000e-01,8.00000000e-01,7.00000000e-01,5.00000000e-01, &
              3.50000000e-01,2.20000000e-01,1.20000000e-01,6.50000000e-02,3.20000000e-02,1.40000000e-02, &
              6.20000000e-03,2.50000000e-03,1.05000000e-03,4.83227951e-06/)

        ! 100 K 3.3
        Naval3=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
                7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
                9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
                1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
        gam3=(/1.00000000e+00,9.80000000e-01,9.00000000e-01,8.50000000e-01,7.00000000e-01,5.00000000e-01, &
              3.00000000e-01,1.50000000e-01,1.00000000e-01,6.00000000e-02,4.00000000e-02,2.30000000e-02, &
              1.45000000e-02,9.00000000e-03,5.00000000e-03,3.00000000e-03,1.50000000e-03,8.00000000e-04, &
              3.80000000e-04,1.80000000e-04,8.30000000e-05,6.81737544e-07/)

        DO ig=1,ngrid

          phi3n4(ig,:)=1.
          phi2n4(ig,:)=1.
          phin3(ig,:)=1.
          phi4(ig,:)=1.
          nd(ig,:)=1.
          na(ig,:)=1.

          !write(*,*) 'TB22 pplev=',pplev(1,:)
          !write(*,*) 'TB22 tlev=',pt(1,:)

          DO l=1,nlayer
           !!calculation of density
           rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))

           !!calculation of column number density for ch4 (cm-2)
           na(ig,l)=vmrch4(ig,l)/100./(mmol(igcm_n2)*1.e-3)* &
                             pplev(ig,l)/g*avog*1.e-4

           !!calculation of atmospheric number density (cm-3) for each layer
           nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)*1.e-3)*avog*1.e-6

           !!calculation of the phi :
           phin3(ig,l)=dphin3*nd(ig,l)
           phi3(ig,l)=dphi3*nd(ig,l)
           phi2n4(ig,l)=dphi2n4*nd(ig,l)
           phi3n4(ig,l)=dphi3n4*nd(ig,l)
           phi4(ig,l)=dphi4*nd(ig,l)
           !! Getting phi4 depending on T and vmr
           vmr=vmrch4(ig,l)/100.
           phi4T(ig,l)=k4zT(pt(ig,l),vmr)/A4*nd(ig,l)
          ENDDO
          !write(*,*) 'TB22 nd=',nd(1,:)
          !write(*,*) 'TB22 vmr=',vmrch4
          !write(*,*) 'TB22 phi4T=',phi4T

          !! interpolation of na values in cm-2 from Strobel values for
          !the 100K model:
          CALL interp_line(Naval2,gam2,22,na(ig,:),gammafi23,nlayer)
          CALL interp_line(Naval3,gam3,22,na(ig,:),gammafi33,nlayer)

          !write(*,*) 'TB22rad na = ',na(1,:)
          !write(*,*) 'TB22rad gammafi = ',gammafi23

          DO l=1,nlayer

           !! Getting ratio depending on temperature and n (cm-2)
           ratio_sq=ratio_square(pt(ig,l),na(ig,l))
           ratio_af=ratio_affine(pt(ig,l),na(ig,l))
           !! Modified gamma depending on T
           gammafi23mod(l)=min(ratio_sq*gammafi23(l),1.)
           gammafi33mod(l)=min(ratio_af*gammafi33(l),1.)

           !! Heating : cf Strobel 1996 and Zalucha 2011
           ! eps 2.3
           eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23mod(l))*     &
                      (beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)*  &
                      phi4T(ig,l)/(phi4T(ig,l)+gammafi23mod(l))*   &
                      (3.*phi2n4(ig,l)+1.)/(1.+phi2n4(ig,l))*(1-beta23)/3.)

           ! eps 3.3
           eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33mod(l))*     &
                      (beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)*  &
                      phi4T(ig,l)/(phi4T(ig,l)+gammafi33mod(l))*   &
                      (1-beta33))

           ! Cooling from Zalucha et al 2013 eq10
           eps76(ig,l)=1./(1.+gammafi23mod(l)/(2.*phi4T(ig,l)))

          ENDDO
          !write(*,*) 'TB22rad eps23 = ',eps23(1,:)
          !write(*,*) 'TB22rad eps33 = ',eps33(1,:)
          !write(*,*) 'TB22rad gamma33 = ',gammafi33mod
          !write(*,*) 'TB22rad eps76 = ',eps76(1,:)

        ENDDO   ! ngrid

      ELSE     !    if (strobel)
         !if strobel = false

         klw = 3.6E-4  ! fit strobel
         !eps_sw1Pa = 0.6 ! fit strobel
         eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5
         !eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25

         ksw = 8.5E-5*(0.9/(eps_sw1Pa-0.1)-1.)
         do l=1, nlayer
           do ig=1, ngrid
              rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
              eps23(ig,l) =  0.1+ 0.9/(1.+ksw/rho(ig,l))
              eps76(ig,l) = 1./(1.+klw/rho(ig,l))
           end do
         end do

      END IF

      end subroutine nlte_ch4

      real function k4zT(temp,vmr)
      implicit none
      real, intent(in) ::  temp,vmr
      k4zT=1.9e-15*exp((temp-240.)/70.)*(1+10.*vmr)
      return
      end function k4zT

      !! Log log affine and square function for ratio adjustment and extrapolation

      real function log_square(x,y1,y2,x1,x2,coef)
      implicit none
      real, intent(in) ::  x,y1,y2,x1,x2,coef
      real aa,bb
      aa=(log(y2)-log(y1))/((log(coef*x2))**4-(log(coef*x1)**4))
      bb=log(y2)-aa*(log(coef*x2))**4
      log_square=exp(aa*(log(coef*x))**4+bb)
      return
      end function log_square

      real function log_affine(x,y1,y2,x1,x2)
      real, intent(in) :: x,y1,y2,x1,x2
      real aa,bb
      aa=(log(y2)-log(y1))/(log(x2)-log(x1))
      bb=log(y2)-aa*log(x2)
      log_affine=exp(aa*log(x)+bb)
      return
      end function log_affine

      !! Ratio of gamma depending on T. Linear approximation
      real function ratio_affine(temp,n)
      implicit none
      real, intent(in) :: temp,n
      real aa,bb,ratio0
      real log_affine
      aa=(1.-0.2)/(100.-40.)
      bb=1.-aa*100.
      ratio0=aa*temp+bb
      ratio_affine=log_affine(n,1.,ratio0,1e15,1e22)
      return
      end function ratio_affine

      !! Ratio of gamma depending on T. Linear approximation
      real function ratio_square(temp,n)
      implicit none
      real, intent(in) :: temp,n
      real aa,bb,ratio0,coef
      real log_square
      coef=1.e-12
      aa=(1.-0.2)/(100.-40.)
      bb=1.-aa*100.
      ratio0=aa*temp+bb
      ratio_square=log_square(n,1.,ratio0,1e15,1e22,coef)
      return
      end function ratio_square