source: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 @ 832

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

      use radinc_h
      use radcommon_h
      use watercommon_h
      use datafile_mod, only: datadir
      use ioipsl_getincom 
      use gases_h
      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
      use aerosol_mod

      USE tracer_h

      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"

!-----------------------------------------------------------------------
!     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(ngrid,nlayermx,L_NSPECTV,naerkind)
      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)

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

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

      !!! this is a local instance of a variable saved in physiq.F and being an argument
      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
      REAL, INTENT(INOUT) :: nueffrad(ngrid,nlayermx,naerkind)

!!     OUTPUT
      REAL dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
      REAL dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
      REAL fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
      REAL fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
      REAL fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
      REAL fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
      REAL fluxtop_dn(ngrid)    ! incident top of atmosphere SW flux (W/m2)
      REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
      REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)

!-----------------------------------------------------------------------
!     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 szangle
      logical global1d
      save szangle,global1d
      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(ngrid,nlayer,nq)
      REAL qsurf(ngrid,nq)       ! 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

      !REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
      !REAL :: QREFir3d(ngrid,nlayermx,naerkind)
      !save QREFvis3d, QREFir3d 
      real, dimension(:,:,:), save, allocatable :: QREFvis3d
      real, dimension(:,:,:), save, allocatable :: 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!!!

      !real ptime, pday
      logical OLRz
      real*8 NFLUXGNDV_nu(L_NSPECTV)

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

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

      real maxrad, minrad
            
      real CBRT
      external CBRT

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

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

         ALLOCATE(QREFvis3d(ngrid,nlayermx,naerkind))
         ALLOCATE(QREFir3d(ngrid,nlayermx,naerkind))

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

!--------------------------------------------------
!     Effective radius and variance of the aerosols
         if(naerkind.gt.4)then
            print*,'Code not general enough to deal with naerkind > 4 yet.'
            call abort
         endif
         call su_aer_radii(ngrid,reffrad,nueffrad)
         
         

!--------------------------------------------------
!     set up correlated k
         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

         OLR_nu(:,:) = 0.
         OSR_nu(:,:) = 0.

         if (ngrid.eq.1) then
           PRINT*, 'Simulate global averaged conditions ?'
           global1d = .false. ! default value
           call getin("global1d",global1d)
           write(*,*) "global1d = ",global1d
           ! Test of incompatibility:
           ! if global1d is true, there should not be any diurnal cycle
           if (global1d.and.diurnal) then
            print*,'if global1d is true, diurnal must be set to false'
            stop
           endif

           if (global1d) then
             PRINT *,'Solar Zenith angle (deg.) ?'
             PRINT *,'(assumed for averaged solar flux S/4)'
             szangle=60.0  ! default value
             call getin("szangle",szangle)
             write(*,*) "szangle = ",szangle
           endif
         endif

         firstcall=.false.   

      end if

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

!--------------------------------------------------
!     Effective radius and variance of the aerosols
      do iaer=1,naerkind

         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
            call co2_reffrad(ngrid,nq,pq,reffrad)
            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
         end if
         if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
            call h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
         endif
         if(iaer.eq.iaero_dust)then
            call dust_reffrad(ngrid,reffrad)
            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
         endif
         if(iaer.eq.iaero_h2so4)then
            call h2so4_reffrad(ngrid,reffrad)
            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
         endif
      end do !iaer=1,naerkind


!     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,ngrid

!=======================================================================
!     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' 
                        print*,'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'
                        print*,'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.d0)
               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
!
            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

!     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 ((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

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

      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)
         !qvar=0.005                   ! completely constant profile (JL)

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

      if(.not.kastprof)then
      ! IMPORTANT: Now convert from kg/kg to mol/mol
         do k=1,L_LEVELS
            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
         end do
      end if

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

         ! initial values equivalent to mugaz
         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'
         print*,'Assumes that the variable gas is H2O!!!'
         print*,'Assumes that there is only one tracer'
         !i_var=igcm_h2o_vap
         i_var=1
         if(nq.gt.1)then
            print*,'Need 1 tracer only to run kcm1d.e' 
            stop
         endif
         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)

         print*,'Warning: reducing qvar in callcorrk.'
         print*,'Temperature profile no longer consistent ', &
                            'with saturated H2O.'
         do k=1,L_LEVELS
            qvar(k) = qvar(k)*satval
         end do

      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*,'level,grid,T',k,ig,tlevrad(k)
            print*,'WARNING, OVERRIDING FOR TEST'
            !call abort
            tlevrad(k)=tgasmax
         endif
      enddo

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


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

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

            fluxtoplanet=0.

            if((ngrid.eq.1).and.(global1d))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,              &
                 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)/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))  &
                *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)))

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


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

!          I am useful, please don`t remove me!
!           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 .and. .not.CLFvarying) then
        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
        IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
        IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
        IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
        IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
        IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
      endif


    end subroutine callcorrk