subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
          albedo,emis,mu0,pplev,pplay,pt,                      & 
          tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
          dtlw,dtsw,fluxsurf_lw,                               &
          fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
          reffrad,tau_col,cloudfrac,totcloudfrac,   &
          clearsky,firstcall,lastcall)


      use radinc_h
      use radcommon_h
      use watercommon_h
      use datafile_mod, only: datadir
      use ioipsl_getincom 

      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)
!
!==================================================================

#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"
#include "gases.h"

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

!     INPUT
      INTEGER icount
      INTEGER ngrid,nlayer
      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
      REAL albedo(ngrid)                    ! SW albedo
      REAL emis(ngrid)                      ! LW emissivity
      REAL pplay(ngrid,nlayermx)            ! pres. level in GCM mid of layer
      REAL pplev(ngrid,nlayermx+1)          ! pres. level at GCM layer boundaries

      REAL pt(ngrid,nlayermx)               ! air temperature (K)
      REAL tsurf(ngrid)                     ! surface temperature (K)
      REAL dist_star,mu0(ngrid)             ! distance star-planet (AU)
      REAL fract(ngrid)                     ! fraction of day

!     Globally varying aerosol optical properties on GCM grid
!     Not needed everywhere so not in radcommon_h
      REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
      REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
      REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)

      REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
      REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
      REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)

      REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
      REAL :: QREFir3d(ngridmx,nlayermx,naerkind)

!      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
!      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...

      REAL reffrad(ngrid,nlayer,naerkind)
      REAL nueffrad(ngrid,nlayer,naerkind)

!     OUTPUT
      REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
      REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
      REAL fluxsurf_lw(ngridmx)   ! incident LW flux to surf (W/m2)
      REAL fluxtop_lw(ngridmx)    ! outgoing LW flux to space (W/m2)
      REAL fluxsurf_sw(ngridmx)   ! incident SW flux to surf (W/m2)
      REAL fluxabs_sw(ngridmx)    ! SW flux absorbed by planet (W/m2)
      REAL fluxtop_dn(ngridmx)    ! incident top of atmosphere SW flux (W/m2)

!-----------------------------------------------------------------------
!     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 tauaero(L_LEVELS+1,naerkind)
      REAL*8 nfluxtopv,nfluxtopi,nfluxtop
      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,albv,acosz

      INTEGER ig,l,k,nw,iaer,irad

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

      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
      REAL pq(ngridmx,nlayer,nq)
      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
      integer nq

!     Local aerosol optical properties for each column on RADIATIVE grid
      real*8  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
      real*8  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
      real*8  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
      real*8  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
      real*8  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
      real*8  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)

      save qxvaer, qsvaer, gvaer
      save qxiaer, qsiaer, giaer
      save QREFvis3d, QREFir3d 

      REAL tau_col(ngrid) ! diagnostic from aeropacity

!     Misc.
      logical firstcall, lastcall, nantest
      real*8  tempv(L_NSPECTV)
      real*8  tempi(L_NSPECTI)
      real*8  temp,temp1,temp2,pweight
      character(len=10) :: tmp1
      character(len=10) :: tmp2

!     for fixed water vapour profiles
      integer i_var
      real RH
      real*8 pq_temp(nlayer)
      real ptemp, Ttemp, qsat

!      real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!

!     for OLR spec
      integer OLRcount
      save OLRcount
      integer OLRcount2
      save OLRcount2
      character(len=2)  :: tempOLR
      character(len=30) :: filenomOLR
      !real ptime, pday
      logical OLRz
      real OLR_nu(ngrid,L_NSPECTI)
      !real GSR_nu(ngrid,L_NSPECTV)
      real OSR_nu(ngrid,L_NSPECTV)
      real*8 NFLUXGNDV_nu(L_NSPECTV)

!     for H2O cloud fraction in aeropacity
      real cloudfrac(ngridmx,nlayermx)
      real totcloudfrac(ngridmx)
      logical clearsky

!     Allow variations in cp with location
      real cpp3D(ngridmx,nlayermx)   ! specific heat capacity at const. pressure

      ! for weird cloud test
      real pqtest(ngridmx,nlayer,nq)

!     are particle radii fixed?
      logical radfixed
      real maxrad, minrad

      real CBRT
      external CBRT

!     included by RW for runaway greenhouse 1D study
      real muvar(ngridmx,nlayermx+1)
      real vtmp(nlayermx)
      REAL*8 muvarrad(L_LEVELS)

      radfixed=.false.

!=======================================================================
!     Initialization on first call

      qxvaer(:,:,:)=0.0
      qsvaer(:,:,:)=0.0
      gvaer(:,:,:) =0.0

      qxiaer(:,:,:)=0.0
      qsiaer(:,:,:)=0.0
      giaer(:,:,:) =0.0

      if(firstcall) then

         call system('rm -f surf_vals_long.out')

!--------------------------------------------------
!     Effective radius and variance of the aerosols

!     CO2 ice:
         DO l=1,nlayer
            DO ig=1,ngrid
               reffrad(ig,l,1)  = 1.e-4
               nueffrad(ig,l,1) = 0.1 
!     these values will change once the microphysics gets to work
!     UNLESS tracer=.false., in which case we should be working with
!     a fixed aerosol layer, and be able to define reffrad in a 
!     .def file. To be improved!
            ENDDO
         ENDDO

!     H2O ice:
         if(naerkind.eq.2)then
            DO l=1,nlayer
               DO ig=1,ngrid
                  reffrad(ig,l,naerkind)  = 1.e-5
                  nueffrad(ig,l,naerkind) = 0.1
               ENDDO  
            ENDDO
         endif

         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
         call getin("corrkdir",corrkdir)
         print*, "corrkdir = ",corrkdir
         write( tmp1, '(i3)' ) L_NSPECTI
         write( tmp2, '(i3)' ) L_NSPECTV
         banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
         banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))

         call sugas_corrk       ! set up gaseous absorption properties
         call setspi            ! basic infrared properties
         call setspv            ! basic visible properties
         call suaer_corrk       ! set up aerosol optical properties

         Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed


         if((igcm_h2o_vap.eq.0) .and. varactive)then
            print*,'varactive in callcorrk but no h2o_vap tracer.'
            stop
         endif

         firstcall=.false.   

      end if

!=======================================================================
!     Initialization on every call    

      do l=1,nlayer
         do ig=1,ngrid
            do iaer=1,naerkind
               nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
            enddo
         enddo
      enddo


      if(kastprof)then
         radfixed=.true.
      endif

      if(radfixed)then
         do l=1,nlayer
            do ig=1,ngrid
               reffrad(ig,l,1)    = 5.e-5 ! CO2 ice
            enddo
         enddo
         print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
         if(naerkind.ge.2)then
            do l=1,nlayer
               do ig=1,ngrid
                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
                  reffrad(ig,l,2) = 5.e-6 ! H2O ice
               enddo
            enddo
            print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
         endif
         if(naerkind.eq.3)then
            do l=1,nlayer
               do ig=1,ngrid
                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
               enddo
            enddo
            print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
         endif
         if(naerkind.gt.3)then
            print*,'Code not general enough to deal with naerkind > 3 yet.'
            call abort
         endif

      else

         maxrad=0.0
         !averad=0.0
         minrad=1.0
         do l=1,nlayer

            !masse = (pplev(ig,l) - pplev(ig,l+1))/g

            do ig=1,ngrid
               if(tracer)then
                  reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
                       (4*Nmix_co2*pi*rho_co2) )
               endif
               reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
               reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)

               !averad = averad + reffrad(ig,l,1)*area(ig)
               maxrad = max(reffrad(ig,l,1),maxrad)
               minrad = min(reffrad(ig,l,1),minrad)
            enddo
         enddo
         if(igcm_co2_ice.gt.0)then
            print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
            print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
         endif

         if((naerkind.ge.2).and.water)then
            maxrad=0.0
            minrad=1.0
            do l=1,nlayer
               do ig=1,ngrid
                  reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
                       (4*Nmix_h2o*pi*rho_ice) )
                  reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
                  reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)

                  maxrad = max(reffrad(ig,l,2),maxrad)
                  minrad = min(reffrad(ig,l,2),minrad)
               enddo
            enddo
            print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
            print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
         endif

         if(naerkind.eq.3)then
            do l=1,nlayer
               do ig=1,ngrid
                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
               enddo
            enddo
         endif

      endif


!     how much light we get
      fluxtoplanet=0
      do nw=1,L_NSPECTV
         stel(nw)=stellarf(nw)/(dist_star**2)
         fluxtoplanet=fluxtoplanet + stel(nw)
      end do

      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
           QIRsQREF3d,omegaIR3d,gIR3d,                             &
           QREFvis3d,QREFir3d)                                     ! get 3D aerosol optical properties

      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
           reffrad,QREFvis3d,QREFir3d,                             & 
           tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths


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

!=======================================================================
!     Transformation of the GCM variables

!-----------------------------------------------------------------------
!     Aerosol optical properties Qext, Qscat and g
!     The transformation in the vertical is the same as for temperature
           
!     shortwave
            do iaer=1,naerkind
               DO nw=1,L_NSPECTV 
                  do l=1,nlayermx

                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
                         *QREFvis3d(ig,nlayermx+1-l,iaer)

                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
                         *QREFvis3d(ig,max(nlayermx-l,1),iaer)

                     qxvaer(2*l,nw,iaer)  = temp1
                     qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                     temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
                     temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)

                     qsvaer(2*l,nw,iaer)  = temp1
                     qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                     temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
                     temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)

                     gvaer(2*l,nw,iaer)  = temp1
                     gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                  end do

                  qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
                  qxvaer(2*nlayermx+1,nw,iaer)=0.

                  qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
                  qsvaer(2*nlayermx+1,nw,iaer)=0.

                  gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
                  gvaer(2*nlayermx+1,nw,iaer)=0.

               end do


!     longwave
               DO nw=1,L_NSPECTI 
                  do l=1,nlayermx

                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
                          *QREFir3d(ig,nlayermx+1-l,iaer)

                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
                          *QREFir3d(ig,max(nlayermx-l,1),iaer)

                     qxiaer(2*l,nw,iaer)  = temp1
                     qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                     temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
                     temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)

                     qsiaer(2*l,nw,iaer)  = temp1
                     qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                     temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
                     temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)

                     giaer(2*l,nw,iaer)  = temp1
                     giaer(2*l+1,nw,iaer)=(temp1+temp2)/2

                  end do

                  qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
                  qxiaer(2*nlayermx+1,nw,iaer)=0.

                  qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
                  qsiaer(2*nlayermx+1,nw,iaer)=0.

                  giaer(1,nw,iaer)=giaer(2,nw,iaer)
                  giaer(2*nlayermx+1,nw,iaer)=0.

               end do
            end do


            ! test / correct for freaky s. s. albedo values
            do iaer=1,naerkind
               do k=1,L_LEVELS+1

                  do nw=1,L_NSPECTV
                     if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
                        print*,'Serious problems with qsvaer values in callcorrk'
                        call abort
                     endif
                     if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
                        qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
                     endif
                  end do

                  do nw=1,L_NSPECTI 
                     if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
                        print*,'Serious problems with qsiaer values in callcorrk'
                        call abort
                     endif
                     if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
                        qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
                     endif
                  end do

               end do
            end do

!-----------------------------------------------------------------------
!     Aerosol optical depths
            
         do iaer=1,naerkind     ! a bug was here           
            do k=0,nlayer-1
               
               pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))

               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)

               tauaero(2*k+2,iaer)=max(temp*pweight,0.0)
               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.0)
!
            end do
            ! boundary conditions
            tauaero(1,iaer)          = tauaero(2,iaer)
            tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
            !tauaero(1,iaer)          = 0.
            !tauaero(L_LEVELS+1,iaer) = 0.
         end do
         !print*,'Note changed tauaero BCs in callcorrk!'

!     Albedo and emissivity
         albi=1-emis(ig)        ! longwave
         albv=albedo(ig)        ! shortwave 

      if(noradsurf.and.(albv.gt.0.0))then
         print*,'For open lower boundary in callcorrk must'
         print*,'have surface albedo set to zero!'
         call abort
      endif

      if(ngridmx.eq.1) then       ! fixed zenith angle 'szangle' in 1D
         acosz = cos(pi*szangle/180.0)
         print*,'acosz=',acosz,', szangle=',szangle
      else
         acosz=mu0(ig)          ! cosine of sun incident angle
      endif


!-----------------------------------------------------------------------
!     Water vapour (to be generalised for other gases eventually)
      
      if(varactive)then

         i_var=igcm_h2o_vap
         do l=1,nlayer
            qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
            ! Average approximation as for temperature...
         end do
         qvar(1)=qvar(2)
         qvar(2*nlayermx+1)=qsurf(ig,i_var)

      elseif(varfixed)then

         do l=1,nlayermx        ! here we will assign fixed water vapour profiles globally
            RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
            if(RH.lt.0.0) RH=0.0
            
            ptemp=pplay(ig,l)
            Ttemp=pt(ig,l)
            call watersat(Ttemp,ptemp,qsat)

            !pq_temp(l) = qsat      ! fully saturated everywhere
            pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
         end do
         
         do l=1,nlayer
            qvar(2*l)   = pq_temp(nlayer+1-l)
            qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
         end do
         qvar(1)=qvar(2)

         ! Lowest layer of atmosphere
         RH = satval * (1 - 0.02) / 0.98
         if(RH.lt.0.0) RH=0.0

         ptemp = pplev(ig,1)
         Ttemp = tsurf(ig)
         call watersat(Ttemp,ptemp,qsat)

         !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
         qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)

      else
         do k=1,L_LEVELS
            qvar(k) = 1.0D-7
         end do
      end if

      ! IMPORTANT: Now convert from kg/kg to mol/mol
      do k=1,L_LEVELS
         qvar(k) = qvar(k)/epsi
      end do



!-----------------------------------------------------------------------
!     kcm mode only
      if(kastprof)then

         DO l=1,nlayer
            muvarrad(2*l)   = mugaz
            muvarrad(2*l+1) = mugaz
         END DO

         do k=1,L_LEVELS
            qvar(k) = 0.0
         end do
         print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
      endif


      if(kastprof.and.(ngasmx.gt.1))then

         DO l=1,nlayer
            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
                                muvar(ig,max(nlayer+1-l,1)))/2
         END DO
      
         muvarrad(1) = muvarrad(2)
         muvarrad(2*nlayermx+1)=muvar(ig,1)

         print*,'Recalculating qvar with VARIABLE epsi for kastprof'
         i_var=igcm_h2o_vap
         do l=1,nlayer
            vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
         end do

         do l=1,nlayer
            qvar(2*l)   = vtmp(nlayer+1-l)
            qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
         end do
         qvar(1)=qvar(2)
         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O


      endif


      ! Keep values inside limits for which we have radiative transfer coefficients
      if(L_REFVAR.gt.1)then ! there was a bug here!
         do k=1,L_LEVELS
            if(qvar(k).lt.wrefvar(1))then
               qvar(k)=wrefvar(1)+1.0e-8
            elseif(qvar(k).gt.wrefvar(L_REFVAR))then
               qvar(k)=wrefvar(L_REFVAR)-1.0e-8
            endif
         end do
      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) = max(pgasmin,0.0001*plevrad(3))

      tlevrad(1) = tlevrad(2)
      tlevrad(2*nlayermx+1)=tsurf(ig)
      
      tmid(1) = tlevrad(2)
      tmid(2) = tlevrad(2)
      pmid(1) = plevrad(2)
      pmid(2) = plevrad(2)
      
      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)


      ! 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, exiting.'
            print*,'WARNING, OVERRIDING FOR TEST'
            !call abort
         elseif(tlevrad(k).gt.tgasmax)then
            print*,'Maximum temperature is outside the radiative'
            print*,'transfer kmatrix bounds, exiting.'
            print*,'WARNING, OVERRIDING FOR TEST'
            !print*, 'T=',pt
            !call abort
         endif
      enddo

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


!-----------------------------------------------------------------------
!     Shortwave

         if(fract(ig) .ge. 1.0e-4) then ! only during daylight!

            fluxtoplanet=0.

            if((ngridmx.eq.1).and.(.not.(diurnal.or.tlocked)))then
               do nw=1,L_NSPECTV
                  stel_fract(nw)= stel(nw) * 0.25 / acosz
                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
                                ! globally averaged = divide by 4
                                ! but we correct for solar zenith angle
               end do
            else
               do nw=1,L_NSPECTV
                  stel_fract(nw)= stel(nw) * fract(ig)
                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
               end do
            endif

            call optcv(dtauv,tauv,taucumv,plevrad,                 &
                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
                 tmid,pmid,taugsurf,qvar,muvarrad)


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

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




!-----------------------------------------------------------------------
!     Longwave

         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
              taugsurfi,qvar,muvarrad)

         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
              wnoi,dwni,cosbi,wbari,gweight,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)=fluxdnv(1)

         fluxtop_lw(ig)  = real(nfluxtopi)
         fluxabs_sw(ig)  = real(-nfluxtopv)
         fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
         fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))


         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*,'aerosol=',aerosol(ig,:,:)
            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)
            end do
            do nw=1,L_NSPECTV 
               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
               OSR_nu(ig,nw)=nfluxoutv_nu(nw)
            end do
         endif

!     Finally, the heating rates
         if(nonideal)then

            DO l=2,L_NLAYRAD
               dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
                    *g/(cpp3D(ig,L_NLAYRAD+1-l)                &
                    *scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
               dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
                    *g/(cpp3D(ig,L_NLAYRAD+1-l)                &
                    *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)           &
                 *g/(cpp3D(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(1)))
            dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
                 *g/(cpp3D(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(1)))

         else

            DO l=2,L_NLAYRAD
               dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
                   *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
               dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
                   *g/(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)           &
                *g/(cpp*scalep*(plevrad(3)-plevrad(1)))
            dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
                *g/(cpp*scalep*(plevrad(3)-plevrad(1)))

         endif

!     ---------------------------------------------------------------
      end do                    ! end of big loop over every GCM column (ig = 1:ngrid)


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

!     IR spectral output, for exoplanet observational comparison
      if(specOLR)then
        if(ngrid.ne.1)then
          !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
          !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
        endif
      endif

      if(lastcall.and.(ngrid.eq.1))then

           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)

           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

    end subroutine callcorrk