subroutine sofie_find_ice(nz,z,ext9,ext10,zlow,t,tmax,
     &   ztop,zmax,zbot,pmc_flag)

c------------------------------------------------------------------------------------------------
c Routine analyzes SOFIE extinction profiles and finds the ice layer top, peak, and bottom altitudes.
c
c Input:
c   nz.........# layers in profile
c   z..........altitude grid, km
c   ext9.......band  9 extinction profile (1/km)
c   ext10......band 10 extinction profile (1/km)
c   zlow.......minimum altitude to consider in looking for PMCs (km, recommended: 75 km)
c   t..........temperature profile (K)
c   tmax.......maximum temperature for it to be ice (K, recommend: 180 K with current warm bias)
c 
c Output:
c   ztop.......ice layer top altitude (km,  zmax = 0 if no ice is found)
c   zmax.......altitude of the peak band 9 extinction (km)
c   zbot.......ice layer bottom altitude (km)
c   pmc_flag... flag set to 1 if PMC   0 if no PMC
c
c Notes:  the routine assumes taht altitude goes from top - down,  i.e. z(1) > z(2)
c
c Source: Mark Hervig, GATS Inc.
c
c Revisions:
c   Mark Hervig, June 17, 2008, cleaned up, trying to make consistent with IDL version.
c
c------------------------------------------------------------------------------------------------

c      IMPLICIT REAL*8 (A-H,O-Z)
      integer pmc_flag
      dimension z(nz), t(nz),ext9(nz), ext10(nz),kice(900)
      
      ztop   = 0.
      zmax   = 0.
      zbot   = 99.
      
      e9max  = 0.
      
      imax   = 0    ! index of Zmax
      itop   = 0    ! index od Ztop
      ibot   = 0    ! index of Zbot
      kmax   = 0    ! index of Zmax
      ktop   = 0    ! index od Ztop
      kbot   = 0    ! index of Zbot
      
      bad    = 0.
      
      enoi   = 1e-7  ! the bit-noise in extinction (1/km)   
      minlay = 7     ! must have this many layers to call it ice
      maxgap = 6     ! don't allow gaps of more than this many points in ice layer
      nice   = 0     ! # of ice layers found
      
      R910_mn = 1.3 ! min band 9/10 ratio that is consistent with ice
      R910_mx = 2.4 ! max band 9/10 ratio that is consistent with ice
      
c-    find altitudes where the band 9/10 extinction ratio is consistent with ice
  
      pmc_flag=1

c      print*, "in sofie find ice"
c      stop 
      do i = 1,nz
      
        R910 = ext9(i) / ext10(i)       
    
        if (R910 .ge. R910_mn .and. R910 .le. R910_mx .and.  
     #      ext9(i) .gt. enoi .and. ext10(i) .gt. enoi .and. 
     #      z(i) .gt. zlow .and. z(i) .lt. 97.0 .and. t(i) .lt. tmax) 
     #     then 
     
          nice = nice + 1
          kice(nice) = i
                                     
          if (ext9(i) .gt. e9max) then  ! extinction peak Z
            zmax  = z(i)
            imax  = i
            kmax  = nice
            e9max = ext9(i)
          endif
          
          if (z(i) .gt. ztop) then  ! ice layer top Z
            ztop = z(i)  
            itop = i
            ktop = nice
          endif
          
          if (z(i) .lt. zbot) then  ! ice layer bottom Z
            zbot = z(i)
            ibot = i
            kbot = nice
          endif
         
        endif
      enddo  ! over altitude 
         
c-   if there are enough layers that looks like ice, go to work,
c    look for continuity from peak up & peak down, 
c    cut off the layer if a gap > maxgap is encountered

      if (nice .ge. minlay) then 
     
        zhi = 0.
        zlo = 0.
     
        do j = kmax,2,-1                   ! search from peak up       
          gap = abs(kice(j) - kice(j-1))

          if (gap .ge. maxgap) then 
            zhi  = z(kice(j))
            ztop = zhi
            itop = kice(j)
            ktop = j
            go to 1
          endif       
        enddo
    1   continue
     
        do j = kmax,nice-1                  ! search from peak down       
          gap = abs( kice(j) - kice(j+1) )

          if (gap .ge. maxgap) then
            zlo  = z(kice(j))
            zbot = zlo
            ibot = kice(j)
            kbot = j
            go to 2
         endif
        enddo
    2   continue
          
       if (zhi .gt. 0 .or. zlow .gt. 0) then
      
         nice = ibot - itop + 1       
         nbad = 0     
       
         do i = itop,ibot     
           R910 = ext9(i) / ext10(i)        
           if (R910 .le. R910_mn .or. R910 .ge. R910_mx) nbad = nbad + 1
         enddo
        
         bad = nbad / nice
       endif
       
      endif
            
      if (nice .lt. minlay .or. bad .gt. 0.20) then
        ztop = 0.
        zmax = 0.
        zbot = 0.
c        kice = 0
        kmax = 0.
        nice = 0
        pmc_flag=0
      endif
      if(pmc_flag .eq. 1 )  print*, 'found ice', ztop, zmax, zbot   

c       stop 
c-    done

      return
      end