source: trunk/LMDZ.MARS/libf/phymars/nirco2abs.F @ 1215

      SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq,
     $     mu0,fract,declin,pdtnirco2)
                                                   
       use tracer_mod, only: igcm_co2, igcm_o
       use comdiurn_h, only: sinlon, coslon, sinlat, coslat
       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   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   ngrid                 number of gridpoint of horizontal grid
c   nlayer                Number of layer
c   dist_sol              sun-Mars distance (AU)
c   mu0(ngrid)          
c   fract(ngrid)          day fraction of the time interval
c   declin                latitude of subslar point
c
c   output:
c   -------
c
c   pdtnirco2(ngrid,nlayer)      Heating rate (K/s)
c
c
c=======================================================================
c
c    0.  Declarations :
c    ------------------
c
!#include "dimensions.h"
!#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
!#include "comdiurn.h"
#include "nirdata.h"
!#include "tracer.h"

c-----------------------------------------------------------------------
c    Input/Output
c    ------------
      integer,intent(in) :: ngrid ! number of (horizontal) grid points
      integer,intent(in) :: nlayer ! number of atmospheric layers
      real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
      real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
      integer,intent(in) :: nq ! number of tracers
      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
      real,intent(in) :: mu0(ngrid) ! solar angle
      real,intent(in) :: fract(ngrid) ! day fraction of the time interval
      real,intent(in) :: declin ! latitude of sub-solar point
      
      real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
c
c    Local variables :
c    -----------------
      INTEGER l,ig, n, nstep,i
      REAL co2heat0, zmu(ngrid)

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

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
c     p0noonlte is a pressure below which non LTE effects are significant.
c     REAL p0nonlte
c     DATA p0nonlte/7.5e-3/
c     SAVE p0nonlte

c     parameters for CO2 heating fit
      real n_a, n_p0, n_b
      parameter (n_a=1.1956475)
      parameter (n_b=1.9628251)
      parameter (n_p0=0.0015888279)

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

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

c     Initialisation
c     --------------
      if (firstcall) then
        if (nircorr.eq.1) then
          ! we will need co2 and o tracers
          ico2=igcm_co2
          if (ico2==0) then
            write(*,*) "nirco2abs error: I need a CO2 tracer"
            write(*,*) "     when running with nircorr==1"
            stop
          endif
          io=igcm_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     co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
      co2heat0=n_a*(1.52/dist_sol)**2/daysec

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

      IF (diurnal) THEN  
         do ig=1,ngrid
            zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35.

            if(nircorr.eq.1) then
               do l=1,nlayer
                  pyy(l)=pplay(ig,l)
               enddo

               call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
               call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
               call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
            endif

            do l=1,nlayer
!           Calculations for the O/CO2 correction
               if(nircorr.eq.1) then
                  cor0=1./(1.+n_p0/pplay(ig,l))**n_b
                  if(pq(ig,l,ico2).gt.1.e-6) then
                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                  else
                     oco2gcm=1.e6
                  endif
                  cociente1=oco2gcm/oldoco2(l)
                  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*sqrt((700.*zmu(ig))/pplay(ig,l))
     &             /(1.+n_p0/pplay(ig,l))**n_b
!           Corrections from tabulation
     $              * cor1(l) * p2011
c           OLD SCHEME (forget et al. 1999)
c    s           co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
c    s          / (1+p0nonlte/pplay(ig,l))
           enddo
         enddo
         

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)
            ztim2=COS(declin)*COS(2.*pi*(zday_int-.5))
            ztim3=-COS(declin)*SIN(2.*pi*(zday_int-.5))
            CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
     s             ztim1,ztim2,ztim3,
     s             mu0_int,fract_int)
            do ig=1,ngrid
               zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35.

               if(nircorr.eq.1) then
                  do l=1,nlayer
                     pyy(l)=pplay(ig,l)
                  enddo
               call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
               call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
               call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
               endif

               do l=1,nlayer
                  if(nircorr.eq.1) then
                     cor0=1./(1.+n_p0/pplay(ig,l))**n_b
                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                     cociente1=oco2gcm/oldoco2(l)
                     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_int(ig).gt.0.) pdtnirco2(ig,l)=
     &                 pdtnirco2(ig,l) + (1/float(nstep))*
     &                 co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
     &                 /(1.+n_p0/pplay(ig,l))**n_b
!                      Corrections from tabulation
     $                 * cor1(l) * p2011
               enddo
            enddo
         end do
      END IF

      return
      end 


      
      subroutine interpnir(escout,p,nlayer,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(nlayer) on pressure grid p(nlayer).
C
      real escout(nlayer),p(nlayer)
      real escin(nl),pin(nl),wm,wp
      integer nl,nlayer,n1,n,nm,np
      do n1=1,nlayer
         if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
            escout(n1) = 0.0
         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