; Write out HAMSR retrieval data in NetCDF format for bufr fortran program
; Also thins the data to the method the user provides. See below for more information
; Andrew Kren
; May 29, 2018
; June 13, 2018: Updated to provide second filter to thin based on time windows

begin

path_hamsr = "/scratch4/BMC/shout/Andrew.Kren/HAMSR/" ; location of HAMSR data

date = tostring(start_d) ; date you want to process, from wrapper ksh script

thin_method = "superob" ; "superob", "random" or "time_dist" distance thinning on the data. time_dist thinning uses a distance, time, and structure thinning method

;;;;;;;;;;; settings for random thinning below ;;;;;;;;;;;;;;;;
retain = 0.5  ; percentage of data to retain after random resampling, 0.8 equals 80% if using random thinning

;;;;;;;;;;; settings for time_dist thinning below ;;;;;;;;;;
temp_threshold = 0.5   ; temperature threshold (degC) to use to thin data if using time_dist thinning
shum_threshold = 0.5   ; specific humidity threshold (g/kg) to use to thin data if using time_dist thinning
pct_threshold = 5.0   ; percent threshold to use to thin data if using time_dist thinning
dist_retain = 2.0      ; value for minimum distance (km) to retain data if using time_dist thinning
time_retain = 0.4      ; value for minimum time distance (minutes) to retain data if using time_dist thinning
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; end of user's settings
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
print("Reading HAMSR data...")

fils = systemfunc("ls " + path_hamsr + "*" +date+"*.nc")  ; file paths
w = addfiles(fils,"r")

; read data
lat = w[:]->lat
lat := lat * 0.001
lon = w[:]->lon
lon := lon *0.001
time = w[:]->time
time@units = "seconds since 2000-01-01 00:00:00.0"
t = short2flt(w[:]->ham_airT)
q = short2flt(w[:]->ham_airQ)
rh = short2flt(w[:]->ham_airRH)
lev = short2flt(w[:]->ham_pres_levels)
lev := lev(0:24)
t@units = "K"
delete(w)
t@_FillValue = -999
q@_FillValue = -999
time@_FillValue = -999

; first filter out dates which are -nan (undefined) in the original dataset
good_data=ind(time.ne."-nan")
t := t(good_data,:)
lat := lat(good_data)
lon := lon(good_data)
time := time(good_data)
q := q(good_data,:)
rh := rh(good_data,:)
; calculate specific humidity
shum = mixhum_ptrh(conform(t,lev,1),t,rh,-2) ; g/kg

; profiles before thinning
num_before = dimsizes(time)

; convert temperature to degC after calculating shum since shum needs K unit
t = t - 273.15 ; convert to celsius
t@units = "degC"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; plot data spatially prior to thinning the data
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

outfile = "hamsr_spatial_map_before" 
wks = gsn_open_wks("png",outfile)
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
res = True
res@pmTickMarkDisplayMode  = "Always"           ; turn on fancy tickmarks
res@mpOutlineOn = True                      ; turn the map outline on
res@mpDataBaseVersion = "MediumRes"         ; Medium resolution
res@mpDataSetName = "Earth..4"
res@gsnMaximize          = True
res@tiMainString         = "HAMSR before thinning " + date

;---Zoom in on United States.
  res@mpMinLatF             = min(lat)-5.0
  res@mpMaxLatF             = max(lat)+5.0
  res@mpMinLonF             = min(lon)-5.0
  res@mpMaxLonF             = max(lon)+5.0
  res@mpCenterLonF          = floor((min(lon)-5.0 + min(lon)+5.0)/2.0)
res@mpFillOn = False
res@gsnRightString = ""
res@gsnLeftString = ""
res@mpGridAndLimbOn = True
res@mpGridLatSpacingF = 10.0
res@mpGridLonSpacingF = 10.0
res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries
res@mpUSStateLineThicknessF = 2.5
res@gsnFrame             = False               ; don't advance frame
res@gsnDraw              = False               ; don't draw plot

plot = gsn_csm_map(wks,res)


gres = True
gres@gsMarkerColor = "Red"   
gres@gsMarkerIndex = 1 ; dots
dum3 = gsn_add_polymarker(wks,plot,lon,lat,gres)


draw(plot)
frame(wks)

delete(wks)
delete(res)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Perform thinning based on user input
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;; random thinning ;;;;;;;;;;;;;;;;;;;;;;
if (thin_method.eq."random") then ; random subsampling for thinning

; generate indices for resampling without replacement to sub-sample the HAMSR data
ntime = dimsizes(time)
iwo = generate_sample_indices(ntime,0) ; indices from entire period
n_retain = toint(retain*ntime)
iwo_retain = iwo(0:n_retain-1)
t := t(iwo_retain,:)
lat := lat(iwo_retain)
lon := lon(iwo_retain)
time := time(iwo_retain)
q := q(iwo_retain,:)
rh := rh(iwo_retain,:)
shum := shum(iwo_retain,:)

end if
;;;;;;;;;;;;;; random thin end ;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;; superob thinning ;;;;;;;;;;;;;;;;;;;;;;
if (thin_method.eq."superob") then
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
; regrid the unstructured data to a triangular mesh grid

    r2d     = 180.0d/(atan(1)*4.0d)  ; Radian to Degree
    lonCell = lon;*r2d
    latCell = lat;*r2d

    Opt                 = True                     ; Regridding optioins

    Opt@SrcFileName     = "HAMSR_" + date + "out.nc"
    Opt@DstFileName     = "HAMSR_" + date + ".nc"
    Opt@WgtFileName     = "HAMSR_" + date + "_wgt.nc"
    Opt@ForceOverwrite  = True
    Opt@SrcESMF         = True

    Opt@SrcGridLat      = latCell                  ; source grid
    Opt@SrcGridLon      = lonCell
    ; determine min/max lat/lons for grid settings
    lon_e = where(lon.lt.0,lon+360.0,lon)

    Opt@DstGridType     = "0.125deg"                ; destination grid
    Opt@DstTitle        = "0.125 grid resolution"
    Opt@DstLLCorner     = (/todouble(min(lat)),todouble(min(lon_e))/)
    Opt@DstURCorner     = (/todouble(max(lat)),todouble(max(lon_e))/) 

    Opt@Debug           = True
    Opt@PrintTimings  = True

    t_regrid = ESMF_regrid(t(0,:),Opt)
    printVarSummary(t_regrid)
    exit

end if
;;;;;;;;;;;; end superob thinning ;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;; time_dist thinning ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (thin_method.eq."time_dist") then  ; time thinning of the data

; reorder arrays since time in original files is non-monotonic
;--- return permutation vector for sorting time
ip = dim_pqsort(t&time,1)

; Use ip vector to reorder arrays. This will also have the effect
; of  reordering the time coordinate array.
;
t := t(ip,:)   ; now is monontically increasing
time := time(ip)
lat := lat(ip)
lon := lon(ip)
q := q(ip,:)
rh := rh(ip,:)
shum := shum(ip,:)

t@_FillValue = -999
q@_FillValue = -999
time@_FillValue = -999
rh@_FillValue =	-999
shum@_FillValue = -999
lat@_FillValue = -999
lon@_FillValue = -999

; make arrays
t_new = new((/dimsizes(time),dimsizes(lev)/),"float",t@_FillValue)
shum_new = t_new
lat_new = new((/dimsizes(time)/),"float",t@_FillValue)
lon_new = lat_new
time_new = new((/dimsizes(time)/),"double",t@_FillValue)

do i=0,dimsizes(time)-1  

;;;;
if (i.eq.0) then ; if first profile then save and set reference profile 
ref_time = time(i)    ; reference time set to first data profile when first iterating through the data
ref_t = t(i,:)          ; reference temp.
ref_q = shum(i,:)       ; reference shum.
ref_lat = lat(i)        ; ref lat
ref_lon = lon(i)        ; ref lon

; save initial profile
t_new(i,:) = t(i,:)
shum_new(i,:) = shum(i,:)
time_new(i) = time(i)
lat_new(i) = lat(i)
lon_new(i) = lon(i)
end if
;;;;

; compare reference profile to next profile

; 1. compute distance first

gcdist = gc_latlon(ref_lat,ref_lon,lat(i),lon(i),2,4) ; km distance

; 2. compute time second

tdist = abs(ref_time - time(i)) / 60.0 ; seconds elapsed converted to minutes

; 3. check structure last, abs mean change and mean % change

t_diff = avg( abs(t(i,:) - ref_t) ) ; degC scalar
shum_diff = avg( abs(shum(i,:) - ref_q) ) ; g/kg scalar

t_pct_diff = 100.0 * ( avg( abs( (t(i,:) - ref_t)/ref_t ) ) ) ; % scalar
q_pct_diff = 100.0 * ( avg( abs( (shum(i,:) - ref_q)/ref_q ) ) ) ; % scalar

; check for new reference profile
if (gcdist.gt.dist_retain.and.tdist.gt.time_retain) then ; satisfies the time and dist thresholds

; check strcture last
if (t_diff.gt.temp_threshold.and.t_pct_diff.gt.pct_threshold .or. shum_diff.gt.shum_threshold.and.q_pct_diff.gt.pct_threshold) then

; set new reference profile if meets all criteria, otherwise will not get here and remain a _FillValue
ref_time = time(i)    ; reference time set to first data profile when first iterating through the data
ref_t = t(i,:)          ; reference temp.
ref_q = shum(i,:)	; reference shum.
ref_lat = lat(i)        ; ref lat
ref_lon = lon(i)

; save data to new arrays
t_new(i,:) = t(i,:)
shum_new(i,:) = shum(i,:)
time_new(i) = time(i)
lat_new(i) = lat(i)
lon_new(i) = lon(i)

end if ; end structure check

end if ; end time-dist check if statement

end do

iwo_retain = ind(.not.ismissing(t_new(:,0))) ; filter out missing data and remake arrays
t := t_new(iwo_retain,:)
lat := lat_new(iwo_retain)
lon := lon_new(iwo_retain)
time := time_new(iwo_retain)
shum := shum_new(iwo_retain,:)

delete(ip)
end if
;;;;;;;;;;;;;;;;;;;;;;;;; time_dist thin end ;;;;;;;;;;;;;;;;;;;;;;;;;;

; profiles after thinning
num_after = dimsizes(time)

print("Number of profiles retained: " + tofloat(num_after)/tofloat(num_before)*100.0 + "%")

; reorder arrays since time in original files is non-monotonic
;--- return permutation vector for sorting time
ip = dim_pqsort(t&time,1)

; Use ip vector to reorder arrays. This will also have the effect
; of  reordering the time coordinate array.
;
t := t(ip,:)   ; now is monontically increasing
time := time(ip)
lat := lat(ip)
lon := lon(ip)
shum := shum(ip,:)

shum@_FillValue = -999
copy_VarMeta(t,shum)
shum@comment = "HAMSR Specific Humidity"
shum@long_name = "HAMSR Specific Humidity"
shum@units = "g/kg"
lat@_FillValue = -999
lon@_FillValue = -999
t@_FillValue = -999
time@_FillValue = -999

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; plot data spatially

outfile = "hamsr_spatial_map_after"                                                                                                                               
wks = gsn_open_wks("png",outfile)
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
res = True
res@pmTickMarkDisplayMode  = "Always"           ; turn on fancy tickmarks
res@mpOutlineOn = True                      ; turn the map outline on
res@mpDataBaseVersion = "MediumRes"         ; Medium resolution
res@mpDataSetName = "Earth..4"
res@gsnMaximize          = True
res@tiMainString         = "HAMSR after thinning " + date

;---Zoom in on United States.
res@mpMinLatF             = min(lat)-5.0
  res@mpMaxLatF             = max(lat)+5.0
  res@mpMinLonF             = min(lon)-5.0
  res@mpMaxLonF             = max(lon)+5.0
  res@mpCenterLonF          = floor((min(lon)-5.0 + min(lon)+5.0)/2.0)
res@mpFillOn = False
res@gsnRightString = ""
res@gsnLeftString = ""
res@mpGridAndLimbOn = True
res@mpGridLatSpacingF = 10.0
res@mpGridLonSpacingF = 10.0
res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries
res@mpUSStateLineThicknessF = 2.5
res@gsnFrame             = False               ; don't advance frame
res@gsnDraw              = False               ; don't draw plot

plot = gsn_csm_map(wks,res)



gres = True
gres@gsMarkerColor = "Blue"
gres@gsMarkerIndex = 1 ; dots
dum3 = gsn_add_polymarker(wks,plot,lon,lat,gres)


draw(plot)
frame(wks)

delete(wks)
delete(res)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;
; Second step of data thinning process
; Go through the assimilation time windows 
; Filter out pts close in proximity
; Keep obs closest to the center of the time window
;;;;;;;;;;;;;;;;;

;;; time-window thinning

; get dates
time_hamsr = cd_calendar(time,0)
year = tointeger(time_hamsr(:,0))
month = tointeger(time_hamsr(:,1))
day = tointeger(time_hamsr(:,2))
hour = tointeger(time_hamsr(:,3))
minute = tointeger(time_hamsr(:,4))
sec = tointeger(time_hamsr(:,5))

; make new arrays for all four assimilation cycles
t_save = new((/dimsizes(time),dimsizes(lev)/),"float",-999)
shum_save = t_save
time_save = new((/dimsizes(time)/),"double",-999)
lat_save = new((/dimsizes(time)/),"float",-999)
lon_save = lat_save

t_temp = t         ; temp arrays for manipulating prior thinned data
shum_temp = shum
lat_temp = lat
lon_temp = lon
time_temp = time

z = 0 ; counter to save the data after thinning, we increment this after each profile and cycle
      ; does not initialize after each cycle

print("Thinning for first 00z timeframe")
;;; 00z one

num_window_profiles = ind(hour.ge.21.or.hour.lt.0) ; 00z

if (.not.all(ismissing(num_window_profiles))) then ; check if missing since could be missing utc times

do i=0,dimsizes(num_window_profiles)-1 ; loop over this window's profiles

if (.not.ismissing(lat_temp(num_window_profiles(i)))) then

; check distance from this ob to all other obs, if within some distance, keep profile closest to center time window

distances = gc_latlon(lat_temp(num_window_profiles(i)),lon_temp(num_window_profiles(i)),lat_temp(num_window_profiles),lon_temp(num_window_profiles),2,4) ; returns km distance

; check distances from this point to all others  ; dist_retain
dist_check = ind(distances.lt.dist_retain.and.distances.ne.0)

if (.not.all(ismissing(dist_check))) then ; then points are nearby

; keep profile closest to central time window

; compute time difference
time_00z = cd_inv_calendar(year(0),month(0),day(0)+1,0,0,0,"seconds since 2000-01-01 00:00:00.0",0)

time_diff = abs(time_00z-time(num_window_profiles(dist_check))) ; abs difference

; select profile closest to central time window, discard the rest
print("selecting profile closest to the central time window...")

t_save(z,:) = t_temp(num_window_profiles(dist_check(minind(time_diff))),:)
shum_save(z,:) = shum_temp(num_window_profiles(dist_check(minind(time_diff))),:)
lat_save(z) = lat_temp(num_window_profiles(dist_check(minind(time_diff))))
lon_save(z) = lon_temp(num_window_profiles(dist_check(minind(time_diff))))
time_save(z) = time_temp(num_window_profiles(dist_check(minind(time_diff))))

; remove data we already went through
t_temp(num_window_profiles(dist_check),:) = t@_FillValue
shum_temp(num_window_profiles(dist_check),:) = shum@_FillValue
lat_temp(num_window_profiles(dist_check)) = lat@_FillValue
lon_temp(num_window_profiles(dist_check)) = lon@_FillValue
time_temp(num_window_profiles(dist_check)) = time@_FillValue

delete(time_diff)
delete(dist_check)
z = z + 1  ; increment z

else ; keep this profile

print("keeping this profile...")
t_save(z,:) = t_temp(num_window_profiles(i),:)
shum_save(z,:) = shum_temp(num_window_profiles(i),:)
lat_save(z) = lat_temp(num_window_profiles(i))
lon_save(z) = lon_temp(num_window_profiles(i))
time_save(z) = time_temp(num_window_profiles(i))

z = z + 1  ; increment z
end if

delete(distances)

end if

end do

end if


print("Thinning for second 00z timeframe")
delete(num_window_profiles)
;;; 00z second

num_window_profiles = ind(hour.ge.0.and.hour.lt.3) ; 00z

if (.not.all(ismissing(num_window_profiles))) then ; check if missing since could be missing utc times

do i=0,dimsizes(num_window_profiles)-1 ; loop over this window's profiles

if (.not.ismissing(lat_temp(num_window_profiles(i)))) then

; check distance from this ob to all other obs, if within some distance, keep profile closest to center time window

distances = gc_latlon(lat_temp(num_window_profiles(i)),lon_temp(num_window_profiles(i)),lat_temp(num_window_profiles),lon_temp(num_window_profiles),2,4) ; returns km distance

; check distances from this point to all others  ; dist_retain
dist_check = ind(distances.lt.dist_retain.and.distances.ne.0)

if (.not.all(ismissing(dist_check))) then ; then points are nearby

; keep profile closest to central time window

; compute time difference
time_00z = cd_inv_calendar(year(0),month(0),day(0),0,0,0,"seconds since 2000-01-01 00:00:00.0",0)

time_diff = abs(time_00z-time(num_window_profiles(dist_check))) ; abs difference

; select profile closest to central time window, discard the rest
print("selecting profile closest to the central time window...")

t_save(z,:) = t_temp(num_window_profiles(dist_check(minind(time_diff))),:)
shum_save(z,:) = shum_temp(num_window_profiles(dist_check(minind(time_diff))),:)
lat_save(z) = lat_temp(num_window_profiles(dist_check(minind(time_diff))))
lon_save(z) = lon_temp(num_window_profiles(dist_check(minind(time_diff))))
time_save(z) = time_temp(num_window_profiles(dist_check(minind(time_diff))))

; remove data we already went through
t_temp(num_window_profiles(dist_check),:) = t@_FillValue
shum_temp(num_window_profiles(dist_check),:) = shum@_FillValue
lat_temp(num_window_profiles(dist_check)) = lat@_FillValue
lon_temp(num_window_profiles(dist_check)) = lon@_FillValue
time_temp(num_window_profiles(dist_check)) = time@_FillValue

delete(time_diff)
delete(dist_check)
z = z + 1  ; increment z

else ; keep this profile

print("keeping this profile...")
t_save(z,:) = t_temp(num_window_profiles(i),:)
shum_save(z,:) = shum_temp(num_window_profiles(i),:)
lat_save(z) = lat_temp(num_window_profiles(i))
lon_save(z) = lon_temp(num_window_profiles(i))
time_save(z) = time_temp(num_window_profiles(i))

z = z + 1  ; increment z
end if

delete(distances)

end if

end do

end if

;;; 06z
print("Thinning for 06z timeframe")
delete(num_window_profiles)

num_window_profiles = ind(hour.ge.3.and.hour.lt.9) ; 06z

if (.not.all(ismissing(num_window_profiles))) then ; check if missing since could be missing utc times

do i=0,dimsizes(num_window_profiles)-1 ; loop over this window's profiles

if (.not.ismissing(lat_temp(num_window_profiles(i)))) then

; check distance from this ob to all other obs, if within some distance, keep profile closest to center time window

distances = gc_latlon(lat_temp(num_window_profiles(i)),lon_temp(num_window_profiles(i)),lat_temp(num_window_profiles),lon_temp(num_window_profiles),2,4) ; returns km distance

; check distances from this point to all others  ; dist_retain
dist_check = ind(distances.lt.dist_retain.and.distances.ne.0)

if (.not.all(ismissing(dist_check))) then ; then points are nearby

; keep profile closest to central time window

; compute time difference
time_06z = cd_inv_calendar(year(0),month(0),day(0),6,0,0,"seconds since 2000-01-01 00:00:00.0",0)

time_diff = abs(time_06z-time(num_window_profiles(dist_check))) ; abs difference

; select profile closest to central time window, discard the rest
print("selecting profile closest to the central time window...")

t_save(z,:) = t_temp(num_window_profiles(dist_check(minind(time_diff))),:)
shum_save(z,:) = shum_temp(num_window_profiles(dist_check(minind(time_diff))),:)
lat_save(z) = lat_temp(num_window_profiles(dist_check(minind(time_diff))))
lon_save(z) = lon_temp(num_window_profiles(dist_check(minind(time_diff))))
time_save(z) = time_temp(num_window_profiles(dist_check(minind(time_diff))))

; remove data we already went through
t_temp(num_window_profiles(dist_check),:) = t@_FillValue
shum_temp(num_window_profiles(dist_check),:) = shum@_FillValue
lat_temp(num_window_profiles(dist_check)) = lat@_FillValue
lon_temp(num_window_profiles(dist_check)) = lon@_FillValue
time_temp(num_window_profiles(dist_check)) = time@_FillValue

delete(time_diff)
delete(dist_check)
z = z + 1  ; increment z

else ; keep this profile

print("keeping this profile...")

t_save(z,:) = t_temp(num_window_profiles(i),:)
shum_save(z,:) = shum_temp(num_window_profiles(i),:)
lat_save(z) = lat_temp(num_window_profiles(i))
lon_save(z) = lon_temp(num_window_profiles(i))
time_save(z) = time_temp(num_window_profiles(i))

z = z + 1  ; increment z
end if

delete(distances)

end if

end do

end if

;;; 12z
print("Thinning for 12z timeframe")
delete(num_window_profiles)

num_window_profiles = ind(hour.ge.9.and.hour.lt.15) ; 12z

if (.not.all(ismissing(num_window_profiles))) then ; check if missing since could be missing utc times

do i=0,dimsizes(num_window_profiles)-1 ; loop over this window's profiles

if (.not.ismissing(lat_temp(num_window_profiles(i)))) then

; check distance from this ob to all other obs, if within some distance, keep profile closest to center time window

distances = gc_latlon(lat_temp(num_window_profiles(i)),lon_temp(num_window_profiles(i)),lat_temp(num_window_profiles),lon_temp(num_window_profiles),2,4) ; returns km distance

; check distances from this point to all others  ; dist_retain
dist_check = ind(distances.lt.dist_retain.and.distances.ne.0)

if (.not.all(ismissing(dist_check))) then ; then points are nearby

; keep profile closest to central time window

; compute time difference
time_12z = cd_inv_calendar(year(0),month(0),day(0),12,0,0,"seconds since 2000-01-01 00:00:00.0",0)

time_diff = abs(time_12z-time(num_window_profiles(dist_check))) ; abs difference

; select profile closest to central time window, discard the rest
print("selecting profile closest to the central time window...")

t_save(z,:) = t_temp(num_window_profiles(dist_check(minind(time_diff))),:)
shum_save(z,:) = shum_temp(num_window_profiles(dist_check(minind(time_diff))),:)
lat_save(z) = lat_temp(num_window_profiles(dist_check(minind(time_diff))))
lon_save(z) = lon_temp(num_window_profiles(dist_check(minind(time_diff))))
time_save(z) = time_temp(num_window_profiles(dist_check(minind(time_diff))))

; remove data we already went through
t_temp(num_window_profiles(dist_check),:) = t@_FillValue
shum_temp(num_window_profiles(dist_check),:) = shum@_FillValue
lat_temp(num_window_profiles(dist_check)) = lat@_FillValue
lon_temp(num_window_profiles(dist_check)) = lon@_FillValue
time_temp(num_window_profiles(dist_check)) = time@_FillValue

delete(time_diff)
delete(dist_check)
z = z + 1  ; increment z

else ; keep this profile

print("keeping this profile...")
t_save(z,:) = t_temp(num_window_profiles(i),:)
shum_save(z,:) = shum_temp(num_window_profiles(i),:)
lat_save(z) = lat_temp(num_window_profiles(i))
lon_save(z) = lon_temp(num_window_profiles(i))
time_save(z) = time_temp(num_window_profiles(i))

z = z + 1  ; increment z
end if

delete(distances)

end if

end do

end if

;;; 18z
print("Thinning for 18z timeframe")
delete(num_window_profiles)

num_window_profiles = ind(hour.ge.15.and.hour.lt.21) ; 18z

if (.not.all(ismissing(num_window_profiles))) then ; check if missing since could be missing utc times

do i=0,dimsizes(num_window_profiles)-1 ; loop over this window's profiles

if (.not.ismissing(lat_temp(num_window_profiles(i)))) then

; check distance from this ob to all other obs, if within some distance, keep profile closest to center time window

distances = gc_latlon(lat_temp(num_window_profiles(i)),lon_temp(num_window_profiles(i)),lat_temp(num_window_profiles),lon_temp(num_window_profiles),2,4) ; returns km distance

; check distances from this point to all others  ; dist_retain
dist_check = ind(distances.lt.dist_retain.and.distances.ne.0)

if (.not.all(ismissing(dist_check))) then ; then points are nearby

; keep profile closest to central time window

; compute time difference
time_18z = cd_inv_calendar(year(0),month(0),day(0),18,0,0,"seconds since 2000-01-01 00:00:00.0",0)

time_diff = abs(time_18z-time(num_window_profiles(dist_check))) ; abs difference

; select profile closest to central time window, discard the rest
print("selecting profile closest to the central time window...")

t_save(z,:) = t_temp(num_window_profiles(dist_check(minind(time_diff))),:)
shum_save(z,:) = shum_temp(num_window_profiles(dist_check(minind(time_diff))),:)
lat_save(z) = lat_temp(num_window_profiles(dist_check(minind(time_diff))))
lon_save(z) = lon_temp(num_window_profiles(dist_check(minind(time_diff))))
time_save(z) = time_temp(num_window_profiles(dist_check(minind(time_diff))))

; remove data we already went through
t_temp(num_window_profiles(dist_check),:) = t@_FillValue
shum_temp(num_window_profiles(dist_check),:) = shum@_FillValue
lat_temp(num_window_profiles(dist_check)) = lat@_FillValue
lon_temp(num_window_profiles(dist_check)) = lon@_FillValue
time_temp(num_window_profiles(dist_check)) = time@_FillValue

delete(time_diff)
delete(dist_check)
z = z + 1  ; increment z

else ; keep this profile

print("keeping this profile...")
t_save(z,:) = t_temp(num_window_profiles(i),:)
shum_save(z,:) = shum_temp(num_window_profiles(i),:)
lat_save(z) = lat_temp(num_window_profiles(i))
lon_save(z) = lon_temp(num_window_profiles(i))
time_save(z) = time_temp(num_window_profiles(i))

z = z + 1  ; increment z
end if

delete(distances)

end if

end do

end if

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; plot data spatially

outfile = "hamsr_spatial_map_after_window"
wks = gsn_open_wks("png",outfile)
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
res = True
res@pmTickMarkDisplayMode  = "Always"           ; turn on fancy tickmarks
res@mpOutlineOn = True                      ; turn the map outline on
res@mpDataBaseVersion = "MediumRes"         ; Medium resolution
res@mpDataSetName = "Earth..4"
res@gsnMaximize          = True
res@tiMainString         = "HAMSR after second thinning " + date

;---Zoom in on United States.
res@mpMinLatF             = min(lat)-5.0
  res@mpMaxLatF             = max(lat)+5.0
  res@mpMinLonF             = min(lon)-5.0
  res@mpMaxLonF             = max(lon)+5.0
  res@mpCenterLonF          = floor((min(lon)-5.0 + min(lon)+5.0)/2.0)
res@mpFillOn = False
res@gsnRightString = ""
res@gsnLeftString = ""
res@mpGridAndLimbOn = True
res@mpGridLatSpacingF = 10.0
res@mpGridLonSpacingF = 10.0
res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries
res@mpUSStateLineThicknessF = 2.5
res@gsnFrame             = False               ; don't advance frame
res@gsnDraw              = False               ; don't draw plot

plot = gsn_csm_map(wks,res)



gres = True
gres@gsMarkerColor = "Green4"
gres@gsMarkerIndex = 1 ; dots
dum3 = gsn_add_polymarker(wks,plot,lon_save,lat_save,gres)


draw(plot)
frame(wks)

delete(wks)
delete(res)

;;; end time-window thinning

; reorder arrays after time-window thinning
delete(iwo_retain)
iwo_retain = ind(.not.ismissing(t_save(:,0))) ; filter out missing data and remake arrays
t := t_save(iwo_retain,:)
lat := lat_save(iwo_retain)
lon := lon_save(iwo_retain)
time := time_save(iwo_retain)
shum := shum_save(iwo_retain,:)


; reorder arrays since time in original files is non-monotonic
;--- return permutation vector for sorting time
delete(ip)
ip = dim_pqsort(t&time,1)

; Use ip vector to reorder arrays. This will also have the effect
; of  reordering the time coordinate array.
;
t := t(ip,:)   ; now is monontically increasing
time := time(ip)
lat := lat(ip)
lon := lon(ip)
shum := shum(ip,:)

; calculate dates
time_hamsr := cd_calendar(time,0)
year := tostring(tointeger(time_hamsr(:,0)))
month := tostring(tointeger(time_hamsr(:,1)))
day := tostring(tointeger(time_hamsr(:,2)))
hour := tostring(tointeger(time_hamsr(:,3)))
minute := tostring(tointeger(time_hamsr(:,4)))
sec := tostring(tointeger(time_hamsr(:,5)))

; fix time to add missing zeros to certain dates and times
month = where(tointeger(month).lt.10,"0"+month,month)
day = where(tointeger(day).lt.10,"0"+day,day)
hour = where(tointeger(hour).lt.10,"0"+hour,hour)
minute = where(tointeger(minute).lt.10,"0"+minute,minute)
sec = where(tointeger(sec).lt.10,"0"+sec,sec)

; combine dates and lat/lon together into one
dates = year ; create array of strings
lat_lon = year
do i=0,dimsizes(year)-1

strs1 = (/year(i),month(i),day(i),hour(i),minute(i),sec(i)/)
temp_date = str_join(strs1,"    ")
dates(i) = temp_date

delete(temp_date)
delete(strs1)

strs1 = (/tostring(lat(i)),tostring(lon(i))/)
temp_date = str_join(strs1,"       ")

lat_lon(i) = temp_date
delete(temp_date)
delete(strs1)
end do

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; write data to text file for fortran later...
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
print("Writing HAMSR retrievals to file...")

output_file	  = "hamsr_tq_values.txt"
system("rm -rf "+ output_file)   ; remove any pre-existing file
opt = True
opt@fout = output_file
opt@row = False ; don't print row numbers

; conform time and lat to t and shum

dates_c = conform_dims(dimsizes(t),dates,0)
lat_lon_c = conform_dims(dimsizes(t),lat_lon,0)

; lat,lon,temp,and shum data
z = 0 ; for levels
do i=219,219;dimsizes(year)-1 
;print("Working on date..." + str_squeeze(dates(i)))
; plot pressure levels only once
if (z.eq.0) then
do k=0,24 ; loop over levels
if (k.eq.0) then
write_table(output_file,"w",[/(/tostring(tointeger(lev(k)))/)/],"%s")
else
write_table(output_file,"a",[/(/tostring(tointeger(lev(k)))/)/],"%s")
end if
end do
z = 1  	
end if

; dates and lat/lon
write_table(output_file,"a",[/(/dates(i)/)/],"%s")
write_table(output_file,"a",[/(/lat_lon(i)/)/],"%s")

; t and q
alist = [/t(i,0),t(i,1),t(i,2),t(i,3),t(i,4),t(i,5),t(i,6),t(i,7),t(i,8),t(i,9),t(i,10),t(i,11),t(i,12),t(i,13),t(i,14),t(i,15),t(i,16),t(i,17),t(i,18),t(i,19),t(i,20),t(i,21),t(i,22),t(i,23),t(i,24)/]
write_table(output_file,"a",alist,"%f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f")

alist = [/shum(i,0),shum(i,1),shum(i,2),shum(i,3),shum(i,4),shum(i,5),shum(i,6),shum(i,7),shum(i,8),shum(i,9),shum(i,10),shum(i,11),shum(i,12),shum(i,13),shum(i,14),shum(i,15),shum(i,16),shum(i,17),shum(i,18),shum(i,19),shum(i,20),shum(i,21),shum(i,22),shum(i,23),shum(i,24)/]
write_table(output_file,"a",alist,"%f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f   %f")

end do



print("Program ended successfully")
end