subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
          tsurf,fract,dist_star,                               &
          dtlw,dtsw,fluxsurf_lw,                               &
          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
          fluxabs_sw,fluxtop_dn,                               &
          OLR_nu,OSR_nu,                                       &
          lastcall)

      use mod_phys_lmdz_para, only : is_master
      use radinc_h
      use radcommon_h
      use gases_h
      USE tracer_h
      use callkeys_mod, only: global1d, szangle
      use comcstfi_mod, only: pi, mugaz, cpp
      use callkeys_mod, only: diurnal,tracer,        &
                              strictboundcorrk,specOLR

      implicit none

!==================================================================
!
!     Purpose
!     -------
!     Solve the radiative transfer using the correlated-k method for
!     the gaseous absorption and the Toon et al. (1989) method for
!     scatttering due to aerosols.
!
!     Authors
!     ------- 
!     Emmanuel 01/2001, Forget 09/2001
!     Robin Wordsworth (2009)
!
!==================================================================

!-----------------------------------------------------------------------
!     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
!     Layer #1 is the layer near the ground. 
!     Layer #nlayer is the layer at the top.
!-----------------------------------------------------------------------


      ! INPUT
      INTEGER,INTENT(IN) :: ngrid                  ! Number of atmospheric columns.
      INTEGER,INTENT(IN) :: nlayer                 ! Number of atmospheric layers.
      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (X/m2).
      INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
      REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
      REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral Short Wavelengths Albedo. By MT2015
      REAL,INTENT(IN) :: emis(ngrid)               ! Long Wave emissivity.
      REAL,INTENT(IN) :: mu0(ngrid)                ! Cosine of sun incident angle.
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)     ! Inter-layer pressure (Pa).
      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
      REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
      REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
      logical,intent(in) :: lastcall               ! Signals last call to physics.
      
      ! OUTPUT
      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)             ! Incident LW flux to surf (W/m2).
      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)             ! Incident SW flux to surf (W/m2)
      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)          ! Absorbed SW flux by the surface (W/m2). By MT2015.
      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)              ! Outgoing LW flux to space (W/m2).
      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)              ! SW flux absorbed by the planet (W/m2).
      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)              ! Incident top of atmosphere SW flux (W/m2).
      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
      
      
!-----------------------------------------------------------------------
!     Declaration of the variables required by correlated-k subroutines
!     Numbered from top to bottom (unlike in the GCM)
!-----------------------------------------------------------------------

      REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
      REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)

      ! Optical values for the optci/cv subroutines
      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)

      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
      REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
      REAL*8 nfluxtopi_nu(L_NSPECTI)                 ! Net band-resolved IR flux at TOA (W/m2).
      REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)         ! For 1D diagnostic.
      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
      REAL*8 albi,acosz
      REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.

      INTEGER ig,l,k,nw

      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)

      logical OLRz
      real*8 NFLUXGNDV_nu(L_NSPECTV)
   
      ! Included by MT for albedo calculations.      
      REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
      REAL*8 surface_stellar_flux   ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.


!=======================================================================
!             I.  Initialization on every call   
!=======================================================================
 

      ! How much light do we get ?
      do nw=1,L_NSPECTV
         stel(nw)=stellarf(nw)/(dist_star**2)
      end do


!-----------------------------------------------------------------------    
      do ig=1,ngrid ! Starting Big Loop over every GCM column
!-----------------------------------------------------------------------


!=======================================================================
!              II.  Transformation of the GCM variables
!=======================================================================


         ! Albedo and Emissivity.
         albi=1-emis(ig)   ! Long Wave.
         DO nw=1,L_NSPECTV ! Short Wave loop.
            albv(nw)=albedo(ig,nw)
         ENDDO

      if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
         acosz = cos(pi*szangle/180.0)
         print*,'acosz=',acosz,', szangle=',szangle
      else
         acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
      endif

!-----------------------------------------------------------------------
!     Pressure and temperature
!-----------------------------------------------------------------------

      DO l=1,nlayer
         plevrad(2*l)   = pplay(ig,nlayer+1-l)/scalep
         plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
         tlevrad(2*l)   = pt(ig,nlayer+1-l)
         tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
      END DO
      
      plevrad(1) = 0.
      plevrad(2) = 0.   !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.

      tlevrad(1) = tlevrad(2)
      tlevrad(2*nlayer+1)=tsurf(ig)
      
      pmid(1) = max(pgasmin,0.0001*plevrad(3))
      pmid(2) =  pmid(1)

      tmid(1) = tlevrad(2)
      tmid(2) = tmid(1)
    
      DO l=1,L_NLAYRAD-1
         tmid(2*l+1) = tlevrad(2*l+1)
         tmid(2*l+2) = tlevrad(2*l+1)
         pmid(2*l+1) = plevrad(2*l+1)
         pmid(2*l+2) = plevrad(2*l+1)
      END DO
      pmid(L_LEVELS) = plevrad(L_LEVELS)
      tmid(L_LEVELS) = tlevrad(L_LEVELS)

!!Alternative interpolation:
!         pmid(3) = pmid(1)
!         pmid(4) = pmid(1) 
!         tmid(3) = tmid(1)
!         tmid(4) = tmid(1)
!      DO l=2,L_NLAYRAD-1
!         tmid(2*l+1) = tlevrad(2*l)
!         tmid(2*l+2) = tlevrad(2*l)
!         pmid(2*l+1) = plevrad(2*l)
!         pmid(2*l+2) = plevrad(2*l)
!      END DO
!      pmid(L_LEVELS) = plevrad(L_LEVELS-1)
!      tmid(L_LEVELS) = tlevrad(L_LEVELS-1)

      ! Test for out-of-bounds pressure.
      if(plevrad(3).lt.pgasmin)then
         print*,'Minimum pressure is outside the radiative'
         print*,'transfer kmatrix bounds, exiting.'
         call abort
      elseif(plevrad(L_LEVELS).gt.pgasmax)then
         print*,'Maximum pressure is outside the radiative'
         print*,'transfer kmatrix bounds, exiting.'
         call abort
      endif

      ! Test for out-of-bounds temperature.
      do k=1,L_LEVELS
         if(tlevrad(k).lt.tgasmin)then
            print*,'Minimum temperature is outside the radiative'
            print*,'transfer kmatrix bounds'
            print*,"k=",k," tlevrad(k)=",tlevrad(k)
            print*,"tgasmin=",tgasmin
            if (strictboundcorrk) then
              call abort
            else
              print*,'***********************************************'
              print*,'we allow model to continue with tlevrad=tgasmin' 
              print*,'  ... we assume we know what you are doing ... '
              print*,'  ... but do not let this happen too often ... '
              print*,'***********************************************'
              !tlevrad(k)=tgasmin
            endif
         elseif(tlevrad(k).gt.tgasmax)then
!            print*,'Maximum temperature is outside the radiative'
!            print*,'transfer kmatrix bounds, exiting.'
!            print*,"k=",k," tlevrad(k)=",tlevrad(k)
!            print*,"tgasmax=",tgasmax
            if (strictboundcorrk) then
              call abort
            else
!              print*,'***********************************************'
!              print*,'we allow model to continue with tlevrad=tgasmax'  
!              print*,'  ... we assume we know what you are doing ... '
!              print*,'  ... but do not let this happen too often ... '
!              print*,'***********************************************'
              !tlevrad(k)=tgasmax
            endif
         endif
      enddo
      do k=1,L_NLAYRAD+1
         if(tmid(k).lt.tgasmin)then
            print*,'Minimum temperature is outside the radiative'
            print*,'transfer kmatrix bounds, exiting.'
            print*,"k=",k," tmid(k)=",tmid(k)
            print*,"tgasmin=",tgasmin
            if (strictboundcorrk) then
              call abort
            else
              print*,'***********************************************'
              print*,'we allow model to continue with tmid=tgasmin'
              print*,'  ... we assume we know what you are doing ... '
              print*,'  ... but do not let this happen too often ... '
              print*,'***********************************************'
              tmid(k)=tgasmin
            endif
         elseif(tmid(k).gt.tgasmax)then
!            print*,'Maximum temperature is outside the radiative'
!            print*,'transfer kmatrix bounds, exiting.'
!            print*,"k=",k," tmid(k)=",tmid(k)
!            print*,"tgasmax=",tgasmax
            if (strictboundcorrk) then
              call abort
            else
!              print*,'***********************************************'
!              print*,'we allow model to continue with tmid=tgasmin'
!              print*,'  ... we assume we know what you are doing ... '
!              print*,'  ... but do not let this happen too often ... '
!              print*,'***********************************************'
              tmid(k)=tgasmax
            endif
         endif
      enddo

!=======================================================================
!          III. Calling the main radiative transfer subroutines
!=======================================================================


         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
         glat_ig=glat(ig)

!-----------------------------------------------------------------------
!        Short Wave Part
!-----------------------------------------------------------------------

         if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
            if((ngrid.eq.1).and.(global1d))then
               do nw=1,L_NSPECTV
                  stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
               end do
            else
               do nw=1,L_NSPECTV
                  stel_fract(nw)= stel(nw) * fract(ig)
               end do
            endif
            
            call optcv(pq(ig,:,1:nmicro),nlayer,plevrad,tmid,pmid,        &
                 dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf)

            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
                 acosz,stel_fract,                                 &
                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,   &
                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)

         else ! During the night, fluxes = 0.
            nfluxtopv       = 0.0d0
            fluxtopvdn      = 0.0d0
            nfluxoutv_nu(:) = 0.0d0
            nfluxgndv_nu(:) = 0.0d0
            do l=1,L_NLAYRAD
               fmnetv(l)=0.0d0
               fluxupv(l)=0.0d0
               fluxdnv(l)=0.0d0
            end do
         end if


         ! Equivalent Albedo Calculation (for OUTPUT). MT2015
         if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.       
            surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))      
            if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
               DO nw=1,L_NSPECTV                  
                  albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
               ENDDO
               albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
               albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
            else
               albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
            endif
         else
            albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
         endif


!-----------------------------------------------------------------------
!        Long Wave Part
!-----------------------------------------------------------------------

         call optci(pq(ig,:,1:nmicro),nlayer,plevrad,tlevrad,tmid,pmid,    &
              dtaui,taucumi,cosbi,wbari,taugsurfi)

         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
              wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         & 
              fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)

!-----------------------------------------------------------------------
!     Transformation of the correlated-k code outputs
!     (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)

!     Flux incident at the top of the atmosphere
         fluxtop_dn(ig)=fluxtopvdn 

         fluxtop_lw(ig)  = real(nfluxtopi)
         fluxabs_sw(ig)  = real(-nfluxtopv)
         fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
         fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
         
!        Flux absorbed by the surface. By MT2015.         
         fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))

         if(fluxtop_dn(ig).lt.0.0)then
            print*,'Achtung! fluxtop_dn has lost the plot!'
            print*,'fluxtop_dn=',fluxtop_dn(ig)
            print*,'acosz=',acosz
            print*,'temp=   ',pt(ig,:)
            print*,'pplay=  ',pplay(ig,:)
            call abort
         endif

!     Spectral output, for exoplanet observational comparison
         if(specOLR)then
            do nw=1,L_NSPECTI 
               OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
            end do
            do nw=1,L_NSPECTV 
               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
               OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
            end do
         endif

!     Finally, the heating rates

         DO l=2,L_NLAYRAD
            dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
            dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
         END DO      

!     These are values at top of atmosphere
         dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
         dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))


!-----------------------------------------------------------------------    
      end do ! End of big loop over every GCM column.
!-----------------------------------------------------------------------



!-----------------------------------------------------------------------
!     Additional diagnostics
!-----------------------------------------------------------------------

      ! IR spectral output, for exoplanet observational comparison
      if(lastcall.and.(ngrid.eq.1))then  ! could disable the 1D output, they are in the diagfi and diagspec... JL12

         print*,'Saving scalar quantities in surf_vals.out...'
         print*,'psurf = ', pplev(1,1),' Pa'
         open(116,file='surf_vals.out')
         write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
                      real(-nfluxtopv),real(nfluxtopi) 
         close(116)


!          USEFUL COMMENT - Do Not Remove.
!
!           if(specOLR)then
!               open(117,file='OLRnu.out')
!               do nw=1,L_NSPECTI
!                  write(117,*) OLR_nu(1,nw)
!               enddo
!               close(117)
!
!               open(127,file='OSRnu.out')
!               do nw=1,L_NSPECTV
!                  write(127,*) OSR_nu(1,nw)
!               enddo
!               close(127)
!           endif

           ! OLR vs altitude: do it as a .txt file.
         OLRz=.false.
         if(OLRz)then
            print*,'saving IR vertical flux for OLRz...'
            open(118,file='OLRz_plevs.out')
            open(119,file='OLRz.out')
            do l=1,L_NLAYRAD
               write(118,*) plevrad(2*l)
               do nw=1,L_NSPECTI
                  write(119,*) fluxupi_nu(l,nw) 
               enddo
            enddo 
            close(118)
            close(119)
         endif

      endif

      ! See physiq.F for explanations about CLFvarying. This is temporary.
      if (lastcall) then
        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
        IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
!$OMP BARRIER
!$OMP MASTER
        IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
        IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
        IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
!$OMP END MASTER
!$OMP BARRIER
      endif


    end subroutine callcorrk