SUBROUTINE nirco2abs(nlon,nlev,nplay,dist_sol,nq,pq,
     $     mu0,fract,pdtnirco2)

       use dimphy
       use geometry_mod, only: longitude_deg, latitude_deg
       use chemparam_mod, only: i_co2, i_o
c       use compo_hedin83_mod2 


       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing heating rate due to
c   absorption by CO2 in the near-infrared
c   This version includes NLTE effects
c
c   (Scheme to be described in Forget et al., JGR, 2003)
c   (old Scheme described in Forget et al., JGR, 1999)
c
c   This version updated with a new functional fit,
c   see NLTE correction-factor of Lopez-Valverde et al (1998)
c   Stephen Lewis 2000
c
c   apr 2019 d.quirino   Improving NLTE params, SOIR/SPICAV Temp comparison 
c   oct 2014 g.gilli     Coupling with photochemical model
C   jan 2014 g.gilli     Revision (following martian non-lte param)    
C   jun 2013 l.salmi     First adaptation to Venus and NIR NLTE param
c   jul 2011 malv+fgg    New corrections for NLTE implemented
c   08/2002 : correction for bug when running with diurnal=F
c
c   author:  Frederic Hourdin 1996 
c   ------
c            Francois Forget 1999 
c
c   input:
c   ----- 
c   nlon                 number of gridpoint of horizontal grid
c   nlev                Number of layer
c   dist_sol              sun-Venus distance (AU)
c   mu0(nlon)          
c   fract(nlon)        day fraction of the time interval
c   declin                latitude of subslar point
c
c   output:
c   -------
c
c   pdtnirco2(nlon,nlev)      Heating rate (K/sec)
c
c
c=======================================================================
c
c    0.  Declarations :
c    ------------------
c

#include "YOMCST.h"
#include "clesphys.h"
c#include "comdiurn.h"
#include "nirdata.h"
c#include "tracer.h"
#include "mmol.h"
c-----------------------------------------------------------------------
c    Input/Output
c    ------------
      integer,intent(in) :: nlon ! number of (horizontal) grid points
      integer,intent(in) :: nlev ! number of atmospheric layers

      real,intent(in) :: nplay(nlon,nlev) ! Pressure
      real,intent(in) :: dist_sol ! Sun-Venus distance (in AU)
      integer,intent(in) :: nq ! number of tracers
      real,intent(in) :: pq(nlon,nlev,nq) ! mass mixing ratio tracers
      real,intent(in) :: mu0(nlon) ! solar angle
      real,intent(in) :: fract(nlon) ! day fraction of the time interval
c      real,intent(in) :: declin ! latitude of sub-solar point
      real :: co2vmr_gcm(nlon,nlev), o3pvmr_gcm(nlon,nlev)
 
      real,intent(out) :: pdtnirco2(nlon,nlev) ! heating rate (K/sec) 

c
c    Local variables :
c    -----------------
      INTEGER l,ig, n, nstep,i
      REAL co2heat0, zmu(nlon)

c     special diurnal=F 
      real mu0_int(nlon),fract_int(nlon),zday_int
      real ztim1,ztim2,ztim3,step

      logical onepeak
      parameter (onepeak=.false.)
c     parameter (onepeak=.true.)
c
c   local saved variables
c   ---------------------
      logical,save :: firstcall=.true.
      integer,save :: ico2=0 ! index of "co2" tracer
      integer,save :: io=0 ! index of "o" tracer

ccc=================================================
cccc     parameters for CO2 heating fit
ccc=================================================

c--------------------------------------------------
c One-peak martian-type fit => Gabriella (2014+)
c--------------------------------------------------
c     n_a  =  heating rate for Venusian day at p0, r0, mu =0 [K day-1]
c     Here p0 = p_cloud top [Pa]
c     n_p0 = is a pressure below which non LTE effects are significant [Pa]
c     n_a Solar heating [K/Eday] at the cloud top, taken from Crisps table     
 
      real n_a, n_p0, n_b, p_ctop

cc "Nominal" values used in Gilli+2017
c       parameter (n_a = 18.13/86400.0)     !c     K/Eday  ---> K/sec   
c       parameter (p_ctop=13.2e2)
c       parameter (n_p0=0.008) 

cc "New" values used to improve SPICAV/SOIR Temperature comparision (D.Quirino)
cc Gilli+2021
       parameter (n_a = 15.92/86400.0)     !c     K/Eday  ---> K/sec   
       parameter (p_ctop=19.85e2)
       parameter (n_p0=0.1)  
       parameter (n_b=1.362)

c    -- NLTE Param v2  --
C       parameter (n_p0=0.01)  
c       parameter (n_b = 1.3)

c--------------------------------------------------
c Multi-peaks Roldan-type fit => Laura (2013)
c New paramaters (Param9*0.5) => Enora (2021)
c--------------------------------------------------
c ENORA FINE TUNING used for VCD 1.1 
c  (fit to fig 12 Roldan-2000)
      real n_coFB, n_aFB, n_bFB, n_p0FB, n_eFB
      real n_coISO, n_aISO, n_bISO, n_p0ISO, n_eISO
      real n_coFH, n_aFH, n_bFH, n_p0FH, n_eFH
      real n_co43, n_a43, n_b43, n_p043, n_e43
      real n_co43b, n_a43b, n_b43b, n_p043b, n_e43b
      real n_conir, n_anir, n_bnir, n_p0nir, n_enir

      parameter (n_coFB=119./86400.0)      !c     K/Eday ---> K/sec
      parameter (n_aFB=0.185)
      parameter (n_bFB=3.7)
      parameter (n_p0FB=2.9e-4)
      parameter (n_eFB=0.76)

      parameter (n_coISO=265./86400.0)      !c     K/Eday ---> K/sec
      parameter (n_aISO=0.313)
      parameter (n_bISO=1.65)
      parameter (n_p0ISO=0.076)
      parameter (n_eISO=0.99)

      parameter (n_coFH=2.5/86400.0)      !c     K/Eday ---> K/sec
      parameter (n_aFH=3.98)
      parameter (n_bFH=2.9)
      parameter (n_p0FH=0.17)
      parameter (n_eFH=2.16)

      parameter (n_co43=55./86400.0)      !c     K/Eday ---> K/sec
      parameter (n_a43=0.625)
      parameter (n_b43=2.6)
      parameter (n_p043=0.043)
      parameter (n_e43=1.654)

!     parameter (n_co43b=100./86400.0)      !c     K/Eday ---> K/sec
! => fine tuning: not affected by the *0.5 below (see ENORA FINE TUNING)
      parameter (n_co43b=200./86400.0)      !c     K/Eday ---> K/sec
      parameter (n_a43b=5.5)
      parameter (n_b43b=2.3)
      parameter (n_p043b=1.)
      parameter (n_e43b=0.4)

      parameter (n_conir=6.5/86400.0)      !c     K/Eday ---> K/sec
      parameter (n_anir=35.65)
      parameter (n_bnir=2.1)
      parameter (n_p0nir=0.046)
      parameter (n_enir=0.9)

      real :: picFB(nlon,nlev), picISO(nlon,nlev), picFH(nlon,nlev)
      real :: pic43(nlon,nlev), pic43b(nlon,nlev), picnir(nlon,nlev)

ccc=================================================

c     Variables added to implement NLTE correction factor (feb 2011)
      real    pyy(nlev)
      real    cor1(nlev),oldoco2(nlev),alfa2(nlev)
      real    p2011,cociente1,merge
      real    cor0,oco2gcm

c----------------------------------------------------------------------  
c     Initialisation
c     --------------
      if (firstcall) then
        if (nircorr.eq.1) then
c          ! we will need co2 and o tracers
          ico2= i_co2
          if (ico2==0) then
            write(*,*) "nirco2abs error: I need a CO2 tracer"
            write(*,*) "     when running with nircorr==1"
           stop
          endif
          io=i_o
          if (io==0) then
            write(*,*) "nirco2abs error: I need an O tracer"
            write(*,*) "     when running with nircorr==1"
            stop
          endif
        endif
        firstcall=.false.
      endif
c     --------------
c     co2heat0 is correction for dist_sol (is 1 for Venus)
      co2heat0=(0.7233/dist_sol)**2

      pdtnirco2(:,:)=0.
c----------------------------------------------------------------------  

c      
c     Simple calcul for a given sun incident angle (if cycle_diurne=T)
c     --------------------------------------------

      IF (cycle_diurne) THEN  

         do ig=1,nlon    
            zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35.

c---------------------------
           if (onepeak) then
c---------------------------
            if(nircorr.eq.1) then
               do l=1,nlev
                  pyy(l)=nplay(ig,l)
               enddo
               call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
               call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
               call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
            endif
            do l=1,nlev
c           Calculations for the O/CO2 correction
               if(nircorr.eq.1) then
                  cor0=1./(1.+n_p0/nplay(ig,l))**n_b
                  if(pq(ig,l,ico2) .gt. 1.e-6) then
                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                     ! handle the rare cases when pq(ig,l,io)<0
                     if (pq(ig,l,io).lt.0) then
                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
     &                            " pq(ig,l,io)=",pq(ig,l,io)
                       oco2gcm=1.e6
                     endif
                  else
                     oco2gcm=1.e6
                  endif
                  cociente1=oco2gcm/oldoco2(l)
                  
c                  WRITE(*,*) "nirco2abs line 211", l, cociente1

                  merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
     $                 (1.-alfa2(l))
                  merge=10**merge
                  p2011=sqrt(merge)*cor0

               else if (nircorr.eq.0) then
                  p2011=1.
                  cor1(l)=1.
               endif

              if(fract(ig).gt.0.) pdtnirco2(ig,l)=
     &             co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
     &             /(1.+n_p0/nplay(ig,l))**n_b
c           Corrections from tabulation
     $              * cor1(l) * p2011
              
            enddo !nlev
c---------------------------
           else  ! multipeak
c---------------------------
            do l=1,nlev
               if(fract(ig).gt.0.) then 
                   picFB(ig,l)=n_coFB
     &              *((n_aFB/nplay(ig,l))**n_eFB)
     &		    *zmu(ig)**0.82
     &             /(1.+n_p0FB/nplay(ig,l))**n_bFB

                   picISO(ig,l)=n_coISO
     &              *((n_aISO/nplay(ig,l))**n_eISO)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0ISO/nplay(ig,l))**n_bISO

                   picFH(ig,l)=n_coFH
     &              *((n_aFH/nplay(ig,l))**n_eFH)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0FH/nplay(ig,l))**n_bFH

                   pic43(ig,l)=n_co43
     &              *((n_a43/nplay(ig,l))**n_e43)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p043/nplay(ig,l))**n_b43

                   pic43b(ig,l)=n_co43b
     &              *((n_a43b/nplay(ig,l))**n_e43b)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p043b/nplay(ig,l))**n_b43b

                   picnir(ig,l)=n_conir
     &              *((n_anir/nplay(ig,l))**n_enir)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0nir/nplay(ig,l))**n_bnir

	           pdtnirco2(ig,l)=co2heat0*
     &           (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
     &           +pic43b(ig,l)+picnir(ig,l))*0.5  ! *0.5 = ENORA FINE TUNING
		  
               endif
            enddo !nlev
c---------------------------
           endif
c---------------------------
         enddo  !nlon


c     Averaging over diurnal cycle (if diurnal=F)
c     -------------------------------------------
c     NIR CO2 abs is slightly non linear. To remove the diurnal
c     cycle, it is better to average the heating rate over 1 day rather
c     than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)

      ELSE      ! if (.not.diurnal) then
         nstep = 20    ! number of integration step /sol
         do n=1,nstep

            zday_int = (n-1)/float(nstep)

            CALL zenang(0.,zday_int,RDAY/nstep,
     &                  latitude_deg,longitude_deg,
     &                  mu0_int,fract_int)

          do ig=1,nlon
               zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35.

c---------------------------
           if (onepeak) then
c---------------------------
            if(nircorr.eq.1) then
               do l=1,nlev
                  pyy(l)=nplay(ig,l)
               enddo
               call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
               call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
               call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
            endif
            do l=1,nlev
c           Calculations for the O/CO2 correction
               if(nircorr.eq.1) then
                  cor0=1./(1.+n_p0/nplay(ig,l))**n_b
                  if(pq(ig,l,ico2) .gt. 1.e-6) then
                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                     ! handle the rare cases when pq(ig,l,io)<0
                     if (pq(ig,l,io).lt.0) then
                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
     &                            " pq(ig,l,io)=",pq(ig,l,io)
                       oco2gcm=1.e6
                     endif
                  else
                     oco2gcm=1.e6
                  endif
                  cociente1=oco2gcm/oldoco2(l)
                  
c                  WRITE(*,*) "nirco2abs line 211", l, cociente1

                  merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
     $                 (1.-alfa2(l))
                  merge=10**merge
                  p2011=sqrt(merge)*cor0

               else if (nircorr.eq.0) then
                  p2011=1.
                  cor1(l)=1.
               endif

              if(fract(ig).gt.0.) pdtnirco2(ig,l)=
     &              pdtnirco2(ig,l) + (1/float(nstep))*
     &             co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
     &             /(1.+n_p0/nplay(ig,l))**n_b
c           Corrections from tabulation
     $              * cor1(l) * p2011 
              
            enddo !nlev
c---------------------------
           else  ! multipeak
c---------------------------
            do l=1,nlev
               if(fract(ig).gt.0.) then 
                   picFB(ig,l)=n_coFB
     &              *((n_aFB/nplay(ig,l))**n_eFB)
     &		    *zmu(ig)**0.82
     &             /(1.+n_p0FB/nplay(ig,l))**n_bFB

                   picISO(ig,l)=n_coISO
     &              *((n_aISO/nplay(ig,l))**n_eISO)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0ISO/nplay(ig,l))**n_bISO

                   picFH(ig,l)=n_coFH
     &              *((n_aFH/nplay(ig,l))**n_eFH)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0FH/nplay(ig,l))**n_bFH

                   pic43(ig,l)=n_co43
     &              *((n_a43/nplay(ig,l))**n_e43)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p043/nplay(ig,l))**n_b43

                   pic43b(ig,l)=n_co43b
     &              *((n_a43b/nplay(ig,l))**n_e43b)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p043b/nplay(ig,l))**n_b43b

                   picnir(ig,l)=n_conir
     &              *((n_anir/nplay(ig,l))**n_enir)
     &		    *zmu(ig)**0.55
     &             /(1.+n_p0nir/nplay(ig,l))**n_bnir

	           pdtnirco2(ig,l)=
     &               pdtnirco2(ig,l)+(1/float(nstep))*co2heat0*
     &           (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
     &           +pic43b(ig,l)+picnir(ig,l))*0.5  ! *0.5 = ENORA FINE TUNING
		  
               endif
            enddo !nlev
c---------------------------
           endif
c---------------------------
          enddo  !nlon
         enddo  !nstep
      
      END IF  ! diurnal cycle

      return
      end 

      
      subroutine interpnir(escout,p,nlev,escin,pin,nl)
C
C subroutine to perform linear interpolation in pressure from 1D profile 
C escin(nl) sampled on pressure grid pin(nl) to profile
C escout(nlev) on pressure grid p(nlev).
C
      real escout(nlev),p(nlev)
      real escin(nl),pin(nl),wm,wp
      integer nl,nlev,n1,n,nm,np
      do n1=1,nlev
         if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
c            escout(n1) = 0.0
            escout(n1) = 1.e-15
         else
            do n = 1,nl-1
               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
                  nm=n
                  np=n+1
                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
                  wp=1.0 - wm
               endif
            enddo
            escout(n1) = escin(nm)*wm + escin(np)*wp
         endif
      enddo
      return
      end