subroutine relayr(delx,xtop,tmp10,tmptop,x101,tmp101,
     1  pnot,xpnot,h,ps)
c ----------------------------------------------------------------
c      subroutine relayr
c
c
c     version 2.0   june 1984
c     algorithm designed by dr. p.k. bhartia ; coded by david lee ,sasc
c
c      purpose -
c
c          find the cumulative ozone amount and ozone-weighted average
c          temperatures at 101 levels from input values for 10 umkehr
c          layers using a spline fit.
c
c
c      method / procedures
c
c          1) define input and output pressure levels
c          2) compute accumulated ozone amounts at input levels
c          3) compute ozone weighted average temperatures above input
c             pressure levels
c          4) perform a spline fit to temperatures and accumulated
c             ozones at input levels
c          5) evaluate spline fit at the lowest 61 pressure levels
c          6) find slope for the accumulated ozone at the 41st. level
c          7) evaluaate accumulated ozone amounts for the top 40 levels
c             from the slope computed from the 41st level
c          8) compute pressures for altitudes at 1 km intervals
c
c
c     calling sequence - call relayr(delx,xtop,tmp10,tmptop,x101,tmp101,
c                        h,ps)
c
c
c     variable     type    i/o   description
c     --------     ----    ---   -----------
c
c     delx(10)      r*8     i    layer ozone amount (d.u) for
c                                10 umkehr layers. (layers 0-9, layer
c                                0 being the bottom half of the atmos-
c                                phere). indexing of the layers are
c                                described in the table below
c
c     xtop          r*8     i    cumulative ozone amount(m-atm-cm)
c                                above umkehr layer 9
c
c     tmp10(10)     r*8     i    average temperatures (degrees kelvin)
c                                for 10 umkehr layers
c
c     tmptop        r*8     i    average temperature for the
c                                atmosphere above umkehr layer 9
c
c     p11(11)       r*8          pressure levels (atm.) defining the 10
c                                umkehr layers. indexing of the pressure
c                                levels are described in the table
c                                below
c
c     x11(11)       r*8          cumulative ozone amount (m-atm-cm)
c                                at pressure levels described above
c                                for p11
c
c     pnot          r*8    i     pressure of reflecting surface
c
c     xpnot         r*8    o     cumulative ozone at pnot
c
c     p101(101)     r*8          pressure (atm.) at 101 atmospheric
c                                levels that define the output layer
c                                ozone amounts and temperatures
c
c     x101(101)     r*8    o     output cumulative ozone amounts
c                                (atm-cm) at levels corresponding
c                                to p101.
c
c
c     tmp101(101)   r*8    o     output ozone weighted temperatures
c                                (deg-kelvin) for atmospheres above
c                                each of the 101 pressure levels.
c
c     z11(11)       r*8          altitude (km) corresponding to p11
c                                defining the umkehr layers
c
c     h(82)         r*8    o     altitude array at 1 km interval from
c                                ground to 82 km
c
c     ps(82)        r*8    o     pressures corresponding to the height
c                                at the h array
c
c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c
c     the table below describes the order of indices for variables
c     in this subroutine
c
c     pressure       umkehr        indices         indices      index
c      (atm.)        layer #     for p11,x11       for delx      for
c                                  & t11ave        & tmp10       z11
c
c      1/1024 -------------------------1--------------------------11
c                       9                             1
c      1/512---------------------------2--------------------------10
c                       8                             2
c      1/256---------------------------3---------------------------9
c                       7                             3
c      1/128---------------------------4---------------------------8
c                       6                             4
c      1/64 ---------------------------5---------------------------7
c                       5                             5
c      1/32 ---------------------------6---------------------------6
c                       4                             6
c      1/16 ---------------------------7---------------------------5
c                       3                             7
c      1/8  ---------------------------8---------------------------4
c                       2                             8
c      1/4  ---------------------------9---------------------------3
c                       1                             9
c      1/2  --------------------------10---------------------------2
c                       0                            10
c       1   --------------------------11---------------------------1
c
c
c    last modified 03/10/95...dave flittner
c    purpose: set earth radius equal to value in common block consts
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.
c--------------------------------------------------------------------
      implicit integer*4(i-n),real*8(a-h,o-z)
      real*8 p11(11),p11log(11),p101(101),p101lg(101),
     1 x11(11),x11log(11),x101(101),x101lg(101),delx(10),
     2 tmp10(10),tmp101(101),t11ave(11)
      dimension z11(11),plginv(11),ps(82),h(82)
      dimension cx(11),bparx(4),ct(11),bpart(4),cz(11),bparz(4)
c
      include "consts.inc"
      include "prints.inc"
c
c*****define presures (atm.) for layer ozone amounts and temperatures.
c     p11(11) is the at the bottom of the atmosphere
c
      do 10 i=1,11
         p11(i)=1.0*2.**(i-11.)
         p11log(i)=dlog(p11(i))
10    continue
c
c*****define pressures (atm.) at 101 levels for output ozone and
c     temperatures  p101(101) is at the bottom of the atmosphere
c
      do 20 i=1,101
         p101(i)=1.0*10.**((i-101.)/20.)
         p101lg(i)=dlog(p101(i))
20    continue
c
c*****compute cumulative ozone at 11 pressure levels & convert to atm-cm
c
      x11(1)=xtop/1000.
      x11log(1)=dlog(x11(1))
      do 110 i=2,11
         x11(i)=x11(i-1)+delx(i-1)/1000.
         x11log(i)=dlog(x11(i))
110   continue
c
c*****compute ozone weighted average temperatures at 11 pressure level
c
      tsum=tmptop*xtop/1000.
      t11ave(1)=tmptop
      do 150 i=2,11
         tsum=tsum+tmp10(i-1)*delx(i-1)/1000.
         t11ave(i)=tsum/x11(i)
150   continue
c
c*****set up boundary value parametes for cumulative ozone spline fit
c
      bparx(1)=0.
      bparx(2)=0.
      bparx(3)=1.
c
c*****f'(nx)=1.707*layer 0 ozone / total ozone
c
      deltax=x11log(11)-x11log(10)
      deltap=p11log(11)-p11log(10)
      fnx=1.707*delx(10)/1000./x11(11)
      bparx(4)=6.0/deltap*(fnx-deltax/deltap)
c
c*****evaluate cumulative ozone at the lowest 61 output preure level
c     using a spline fit to point at input pressure levels
c
      nx=11
      ic=10
      l1=41
      l2=101
      m=l2-l1+1
      call spline(p11log,x11log,nx,1.d30,fnx,cx)
      do 175 i = l1,l2
         call splint(p11log,x11log,cx,nx,p101lg(i),x101lg(i))
175   continue
c
c***** evaluate cumulative ozone at pnot
c
      xpnot=x11(11)
      if (pnot.ge.p11(11)) go to 180
      p0log=dlog(pnot)
      call splint(p11log,x11log,cx,nx,p0log,xp0log)
c     if(ier.ne.0) go to 900
      xpnot=dexp(xp0log)
  180 continue
      if (jprint(6).eq.1) write(33,2000)  xpnot
c
c*****evaluate cumultive ozone at uppermost 40 levels based on
c     slope computed from the 41st and 42nd levels
c
      slope=(x101lg(l1)-x101lg(l1+1))/(p101lg(l1)-p101lg(l1+1))
c
      l1m1=l1-1
      do 200 l=1,l1m1
      x101lg(l)=x101lg(l1)+slope*(p101lg(l)-p101lg(l1))
  200 continue
c
c***** convert log of cumulative ozone to cumulative ozone
c
      do 250 l=1,101
      x101(l)=dexp(x101lg(l))
  250 continue
c
c*****evaluate ozone-weighted average temperatures for atmospheres
c     above lowest 61 output pressure levels using a spline fit to
c     points at input pressure levels
c
      bpart(1)=0.
      bpart(2)=0.
      bpart(3)=0.
      bpart(4)=0.
      call spline(p11log,t11ave,nx,1.d30,1.d30,ct)
      do 275 i = l1,l2
         call splint(p11log,t11ave,ct,nx,p101lg(i),tmp101(i))
275   continue
c
c*****evaluate weighted temperatures in the upper most 40 layers
c     based on constant lapse rate of -1.5 degrees per layer)
195   do 300 l=1,l1m1
         tmp101(l)=tmp101(l1)-1.5*(l1-l)
 300  continue
c
c*****compute heights z11 (km.) corresponding to the 11 input pressure
c     levels p11 (in inverse order) using the equations :
c     dz/dlogp=n*r*t/g
c     g=g0*(r0/(r0+z))**2
c
c     constant is g0/(n*r); z11(1) is at bottom of atmosphere
c
      const=1.0/(-34.163)
      r0=r
      dzhalf=2.5
      z11(1)=0.
      plginv(1)=p11log(11)
      do 400 i=2,11
         j=12-i
         plginv(i)=p11log(j)
         dlogp=plginv(i)-plginv(i-1)
         z=z11(i-1)+dzhalf
         t=tmp10(j)
         dz=const*((r0+z)/r0)**2*t*dlogp
         z11(i)=z11(i-1)+dz
         dzhalf=dz/2.0
400   continue
c
c*****perform a spline fit to log(p) and z
c
      bparz(1)=0.
      bparz(2)=0.
      bparz(3)=0.
      bparz(4)=0.
c     if(ier.ne.0) go to 900
      call spline(z11,plginv,nx,1.d30,1.d30,cz)
c
c*****set up altitude arrays at 1 km constant intervals
c
      do 320 i=1,82
         h(i)=float(i-1)*1.0
 320  continue
c
c*****evaluate log of pressure at 1 km intervals from the ground to
c     the top of umkehr layer 9 based on spline fit
c
      m=z11(11)+1.
c     if(ier.ne.0) go to 900
      do 345 i = 1,m
         call splint(z11,plginv,cz,nx,h(i),ps(i))
345   continue
c
c*****evaluate log of pressues above umkehr layer 9 based on constant
c     scale height
c
      slope=(ps(m)-ps(m-1))/(h(m)-h(m-1))
      mp1=m+1
      do 450 i=mp1,82
      ps(i)=ps(m)+slope*(h(i)-h(m))
 450  continue
c
c***** convert log of pressure to pressure
c
      do 460 i=1,82
      ps(i)=dexp(ps(i))
  460 continue
c
c*****  dump out variout input & ouput data
c
      if (jprint(6).eq.1)
     1 write(33,6200) tmptop,xtop,(i,p11(i),x11(i),t11ave(i),z11(12-i),
     2 tmp10(i),delx(i),i=1,10),(i,p11(i),x11(i),t11ave(i),z11(12-i),
     3 i=11,11)
      if (jprint(6) .eq. 1) write(33,6210) (h(i),ps(i),i=1,82)
  900 continue
      return
2000  format(1x,'xpnot=',f10.5//)
6200  format(//' ozone,temperature & other variables in umkehr layers:'/
     1 /t28,'weighted',t51,'layer',t61,'layer'/1x,'index',t13,'p11',t23,
     2 'x11',t32,'temp',t43,'z11',t52,'temp',t61,'ozone'
     3 /t46,2f10.2/10(i3,2x,e10.4,f10.3,2f10.2/t46,2f10.2/),
     4 i3,2x,e10.4,f10.3,2f10.2)
6210  format(/'    pressures for first 82 km altitude :'/5(8x,'h',9x,
     1 'ps',2x)/17(5(f10.1,e12.4)/))
      end