subroutine stnp81(ilat, estozn, lprint, stnpro)
c
c***********************************************************************
c
c     stnp81
c
c     march, 2001  by charlie wellemeyer SSAI
c
c     purpose
c        determine a standard profile associated with a specific total
c        ozone and a single latitude band.
c
c     method
c        interpolates to input total ozone within a single
c        latitude band.  pure low latitude profiles
c        are used for latitudes less than 30 degrees and 
c        pure low latitude profiles are used for latitudespure
c        high latitude profiles are used for latitudes greater
c        than 75 degrees. all other profiles are interpolated in
c        latitude with the mid latitude profiles, which are
c        representative of 45 degrees latitude.
c
c  low to mid latitude interpolations
c
c         1    2    3    4    5    6    7    8    9   10
c low              225  275  325
c mid              225  275  325  375  425  475  525  575
c high   125  175  225  275  325  375  425  475  525  575
c
c
c     modification
c        none
c
c     variables
c       name      type  i/o   description
c       ----      ----  ---   ------------
c
c      arguments
c       ilat       i*4    i    latitude index
c       estozn     r*8    i    total ozone interpolation target
c       lprint(30) l*4    i    logical switches for debug printout
c       stnpro(81) r*8    o    interpolated standard ozone profile
c
c      internal
c       indoz      i*4         latitude index for lower interpolation node
c       indm1      i*4         latitude index for upper interpolation node
c
c***********************************************************************
c
c     implicit none
c
c     -- input parameters
c
      integer ilat
      real estozn
      logical lprint(30)
c
c  internal parameters
c
      real stnpr1(81), stnpr2(81), ozind(21), tot,
     1     ozone, abslat
      real thrty/30.0/,frtyfv/45.0/,sixty/60.0/
      integer indl1, indl2, indm1, indm2, indh1, indh2, indoz, ilev, i
      real tot1, tot2, total
c
c     output parameters
c
      real stnpro(81)
c
      data ozind/225.0,275.0,325.0,225.0,275.0,325.0,375.0,
     1           425.0,475.0,525.0,575.0,125.0,175.0,225.0,
     2           275.0,325.0,375.0,425.0,475.0,525.0,575.0/
c
      include 'prof81.com'
c
c  estimate standard profile for given ozone and latitude band
c
      ozone=estozn
c
      if(ilat.eq.2) then 
        indm1 = 4
        do 200 indoz=4,10                     
          if(ozone.gt.ozind(indoz)) indm1 = indoz                   
  200   continue                                                  
        indm2 = indm1 + 1                                      
c                                                           
        tot = 0.0
        tot1 = 0.0
        tot2 = 0.0
        do 300 ilev=1,81                                  
          stnpr1(ilev) = terp(ozone,ozind(indm1),ozind(indm2),       
     1         ozp81(ilev,indm1),ozp81(ilev,indm2))              
          tot=tot+stnpr1(ilev)
          tot1=tot1+ozp81(ilev,indm1)
          tot2=tot2+ozp81(ilev,indm2)
  300   continue                                                 
c
        if(lprint(25)) write(6,7010) indm1, indm2, tot, tot1,
     1     tot2, stnpr2
 7010   format('1 stnprf:  indm1, indm2=',2i5,'  toz =',3f8.2,
     1     '   mid-lat prof:',/,9(10f7.2/))
c
      endif                                                    
c                                                            
c  low to mid latitude interpolations                     
c                                                       
c         1    2    3    4    5    6    7    8    9   10
c low              225  275  325 
c mid              225  275  325  375  425  475  525  575
c high   125  175  225  275  325  375  425  475  525  575
c
      if(ilat.eq.1) then                          
        indl1 = 1
        do 350 indoz=1,2                     
          if(ozone.gt.ozind(indoz)) indl1 = indoz                    
  350   continue                                                  
        indl2 = indl1 + 1                                      
c                                   
        tot = 0.0
        do 400 ilev=1,81         
          stnpr1(ilev) = terp(ozone,ozind(indl1),ozind(indl2),       
     1         ozp81(ilev,indl1),ozp81(ilev,indl2))             
          tot=tot+stnpr1(ilev)
  400   continue                                                
c
        if(lprint(25)) write(6,7020) indl1, indl2, tot, stnpr1
 7020   format('2 stnprf:  indl1, indl2=',2i5,'   toz =',f8.2,
     1     '   low-lat prof',/,9(10f7.2/))
c
      endif
c                                                        
c  mid to high latitude interpolations                  
c                                                      
c         1    2    3    4    5    6    7    8    9   10
c low              225  275  325
c mid              225  275  325  375  425  475  525  575
c high   125  175  225  275  325  375  425  475  525  575
c
      if(ilat.eq.3) then                                            
        indh1 = 12
        do 600 indoz=12,20                         
          if(ozone.gt.ozind(indoz)) indh1 = indoz
  600   continue                              
        indh2 = indh1 + 1                
c                                     
        tot = 0.0
        do 700 ilev=1,81          
          stnpr1(ilev) = terp(ozone,ozind(indh1),ozind(indh2),        
     1         ozp81(ilev,indh1),ozp81(ilev,indh2))              
          tot=tot+stnpr1(ilev)
          tot1=tot1+ozp81(ilev,indh1)
          tot2=tot2+ozp81(ilev,indh2)
c     if(lprint(25)) write(6,'(6f10.2)') stnpr1(ilev),ozone,
c    1     ozind(indh1),ozind(indh2),ozp81(ilev,indh1),
c    2     ozp81(ilev,indh2)
  700   continue                                                  
c
        if(lprint(25)) write(6,7030) indh1, indh2, tot, tot1,
     1     tot2, stnpr1
 7030   format('3 stnprf:  indh1, indh2=',2i5,'   toz =',3f8.2,
     1     '   hi-lat prof',/,9(10f7.2/))
c                                                       
      endif                                            
c
      total = 0.0
      do 900 i=1,81
        stnpro(i) = stnpr1(82-i)
        total = total + stnpro(i)
  900 continue
      if(abs(total-ozone).gt.0.1) then
        print*,'non-linearity error in stnprf'
        print*,total,ozone,abs(total-ozone)
      endif
c
      if(lprint(25)) write(6,7050) total, stnpro
 7050 format('5 stnprf: Total ozone =',f8.2,
     1     '   retrieved standard profile:',/,9(10f7.2/))
c
  999 continue
      return                                                         
c                                                                  
c  formats                                                       
c                                                             
      end                                                   
      function terp(x,x1,x2,y1,y2)
c
      real x,x1,x2,y1,y2
c                           
      terp = (x - x1) / (x2 - x1) * (y2 - y1) + y1 
c                                          
      return  
      end