C**Revisions:
C**Added MTSAT1 logic 5/2006
c23456789112345678921234567893123456789412345678951234567896123456789712
c***************************************************************************
      subroutine InitCKDist_ckm(FileNames,DataType,psc3temp,pSolZenMax,
     &                      pTwilit,i85modis)
c******INPUTS
c      FileNames(1)   =   virs.corrk.coefs.plus.n2o
c                         goes.corrk.coefs
c                         modis.corrk.coefs
c                         mtsat.corrk.coefs
c      FileNames(2)   =   SkinTderiv.ch4.allvz.dy.dat 
c      FileNames(3)   =   LNPWderiv.ch4.allvz.dy.dat
c      DataType = 0=AVHRR/N12,14 1=AVHRR/N15-17, 2=VIRS, 3=TERRA, 6=AQUA, 11=G10TSA, 
c                 21=G7, 22=G8, 23=G10, 24=G12, 25=G09, 26=G11, 40=MSG1, 50=MTSAT-1R
c      psc3temp = instrument's near IR solar constant (3.7-3.9um) in K units
c                 Setting to <=0.0 will use defaults in AssignCommonBlockConst_ckm.
c      pSolZenMax = maximum angle for daytime processing (82.0 deg)
c                   Not currently used.
c      pTwilit = maximum angle for twilight processing (87.5)
c      i85modis = 1 For modis 8.5um band (stratospheric ozone)
c               = 0 Do not do modis 8.5um band
c****table of near-ir solar constants (psc3temp)
c       347.00 = noaa10 and used for TOGA too
c       357.37 = noaa?, kansas, 1991 from dave (and apr94)
c       362.00 = noaa9 - don't know when
c       352.00 = noaa14
c       347.00 = noaa9 1986 global runs for 4/95 CERES STM
c       342.00 = noaa9 1986 global runs for 11/97 according to Minnis
c       359.64 = noaa9 1986 as of 11/97 according to Arduini
c       344.60 = GOES8 according to Arduini
c       355.40 = VIRS 1/98 according to Sunny
c       353.21 = MODIS 10/00 according to Sunny

      include './corrk.h'
      IMPLICIT NONE
      
      integer DataType,mlv,i85flag,i85modis

      real psc3temp,pSolZenMax,pTwilit
      real sc3temp,bsun,solzenmax,twilit
      real p0(MAXCKLEV),t0(MAXCKLEV),u0(MAXCKLEV)
      real z0(MAXCKLEV),rh0(MAXCKLEV),o30(MAXCKLEV)
      character*200 FileNames(10)

      COMMON /ktemps/ mlv,p0,t0,z0,u0,rh0,o30,i85flag
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit

      call AssignCommonBlockConst_ckm(DataType,pSolZenMax,pTwilit)
      call Read_Corrk_ckm (FileNames, DataType)

      if(psc3temp .gt. 0.0)then  !user-input near-ir solar const
        sc3temp=psc3temp
      endif

      i85flag=i85modis

      return
      end
c*****************************************************************************
      subroutine MoveBTReg_ckm(iupdown,iradio,SBT,SRef,klevs,sza,wave,
     &                         emit,chi,FBT,FRef)
c*****************************************************************************
c      Purpose:  Move surface skin T or reflectance to TOA using the
c                number of levels input by the user (klevs).  
c                Also moves TOA BT to the sfc.
c ***** INPUT
c     SBT = starting BT(K)
c     SRef = starting reflectance of near-ir channel to move to TOA/Sfc
c            Set to -1, or <0 if it's not being used.  Currently used
c            for Bob Arduini's bidir snow reflectance (3.7um) model.
c     sza = solar zenith angle (deg)
c     wave = wavelength (um) (see AssignCommonBlockConst_ckm)
c     emit = surface emittance for the wavelength.
c     chi = bidirectional model; used for near-ir channel in daytime
c     klevs = number of levels
c     iupdown=solve up(1) or down(2)
c     iradio=radiometric (1) or true Tskin (2) for down direction.

c ***** OUTPUT
c     FBT = Final Brightness Temperature
c     FRef = Final reflectance of near-ir channel.

      IMPLICIT NONE

      INTEGER  klevs, icorrK, laybot, laytop, ichan, iupdown, iradio
      REAL  sza, wave
      REAL  solzenmax, bsun,twilit, sc3temp 
      REAL  emit, Tfinal, SBT, TOAout
      REAL  SRef
      REAL  chi, B3SurfRad, SurfRef, FBT, FRef

      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
          
      icorrK=1
      laybot = 1        ! bottom layer
      laytop = klevs-1  ! top layer (layers=levels-1)

      if(wave .lt. 5.0)ichan=1
      if(wave .ge. 10.0 .and. wave .le. 11.4)ichan=2
      if(wave .ge. 11.8 .and. wave .le. 13.5)ichan=3
      if(wave .ge. 8.0 .and. wave .le. 9.0)ichan=4


c*** UseCorrKTempE_ckm: Use emissivity
c*** UseCorrKTemp_ckm: Assumes emissivity is 1

      if(iupdown .eq. 2 .and. iradio .eq. 1)emit=1.0

      if(wave .lt. 5.0 .and. sza .le. twilit)then !daytime solar+ir part
        icorrK=1
        if(iupdown .eq. 1)then
          call GetT3TOA_ckm(SBT,emit,chi,sza,icorrK,B3SurfRad,SurfRef,FBT)
        else         
          Call UseCorrKTempE_ckm(SBT,laybot,laytop,iupdown,ichan,FBT,emit)  
        endif
      else            
        Call UseCorrKTempE_ckm(SBT,laybot,laytop,iupdown,ichan,FBT,emit)  
      endif


      if(wave .lt. 5.0)then !Reflectance
        if(sza .le. twilit .and. SRef .gt. 0.0)then !daytime solar part
          call UseCorrKRef_ckm(SRef,laybot,laytop,iupdown,FRef)   
        elseif(sza .gt. twilit)then
          FRef=0.0
        else
          FRef=-0.999
        endif       
      endif
 200  format(f8.2,4i5,2f8.2)

      return
      end

c***********************************************************************
      SUBROUTINE Read_Corrk_ckm (FileNames, DataType)
c*********************************************************************** 
c  !F77
c  Name: Subroutine: Read_Corrk_ckm
c  
c  Purpose:
c  Read data files for Corrk algorithm
c    
c  !Input Parameters:
c  character*200 FileNames(10)
c  integer DataType
c     FileNames(1)   =    virs.corrk.coefs.plus.n2o (AVHRR, VIRS)
c                         goes.corrk.coefs
c                         modis.corrk.coefs
c     FileNames(2)   =    SkinTderiv.ch4.allvz.dy.dat 
c     FileNames(3)   =    LNPWderiv.ch4.allvz.dy.dat
c     FileNames(4)   =    sno3.75.bin
c   
c  !Output parameters:
c  none
c  
c  External Files:
c  None
c   
c  Called Modules: 
c  none
c  
c  Local Variables:
c  see declaration
c
c  !Revision History:
c  6/2003 : Initial Version
c  
c  !END
c***********************************************************************

      IMPLICIT NONE

c *** input parameters ***
      integer DataType
      character*200 FileNames(10)
      integer i,j

c *** output parameters ***
c     none


c *** local parameters ***

      integer num_pressure,max_coefs,iorder
      integer num_segment_3 
      integer i_goes_segment, i_modis_segment, num_seg 
      integer icoef,kpres,iseg,jorder,nsza,nvza,nraz

      parameter(nsza=21, nvza=21, nraz=24)
      real sza(nsza),vza(nvza),raz(nraz)
      real alb(nsza),dfalb
      real refl(nsza,nvza,nraz),brdf(nsza,nvza,nraz)
     
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      parameter (num_segment_3 = 5)
      parameter (i_goes_segment = 1)
      parameter (i_modis_segment = 1)

      integer n_h2o_3(num_segment_3),n_ch4_3(num_segment_3),n_n2o_3(num_segment_3),
     +  n_h2o_4,n_co2_4,n_h2o_5,n_co2_5

      real f_h2o_3(max_coefs,num_segment_3), f_h2o_4(max_coefs), f_h2o_5(max_coefs)
      real f_ch4_3(max_coefs,num_segment_3), f_co2_4(max_coefs), f_co2_5(max_coefs)
      real f_n2o_3(max_coefs,num_segment_3)
      real k_h2o_3(max_coefs,num_segment_3), k_h2o_4(max_coefs), k_h2o_5(max_coefs)
      real k_ch4_3(max_coefs,num_segment_3), k_co2_4(max_coefs), k_co2_5(max_coefs)
      real k_n2o_3(max_coefs,num_segment_3)
      real coefk_h2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_h2o_4(max_coefs,iorder,num_pressure)
      real coefk_h2o_5(max_coefs,iorder,num_pressure)
      real coefk_ch4_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_n2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_4(max_coefs,iorder,num_pressure)
      real coefk_co2_5(max_coefs,iorder,num_pressure)

      real dtoadskin(30,20,6),dtoadlnpw(30,20,6)
      

c *** common blocks ***   
      common/corrkcoef/ 
     +  n_h2o_3, f_h2o_3, k_h2o_3, coefk_h2o_3,	    ! from virs.corrk.coefs.plus.n2o,
     +  n_ch4_3, f_ch4_3, k_ch4_3, coefk_ch4_3,     ! mtsat.corrk.coefs
     +  n_n2o_3, f_n2o_3, k_n2o_3, coefk_n2o_3,
     +  n_h2o_4, f_h2o_4, k_h2o_4, coefk_h2o_4,
     +  n_co2_4, f_co2_4, k_co2_4, coefk_co2_4,
     +  n_h2o_5, f_h2o_5, k_h2o_5, coefk_h2o_5,
     +  n_co2_5, f_co2_5, k_co2_5, coefk_co2_5

      common/derivs/dtoadskin,dtoadlnpw             ! from SkinTderiv.ch4.allvz.dy.dat or
                                                    !      LNPWderiv.ch4.allvz.dy.dat

      common /snowmod/sza,vza,raz,refl,alb,dfalb,brdf  !snow 3.7um model
 

c *** read in corrk coeffs ***
      if(DataType .ne. 50)then
        open(2, file = FileNames(1), form='formatted')
      else
        open(2, file = FileNames(1))
      endif

      if(DataType .lt. 3 .or. DataType .eq. 50) then     
        num_seg = num_segment_3
      elseif (DataType .eq. 3 .or. DataType .eq. 6) then
        num_seg = i_modis_segment
      else
        num_seg = i_goes_segment
      endif

      if(DataType .ne. 50)then  !not mtsat1

      do iseg = 1, num_seg
        read(2,'(/,1x,i5)') n_h2o_3(iseg)                              ! Channel 3 water vapor
        read(2,'(1x,6f13.8)') (f_h2o_3(icoef,iseg),icoef=1,max_coefs)
        read(2,'(1x,6f13.8)') (k_h2o_3(icoef,iseg),icoef=1,max_coefs)
        read(2,'(1x,6e13.5)') (((coefk_h2o_3(icoef,jorder,kpres,iseg),
     +   icoef=1,max_coefs),jorder=1,iorder), kpres=1,num_pressure)
        if (DataType.ne.3.and.DataType.ne.6) then
          read(2,'(/,1x,i5)') n_ch4_3(iseg)                              ! Channel 3 methane
          read(2,'(1x,6f13.8)') (f_ch4_3(icoef,iseg),icoef=1,max_coefs)
          read(2,'(1x,6f13.8)') (k_ch4_3(icoef,iseg),icoef=1,max_coefs)
          read(2,'(1x,6e13.5)') (((coefk_ch4_3(icoef,jorder,kpres,iseg), 
     +      icoef=1,max_coefs),jorder=1,iorder), kpres=1,num_pressure)
          read(2,'(/,1x,i5)') n_n2o_3(iseg)                              ! Channel 3 N2O
          read(2,'(1x,6f13.8)') (f_n2o_3(icoef,iseg),icoef=1,max_coefs)
          read(2,'(1x,6f13.8)') (k_n2o_3(icoef,iseg),icoef=1,max_coefs)
          read(2,'(1x,6e13.5)') (((coefk_n2o_3(icoef,jorder,kpres,iseg), 
     +      icoef=1,max_coefs),jorder=1,iorder), kpres=1,num_pressure)
        end if
      end do
      read(2,'(/,1x,i5)') n_h2o_4                                  ! Channel 4 water vapor
      read(2,'(1x,6f13.8)') f_h2o_4
      read(2,'(1x,6f13.8)') k_h2o_4
      read(2,'(1x,6e13.5)') coefk_h2o_4
      if (DataType .lt. 3) then
        read(2,'(/,1x,i5)') n_co2_4                      ! Channel 4 carbon dioxide
        read(2,'(1x,6f13.8)') f_co2_4
        read(2,'(1x,6f13.8)') k_co2_4
        read(2,'(1x,6e13.5)') coefk_co2_4
      end if
      read(2,'(/,1x,i5)') n_h2o_5                         ! Channel 5 water vapor
      read(2,'(1x,6f13.8)') f_h2o_5
      read(2,'(1x,6f13.8)') k_h2o_5
      read(2,'(1x,6e13.5)') coefk_h2o_5
      close(2)

      else  !read mtsat1

c*** 3.75um channel        
        do iseg=1,num_seg
          do icoef=1,max_coefs
            f_h2o_3(icoef,iseg)=1.0
            k_h2o_3(icoef,iseg)=0.0
            f_ch4_3(icoef,iseg)=1.0
            k_ch4_3(icoef,iseg)=0.0
            f_n2o_3(icoef,iseg)=1.0
            k_n2o_3(icoef,iseg)=0.0
            do jorder=1,iorder
            do kpres=1,num_pressure
              coefk_h2o_3(icoef,jorder,kpres,iseg)=0.0
              coefk_ch4_3(icoef,jorder,kpres,iseg)=0.0
              coefk_n2o_3(icoef,jorder,kpres,iseg)=0.0
            enddo
            enddo
          enddo
          read(2,*)
          read(2,*)
          read(2,*)n_h2o_3(iseg)
          read(2,*)(f_h2o_3(icoef,iseg),icoef=1,n_h2o_3(iseg))
          read(2,*)(k_h2o_3(icoef,iseg),icoef=1,n_h2o_3(iseg))
          do icoef=1,n_h2o_3(iseg)
          do jorder=1,iorder
            read(2,*)(coefk_h2o_3(icoef,jorder,kpres,iseg),kpres=1,num_pressure)
          enddo
          enddo	  
          read(2,*)
          read(2,*)
          read(2,*)n_ch4_3(iseg)
          read(2,*)(f_ch4_3(icoef,iseg),icoef=1,n_ch4_3(iseg))
          read(2,*)(k_ch4_3(icoef,iseg),icoef=1,n_ch4_3(iseg))
          do icoef=1,n_ch4_3(iseg)
          do jorder=1,iorder
            read(2,*)(coefk_ch4_3(icoef,jorder,kpres,iseg),kpres=1,num_pressure)
          enddo
          enddo
          read(2,*)
          read(2,*)
          read(2,*)n_n2o_3(iseg)
	  if(n_n2o_3(iseg) .gt. 0)then
            read(2,*)(f_n2o_3(icoef,iseg),icoef=1,n_n2o_3(iseg))
            read(2,*)(k_n2o_3(icoef,iseg),icoef=1,n_n2o_3(iseg))
	  endif
          do icoef=1,n_n2o_3(iseg)
          do jorder=1,iorder
            read(2,*)(coefk_n2o_3(icoef,jorder,kpres,iseg),kpres=1,num_pressure)
          enddo
          enddo
          if(n_h2o_3(iseg) .eq. 0)n_h2o_3(iseg)=1
          if(n_ch4_3(iseg) .eq. 0)n_ch4_3(iseg)=1
          if(n_n2o_3(iseg) .eq. 0)n_n2o_3(iseg)=1
        enddo

c*** 11um channel
          do icoef=1,max_coefs
            f_h2o_4(icoef)=1.0
            k_h2o_4(icoef)=0.0
            f_co2_4(icoef)=1.0
            k_co2_4(icoef)=0.0            
            do jorder=1,iorder
            do kpres=1,num_pressure
              coefk_h2o_4(icoef,jorder,kpres)=0.0
              coefk_co2_4(icoef,jorder,kpres)=0.0
            enddo
            enddo
          enddo
          read(2,*)
          read(2,*)
          read(2,*)n_h2o_4
          read(2,*)(f_h2o_4(icoef),icoef=1,n_h2o_4)
          read(2,*)(k_h2o_4(icoef),icoef=1,n_h2o_4)
          do icoef=1,n_h2o_4
          do jorder=1,iorder
            read(2,*)(coefk_h2o_4(icoef,jorder,kpres),kpres=1,num_pressure)
          enddo
          enddo
          read(2,*)
          read(2,*)
          read(2,*)n_co2_4
          read(2,*)(f_co2_4(icoef),icoef=1,n_co2_4)
          read(2,*)(k_co2_4(icoef),icoef=1,n_co2_4)
          do icoef=1,n_co2_4
          do jorder=1,iorder
            read(2,*)(coefk_co2_4(icoef,jorder,kpres),kpres=1,num_pressure)
          enddo
          enddo
	  
c*** 12um channel
          do icoef=1,max_coefs
            f_h2o_5(icoef)=1.0
            k_h2o_5(icoef)=0.0
            f_co2_5(icoef)=1.0
            k_co2_5(icoef)=0.0            
            do jorder=1,iorder
            do kpres=1,num_pressure
              coefk_h2o_5(icoef,jorder,kpres)=0.0
              coefk_co2_5(icoef,jorder,kpres)=0.0
            enddo
            enddo
          enddo
	  
          read(2,*)
          read(2,*)
          read(2,*)n_h2o_5
          read(2,*)(f_h2o_5(icoef),icoef=1,n_h2o_5)
          read(2,*)(k_h2o_5(icoef),icoef=1,n_h2o_5)
          do icoef=1,n_h2o_5
          do jorder=1,iorder
            read(2,*)(coefk_h2o_5(icoef,jorder,kpres),kpres=1,num_pressure)
          enddo
          enddo
	 	  
          read(2,*)
          read(2,*)
          read(2,*)n_co2_5
          read(2,*)(f_co2_5(icoef),icoef=1,n_co2_5)
          read(2,*)(k_co2_5(icoef),icoef=1,n_co2_5)
          do icoef=1,n_co2_5
          do jorder=1,iorder
            read(2,*)(coefk_co2_5(icoef,jorder,kpres),kpres=1,num_pressure)
          enddo
          enddo
	  close(2)
      endif  !end reading ck data
      

c *** read in Tskin corrections for soundings ***
      open(2, file = FileNames(2), form='formatted')
      read(2,1000) dtoadskin
      close(2)

c *** read in precitable water corrections for soundings ***
      open(2, file = FileNames(3), form='formatted')
      read(2,1000) dtoadlnpw
      close(2)

c *** read in 3.7um snow models (from Bob Arduini)
        open(2,file=Filenames(4), form='unformatted')
        read(2) sza, vza, raz
        read(2) refl
        read(2) alb, dfalb
        read(2) brdf
        close(2)


 1000 format(30(f8.4,1x)) 
 
      return
      end

c*********************************************************************************
      subroutine snow_model37_ckm(vz,sz,rz,ref0,esnow,chi3x)
c*********************************************************************************
c ***** this subroutine provides bob arduini's reflectances
c       for snow/ice from solar zenith, viewing, azimuth angles
c       Provided by performing a three way interpolation.
c       This version is for Sep 2001 (higher res angles; uses
c       Bob's diff alb to get emiss)

      IMPLICIT NONE

c ***** variables
c       sz = solar zenith angle in degrees  0.=<sz<=90.      
c       vz = viewing zenith angle in degrees 0.=<vz<=90.
c       rz = relative azimuth angle in degrees, forward scat=0deg 0.=<rz<=180.
c       sza = solar zenith angle bin centers in degrees.
c       vza = viewing angle bin centers in degrees.
c       raz = relative azimuth angle bin centers in degrees.
c       d1,d2 = interpolated viewing angle factors
c       chi3x = bidir model over snow at 3.7um
c       esnow = snow emmitance; based on diffuse alb
c       ref0 = 3.7um snow reflectance
c       a() = interpolated azimuth angle factors
c       c() = interpolated solar angle factors    

      integer*4 iaz(2),ivz(2),isz(2)
      real*4 c(2),a(2),ref0,dtor
      save 
      data dtor / 0.01745329252 /  

      integer nsza,nvza,nraz,i,ksz,kaz
      parameter(nsza=21, nvza=21, nraz=24)

      real sza(nsza),vza(nvza),raz(nraz)
      real alb(nsza),dfalb
      real refl(nsza,nvza,nraz),brdf(nsza,nvza,nraz)
      real vz,sz,rz,esnow,chi3x,d1,d2,slop,cept

      common /snowmod/sza,vza,raz,refl,alb,dfalb,brdf  !3.7um snow model data

c  compute snow emittance 
      esnow=1.0-dfalb

      if(sz .ge. 90.0)then
        ref0=0.0
        chi3x=1.0
        return
      endif

c ***** set up angular bins
      do i=1,nraz-1
        if(rz .ge. raz(i) .and. rz .lt. raz(i+1))then     
          iaz(1) = i
          iaz(2) = i+1
        endif
      enddo
      if(rz .eq. 180.0)then
        iaz(1) = nraz
        iaz(2) = nraz
      endif

      do i=1,nvza-1
        if(vz .ge. vza(i) .and. vz .lt. vza(i+1))then     
          ivz(1) = i
          ivz(2) = i+1
        endif
      enddo
      
      if(vz .eq. 90.0)then
        ivz(1) = nvza
        ivz(2) = nvza
      endif

      do i=1,nsza-1
        if(sz .ge. sza(i) .and. sz .lt. sza(i+1))then     
          isz(1) = i
          isz(2) = i+1
        endif
      enddo
     
      if(sz .eq. 90.0)then
        isz(1) = nsza
        isz(2) = nsza
      endif

c ***** loop through for interp
      
         do 20 ksz = 1,2
            c(ksz) = 0.
            do 30 kaz = 1,2     
              d1 = refl(isz(ksz),ivz(1),iaz(kaz))
              d2 = refl(isz(ksz),ivz(2),iaz(kaz))          
               if(ivz(1) .ne. ivz(2)) then
                  call lineq_ckm(vza(ivz(1)),vza(ivz(2)),d1,d2,slop,cept)
                  a(kaz) = slop*vz + cept
               else
                  a(kaz) = d1
               endif
 30         continue
            if(iaz(1) .ne. iaz(2)) then
               call lineq_ckm(raz(iaz(1)),raz(iaz(2)),a(1),a(2),slop,cept)
               c(ksz) = slop*rz + cept
            else
               c(ksz) = a(1)
            endif
 20      continue     
          if(isz(1) .ne. isz(2)) then
             call lineq_ckm(sza(isz(1)),sza(isz(2)),c(1),c(2),slop,cept)
             ref0 = slop*sz + cept
          else
             ref0 = c(1)
          endif
        
          if(ref0.lt.0)then
            print *,'ref0<0 ',ref0,sz,vz,rz
          endif

c ***** loop through for interp
      
         do 40 ksz = 1,2
            c(ksz) = 0.
            do 50 kaz = 1,2     
              d1 = brdf(isz(ksz),ivz(1),iaz(kaz))
              d2 = brdf(isz(ksz),ivz(2),iaz(kaz))  
               if(ivz(1) .ne. ivz(2)) then
                  call lineq_ckm(vza(ivz(1)),vza(ivz(2)),d1,d2,slop,cept)
                  a(kaz) = slop*vz + cept
               else
                  a(kaz) = d1
               endif
 50         continue
            if(iaz(1) .ne. iaz(2)) then
               call lineq_ckm(raz(iaz(1)),raz(iaz(2)),a(1),a(2),slop,cept)
               c(ksz) = slop*rz + cept
            else
               c(ksz) = a(1)
            endif
 40      continue     
          if(isz(1) .ne. isz(2)) then
             call lineq_ckm(sza(isz(1)),sza(isz(2)),c(1),c(2),slop,cept)
             chi3x = slop*sz + cept
          else
             chi3x = c(1)
          endif
        
          if(chi3x.lt.0)then
            print *,'chi3x<0 ',chi3x,sz,vz,rz
          endif


 100  format(3f7.3)

      return
      end
c*****************************************************************************
      subroutine lineq_ckm(x1,x2,y1,y2,a,b) 

      IMPLICIT NONE

      real x1,x2,y1,y2,a,b

      if((x2-x1).eq.0) then
         print *,'lineq_ckm x1=x2',x1,x2,y1,y2
         a = 0.
         b = (y1+y2)*0.5
      else
         a = (y2-y1)/(x2-x1) 
         b = y1 - a * x1 
      endif

      return
      end 
c*****************************************************************************
      Subroutine AssignCommonBlockConst_ckm (DataType,pSolZenMax,pTwilit)
c    !************************************************************************
c    ! !f77
c    ! Name: Subroutine AssignCommonBlockConst_ckm
c    !
c    ! Purpose: assign several constants to common blocks :
c    !                   common /radiat/ c
c    !
c    ! !Input Parameters:
c    ! integer : DataType
c    !
c    ! !Output parameters:
c    ! none
c    !
c    ! Local Variables:
c    ! See declaration
c    !
c    ! Comments: 
c    !
c    ! !Revision History:
c    ! 06/25/98 : Initial delivery of software.
c    ! Sunny Sun-Mack - SAIC (s.sun-mack@larc.nasa.gov)
c    !
c    ! !END
c    !*************************************************************************

      include './corrk.h'
      implicit none

c     *** input parameters ***
      integer DataType
      real    pSolZenMax
      real    pTwilit

c     *** local parameters ***
      integer MaxChanCK
      parameter (MaxChanCK=4)
      integer i, num_pressure, num_segment_3, initsk, initpw
      parameter (num_segment_3 = 5)
      parameter (num_pressure = 19)
      real    wave(MaxChanCK), c(MaxChanCK,2), c2(num_segment_3,2)
      real    sc3temp, bsun, solzenmax, twilit
      real stp(num_pressure)
      integer iskin(30),ipw(20)

      integer nearIRchan
      integer*2   fil
      real        xfil
      real ramda_3(num_segment_3), delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)

      common /radiat/ c
      common /radiat2/ c2
      common /radiat3/ wave
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
      common /standpres/stp
      common /adjuststuff/iskin,ipw

      COMMON /nearIR/ nearIRchan
      common /fills/   xfil,fil
      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg
 
      fil = -9
      xfil = -9.0
c
c **  030500: ramda_3, delw_3 and response_3 will need to be changed for each satellite
c             However, for now they are hardwired to VIRS, all that is available
c             Also, bbnu_5seg and bsun_seg are defined later in fillatmos_ckm
c
c**delw is the wavelength interval
      if(DataType .eq. 50)then  !MTSAT
       ramda_3(1)=3.9526
       ramda_3(2)=3.8536
       ramda_3(3)=3.7488
       ramda_3(4)=3.6496 
       ramda_3(5)=3.5746
       delw_3(1)=0.0937
       delw_3(2)=0.1040
       delw_3(3)=0.1054
       delw_3(4)=0.0933
       delw_3(5)=0.0575
       response_3(1)=0.7990
       response_3(2)=0.9406
       response_3(3)=0.8760
       response_3(4)=0.7204
       response_3(5)=0.6197
      else  !VIRS, AVHRR
       ramda_3(1)=3.9604
       ramda_3(2)=3.8536
       ramda_3(3)=3.7488
       ramda_3(4)=3.6496 
       ramda_3(5)=3.5492
       delw_3(1)=0.1099
       delw_3(2)=0.1040
       delw_3(3)=0.1054
       delw_3(4)=0.0933
       delw_3(5)=0.1071
       response_3(1)=0.24603
       response_3(2)=0.97803
       response_3(3)=0.96178
       response_3(4)=0.96992
       response_3(5)=0.62487
      endif

c**   choose your instrument's near IR solar constant T in K units
c       347.00 = noaa10 and used for TOGA too
c       357.37 = noaa?, kansas, 1991 from dave (and apr94)
c       362.00 = noaa9 - don't know when
c       352.00 = noaa14
c       347.00 = noaa9 1986 global runs for 4/95 CERES STM
c       342.00 = noaa9 1986 global runs for 11/97 according to Minnis
c       359.64 = noaa9 1986 as of 11/97 according to Arduini
c       344.60 = GOES8 according to Arduini
c       355.40 = VIRS 1/98 according to Sunny
c       354.898 = VIRS 6/02 according to Sunny
c       353.21 = MODIS 10/00 according to Sunny
c       354.741 = MODIS 6/02 according to Sunny

      if ( DataType .eq. 0 ) then	! AVHRR/N12,14
        wave(1)= 3.75
        wave(2)=10.80
        wave(3)=12.00
        wave(4)=-9.00
        sc3temp=352.0              ! **** Check appropriate NOAA satellite ****
      elseif ( DataType .eq. 1 ) then	! AVHRR/N15-17
        wave(1)= 3.72
        wave(2)=10.80
        wave(3)=12.00
        wave(4)=-9.00
        sc3temp=360.07              ! **** Check appropriate NOAA satellite ****
      else if ( DataType .eq. 2 ) then  ! VIRS
        wave(1)= 3.78
        wave(2)=10.75
        wave(3)=11.95
        wave(4)=-9.00 
        sc3temp=354.898
      else if ( DataType .eq. 3 ) then  ! MODIS(Terra)
        wave(1)= 3.79
        wave(2)=11.03
        wave(3)=12.02
        wave(4)= 8.55
        sc3temp=354.81
      else if ( DataType .eq. 6 ) then  ! MODIS(Aqua)
        wave(1)= 3.79
        wave(2)=11.03
        wave(3)=12.02
        wave(4)= 8.55
        sc3temp=355.13
      else if ( DataType .eq. 11 ) then  ! use GOES10 for TSA
        wave(1)= 3.90
        wave(2)=10.80
        wave(3)=11.90
        wave(4)=-9.00 
        sc3temp=354.741
      else if ( DataType .eq. 21 ) then  ! GOES7
        wave(1)= 3.94
        wave(2)=11.20
        wave(3)=12.66
        wave(4)=-9.00 
        sc3temp=344.60
      else if ( DataType .eq. 22 ) then  ! GOES8
        wave(1)= 3.90
        wave(2)=10.80
        wave(3)=11.90
        wave(4)=-9.00 
        sc3temp=344.60
      else if (DataType .eq. 23) then  ! GOES10
        wave(1)= 3.92
        wave(2)=10.80
        wave(3)=11.90
        wave(4)=-9.00 
        sc3temp=344.16
      else if (DataType .eq. 26) then  ! GOES11
        wave(1)= 3.91
        wave(2)=10.80
        wave(3)=11.90
        wave(4)=-9.00 
        sc3temp=344.89
      else if ( DataType .eq. 24 ) then ! GOES12 
        wave(1)= 3.90
        wave(2)=10.80
        wave(3)=13.30   !added 12/2005
        wave(4)=-9.00 
        sc3temp=345.11
      else if ( DataType .eq. 25 ) then ! GOES09 
        wave(1)= 3.91
        wave(2)=10.80
        wave(3)=11.90
        wave(4)=-9.00 
        sc3temp=344.47
      else if ( DataType .eq. 40 ) then ! MSG1
        wave(1)= 3.92
        wave(2)=10.80
        wave(3)=12.00
        wave(4)=-9.00 
        sc3temp=343.95
      else if ( DataType .eq. 50 ) then ! MTSAT-1R  !added 12/2005
        wave(1)= 3.75
        wave(2)=10.80
        wave(3)=12.00
        wave(4)=-9.00 
        sc3temp=355.4
      else
        print *,'Corrk Error: Unrecognized DataType...Exiting'
        stop
      endif

c**   nearIRchan is used in Init_Cloud as well as simulite and other subroutines 
c     related to the emittance parameterizations 

      if (DataType.le.6) then
        nearIRchan = 3                                 ! AVHRR, MODIS or VIRS channel
      elseif(DataType .eq. 50)then
        nearIRchan = 5
      else
        nearIRchan = 2                                 ! GOES channel
      end if
      
     
      solzenmax = pSolZenMax

      twilit = pTwilit
c
c Calculate Planck constants
c
      do i = 1,MaxChanCK
        c(i,1)=1.191044e8/(wave(i)*wave(i)*wave(i)*wave(i)*wave(i))
        c(i,2)=1.438769e4/wave(i)
      enddo
c
c Calculate Planck constants for 3.7um subinterval wavelengths - used in fillatmos_ckm
c
      do i = 1,num_segment_3
        c2(i,1)=1.191044e8/(ramda_3(i)*ramda_3(i)*ramda_3(i)*ramda_3(i)*ramda_3(i))
        c2(i,2)=1.438769e4/ramda_3(i)
      enddo
c
c**   begin assignments for VIST subroutine adjust and associated subroutines
c
      stp(1)= 0.251
      stp(2)= 0.398
      stp(3)= 0.631
      stp(4)= 1.000
      stp(5)= 1.58
      stp(6)= 2.51 
      stp(7)= 3.98
      stp(8)= 6.31
      stp(9)= 10.0
      stp(10)= 15.8
      stp(11)= 25.1
      stp(12)= 39.8
      stp(13)= 63.1
      stp(14)= 100.0
      stp(15)= 158.0
      stp(16)= 251.0
      stp(17)= 398.0
      stp(18)= 631.0
      stp(19)= 1000.0

      initsk=195
      do i=1,30
        initsk=initsk+5
        iskin(i)=initsk
      enddo

c**   pw range=0.04-6.327; scale by 1000
      ipw(1)=40
      initpw=0
      do i=2,20
        initpw=initpw+333
        ipw(i)=initpw
      enddo      
c
c**   end assignments for VIST subroutine adjust and associated subroutines
c
      return
      end 
      

c*****************************************************************************
      Subroutine AssignCommonBlocks_ckm (    MOANumLevels,  
     &                    MOAPressure,   MOATemperature,
     &                    MOAHeight,     MOARHumidity,  MOAOzone,                   
     &                    MeanRelHum,       MeanRadiance,
     &                    SurfaceRadiance,  OpticalDepth,
     &                    AttenuationRatio, TotalReflectRat,
     &                    SolarRadiance,    RatioGoingUp)

c    !*************************************************************************
c    ! !f77
c    ! Name: Subroutine AssignCommonBlocks_ckm
c    !
c    ! Purpose: Put CorrK results to common blocks :
c    !             COMMON /ktemps/ mlv,p0,t0,z0,u0,rh0,o30,i85flag
c    !             COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4
c    !             COMMON /krads/ bnu,bnu0
c    !             COMMON /ktaus/ attenuate,ratio
c    !             COMMON /constants/ sc3temp,bsun,solzenmax,twilit
c    !
c    ! !Input Parameters:
c    ! See declaration
c    !
c    ! !Output parameters:
c    ! See declaration
c    !
c    ! Local Variables:
c    ! See declaration
c    !
c    ! Comments: 
c    !
c    ! !Revision History:
c    ! 10/20/97 : Initial delivery of software.
c    ! Sunny Sun-Mack - SAIC (s.sun-mack@larc.nasa.gov)
c    !
c    ! !END
c    !*************************************************************************

      include './corrk.h'
      implicit none

      integer MaxChanCK

      PARAMETER (MaxChanCK = 4)

c  *** input parameters ***
      integer  MOANumLevels
      real     MOAPressure     (MAXCKLEV)
      real     MOATemperature  (MAXCKLEV)
      real     MOAHeight       (MAXCKLEV)
      real     MOARHumidity    (MAXCKLEV)
      real     MOAOzone        (MAXCKLEV)
      real     MeanRelHum      (MAXCKLEV)
      real     MeanRadiance    (MAXCKLEV, MaxChanCK)
      real     SurfaceRadiance (MaxChanCK)
      real     OpticalDepth    (MAXCKLEV, MaxChanCK)
      real     AttenuationRatio(MAXCKLEV)
      real     TotalReflectRat
      real     SolarRadiance   (MAXCKLEV)
      real     RatioGoingUp    (MAXCKLEV)


c  *** local parameter ***
      integer  i, j, l, ic
      integer  mlv, i85flag
      real     sc3temp, bsun, solzenmax, twilit

      real p0(MAXCKLEV),t0(MAXCKLEV),u0(MAXCKLEV)
      real z0(MAXCKLEV),rh0(MAXCKLEV),o30(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV),o3lay(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real attenuate(MAXCKLEV,MaxChanCK),ratio(MAXCKLEV),SolarRad(MAXCKLEV)
      real TotalRefRat, ratio_up(MAXCKLEV)
          
      real absorb_bnu(MAXCKLEV,MaxChanCK),absorb_totaltau(MAXCKLEV,MaxChanCK)
      real interCoef(MAXCKLEV,MaxChanCK)                  
      real absorb_slope(MAXCKLEV,MaxChanCK,MAXCKLEV)       
      real absorb_intercept(MAXCKLEV,MaxChanCK,MAXCKLEV)   
      real xslope, xintercept                  
      real planck_ckm, planckiold_ckm

      COMMON /absorb_up/ absorb_slope,absorb_intercept 
      COMMON /ktemps/    mlv,p0,t0,z0,u0,rh0,o30,i85flag
      COMMON /kprofs/    play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/     bnu,bnu0
      COMMON /ktaus/     attenuate,ratio, TotalRefRat, SolarRad, ratio_up
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit

c *** assign common block /ktemps/ ***
      mlv = MOANumLevels
      do i = 1, MOANumLevels
         p0(i)  = MOAPressure    (i)
         t0(i)  = MOATemperature (i)
         z0(i)  = MOAHeight      (i)
         rh0(i) = MOARHumidity   (i)
         o30(i) = MOAOzone       (i)
      end do


c *** assign common block /kprofs/ ***
      do i = 1, MOANumLevels-1
         ulay(i)    = MeanRelHum    (i)
      end do


c *** assign common block /krads/ ***
      do j = 1, MaxChanCK
         do i = 1, MOANumLevels-1
            bnu(i,j) = MeanRadiance(i, j) 
         end do
         bnu0(j) = SurfaceRadiance (j)
      end do


c *** assign common block /ktaus/ ***
      do j = 1, MaxChanCK
         do i = 1, MOANumLevels-1
            attenuate(i,j) = OpticalDepth (i,j) 
         end do
      end do

      do i = 1, MOANumLevels-1
         ratio(i)    = AttenuationRatio (i) 
         SolarRad(i) = SolarRadiance(i)
         ratio_up(i) = RatioGoingup(i)
      end do
      
      TotalRefRat = TotalReflectRat

c *** convert bnu to Garber's bnu 
c *** absorb_num and absorb_totaltau

      do i=1,MOANumLevels-1
        do j=1,MaxChanCK
c         absorb_bnu(i,j) = planck_ckm(j,planckiold_ckm(j,bnu(i,j)))
          absorb_bnu(i,j) = bnu(i,j)
          absorb_totaltau(i,j)=attenuate(i,j)
        end do
      end do


c *** calculate/assign common block /absorb_up/ ***
      do i = 1,MaxChanCK
         do j = 1,MOANumLevels-1
            interCoef(j,i)=absorb_bnu(j,i)*(1.0-absorb_totaltau(j,i))
         end do
      end do

      do l = 1, MOANumLevels-1
         do i = 1, MaxChanCK
            xSlope     = 1.0
            xIntercept = 0.0
            do j = l, MOANumLevels-1
               xSlope     = xSlope * absorb_totaltau(j,i)
               xIntercept = xIntercept * absorb_totaltau(j,i) 
     &                      + interCoef (j,i)
               absorb_slope     (j,i,l) = xSlope
               absorb_intercept (j,i,l) = xIntercept
            end do
         end do
      end do

      return
      end 
c****************************************************************************
      subroutine GetCorrK_ckm(DataType,    iupdown,     dis,             !inputs 
     &                    vz,               sz,             SBT,
     &                    emit,             MOANumLevels,
     &                    MOAPressure,      MOATemperature,
     &                    MOAHeight,        MOARHumidity,   MOAOzone,
     &                    MeanRelHum,       MeanRadiance,              !outputs
     &                    SurfaceRadiance,  OpticalDepth, 
     &                    AttenuationRatio, TotalReflectRat,
     &                    SolarRadiance,    RatioGoingUp)

c    !*************************************************************************
c    ! !f77
c    ! Name: Subroutine GetCorrK_ckm
c    !
c    ! Purpose: Calculate CorrK and return CorrK information
c    !
c    ! !Input Parameters:
c    ! See declaration
c    !
c    ! !Output parameters:
c    ! See declaration
c    !
c    ! Local Variables:
c    ! See declaration
c    !
c    ! Comments:
c    !
c    ! !Revision History:
c    ! 10/20/97 : Initial delivery of software.
c    ! Sunny Sun-Mack - SAIC (s.sun-mack@larc.nasa.gov)
c    !
c    ! !END
c    !*************************************************************************

      include './corrk.h'
      implicit none

      integer MaxChanCK
      PARAMETER (MaxChanCK = 4)

c  *** input parameters ***
      real     vz
      real     sz
      real     emit(MAXCHAN)

      integer  MOANumLevels
      real     MOAPressure     (MAXCKLEV)
      real     MOATemperature  (MAXCKLEV)
      real     MOAHeight       (MAXCKLEV)
      real     MOARHumidity    (MAXCKLEV)
      real     MOAOzone        (MAXCKLEV)


c   *** output parameters ***
      real     MeanRelHum      (MAXCKLEV)
      real     MeanRadiance    (MAXCKLEV, MaxChanCK)
      real     SurfaceRadiance (MaxChanCK)
      real     OpticalDepth    (MAXCKLEV, MaxChanCK)
      real     AttenuationRatio(MAXCKLEV)
      real     TotalReflectRat
      real     SolarRadiance   (MAXCKLEV)
      real     RatioGoingUp   (MAXCKLEV)


c  *** local parameter ***
      integer  i, j
      integer  mlv, i85flag
      real     dis, SBT, SkinT, precip_water
      integer  DataType, iupdown

      real p0(MAXCKLEV),t0(MAXCKLEV),u0(MAXCKLEV)
      real z0(MAXCKLEV),rh0(MAXCKLEV),o30(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV),bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)

      real attenuate(MAXCKLEV,MaxChanCK),ratio(MAXCKLEV),SolarRad(MAXCKLEV)
      real TotalRefRat,ratio_up(MAXCKLEV)

      COMMON /ktemps/ mlv,p0,t0,z0,u0,rh0,o30,i85flag

      COMMON /kprofs/play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0
      COMMON /ktaus/ attenuate,ratio, TotalRefRat, SolarRad, ratio_up
      COMMON /SkinTkeep/ SkinT, precip_water
c
c       The following routines are called whenever the sounding is changed.
c       'fillatmos_ckm' calculates layer mean T and RH from the profile
c       'CorrKsetup_ckm' calls the correlated K routines from Kratz and sets up
c                    several arrays containing layer mean transmission info
c                    that are used to calculate Tskin and surface emittance


      mlv = MOANumLevels
      do i = 1, MOANumLevels
         p0(i) = MOAPressure (i)
         t0(i) = MOATemperature (i)
         z0(i) = MOAHeight(i)
         rh0(i) = MOARHumidity(i)
         o30(i) = MOAOzone(i)
      end do


      if(iupdown .eq. 1)then
        SkinT=SBT
      else
        SkinT=0.0
      endif

      call fillatmos_ckm

      if(DataType .le. 2) then
        call CorrKsetup_ckm(vz,sz,dis,emit(1),emit(2),emit(3))
c	print*,'Called CorrK setup for AVHRR,VIRS'
      elseif (DataType .eq. 3 .or. DataType .eq. 6) then
        call CorrKsetup_modis_ckm(vz,sz,dis,emit(1),emit(2),emit(3))
c	print*,'Called CorrK setup for MODIS'
      elseif(DataType .eq. 50)then
        call CorrKsetup_mtsat_ckm(vz,sz,dis,emit(1),emit(2),emit(3))
c	print*,'Called CorrK setup for MTSAT'
      else
        call CorrKsetup_goes_ckm(vz,sz,dis,emit(1),emit(2),emit(3))
c	print*,'Called CorrK setup for GOES'
      end if

      do i = 1, MOANumLevels-1
         MeanRelHum     (i) = ulay(i)
c        MeanPressure   (i) = play(i)
c        MeanTemperature(i) = tlay(i)
c        DeltaPressure  (i) = delplay(i)
c        DeltaHeight    (i) = delzlay(i)
c        CO2Density     (i) = uco2(i)
c        CH4Density     (i) = uch4(i)
c        CFCDensity     (i) = ucfc(i)
c        N2ODensity     (i) = un2o(i)
      end do


      do j = 1, MaxChanCK
         do i = 1, MOANumLevels-1
            MeanRadiance(i, j) = bnu(i,j)
         end do
         SurfaceRadiance (j)   = bnu0(j)
      end do


      do j = 1, MaxChanCK
         do i = 1, MOANumLevels-1
            OpticalDepth (i,j) = attenuate(i,j)
         end do
      end do

      do i = 1, MOANumLevels-1
         AttenuationRatio(i) = ratio(i)
         SolarRadiance(i)    = SolarRad(i)
         RatioGoingUp(i)     = ratio_up(i)
      end do

      TotalReflectRat = TotalRefRat

c     SConstInRadCh3 = bsun

      return
      end
c*****************************************************************************
      subroutine SetCorrkReg_ckm(DataType,  iupdown, discf,                
     &                    vza,           sza,         SBT,
     &                    emit,          NumLevels,
     &                    Pressure,      Temperature,
     &                    Height,        RHumidity,   Ozone)
c*******************************************************************************

c     This subroutine calls GetCorrK_ckm and AssignCommonBlocks_ckm, it should be
c     called for each region or profile.

      include './corrk.h'
      IMPLICIT NONE

      integer DataType,NumLevels,MaxChanCK,iupdown

      parameter (MaxChanCK=4)

      real    Height(MAXCKLEV), Pressure(MAXCKLEV),
     &        Temperature(MAXCKLEV), RHumidity(MAXCKLEV), Ozone(MAXCKLEV)
      real     MeanRelHum      (MAXCKLEV)
      real     MeanRadiance    (MAXCKLEV, MaxChanCK)
      real     SurfaceRadiance (MaxChanCK)
      real     OpticalDepth    (MAXCKLEV, MaxChanCK)
      real     AttenuationRatio(MAXCKLEV)
      real     TotalReflectRat
      real     SolarRadiance   (MAXCKLEV)
      real     RatioGoingUp   (MAXCKLEV)
      real  emit(MAXCHAN)
      real  SBT,discf,vza,sza

       call GetCorrK_ckm(DataType,   iupdown,     discf,                  !inputs 
     &                vza,            sza,           SBT,
     &                emit,           NumLevels,
     &                Pressure,         Temperature,
     &                Height,           RHumidity,     Ozone,
     &                MeanRelHum,       MeanRadiance,          !outputs
     &                SurfaceRadiance,  OpticalDepth, 
     &                AttenuationRatio, TotalReflectRat,
     &                SolarRadiance,    RatioGoingUp)


       call AssignCommonBlocks_ckm (NumLevels,  
     &                    Pressure,   Temperature,
     &                    Height,     RHumidity,  Ozone,                   
     &                    MeanRelHum,       MeanRadiance,
     &                    SurfaceRadiance,  OpticalDepth,
     &                    AttenuationRatio, TotalReflectRat,
     &                    SolarRadiance,    RatioGoingUp)


      return
      end
c*****************************************************************************
      subroutine GetBTCorrk_ckm(DataType, iupdown, discf,                 
     &                    vza,  sza,   aza,  SBT, iradio, waved,
     &                    emit, emsnow, chiopen, NumLevels, nch,
     &                    Pressure,      Temperature,
     &                    Height,        RHumidity,   Ozone,
     &                    snowflg, icep, wave, FBT, SRef, FRef)
c*******************************************************************************

c     This subroutine calls SetCorrkReg_ckm and GetBTReg_ckm to obtain T3.7, T11, T12,
c     R3.7 at TOA.  It should be called for each region or profile.
c**INPUTS:
c      DataType = 0,1=AVHRR, 2=VIRS, 3=TERRA, 6=AQUA, 11=G10TSA, 
c                 21=G07, 22=G08, 23=G10, 24=G12, 25=G09, 26=G11, 40=MSG1, 50=MTSAT-1R
c      discf = earth-sun distance correction factor
c      vza = view zenith angle
c      sza = solar zenith angle
c      aza = relative azimuth angle, forward scatter=0deg
c      SBT = starting brightness temperature at sfc or TOA
c      emit = surface emittance for each wavelength (non-snow)
c      emsnow = The snow/ice emittance for each wavelength
c      NumLevels = #levels for corr-k (min of 8); level #1 is at the surface
c      nch = number of channels to process (1-4)
c      Pressure = pressure profile (mb)
c      Temperature = temperature profile (K)
c      Height = Height profile (m)
c      RHumidity = relative humidity profile (%)
c      Ozone = Ozone profile (g/g)
c      wave = wavelengths for TOA BT (3.7,11,12um)
c      chiopen = bidirectional refl factor for near-ir solar part for non-snow.
c      snowflg=snow flag from CERES maps (-1, 0=no snow, or 1=snow)
c      icep=ice percent from CERES maps (-1, 0-100)
c      SRef=Starting 3.0-4.0um reflectance. Set to -1 to use BobA snow model.
c      iupdown 1=solve upward, 2=solve downward
c      iradio=radiometric tsfc from toa; emit=1.0 (1) or true Tskin from TOA (2)
c      waved=wavelength to process, for downward direction only.
c      
c**Outputs:
c      FBT = Top-of-atmosphere or surface brightness temperatures for the input wavelengths
c      FRef = Top-of-atmosphere or surface Reflectance of near-ir band


      include './corrk.h'
      IMPLICIT NONE

      integer DataType,NumLevels,MaxChanCK,nch,i,ich,snowflg,iupdown
      integer iradio

      parameter (MaxChanCK=4)

      real  Height(MAXCKLEV),Height_km(MAXCKLEV),Pressure(MAXCKLEV),
     &      Temperature(MAXCKLEV), RHumidity(MAXCKLEV), Ozone(MAXCKLEV)
      real  emit(MAXCHAN),wave(MaxChanCK),emsnow(MaxChanCK)
      real  SBT,discf,FBT(MaxChanCK),chi3sno,refsfcck,waved
      real  vza,sza,aza,chi,chiopen,SRef,SRefm,FRef,icep,snoicep

      
      do i=1,NumLevels
        Height_km(i)=Height(i)/1000.0
      enddo


c***combine snow and ice maps into one
      snoicep=icep/100.0
      if(snowflg .eq. 1)then
        snoicep=1.0
      endif

      SRefm=-9.0

      call snow_model37_ckm(vza,sza,aza,SRefm,emsnow(1),chi3sno)

c      print*,'snow 37 emit:',emsnow(1)

      
      if(snoicep .eq. 0.0)then
        chi=chiopen
      else      
        do i=1,nch
          emit(i)=(emsnow(i)*snoicep) + (emit(i)*(1-snoicep))
        enddo
        chi=(chi3sno*snoicep) + (chiopen*(1-snoicep))
      endif

      if(SRef .le. 0.0 .and. iupdown .eq. 1)SRef=SRefm  !take from BobA snow model.
      
      
      call SetCorrkReg_ckm(DataType, iupdown, discf,               
     &                 vza, sza, SBT,
     &                 emit,NumLevels,
     &                 Pressure, Temperature,
     &                 Height_km, RHumidity, Ozone)


c***Get BT 3.7,11,12um, Ref3.7
      do ich=1,nch
        if(iupdown .eq. 1 .or. (iupdown .eq. 2 .and. waved .eq. wave(ich)))then
          call MoveBTReg_ckm(iupdown,iradio,SBT,SRef,NumLevels,sza,wave(ich),
     &                     emit(ich),chi,FBT(ich),FRef)
          if(FBT(ich) .ge. 100.0 .and. FBT(ich) .le. 350.0)then
          else  !bounds or nan error
            write(47,100) (Height_km(i), Pressure(i),
     &       Temperature(i), RHumidity(i),i=1,NumLevels)
           write(47,*)ich,Height_km(1),SBT,FBT(ich),NumLevels
           write(47,*)
           print 101,'Corrk error:ich,SBT,FBT(1-3),sza,vza:',ich,SBT,FBT(ich),sza,vza
           FBT(ich)=-99.0
c           stop
          endif
        else
          FBT(ich)=-99.0
          FRef=-0.999
        endif
      enddo

100   format(f9.4,3f8.2)
101   format(a35,i4,4f8.2)
      return
      end
c******************************************************************************  
      Subroutine fillatmos_ckm
c
c	This subroutine calculates mean layer values of:  temperature (K)
c                                                     pressure (atm)
c                                                     specific humidity (g/cm^2)
c                                                     CO2 (g/cm^2)
c                                                     CH4 (g/cm^2)
c                                                     N2O (g/cm^2)
c                                                     CFC-F12 & CFC-F11 (g/cm^2)
c
c	These are needed to use Dave Kratz's correlated K's
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      parameter (num_segment_3 = 5)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV),bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real ramda_3(num_segment_3),delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0

      COMMON /SkinTkeep/ SkinT, precip_water

      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg
     
      solar_temp = 5710.0

      sfcP = 1013.25	! This is used to convert millibars to atmospheres

      precip_water=0.0

c** uf12, uf11 are the cfc's
      do m = 1,klevs - 1
        tlay(m) = (tk(m) + tk(m+1)) / 2.0
        play(m) = (pk(m) + pk(m+1)) / 2.0 / sfcP
        delplay(m) = (pk(m) - pk(m+1))
	delzlay(m) = zk(m+1) - zk(m)
        zlay=(zk(m+1)+zk(m))/2.0
        uco2(m) = play(m) * 1.01325e+06 /
     *          (1.3805e-16 * tlay(m)) * 3.50e-04 * delzlay(m) *
     *          1.0e+05 / 6.023e+23 * 44.00995
        uch4(m) = play(m) * 1.01325e+06 /
     *          (1.3805e-16 * tlay(m)) * 1.75e-06 * delzlay(m) *
     *          1.0e+05 / 6.023e+23 * 16.04303
        un2o(m) = play(m) * 1.01325e+06 / (1.3805e-16 * tlay(m)) *
     *          2.80e-07 * delzlay(m) * 1.0e+05 / 6.023e+23 * 44.0128     
        uf12(m)=play(m)*1.01325e+06/
     &         (1.3805e-16*tlay(m))*0.503e-09*delzlay(m)*
     &         1.0e+05/6.023e+23*120.9139
        uf11(m)=play(m)*1.01325e+06/
     &         (1.3805e-16*tlay(m))*0.268e-09*delzlay(m)*
     &         1.0e+05/6.023e+23*137.3685

c       print*,m+1,zk(m+1),uf11(m)

       if(zk(m+1) .gt. 25.0)then
         uf12(m)=uf12(m)*exp((-7.5-(zk(m+1)-25.0))/5.0)
         uf11(m)=uf11(m)*exp((-7.5-(zk(m+1)-25.0))/3.0)
       elseif(zk(m+1) .ge. 15.5)then
         uf12(m)=uf12(m)*exp((15.5-zk(m+1))/5.0)
         uf11(m)=uf11(m)*exp((15.5-zk(m+1))/3.0)
       endif

c       print*,m+1,zk(m+1),uf11(m)
 
        bnu(m,1) = planck_ckm(1,tlay(m))
        bnu(m,2) = planck_ckm(2,tlay(m))
        bnu(m,3) = planck_ckm(3,tlay(m))
        bnu(m,4) = planck_ckm(4,tlay(m))

        do iseg = 1, num_segment_3
          bnu_5seg(m,iseg) = planck2_ckm(iseg,tlay(m))
        end do

        if(uk(m) .gt. 0.) then 
          ulay(m) = delzlay(m) * (uk(m) + uk(m+1)) / 2.0
        else if(rhk(m) .gt. 0.)  then
c         qave = (q0(m) + q0(m+1)) / 2.0
          rhave = (rhk(m) + rhk(m+1)) / 2.0
          qave = SPEC(tlay(m),play(m)*sfcP,rhave)
          ulay(m) = 1.02 * delplay(m) * qave  !PW=(1/g)*dP*q; fact of 10 from mb-Pa, kg/m2 to g/cm2
c       else if(q0(m) .gt. 0.)  then
c         qave = (q0(m) + q0(m+1)) / 2.0
c         rhave = RELHUM(tlay(m),play(m)*sfcP,qave)
c         ulay(m) = 1.02 * delplay(m) * qave
        else
          print*,'error with moisture sounding',m,rhk(m)
          stop2
        end if
        precip_water=precip_water+ulay(m)


***compute ozone for each layer by changing level 
c        values (g/g) to Kratz world (g/(cm*cm)/km)
        if(i85flag .eq. 1)then
        if (o3k(m) .gt. 0.) then
          o3ave = (o3k(m) + o3k(m+1)) / 2.0
          o3lay(m) = 1.02 * delplay(m) * o3ave
        else
          o3lay(m)=0.0
        end if
        endif

      end do

c     do m = 1,klevs - 1
c       if (tlay(m).lt.tlay(m+1)) then
c          mmmm=1
c       end if
c     end do
      
      if (SkinT.eq.0.0) then
        SkinT = tk(1)		  ! Use surface T if no Skin T
      end if
	  
      bnu0(1) = planck_ckm(1,SkinT)	  ! Define surface radiance using SkinT
      bnu0(2) = planck_ckm(2,SkinT)
      bnu0(3) = planck_ckm(3,SkinT)
      bnu0(4) = planck_ckm(4,SkinT)

c     sigma=5.67e-05
c     sconst=1.368e+06
c     fsun=sigma*solar_temp**4

      do iseg = 1, num_segment_3
        bsun_seg(iseg) = planck2_ckm(iseg,solar_temp)
c       bsun_seg(iseg) = bsun_seg(iseg)/fsun*sconst
      end do
  
      return
      end
c********************************************************************************	  
      Subroutine FindCloudLayer_ckm(Tcloud,icloudlayer)
c
c	This subroutine finds the layer in which the "cloud temperature" resides.
c	This is done for purposes of applying the correlated Ks.
c
c	Our basic assumption is that the the cloud is infinitesimally thin at the level
c	of the cloud temperature.  We also assume that all of the humidity in that layer
c	is either inside or below the cloud.  Therefore, we include the cloud layer in
c	gaseous absorption calculations from below, but stop at the next layer above the
c	cloud when making calculations above the cloud.
c
c	In general, when when calling UseCorrKRatio_ckm, you should use either:
c			laylow = icloudlayer + 1
c		or	layhigh = icloudlayer
c
c	A value of Tcloud = 0. will return icloudlayer = top layer of atmosphere (klevs - 1)
c
c*************************************************************************
c        Begin                CHANGE  5/2000
c*************************************************************************
c
c   It has been decided by the powers that be (Pat Minnis) that we should include the
c   absorption from the layer that the cloud resides in the calculations ABOVE the cloud.
c   This should result in slightly reduced particle sizes & account for molecular
c   absorption within the cloud.  So we are now assuming that the cloud is at the bottom
c   of the layer, instead of the top.
c
c   What this means is that in general when we call UseCorrK_ckm:
c     laylow = icloudlayer
c     layhigh = icloudlayer - 1
c   instead of the convention described above.
c
c   However, it would be a pain to change in all the places where we call UseCorrK_ckm,
c   so the easier solution is to change the definition of icloudlayer so that it's
c   value is decreased by 1.
c
c   This change has been made to the final executable line of this subroutine.  Since
c   laylow & layhigh are compared at the start of UseCorrK_ckm, the chance that we end up
c   with a cloudlayer = 0 will not cause problems.
c
c*************************************************************************
c        End                  CHANGE  5/2000
c*************************************************************************
      
      include './corrk.h'
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      REAL Tcloud
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      
      if(Tcloud .gt. 0.) then
        mlev = 2
        do while(Tcloud .lt. tk(mlev) .and. mlev .le. klevs - 1)
          mlev = mlev + 1
        end do
        icloudlayer = mlev - 1
        icloudlayer = icloudlayer  - 1
      else
        icloudlayer = klevs - 1
      end if
  
c     icloudlayer = icloudlayer  - 1                  !correction 10/2000

      return
      end
c************************************************************************************	  
      Subroutine CorrKsetup_ckm(vz,sz,dis,surfemit3,surfemit4,surfemit5)
c
c	This subroutine sets up the correlated Ks.
c
c	Kratz's code is used to calculate upward flux at each level.  I then solve for
c	  a total optical depth for each layer.  These are then used in routine "UseCorrK_ckm".
c     These optical depths are only valid for the given angle.
c
c   The major changes in this version:
c      1)  The surface emittances are used (these are now input parameters)
c      2)  Channel 3 surface reflectance = 1. - surfemit3 instead of a constant value of 0.10
c      3)  For channel 3, we now solve downward first.  We then use the solar radiance at
c          the surface for bref, instead of using the radiance at the TOA
c      4)  We now use subroutine ck for all correlated k calculations (except for those
c          that involve CFCs).  The inputs to subroutine ck have changed to include the 
c          coefficients.  These coefficients are now read from an ancillary data file during
c          the first call to CorrKsetup_ckm
c      5)  Kratz's latest edition of correlated ks for Channel 3 has the channel divided into 5 
c          spectral intervals.  Comparisons with the old version shows only a 1-2% improvement.
c          I've decided to stick with the old version.  The new version would have required 
c          separate bnu arrays for each spectral interval & the run-time of the Channel 3 
c          portion of the code would have increased by a factor of 5.
c      6)  Kratz has added absorption by CFCs to channels 4 & 5.  This changes the COMMON 
c          block /kprofs/ which now includes the array "ucfc".  This array is filled in 
c          subroutine fillatmos_ckm.  The correlated k's are calculated in subroutines,
c          ck_cfc_4_ckm & ck_cfc_5_ckm.
c      7)  The window is now modeled using the Clough continuum, instead of Roberts.  This
c          uses a parameterization developed by Fred Rose.  The old Subroutine "window" has been
c          replaced by two channel-specific routines, window_4_ckm & window_5_ckm.
c      8)  We no longer calculate the total optical depth (totaltau).  This was an unnecessary 
c          and time-consuming step that I did originally to check if the numbers looked correct.  
c          This led to some of the problems that we had with NAN results when we tried to take the
c          logarithm of 0. or a negative number.
c          We now just calculate the attenuation.  This causes a change in Subroutines UseCorrK_ckm  & GetSurtTemp
c          since we will no longer need to take the EXP of the optical depth.  It also changes the common
c          block /ktaus/.
c      9)  The layer indexing in the code was changed slightly from Kratz's original.
c          Now, level 1 = surface, level 2 = the first level above the surface, etc. up to level klevs
c          And, layer 1 = the layer between levels 1 & 2, etc. up to layer (klevs - 1) = nlay
c          Kratz's old numbering system was different in the upward and downward directions.
c     10)  I've changed the cos(sz) & cos(vz) into sec(sz) & sec(vz) 
c     11)  I now check to see if we are dividing by 0. when calculating the layer attenuation.
c          This should eliminate NAN problems when there is a temperature inversion.
c  
c
c	It is assumed that Subroutine fillatmos_ckm has been called to fill in values of bnu for
c    each atmospheric layer and bnu0 for the surface radiance.
c
c  INPUT PARAMETERS:	
c		vz          = viewing zenith
c		sz          = solar zenith angle
c       surfemit3   = VIRS Channel 3 (3.75 µm) surface emittance
c       surfemit4   = VIRS Channel 4 (10.8 µm) surface emittance
c       surfemit5   = VIRS Channel 5 (11.9 µm) surface emittance
c       bnu(ilay)   = emitted radiance from layer ilay
c       bnu0(ichan) = emitted radiance for the surface for channel ichan (should be planck(SkinTemp))
c                     where ichan = 1 for Channel 3
c                           ichan = 2 for Channel 4
c                           ichan = 3 for Channel 5
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      parameter (num_segment_3 = 5)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV),bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real dis

      real sx4(MAXCKLEV),tsx(MAXCKLEV)
c     real sx3(MAXCKLEV)
      real sxz3(MAXCKLEV)
      real tsxb(MAXCKLEV)
c     real tsxa(MAXCKLEV)
      real tsyb(MAXCKLEV),tsya(MAXCKLEV)
      real sy3(MAXCKLEV),sz3(MAXCKLEV)
      real sx5(MAXCKLEV)

      real tsxb_sum(MAXCKLEV),weight_sum,bnu_sum(MAXCKLEV)
      real tsyb_sum(MAXCKLEV),tsya_sum(MAXCKLEV)
      real ramda_3(num_segment_3), delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)

      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      integer n_h2o_3(num_segment_3),n_ch4_3(num_segment_3),n_h2o_4,n_co2_4,n_h2o_5
      integer n_n2o_3(num_segment_3),n_co2_5
      real f_h2o_3(max_coefs,num_segment_3), f_h2o_4(max_coefs), f_h2o_5(max_coefs)
      real f_ch4_3(max_coefs,num_segment_3), f_co2_4(max_coefs), f_co2_5(max_coefs)
      real f_n2o_3(max_coefs,num_segment_3)
      real f_cfc_4(max_coefs), f_cfc_5(max_coefs)
      real k_h2o_3(max_coefs,num_segment_3), k_h2o_4(max_coefs), k_h2o_5(max_coefs)
      real k_ch4_3(max_coefs,num_segment_3), k_co2_4(max_coefs), k_co2_5(max_coefs)
      real k_n2o_3(max_coefs,num_segment_3)
      real coefk_h2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_h2o_4(max_coefs,iorder,num_pressure)
      real coefk_h2o_5(max_coefs,iorder,num_pressure)
      real coefk_ch4_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_n2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_4(max_coefs,iorder,num_pressure)
      real coefk_co2_5(max_coefs,iorder,num_pressure)
      real tau_h2o_3(max_coefs,MAXCKLEV)
      real tau_h2o_4(max_coefs,MAXCKLEV)
      real tau_h2o_5(max_coefs,MAXCKLEV)
      real tau_ch4_3(max_coefs,MAXCKLEV)
      real tau_n2o_3(max_coefs,MAXCKLEV)
      real tau_co2_4(max_coefs,MAXCKLEV)
      real tau_cfc_4(max_coefs,MAXCKLEV)
      real tau_cfc_5(max_coefs,MAXCKLEV)
      real tsxb_up(MAXCKLEV), tsxb_up_sum(MAXCKLEV)

      common/corrkcoef/ n_h2o_3, f_h2o_3, k_h2o_3, coefk_h2o_3,
     +  n_ch4_3, f_ch4_3, k_ch4_3, coefk_ch4_3,
     +  n_n2o_3, f_n2o_3, k_n2o_3, coefk_n2o_3,
     +  n_h2o_4, f_h2o_4, k_h2o_4, coefk_h2o_4,
     +  n_co2_4, f_co2_4, k_co2_4, coefk_co2_4,
     +  n_h2o_5, f_h2o_5, k_h2o_5, coefk_h2o_5,
     +  n_co2_5, f_co2_5, k_co2_5, coefk_co2_5


      real tauf_3(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_4(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_5(max_coefs,max_coefs,MAXCKLEV)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),
     +  solarrad(MAXCKLEV),ratio_up(MAXCKLEV)
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg

      nlay = klevs - 1			! nlay = the number of layers (= number of levels - 1)
      secvz = 1. / COS(.01745329*vz)
      if(sz .lt. twilit) then
        secsz = 1. / COS(.01745329*sz)
        cossz = COS(.01745329*sz)
      else
        cossz = 0.
        secsz = 0.
      end if
c
      ichan = 1
      bsun = planck_ckm(ichan,sc3temp)
      bsun = bsun * dis

	  
C   -----------------------------------------------------------------
C    First, we do channel 3 downward in the solar direction and
C    upward in the viewing direction
C   -----------------------------------------------------------------
	  
      surfref = 1. - surfemit3
      if (surfref.le.0.02) surfref=0.02
c     bref = bsun * surfref * cossz
c**	  write(28,*) bref,surfref

c     weight_sum = 0.0
      do m = 1,klevs
        tsxb_sum(m) = 0.0
        tsya_sum(m) = 0.0
        tsyb_sum(m) = 0.0
        bnu_sum(m) = 0.0
        tsxb_up_sum(m) = 0.0
      end do

      do iseg = 1,num_segment_3
 
      do m = 1,klevs
        tsx(m) = 0.0
c       tsxa(m) = 0.0
        tsxb(m) = 0.0
        tsya(m) = 0.0
        tsyb(m) = 0.0
        tsxb_up(m) = 0.0
      end do
	  
      tb = 0.
	   
c     back = bnu0(1) * surfemit3

      call ck_ckm(n_h2o_3(iseg),nlay,ulay,f_h2o_3(1,iseg),k_h2o_3(1,iseg),
     +  coefk_h2o_3(1,1,1,iseg),play,tlay,tau_h2o_3)
      call ck_ckm(n_ch4_3(iseg),nlay,uch4,f_ch4_3(1,iseg),k_ch4_3(1,iseg),
     +  coefk_ch4_3(1,1,1,iseg),play,tlay,tau_ch4_3)
      call ck_ckm(n_n2o_3(iseg),nlay,un2o,f_n2o_3(1,iseg),k_n2o_3(1,iseg),
     +  coefk_n2o_3(1,1,1,iseg),play,tlay,tau_n2o_3)

      do m = 1, nlay
        do i = 1, n_h2o_3(iseg)
          do ia = 1, n_ch4_3(iseg)
            do ib = 1, n_n2o_3(iseg)
              tauf_3(i,ia,ib, m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m) + tau_n2o_3(ib,m)
c             tauf_3(i,ia,m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m)  
            end do
          end do
        end do
      end do	   
 
      sy3(klevs) = 0.
      sz3(klevs) = bsun_seg(iseg) * cossz
      do i = 1, n_h2o_3(iseg)
        do ia = 1, n_ch4_3(iseg)
          do ib = 1, n_n2o_3(iseg)
            tsya(klevs) = tsya(klevs) + sy3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            tsyb(klevs) = tsyb(klevs) + sz3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            do m = klevs-1, 1, -1			! Solve Downward emitted (vz)
              sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +          bnu_5seg(m,iseg) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              tsya(m) = tsya(m) + sy3(m) * f_h2o_3(i,iseg) *
     +          f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)

              if(cossz .gt. 0.) then                     ! Do the downward solar component
                sz3(m) = sz3(m+1) * exp(-tauf_3(i,ia,ib,m) * secsz)
                tsyb(m) = tsyb(m) + sz3(m) * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)         

                sxz3_up = sz3(m)
                do m_up = m+1, klevs                      ! Do the reflected upward solar component (vz)
                  sxz3_up = sxz3_up * exp(-tauf_3(i,ia,ib,m_up-1) * secvz)
                end do
                tsxb_up(m) = tsxb_up(m) + sxz3_up * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
              end if
            end do

            if(cossz .le. 0.0) then   !added 12/2005
               sz3(1) = 0.0
            endif

            sxz3(1) = sz3(1) * surfref
            tsxb(1) = tsxb(1) + sxz3(1) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +        f_n2o_3(ib,iseg)
            do m = 2, klevs
              if(cossz .gt. 0.) then                       ! Do the upward solar component (vz)
                sxz3(m) = sxz3(m-1) * exp(-tauf_3(i,ia,ib,m-1) * secvz)
                tsxb(m) = tsxb(m) + sxz3(m) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +            f_n2o_3(ib,iseg) 
              end if
            end do
          end do
        end do
      end do

c     weight_sum = weight_sum + delw_3(iseg) * response_3(iseg)
      do m = 1,klevs
c       weight_sum(m) = weight_sum(m) + delw_3(iseg) * response_3(iseg)
        tsxb_sum(m) = tsxb_sum(m) + tsxb(m) * delw_3(iseg) * response_3(iseg)
        tsya_sum(m) = tsya_sum(m) + tsya(m) * delw_3(iseg) * response_3(iseg)
        tsyb_sum(m) = tsyb_sum(m) + tsyb(m) * delw_3(iseg) * response_3(iseg)
        tsxb_up_sum(m) = tsxb_up_sum(m) + tsxb_up(m) * delw_3(iseg) *
     +    response_3(iseg)
        bnu_sum(m) = bnu_sum(m) + bnu_5seg(m,iseg) * delw_3(iseg) * response_3(iseg)
      end do
      end do

	  
c  Now let's calculate the attenuation ratios for each layer for the 
c  reflected part and a composite optical depth for the emission part

      solarrad(klevs) = bsun * cossz
      do m = nlay,1,-1
        if(tsya_sum(m) - bnu_sum(m) .ne. 0.) then
          attenuate(m,1) = (tsya_sum(m) - bnu_sum(m)) / 
     +      (tsya_sum(m+1) - bnu_sum(m))
        else
          attenuate(m,1) = 0.
        end if
        if(cossz .gt. 0.) then
          ratiorefl(m) = tsxb_sum(m+1) / tsxb_sum(m)
          ratio_up(m) = tsxb_up_sum(m) / tsyb_sum(m)
          solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1)
c         solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1) +   ! Include downward emission
c    +      tsya_sum(m)/weight_sum                  
        else
          ratiorefl(m) = 1.
          solarrad(m) = 0.
        end if
      end do
      if(cossz .gt. 0.) then
        TotalRefRat  = tsxb_sum(nlay + 1) / tsxb_sum(1)
      else
        TotalRefRat  = 1.
      end if

      Radorig = 0.
      do m = nlay, 1, -1                                         ! Correctly add the downwelling
c       solarrad(m) = solarrad(m) + tsya_sum(m) / weight_sum
        iupdown = 2
        ichan = 1
        laylow = m
        layhigh = m
        call UseCorrK_down_ckm(RadOrig,laylow,layhigh,iupdown,ichan,vz,Radfinal)
        solarrad(m) = solarrad(m) + Radfinal
        RadOrig = Radfinal
      end do
 
C   -----------------------------------------------------------------
C    Next, we do channel 4 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.
      
c     back = bnu0(2) * surfemit4

      call ck_ckm(n_h2o_4,nlay,ulay,f_h2o_4,k_h2o_4,
     +  coefk_h2o_4,play,tlay,tau_h2o_4)
      call ck_ckm(n_co2_4,nlay,uco2,f_co2_4,k_co2_4,
     +  coefk_co2_4,play,tlay,tau_co2_4)
      call ck_cfc_4_ckm(n_cfc_4,nlay,uf12,f_cfc_4,
     +  play,tlay,tau_cfc_4)

      do m = 1, nlay
        call window_4_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i = 1, n_h2o_4
          do ia = 1, n_co2_4
            do ib = 1, n_cfc_4
              tauf_4(i,ia,ib,m) = tau_h2o_4(i,m) + twin + 
     +          tau_co2_4(ia,m) + tau_cfc_4(ib,m)
            end do
          end do
        end do      
      end do
 
      sx4(klevs) = 0.0
      do i = 1, n_h2o_4
        do ia = 1, n_co2_4
          do ib = 1, n_cfc_4
            do m = klevs-1,1,-1
              sx4(m) = sx4(m+1) * exp(-tauf_4(i,ia,ib,m) * 
     *          secvz) + bnu(m,2) * (1.0 - 
     *          exp(-tauf_4(i,ia,ib,m) * secvz))
              tsx(m) = tsx(m) + sx4(m) * f_h2o_4(i) * 
     *          f_co2_4(ia) * f_cfc_4(ib)
            end do
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth
	  
      do m = 1,nlay
        if(tsx(m+1) - bnu(m,2) .ne. 0.) then
	  attenuate(m,2) = (tsx(m) - bnu(m,2)) / (tsx(m+1) - 
     *      bnu(m,2))
        else
          attenuate(m,2) = 0.
        end if
c       if (attenuate(m,2).lt.0.0.or.attenuate(m,2).gt.2.0) then
c         mmmm=1
c       end if
      end do
	   
C   -----------------------------------------------------------------
C    Now, we do channel 5 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.

c     back = bnu0(3) * surfemit5

      call ck_ckm(n_h2o_5,nlay,ulay,f_h2o_5,k_h2o_5,
     +  coefk_h2o_5,play,tlay,tau_h2o_5)
      call ck_cfc_5_ckm(n_cfc_5,nlay,uf11,f_cfc_5,
     +  play,tlay,tau_cfc_5)

      do m=1, nlay
        call window_5_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i=1, n_h2o_5
          do ib=1, n_cfc_5
            tauf_5(i,ib,m) = tau_h2o_5(i,m) + twin + 
     +        tau_cfc_5(ib,m)
          end do
        end do      
      end do
      
 
      sx5(klevs) = 0.0
      do i = 1, n_h2o_5
        do ib = 1, n_cfc_5
          do m = klevs-1,1,-1
            sx5(m) = sx5(m+1) * exp(-tauf_5(i,ib,m) * 
     *        secvz) + bnu(m,3) * (1.0 - 
     *        exp(-tauf_5(i,ib,m) * secvz))
            tsx(m) = tsx(m) + sx5(m) * f_h2o_5(i) * f_cfc_5(ib)
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth

      do m = 1,nlay
        if(tsx(m+1) - bnu(m,3) .ne. 0.) then
	  attenuate(m,3) = (tsx(m) - bnu(m,3)) / (tsx(m+1) - 
     +      bnu(m,3))
        else
          attenuate(m,3) = 0.
        end if
      end do
      
 100  format(3i10,5f10.3)
	  
      return 
      end


c***************************************************************************************** 
      Subroutine CorrKsetup_mtsat_ckm(vz,sz,dis,surfemit3,surfemit4,surfemit5)
C
C**created 5/2006
c
c	This subroutine sets up the correlated Ks.
c
c	Kratz's code is used to calculate upward flux at each level.  I then solve for
c	  a total optical depth for each layer.  These are then used in routine "UseCorrK_ckm".
c     These optical depths are only valid for the given angle.
c
c   The major changes in this version:
c      1)  The surface emittances are used (these are now input parameters)
c      2)  Channel 3 surface reflectance = 1. - surfemit3 instead of a constant value of 0.10
c      3)  For channel 3, we now solve downward first.  We then use the solar radiance at
c          the surface for bref, instead of using the radiance at the TOA
c      4)  We now use subroutine ck for all correlated k calculations (except for those
c          that involve CFCs).  The inputs to subroutine ck have changed to include the 
c          coefficients.  These coefficients are now read from an ancillary data file during
c          the first call to CorrKsetup_ckm
c      5)  Kratz's latest edition of correlated ks for Channel 3 has the channel divided into 5 
c          spectral intervals.  Comparisons with the old version shows only a 1-2% improvement.
c          I've decided to stick with the old version.  The new version would have required 
c          separate bnu arrays for each spectral interval & the run-time of the Channel 3 
c          portion of the code would have increased by a factor of 5.
c      6)  Kratz has added absorption by CFCs to channels 4 & 5.  This changes the COMMON 
c          block /kprofs/ which now includes the array "ucfc".  This array is filled in 
c          subroutine fillatmos_ckm.  The correlated k's are calculated in subroutines,
c          ckf12_mtsat_ckm & ckf11_mtsat_ckm.
c      7)  The window is now modeled using the Clough continuum, instead of Roberts.  This
c          uses a parameterization developed by Fred Rose.  The old Subroutine "window" has been
c          replaced by two channel-specific routines, window_4_ckm & window_5_ckm.
c      8)  We no longer calculate the total optical depth (totaltau).  This was an unnecessary 
c          and time-consuming step that I did originally to check if the numbers looked correct.  
c          This led to some of the problems that we had with NAN results when we tried to take the
c          logarithm of 0. or a negative number.
c          We now just calculate the attenuation.  This causes a change in Subroutines UseCorrK_ckm  & GetSurtTemp
c          since we will no longer need to take the EXP of the optical depth.  It also changes the common
c          block /ktaus/.
c      9)  The layer indexing in the code was changed slightly from Kratz's original.
c          Now, level 1 = surface, level 2 = the first level above the surface, etc. up to level klevs
c          And, layer 1 = the layer between levels 1 & 2, etc. up to layer (klevs - 1) = nlay
c          Kratz's old numbering system was different in the upward and downward directions.
c     10)  I've changed the cos(sz) & cos(vz) into sec(sz) & sec(vz) 
c     11)  I now check to see if we are dividing by 0. when calculating the layer attenuation.
c          This should eliminate NAN problems when there is a temperature inversion.
c  
c
c	It is assumed that Subroutine fillatmos_ckm has been called to fill in values of bnu for
c    each atmospheric layer and bnu0 for the surface radiance.
c
c  INPUT PARAMETERS:	
c		vz          = viewing zenith
c		sz          = solar zenith angle
c       surfemit3   = VIRS Channel 3 (3.75 µm) surface emittance
c       surfemit4   = VIRS Channel 4 (10.8 µm) surface emittance
c       surfemit5   = VIRS Channel 5 (11.9 µm) surface emittance
c       bnu(ilay)   = emitted radiance from layer ilay
c       bnu0(ichan) = emitted radiance for the surface for channel ichan (should be planck(SkinTemp))
c                     where ichan = 1 for Channel 3
c                           ichan = 2 for Channel 4
c                           ichan = 3 for Channel 5
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      parameter (num_segment_3 = 5)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV),bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real dis

      real sx4(MAXCKLEV),tsx(MAXCKLEV)
c     real sx3(MAXCKLEV)
      real sxz3(MAXCKLEV)
      real tsxb(MAXCKLEV)
c     real tsxa(MAXCKLEV)
      real tsyb(MAXCKLEV),tsya(MAXCKLEV)
      real sy3(MAXCKLEV),sz3(MAXCKLEV)
      real sx5(MAXCKLEV)

      real tsxb_sum(MAXCKLEV),weight_sum,bnu_sum(MAXCKLEV)
      real tsyb_sum(MAXCKLEV),tsya_sum(MAXCKLEV)
      real ramda_3(num_segment_3), delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)

      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      integer n_h2o_3(num_segment_3),n_ch4_3(num_segment_3),n_h2o_4,n_co2_4,n_h2o_5
      integer n_n2o_3(num_segment_3),n_co2_5
      real f_h2o_3(max_coefs,num_segment_3), f_h2o_4(max_coefs), f_h2o_5(max_coefs)
      real f_ch4_3(max_coefs,num_segment_3), f_co2_4(max_coefs), f_co2_5(max_coefs)
      real f_n2o_3(max_coefs,num_segment_3)
      real f_cfc_4(max_coefs), f_cfc_5(max_coefs)
      real k_h2o_3(max_coefs,num_segment_3), k_h2o_4(max_coefs), k_h2o_5(max_coefs)
      real k_ch4_3(max_coefs,num_segment_3), k_co2_4(max_coefs), k_co2_5(max_coefs)
      real k_n2o_3(max_coefs,num_segment_3)
      real coefk_h2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_h2o_4(max_coefs,iorder,num_pressure)
      real coefk_h2o_5(max_coefs,iorder,num_pressure)
      real coefk_ch4_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_n2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_4(max_coefs,iorder,num_pressure)
      real coefk_co2_5(max_coefs,iorder,num_pressure)
      real tau_h2o_3(max_coefs,MAXCKLEV)
      real tau_h2o_4(max_coefs,MAXCKLEV)
      real tau_h2o_5(max_coefs,MAXCKLEV)
      real tau_ch4_3(max_coefs,MAXCKLEV)
      real tau_n2o_3(max_coefs,MAXCKLEV)
      real tau_co2_4(max_coefs,MAXCKLEV),tau_co2_5(max_coefs,MAXCKLEV)
      real tau_cfc_4(max_coefs,MAXCKLEV)
      real tau_cfc_5(max_coefs,MAXCKLEV)
      real tsxb_up(MAXCKLEV), tsxb_up_sum(MAXCKLEV)
      real taufact(3)

      common/corrkcoef/ n_h2o_3, f_h2o_3, k_h2o_3, coefk_h2o_3,
     +  n_ch4_3, f_ch4_3, k_ch4_3, coefk_ch4_3,
     +  n_n2o_3, f_n2o_3, k_n2o_3, coefk_n2o_3,
     +  n_h2o_4, f_h2o_4, k_h2o_4, coefk_h2o_4,
     +  n_co2_4, f_co2_4, k_co2_4, coefk_co2_4,
     +  n_h2o_5, f_h2o_5, k_h2o_5, coefk_h2o_5,
     +  n_co2_5, f_co2_5, k_co2_5, coefk_co2_5

      real tauf_3(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_4(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_5(max_coefs,max_coefs,MAXCKLEV)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),
     +  solarrad(MAXCKLEV),ratio_up(MAXCKLEV)
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg


      nlay = klevs - 1			! nlay = the number of layers (= number of levels - 1)
      secvz = 1. / COS(.01745329*vz)
      if(sz .lt. twilit) then
        secsz = 1. / COS(.01745329*sz)
        cossz = COS(.01745329*sz)
      else
        cossz = 0.
        secsz = 0.
      end if
c
      ichan = 1
      bsun = planck_ckm(ichan,sc3temp)
      bsun = bsun * dis

	  
C   -----------------------------------------------------------------
C    First, we do channel 5/3.75um downward in the solar direction and
C    upward in the viewing direction
C   -----------------------------------------------------------------
	  
      surfref = 1. - surfemit3
      if (surfref.le.0.02) surfref=0.02

      do m = 1,klevs
        tsxb_sum(m) = 0.0
        tsya_sum(m) = 0.0
        tsyb_sum(m) = 0.0
        bnu_sum(m) = 0.0
        tsxb_up_sum(m) = 0.0
      end do

      do iseg = 1,num_segment_3
 
      do m = 1,klevs
        tsx(m) = 0.0
        tsxb(m) = 0.0
        tsya(m) = 0.0
        tsyb(m) = 0.0
        tsxb_up(m) = 0.0
      end do
	  
      tb = 0.
	   
c     back = bnu0(1) * surfemit3
  
      if(iseg .eq. 1)then
        taufact(1)=0.869372
        taufact(2)=0.917775
      else
        taufact(1)=1.0
        taufact(2)=1.0
      endif
      taufact(3)=1.0
      
      call ck_mtsat_ckm(n_h2o_3(iseg),nlay,ulay,f_h2o_3(1,iseg),k_h2o_3(1,iseg),
     +  coefk_h2o_3(1,1,1,iseg),play,tlay,taufact(1),tau_h2o_3)      
      call ck_mtsat_ckm(n_ch4_3(iseg),nlay,uch4,f_ch4_3(1,iseg),k_ch4_3(1,iseg),
     +  coefk_ch4_3(1,1,1,iseg),play,tlay,taufact(2),tau_ch4_3)     
      call ck_mtsat_ckm(n_n2o_3(iseg),nlay,un2o,f_n2o_3(1,iseg),k_n2o_3(1,iseg),
     +  coefk_n2o_3(1,1,1,iseg),play,tlay,taufact(3),tau_n2o_3)

      do m = 1, nlay
        do i = 1, n_h2o_3(iseg)
          do ia = 1, n_ch4_3(iseg)
            do ib = 1, n_n2o_3(iseg)
              tauf_3(i,ia,ib, m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m) + tau_n2o_3(ib,m) 
c	      print*,iseg,tau_h2o_3(i,m),tau_ch4_3(ia,m),tau_n2o_3(ib,m)
            end do
          end do
        end do
      end do	   
 
      sy3(klevs) = 0.
      sz3(klevs) = bsun_seg(iseg) * cossz
      do i = 1, n_h2o_3(iseg)
        do ia = 1, n_ch4_3(iseg)
          do ib = 1, n_n2o_3(iseg)
            tsya(klevs) = tsya(klevs) + sy3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            tsyb(klevs) = tsyb(klevs) + sz3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            do m = klevs-1, 1, -1			! Solve Downward emitted (vz)
              sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +          bnu_5seg(m,iseg) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              tsya(m) = tsya(m) + sy3(m) * f_h2o_3(i,iseg) *
     +          f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)

              if(cossz .gt. 0.) then                     ! Do the downward solar component
                sz3(m) = sz3(m+1) * exp(-tauf_3(i,ia,ib,m) * secsz)
                tsyb(m) = tsyb(m) + sz3(m) * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)         

                sxz3_up = sz3(m)
                do m_up = m+1, klevs                      ! Do the reflected upward solar component (vz)
                  sxz3_up = sxz3_up * exp(-tauf_3(i,ia,ib,m_up-1) * secvz)
                end do
                tsxb_up(m) = tsxb_up(m) + sxz3_up * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
              end if
            end do

            if(cossz .le. 0.0) then   !added 12/2005
               sz3(1) = 0.0
            endif

            sxz3(1) = sz3(1) * surfref
            tsxb(1) = tsxb(1) + sxz3(1) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +        f_n2o_3(ib,iseg)
            do m = 2, klevs
              if(cossz .gt. 0.) then                       ! Do the upward solar component (vz)
                sxz3(m) = sxz3(m-1) * exp(-tauf_3(i,ia,ib,m-1) * secvz)
                tsxb(m) = tsxb(m) + sxz3(m) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +            f_n2o_3(ib,iseg) 
              end if
            end do
          end do
        end do
      end do

      do m = 1,klevs
        tsxb_sum(m) = tsxb_sum(m) + tsxb(m) * delw_3(iseg) * response_3(iseg)
        tsya_sum(m) = tsya_sum(m) + tsya(m) * delw_3(iseg) * response_3(iseg)
        tsyb_sum(m) = tsyb_sum(m) + tsyb(m) * delw_3(iseg) * response_3(iseg)
        tsxb_up_sum(m) = tsxb_up_sum(m) + tsxb_up(m) * delw_3(iseg) *
     +    response_3(iseg)
        bnu_sum(m) = bnu_sum(m) + bnu_5seg(m,iseg) * delw_3(iseg) * response_3(iseg)
      end do
      end do

	  
c  Now let's calculate the attenuation ratios for each layer for the 
c  reflected part and a composite optical depth for the emission part

      solarrad(klevs) = bsun * cossz
      do m = nlay,1,-1
        if(tsya_sum(m) - bnu_sum(m) .ne. 0.) then
          attenuate(m,1) = (tsya_sum(m) - bnu_sum(m)) / 
     +      (tsya_sum(m+1) - bnu_sum(m))
        else
          attenuate(m,1) = 0.
        end if
        if(cossz .gt. 0.) then
          ratiorefl(m) = tsxb_sum(m+1) / tsxb_sum(m)
          ratio_up(m) = tsxb_up_sum(m) / tsyb_sum(m)
          solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1)              
        else
          ratiorefl(m) = 1.
          solarrad(m) = 0.
        end if
      end do
      if(cossz .gt. 0.) then
        TotalRefRat  = tsxb_sum(nlay + 1) / tsxb_sum(1)
      else
        TotalRefRat  = 1.
      end if

      Radorig = 0.
      do m = nlay, 1, -1                                         ! Correctly add the downwelling
        iupdown = 2
        ichan = 1
        laylow = m
        layhigh = m
        call UseCorrK_down_ckm(RadOrig,laylow,layhigh,iupdown,ichan,vz,Radfinal)
        solarrad(m) = solarrad(m) + Radfinal
        RadOrig = Radfinal
      end do
 
C   -----------------------------------------------------------------
C    Next, we do channel 2/10.8um for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.
      
c     back = bnu0(2) * surfemit4

      taufact(1)=1.0
      taufact(2)=1.0

      call ck_mtsat_ckm(n_h2o_4,nlay,ulay,f_h2o_4,k_h2o_4,coefk_h2o_4,play,
     &            tlay,taufact(1),tau_h2o_4)
      call ck_mtsat_ckm(n_co2_4,nlay,uco2,f_co2_4,k_co2_4,coefk_co2_4,play,
     &            tlay,taufact(2),tau_co2_4)
      call ckf12_mtsat_ckm(n_cfc_4,nlay,uf12,f_cfc_4,play,tlay,tau_cfc_4)

      do m = 1, nlay
        call window_mtsat_ch2_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i = 1, n_h2o_4
          do ia = 1, n_co2_4
            do ib = 1, n_cfc_4
              tauf_4(i,ia,ib,m) = tau_h2o_4(i,m) + twin + tau_co2_4(ia,m) + 
     &                           tau_cfc_4(ib,m)
c              print*,'10.8um:',m,tau_h2o_4(i,m),tau_co2_4(ia,m),tau_cfc_4(ib,m),twin
            end do
          end do
        end do      
      end do
 
      sx4(klevs) = 0.0
      do i = 1, n_h2o_4
        do ia = 1, n_co2_4
          do ib = 1, n_cfc_4
            do m = klevs-1,1,-1
              sx4(m) = sx4(m+1) * exp(-tauf_4(i,ia,ib,m) * secvz) + bnu(m,2) *
     &                (1.0 - exp(-tauf_4(i,ia,ib,m) * secvz))
            tsx(m) = tsx(m) + sx4(m) * f_h2o_4(i) * f_co2_4(ia) * f_cfc_4(ib)
            end do
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth
	  
      do m = 1,nlay
        if(tsx(m+1) - bnu(m,2) .ne. 0.) then
	  attenuate(m,2) = (tsx(m) - bnu(m,2)) / (tsx(m+1) - bnu(m,2))
c	  print*,m,attenuate(m,2)
        else
          attenuate(m,2) = 0.
        end if
      end do
	   
C   -----------------------------------------------------------------
C    Now, we do channel 3/12um for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.

c     back = bnu0(3) * surfemit5

     
      taufact(1)=1.0
      taufact(2)=0.905499
     
      call ck_mtsat_ckm(n_h2o_5,nlay,ulay,f_h2o_5,k_h2o_5,coefk_h2o_5,play,
     &                  tlay,taufact(1),tau_h2o_5)
      call ck_mtsat_ckm(n_co2_5,nlay,uco2,f_co2_5,k_co2_5,coefk_co2_5,play,
     &                 tlay,taufact(2),tau_co2_5)
      call ckf11_mtsat_ckm(n_cfc_5,nlay,uf11,f_cfc_5,play,tlay,tau_cfc_5)

      do m=1, nlay
        call window_mtsat_ch3_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i=1, n_h2o_5
          do ia = 1, n_co2_5
            do ib = 1, n_cfc_5
              tauf_5(i,ib,m) = tau_h2o_5(i,m) + twin + tau_co2_5(ia,m) + 
     &                        tau_cfc_5(ib,m)
c              print*,'12.0um:',i,tau_h2o_5(i,m),tau_co2_5(ia,m),tau_cfc_5(ib,m),twin
            end do
          enddo
        end do      
      end do
 
      sx5(klevs) = 0.0
      do i = 1, n_h2o_5
        do ia = 1, n_co2_5
          do ib = 1, n_cfc_5
            do m = klevs-1,1,-1
              sx5(m) = sx5(m+1) * exp(-tauf_5(i,ib,m) * secvz) + bnu(m,3) *
     &                (1.0 - exp(-tauf_5(i,ib,m) * secvz))
              tsx(m) = tsx(m) + sx5(m) * f_h2o_5(i) * f_co2_5(ia) * f_cfc_5(ib)
            enddo
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth

      do m = 1,nlay
        if(tsx(m+1) - bnu(m,3) .ne. 0.) then
	  attenuate(m,3) = (tsx(m) - bnu(m,3)) / (tsx(m+1) - 
     +      bnu(m,3))
        else
          attenuate(m,3) = 0.
        end if
      end do
	  
      return 
      end
c***********************************************************************************************
	   
      Subroutine UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Radfinal)
c
c	This subroutine uses the top-to-bottom of layer radiance ratios and composite
c   optical depths calculated in CorrKsetup_ckm to correct radiances for gaseous absorption.
c   These are valid only for the viewing zenith angles at which they were calculated.
c		
c		RadOrig = Starting radiance (to be corrected)
c		Radfinal = Final, corrected radiance
c		laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c		layhigh = highest layer of atmosphere to be used
c		iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c               ichan = channel = 3 for 12.0 µm
c                           = 2 for 10.7 µm
c                           = 1 for  3.75 µm
c	
c	In general, when calling UseCorrKRatio_ckm when dealing with clouds, you should use either:
c			laylow = icloudlayer + 1
c		or	layhigh = icloudlayer
c	 where	icloudlayer is derived from subroutine FindLayerCloud
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)
      real ratio_up(MAXCKLEV)
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /krads/ bnu,bnu0
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
	  
	  
      if(layhigh .ge. laylow) then
        if(iupdown .eq. 1) then		! iupdown = 1 means solve upward (add gas absorption)
	  
          do m = laylow,layhigh
            tran = attenuate(m,ichan)
            if(tran .le. 0.05)print*,'Warning, corrk up opaque layer (m,ichan):',m,ichan
            if(m .eq. laylow) then
              tsy = RadOrig * tran + (1. - tran) * bnu(m,ichan)
	    else
              tsy = tsy * tran + (1. - tran) * bnu(m,ichan)
            end if
c            print*,m,tsy
          end do
          Radfinal = tsy
         	  
        else if(iupdown .eq. 2) then	! iupdown = 2 means solve downward (remove gas absorption)
	    
	  do m = layhigh,laylow,-1
	    tran = attenuate(m,ichan)
            if(tran .le. 0.05)print*,'Warning, corrk down opaque layer (m,ichan):',m,ichan
	    if(m .eq. layhigh) then
              if(tran .gt. 0.05)then
	        tsy = (RadOrig - (1. - tran) * bnu(m,ichan)) / tran
              else
                tsy = RadOrig
              endif
	    else
              if(tran .gt. 0.05)then
                tsy = (tsy - (1. - tran) * bnu(m,ichan)) / tran             
              endif
	    end if
c            print*,m,tsy
          end do
          
          Radfinal = tsy
          		  
        end if
	  
      else
	Radfinal = RadOrig
      end if
	  
      return
      end
C********************************************************************************************
	  
      Subroutine UseCorrKTemp_ckm(TOrig,laylow,layhigh,iupdown,ichan,Tfinal)
c
c	This subroutine uses the atmospheric optical depths calculated in CorrKsetup_ckm to correct
c	radiances for gaseous absorption.  These optical depths are valid only for the viewing
c	zenith angles at which they were calculated
c		
c		TOrig = Starting temperature (to be corrected)
c		Tfinal = Final, corrected temperature
c		laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c		layhigh = highest layer of atmosphere to be used
c		iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c		ichan = channel = 3 for 12.0µ 
c					    = 2 for 10.7µ 
c                       =1 for 3.9µ
c
      RadOrig = planck_ckm(ichan,TOrig)
	  
      Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Radfinal)
	  
      Tfinal = plancki_ckm(ichan,Radfinal)
	  
      return
      end

      Subroutine UseCorrKTempE_ckm(TOrig,laylow,layhigh,iupdown,
     +  ichan,Tfinal,SurfEmit)
c
c	This subroutine uses the atmospheric optical depths calculated in CorrKsetup_ckm to correct
c	radiances for gaseous absorption.  These optical depths are valid only for the viewing
c	zenith angles at which they were calculated
c		
c		TOrig = Starting temperature (to be corrected)
c		Tfinal = Final, corrected temperature
c		laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c		layhigh = highest layer of atmosphere to be used
c		iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c		ichan = channel = 3 for 12.0µ 
c					    = 2 for 10.7µ 
c                       =1 for 3.9µ
c
	  
c     iwr = 0

      if(iupdown .eq. 1)then
        RadOrig = planck_ckm(ichan,TOrig) * SurfEmit	  
        Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Radfinal)
	Tfinal = plancki_ckm(ichan,Radfinal)
      elseif(iupdown .eq. 2)then
        RadOrig = planck_ckm(ichan,TOrig)  
        Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Radfinal)
        Radfinal=Radfinal/SurfEmit
	Tfinal = plancki_ckm(ichan,Radfinal)
      else
        print*,'iupdown flag must be 1 or 2'
        stop
      endif
	  
      return
      end
	  
      Subroutine GetSurfTemp_ckm(T4toa,Tsurf,sza,vz,dis,xp,
     +  SurfEmit3,SurfEmit4,SurfEmit5,DataType)
c
c	This subroutine derives a surface skin temperature (Tsurf) using the observed
c	IR window TOA clear-sky temperature (T4toa)
c
c	Once Tsurf is derived, this value should be used as tk(1) and the values of
c	bnu0(i), i=1,3 need to be changed for consistency.
c
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)

      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0

      real SkinT, precip_water
      COMMON /SkinTkeep/ SkinT, precip_water

      REAL SurfEmit3, SurfEmit4, SurfEmit5

      integer DataType

      Ttop = 0.
      Call FindCloudLayer_ckm(Ttop,laytop)
      laylow = 1
      layhigh = laytop
      iupdown = 1
      ichan = 2
c       SkinT= moa skin
      Call UseCorrKTempE_ckm(SkinT,laylow,layhigh,iupdown,ichan,T4model,
     +  SurfEmit4)

      if(ABS(T4toa - T4model) .gt. 0.2) then	! Don't bother if the difference is small

        if(XP .lt. 0.75)then		! Use a small sigma value for water regions
          sigtsk = 1.0
        else
          sigtsk = 5.0			! Use a larger sigma value for water regions
        end if
        siglnpw = 0.1
        sigtoa = 0.5

        if(vz .le. 7.5) then
          ivzbin=1 
        else if(vz .le. 22.5) then
          ivzbin=2
        else if(vz .le. 37.5) then
          ivzbin=3
        else if(vz .le. 52.5) then
          ivzbin=4
        else if(vz .le. 67.5) then
          ivzbin=5
        else if(vz .le. 82.5) then
          ivzbin=6
        else
          print*,'vz not in range 0.0-82.5...stopping'
          stop
        endif

        call adjust_ckm(ivzbin,T4toa,T4model,sigtsk,siglnpw,sigtoa,dskin,
     +    pwratio)

        SkinTnew = SkinT - dskin

        bnu0(1) = planck_ckm(1,SkinT)
        bnu0(2) = planck_ckm(2,SkinT)
        bnu0(3) = planck_ckm(3,SkinT)
        bnu0(4) = planck_ckm(4,SkinT)

        do m = 1, klevs-1
          ulay(m) = ulay(m) * pwratio
        end do
        if (DataType.le.2) then
          call CorrKsetup_ckm(vz,sza,dis,surfemit3,surfemit4,surfemit5)
        else if (DataType.eq.3.or.DataType.eq.6) then
          call CorrKsetup_modis_ckm(vz,sza,dis,surfemit3,surfemit4,surfemit5)
	elseif(DataType .eq. 50)then
	  call CorrKsetup_mtsat_ckm(vz,sza,dis,surfemit3,surfemit4,surfemit5)
        else
          call CorrKsetup_goes_ckm(vz,sza,dis,surfemit3,surfemit4,surfemit5)
        end if
      else
        SkinTnew=SkinT
      end if

      Tsurf = SkinTnew

      return
      end

      Subroutine GetTcloud_ckm(Tcloudin,Tcloudout,Tskin,T4,Emit,icorrK,laytop,
     +                     cloudamount,SurfEmit,btrop,Rad4clearTOA)

      c = cloudamount

      if(icorrK .eq. 1) then
        Call FindCloudLayer_ckm(tcloudin,icloudlayer)
        laylow = icloudlayer + 1
        layhigh = laytop
        iupdown = 2
        ichan = 2
        Call UseCorrKTemp_ckm(T4,laylow,layhigh,iupdown,ichan,T4Corr)
        laylow = 1
        layhigh = icloudlayer
        iupdown = 1
        ichan = 2
        Call UseCorrKTempE_ckm(Tskin,laylow,layhigh,iupdown,ichan,
     +    TC4Bot,SurfEmit)
        if (TC4Bot.lt.T4Corr) T4Corr=TC4Bot
        b4TOA = planck_ckm(2,T4Corr)
        b4Bot = planck_ckm(2,TC4Bot)
c       write(44,*) ' iteration = ',jiter,' Corrected T4 = ',T4Corr,' Cloud Bottom T4 = ',TC4Bot
      else
        b4TOA = planck_ckm(2,T4)
        b4Bot = planck_ckm(2,Tskin)
      endif
c     Bcloud = (b4TOA - b4Bot) / Emit + b4Bot
      Bcloud = (b4TOA - c*(1.-Emit)*b4Bot - (1.-c)*
     +  Rad4clearTOA) / (c*Emit)
c
      IF(Bcloud .lt. btrop) then
        Bcloud = btrop
      end if
      Tcloudout = plancki_ckm(2,Bcloud)
c     if(Tcloudout .lt. TROPT) Tcloudout = TROPT

      return
      end

c      Subroutine GetTrueEmit_ckm(Tcloud,Tskin,tau,cvz,icorrK,ir,
c     +  SurfEmit,DataType,emit4)
c
c      integer DataType
c
c      emitlimit = 4.0
c
c      if(icorrK .eq. 1) then
c        Call FindCloudLayer_ckm(tcloud,icloudlayer)
c        laylow = 1
c        layhigh = icloudlayer
c        iupdown = 1
c        ichan = 2
c        Call UseCorrKTempE_ckm(Tskin,laylow,layhigh,iupdown,ichan,
c     +    TC4Bot,SurfEmit)
c      else
c        TC4Bot = Tskin
c      endif
c
c      ichan = 4
c      call zcompcala(ichan,ir,cvz,tau,TC4bot,Tcloud,DataType,emit4)
c
c      if(emit4 .gt. emitlimit) emit4 = emitlimit
c
c      return
c      end

      subroutine GetClearTOA_ckm(TC3,TC4,TC5,Rad3clearTOA,Rad4clearTOA,
     +  Rad5clearTOA)

      ichan = 1
      Rad3clearTOA = planck_ckm(ichan,TC3)
      ichan = 2
      Rad4clearTOA = planck_ckm(ichan,TC4)
      ichan = 3
      Rad5clearTOA = planck_ckm(ichan,TC5)

      return
      end

      Subroutine UseCorrKSolar_ckm(RadOrig,laylow,layhigh,iupdown,
     +  Radfinal)
c
c	This subroutine uses the top-to-bottom of layer radiance ratios calculated in CorrKsetup_ckm
c	to correct radiances for gaseous absorption.  These ratios are valid only for  
c	angles at which they were calculated (sza for downward, vza for upward)
c		
c		RadOrig = Starting radiance (to be corrected)
c		Radfinal = Final, corrected radiance
c		laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c		layhigh = highest layer of atmosphere to be used
c		iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c		ichan = 1 (3.9µ)
c
c  NOTES:
c	Solving "upward" (iupdown = 1) is used to derive the UPWARD TOA 
c	radiance from the reflected solar radiation.
c
c	The DOWNWARD solar component has been solved for in CorrKSetup_ckm.  We can
c	simply use the radiance value at the level desired.	 (iupdown = 2)   
c
c	Solar radiance that has been reflected from a cloud or the surface will
c	be attenuated at a different rate from the emitted portion.  The
c	attenuation ratio is independent of the magnitude of the solar
c	component - however, it must be handled separately from the emitted
c	component.  
c
c	These ratios are used in a similar manner to those in subroutine
c	UseCorrKRatio_ckm.  One difference is that they cannot be used to derive
c	the surface temperature from a TOA radiance.  In order to do this,
c	one must be able to separate the solar component from the emitted
c	component and one can't know this beforehand.
c
c	When iupdown = 2, the only necessary input argument is laylow
c	
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)
      real ratio_up(MAXCKLEV)

      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /krads/ bnu,bnu0
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
      
      if(iupdown .eq. 1) then				! iupdown = 1 means solve upward (add gas absorption)
        
        if(layhigh .ge. laylow) then
          if(laylow .eq. 1 .and. layhigh .eq. klevs - 1) then
            Radfinal = RadOrig * TotalRefRat
          else if (layhigh .eq. klevs-1) then
            Radfinal = RadOrig * ratio_up(laylow)
          else
            do m = laylow,layhigh
              if(m.eq.laylow) then
                tsy = RadOrig * ratiorefl(m)
              else
                tsy = tsy * ratiorefl(m)
              end if
            end do
            Radfinal = tsy
          end if
	else
	  Radfinal = RadOrig
	end if
	  
      else if(iupdown .eq. 2) then			! iupdown = 2 means downward (also add gas absorption)

c ******     begin Heck fix for tcloud < tk(klevs) ***** 4/15/98
         if (laylow.eq.klevs) then
            Radfinal = solarrad(laylow-1)
         else
c ******    end Heck fix for tcloud < tk(klevs)    ***** 4/15/98 
 
           Radfinal = solarrad(laylow)
         end if
      end if
	  
      return
      end

      Subroutine UseCorrKRef_ckm(RefOrig,laylow,layhigh,iupdown,RefFinal)
c
c       This subroutine is intended for use with Bob Arduini 3.7um snow models to take the
c       clear-sky surface values to TOA.
c	This subroutine uses the top-to-bottom of layer radiance ratios calculated in CorrKsetup_ckm 
c	to correct radiances for gaseous absorption.  These ratios are valid only for  
c	angles at which they were calculated (sza for downward, vza for upward)
c		
c		RefOrig = Starting reflectance (to be corrected)
c		RefFinal = Final, corrected reflectance
c		laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c		layhigh = highest layer of atmosphere to be used
c		iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c		ichan = 1 (3.9µ)
c
c  NOTES:
c	Solving "upward" (iupdown = 1) is used to derive the UPWARD TOA 
c	radiance from the reflected solar radiation.
c
c	Solar radiance that has been reflected from a cloud or the surface will
c	be attenuated at a different rate from the emitted portion.  The
c	attenuation ratio is independent of the magnitude of the solar
c	component - however, it must be handled separately from the emitted
c	component.  
c
c	These ratios are used in a similar manner to those in subroutine
c	UseCorrKRatio_ckm.  One difference is that they cannot be used to derive
c	the surface temperature from a TOA radiance.  In order to do this,
c	one must be able to separate the solar component from the emitted
c	component and one can't know this beforehand.
c	
	  
      include './corrk.h'
      integer MaxChanCK,iupdown
      parameter (MaxChanCK=4)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)
      real ratio_up(MAXCKLEV)

      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /krads/ bnu,bnu0
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up

      
      if(iupdown .eq. 1) then	! iupdown = 1 means solve upward (add gas absorption)

        RadOrig = RefOrig * solarrad(klevs)
        
        if(layhigh .ge. laylow) then
          if(laylow .eq. 1 .and. layhigh .eq. klevs - 1) then
            Radfinal = RadOrig * TotalRefRat
          elseif (layhigh .eq. klevs-1) then
            Radfinal = RadOrig * ratio_up(laylow)
          else
            do m = laylow,layhigh
              if(m.eq.laylow) then
                tsy = RadOrig * ratiorefl(m)
              else
                tsy = tsy * ratiorefl(m)
              end if
            end do
            Radfinal = tsy
          endif
	else
	  Radfinal = RadOrig
	end if

        RefFinal=Radfinal/solarrad(klevs)


c        print '(4f7.4)', RefOrig,RefFinal,TotalRefRat,solarrad(klevs)
	  
      elseif(iupdown .eq. 2) then  ! iupdown = 2 means solve downward

        RadOrig = RefOrig * solarrad(klevs)

        if(layhigh .ge. laylow) then
          if(laylow .eq. 1 .and. layhigh .eq. klevs - 1) then
            Radfinal = RadOrig * 1.0/TotalRefRat
          elseif (layhigh .eq. klevs-1) then
            Radfinal = RadOrig * 1.0/ratio_up(laylow)
          else
            do m = laylow,layhigh
              if(m.eq.laylow) then
                tsy = RadOrig * 1.0/ratiorefl(m)
              else
                tsy = tsy * 1.0/ratiorefl(m)
              end if
            end do
            Radfinal = tsy
          endif
	else
	  Radfinal = RadOrig
	end if 

        RefFinal=Radfinal/solarrad(klevs)

      end if
	  
      return
      end
	  
      Subroutine GetSurfRefRad_ckm(Tskin,T3TOA,sz,icorrK,B3SurfRad,
     +  SurfRef,SurfEmit)
c
c   This subroutine calculates the amount of radiance reflected by
c   the surface.  It assumes that all of the difference between the
c   channel 3 and channel 4 clear temperatures is due to reflection
c   (i.e. the channel 3 and 4 emittances are both = 1)
c
c        Tskin                  = the surface radiating skin temperature (usually = T4clear)
c        Rad3surf               = BB3(Tskin)
c        T3TOA                  = the channel 3 TOA clear temperature
c        Rad3TOA                = BB3(T3TOA)
c        Rad3skin               = Rad3TOA corrected to the surface
c        Rad3surfTOA            = Rad3surf corrected to the TOA
c        sz                     = solar zenith angle
c        B3SurfRad              = surface-reflected solar radiance
c        SurfRef                = Channel 3 surface reflectance
c
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)
      real ratio_up(MAXCKLEV)
       
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
       
      Rad3surf = planck_ckm(1,Tskin)*SurfEmit
      Rad3TOA = planck_ckm(1,T3TOA)
       
      if(icorrK .eq. 0) then
        B3SurfRad = Rad3TOA - Rad3surf
        SurfRef = B3SurfRad / (bsun * COS(.01745329*sz))
      else if(icorrK .eq. 1) then
        Ttop = 0.
        Call FindCloudLayer_ckm(Ttop,laytop)
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 1
c       Rad3Surf = Rad3Surf*SurfEmit
        Call UseCorrK_ckm(Rad3TOA,laylow,layhigh,iupdown,ichan,Rad3skin)
        iupdown = 1
        Call UseCorrK_ckm(Rad3surf,laylow,layhigh,iupdown,ichan,Rad3surfTOA)
       
        B3SurfRad = (Rad3TOA - Rad3surfTOA) / TotalRefRat
        SurfRef = B3SurfRad / solarrad(1)
c       SurfEmit = Rad3skin / Rad3surf
      end if
 
      if(B3SurfRad .le. 0.) then
        B3SurfRad = 0.
        SurfRef = 0.
      end if
 
      return
      end

      Subroutine GetT3TOA_ckm(Tskin,SurfEmit,chi,sz,icorrK,B3SurfRad,
     +  SurfRef,T3TOA)
c
c   This subroutine calculates the TOA channel 3 temperature assuming
c	that we know the surface emittance and the skin temperature.  The surface 
c	reflectance is computed using albedo = (1 - emittance) and the visible BDRM.
c
c	INPUT:
c        Tskin 			= the surface radiating skin temperature (usually = T4clear)
c        sz 			= solar zenith angle
c        SurfEmit		= Channel 3 surface emittance
c        chi            = surface clear anisotropic factor
c        icorrK         = Correlated K flag
c	OUTPUT:
c        B3SurfRad 		= surface-reflected solar radiance
c        SurfRef 		= Channel 3 surface reflectance
c        T3TOA 			= the channel 3 TOA clear temperature
c	LOCAL:
c        Rad3surf		= BB3(Tskin) * surface emittance
c        Rad3TOA 		= BB3(T3TOA)
c        Rad3surfTOA 	= Rad3surf corrected to the TOA
c
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)	
      real ratio_up(MAXCKLEV)
  
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit

	  
      Rad3surf = planck_ckm(1,Tskin) * SurfEmit
      SurfRef = (1. - SurfEmit) * chi
      if(SurfRef .lt. 0.) SurfRef = 0.
  
      if(icorrK .eq. 0) then
	B3SurfRad =  SurfRef * (bsun * COS(.01745329*sz))
	Rad3TOA = B3SurfRad + Rad3surf
      else if(icorrK .eq. 1) then
        Ttop = 0.
        Call FindCloudLayer_ckm(Ttop,laytop)
        laylow = 1
        layhigh = laytop
        iupdown = 1
        ichan = 1
c        Call UseCorrK_ckm(Rad3surf,laylow,layhigh,iupdown,ichan,Rad3surfTOA)
	
	B3SurfRad = SurfRef * solarrad(1)
	
c**	Rad3TOA = Rad3surfTOA + B3SurfRad * TotalRefRat	  ! Method #1
	
	Radtest = Rad3surf + B3SurfRad			  ! Method #2
        Call UseCorrK_ckm(Radtest,laylow,layhigh,iupdown,ichan,Rad3testTOA)
	Rad3TOA = Rad3testTOA
      end if
  
      T3TOA = plancki_ckm(1,Rad3TOA)

  
      return
      end
	  
      Subroutine GetEmit3Day_ckm(Tskin,T3TOA,chi,sz,icorrK,B3SurfRad,
     +  SurfRef,SurfEmit)
c
c   This subroutine calculates the channel 3 emittance in the daytime using
c	the emittance/reflectance relation assumed in GetT3TOA_ckm
c
c	INPUT:
c        Tskin 			= the surface radiating skin temperature (usually = T4clear)
c        sz 			= solar zenith angle
c        T3TOA 			= the channel 3 TOA clear temperature
c        chi            = surface clear anisotropic factor
c        icorrK         = Correlated K flag
c	OUTPUT:
c        SurfEmit		= Channel 3 surface emittance
c        B3SurfRad 		= surface-reflected solar radiance
c        SurfRef 		= Channel 3 surface reflectance
c	LOCAL:
c        Rad3surf		= BB3(Tskin) 
c        Rad3TOA 		= BB3(T3TOA)
c        Rad3surfTOA 	= Rad3surf corrected to the TOA
c
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)	  
      real ratio_up(MAXCKLEV)

      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
	  
      Rad3surf = planck_ckm(1,Tskin)
      Rad3TOA = planck_ckm(1,T3TOA)
  
      if(icorrK .eq. 0) then
        chisol = chi * bsun * COS(.01745329*sz)
	SurfEmit = (Rad3TOA - chisol) / (Rad3surf - chisol)
	SurfRef = (1. - SurfEmit) * chi
	B3SurfRad = (1. - SurfEmit) * chisol
      else if(icorrK .eq. 1) then
        Ttop = 0.
        Call FindCloudLayer_ckm(Ttop,laytop)
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 1
        Call UseCorrK_ckm(Rad3TOA,laylow,layhigh,iupdown,ichan,Rad3skin)	
	
	chisol = chi * solarrad(1)
	SurfEmit = (Rad3skin - chisol) / (Rad3surf - chisol)
	SurfRef = (1. - SurfEmit) * chi
	B3SurfRad = (1. - SurfEmit) * chisol
      end if
  
      if(B3SurfRad .le. 0.) then
        B3SurfRad = 0.
	SurfRef = 0.
c     else
c       SurfEmit = 1.
      end if
  
      return
      end

      Subroutine UseCorrK_down_ckm(RadOrig,laylow,layhigh,iupdown,ichan,vz,Radfinal)
c
c       This subroutine uses the atmospheric attenuation factors calculated in CorrKsetup_ckm
c   to calculate downwelling radiation.
c
c   REMEMBER: the original corr-k routines all deal with UPWELLING radiation.  The
c             array iupdown is not referring to upwelling vs. downwelling radiation, but
c             refers to whether we are adding or removing gaseous absorption as we move
c             upwards or down through the atmosphere.
c
c               RadOrig = Starting radiance (to be corrected)
c               Radfinal = Final, corrected radiance
c               laylow = lowest layer of atmosphere to be used (layer 1 is at surface)
c               layhigh = highest layer of atmosphere to be used
c               iupdown = direction of calculations = 1 means solve upward (add gas absorption)
c                                           = 2 means solve downward (remove gas absorption)
c               ichan = channel = 1 for 12.0 µm
c                                           = 2 for 10.7 µm
c                       = 3 for  3.75 µm
c       vz = the viewing zenith angle at which the atmospheric absortion was calculated
c            This is modified to adjust to an angle of 53°, which is used to calculate
c            the diffuse radiation
c
c       In general, when calling UseCorrKRatio_ckm when dealing with clouds, you should use either:
c                       laylow = icloudlayer + 1
c               or      layhigh = icloudlayer
c        where  icloudlayer is derived from subroutine FindLayerCloud
c

      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV),zk(MAXCKLEV),rhk(MAXCKLEV)
      real o3k(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV)
      real ratio_up(MAXCKLEV)
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up

      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /krads/ bnu,bnu0

      angle_factor = cos(.01745329*vz) / cos(.01745329*53.)

      if(layhigh .ge. laylow) then
        if(iupdown .eq. 2) then             ! iupdown = 2 means solve downward (add gas absorption)

          do m = layhigh,laylow,-1
            tran = attenuate(m,ichan)**angle_factor
            if(tran .le. 0.05)print*,'Warning, corrk dnwell down opaque layer (m,ichan):',m,ichan
            if(m.eq.layhigh) then
              tsy = RadOrig * tran + (1. - tran) * bnu(m,ichan)
            else
              tsy = tsy * tran + (1. - tran) * bnu(m,ichan)
            end if
          end do
          Radfinal = tsy
        else if(iupdown .eq. 1) then         ! iupdown = 1 means solve upward (remove gas absorption)

          do m = laylow,layhigh
            tran = attenuate(m,ichan)**angle_factor
            if(tran .le. 0.05)print*,'Warning, corrk dnwell up opaque layer (m,ichan):',m,ichan
            if(m.eq.laylow) then
              if(tran .gt. 0.05)then
                tsy = (RadOrig - (1. - tran) * bnu(m,ichan)) / tran
              else
                tsy=RadOrig
              endif
            else
              if(tran .gt. 0.05)then
                tsy = (tsy - (1. - tran) * bnu(m,ichan)) / tran
              endif
            end if
          end do
          Radfinal = tsy
        end if
      else
        Radfinal = RadOrig
      end if
c     write(26,*) laylow,layhigh,iupdown,ichan
c     write(26,*) RadOrig,Radfinal
      return
      end
  
c****************************************************************************************
      subroutine ck_ckm(ni,mlv,u,f,k,coefk,p,t,tau)
c
c    This subroutine will calculate the optical depths (tau) for 
c    any of Dave Kratz's standard correlated k fits.  It replaces the
c    old routines like ckco2 & ckch4, etc.  Instead of the endless
c    data statements & repeated code, we now use one subroutine that
c    gets the correlated k coefficients as an input argument.
c
c     The basic equation used to derive the k's for a given
c     temperature is:
c       k(T) = k(250) + k' * (T - 250) + k'' * (T - 250)^2
c
c  INPUT:
c     k(i)     = the k-coefficient for the ith subinterval (cm^2/gm).
c     f(i)     = the weight (the percentage of the wavenumber interval) of
c                the ith subinterval.  Note that the summation of f(i) over all
c                i's must equal 1.
c     ni       = the number of subintervals.
c     mlv      = the number of layers.
c     p        = the layer pressure (atm), 
c     t        = the layer temperature (K), 
c     u        = the gas column abundance (gm/cm^2).
c     coefk(i,jt,jp) = the temperature correction coefficients: k, k' 
c                      and k'' specified at 19 pressures which are used 
c                      to calculate the k-coefficients.
c
c  OTHER KEY VARIABLES:
c     pmb      = the layer pressure in hPa (mb).
c     stp(jp)  = the 19 standard pressures (hPa)
c     fkg(i,m) = the interpolated pressure and temperature correction
c                coefficient from coefk(i,jt,jp).
c
c  OUTPUT:
c     tau(i,m) = the optical depth for subinterval i and layer m.
c
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder, num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)

      common /standpres/stp

      do i = 1, ni
        do m = 1, mlv
          ml = 1
          tdif = t(m) - 250.0
          pmb = p(m)*1013.25
          if(pmb .lt. stp(1))then
            x1 = coefk(i,1,1) + (coefk(i,2,1) + coefk(i,3,1) * tdif)
     +        * tdif
            fkg(i,m) = x1 * pmb / stp(1)
          else if (pmb .ge. stp(19)) then
            x1 = coefk(i,1,18) + (coefk(i,2,18) + coefk(i,3,18) *
     +        tdif) * tdif
            x2 = coefk(i,1,19) + (coefk(i,2,19) + coefk(i,3,19) * 
     +        tdif) * tdif
            fkg(i,m) = x1 + (x2-x1) / (stp(19)-stp(18)) * (pmb - 
     +        stp(18))
          else
            do while(pmb .ge. stp(ml))
              ml=ml+1
            end do
            x1 = coefk(i,1,ml-1) + (coefk(i,2,ml-1) + coefk(i,3,ml-1) * 
     +        tdif) * tdif
            x2 = coefk(i,1,ml) + (coefk(i,2,ml) + coefk(i,3,ml) * 
     +        tdif) * tdif
            fkg(i,m) = x1 + (x2-x1) / (stp(ml) - stp(ml-1)) * (pmb - 
     +        stp(ml-1))
          end if
        end do
      end do
      do i = 1, ni
        do m = 1, mlv
          tau(i,m) = k(i) * u(m) * fkg(i,m)
c         write(6,*)m,i,u(m),k(i),fkg(i,m),tau(i,m)
        end do
      end do
      return
      end
c***************************************************************************************
      subroutine ck_mtsat_ckm(ni,mlv,u,f,k,coefk,p,t,fact,tau)
c
c    This subroutine will calculate the optical depths (tau) for 
c    any of Dave Kratz's standard correlated k fits.  It replaces the
c    old routines like ckco2 & ckch4, etc.  Instead of the endless
c    data statements & repeated code, we now use one subroutine that
c    gets the correlated k coefficients as an input argument.
c
c     The basic equation used to derive the k's for a given
c     temperature is:
c       k(T) = k(250) + k' * (T - 250) + k'' * (T - 250)^2
c
c  INPUT:
c     k(i)     = the k-coefficient for the ith subinterval (cm^2/gm).
c     f(i)     = the weight (the percentage of the wavenumber interval) of
c                the ith subinterval.  Note that the summation of f(i) over all
c                i's must equal 1.
c     ni       = the number of subintervals.
c     mlv      = the number of layers.
c     p        = the layer pressure (atm), 
c     t        = the layer temperature (K), 
c     u        = the gas column abundance (gm/cm^2).
c     fact     = factor to multily tau by
c     coefk(i,jt,jp) = the temperature correction coefficients: k, k' 
c                      and k'' specified at 19 pressures which are used 
c                      to calculate the k-coefficients.
c
c  OTHER KEY VARIABLES:
c     pmb      = the layer pressure in hPa (mb).
c     stp(jp)  = the 19 standard pressures (hPa)
c     fkg(i,m) = the interpolated pressure and temperature correction
c                coefficient from coefk(i,jt,jp).
c
c  OUTPUT:
c     tau(i,m) = the optical depth for subinterval i and layer m.
c
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder, num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)

      common /standpres/stp

      do i = 1, ni
        do m = 1, mlv
          ml = 1
          tdif = t(m) - 250.0
          pmb = p(m)*1013.25
          if(pmb .lt. stp(1))then
            x1 = coefk(i,1,1) + (coefk(i,2,1) + coefk(i,3,1) * tdif)
     +        * tdif
            fkg(i,m)=x1
          else if (pmb .ge. stp(19)) then
            x1 = coefk(i,1,18) + (coefk(i,2,18) + coefk(i,3,18) *
     +        tdif) * tdif
            x2 = coefk(i,1,19) + (coefk(i,2,19) + coefk(i,3,19) * 
     +        tdif) * tdif
            fkg(i,m) = x1 + (x2-x1) / (stp(19)-stp(18)) * (pmb - 
     +        stp(18))
          else
            do while(pmb .ge. stp(ml))
              ml=ml+1
            end do
            x1 = coefk(i,1,ml-1) + (coefk(i,2,ml-1) + coefk(i,3,ml-1) * 
     +        tdif) * tdif
            x2 = coefk(i,1,ml) + (coefk(i,2,ml) + coefk(i,3,ml) * 
     +        tdif) * tdif
            fkg(i,m) = x1 + (x2-x1) / (stp(ml) - stp(ml-1)) * (pmb - 
     +        stp(ml-1))
          end if
        end do
      end do
      do i = 1, ni
        do m = 1, mlv
          tau(i,m) = fact * k(i) * u(m) * fkg(i,m)
c         write(6,*)m,i,u(m),k(i),fkg(i,m),tau(i,m)
        end do
      end do
      return
      end
c**************************************************************************************
      subroutine ck_cfc_4_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=2
c     f(1)=ratio of spectral interval used vs. cfc spectral interval
c      cor=1.000000
      cor=0.998434
      f(1)=0.77778*cor
      f(2)=1.0-f(1)
      do m=1,mlv
        k(1)=3814.461*(1.0000-2.3034e-04*(t(m)-250.0)+
     *       5.1926e-07*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
        tau(2,m)=0.0
c       write(6,*)m,u(m),k(1),tau(1,m)
      end do
      return
      end
c *********************************************************************
      subroutine ckf12_mtsat_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      ni=2
c     f(1)=ratio of spectral interval used vs. f12 spectral interval
c      cor=1.000000
      cor=0.987809
      f(1)=0.87500*cor
      f(2)=1.0-f(1)
      do m=1,mlv
        k(1)=3814.461*(1.0000-2.3034e-04*(t(m)-250.0)+
     *       5.1926e-07*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
        tau(2,m)=0.0
c       write(6,*)m,u(m),k(1),tau(1,m)
      end do
      return
      end
c****************************************************************************************
      subroutine window_4_ckm (amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 880 to 970 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)

      aa(1)=9.292
      aa(2)=0.9999
      aa(3)=-0.02433
      aa(4)=0.01961
      aa(5)=-0.7897
      aa(6)=3.467e-05
      aa(7)=0.9923
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

      twin = exp ( tau_log)
c
      return
      end
c****************************************************************************************
      subroutine window_mtsat_ch2_ckm (amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.4 continuum over 880 to 970 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.4 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
c     data aa/9.292,0.9999,-0.02433,0.01961,-0.7897,3.467e-05,0.9923/
      data aa/9.322,1.0000,-0.02439,0.01067,-0.6760,-2.685e-04,0.9901/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

        twin = exp ( tau_log)
c
        return
        end
c*********************************************************************************************
      subroutine ck_cfc_5_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=2
c     f(1)=ratio of spectral interval used vs. cfc spectral interval
c     correction for spectral interval offset
c      cor=1.000000
      cor=0.945681
      f(1)=0.625*cor
      f(2)=1.0-f(1)
      do m=1,mlv
        k(1)=5611.777*(1.0000+1.2661e-03*(t(m)-250.0)+
     *       3.5571e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
        tau(2,m)=0.0
c       write(6,*)m,u(m),k(1),tau(1,m)
      end do
      return
      end
c*******************************************************************************************
      subroutine ckf11_mtsat_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      ni=1
      cor11=1.02304
      f(1)=1.0
      do m=1,mlv
        fk11=cor11*5151.471*(1.0000-2.3775e-04*(t(m)-250.0)-
     *       4.1043e-06*(t(m)-250.0)**2)
        tau(1,m)=fk11*u(m)
      end do
      return
      end
c******************************************************************************************
      subroutine window_5_ckm(amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 800 to 880 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)

      aa(1)=9.646
      aa(2)=1.000
      aa(3)=-0.02449
      aa(4)=0.03224
      aa(5)=-0.6462
      aa(6)=-2.144e-04
      aa(7)=0.9806
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

      twin = exp ( tau_log)
c
      return
      end

c***************************************************************************
      subroutine window_mtsat_ch3_ckm(amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.4 continuum over 800 to 860 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.4 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
c
c     data aa/9.199,1.006,-0.02386,0.10630, 0.8928,-2.058e-03,0.9330/
      data aa/9.499,1.003,-0.02417,0.03857, 0.7636,-1.042e-03,0.9565/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

        twin = exp ( tau_log )
c
        return
        end
c****************************************************************************
      subroutine adjust_ckm(ivzbin,TOBS,TMODEL,sigtsk,siglnpw,sigtoa,
     +  dskin,pwratio)

c     This routine will adjust or correct a sounding so that the
c     satellite observed channel 4 brightness temperature agrees
c     with the model temperature.

      include './corrk.h'
      parameter(mflx=1,nf=1, mcldc=1,ncc=1, nvar=2,ntv=2,iprt=0)
      real p0(MAXCKLEV),t0(MAXCKLEV),u0(MAXCKLEV),z0(MAXCKLEV)
      real rh0(MAXCKLEV),o30(MAXCKLEV)
      real dfx(mflx),tol(mflx),flx(0:mcldc,mflx)
      real dfdv(0:mcldc,nvar,mflx)
      real org(0:mcldc),sigs(0:mcldc),adj(0:mcldc)
      real sigv(0:mcldc,nvar),adjv(0:mcldc,nvar),badj(nvar)
      real lnpwder

      COMMON /ktemps/ mlv,p0,t0,z0,u0,rh0,o30,i85flag
      COMMON /SkinTkeep/ tskin, pw

      call getderiv_ckm(tskin,pw,ivzbin,tder,lnpwder)

c This tune_xxx_ckm algorithm is a general algorithm that includes Options
c To adjust CLOUD FRACTION , For the CLEAR SKY ONLY SET UP arrays
c must be dimensioned to have a clear and cloud component hence 
c (0:mcldc) where mcldc =1 and ncc=1
c In this set up there are 2 tuning variables Skint , lnPw hence nvar=2 and ntv=2
c Only LW Flux is being tuned Hence mflx=1 and nf=1

cc Apriori user assigned uncertanties "sigmas"
      sigs(0)=0.000001  !! Cloud FRACTION SIGMA CLEAR SKY
      sigs(1)=0.000001  !! Cloud FRACTION SIGMA CLOUDY SKY

      sigv(0,1) = sigtsk !!SKINT SIG 
      sigv(0,2) = siglnpw !!LNPW SIG  !0.05 c
      sigv(1,1) = 0.0     !!SKINT SIG 
      sigv(1,2) = 0.0     !!LNPW SIG  !0.05 c
      tol(1) 	  = sigtoa !!TOAT Tolerance (K)

cc PRECOMPUTED DERIVATIVES From Look-up Table 
      dfdv(0,1,1) = tder !! <<-- SKIN TEMP DERIVATIVE { ( TOAT(t2)-TOAT(t1) )/ (t2-t1) }
      dfdv(0,2,1) = lnpwder  !!  <<-- lnPW DERIVATIVE { (TOAT(pw2)-TOAT(pw1) )/ (lnPW2-lnPw1) }
      dfdv(1,1,1) = 0.0 !! <<-- SKIN TEMP DERIVATIVE { ( TOAT(t2)-TOAT(t1) )/ (t2-t1) }
      dfdv(1,2,1) = 0.0  !!  <<-- lnPW DERIVATIVE { (TOAT(pw2)-TOAT(pw1) )/ (lnPW2-lnPw1) }

      org(0)=1.0 	!!! ALL CLEAR SKY 
      org(1)=0.0 	!!! ALL CLEAR SKY 
      dfx(1)= TMODEL-TOBS 	!!! OBS-MODEL ERROR

      flx(0,1)= TMODEL 	!!! MODEL CSOLR
      flx(1,1)= TMODEL 	!!! FAKE IT OUT NO CLD FRAC DERIVATIVE

      call tune_xxx_ckm(mflx,nf,mcldc,ncc,org,dfx,tol,flx,sigs,adj,
     1              nvar,ntv,dfdv,sigv,adjv,badj,iprt,iok)
      dskin=badj(1)
      dlnpw=badj(2)

      pwc=exp(alog(pw)-dlnpw)  !correct bad profile based on derivative
      pwratio=pwc/pw  !find ratio to correct bad profile

      return
      end

      subroutine tune_xxx_ckm(mflx,nf,mcldc,ncc,org,dfx,tol,flx,sigs,adj,
     &  nvar,ntv,dfdv,sigv,adjv,badj,iprt,iok)

      parameter(mv=20,ret=0.0001)
      real dfx(mflx),tol(mflx)
      real flx(0:mcldc,mflx)
      real dfdv(0:mcldc,nvar,mflx)
      real org(0:mcldc),sigs(0:mcldc),adj(0:mcldc)
      real su(0:mv),u(0:mv),wt(0:mv),ce(0:mv)
      real at(0:mv,mv)
      real BBB(mv),te(mv),uv(mv),bsig(mv),r(mv)
      real atv(mv,mv),aaa(mv,mv)
      real adjv(0:mcldc,nvar)
      real a(mv+1)
      real badj(nvar)
      real sigv(0:mcldc,nvar)
      integer INDX(mv)
      integer iprt,iok

! INPUT!
! MCLDC - Max dimension for # of Cloud Conditions
! NCC  - Actual # of Cloud Conditions ond used in tuning
! ORG(0:mcldc) - unadjusted cloud fraction
! DF(mflx)   - Inital Flux error ( OBS - Modeled)
! TOL(mflx)  - 1_sigma tolerance for observed flux
! FLX(0:mcldc,mflx)   - Array holding Model Fluxes for each Cloud Condition
! sigs(0:mcldc) -1 sigma uncertanty for Cloud Cond. Fractions (ABSOLUTE)
		  
! NVAR - Max dimension for # of Tuning Variables
! NTV - actual # of tuning variables
! DFDV(0:mcldc,nvar,mflx) -Partial Derivatives for each tuning variable for each
!			Cloud Condition for each flux
!  sigv(0:mcldc,nvar)  - 1 sigma uncertanty for tuning variable
! iprt - turn printout of statistics 1=on 0=off

! OUTPUT!
! adj(0:mcldc) - sugested adjustments to be applied to Cloud Fractions
! adjv(0:mcldc,nvar) - sugested adjustments to be applied to tuning variables
! badj(nvar) - adjustment over entire footprint.

	if(mcldc.gt.mv)    stop '1 redimension tune_xxx_ckm'
	if(nvar .gt.mv)    stop '2 redimension tune_xxx_ckm'
	if(ntv  .gt.mv)    stop '3 redimension tune_xxx_ckm'
	if(ntv .gt. nvar)  stop '4 redimension tune_xxx_ckm'
	if(ncc .gt. mcldc) stop '5 redimension tune_xxx_ckm'

	ipass=0
	iok=0
1000	continue
! --- Initalize variables
	do i=1,mv
          a(i)  =0
          bbb(i)=0
	indx(i)	=0
	te(i) 	=0
c       sc(i)	=0
c       sa(i)	=0

	 do j=1,mv
	 aaa(i,j)=0.0
	 enddo
	enddo

	do  i=0,mcldc
	adj(i)	=0.0
	ce(i)	=0.0

	do j=1,nvar
	adjv(i,j)=0.0
	enddo
	enddo

	do  j=1,nvar
	 do k=1,nf
	 atv(j,k)=0.0
	 enddo
	bsig(j)  = 0.0	
	badj(j)	 = 0.0
	enddo

! --- Sum Flux @ Cloud Condition
	do ic =0,ncc
	 do k=1,nf
	 energy = org(ic)*flx(ic,k)
	 te(k) =te(k) + energy
	 ce(ic)=ce(ic)+ energy
	 enddo !k
	enddo !ic

	teall=0
	do k=1,nf
	teall=teall+te(k) 
	enddo
! --- Set up for solution

	do k=1,nf
	a(k)=-2.0*dfx(k)
	enddo
	a(nf+1)=0

	do i=0,ncc  !!! Cloud Conditions

	    if(ipass.eq.0) then
	        su(i) =  abs( sigs(i) )              !!! ALL ABSOLUTE
		if(i.ge.1 .and. org(i).le.sigs(i) ) su(i)= org(i)
	    elseif(ipass.eq.1 )then
 		su(0) =  abs( sigs(0) )              !!! CS ABSOLUTE       
	        if(i.ge.1)then
		su(i) =  abs( sigs(i) * org(i) )     !!! FS RELATIVE
		endif
	    elseif(ipass.eq.2) then
	        su(i) =  abs( sigs(i) * org(i) )     !!! ALL RELATIVE
	    elseif(ipass.ge.3 .and. ipass.le.10 ) then
		su(i) =  abs( sigs(i)*org(i))* 4.5/real(ipass**2) !REDUCE
	    else 
	        stop ' subroutine tune_xxx_ckm error ipass'
	    endif

	do k=1,nf
	at(i,k)= flx(i,k) * su(i)
	enddo

	wt(i) = ce(i) / teall  ! ENERGY FRAC for each cloud condition...
!	wt(i) =org(i)
	do j=1,ntv   !!! Bulk Sigma
	bsig(j)  =bsig(j) + sigv(i,j)  *wt(i) 
	enddo ! j

	enddo !i
! --- 
	do i=0,ncc   !!! Cloud Condition
	do j=1,ntv   !!! Property variable
	do k=1,nf    !! FLUX
	atv(j,k) = atv(j,k)+  dfdv(i,j,k) * bsig(j) * org(i) 
	enddo !k
	enddo !j
	enddo !i

! --- SOLVE !!!

	do ij=1,nf
	aaa(ij,ij)=  tol(ij)*tol(ij)
	enddo

	do ic=0,ncc !!! CLOUD CONDITIONS

	do j=1,nf
	do i=1,nf
	aaa(i,j)=aaa(i,j)+ at(ic,i)*at(ic,j)
	enddo !i
	enddo !j

	do ij= 1 , nf
	aaa(nf+1,ij)=aaa(nf+1,ij)+ at(ic,ij)*su(ic)
	aaa(ij,nf+1)=aaa(nf+1,ij)
	enddo !ij

	aaa(nf+1,nf+1)= aaa(nf+1,nf+1)+su(ic)*su(ic)
	enddo !! ic ; CLOUD CONDITION
! ---
	do iv=1,ntv  !! CLOUD PROPERTY TERMS
	do j=1,nf
	do i=1,nf
	aaa(i,j) = aaa(i,j) + atv(iv,i)*atv(iv,j)
	enddo !i
	enddo !j

	enddo  !! CLOUD PROPERTY TERMS

! --- Solve Matrix

	do k=1,nf+1
	bbb(k)=a(k)
	enddo

	call LUDCMP_ckm(AAA,NF+1,MV,INDX,DDD)
	call LUBKSB_ckm(AAA,NF+1,MV,INDX,bbb)

! --- Solve for adjustments
	do j=1,ntv !!! ADJUST PROPERTY Variables
	    uuu=0.0
	    do k=1,nf
	    uuu=uuu+ atv(j,k)*bbb(k)
	    enddo
	uv(j)= -0.5* uuu 
	badj(j)= uv(j)*bsig(j) !! Bulk Footprint Adjustment

! -- Individulal cloud condition property adjustment
	 do i=0,ncc
	     if(org(i) .gt.0) then
	       adjv( i,j ) = uv(j) *  sigv(i,j)* wt(i)/org(i)
	     else
	       adjv( i,j ) =0.0
	     endif
	enddo !i Cloud Condition
	enddo !j Property variable


	sumadj=0
	do i=0,ncc !!! ADJUST CLOUD FRACTIONS

	    uuu=0
	    do k=1,nf
	    uuu = uuu + at(i,k)*bbb(k)
	    enddo
	uuu=uuu+ su(i) * bbb(nf+1)

	  u( i ) = -0.5 * uuu
	adj( i ) = u( i )*su( i )
	sumadj = sumadj+adj(i)
	enddo

! --- Check for negative cloud fraction
	do ic =0,ncc
	 if(org(ic)+adj(ic) .lt.0  .or. org(ic)+adj(ic) .gt. 1.0+ret ) then

	 ipass=ipass+1

	 go to 1000
	 endif
	enddo
! --- Find Flux Allowances
	do k=1,nf
	R(k) =0.5 *tol(k)*bbb(k)
c       ef(k)=R(k)*tol(k)
	if(abs( R(k) ).gt.3.0 ) iok=iok+100  !! > 3 sigma FLUX
	enddo

	do ic=0,ncc
	if(abs( u(ic) ).gt.3.0 ) iok=iok+ 10  !! > 3 sigma CLOUD FRAC
	enddo
	
	do j=1,ntv
	if(abs( uv(j) ).gt.3.0 ) iok=iok+  1 !! > 3 sigma Property Variable
	enddo
! --- ERROR CHECK
	s1=0
	s2=0
	s3=0
	do ic =0,ncc
	s1=s1+org(ic)
	s2=s2+adj(ic)
	s3=s3+org(ic)+adj(ic)
	enddo

	if( abs(s1-1.00) .gt. ret .or. 
     &    abs(s2).gt.ret .or. abs(s3-1.00) .gt. ret) then
	print*,sumadj,' =Sum of Cloud fraction adjustments'
	print *, s1,s2,s3
	stop ' ERROR CLD FRAC ADJUST'
	endif
! --- ERROR CHECK

	return
	end
!========================================================================
      SUBROUTINE LUBKSB_ckm(A,N,NP,INDX,B)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END
!=============================================================================
      SUBROUTINE LUDCMP_ckm(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) STOP 'Singular matrix.'
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END

      Subroutine getderiv_ckm(SkinTin,PWin,ivzbin,tder,lnpwder)

c     This routine finds the derivative corresponding to
c     the nearest 5deg temp, .333pw box for the given skin temp
c     and precip water (no interpolation for in between 5deg/0.33
c     increments).
      
      integer iskin(30),ipw(20)
      real dtoadskin(30,20,6),dtoadlnpw(30,20,6),lnpwder   

      common/derivs/dtoadskin,dtoadlnpw
      common /adjuststuff/iskin,ipw

      smalldif=100.0
      do i=1,30
        skintdif=abs(SkinTin-iskin(i))
        if(skintdif .lt. smalldif)then
          smalldif=skintdif
          indexsk=i
        endif
      enddo

      smalldif=7.0*1000
      ipwin=int(PWin*1000)
      do i=1,20
        pwdif=abs(ipwin-ipw(i))
        if(pwdif .lt. smalldif)then
          smalldif=pwdif
          indexpw=i
        endif       
      enddo

      indexpw=21-indexpw
      tder=dtoadskin(indexsk,indexpw,ivzbin)
      lnpwder=dtoadlnpw(indexsk,indexpw,ivzbin)

      return
      end
c
c-----------------------------------------------------------------------
c=======================================================================
c-----------------------------------------------------------------------
c
      real function planck_ckm(ic,temp)
c
c-----------------------------------------------------------------------
c
c variable     description
c --------     ---------------------------------------------------------
c ic           channel number
c temp         temperature, K
c
c aplanck      blackbody radiance, (W/m/m)/micron/sr
c
c-----------------------------------------------------------------------
c
c Allocate memory
c
c   Parameters
c
      parameter (nc=4)
c
c   Arrays
c
      real c(nc,2)
c
c   Scalars
c
      integer ic,DataType
      real temp
c
c   Common blocks
c
      common /radiat/ c
      common/DT/DataType
c
c-----------------------------------------------------------------------
c
c Calculate Planck function
c
      temp_adj=temp
      if (ic.eq.1) then
        if (DataType.eq.2) then
          temp_adj = (temp_adj + 1.132)/1.004102
        else if (DataType.eq.3) then
          temp_adj = (temp_adj + 0.361)/1.001392
        end if
      end if

      planck_ckm=c(ic,1)/(exp(c(ic,2)/temp_adj)-1.)
      return
c
c-----------------------------------------------------------------------
c
c Finem respice
c
      end
c
c-----------------------------------------------------------------------
c=======================================================================
c-----------------------------------------------------------------------
c
      real function plancki_ckm(ic,rad)
c
c-----------------------------------------------------------------------
c
c variable     description
c --------     ---------------------------------------------------------
c ic           channel number
c rad          blackbody radiance, (W/m/m)/micron/sr
c
c aplancki     temperature, K
c
c-----------------------------------------------------------------------
c
c Allocate memory
c
c   Parameters
c
      parameter (nc=4)
c
c   Arrays
c
      real c(nc,2)
c
c   Scalars
c
      integer ic,DataType
      real rad,test
c
c   Common blocks
c
      common /radiat/ c
      common/DT/DataType
c
c-----------------------------------------------------------------------
c
c Calculate inverse Planck function
c
      test=1.+c(ic,1)/rad
      if (test.le.0.) then
         plancki_ckm=0.
      else
         plancki_ckm=c(ic,2)/alog(test)
      endif

      temp_adj=plancki_ckm
      if (ic.eq.1) then
        if (DataType.eq.2) then
          temp_adj = temp_adj * 1.004102 - 1.132
        else if (DataType.eq.3) then
          temp_adj = temp_adj * 1.001392 - 0.361
        end if
      end if

      plancki_ckm=temp_adj

      return
c
c-----------------------------------------------------------------------
c
c Finem respice
c
      end
      
c
c-----------------------------------------------------------------------
c=======================================================================
c-----------------------------------------------------------------------
c
      real function planck2_ckm(ic,temp)
c
c-----------------------------------------------------------------------
c
c variable     description
c --------     ---------------------------------------------------------
c ic           channel number
c temp         temperature, K
c
c aplanck      blackbody radiance, (W/m/m)/micron/sr
c
c-----------------------------------------------------------------------
c
c Allocate memory
c
c   Parameters
c
      parameter (nc=5)
c
c   Arrays
c
      real c(nc,2)
c
c   Scalars
c
      integer ic,DataType
      real temp
c
c   Common blocks
c
      common /radiat2/ c
      common/DT/DataType
c
c-----------------------------------------------------------------------
c
c Calculate Planck function
c
      temp_adj=temp
      if (ic.eq.1) then
        if (DataType.eq.2) then
          temp_adj = (temp_adj + 1.132)/1.004102
        else if (DataType.eq.3) then
          temp_adj = (temp_adj + 0.361)/1.001392
        end if
      end if

      planck2_ckm=c(ic,1)/(exp(c(ic,2)/temp_adj)-1.)
      return

c Finem respice
c
      end

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

      Subroutine CorrKsetup_goes_ckm(vz,sz,dis,surfemit3,surfemit4,surfemit5)
c
c	This subroutine sets up the correlated Ks.
c
c	Kratz's code is used to calculate upward flux at each level.  I then solve for
c	  a total optical depth for each layer.  These are then used in routine "UseCorrK_ckm".
c     These optical depths are only valid for the given angle.
c
c   The major changes in this version:
c      1)  The surface emittances are used (these are now input parameters)
c      2)  Channel 3 surface reflectance = 1. - surfemit3 instead of a constant value of 0.10
c      3)  For channel 3, we now solve downward first.  We then use the solar radiance at
c          the surface for bref, instead of using the radiance at the TOA
c      4)  We now use subroutine ck for all correlated k calculations (except for those
c          that involve CFCs).  The inputs to subroutine ck have changed to include the 
c          coefficients.  These coefficients are now read from an ancillary data file during
c          the first call to CorrKsetup_ckm
c      5)  Kratz's latest edition of correlated ks for Channel 3 has the channel divided into 5 
c          spectral intervals.  Comparisons with the old version shows only a 1-2% improvement.
c          I've decided to stick with the old version.  The new version would have required 
c          separate bnu arrays for each spectral interval & the run-time of the Channel 3 
c          portion of the code would have increased by a factor of 5.
c      6)  Kratz has added absorption by CFCs to channels 4 & 5.  This changes the COMMON 
c          block /kprofs/ which now includes the array "ucfc".  This array is filled in 
c          subroutine fillatmos_ckm.  The correlated k's are calculated in subroutines,
c          ck_cfc_4_ckm & ck_cfc_5_ckm.
c      7)  The window is now modeled using the Clough continuum, instead of Roberts.  This
c          uses a parameterization developed by Fred Rose.  The old Subroutine "window" has been
c          replaced by two channel-specific routines, window_4_ckm & window_5_ckm.
c      8)  We no longer calculate the total optical depth (totaltau).  This was an unnecessary 
c          and time-consuming step that I did originally to check if the numbers looked correct.  
c          This led to some of the problems that we had with NAN results when we tried to take the
c          logarithm of 0. or a negative number.
c          We now just calculate the attenuation.  This causes a change in Subroutines UseCorrK_ckm  & GetSurtTemp
c          since we will no longer need to take the EXP of the optical depth.  It also changes the common
c          block /ktaus/.
c      9)  The layer indexing in the code was changed slightly from Kratz's original.
c          Now, level 1 = surface, level 2 = the first level above the surface, etc. up to level klevs
c          And, layer 1 = the layer between levels 1 & 2, etc. up to layer (klevs - 1) = nlay
c          Kratz's old numbering system was different in the upward and downward directions.
c     10)  I've changed the cos(sz) & cos(vz) into sec(sz) & sec(vz) 
c     11)  I now check to see if we are dividing by 0. when calculating the layer attenuation.
c          This should eliminate NAN problems when there is a temperature inversion.
c  
c
c	It is assumed that Subroutine fillatmos_ckm has been called to fill in values of bnu for
c    each atmospheric layer and bnu0 for the surface radiance.
c
c  INPUT PARAMETERS:	
c		vz          = viewing zenith
c		sz          = solar zenith angle
c       surfemit3   = VIRS Channel 3 (3.75 µm) surface emittance
c       surfemit4   = VIRS Channel 4 (10.8 µm) surface emittance
c       surfemit5   = VIRS Channel 5 (11.9 µm) surface emittance
c       bnu(ilay)   = emitted radiance from layer ilay
c       bnu0(ichan) = emitted radiance for the surface for channel ichan (should be planck(SkinTemp))
c                     where ichan = 1 for Channel 3
c                           ichan = 2 for Channel 4
c                           ichan = 3 for Channel 5
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      parameter (num_segment_3 = 5)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real dis

      real sx4(MAXCKLEV),tsx(MAXCKLEV)
c     real sx3(MAXCKLEV)
      real sxz3(MAXCKLEV)
      real tsxb(MAXCKLEV)
c     real tsxa(MAXCKLEV)
      real tsyb(MAXCKLEV),tsya(MAXCKLEV)
      real sy3(MAXCKLEV),sz3(MAXCKLEV)
      real sx5(MAXCKLEV)

      real tsxb_sum(MAXCKLEV),weight_sum,bnu_sum(MAXCKLEV)
      real tsyb_sum(MAXCKLEV),tsya_sum(MAXCKLEV)
      real ramda_3(num_segment_3), delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)

      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      integer n_h2o_3(num_segment_3),n_ch4_3(num_segment_3),n_h2o_4,n_co2_4,n_h2o_5
      integer n_n2o_3(num_segment_3)
      real f_h2o_3(max_coefs,num_segment_3), f_h2o_4(max_coefs), f_h2o_5(max_coefs)
      real f_ch4_3(max_coefs,num_segment_3), f_co2_4(max_coefs)
      real f_n2o_3(max_coefs,num_segment_3)
      real f_o3_4(max_coefs), f_co2_5(max_coefs)
      real k_h2o_3(max_coefs,num_segment_3), k_h2o_4(max_coefs), k_h2o_5(max_coefs)
      real k_ch4_3(max_coefs,num_segment_3), k_co2_4(max_coefs)
      real k_n2o_3(max_coefs,num_segment_3)
      real k_o3_4(max_coefs), k_co2_5(max_coefs)
      real coefk_h2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_h2o_4(max_coefs,iorder,num_pressure)
      real coefk_h2o_5(max_coefs,iorder,num_pressure)
      real coefk_ch4_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_n2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_4(max_coefs,iorder,num_pressure)
      real coefk_o3_4(max_coefs,iorder,num_pressure)
      real coefk_co2_5(max_coefs,iorder,num_pressure)
      real tau_h2o_3(max_coefs,MAXCKLEV)
      real tau_h2o_4(max_coefs,MAXCKLEV)
      real tau_h2o_5(max_coefs,MAXCKLEV)
      real tau_ch4_3(max_coefs,MAXCKLEV)
      real tau_n2o_3(max_coefs,MAXCKLEV)
      real tau_co2_4(max_coefs,MAXCKLEV)
      real tau_o3_4(max_coefs,MAXCKLEV)
      real tau_co2_5(max_coefs,MAXCKLEV)
      real tsxb_up(MAXCKLEV), tsxb_up_sum(MAXCKLEV)

      common/corrkcoef/ n_h2o_3, f_h2o_3, k_h2o_3, coefk_h2o_3,
     +  n_ch4_3, f_ch4_3, k_ch4_3, coefk_ch4_3,
     +  n_n2o_3, f_n2o_3, k_n2o_3, coefk_n2o_3,
     +  n_h2o_4, f_h2o_4, k_h2o_4, coefk_h2o_4,
     +  n_co2_4, f_co2_4, k_co2_4, coefk_co2_4,
     +  n_h2o_5, f_h2o_5, k_h2o_5, coefk_h2o_5,
     +  n_co2_5, f_co2_5, k_co2_5, coefk_co2_5

      real tauf_3(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_4(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_5(max_coefs,max_coefs,MAXCKLEV)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),
     +  solarrad(MAXCKLEV),ratio_up(MAXCKLEV)
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg

      i_goes_segment=1

      nlay = klevs - 1               ! nlay = the number of layers (= number of levels - 1)
      secvz = 1. / COS(.01745329*vz)
      if(sz .lt. twilit) then
        secsz = 1. / COS(.01745329*sz)
        cossz = COS(.01745329*sz)
      else
        cossz = 0.
        secsz = 0.
      end if
c
      ichan = 1
      bsun = planck_ckm(ichan,sc3temp)
      bsun = bsun * dis
	  
C   -----------------------------------------------------------------
C    First, we do channel 3 downward in the solar direction and
C    upward in the viewing direction
C   -----------------------------------------------------------------
	  
      surfref = 1. - surfemit3
      if (surfref.le.0.02) surfref=0.02
c     bref = bsun * surfref * cossz
c**	  write(28,*) bref,surfref

c     weight_sum = 0.0
      do m = 1,klevs
        tsxb_sum(m) = 0.0
        tsya_sum(m) = 0.0
        tsyb_sum(m) = 0.0
        bnu_sum(m) = 0.0
        tsxb_up_sum(m) = 0.0
      end do

      do iseg = 1,i_goes_segment
 
      do m = 1,klevs
        tsx(m) = 0.0
c       tsxa(m) = 0.0
        tsxb(m) = 0.0
        tsya(m) = 0.0
        tsyb(m) = 0.0
        tsxb_up(m) = 0.0
      end do
	  
      tb = 0.
	   
c     back = bnu0(1) * surfemit3

      call ck_ckm(n_h2o_3(iseg),nlay,ulay,f_h2o_3(1,iseg),k_h2o_3(1,iseg),
     +  coefk_h2o_3(1,1,1,iseg),play,tlay,tau_h2o_3)
      call ck_ckm(n_ch4_3(iseg),nlay,uch4,f_ch4_3(1,iseg),k_ch4_3(1,iseg),
     +  coefk_ch4_3(1,1,1,iseg),play,tlay,tau_ch4_3)
      call ck_ckm(n_n2o_3(iseg),nlay,un2o,f_n2o_3(1,iseg),k_n2o_3(1,iseg),
     +  coefk_n2o_3(1,1,1,iseg),play,tlay,tau_n2o_3)

      do m = 1, nlay
        do i = 1, n_h2o_3(iseg)
          do ia = 1, n_ch4_3(iseg)
            do ib = 1, n_n2o_3(iseg)
              tauf_3(i,ia,ib, m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m) + tau_n2o_3(ib,m)
c             tauf_3(i,ia,m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m)  
            end do
          end do
        end do
      end do	   
 
      sy3(klevs) = 0.
      if (i_goes_segment .eq. 1) then
        sz3(klevs) = bsun * cossz
      else
        sz3(klevs) = bsun_seg(iseg) * cossz
      end if
      do i = 1, n_h2o_3(iseg)
        do ia = 1, n_ch4_3(iseg)
          do ib = 1, n_n2o_3(iseg)
            tsya(klevs) = tsya(klevs) + sy3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            tsyb(klevs) = tsyb(klevs) + sz3(klevs) * f_h2o_3(i,iseg) *
     +        f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
            do m = klevs-1, 1, -1			! Solve Downward emitted (vz)
              if (i_goes_segment .eq. 1) then
                sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +            bnu(m,1) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              else
                sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +            bnu_5seg(m,iseg) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              end if
              tsya(m) = tsya(m) + sy3(m) * f_h2o_3(i,iseg) *
     +          f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)

              if(cossz .gt. 0.) then                     ! Do the downward solar component
                sz3(m) = sz3(m+1) * exp(-tauf_3(i,ia,ib,m) * secsz)
                tsyb(m) = tsyb(m) + sz3(m) * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)         

                sxz3_up = sz3(m)
                do m_up = m+1, klevs                      ! Do the reflected upward solar component (vz)
                  sxz3_up = sxz3_up * exp(-tauf_3(i,ia,ib,m_up-1) * secvz)
                end do
                tsxb_up(m) = tsxb_up(m) + sxz3_up * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_n2o_3(ib,iseg)
              end if
            end do

            if(cossz .le. 0.0) then  !added 12/2005
              sz3(1) = 0.0
            endif

            sxz3(1) = sz3(1) * surfref
            tsxb(1) = tsxb(1) + sxz3(1) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +        f_n2o_3(ib,iseg)
            do m = 2, klevs
              if(cossz .gt. 0.) then                       ! Do the upward solar component (vz)
                sxz3(m) = sxz3(m-1) * exp(-tauf_3(i,ia,ib,m-1) * secvz)
                tsxb(m) = tsxb(m) + sxz3(m) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) *
     +            f_n2o_3(ib,iseg) 
              end if
            end do
          end do
        end do
      end do

c     weight_sum = weight_sum + delw_3(iseg) * response_3(iseg)
      do m = 1,klevs
c       weight_sum(m) = weight_sum(m) + delw_3(iseg) * response_3(iseg)
        tsxb_sum(m) = tsxb_sum(m) + tsxb(m) * delw_3(iseg) * response_3(iseg)
        tsya_sum(m) = tsya_sum(m) + tsya(m) * delw_3(iseg) * response_3(iseg)
        tsyb_sum(m) = tsyb_sum(m) + tsyb(m) * delw_3(iseg) * response_3(iseg)
        tsxb_up_sum(m) = tsxb_up_sum(m) + tsxb_up(m) * delw_3(iseg) *
     +    response_3(iseg)
        if (i_goes_segment .eq. 1) then
          bnu_sum(m) = bnu_sum(m) + bnu(m,1) * delw_3(iseg) * response_3(iseg)
        else
          bnu_sum(m) = bnu_sum(m) + bnu_5seg(m,iseg) * delw_3(iseg) * response_3(iseg)
        end if
      end do
      end do

	  
c  Now let's calculate the attenuation ratios for each layer for the 
c  reflected part and a composite optical depth for the emission part

      solarrad(klevs) = bsun * cossz
      do m = nlay,1,-1
        if(tsya_sum(m) - bnu_sum(m) .ne. 0.) then
          attenuate(m,1) = (tsya_sum(m) - bnu_sum(m)) / 
     +      (tsya_sum(m+1) - bnu_sum(m))
        else
          attenuate(m,1) = 0.
        end if
        if(cossz .gt. 0.) then
          ratiorefl(m) = tsxb_sum(m+1) / tsxb_sum(m)
          ratio_up(m) = tsxb_up_sum(m) / tsyb_sum(m)
          solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1)
c         solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1) +      ! Include downward emission
c    +      tsya_sum(m)/weight_sum                  
        else
          ratiorefl(m) = 1.
          solarrad(m) = 0.
        end if
      end do
      if(cossz .gt. 0.) then
        TotalRefRat  = tsxb_sum(nlay + 1) / tsxb_sum(1)
      else
        TotalRefRat  = 1.
      end if

      Radorig = 0.
      do m = nlay, 1, -1                                         ! Correctly add the downwelling
c       solarrad(m) = solarrad(m) + tsya_sum(m) / weight_sum
        iupdown = 2
        ichan = 1
        laylow = m
        layhigh = m
        call UseCorrK_down_ckm(RadOrig,laylow,layhigh,iupdown,ichan,vz,Radfinal)
        solarrad(m) = solarrad(m) + Radfinal
        RadOrig = Radfinal
      end do
 
C   -----------------------------------------------------------------
C    Next, we do channel 4 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.
      
c     back = bnu0(2) * surfemit4

      call ck_ckm(n_h2o_4,nlay,ulay,f_h2o_4,k_h2o_4,
     +  coefk_h2o_4,play,tlay,tau_h2o_4)
      call ckco2_goes_chan4_ckm(n_co2_4,nlay,uco2,f_co2_4,
     +  play,tlay,tau_co2_4)
c     call cko3_goes_chan4_ckm(n_o3_4,nlay,uo3,f_o3_4,play,tlay,tau_o3_4)  ! Don't use O3
      n_o3_4 = 1                              ! Don't use O3
      f_o3_4(1) = 1.0                         ! Don't use O3

      do m = 1, nlay
        call window_goes_chan4_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i = 1, n_h2o_4
          do ia = 1, n_co2_4
            do ib = 1, n_o3_4
	      tau_o3_4(ib,m)=0.0   ! Don't use O3
              tauf_4(i,ia,ib,m) = tau_h2o_4(i,m) + twin + tau_co2_4(ia,m) +
     &          tau_o3_4(ib,m)
c                print*,'10.8um:',m,tau_h2o_4(i,m),tau_co2_4(ia,m),twin
            end do
          end do
        end do      
      end do
 
      sx4(klevs) = 0.0
      do i = 1, n_h2o_4
        do ia = 1, n_co2_4
          do ib = 1, n_o3_4
            do m = klevs-1,1,-1
              sx4(m) = sx4(m+1) * exp(-tauf_4(i,ia,ib,m) * 
     *          secvz) + bnu(m,2) * (1.0 - 
     *          exp(-tauf_4(i,ia,ib,m) * secvz))
              tsx(m) = tsx(m) + sx4(m) * f_h2o_4(i) * 
     *          f_co2_4(ia) * f_o3_4(ib)
            end do
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth
	  
      do m = 1,nlay
        if(tsx(m+1) - bnu(m,2) .ne. 0.) then
	  attenuate(m,2) = (tsx(m) - bnu(m,2)) / (tsx(m+1) - 
     *      bnu(m,2))
c          print*,m,attenuate(m,2)
        else
          attenuate(m,2) = 0.
        end if
c       if (attenuate(m,2).lt.0.0.or.attenuate(m,2).gt.2.0) then
c         mmmm=1
c       end if
      end do
	   
C   -----------------------------------------------------------------
C    Now, we do channel 5 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.

c     back = bnu0(3) * surfemit5

    
      call ck_ckm(n_h2o_5,nlay,ulay,f_h2o_5,k_h2o_5,
     +  coefk_h2o_5,play,tlay,tau_h2o_5)
      call ckco2_goes_chan5_ckm(n_co2_5,nlay,uco2,f_co2_5,play,
     +  tlay,tau_co2_5)

      
      do m=1, nlay
        call window_goes_chan5_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i=1, n_h2o_5
          do ib=1, n_co2_5
            tauf_5(i,ib,m) = tau_h2o_5(i,m) + twin + tau_co2_5(ib,m)
c	     print*,'12.0um:',i,tau_h2o_5(i,m),tau_co2_5(ib,m),twin
          end do
        end do      
      end do
 
     
      sx5(klevs) = 0.0
      do i = 1, n_h2o_5
        do ib = 1, n_co2_5
          do m = klevs-1,1,-1
            sx5(m) = sx5(m+1) * exp(-tauf_5(i,ib,m) * 
     *        secvz) + bnu(m,3) * (1.0 - 
     *        exp(-tauf_5(i,ib,m) * secvz))
            tsx(m) = tsx(m) + sx5(m) * f_h2o_5(i) * f_co2_5(ib)
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth

      do m = 1,nlay
        if(tsx(m+1) - bnu(m,3) .ne. 0.) then
	  attenuate(m,3) = (tsx(m) - bnu(m,3)) / (tsx(m+1) - 
     +      bnu(m,3))
        else
          attenuate(m,3) = 0.
        end if
      end do
	  
 100  format(3i10,5f10.3)
 
      return 
      end
      
c *********************************************************************
      subroutine ckco2_goes_chan4_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.827729*0.01929046*(1.0000+3.5973e-02*(t(m)-250.0)+
     *       4.08281e-04*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
         
c *********************************************************************
      subroutine cko3_goes_chan4_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.630813*24.49382*(1.0000+1.6825e-02*(t(m)-250.0)+
     *       8.7256e-05*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end

c**********************************************************************
      subroutine window_goes_chan4_ckm(amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 870 to 1000 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
      data aa/9.227,0.9998,-0.02426,0.01972,-0.7891,5.843e-05,0.9922/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

      twin = exp ( tau_log)
c
      return
      end
               
c *********************************************************************
      subroutine ckco2_goes_chan5_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.705998*0.02862606*(1.0000+3.2154e-02*(t(m)-250.0)+
     *       3.3047e-04*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
         
c**********************************************************************
      subroutine window_goes_chan5_ckm (amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 790 to 880 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
      data aa/9.651,1.000,-0.02445,0.03456,-0.6236,-2.403e-04,0.9788/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
          tau_log   = aa(1)              +
     *                aa(2) * log(amnt)  +
     *                aa(3) * temp       +
     *                aa(4) * patm       +
     *                aa(5) * ph2o       +
     *                aa(6) * amnt       +
     *                aa(7) * log(ph2o)

      twin = exp ( tau_log)
c
      return
      end
c*****************************************************************************************
      Subroutine CorrKsetup_modis_ckm(vz,sz,dis,surfemit3,surfemit4,surfemit5)
c
c	This subroutine sets up the correlated Ks.
c
c	Kratz's code is used to calculate upward flux at each level.  I then solve for
c	  a total optical depth for each layer.  These are then used in routine "UseCorrK_ckm".
c     These optical depths are only valid for the given angle.
c
c   The major changes in this version:
c      1)  The surface emittances are used (these are now input parameters)
c      2)  Channel 3 surface reflectance = 1. - surfemit3 instead of a constant value of 0.10
c      3)  For channel 3, we now solve downward first.  We then use the solar radiance at
c          the surface for bref, instead of using the radiance at the TOA
c      4)  We now use subroutine ck for all correlated k calculations (except for those
c          that involve CFCs).  The inputs to subroutine ck have changed to include the 
c          coefficients.  These coefficients are now read from an ancillary data file during
c          the first call to CorrKsetup_ckm
c      5)  Kratz's latest edition of correlated ks for Channel 3 has the channel divided into 5 
c          spectral intervals.  Comparisons with the old version shows only a 1-2% improvement.
c          I've decided to stick with the old version.  The new version would have required 
c          separate bnu arrays for each spectral interval & the run-time of the Channel 3 
c          portion of the code would have increased by a factor of 5.
c      6)  Kratz has added absorption by CFCs to channels 4 & 5.  This changes the COMMON 
c          block /kprofs/ which now includes the array "ucfc".  This array is filled in 
c          subroutine fillatmos_ckm.  The correlated k's are calculated in subroutines,
c          ck_cfc_4_ckm & ck_cfc_5_ckm.
c      7)  The window is now modeled using the Clough continuum, instead of Roberts.  This
c          uses a parameterization developed by Fred Rose.  The old Subroutine "window" has been
c          replaced by two channel-specific routines, window_4_ckm & window_5_ckm.
c      8)  We no longer calculate the total optical depth (totaltau).  This was an unnecessary 
c          and time-consuming step that I did originally to check if the numbers looked correct.  
c          This led to some of the problems that we had with NAN results when we tried to take the
c          logarithm of 0. or a negative number.
c          We now just calculate the attenuation.  This causes a change in Subroutines UseCorrK_ckm  & GetSurtTemp
c          since we will no longer need to take the EXP of the optical depth.  It also changes the common
c          block /ktaus/.
c      9)  The layer indexing in the code was changed slightly from Kratz's original.
c          Now, level 1 = surface, level 2 = the first level above the surface, etc. up to level klevs
c          And, layer 1 = the layer between levels 1 & 2, etc. up to layer (klevs - 1) = nlay
c          Kratz's old numbering system was different in the upward and downward directions.
c     10)  I've changed the cos(sz) & cos(vz) into sec(sz) & sec(vz) 
c     11)  I now check to see if we are dividing by 0. when calculating the layer attenuation.
c          This should eliminate NAN problems when there is a temperature inversion.
c  
c
c	It is assumed that Subroutine fillatmos_ckm has been called to fill in values of bnu for
c    each atmospheric layer and bnu0 for the surface radiance.
c
c  INPUT PARAMETERS:	
c		vz          = viewing zenith
c		sz          = solar zenith angle
c       surfemit3   = VIRS Channel 3 (3.75 µm) surface emittance
c       surfemit4   = VIRS Channel 4 (10.8 µm) surface emittance
c       surfemit5   = VIRS Channel 5 (11.9 µm) surface emittance
c       bnu(ilay)   = emitted radiance from layer ilay
c       bnu0(ichan) = emitted radiance for the surface for channel ichan (should be planck(SkinTemp))
c                     where ichan = 1 for Channel 3
c                           ichan = 2 for Channel 4
c                           ichan = 3 for Channel 5
c
	  
      include './corrk.h'
      integer MaxChanCK
      parameter (MaxChanCK=4)
      parameter (num_segment_3 = 5)
      real pk(MAXCKLEV),tk(MAXCKLEV),uk(MAXCKLEV)
      real zk(MAXCKLEV),rhk(MAXCKLEV),o3k(MAXCKLEV)
      real play(MAXCKLEV),tlay(MAXCKLEV)
      real delplay(MAXCKLEV),delzlay(MAXCKLEV)
      real ulay(MAXCKLEV),uco2(MAXCKLEV),uch4(MAXCKLEV)
      real uf12(MAXCKLEV),uf11(MAXCKLEV),un2o(MAXCKLEV)
      real o3lay(MAXCKLEV)
      real bnu(MAXCKLEV,MaxChanCK),bnu0(MaxChanCK)
      real dis

      real sx4(MAXCKLEV),tsx(MAXCKLEV)
c     real sx3(MAXCKLEV)
      real sxz3(MAXCKLEV)
      real tsxb(MAXCKLEV)
c     real tsxa(MAXCKLEV)
      real tsyb(MAXCKLEV),tsya(MAXCKLEV)
      real sy3(MAXCKLEV),sz3(MAXCKLEV)
      real sx5(MAXCKLEV)

      real tsxb_sum(MAXCKLEV),weight_sum,bnu_sum(MAXCKLEV)
      real tsyb_sum(MAXCKLEV),tsya_sum(MAXCKLEV)
      real ramda_3(num_segment_3), delw_3(num_segment_3), response_3(num_segment_3)
      real bnu_5seg(MAXCKLEV,num_segment_3), bsun_seg(num_segment_3)

      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)

      integer n_h2o_3(num_segment_3),n_ch4_3(num_segment_3),n_h2o_4,n_co2_4,n_h2o_5
      integer n_co2_3(num_segment_3)
      real f_h2o_3(max_coefs,num_segment_3), f_h2o_4(max_coefs), f_h2o_5(max_coefs)
      real f_ch4_3(max_coefs,num_segment_3), f_co2_4(max_coefs)
      real f_co2_3(max_coefs,num_segment_3)
      real f_o3_4(max_coefs), f_co2_5(max_coefs)
      real k_h2o_3(max_coefs,num_segment_3), k_h2o_4(max_coefs), k_h2o_5(max_coefs)
      real k_ch4_3(max_coefs,num_segment_3), k_co2_4(max_coefs)
      real k_co2_3(max_coefs,num_segment_3)
      real k_o3_4(max_coefs), k_co2_5(max_coefs)
      real coefk_h2o_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_h2o_4(max_coefs,iorder,num_pressure)
      real coefk_h2o_5(max_coefs,iorder,num_pressure)
      real coefk_ch4_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_3(max_coefs,iorder,num_pressure,num_segment_3)
      real coefk_co2_4(max_coefs,iorder,num_pressure)
      real coefk_o3_4(max_coefs,iorder,num_pressure)
      real coefk_co2_5(max_coefs,iorder,num_pressure)
      real tau_h2o_3(max_coefs,MAXCKLEV)
      real tau_h2o_4(max_coefs,MAXCKLEV)
      real tau_h2o_5(max_coefs,MAXCKLEV)
      real tau_ch4_3(max_coefs,MAXCKLEV)
      real tau_co2_3(max_coefs,MAXCKLEV)
      real tau_co2_4(max_coefs,MAXCKLEV)
      real tau_o3_4(max_coefs,MAXCKLEV)
      real tau_co2_5(max_coefs,MAXCKLEV)
      real tsxb_up(MAXCKLEV), tsxb_up_sum(MAXCKLEV)

      integer n_h2o_85, n_ch4_85, n_o3_85, n_n2o_85
      real f_h2o_85(max_coefs), f_ch4_85(max_coefs), f_o3_85(max_coefs), 
     +  f_n2o_85(max_coefs)
      real sx85(MAXCKLEV)
      real tau_h2o_85(max_coefs,MAXCKLEV)
      real tau_ch4_85(max_coefs,MAXCKLEV)
      real tau_o3_85(max_coefs,MAXCKLEV)
      real tau_n2o_85(max_coefs,MAXCKLEV)
      real tauf_85(max_coefs,max_coefs,max_coefs,max_coefs,MAXCKLEV)
c     real attenuate_85(MAXCKLEV), bnu_85(MAXCKLEV), bnu0_85
c     common/ktau_85/attenuate_85, bnu_85, bnu0_85
***
***  WARNING, changed the 3rd line of corrkcoef common block
***   to account for differences between goes and modis.
***   Used n_co2_3, etc., in place of n_n2o_3, etc.
***   n2o is used for goes and for avhrr & virs.
***
      common/corrkcoef/ n_h2o_3, f_h2o_3, k_h2o_3, coefk_h2o_3,
     +  n_ch4_3, f_ch4_3, k_ch4_3, coefk_ch4_3,
     +  n_co2_3, f_co2_3, k_co2_3, coefk_co2_3,
     +  n_h2o_4, f_h2o_4, k_h2o_4, coefk_h2o_4,
     +  n_co2_4, f_co2_4, k_co2_4, coefk_co2_4,
     +  n_h2o_5, f_h2o_5, k_h2o_5, coefk_h2o_5,
     +  n_co2_5, f_co2_5, k_co2_5, coefk_co2_5

      real tauf_3(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_4(max_coefs,max_coefs,max_coefs,MAXCKLEV)
      real tauf_5(max_coefs,max_coefs,MAXCKLEV)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),
     +  solarrad(MAXCKLEV),ratio_up(MAXCKLEV)
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up
	  
      COMMON /ktemps/ klevs,pk,tk,zk,uk,rhk,o3k,i85flag
      COMMON /kprofs/ play,tlay,delplay,delzlay,ulay,uco2,uch4,uf12,uf11,un2o,o3lay
      COMMON /krads/ bnu,bnu0
      COMMON /constants/ sc3temp,bsun,solzenmax,twilit
      common /Chan3_5segments/ ramda_3, delw_3, response_3, bnu_5seg, bsun_seg

      i_modis_segment=1

      nlay = klevs - 1               ! nlay = the number of layers (= number of levels - 1)
      secvz = 1. / COS(.01745329*vz)
      if(sz .lt. twilit) then
        secsz = 1. / COS(.01745329*sz)
        cossz = COS(.01745329*sz)
      else
        cossz = 0.
        secsz = 0.
      end if
c
      ichan = 1
      bsun = planck_ckm(ichan,sc3temp)
      bsun = bsun * dis
	  
C   -----------------------------------------------------------------
C    First, we do channel 20 downward in the solar direction and
C    upward in the viewing direction
C   -----------------------------------------------------------------
	  
      surfref = 1. - surfemit3
      if (surfref.le.0.02) surfref=0.02
c     bref = bsun * surfref * cossz

c     weight_sum = 0.0
      do m = 1,klevs
        tsxb_sum(m) = 0.0
        tsya_sum(m) = 0.0
        tsyb_sum(m) = 0.0
        bnu_sum(m) = 0.0
        tsxb_up_sum(m) = 0.0
      end do

      do iseg = 1,i_modis_segment
 
      do m = 1,klevs
        tsx(m) = 0.0
c       tsxa(m) = 0.0
        tsxb(m) = 0.0
        tsya(m) = 0.0
        tsyb(m) = 0.0
        tsxb_up(m) = 0.0
      end do
	  
      tb = 0.
	   
c     back = bnu0(1) * surfemit3

      call ck_ckm(n_h2o_3(iseg),nlay,ulay,f_h2o_3(1,iseg),k_h2o_3(1,iseg),
     +  coefk_h2o_3(1,1,1,iseg),play,tlay,tau_h2o_3)
      call ckco2_modis_chan20_ckm(n_co2_3(iseg),nlay,uco2,f_co2_3(1,iseg),
     +  play,tlay,tau_co2_3)
      call ckch4_modis_chan20_ckm(n_ch4_3(iseg),nlay,uch4,f_ch4_3(1,iseg),
     +  play,tlay,tau_ch4_3)
      do m = 1, nlay
        do i = 1, n_h2o_3(iseg)
          do ia = 1, n_ch4_3(iseg)
            do ib = 1, n_co2_3(iseg)
              tauf_3(i,ia,ib, m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m) + 
     +          tau_co2_3(ib,m)                                            ! changed 10/2000
c****         tauf_3(i,ia,ib, m) = tau_h2o_3(i,m) + tau_ch4_3(ia,m)        ! to test w/o co2
            end do
          end do
        end do
      end do	   
 
      sy3(klevs) = 0.
      if (i_modis_segment .eq. 1) then
        sz3(klevs) = bsun * cossz
      else
        sz3(klevs) = bsun_seg(iseg) * cossz
      end if
      do i = 1, n_h2o_3(iseg)
        do ia = 1, n_ch4_3(iseg)
          do ib = 1, n_co2_3(iseg)
            tsya(klevs) = tsya(klevs) + sy3(klevs) * f_h2o_3(i,iseg) * 
     +        f_ch4_3(ia,iseg) * f_co2_3(ib,iseg)
            tsyb(klevs) = tsyb(klevs) + sz3(klevs) * f_h2o_3(i,iseg) * 
     +        f_ch4_3(ia,iseg) * f_co2_3(ib,iseg)
            do m = klevs-1, 1, -1			! Solve Downward emitted (vz)
              if (i_modis_segment .eq. 1) then
                sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +            bnu(m,1) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              else
                sy3(m) = sy3(m+1) * exp(-tauf_3(i,ia,ib,m) * secvz)+
     +            bnu_5seg(m,iseg) * (1.0 - exp(-tauf_3(i,ia,ib,m) * secvz))
              end if
              tsya(m) = tsya(m) + sy3(m) * f_h2o_3(i,iseg) * 
     +          f_ch4_3(ia,iseg) * f_co2_3(ib,iseg)             ! changed 10/2000

              if(cossz .gt. 0.) then                       ! Do the downward solar component
                sz3(m) = sz3(m+1) * exp(-tauf_3(i,ia,ib,m) * secsz)
                tsyb(m) = tsyb(m) + sz3(m) * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_co2_3(ib,iseg)           ! changed 10/2000

                sxz3_up = sz3(m)
                do m_up = m+1, klevs                      ! Do the reflected upward solar component (vz)
                  sxz3_up = sxz3_up * exp(-tauf_3(i,ia,ib,m_up-1) * secvz)
                end do
                tsxb_up(m) = tsxb_up(m) + sxz3_up * f_h2o_3(i,iseg) * 
     +            f_ch4_3(ia,iseg) * f_co2_3(ib,iseg)             ! added 10/2000
              end if
            end do

            if(cossz .le. 0.0) then !added 12/2005
              sz3(1) = 0.0
            endif

            sxz3(1) = sz3(1) * surfref
            tsxb(1) = tsxb(1) + sxz3(1) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) * 
     +        f_co2_3(ib,iseg)                       ! changed 10/2000
            do m = 2, klevs
              if(cossz .gt. 0.) then                       ! Do the upward solar component (vz)
                sxz3(m) = sxz3(m-1) * exp(-tauf_3(i,ia,ib,m-1) * secvz)
                tsxb(m) = tsxb(m) + sxz3(m) * f_h2o_3(i,iseg) * f_ch4_3(ia,iseg) * 
     +            f_co2_3(ib,iseg)           ! changed 10/2000
              end if
            end do
          end do
        end do
      end do

c     weight_sum = weight_sum + delw_3(iseg) * response_3(iseg)
      do m = 1,klevs
c       weight_sum(m) = weight_sum(m) + delw_3(iseg) * response_3(iseg)
        tsxb_sum(m) = tsxb_sum(m) + tsxb(m) * delw_3(iseg) * response_3(iseg)
        tsya_sum(m) = tsya_sum(m) + tsya(m) * delw_3(iseg) * response_3(iseg)
        tsyb_sum(m) = tsyb_sum(m) + tsyb(m) * delw_3(iseg) * response_3(iseg)
        tsxb_up_sum(m) = tsxb_up_sum(m) + tsxb_up(m) * delw_3(iseg) *
     +    response_3(iseg)
        if (i_modis_segment .eq. 1) then
          bnu_sum(m) = bnu_sum(m) + bnu(m,1) * delw_3(iseg) * response_3(iseg)
        else
          bnu_sum(m) = bnu_sum(m) + bnu_5seg(m,iseg) * delw_3(iseg) * response_3(iseg)
        end if
      end do
      end do

	  
c  Now let's calculate the attenuation ratios for each layer for the 
c  reflected part and a composite optical depth for the emission part

      solarrad(klevs) = bsun * cossz
      do m = nlay,1,-1
        if(tsya_sum(m) - bnu_sum(m) .ne. 0.) then
          attenuate(m,1) = (tsya_sum(m) - bnu_sum(m)) / 
     +      (tsya_sum(m+1) - bnu_sum(m))
        else
          attenuate(m,1) = 0.
        end if
        if(cossz .gt. 0.) then
          ratiorefl(m) = tsxb_sum(m+1) / tsxb_sum(m)
          ratio_up(m) = tsxb_up_sum(m) / tsyb_sum(m)
          solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1)
c         solarrad(m) = solarrad(m+1) * tsyb_sum(m) / tsyb_sum(m+1) +      ! Include downward emission
c    +      tsya_sum(m)/weight_sum                  
        else
          ratiorefl(m) = 1.
          solarrad(m) = 0.
        end if
      end do
      if(cossz .gt. 0.) then
        TotalRefRat  = tsxb_sum(nlay + 1) / tsxb_sum(1)
      else
        TotalRefRat  = 1.
      end if

      Radorig = 0.
      do m = nlay, 1, -1                                         ! Correctly add the downwelling
c       solarrad(m) = solarrad(m) + tsya_sum(m) / weight_sum
        iupdown = 2
        ichan = 1
        laylow = m
        layhigh = m
        call UseCorrK_down_ckm(RadOrig,laylow,layhigh,iupdown,ichan,vz,Radfinal)
        solarrad(m) = solarrad(m) + Radfinal
        RadOrig = Radfinal
      end do
 
C   -----------------------------------------------------------------
C    Next, we do channel 31 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.
      
c     back = bnu0(2) * surfemit4

      call ck_ckm(n_h2o_4,nlay,ulay,f_h2o_4,k_h2o_4,
     +  coefk_h2o_4,play,tlay,tau_h2o_4)
      call ckco2_modis_chan31_ckm(n_co2_4,nlay,uco2,f_co2_4,
     +  play,tlay,tau_co2_4)
c     call cko3_goes_chan4_ckm(n_o3_4,nlay,uo3,f_o3_4,play,tlay,tau_o3_4)  ! Don't use O3
      n_o3_4 = 1                              ! Don't use O3
      f_o3_4(1) = 1.0                         ! Don't use O3

      do m = 1, nlay
        call window_modis_chan31_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i = 1, n_h2o_4
          do ia = 1, n_co2_4
            do ib = 1, n_o3_4
	      tau_o3_4(ib,m)=0.0 ! Don't use O3
              tauf_4(i,ia,ib,m) = tau_h2o_4(i,m) + twin + tau_co2_4(ia,m) +
     &           tau_o3_4(ib,m)
            end do
          end do
        end do      
      end do
 
      sx4(klevs) = 0.0
      do i = 1, n_h2o_4
        do ia = 1, n_co2_4
          do ib = 1, n_o3_4
            do m = klevs-1,1,-1
              sx4(m) = sx4(m+1) * exp(-tauf_4(i,ia,ib,m) * 
     *          secvz) + bnu(m,2) * (1.0 - 
     *          exp(-tauf_4(i,ia,ib,m) * secvz))
              tsx(m) = tsx(m) + sx4(m) * f_h2o_4(i) * 
     *          f_co2_4(ia) * f_o3_4(ib)
            end do
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth
	  
      do m = 1,nlay
        if(tsx(m+1) - bnu(m,2) .ne. 0.) then
	  attenuate(m,2) = (tsx(m) - bnu(m,2)) / (tsx(m+1) - 
     *      bnu(m,2))
        else
          attenuate(m,2) = 0.
        end if
c       print 88,m,tsx(m),attenuate(m,2)
c       if (attenuate(m,2).lt.0.0.or.attenuate(m,2).gt.2.0) then
c         mmmm=1
c       end if
      end do
	   
C   -----------------------------------------------------------------
C    Now, we do channel 32 for the viewing direction only
C   -----------------------------------------------------------------
	  
      do m = 1,klevs
        tsx(m) = 0.0
      end do
	  
      tb = 0.

c     back = bnu0(3) * surfemit5

      call ck_ckm(n_h2o_5,nlay,ulay,f_h2o_5,k_h2o_5,
     +  coefk_h2o_5,play,tlay,tau_h2o_5)
      call ckco2_modis_chan32_ckm(n_co2_5,nlay,uco2,f_co2_5,play,
     +  tlay,tau_co2_5)

      do m=1, nlay
        call window_modis_chan32_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i=1, n_h2o_5
          do ib=1, n_co2_5
            tauf_5(i,ib,m) = tau_h2o_5(i,m) + twin + tau_co2_5(ib,m)
          end do
        end do      
      end do
 
      sx5(klevs) = 0.0
      do i = 1, n_h2o_5
        do ib = 1, n_co2_5
          do m = klevs-1,1,-1
            sx5(m) = sx5(m+1) * exp(-tauf_5(i,ib,m) * 
     *        secvz) + bnu(m,3) * (1.0 - 
     *        exp(-tauf_5(i,ib,m) * secvz))
            tsx(m) = tsx(m) + sx5(m) * f_h2o_5(i) * f_co2_5(ib)
          end do
        end do
      end do
	  
c  Now let's calculate a composite optical depth

      do m = 1,nlay
        if(tsx(m+1) - bnu(m,3) .ne. 0.) then
	  attenuate(m,3) = (tsx(m) - bnu(m,3)) / (tsx(m+1) - 
     +      bnu(m,3))
        else
          attenuate(m,3) = 0.
        end if
c        print*,m,tsx(m),tsx(m+1),bnu(m,3),attenuate(m,3)
      end do

      if (i85flag.eq.1) then

C   -----------------------------------------------------------------
C    Next, we do channel 29 (8.5 µm) for the viewing direction only
C   -----------------------------------------------------------------

      do m = 1,klevs
        tsx(m) = 0.0
      end do

      tb = 0.

      back = bnu0(4)
                       
      call ck_85_ckm(n_h2o_85,nlay,ulay,delw,f_h2o_85,play,tlay,tau_h2o_85)
      call cko3_85_ckm(n_o3_85,nlay,o3lay,delw,f_o3_85,play,tlay,tau_o3_85)         
      call ckch4_85_ckm(n_ch4_85,nlay,uch4,delw,f_ch4_85,play,tlay,tau_ch4_85)
      call ckn2o_85_ckm(n_n2o_85,nlay,un2o,delw,f_n2o_85,play,tlay,tau_n2o_85)

      n_o3_85 = 1

      f_o3_85(1) = 1.                                                         

      n_n2o_85 = 1

      f_n2o_85(1) = 1.

      n_ch4_85 = 1

      f_ch4_85(1) = 1.


      do m = 1, nlay
        call window_modis_85_ckm(ulay(m),twin,play(m),tlay(m),delzlay(m))
        do i = 1, n_h2o_85
          do ia = 1, n_o3_85
            do ib = 1, n_ch4_85
              do ic = 1, n_n2o_85
                tauf_85(i,ia,ib,ic,m) = tau_h2o_85(i,m) + twin + tau_ch4_85(ib,m) +
     +            tau_n2o_85(ic,m)  
                tauf_85(i,ia,ib,ic,m) = tau_h2o_85(i,m) + twin + tau_o3_85(ia,m) + 
     +            tau_ch4_85(ib,m) + tau_n2o_85(ic,m)                          
              end do
            end do
          end do
        end do    
      end do

      do i = 1, n_h2o_85
        do ia = 1, n_o3_85
          do ib = 1, n_ch4_85
            do ic = 1, n_n2o_85
              do m = klevs, 1, -1
                if(m.eq.klevs)then
                  sx85(m) = 0.
                else
                  sx85(m) = sx85(m+1) * exp(-tauf_85(i,ia,ib,ic,m) * secvz) +
     *                   bnu(m,4) * (1.0 - exp(-tauf_85(i,ia,ib,ic,m) * secvz))
                end if
                tsx(m) = tsx(m) + sx85(m) * f_h2o_85(i) * f_o3_85(ia) * 
     +            f_ch4_85(ib) * f_n2o_85(ic)
              end do
            end do
          end do
        end do
      end do

c  Now let's calculate a composite optical depth

      do m = 1,nlay
        if(tsx(m+1) - bnu(m,4) .ne. 0.) then
          attenuate(m,4) = (tsx(m) - bnu(m,4)) / (tsx(m+1) - bnu(m,4))
        else
          attenuate(m,4) = 0.
        end if
c       print 88,m,tsx(m),attenuate(m,4)
c       write(28,'(10h chan 29   ,i10,3f15.1,f10.6)') m, tsx(m), tsx(m+1), 
c    +    bnu_85(m), attenuate_85(m)
      end do

      end if

 88   format(i6,2f9.3)

      return 
      end

c *********************************************************************
c     co2 2605--2735 cm^{-1}
      subroutine ckco2_modis_chan20_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.0017949128*(1.0000-1.2487e-04*(t(m)-250.0)+
     *    3.4136e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
         
c *********************************************************************
c     ch4 2605--2735 cm^{-1}
      subroutine ckch4_modis_chan20_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=14.524085*(1.0000+1.9377e-03*(t(m)-250.0)-
     *    3.0119e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end

c *********************************************************************
      subroutine ckco2_modis_chan31_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.912609*0.01047837*(1.0000+4.1069e-02*(t(m)-250.0)+
     *    5.1759e-04*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
         
c**********************************************************************
      subroutine window_modis_chan31_ckm(amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 835 to 980 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
c
      data aa/9.374,0.9999,-0.02441,0.01989,-0.7838,4.191e-05,0.9918/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
      tau_log   = aa(1)              +
     *            aa(2) * log(amnt)  +
     *            aa(3) * temp       +
     *            aa(4) * patm       +
     *            aa(5) * ph2o       +
     *            aa(6) * amnt       +
     *            aa(7) * log(ph2o)

      twin = exp ( tau_log )
c
      return
      end
         
c *********************************************************************
      subroutine ckco2_modis_chan32_ckm(ni,mlv,u,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)

      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.942022*0.00534925618*(1.0000+4.0490e-02*(t(m)-250.0)+
     *    5.0470e-04*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end

c**********************************************************************
      subroutine window_modis_chan32_ckm(amnt,twin,patm,temp,fac)
c     Parameterization of CKD_2.1 continuum over 835 to 980 cm-1 band
c     Fit inaccurate for amnt < 0.01 g/cm*2;
c          however optical depths < 0.001.
c     INPUT:
c     amnt = h2O layer amount (g/cm**2)
c     patm= pressure (atm)
c     temp = temperature (k)
c     OUTPUT:
c     twin = parameterized CKD_2.1 optical depth over band
c
      parameter (ncoef=7)
      real aa(ncoef)
c
      data aa/9.656,1.000,-0.02453,0.03056,-0.6620,-1.826e-04,0.9820/
c
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
c
      tau_log   = aa(1)              +
     *            aa(2) * log(amnt)  +
     *            aa(3) * temp       +
     *            aa(4) * patm       +
     *            aa(5) * ph2o       +
     *            aa(6) * amnt       +
     *            aa(7) * log(ph2o)

      twin = exp ( tau_log )
c
      return
      end

c *********************************************************************
c     h2o 1140-1200 cm^{-1}
      subroutine ck_85_ckm(ni,mlv,u,delw,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder,num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)

      common /standpres/stp
       
      ni=5
      f(1)=0.624325
      f(2)=0.253393
      f(3)=0.089834
      f(4)=0.026964
      f(5)=0.005484
      k(1)=0.0072755
      do i=2,ni
        k(i)=8.0*k(i-1)
      end do
      data ( (coefk(1,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,
     *0.00000,0.00006,0.00065,0.00265,0.00730,0.01769,
     *0.04487,0.11014,0.24011,0.44539,0.77830,1.26324,
     *1.93516,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,1.275e-06,1.470e-05,5.463e-05,1.524e-04,3.462e-04,
     *8.261e-04,1.945e-03,4.299e-03,7.952e-03,1.399e-02,2.282e-02,
     *3.441e-02,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
     *0.000e+00,5.625e-09,1.125e-07,3.906e-07,1.178e-06,2.606e-06,
     *5.284e-06,1.178e-05,2.511e-05,4.654e-05,7.976e-05,1.279e-04,
     *1.867e-04/
      data ( (coefk(2,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.00000,0.00000,0.00000,0.00000,0.00001,0.00021,
     *0.00100,0.00301,0.00822,0.02028,0.04442,0.09274,
     *0.16324,0.26366,0.39942,0.58242,0.84505,1.18790,
     *1.68415,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,2.125e-07,4.725e-06,
     *2.071e-05,6.063e-05,1.574e-04,3.879e-04,8.413e-04,1.719e-03,
     *2.968e-03,4.783e-03,7.230e-03,1.045e-02,1.519e-02,2.140e-02,
     *2.962e-02,
     *0.000e+00,0.000e+00,0.000e+00,0.000e+00,3.125e-10,3.563e-08,
     *1.459e-07,4.406e-07,1.066e-06,2.516e-06,5.706e-06,1.075e-05,
     *1.748e-05,2.754e-05,4.141e-05,5.808e-05,8.307e-05,1.204e-04,
     *1.526e-04/
      data ( (coefk(3,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.00010,0.00012,0.00032,0.00085,0.00193,0.00458,
     *0.01031,0.02119,0.04015,0.06639,0.10163,0.14948,
     *0.21695,0.30768,0.42844,0.59514,0.84740,1.17004,
     *1.48726,
     *6.413e-06,7.337e-06,1.126e-05,2.530e-05,5.484e-05,1.101e-04,
     *2.316e-04,4.420e-04,7.960e-04,1.290e-03,1.956e-03,2.831e-03,
     *4.046e-03,5.697e-03,7.895e-03,1.074e-02,1.516e-02,2.099e-02,
     *2.736e-02,
     *1.053e-07,1.197e-07,1.322e-07,2.650e-07,5.916e-07,1.016e-06,
     *1.897e-06,3.206e-06,5.294e-06,8.391e-06,1.266e-05,1.748e-05,
     *2.388e-05,3.295e-05,4.598e-05,5.965e-05,7.798e-05,1.001e-04,
     *1.449e-04/
      data ( (coefk(4,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.02081,0.02124,0.02440,0.03128,0.03987,0.05173,
     *0.06697,0.08747,0.11456,0.15100,0.20218,0.27745,
     *0.38290,0.51788,0.65919,0.78981,0.93236,1.06208,
     *1.13645,
     *6.691e-04,6.936e-04,7.050e-04,8.089e-04,9.665e-04,1.168e-03,
     *1.433e-03,1.786e-03,2.247e-03,2.882e-03,3.805e-03,5.040e-03,
     *6.803e-03,9.243e-03,1.208e-02,1.451e-02,1.720e-02,1.948e-02,
     *2.062e-02,
     *7.641e-06,8.097e-06,7.569e-06,7.547e-06,8.419e-06,9.331e-06,
     *1.065e-05,1.229e-05,1.459e-05,1.789e-05,2.317e-05,2.783e-05,
     *3.367e-05,4.389e-05,6.509e-05,7.864e-05,9.405e-05,1.047e-04,
     *9.958e-05/
      data ( (coefk(5,jt,jp), jp = 1, 19), jt = 1, 3)/
     *1.76643,1.76732,1.77176,1.77534,1.78324,1.79220,
     *1.80674,1.82530,1.84127,1.84779,1.83103,1.77690,
     *1.69182,1.57332,1.43992,1.27381,1.10369,0.89454,
     *0.69900,
     *3.494e-02,3.492e-02,3.490e-02,3.487e-02,3.489e-02,3.491e-02,
     *3.488e-02,3.491e-02,3.484e-02,3.469e-02,3.423e-02,3.325e-02,
     *3.176e-02,2.964e-02,2.708e-02,2.392e-02,2.113e-02,1.741e-02,
     *1.365e-02,
     *1.993e-04,1.990e-04,1.970e-04,1.965e-04,1.955e-04,1.953e-04,
     *1.920e-04,1.894e-04,1.848e-04,1.809e-04,1.761e-04,1.714e-04,
     *1.643e-04,1.556e-04,1.381e-04,1.207e-04,1.094e-04,9.372e-05,
     *7.355e-05/

      do i=1,ni
      do m=1,mlv
        ml=1
        pmb=p(m)*1013.25
        if(pmb.lt.stp(1))then
          x1=coefk(i,1,1)+coefk(i,2,1)*(t(m)-250.0)
     *    +coefk(i,3,1)*(t(m)-250.0)**2
          fkg(i,m)=x1*pmb/stp(1)
        else if (pmb.gt.stp(19)) then
          x1=coefk(i,1,18)+coefk(i,2,18)*(t(m)-250.0)
     *    +coefk(i,3,18)*(t(m)-250.0)**2
          x2=coefk(i,1,19)+coefk(i,2,19)*(t(m)-250.0)
     *    +coefk(i,3,19)*(t(m)-250.0)**2
          fkg(i,m)=x1+(x2-x1)/(stp(19)-stp(18))
     *   *(pmb-stp(18))
        else
          do while(pmb.ge.stp(ml))
            ml=ml+1
          end do
          x1=coefk(i,1,ml-1)+coefk(i,2,ml-1)*(t(m)-250.0)
     *    + coefk(i,3,ml-1)*(t(m)-250.0)**2
          x2=coefk(i,1,ml)+coefk(i,2,ml)*(t(m)-250.0)
     *    + coefk(i,3,ml)*(t(m)-250.0)**2
          fkg(i,m)=x1+(x2-x1)/(stp(ml)-stp(ml-1))
     *    *(pmb-stp(ml-1))
        end if
      end do
      end do
      do i=1,ni
        do m=1,mlv
          tau(i,m)=k(i)*u(m)*fkg(i,m)
        end do
      end do
      return
      end
c**********************************************************************
      subroutine window_modis_85_ckm(u,twin,pres,temp,fac)
      omega = 1170.0
      r=8.314d+07
c     amnt=u/6.023e+23*2.687e+19*18.01534
      amnt=u
c     write(6,*)u,amnt
c     a1=4.18+5577.8*exp(-.00787*omega)
      a1=4.18+7815.6*exp(-.0083*omega)
      a1=a1*exp(1800.*(1./temp-1./296.))
      a2=0.0008*a1
      ph2o=amnt/(fac*1.0d+05)*r*temp/18.01534/1.01325d+06
      tau=amnt*(a1*ph2o+a2*(pres-ph2o))
c     tau=amnt*(a1*ph2o)
      if(tau.gt.88)then
        twin=88.0
      else
        twin=tau
      end if
      return
      end
c *********************************************************************
c     o3 1140--1200 cm^{-1}
      subroutine cko3_85_ckm(ni,mlv,u,delw,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder,num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)
       
      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=15.60821*(1.0000+4.6204e-03*(t(m)-250.0)+
     *    8.7844e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
c *********************************************************************
c     ch4 1140--1200 cm^{-1}
      subroutine ckch4_85_ckm(ni,mlv,u,delw,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder,num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)
       
      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=3.44205*(1.0000+1.8357e-02*(t(m)-250.0)+
     *    1.4238e-04*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
c *********************************************************************
c     n2o 1140--1200 cm^{-1}
      subroutine ckn2o_85_ckm(ni,mlv,u,delw,f,p,t,tau)
      include './corrk.h'
      parameter (num_pressure = 19)
      parameter (max_coefs = 6)
      parameter (iorder = 3)
      real f(max_coefs),k(max_coefs),tau(max_coefs,MAXCKLEV)
      real u(MAXCKLEV), t(MAXCKLEV), p(MAXCKLEV)
      real coefk(max_coefs, iorder,num_pressure)
      real stp(num_pressure), fkg(max_coefs,MAXCKLEV)
       
      ni=1
      f(1)=1.00000
      do m=1,mlv
        k(1)=0.93312*69.15198736*(1.0000+7.9894e-04*(t(m)-250.0)+
     *    2.6735e-06*(t(m)-250.0)**2)
        tau(1,m)=k(1)*u(m)
      end do
      return
      end
      
      Subroutine emit_from_eprime_ckm(T3TOA,T4TOA,T5TOA,T85TOA,chi,sz,vz,iuse_downwelling,e_prime,
     +                            SurfEmit3, SurfEmit4,SurfEmit5,SurfEmit85,Tskin_new)
c
c     This subroutine calculates the channel 3 and channel 4 surface emittance in the daytime using
c	     the the nighttime emittance ratio (e_prime).  It assumes that CorrKsetup_ckm has already
c        been called.
c     This version includes formalae for calculating the surface emittance that take into
c        account the reflected down_welling radiation.
c     The downwelling terms can be turned off using iuse_downwelling = 0.  This causes the 
c        subroutine to revert to our old formulation for determining surface emittance.
c
c	INPUT:
c        T3TOA 			  = the channel 3 TOA clear temperature (satellite observation)
c        T4TOA 			  = the channel 4 TOA clear temperature (satellite observation)
c        chi              = surface clear anisotropic factor
c        sz 			  = solar zenith angle
c        vz 			  = viewing zenith angle
c        iuse_downwelling = Downwelling radiation option flag
c                         = 1 means thatt we include downwelling radiation in our calculations
c                         = 0 means thatt we calculate surfemit the old way (no downwelling)
c        e_prime          = the nighttime channel3/channel 4 surface emittance ratio
c                           (e_prime is calculated slightly differently for the downwelling &
c                            no-downwelling cases.  Make sure that you use the proper version)
c	OUTPUT:
c        SurfEmit3        = Channel 3 surface emittance
c        SurfEmit4        = Channel 4 surface emittance
c        Tskin_new        = Derived skin temperature 
c	LOCAL:
c        Rad4surf_down    = Channel 4 downwelling radiation
c        Rad4surf_up_total= Channel 4 total upwelling radiation at the surface
c        Rad4surf_up      = Channel 4 total upwelling radiation at the surface from emission
c                           Rad4surf_up = Rad4surf_up_total if SurfEmit4 = 1 or if 
c                           Rad4surf_down = 0.
c        T4surf           = T4TOA after molecular absorption is removed.  This contains contributions
c                           Tskin, SurfEmit4, and Rad4surf_down
c
c        Rad3surf_down    = Channel 3 downwelling radiation
c        Rad3surf_up      = Channel 3 total upwelling radiation at the surface 
c        Rad3skin_up      = Channel 3 radiation for a blackbody temperature of T4surf
c
      include './corrk.h'

      integer MaxChanCK
      parameter (MaxChanCK=4)

      real attenuate(MAXCKLEV,MaxChanCK),ratiorefl(MAXCKLEV),solarrad(MAXCKLEV),
     +  ratio_up(MAXCKLEV)  
      COMMON /ktaus/ attenuate,ratiorefl,TotalRefRat,solarrad,ratio_up

      logical iuse_downwelling
c
c   Call FindCloudLayer_ckm with T = 0. to get index of top layer of sounding
c		
      Ttop = 0.
      Call FindCloudLayer_ckm(Ttop,laytop)
c
c   Perform calculations with Channel 4 first.  The steps:
c    1)  Define downwelling radiation using correlated k's
c    2)  Derive T4surf by removing moecular absorption from T4TOA using correlated k's
c
      ichan = 2
      if(iuse_downwelling) then
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 2
        RadOrig_down = 0.
        Call UseCorrK_down_ckm(RadOrig_down,laylow,layhigh,iupdown,ichan,vz,Rad4surf_down)
      end if
      iwr = 0
      RadOrig = planck_ckm(ichan,T4TOA)
      laylow = 1
      layhigh = laytop
      iupdown = 2
      ichan = 2
      Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Rad4surf_up_total)
c **** Rad4surf_up = (Rad4surf_up_total - Rad4surf_down) / surfemit4 + Rad4surf_down
      Rad4surf_up = Rad4surf_up_total
      T4surf = plancki_ckm(ichan,Rad4surf_up)
c
c  Now for channel 5
c
      ichan = 3
      if(iuse_downwelling) then
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 3
        RadOrig_down = 0.
        Call UseCorrK_down_ckm(RadOrig_down,laylow,layhigh,iupdown,ichan,vz,Rad5surf_down)
      end if
      ichan = 3
      RadOrig = planck_ckm(ichan,T5TOA)
      laylow = 1
      layhigh = laytop
      iupdown = 2
      ichan = 3
      Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Rad5surf_up_total)

c
c  Now for 8.5 µm channel 
c
      if(iuse_downwelling) then
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 4
        RadOrig_down = 0.
        Call UseCorrK_down_ckm(RadOrig_down,laylow,layhigh,iupdown,ichan,vz,Rad85surf_down)
      end if
      ichan = 4
      RadOrig = planck_ckm(ichan,T85TOA)
      laylow = 1
      layhigh = laytop
      iupdown = 2
      ichan = 4
      Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Rad85surf_up_total)

c
c   Perform calculations with Channel 3.  The steps:
c    1)  Define downwelling radiation using correlated k's
c    2)  Derive T3surf by removing moecular absorption from T3TOA using correlated k's
c    3)  Derive surfemit3
c    4)  Derive Tskin_new = the skin temperature
c    5)  Use Tskin_new to derive Surfemit4
c	       	  
      ichan = 1
      if(iuse_downwelling) then
        laylow = 1
        layhigh = laytop
        iupdown = 2
        ichan = 1
        RadOrig_down = 0.
        Call UseCorrK_down_ckm(RadOrig_down,laylow,layhigh,iupdown,ichan,vz,Rad3surf_down)
      end if		  
      RadOrig = planck_ckm(ichan,T3TOA)
      laylow = 1
      layhigh = laytop
      iupdown = 2
      ichan = 1
      Call UseCorrK_ckm(RadOrig,laylow,layhigh,iupdown,ichan,Rad3surf_up)
      T3surf = plancki_ckm(ichan,Rad3surf_up)
      Rad3skin_up = planck_ckm(ichan,T4surf)
		  		  
      if(iuse_downwelling) then
c ****  surfemit3_new = 1.0 - (Rad3surf_up - e_prime*Rad3skin_up) / 
c ****+                        (chi*solarrad(1) + Rad3surf_down)
        surfemit3_new = 1.0 - (Rad3surf_up - e_prime*Rad3skin_up - 
     +                        (1. - e_prime)*Rad3surf_down) / (chi*solarrad(1))

        Rad3skin_new = (Rad3surf_up - (1.0 - surfemit3_new)*(chi*solarrad(1) + 
     +                  Rad3surf_down)) / surfemit3_new
        ichan = 1
        Tskin_new = plancki_ckm(ichan,Rad3skin_new)
        ichan = 2		  
        Rad4skin_new = planck_ckm(ichan,Tskin_new)
        surfemit4_new = (Rad4surf_up_total - Rad4surf_down) / 
     +                  (Rad4skin_new - Rad4surf_down)
        ichan = 3		  
        Rad5skin_new = planck_ckm(ichan,Tskin_new)
        surfemit5_new = (Rad5surf_up_total - Rad5surf_down) / 
     +                  (Rad5skin_new - Rad5surf_down)
        ichan = 4		  
        Rad85skin_new = planck_ckm(ichan,Tskin_new)
        surfemit85_new = (Rad85surf_up_total - Rad85surf_down) / 
     +                  (Rad85skin_new - Rad85surf_down)
      else
        surfemit3_new = 1.0 - (Rad3surf_up - e_prime*Rad3skin_up) / 
     +                        (chi*solarrad(1))
        Rad3skin_new = (Rad3surf_up - (1.0 - surfemit3_new)*(chi*solarrad(1))) / 
     +                  surfemit3_new
        ichan = 1
        Tskin_new = plancki_ckm(ichan,Rad3skin_new)
        ichan = 2		  
        Rad4skin_new = planck_ckm(ichan,Tskin_new)
        surfemit4_new = (Rad4surf_up_total) / (Rad4skin_new)
        ichan = 3		  
        Rad5skin_new = planck_ckm(ichan,Tskin_new)
        surfemit5_new = (Rad5surf_up_total) / (Rad5skin_new)
        ichan = 4		  
        Rad85skin_new = planck_ckm(ichan,Tskin_new)
        surfemit85_new = (Rad85surf_up_total) / (Rad85skin_new)
      end if

      SurfEmit3 = surfemit3_new
      SurfEmit4 = surfemit4_new
      SurfEmit5 = surfemit5_new
      SurfEmit85 = surfemit85_new
      
      return
      end
	  
