c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/getgascoe.f,v 2.22 2000/08/30 16:38:09 flittner Exp $ subroutine getgascoe(wavelen) c c*********************************************************************** c cccc c subroutine getgascoe c c by: dave flittner c c purpose- c c Load the other (ngas-1) absorbing gas absorption coefficients. c The data base for the absorption coefficient of each gas is c read from file fgasalfa(igas). The is value computed either by: c c 1. Averaging over a slit with fwhm=bw (use 11 A for toms. c The slit type is set by itype=0 for rectangular or itype=1 c for triangular.) c or c 2. Linear interpolation to the desired wavelength. c c c input/output variables - c c variables type i/o description c --------- ---- --- ----------- c formal parameters c wavelen r*8 i wavelength in Angstroms c part of common 'cgas2' c alphagas r*8 o gaseous absorption coefficient c Units=[cm*cm] c c external function slitavg is called for averaging over a slit. c c last mod. March 13, 1998...def c purpose: Do not extrapolate to wavelength outside data set. Also, c do not bandaverage. cccc c*********************************************************************** implicit none include 'parameter.inc' include 'consts.cmn' c passed real*8 wavelen include 'cgas2.cmn' c c local c integer np,mw,nw,itype,mode,iunt,igas parameter (mw=20388) real*8 w(mw),coe(mw),bw,slitavg external slitavg logical sort,monoton,lbandavg c c use bandavgeraged values c lbandavg=.FALSE. c c now get the absorption coef. for this wavelength c call get_lun(iunt) do igas=1,ngas-1 open(iunt,file=fgasalfa(igas),status='old') nw=1 30 read(iunt,*,end=40)w(nw),coe(nw) nw=nw+1 goto 30 40 close(iunt) nw=nw-1 if(lbandavg)then itype=1 bw=11.0d0 alphagas(igas)=slitavg(w,coe,nw,wavelen,bw,itype) else c c now interpolate to the desired wavelength c c do not extrapolate if((wavelen.lt.w(1)).or.(wavelen.gt.w(nw)))then alphagas(igas)=0.0d0 c write(6,*)'getgascoe. wavelength not in ', c &'absorption coef. data set for gas ',igas+1 else sort=.false. c use linear interpolation mode=4 np=1 monoton=.true. call smooth(w,coe,nw,nw,sort,wavelen,alphagas(igas), &np,np,mode,monoton) endif endif c the following is if the input coes are multiplied by 1.0e20 c alphagas(igas)=alphagas(igas)*1.0d-20 c write(6,1000)igas+1,alphagas(igas),alphagas(igas)*Nlosch enddo 1000 format(' getgas. abs. coe. for gas ',i2,' = ',1pe12.4, &' [cm^2] or ',1pe12.4,' [1/atm-cm]') return end c*********************************************************************** subroutine smooth(x,y,num,no,sort,ix,iy,n,np,mode, &monoton) c c routine to interpolate use local derivative, resultant c line is 'smooth', much like a spline fit. C C Reference: Hiroshi Akima, "A new method of interpolation and C smooth curve fitting based on local procedures". C Journal of the Association for Computing Machinery, C Vol. 17, No. 4, October 1970, pp. 589-602. c c Feb. 1, 1990 c c np=dimensions of interpolated arrays c no=dimensions of data arrays c x,y=data arrays c ix,iy=interpolated arrays c num=number of data pairs c n=number of interpolated pairs c implicit none logical sort,linear,up,linend,noover,nounder,monoton,prt integer*4 num,no,n,np,i,j,k,ibegin(2),iend(2),ij,KHI,KLO integer*4 itwosame,itmp,mode real*8 ix(np),iy(np),x(no),y(no),m(4),t(2),p(4) real*8 xtmp,ytmp,dx,dy,dix integer modeo c c Dec. 10, 1991 c If mode=0, perform cubic interpolation in x-y space. c If mode=1, perform cubic interpolation in ln(x)-y space. c If mode=2, perform cubic interpolation in x-ln(y) space. c If mode=3, perform cubic interpolation in ln(x)-ln(y) space. c If mode=4, perform linear interpolation instead of cubic in x-y space. c If mode=5, perform linear interpolation instead of cubic in ln(x)-ln(y) space. c If mode=6, perform linear interpolation instead of cubic in x-ln(y) space. c c if linend = TRUE, points beyond the ends will be computed using linear c function instead of quadratic. prt=.FALSE. linend=.FALSE. itwosame=0 c noover=monoton nounder=monoton modeo=mode c check to see if two or more values of x(I) are the same c also sort values in ascending order. if(sort)then do i=1,num xtmp=x(i) ytmp=y(i) itmp=i do j=1+i,num if(x(j).lt.xtmp)then itmp=j xtmp=x(j) ytmp=y(j) end if enddo if(itmp.gt.i) then x(itmp)=x(i) x(i)=xtmp y(itmp)=y(i) y(i)=ytmp end if if(i.ne.1)then if(x(i).eq.x(i-1))then print*,'2 abcissa points are the same',x(i) itwosame=1 end if end if enddo if (itwosame.eq.1)goto 1000 end if c check to is if x values are in ascending order. if(x(1).lt.x(num))then up=.TRUE. else up=.FALSE. endif c loop for the number of interpolation points if(up)then do 70 i=1,n mode=modeo c find data pairs that bracket the interpolation point if(up)then KLO=1 KHI=NUM else KLO=NUM KHI=1 endif do while(abs(KHI-KLO).GT.1) K=(KHI+KLO)/2 IF(X(K).EQ.IX(i))then IY(i)=Y(K) GOTO 60 ELSE IF(X(K).GT.IX(i))THEN KHI=K ELSE KLO=K ENDIF enddo if(mode.eq.4) then iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 60 endif do 35 ij=1,2 k=KLO+ij-1 if(((k).ge.3).and.((k).le.(num-2)))then ibegin(ij)=1 iend(ij)=4 else if(((k).eq.1))then ibegin(ij)=3 iend(ij)=4 else if((k).eq.2)then ibegin(ij)=2 iend(ij)=4 else if((k).eq.num-1)then ibegin(ij)=1 iend(ij)=3 else if((k).eq.num)then ibegin(ij)=1 iend(ij)=2 end if 35 continue if(mode.eq.0)then linear=.TRUE. goto 44 else if(mode.eq.1)then if(ix(i).le.0.0d0)then mode=0 linear=.TRUE. goto 44 end if do 40 j=KLO-3+ibegin(1),KLO+1-2+iend(2) if(x(j).le.0.0d0) then mode=0 linear=.TRUE. goto 44 end if 40 continue else if(mode.eq.2)then do 41 j=KLO-3+ibegin(1),KLO+1-2+iend(2) if(y(j).le.0.0d0) then mode=0 linear=.TRUE. goto 44 end if 41 continue else if(mode.eq.3)then c check to see if any abssisa is less than or equal to zero do 42 j=KLO-3+ibegin(1),KLO+1-2+iend(2) if(x(j).le.0.0d0) then mode=2 end if 42 continue do 43 j=KLO-3+ibegin(1),KLO+1-2+iend(2) if(y(j).le.0.0d0) then if(mode.eq.2) then linear=.TRUE. mode=0 else mode=1 endif goto 44 end if 43 continue else if(mode.eq.5)then if(x(KLO+1).le.0.0d0)then mode=4 else if(x(KLO).le.0.0d0)then mode=4 else if(y(KLO+1).le.0.0d0)then mode=4 else if(y(KLO).le.0.0d0)then mode=4 else dx=log(x(KLO+1))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=log(y(KLO+1))-log(y(KLO)) p(1)=log(y(KLO)) if(dx.ne.0.0d0)then p(2)=dy/dx else p(2)=0.0d0 endif iy(i)=p(1)+p(2)*(dix) iy(i)=exp(iy(i)) goto 60 endif iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 60 else if(mode.eq.6)then if(y(KLO+1).le.0.0d0)then mode=4 else if(y(KLO).le.0.0d0)then mode=4 else dx=(x(KLO+1))-(x(KLO)) dix=(ix(i))-(x(KLO)) dy=log(y(KLO+1))-log(y(KLO)) p(1)=log(y(KLO)) if(dx.ne.0.0d0)then p(2)=dy/dx else p(2)=0.0d0 endif iy(i)=p(1)+p(2)*(dix) iy(i)=exp(iy(i)) goto 60 endif iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 60 endif linear=.FALSE. 44 do 50 ij=1,2 k=KLO+ij-1 if(linear.or.mode.eq.0)then do 45 j=ibegin(ij),iend(ij) m(j)=(y(k+j-2)-y(k-3+j))/(x(k+j-2)-x(k+j-3)) 45 continue else if(mode.eq.1)then do 46 j=ibegin(ij),iend(ij) m(j)=((y(k+j-2))-(y(k-3+j)))/ & (log(x(k+j-2))-log(x(k+j-3))) 46 continue else if(mode.eq.2)then do 47 j=ibegin(ij),iend(ij) m(j)=(log(y(k+j-2))-log(y(k-3+j)))/ & ((x(k+j-2))-(x(k+j-3))) 47 continue else if(mode.eq.3)then do 48 j=ibegin(ij),iend(ij) m(j)=(log(y(k+j-2))-log(y(k-3+j)))/ & (log(x(k+j-2))-log(x(k+j-3))) 48 continue endif if(ibegin(ij).eq.3)then if(linend)then m(2)=m(3) m(1)=m(3) else m(2)=2.0d0*m(3)-m(4) m(1)=3.0d0*m(3)-2.0d0*m(4) endif else if(ibegin(ij).eq.2)then c m(1)=m(2)+m(3)-m(4) if(linend)then m(1)=m(2) else m(1)=2.0d0*m(2)-m(3) endif else if(iend(ij).eq.2)then if(linend)then m(3)=m(2) m(4)=m(2) else m(3)=2.0d0*m(2)-m(1) m(4)=3.0d0*m(2)-2.0d0*m(1) endif else if(iend(ij).eq.3)then if(linend)then m(4)=m(3) else m(4)=2.0d0*m(3)-m(2) endif end if c if points are equall if((m(1).eq.m(2)).and.(m(3).eq.m(4))) then t(ij)=0.5d0*(m(2)+m(3)) else t(ij)=(abs(m(4)-m(3))*m(2)+abs(m(2)-m(1))*m(3)) * /(abs(m(4)-m(3))+abs(m(2)-m(1))) end if 50 continue if(linear.or.mode.eq.0)then dx=x(KLO+1)-x(KLO) dix=ix(i)-x(KLO) dy=y(KLO+1)-y(KLO) p(1)=y(KLO) else if(mode.eq.1)then dx=log(x(KHI))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=y(KHI)-y(KLO) p(1)=y(KLO) else if(mode.eq.2)then dx=(x(KLO+1))-(x(KLO)) dix=(ix(i))-(x(KLO)) dy=log(y(KLO+1))-log(y(KLO)) p(1)=log(y(KLO)) else if(mode.eq.3)then dx=log(x(KLO+1))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=log(y(KLO+1))-log(y(KLO)) p(1)=log(y(KLO)) end if p(2)=t(1) p(3)=(3.0d0*dy/dx-2.0d0*t(1)-t(2))/dx p(4)=(t(1)+t(2)-2.0d0*dy/dx)/(dx*dx) iy(i)=p(1)+p(2)*(dix)+p(3)*(dix*dix)+p(4)*(dix*dix*dix) if(mode.ge.2)then iy(i)=exp(iy(i)) end if 59 if(nounder)then if((iy(i).lt.y(KLO)).and.(iy(i).lt.y(KHI)))then if(prt)then write(6,*)'under shoots between',KLO,KHI write(6,*)'t1',t(1),'t2',t(2) write(6,*)'slope=',p(3)*2.0+6.0*p(4)*ix(i) c will use linear write(6,*)'old',ix(i),iy(i) endif iy(i)=p(1)+dy*dix/dx if(mode.ge.2)then iy(i)=exp(iy(i)) end if if(prt)then write(6,*)x(KLO),y(KLO) write(6,*)x(KHI),y(KHI) write(6,*)'new',ix(i),iy(i) endif endif endif if(noover)then if((iy(i).gt.y(KLO)).and.(iy(i).gt.y(KHI)))then if(prt)then write(6,*)'going up' write(6,*)'over shoots between',KLO,KHI write(6,*)'t1',t(1),'t2',t(2) write(6,*)'slope=',p(3)*2.0+6.0*p(4)*ix(i) write(6,*)'iy=',iy(i) endif c will use linear if(dx.eq.0)then iy(i)=p(1) else iy(i)=p(1)+dy*dix/dx endif if(mode.ge.2)then iy(i)=exp(iy(i)) end if if(prt)then write(6,*)'p(1)',p(1),'dy=',dy write(6,*)'dix=',dix,dx write(6,*)'iy=',iy(i) write(6,*)'ix=',ix(i),'i=',i,'mode=',mode write(6,*)'ibegin=',ibegin(1),ibegin(2) write(6,*)'iend=',iend(1),iend(2) do k=1,num write(6,*)'x=',x(k),' y=',y(k) enddo endif endif endif 60 continue 70 continue else do 900 i=1,n mode=modeo c find data pairs that bracket the interpolation point KLO=NUM KHI=1 75 IF (KLO-KHI.GT.1) THEN K=(KHI+KLO)/2 IF(X(K).EQ.IX(i))then IY(i)=Y(K) GOTO 800 ELSE IF(X(K).GT.IX(i))THEN KHI=K ELSE KLO=K ENDIF GOTO 75 ENDIF if(mode.eq.4) then iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 900 endif do 80 ij=1,2 k=KLO-ij+1 if(((k).ge.3).and.((k).le.(num-2)))then ibegin(ij)=1 iend(ij)=4 else if(((k).eq.1))then ibegin(ij)=1 iend(ij)=2 else if((k).eq.2)then ibegin(ij)=1 iend(ij)=3 else if((k).eq.num-1)then ibegin(ij)=2 iend(ij)=4 else if((k).eq.num)then ibegin(ij)=3 iend(ij)=4 end if 80 continue if(mode.eq.0)then linear=.TRUE. goto 84 else if(mode.eq.1)then if(ix(i).le.0.0d0)then linear=.TRUE. mode=0 goto 84 else do 81 j=KLO-1+2-iend(2),KLO+3-ibegin(1) if(x(j).le.0.0d0) then linear=.TRUE. mode=0 goto 84 end if 81 continue end if else if(mode.eq.2)then do 82 j=KLO-1+2-iend(2),KLO+3-ibegin(1) if(y(j).le.0.0d0) then linear=.TRUE. mode=0 goto 84 end if 82 continue else if(mode.eq.3)then if(ix(i).le.0.0d0)then linear=.TRUE. mode=0 goto 84 endif do j=KLO-1+2-iend(2),KLO+3-ibegin(1) if(x(j).le.0.0d0) then mode=2 end if enddo do 83 j=KLO-1+2-iend(2),KLO+3-ibegin(1) if(y(j).le.0.0d0) then if(mode.eq.2) then linear=.TRUE. mode=0 else mode=1 endif goto 84 end if 83 continue else if(mode.eq.5)then if(x(KLO-1).le.0.0d0)then mode=4 else if(x(KLO).le.0.0d0)then mode=4 else if(y(KLO-1).le.0.0d0)then mode=4 else if(y(KLO).le.0.0d0)then mode=4 else dx=log(x(KLO-1))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=log(y(KLO-1))-log(y(KLO)) p(1)=log(y(KLO)) if(dx.ne.0.0d0)then p(2)=dy/dx else p(2)=0.0d0 endif iy(i)=p(1)+p(2)*(dix) iy(i)=exp(iy(i)) goto 800 endif iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 800 else if(mode.eq.6)then if(y(KLO-1).le.0.0d0)then mode=4 else if(y(KLO).le.0.0d0)then mode=4 else dx=(x(KLO-1))-(x(KLO)) dix=(ix(i))-(x(KLO)) dy=log(y(KLO-1))-log(y(KLO)) p(1)=log(y(KLO)) if(dx.ne.0.0d0)then p(2)=dy/dx else p(2)=0.0d0 endif iy(i)=p(1)+p(2)*(dix) iy(i)=exp(iy(i)) goto 800 endif iy(i)=y(KLO)+(y(KHI)-y(KLO))* & (ix(i)-x(KLO))/(x(KHI)-x(KLO)) goto 800 endif linear=.FALSE. 84 do 90 ij=1,2 k=KLO-ij+1 if(linear.or.mode.eq.0)then do 85 j=ibegin(ij),iend(ij) m(j)=(y(k+2-j)-y(k+3-j))/ & (x(k+2-j)-x(k+3-j)) 85 continue else if(mode.eq.1)then do 86 j=ibegin(ij),iend(ij) m(j)=((y(k+2-j))-(y(k+3-j)))/ & (log(x(k+2-j))-log(x(k+3-j))) 86 continue else if(mode.eq.2)then do 87 j=ibegin(ij),iend(ij) m(j)=(log(y(k+2-j))-log(y(k+3-j)))/ & ((x(k+2-j))-(x(k+3-j))) 87 continue else if(mode.eq.3)then do 88 j=ibegin(ij),iend(ij) m(j)=(log(y(k+2-j))-log(y(k+3-j)))/ & (log(x(k+2-j))-log(x(k+3-j))) 88 continue endif if(ibegin(ij).eq.3)then m(2)=2*m(3)-m(4) m(1)=m(2)+m(3)-m(4) else if(ibegin(ij).eq.2)then m(1)=m(2)+m(3)-m(4) else if(iend(ij).eq.2)then m(3)=2*m(2)-m(1) m(4)=m(3)+m(2)-m(1) else if(iend(ij).eq.3)then m(4)=m(3)+m(2)-m(1) end if if((m(1).eq.m(2)).and.(m(3).eq.m(4))) then t(ij)=0.5*(m(2)+m(3)) else t(ij)=(abs(m(4)-m(3))*m(2)+abs(m(2)-m(1))*m(3)) * /(abs(m(4)-m(3))+abs(m(2)-m(1))) end if 90 continue if(linear.or.mode.eq.0)then dx=x(KLO-1)-x(KLO) dix=ix(i)-x(KLO) dy=y(KLO-1)-y(KLO) p(1)=y(KLO) else if(mode.eq.1)then dx=log(x(KLO-1))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=y(KLO-1)-y(KLO) p(1)=y(KLO) else if(mode.eq.2)then dx=(x(KLO-1))-(x(KLO)) dix=(ix(i))-(x(KLO)) dy=log(y(KLO-1))-log(y(KLO)) p(1)=log(y(KLO)) else if(mode.eq.3)then dx=log(x(KLO-1))-log(x(KLO)) dix=log(ix(i))-log(x(KLO)) dy=log(y(KLO-1))-log(y(KLO)) p(1)=log(y(KLO)) end if p(2)=t(1) p(3)=(3.0d0*dy/dx-2.0d0*t(1)-t(2))/dx p(4)=(t(1)+t(2)-2.0d0*(dy/dx))/(dx*dx) iy(i)=p(1)+p(2)*(dix)+p(3)*(dix*dix)+p(4)*(dix*dix*dix) if(mode.ge.2)then iy(i)=exp(iy(i)) end if if(nounder)then if((iy(i).lt.y(KLO)).and.(iy(i).lt.y(KHI)))then if(prt)then write(6,*)'under shoots between',KLO,KHI write(6,*)'t1',t(1),'t2',t(2) write(6,*)'slope=',p(3)*2.0+6.0*p(4)*ix(i) write(6,*)'old',ix(i),iy(i) endif c will use linear iy(i)=p(1)+dy*dix/dx if(mode.ge.2)then iy(i)=exp(iy(i)) end if if(prt)then write(6,*)x(KLO),y(KLO),KLO write(6,*)x(KHI),y(KHI),KHI write(6,*)ix(i),iy(i) write(6,*)'p(1)',p(1),dy write(6,*)'dix',dix,dx endif endif endif if(noover)then if((iy(i).gt.y(KLO)).and.(iy(i).gt.y(KHI)))then if(prt)then write(6,*)'over shoots between',KLO,KHI write(6,*)'t1',t(1),'t2',t(2) write(6,*)'slope=',p(3)*2.0+6.0*p(4)*ix(i) endif c will use linear if(dx.eq.0)then iy(i)=p(1) else iy(i)=p(1)+dy*dix/dx endif if(mode.ge.2)then iy(i)=exp(iy(i)) end if if(prt)then write(6,*)x(KLO),y(KLO) write(6,*)x(KHI),y(KHI),iy(i) endif endif endif 800 continue 900 continue endif 1000 return end c***********************************************************************