program sigt
c
c TEST VERSION (21JUN96)
c test of sig to pressure interp code using specific data extracted
c  with wgrib from the sanl files.  The code expects an f77 ieee array of
c  real values for 28 sigma levels followed by the surface pressure.
c
      parameter (id=192,jd=94,lv=28,lp=3)
      dimension dsfcp(id,jd),xsig(id,jd,lv),xprs(id,jd,lp),dprs(lp)
      dimension asfcp(id,jd),aprs(lp),dsig(lv),asig(lv+1)
      data dprs/850.,700.,10./
      data dsig/.9950,.9821,.9644,.9425,.9159,.8838,.8458,.8014,.7508,
     2          .6943,.6329,.5681,.5017,.4357,.3720,.3125,.2582,.2101,
     3          .1682,.1326,.1028,.0782,.0580,.0418,.0288,.0183,.0101,
     4          .0027/
      call sigint(dsig,asig,lv,dprs,aprs,lp)
      do k=1,28
	read(11)((xsig(i,j,k),i=1,id),j=1,jd)
	write(*,1001)k,xsig(96,47,k)
 1001   format(' read ',i5,f12.2)
      enddo
      read(11)dsfcp
	write(*,1001)99,dsfcp(96,47)
      write(*,*)' Call sigsfp '
      call sigsfp(dsfcp,asfcp,id,jd)
c
      write(*,*)' Call sig2prs'
      call sig2prs(xsig,xprs,id,jd,lv,lp,asfcp,asig,aprs,lp)
      open(21,access='direct',recl=id*jd*4)
      do lw=1,lp
	write(21,rec=lw)((xprs(i,j,lw),i=1,id),j=1,jd)
      enddo
      close(21)
      end
c
      subroutine sig2prs(xsig,xprs,ip,jp,ls,lp,asfcp,asig,aprs,nprs)
c
c Interpolate data fields from sigma levels to pressure levels.
c  (initialization calls may be necessary)
c
c on input
c       xsig - full stack of sigma layer data fields
c       ip,jp,ls - dimensions of xsig array
c       asfcp    - log of surface pressure at the same time as the input
c                    xsig data (can be computed from the surface pressure
c                    field itself (sfcp) by calling sigsfp.
c       asig     - log of the sigma level values which correspond to the
c                    levels in the data array xsig, can be computed from
c                    the  sigma level values by calling sigint.
c       aprs     - log of the pressure values for the levels desired in
c                    the output grids, can be computer from pressure
c                    values by calling singint.
c
c       nprs     - number of levels in the desired output array, must be
c                     equal to or less than lp
c       lp       - the k dimension of the pressure level output data array
c                   (i and j dimensions are ip and jp)
c on output
c       xprs     - data values at the pressure levels specified by prs,
c                  The interpolation is simple linear in log P.  Any points
c                  below the lowest sigma level are set to the value of the
c                  lowest sigma level value.  Any points above the highest
c                  sigma level are set to the highest sigma level value.
c
      dimension xprs(ip,jp,lp),dprs(nprs),aprs(nprs)
      dimension sfcp(ip,jp),xsig(ip,jp,ls),dsig(ls),asig(ls)
      dimension asfcp(ip,jp),asig(ls+1)
      do j=1,jp
	do i=1,ip
	  lb=0
	  do n=1,nprs
	    axp=aprs(n)-asfcp(i,j)
c     write(*,*)n,i,j,axp,aprs(n),asfcp(i,j),lb,asig(lb+1)
c     write(*,*)n,lb,exp(axp),exp(asig(lb)),i,j,n,xprs(i,j,n)
   10       if(axp.lt.asig(lb+1)) then
	      lb=lb+1
	      go to 10
	    endif
	    if(lb.le.0) then
	      xprs(i,j,n)=xsig(i,j,1)
	    elseif(lb.ge.ls) then
	      xprs(i,j,n)=xsig(i,j,ls)
	    else
	      rr = (axp-asig(lb))/(asig(lb+1)-asig(lb))
	      xprs(i,j,n)=rr*xsig(i,j,lb+1)+(1.-rr)*xsig(i,j,lb)
	    endif
	  enddo
	enddo
      enddo
      return
c
      entry sigsfp(sfcp,asfcp,ip,jp)
c compute log of surface pressure (asfcp) from the surface pressure
c  field dimensioned ip x jp. Units of sfcp are Pascals as is
c  normally the case in the NCEP Reanalysis data.
      do j=1,jp
       do i=1,ip
	 asfcp(i,j)=log(sfcp(i,j))
       enddo
      enddo
      return
c
      entry sigint(dsig,asig,ls,dprs,aprs,nprs)
c compute log of sigma levels (asig) from sigma values (dsig),
c   for ls levels.
c compute log of pressure surfaces (aprs) from the array of
c   desired pressure surfaces (units mb (hectopascals)) for
c   nprs pressure levels.
      do i=1,ls
	asig(i)=log(dsig(i))
      enddo
      asig(ls+1)=-999999.
      do i=1,nprs
	aprs(i)=log(dprs(i)*100.)
      enddo
      return
      end