pro tdg2_nc_python,year,month
; ****************************************************************
; This procedure will take Kwaj gridded data with 1km H res and 1 km
; height and will assign a label corresponding to the 0.1 deg box number
; from a 0.1 deg grid.  Label will range from 0 to 528 depending on 
; resulting X,Y position.
;
; Write netcdf file with all relevant info.
; Last modified: 1/19/2021 for inclusion of 0.0 rates if no rainy pixels.
; Modified May 2022 for this version to read python-gridded files at 1km level.
; 
; David Marks
; SSAI @ WFF
; ********************************************************************


  year  = string(year,format='(I4.4)')
  month = string(month,format='(I2.2)')

; Directory structures
  in_dir = '/gpmraid/gpmarchive/Radar/KPOL/Gridded/' + year + '/' + month + '/temp/'
  out_dir = '/gpmraid/gpmarchive/Radar/KPOL/Gridded/IMERG_grid/' + year + '/'
  spawn,'mkdir -p ' + out_dir

; Find data files -- these are the 1km H-res gridded files
  wc = in_dir + '*nc.gz'
  files = file_search(wc,COUNT=nf)
  if (nf eq 0) then begin
        print,'No NetCDF files found in ' + in_dir
        stop
  endif
;
;
; ** Files were found....read nc into IDL structure
  for ifiles=0L,nf-1 do begin
        cf = files[ifiles]
        flag = uncomp_file(cf,file)
        if(flag ne 'OK') then begin
            print,flag
            goto,next_file
        endif
        print,'<-- ' + file
	ncdf_read,orig_grid,filename=file,/all

; Set up our 0.1 deg grid centered on Kwaj radar location or IMERG 0.1 deg center
; get_grids,dx,dy,dz,nx,ny,nz,center_lat,center_lon,xgrid,ygrid,zgrid
; 	get_deg_grids,0.1,0.1,1,23,23,1,8.718,167.732,xgrid,ygrid,zgrid; KPOL radar location
 	get_deg_grids,0.1,0.1,1,23,23,1,8.75,167.75,xgrid,ygrid,zgrid; KWAJ -- modified center to match IMERG 0.1 deg center per J. Wang
;	get_deg_grids,0.1,0.1,1,23,23,1,38.25,-75.35,xgrid,ygrid,zgrid; NPOL-Newark -- modified center to match IMERG 0.1 deg center per J. Wang

; ** Read lon/lat for each 1km data point and find closest xgrid / ygrid position
; ** This is set for a 301 x 301 grid
; ** Initialize arrays for each file
	diff_lon = fltarr(23)
	diff_lat = fltarr(23)
	point_label = intarr(301,301)
	point_label[*,*] = -99; initalize
	rc_1d = [ ]
	pt_1d = [ ]
        total_points = intarr(23,23)
        rain_points = intarr(23,23)
        rc_uncond_mean = fltarr(23,23); unconditional
        rc_cond_mean = fltarr(23,23); conditional

	for iposn = 0,300 do begin
	    for jposn = 0,300 do begin
;		my_lon = orig_grid.lon0[iposn,jposn]; for radx gridded file
;		my_lat = orig_grid.lat0[iposn,jposn]
		my_lon = orig_grid.POINT_LONGITUDE[iposn,jposn,1]; for pyart gridded file
		my_lat = orig_grid.POINT_LATITUDE[iposn,jposn,1]
		diff_lon = abs(xgrid-my_lon)
		diff_lat = abs(ygrid-my_lat)
		if min(diff_lon) gt 0.05 or min(diff_lat gt 0.05) then begin; outside of our grid
		    point_label[iposn,jposn] = -99
		    goto,next_point
		endif
		xposn = where(diff_lon eq min(diff_lon))
		xposn = xposn[0]
		yposn = where(diff_lat eq min(diff_lat))
		yposn = yposn[0]
; ** Compute the 0.1 deg box label for this point
		point_label[iposn,jposn] = xposn + (yposn * 23)
; ** Put the point label and rc data for that point in 1d arrays for easier indexing
		pt_1d = [pt_1d,point_label[iposn,jposn]]
		rc_1d = [rc_1d,orig_grid.rc[iposn,jposn,1]]
next_point:
	    endfor
	endfor

; ** 0.1 deg box labels are complete for all points
; ** Compute rr averages (conditional and uncond.) and place in arrays.
; ** There are 529 0.1 deg boxes in our 23 x 23 array.


  	for ibox = 0,528 do begin
	    q = where(pt_1d eq ibox,cq)
	    total_points[ibox] = cq; # total points in this box
	    box_data = rc_1d(q)
	    valid = where(box_data gt 0.0,cvalid)
	    rain_points[ibox] = cvalid; # rainy points in this box
	    if (cvalid eq 0) then begin
		rc_uncond_mean[ibox] = 0.0; no rain points
		rc_cond_mean[ibox] = 0.0
	    endif else begin; we have rates > 0.0
		my_valid = box_data(valid); gt 0.0 rates
		rc_cond_mean[ibox] = mean(my_valid)

		uncond = where(box_data lt 0.0,c_uncond)
   		if (c_uncond gt 0) then box_data(uncond) = 0.0
		rc_uncond_mean[ibox] = mean(box_data)
	    endelse
	endfor; 0.1 deg boxes

; ** parse filename and write data 
  a = strsplit(file,'_',/EXTRACT)
  b = strsplit(a[3],'.',/EXTRACT)
  time = b[0]
  out_file = out_dir + a[0] + '_' + a[1] + a[2] + '_' + time + '_IMERG_GRID.nc'
  print
  print,'    -->' + out_file
  print

; Set up file and handler for netcdf output
  n_lon = n_elements(xgrid); same for lat ygrid
  n_lat = n_elements(ygrid);
  n_ht = n_elements(zgrid)

  fid=ncdf_create(out_file,/CLOBBER)
  d=indgen(3); number of dimensions

; Set Diminsions
  d[0]=ncdf_dimdef(fid,'lon_elements',n_lon)
  d[1]=ncdf_dimdef(fid,'lat_elements',n_lat)
  d[2]=ncdf_dimdef(fid,'ht_elements',n_ht)

; Define variables
  var_id=ncdf_vardef(fid,'Longitude',d[0],/float)
  var_id=ncdf_vardef(fid,'Latitude',d[1],/float)
  var_id=ncdf_vardef(fid,'Height',d[2],/float)
  var_id=ncdf_vardef(fid,'Total_points',[d[0],d[1]],/long)
  var_id=ncdf_vardef(fid,'Rain_points',[d[0],d[1]],/long)
  var_id=ncdf_vardef(fid,'RC_rate_uncond',[d[0],d[1]],/float)
  var_id=ncdf_vardef(fid,'RC_rate_cond',[d[0],d[1]],/float)
; var_id=ncdf_vardef(fid,'_FillValue',1,/float)

; Change modes to data mode
  ncdf_control,fid,/ENDEF

; Write data
  ncdf_varput,fid,'Longitude',xgrid
  ncdf_varput,fid,'Latitude',ygrid
  ncdf_varput,fid,'Height',zgrid
  ncdf_varput,fid,'Total_points',total_points
  ncdf_varput,fid,'Rain_points',rain_points
  ncdf_varput,fid,'RC_rate_uncond',rc_uncond_mean
  ncdf_varput,fid,'RC_rate_cond',rc_cond_mean

; Close file and release handler
  ncdf_close,fid
  print,' > Writing complete....Closing...'
  print
	    
	
next_file:
  spawn,'rm ' + './'+file; remove file copy
  endfor

  stop
  end