c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/std_relayr.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine std_relayr(pnot,xpnot)
c ----------------------------------------------------------------
c      subroutine std_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(h,ps)
c
c
c     variable     type    i/o   description
c     --------     ----    ---   -----------
c via common input.cmn
c     xprf(11)      r*8     i    layer ozone amount (d.u) for
c                                input layers. (layer
c                                1 being the top of the atmos-
c                                phere). indexing of the layers are
c                                described in the table below
c via common input.cmn
c     tmpprf(11)    r*8     i    average temperatures (degrees kelvin)
c                                for input layers
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     x_101(101)     r*8    o     output cumulative ozone amounts
c                                (atm-cm) at levels corresponding
c                                to p101.
c
c
c     t_101(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     hhold(82)  r*8    o through atm_101.cmn: altitude array at 1 km 
c                                interval from ground to 82 km
c
c     pshold(82) r*8    o through atm_101.cmn: pressures corresponding 
c                                to the height 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 xprf      for
c                                  & t11ave        & tmpprt      z11
c                                                     1
c      1/1024 -------------------------1--------------------------11
c                       9                             2
c      1/512---------------------------2--------------------------10
c                       8                             3
c      1/256---------------------------3---------------------------9
c                       7                             4
c      1/128---------------------------4---------------------------8
c                       6                             5
c      1/64 ---------------------------5---------------------------7
c                       5                             6
c      1/32 ---------------------------6---------------------------6
c                       4                             7
c      1/16 ---------------------------7---------------------------5
c                       3                             8
c      1/8  ---------------------------8---------------------------4
c                       2                             9
c      1/4  ---------------------------9---------------------------3
c                       1                            10
c      1/2  --------------------------10---------------------------2
c                       0                            11
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
c     last modified July, 1998...def
c     purpose: allow for use of input profile with elevated surface.
c     This was used in the flux intercomparison.
c
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn.
c
c     last mod. Oct. 19, 1998...def
c     purpose: tidy-up user input profile option.  Use us76 atm
c     for t-p-z above and below the input profile with adjustments
c     at boundaries to match input profile. Will not use o3 weighted
c     temperatures.  Instead will use average of input temps for 
c     temp. to use for the abs coe calc.
c
c     last mod. Oct. 19, 1998...def
c     purpose: also save cumo3 so don't have to read the user profile
c     each time.
c
c     last mod. Nov. 20, 1998...def
c     purpose: fix use of us76 atmo for user profile. Limit surface 
c     pressure to be no greater than 1.0 atm.
c
c last mod. Nov. 23, 1998...def
c purpose: check to see that the cumulative ozone decreases with altitude.
c
c last mod. May 6, 1999...def
c purpose: pass pshold and hhold through contrl.cmn instead of command line.
c
c last mod. Nov. 2, 1999...def
c purpose: use .not. instead of .eqv. in logical testing
c
c last mod. Nov. 12, 1999...def
c purpose: use isothermal atmo. below the lowest point in the user profile 
c and not the standard atmo table.
c
c last mod. May 16, 2000...def
c purpose: add switch for using ozone weighted temperatures with the
c user profile.
c
c ----------------------------------------------------------------
      implicit none
c passed in
      real*8 pnot
c passed out
      real*8 xpnot
c commons
c passed out through atm_101.cmn are pshold and hhold
      include 'parameter.inc'
      include 'consts.cmn'
      include 'switches.cmn'
      include 'input.cmn'
      include 'atm_11.cmn'
      include 'atm_101.cmn'
c local
      integer l1,l2,m,i,l1m1,l,j,mp1
      real*8 p101(101),p101lg(101),twgt(n_101)
      real*8 dlogp,t,const,tbar
      real*8 deltax,deltap,fnx,p0log,xp0log,tsum,slope,r0,dzhalf,z,dz
      integer mode
      real*8 y21,y2n
      real*8 ptmp,xtmp1,xtmp2,htmp1,htmp2,xtmp3,xtmp4,no3_2,no3_1,h_no3
c
      if (jprint(6).eq.1) write(33,1000)
c
c*****define pressures (atm.) at 101 levels for output ozone and
c     temperatures  p101(101) is at the bottom of the atmosphere
c
      if(lv7tab)then
        do 20 i=1,n_101
           p101(i)=1.0*10.**((i-101.)/20.)
           p101lg(i)=dlog(p101(i))
20      continue
        do i=1,n_101
           p_101(i)=p101(i)
           p_101lg(i)=p101lg(i)
        enddo
      else
        call cp101()
      endif
c
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,n_p_lay
           if(prf_type.eq.1)then
             p11(i)=p_lay(i)
           else
             p11(i)=1.0*2.**(i-11.)
             p_lay(i)=p11(i)
           endif
           p11log(i)=dlog(p11(i))
10      continue
c
c*****compute cumulative ozone at 11 pressure levels & convert to atm-cm
c
        x11(1)=xprf(1)/1000.
        x11log(1)=dlog(x11(1))
        do 110 i=2,n_p_lay
           x11(i)=x11(i-1)+xprf(i)/1000.
           x11log(i)=dlog(x11(i))
110     continue
c
c*****compute ozone weighted average temperatures at 11 pressure level
c
        tsum=tmpprf(1)*xprf(1)/1000.
        t11ave(1)=tmpprf(1)
        do 150 i=2,n_p_lay
           tsum=tsum+tmpprf(i)*xprf(i)/1000.
           t11ave(i)=tsum/x11(i)
150     continue
c
c*****f'(nx)=1.707*layer 0 ozone / total ozone
c
        deltax=x11log(n_p_lay)-x11log(n_p_lay-1)
        deltap=p11log(n_p_lay)-p11log(n_p_lay-1)
        fnx=1.707*xprf(n_p_lay)/1000./x11(n_p_lay)
        if(prf_type.eq.0)then
c
c*****evaluate cumulative ozone at the lowest 61 output preure level
c     using a spline fit to point at input pressure levels
c
          l1=41
          l2=101
          m=l2-l1+1
          mode=0
          call spline(p11log,x11log,n_p_lay,1.d30,fnx,cx,mode)
          do 175 i = l2,l1,-1
             call splint(p11log,x11log,cx,n_p_lay,p_101lg(i),x_101lg(i),
     &mode)
175       continue
c
c*****evaluate cumultive ozone at uppermost 40 levels based on
c     slope computed from the 41st and 42nd levels
c
          slope=(x_101lg(l1)-x_101lg(l1+1))/(p_101lg(l1)-p_101lg(l1+1))
          l1m1=l1-1
          do 200 l=l1m1,1,-1
            x_101lg(l)=x_101lg(l1)+slope*(p_101lg(l)-p_101lg(l1))
  200     continue
        else
c
c use C,sigma for the top layer based upon the scale height
c from the layer below the top
c
          slope=(x11log(1)-x11log(2))/(p11log(1)-p11log(2))
          mode=0
          call spline(p11log,x11log,n_p_lay,1.d30,fnx,cx,mode)
          do i = 1,n_101
             if(p101lg(i).lt.p11log(1))then
               x_101lg(i)=x11log(1)+slope*(p101lg(i)-p11log(1))
             else
               call splint(p11log,x11log,cx,n_p_lay,p101lg(i),
     &x_101lg(i),mode)
             endif
          enddo
        endif
c
c***** evaluate cumulative ozone at pnot
c
        xpnot=x11(n_p_lay)
        if (pnot.lt.p11(n_p_lay))then
          p0log=dlog(pnot)
          call splint(p11log,x11log,cx,n_p_lay,p0log,xp0log,mode)
          xpnot=dexp(xp0log)
        endif
      if (jprint(6).eq.1) write(33,2000)  xpnot
c
c***** convert log of cumulative ozone to cumulative ozone
c
      do 250 l=1,n_101
      x_101(l)=dexp(x_101lg(l))
  250 continue
c
c*****set up altitude arrays at 1 km constant intervals from 0-81
c
      do 320 i=1,82
         hhold(i)=float(i-1)*1.0
 320  continue
      const=1.0/(-34.163)
      r0=r
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
        dzhalf=2.5
        z11(1)=0.
        plginv(1)=p11log(n_p_lay)
        do 400 i=2,n_p_lay
           j=n_p_lay+1-i
           plginv(i)=p11log(j)
           dlogp=plginv(i)-plginv(i-1)
           z=z11(i-1)+dzhalf
           t=tmpprf(j+1)
           dz=const*((r0+z)/r0)**2*t*dlogp
           plginv(i)=p11log(n_p_lay+1-i)
           z11(i)=z11(i-1)+dz
c for more rigerious calc. of z11 use the following two lines
c           dz2=-1./(dlogp*const*t/r0/r0+1./(r0+z11(i-1)))-r0-z11(i-1)
c           z11(i)=z11(i-1)+dz2
           dzhalf=dz/2.0
400     continue
        call set_p11_z11_splines()
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(n_p_lay)+1.
        if(m.gt.82)m=82
c
c*****evaluate log of pressues above umkehr layer 9 based on constant
c     scale height
c
        do 345 i = 1,82
           if(hhold(i).gt.z11(n_p_lay))then
             slope=(pshold(m)-pshold(m-1))/(hhold(m)-hhold(m-1))
             pshold(i)=pshold(m)+slope*(hhold(i)-hhold(m))
           else
             call splint(z11,plginv,cz,n_p_lay,hhold(i),pshold(i),mode)
           endif
345     continue
c
c***** convert log of pressure to pressure
c
      do 460 i=1,82
         pshold(i)=dexp(pshold(i))
  460 continue
        call tfromzp()
        if(gc_type.gt.1)call cmp_gc()
        if(lwgttmp)then
          if(lwgt11)then
c the standard method
            call cmptwgt11(p11log,t11ave,n_p_lay,p_101lg,twgt,n_101)
          else
            call cmptwgt101(p_101lg,t_101,x_101lg,n_101,twgt)
          endif
          do i=1,n_101
             t_101(i)=twgt(i)
          enddo
        endif
c
c*****  dump out variout input & ouput data
c
      if (jprint(6).eq.1)then
         write(33,6201)
         write(33,6202) tmpprf(1),xprf(1)
         write(33,6203) (i,p11(i),x11(i),t11ave(i),
     & z11(n_p_lay+1-i),tmpprf(i+1),xprf(i+1),i=1,n_p_lay-1)
         write(33,6204) (i,p11(i),x11(i),t11ave(i),z11(n_p_lay+1-i),
     & i=n_p_lay,n_p_lay)
      endif
c     1 write(33,6200) tmpprf(1),xprf(1),(i,p11(i),x11(i),t11ave(i),
c     2 z11(n_p_lay+1-i),tmpprf(i+1),xprf(i+1),i=1,n_p_lay-1),
c     3 (i,p11(i),x11(i),t11ave(i),z11(n_p_lay+1-i),i=n_p_lay,n_p_lay)
c      if (jprint(6) .eq. 1) write(33,6210) (hhold(i),pshold(i),i=1,82)
      if (jprint(6) .eq. 1) then
         write(33,*)'p_101lg(i),t_101(i,1),x_101lg(i),x_101(i)'
         do i=1,101
            write(33,'(1P4E12.4)')p_101lg(i),t_101(i),x_101lg(i),
     &x_101(i)
         enddo
         write(33,*)'relayr. m=',m
         write(33,*)'hhold(i),no3,pshold(i),tbar'
         write(33,*)' [km]  [#/cc]  [atm]   [K]'
c for the first point (z=0) we will get the o3 layer amount and the
c value of no3(hhold(2)) and assume that have a scale height in the lowest layer
            i=2
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(p_101lg,x_101lg,n_101,y21,y2n,y2_101,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp2,mode)
            ptmp=log(pshold(i))
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp3,mode)
            ptmp=log(pshold(1))
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp4,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2_101,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp2,mode)
      no3_2=(exp(xtmp2)-exp(xtmp1))/(htmp1-htmp2)
            H_no3=0.5/((exp(xtmp4)-exp(xtmp3))/no3_2-1.0d0)
      no3_2=no3_2*2.68684e14
      no3_1=no3_2*exp(1.0d0/H_no3)
c 
         do i=1,81
            tbar=(hhold(i+1)-hhold(i))/dlog(pshold(i+1)/pshold(i))/
     &(const*(1.0+0.5*(hhold(i)+hhold(i+1))/r0)**2)
c also compute the o3 concentration
            if(i.eq.1)then
              write(33,'(1P4E12.4)')hhold(i),no3_1,
     &pshold(i),tbar
            else              
              mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(p_101lg,x_101lg,n_101,y21,y2n,y2_101,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp2,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2_101,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp2,mode)
            write(33,'(1P4E12.4)')hhold(i),
     &(exp(xtmp2)-exp(xtmp1))/(htmp1-htmp2)*2.68684e14,
     &pshold(i),tbar
            endif
         enddo
c now for the last point
            i=82
c also compute the o3 concentration
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(p_101lg,x_101lg,n_101,y21,y2n,y2_101,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2_101,n_101,ptmp,xtmp2,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2_101,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2_101,82,ptmp,htmp2,mode)
            write(33,'(1P4E12.4)')hhold(i),
     &(exp(xtmp2)-exp(xtmp1))/(htmp1-htmp2)*2.68684e14,
     &pshold(i),tbar
      endif
  900 continue
      return
1000  format('relayr. debugging output--------')
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)
6201  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 '          [atm]   [atm-cm]    [K]        [km]      [K]       ',
     4 '[DU]')
6202  format(t46,2f10.2)
6203  format((i3,2x,e10.4,f10.3,2f10.2/t46,2f10.2))
6204  format(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
c