subroutine short(rh,ntau,xu,wu,den,col,prs,hast)

c	***************************************************************
c	   This program calculates heat rates according to MecKay
c	 (Titan paper).      
c	   Modification to calculate the heating rate usinng c-k algrium
c	instead of taking subtractions from the equavalent width.
c						Y. Wang 5/1993
c	***************************************************************

        implicit double precision (a-h,o-z)
        include "parameters.dat"
	character fname*60

	real timeused,time(2)
	real dtime 
	data ch4m /2.67d-23/
	data rs /9.539d0/
c	-----------------------------------------------------------------
c			Input arrays
c	-----------------------------------------------------------------
	dimension xu(3),wu(3),den(nd1),col(nd1),prs(nd1),dcol(nd1)
c	-----------------------------------------------------------------
c			Output Arrays
c	-----------------------------------------------------------------
     &	,hast(nd1), hack(nd1)
c	-----------------------------------------------------------------
c			Internal Arrays
c	-----------------------------------------------------------------
     &	,wmid(50),weight(4,50),zk(4,50),h(nd1),op(4,nd1,50),zkn(4,nd1,50)
     &  , solfx(50),factor(50),trans(nd1,50),qw(nd1),dwmid(50)
     &  , w1(50),w2(50)

c 	write(*,*)xu(1),wu(1) 
	write(*,*)'in short,ntau=',ntau 
	do n = 1,ntau
c 	   write(*,*)'in short',den(n),col(n)
	end do
	timeused = dtime(time)

	call stread(wmid,weight,zk,solfx,factor,ndat)

c	write(*,*)'stread get called'

	do n = 1,ntau
	   h(n) = zero 
	   hast(n) = zero
	end do

	do n = 2,ntau
	    dcol(n) = col(n) - col(n-1)
c	    write(*,*)dcol(n)
	end do
	
c	write(*,*)prs(1),factor(15),zk(1,20),ch4m
	do n = 1,ntau
c	    write(*,*)col(n),prs(n),den(n)
	end do

	do i = 1,ndat
	   do n = 1,ntau
	      do j = 1,4
                 zkn(j,n,i) = (prs(n)/1.d6)**factor(i)*zk(j,i)
                 op(j,n,i) = (prs(n)/1.d6)**factor(i)*zk(j,i)*ch4m
	      end do
	   end do
	end do

c	-------------------------------------------------------------
c	       Changing the wavelengthes from um to nm
c	-------------------------------------------------------------

	do i = 1,ndat
	    wmid(i) = 1.d3*wmid(i)
	end do
c	-------------------------------------------------------------

c	write(*,*)op(4,21,4),factor(4),zk(4,4)
	do nu = 1,3
c	  write(*,*)'in short',wu(nu),xu(nu)
	end do

c	do i = 1,ndat
c          solfx(i)=four*solfx(i)*(rs**2/rh**2)
c          do n = 1,ntau
c	      qw(n) = zero 
c	      do j = 1,4
c		 trans(n,i) = zero
c	         do ni = 1,3 
c     		    dtau = dble(-zkn(j,n,i)*ch4m*col(n)/xu(ni))
c		    trans(n,i) = trans(n,i) + weight(j,i)*dexp(dtau)
c	            qw(n) = qw(n) + wu(ni)*zkn(j,n,i)*ch4m
c     &                            *weight(j,i)*dexp(dtau)
c		 end do
c	      end do 
c	          
c	      hast(n) = hast(n) + half*solfx(i)*qw(n)*den(n)
c	  end do 
c	end do

c	----------------------------------------------------------
c	     The numbers for the solar flux listed in McKay's 
c	Table V (short.dat) were divided by four to represent the
c	globally averaged conditions. See pp48 in McKay et al. 
c	(1989) Titan paper.
c	---------------------------------------------------------
c	write(*,103)
c	read(*,*)fname
c	open(39,file=fname)

	w1(1) = 3.d2
	w2(1) = 3.5d2
	wmid(ndat+1) = 1.742d3
	do l = 2,ndat
	   w1(l) = half*(wmid(l)+wmid(l-1))
	   w2(l) = half*(wmid(l)+wmid(l+1))
	end do

	do l = 2,ndat
	    dwmid(l) = wmid(l+1)-wmid(l-1)
	end do
	dwmid(1) = 5.d1

	do l = 1,ndat
           solfx(l)=four*solfx(l)*(rs**2/rh**2)
	   write(*,*)'in short,solfx=',solfx(l)
	   do nn = 1,ntau 
 	       sumnq = zero 
	       qw(nn) = zero 
	       do j = 1,4
	         sumnu = zero
	         do nu = 1,3
	           s1 = -op(j,1,l) 
	           sum = -op(j,1,l)*col(1)
	           do m = 2,nn
	             s2 = -op(j,m,l) 
	             sum = sum +half*(s2+s1)*dcol(m)
	             s1 = s2
	           end do
	           sumnu = sumnu + dexp(sum/xu(nu))*wu(nu)
                 end do 
 	         sumnq = sumnq+weight(j,l)*dexp(sum) 
	         qw(nn) = qw(nn) + op(j,nn,l)*weight(j,l)*sumnu
	       end do  
	       if(sumnq .ge. one) then
	           trans(nn,l) = one 
	       else
	           trans(nn,l)  =  sumnq
	       end if 

c	       hack(nn) = half*solfx(l)*qw(nn)
	       hack(nn) = qw(nn)
	   end do

	   do n = 1, ntau 
c	       hast(n) = hast(n) + hack(n)*den(n)
	       hast(n) = hast(n) + hack(n)
	   end do 
c	   if(l.eq.11 .or. l.eq.12 .or. l.eq.13 .or. l.eq.14) then
c	      write(39,104)wmid(l),w1(l),w2(l),hast(1)
c	      do n = 2,ntau
c	          write(39,105)hast(n)
c	      end do
c	   end if
	end do

c	write(*,101)
c	read(*,*)fname
c	open(29,file=fname)
c	do l = 1,ndat
c	    write(29,102) wmid(l),trans(1,l)
c	end do
c	close(29)  
	
	timeused = dtime(time)
	write(40,901)timeused

 101	format(' Enter filename for transmissions ->',$)
 102	format(1p2e12.2)
 103	format(' Enter filename for heating rates vs. wavelength >',$)
 104	format(1p4e12.2)
 105	format(36x,1pe12.2)
 
 901	format(//,'time used in  short',/' took',f10.3,'s'/)

         return
        end

c***************************************************************

	subroutine stread(wmid,weight,zk,solfx,factor,ndat)

        implicit double precision (a-h,o-z)
     	dimension  wmid(50),weight(4,50),zk(4,50),solfx(50),factor(50)

	open(19,file = 'short.dat',status='old')

	
	ndat = 20
	do i = 1,ndat 
           read(19,*)wmid(i),(weight(j,i),j=1,4)
           read(19,*)(zk(j,i),j=1,4)
	end do
	do i = 1,ndat
	   read(19,*)solfx(i),factor(i)
	end do
	close(19) 
c 	write(*,*)'in stread'
c  	do i = 1,ndat
c  	   write(*,101)wmid(i),(weight(j,i),j=1,4)
c  	   write(*,102)(zk(j,i),j=1,4)
c 	   write(*,103)solfx(i) 
c  	end do
  101	format(1x,5f7.3)
  102	format(8x,1p4e11.3)
  103	format(f8.2) 

	return
	end