;*****************************************************
;This extracts dates satistifying multiple conditions
;Lyndon Mark Olaguera
;Created: May 27, 2019
;*****************************************************

begin
;************************************************
; variable and file handling
;************************************************
  f  = addfile("test_file_olr_wind_deform.nc","r")       ; open netcdf file

  olr    = f->olr(:,:,:)                  
  olr   := short2flt(olr)                 ; float
  u      = f->uwnd(:,:,:)                 ; float 
  stretch = f->stretch_sph(:,:,:)        ; float

;**************************************************
;Extract over a region
;**************************************************
  latS = 13.0
  latN = 22.0
  lonL = 120.0
  lonR = 122.5

  olr_region=olr(:,{latS:latN},{lonL:lonR})     ;select subregion
  printVarSummary(olr_region)
  printMinMax(olr_region,0)

  u_region=u(:,{latS:latN},{lonL:lonR})
  stretch_region=stretch(:,{latS:latN},{lonL:lonR})

;**************************************************
;Find date when conditions are satisfied
;**************************************************
  olr_threshold = 240 
  olr_thres  = dim_min_n(olr_region,(/1,2/))   ;min at each time step
 
  uwnd_threshold = 0
  uwnd_thres  = dim_min_n(u_region,(/1,2/))   ;min at each time step

  stretch_threshold = 0
  stretch_thres  = dim_max_n(stretch_region,(/1,2/))   ;max at each time step


  n_thres  = ind((olr_thres .le. olr_threshold) .and. (uwnd_thres .le. uwnd_threshold) .and. (stretch_thres .ge. stretch_threshold))    ;time indices
 
  if(ismissing(n_thres(0))) then
     print("++++++++++++++")
     print(" NO values exceed threshold="+threshold)
     print("++++++++++++++")
  exit
  end if
  
;***************************************************
;print
;***************************************************

  ymd       = cd_calendar(u&time, -2)            ;-- all dates
  print(n_thres+"  "+ymd(n_thres))

;***********************************************
; save to a netcdf file
;***********************************************
  diro = "./"
  filo = "olr.uwnd.stretch.threshold.nc"
  ptho = diro+filo
  system("/bin/rm -f "+ptho)
  ncdf = addfile(ptho,"c")

  fAtt = True
  fAtt@title         = "olr_uwnd_stretch >= threshold="+threshold 
  fAtt@source_file   = "daily winds from NCEP-NCAR Reanaysis"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc("date")
  fileattdef(ncdf,fAtt)            ; copy file attributes
  filedimdef(ncdf,"time",-1,True)

  ;ncdf->vr = vr(nt_thres,:,:) ;subset dates


end