undef("first_occurence_ind")
function first_occurence_ind(T:numeric, testValue[1]:numeric)
;
; Must check for special cases:
;      [1] all values.le.testValue  [ Cold regions ]
;      [2] all values.ge.testValue  [ Warm regions ]
;
; +++++ ncl-talk question ++++++++++++++
; 5 Dec 2019
; I'm dealing with the three-dimensional grid temperature data (ordered by time, lat, lon). 
; I want to confirm the first day (position in the time dimension) in each lat/lon that 
; the temperature is greater than 0. It's hard for me to think about how to write the "if" 
; *statement. Are there some relevant statistics function in NCL? 
; ++++++++++++++++++++++++++++++++++++++
;
local t, t1,it, dimt, rankt, ntim, nsta, nlat, mlon  
begin
  t  = T                                                ; local [deleted upon exit]
  if (.not.isatt(t,"_FillValue"))
      t@_FillValue = 1e20
  end if

  dimt = dimsizes(t)
  rankt= dimsizes(dimt)
  ntim = dimt(0)

  if (rankt.eq.1) then                                   ; T[*]  => T(ntim)
      it = new(1, "integer", -999)
      it@long_name = "temporal index for Cross-Over"
      if (all(t.le.testValue) .or. all(t.ge.testValue)) then
          return(it)
      end if
      t  = where(t.le.testValue, t@_FillValue, t)
      t  = where(.not.ismissing(t),testValue,t)          ; all 't' = testValue
      it = dim_minind(t,0)                            ; indices of 1st occurence of testValue
      return(it)
  elseif (rankt.eq.2) then                               ; T[*][*] => T(time,nStations)
      nsta = dimt(1)                                
      it   = new(nsta, "integer", -999)
      it@long_name = "temporal index for Cross-Over"
      copy_VarCoords(t(0,:),it)
      if (all(t.le.testValue) .or. all(t.ge.testValue)) then
          return(it)
      else
          do nst=0,nsta-1
             t1 = t(:,nst)
             if (.not.(all(t1.le.testValue) .or. all(t1.ge.testValue))) then
                 t1 = where(t1.le.testValue, t1@_FillValue, t1)
                 t1 = where(.not.ismissing(t1),testValue,t1)      
                 it(nst) = dim_minind(t1,0)                ; indices of 1st occurence of testValue
             end if
          end do
          return(it)
      end if
  elseif (rankt.eq.3) then                               ; T[*][*][*] => T(time,lat,lon)
      nlat = dimt(1)
      mlon = dimt(2)
      it   = new((/nlat,mlon/), "integer", -999)
      copy_VarCoords(t(0,:,:),it)
      it@long_name = "temporal index for Cross-Over"

      if (.not.(all(t.le.testValue) .or. all(t.ge.testValue))) then
          do nl=0,nlat-1
            do ml=0,mlon-1
               t1 = t(:,nl,ml)
               if (.not.(all(t1.le.testValue) .or. all(t1.ge.testValue))) then
                   t1 = where(t1.le.testValue, t1@_FillValue, t1)
                   t1 = where(.not.ismissing(t1),testValue,t1)      
                   it(nl,ml) = dim_minind(t1,0)                ; indices of 1st occurence of testValue
               end if
            end do
          end do
          return(it)
      end if
  end if                    
end