;=====================================================
;Script for getting higher order derivatives using SPH
;Get dates that satisfy certain conditions
;
;Script provided by Dr. Dennis Shea 
;====================================================
begin
 f = addfile("mslp_jan_2017.nc","r")
 p = f->mslp           ; ordered North-to-South

 if (p&lat(0) .gt. p&lat(1)) then
     p = p(:,::-1,:)      ; ordered South-to-North [required by spherical harmonics]
     print("pressure array reordered South-to-North")
     print("----------")
 end if

 nmsg = num(ismissing(p))   ; count missing values 
 if (nmsg.gt.0) then
     print("EXIT: missing values encountered: nmsg="+nmsg)
     exit
 end if

 printVarSummary(p)
 printMinMax(p,0)
 print("----------")


; Compute Laplacian
; gradsf is a procedure which requires that the return arrays be preallocated
; create arrays to hold output, same size and type as input

 dimp = dimsizes(p)
 dpdx = new(dimp, typeof(p), p@_FillValue) ; dp/dlon => dp/dx
 dpdy = new(dimp, typeof(p), p@_FillValue) ; dp/dlat => dp/dy 

; procedure gradsf will return dp/dx, dp/dy in the preallocated arrays

 gradsf (p,  dpdx, dpdy) 

 dpdx@long_name = "Zonal Pressure gradient: dpdx"
 dpdx@units     = "Pa/m"                        
 copy_VarCoords(p, dpdx)

 dpdy@long_name = "Meridionalgradient: dpdy"
 dpdy@units     = "Pa/m"  
 copy_VarCoords(p, dpdy)


 work = new(dimp, typeof(p), p@_FillValue) ; placeholder for temporary values
 px2  = new(dimp, typeof(p), p@_FillValue) ; d2(px)/dx2 
 py2  = new(dimp, typeof(p), p@_FillValue) ; d2(py)/dy2 

 gradsf (dpdx,  px2, work) 
 gradsf (dpdy, work, py2 ) 

 px2@long_name = "d2(p)/d2lon"
 px2@units = "..."
 copy_VarCoords(p, px2)

 py2@long_name = "d2(p)/d2lat"
 py2@units = "..."
 copy_VarCoords(p, py2)
 print("----------")

 printVarSummary(px2)
 print("----------")
 printVarSummary(py2)
 print("----------")
 
;=======================
;
;=======================
 px2py2 = px2*py2
 px2py2@long_name = "product: [d2(p)/d2lon]*[d2(p)/d2lat]"
 px2py2@units     = "..."
 copy_VarCoords(p, px2py2)
 printVarSummary(px2py2)

;=====================================
;Get Dates
;dpdx = 0 ; dpdY = 0; px2py2 < 0
;H Bluestein; Col characteristics pp.43
;
;The gridpoint must satisfy all three
;======================================

;Create a subregion

  latS = 5.0
  latN = 30.0
  lonL = 120.0
  lonR = 160.0

  dpdx_region=dpdx(:,{latS:latN},{lonL:lonR})  
  dpdy_region=dpdy(:,{latS:latN},{lonL:lonR}) 
  px2py2_region=px2py2(:,{latS:latN},{lonL:lonR}) 

;======================
;detect presence of col
;======================
  dpdx_region_1d = ndtooned(dpdx_region)
  dpdy_region_1d = ndtooned(dpdy_region)
  px2py2_region_1d = ndtooned(px2py2_region)

  ncol  = ind_resolve((ind(dpdx_region_1d.eq. 0) .and. ind(dpdy_region_1d .eq. 0) .and. ind(px2py2_region_1d.lt. 0)),dimsizes(dpdx_region))
  
  printVarSummary(ncol)

  print(ncol)
  print("=====")

  if(ismissing(ncol(0))) then
     print("++++++++++++++")
     print(" NO values exceed threshold="+threshold)
     print("++++++++++++++")
  exit
  end if

;***************************************************
;print dates
;***************************************************

      TIME  = dpdx_region&time             ; all times
      YMD  = cd_calendar(TIME,-2)          ; all yyyy-mm-dd
      ymd   = YMD(ncol)                 ; dates that satisfy criteria


     undef("add_hyphen")
      function add_hyphen(yyyymmdd[*]:integer)
      local yyyy, mmdd, mm, dd
      begin
      yyyy    = yyyymmdd/10000
      mmdd    = yyyymmdd-(yyyy*10000)
      mm      = mmdd/100
      dd      = mmdd-(mm*100)
      return(sprinti("%04d",yyyy) + "-" + sprinti("%02d",mm) + "-" + sprinti("%02d",dd))
      end

      yymd = add_hyphen(ymd)
      print(yymd)
      asciiwrite("all_withrain_"+year+".txt",yymd)

end