c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/user_relayr.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine user_relayr(pnot,xpnot)
c ----------------------------------------------------------------
c--------------------------------------------------------------------
      implicit none
c commons
c passed out through atm_101.cmn are pshold and hhold
      include 'parameter.inc'
      include 'consts.cmn'
      include 'switches.cmn'
      include 'us76atm.cmn'
      include 'atm_101.cmn'
c passed in
      real*8 pnot
c passed out
      real*8 xpnot
c local
      integer i,l,j
      real*8 const,tbar
      real*8 p0log,slope,r0
      integer muser,nuser,mode,mode0
      parameter (muser=135)
      real*8 zuser(muser),puser(muser),o3user(muser),tuser(muser)
      real*8 cumo3(muser),y2user(muser),y21,y2n,junk
      real*8 toto3,H_cumo3,factor,twgt(muser)
      real*8 ptmp,xtmp1,xtmp2,htmp1,htmp2,xtmp3,xtmp4,no3_2,no3_1,h_no3
      save puser,zuser,tuser,o3user,toto3,nuser,cumo3,
     &twgt
c
      if (jprint(6).eq.1) write(33,1000)
c
c*****define pressures (atm.) at 101 levels for output ozone and
c     temperatures  p_101(101) is at the bottom of the atmosphere
c
      call cp101()
c input the user profile from the file user.prf
        call rdprof(puser,zuser,tuser,o3user,
     &toto3,muser,nuser)
        if(puser(nuser).lt.pnot)then
          write(6,*)'relayr. WARNING!!! WARNING!!!'
          write(6,*)'surface pressure of the user profile ',
     &'(puser(nuser)) is less than the'
          write(6,*)'surface pressure (pnot) in the PROF.DAT file.'
          write(6,*)'Should have pnot less than or equal to',
     &' puser(nuser)'
          write(6,*)'pnot = ',pnot
          write(6,*)'puser(nuser)= ',puser(nuser)
        endif
c
        cumo3(1)=0.0
        write(6,*)'level puser zuser cumo3'
        write(6,'(i3,3e12.4)')1,puser(1),zuser(1),cumo3(1)
        do i=2,nuser
           cumo3(i)=cumo3(i-1)+
     &0.5*(zuser(i-1)-zuser(i))*(o3user(i-1)+o3user(i))
           write(6,'(i3,3e12.4)')i,puser(i),zuser(i),cumo3(i)
        enddo
c check to see that input profile is not more than the input
c total amount
        write(6,2990)toto3,cumo3(nuser)
2990    format('the input total column amount was ',e14.6,
     &' and the integrated amount was ',e14.6)
        if(toto3.le.cumo3(nuser))then
          factor=toto3/cumo3(nuser)
c norm the profile
          do i=2,nuser
             cumo3(i)=cumo3(i)*factor
          enddo
c use scale height for the first point
          H_cumo3=(zuser(2)-zuser(3))/log(cumo3(3)/cumo3(2))
c check to see that the scale height is greater than 0.  If it is not,
c then set to 5 km.
          if(H_cumo3.le.0.)then
            write(6,*)'relayr. H_cumo3 =',H_cumo3
            write(6,*)'relayr. will set to 5 km'
            H_cumo3=5.0
          endif
          cumo3(1)=cumo3(2)*exp((zuser(2)-zuser(1))/H_cumo3)
          write(6,3000)factor,toto3
        else
          cumo3(1)=toto3-cumo3(nuser)
          do i=2,nuser
             cumo3(i)=cumo3(i)+cumo3(1)
          enddo
        endif
        do i=1,nuser
           puser(i)=log(puser(i))
           cumo3(i)=log(cumo3(i))
        enddo
c spline from the user profile
        y21=1.0d30
        y2n=1.0d30
        mode=0
        call spline(puser,cumo3,nuser,y21,y2n,y2user,mode)
        do i = 1,n_101
           if(p_101lg(i).lt.puser(1))then
c above the top of the input user profile
c assume ln(cumo3(ln(p)))=linear function using the average
c slope from the last 2 points
             slope=0.5*((cumo3(2)-cumo3(1))/(puser(2)-puser(1))+
     &(cumo3(3)-cumo3(2))/(puser(3)-puser(2)))
             x_101lg(i)=cumo3(1)+slope*(p_101lg(i)-puser(1))
           else if(p_101lg(i).gt.puser(nuser))then
c below the bottom of the input user profile
c extrapolate using constant density taken from the average of the last
c 10 data points of the profile and a pressure scale height of 7 km
             x_101lg(i)=0.0
             do j=nuser-9,nuser
                x_101lg(i)=x_101lg(i)+o3user(j)*0.1
             enddo
             x_101lg(i)=log(-7.0*(puser(nuser)-p_101lg(i))*x_101lg(i)
     &+exp(cumo3(nuser)))
           else
             call splint(puser,cumo3,y2user,nuser,p_101lg(i),
     &x_101lg(i),mode)
c check to see that x_101 increases with decreasing p_101
             if(i.gt.1)then
               if(x_101lg(i).lt.x_101lg(i-1))then
c use linear for this point
c                 write(6,*)'relayr. x dec. ',i,exp(x_101lg(i-1)),
c     &exp(x_101lg(i))
                 mode0=mode
                 mode=4
                 call splint(puser,cumo3,y2user,nuser,p_101lg(i),
     &x_101lg(i),mode)
c set back to cubic
                 mode=mode0
               endif
             endif
           endif
c           write(6,*)'relayr.p x',i,exp(p_101lg(i)),exp(x_101lg(i))
        enddo
        p0log=dlog(pnot)
        xpnot=dexp(cumo3(nuser))
        if(p0log.ne.puser(nuser))then
          mode=0
          y21=1.0d30
          y2n=1.0d30
          call spline(p_101lg,x_101lg,n_101,y21,y2n,y2user,mode)
          call splint(p_101lg,x_101lg,y2user,n_101,p0log,xpnot,mode)
          xpnot=dexp(xpnot)
        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
        y21=1.0d30
        y2n=1.0d30
        do i = 1,82
           if(hhold(i).gt.zuser(1))then
c altitude above the user input profile
c interpolate in standard atmosphere tables and match the pressure
c at the top of the user input profile
             mode=0
             call spline(z76,p76,82,y21,y2n,y2user,mode)
             call splint(z76,p76,y2user,82,zuser(1),junk,mode)
             call splint(z76,p76,y2user,82,hhold(i),pshold(i),mode)
             pshold(i)=pshold(i)+(puser(1)-junk)
           else if(hhold(i).lt.zuser(nuser))then
c altitude below the user input profile
c use constant pressure scale height between nuser point and
c the surface which will be set at 1.0 atm.
             pshold(i)=puser(nuser)*(1.0-
     &(zuser(nuser)-hhold(i))/zuser(nuser))
c for now limit the pressure to be at most 1.0atm
             if(pshold(i).gt.0.0d0)pshold(i)=0.0d0
           else
             mode=0
             call spline(zuser,puser,nuser,y21,y2n,y2user,mode)
             call splint(zuser,puser,y2user,nuser,hhold(i),pshold(i),
     &mode)
           endif
c           write(6,*)'relayr. hhold pshold',hhold(i),exp(pshold(i))
        enddo
c compute the altitude of each one of the pressure levels
        y21=1.0d30
        y2n=1.0d30
        i=n_101
        do while (i.gt.0)
           if(p_101lg(i).gt.puser(nuser))then
c altitude below the user input profile
c user constant scale height between nuser point and the
c surface which will be set at 1.0 atm
             alt_101(i)=p_101lg(i)/puser(nuser)*zuser(nuser)
           else if(p_101lg(i).lt.puser(1))then
             do j=i,1,-1
c altitude above the user input profile
c interpolate in standard atmosphere tables and match the pressure
c at the top of the user input profile
                mode=0
                call spline(p76,z76,82,y21,y2n,y2user,mode)
                call splint(p76,z76,y2user,82,puser(1),junk,mode)
                call splint(p76,z76,y2user,82,p_101lg(i),alt_101(i),
     &mode)
                alt_101(i)=zuser(1)+(alt_101(i)-junk)
             enddo
             goto 459
           else
             mode=0
             call spline(puser,zuser,nuser,y21,y2n,y2user,mode)
             call splint(puser,zuser,y2user,nuser,p_101lg(i),alt_101(i),
     &mode)
           endif
c           write(6,*)'user_relayr. alt_101 p_101',alt_101(i),exp(p_101lg(i))
           i=i-1
        enddo
459     if(gc_type.gt.1)call cmp_gc()
c
c***** convert log of pressure to pressure
c
      do 460 i=1,82
         pshold(i)=dexp(pshold(i))
  460 continue
c
c create a temperature profile at the 101 pressure levels
c this temperature profile will then be used either directly to compute
c the absorption coe. for O3, or be used in the O3 weighting scheme
c to compute an O3 weighted temperature profile that will be used to 
c compute the O3 abs. coe.
c Will assume T(alog(p))=linear function
c
        mode=4
        y21=1.0d30
        y2n=1.0d30
        do i = 1,n_101
           if(p_101lg(i).lt.puser(1))then
c interpolate in standard atmosphere tables
             call spline(p76,t76,82,y21,y2n,y2user,mode)
             call splint(p76,t76,y2user,82,p_101lg(i),t_101(i),mode)
c shift standard atmo to match input profile at last point
             call splint(p76,t76,y2user,82,puser(1),junk,mode)
             t_101(i)=t_101(i)+tuser(1)-junk
           else if(p_101lg(i).gt.puser(nuser))then
c interpolate in standard atmosphere tables
             call spline(p76,t76,82,y21,y2n,y2user,mode)
             call splint(p76,t76,y2user,82,p_101lg(i),t_101(i),mode)
c shift standard atmo to match input profile at last point
             call splint(p76,t76,y2user,82,puser(nuser),junk,mode)
             t_101(i)=t_101(i)+tuser(nuser)-junk
           else
             call spline(puser,tuser,nuser,y21,y2n,y2user,mode)
             call splint(puser,tuser,y2user,nuser,p_101lg(i),
     &t_101(i),mode)
           endif
        enddo
c
        if(lwgttmp)then
c weight the temperature profile using the cummulative o3 profile
          call cmptwgt101(p_101lg,t_101,x_101lg,n_101,twgt)
          do i=1,n_101
             t_101(i)=twgt(i)
          enddo
c average the temp for each layer
          do i=n_101,2,-1
             t_101(i)=0.5d0*(t_101(i)+t_101(i-1))
          enddo
        endif
c
c*****  dump out variout input & ouput data
c
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),x_101lg(i),x_101(i)'
         do i=1,n_101
            write(33,'(1P4E12.4)')p_101lg(i),t_101(i),x_101lg(i),
     &x_101(i)
         enddo
         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,y2user,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp2,mode)
            ptmp=log(pshold(i))
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp3,mode)
            ptmp=log(pshold(1))
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp4,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2user,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2user,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2user,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,y2user,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp2,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2user,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2user,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2user,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,y2user,mode)
            ptmp=log(pshold(i)*0.999)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp1,mode)
            ptmp=log(pshold(i)*1.001)
            call splint(p_101lg,x_101lg,y2user,n_101,ptmp,xtmp2,mode)
            mode=0
            y21=1.0d30
            y2n=1.0d30
            call spline(pshold,hhold,82,y21,y2n,y2user,mode)
            ptmp=(pshold(i)*0.999)
            call splint(pshold,hhold,y2user,82,ptmp,htmp1,mode)
            ptmp=(pshold(i)*1.001)
            call splint(pshold,hhold,y2user,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//)
3000  format(1x,'relayr. The user profile must be adjusted by a',
     &' factor of ',1PE14.6,/,'  to match the input total column',
     &' ozone amount of ',0Pf8.6)
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
c