; Program to read the diagnostic radiance netcdf files for user's specified 
; satellites and time periods
; Plots all satellite passes and writes data to text files

; Andrew Kren
; May 18, 2016

begin

;; Insert the satellite passes you want to plot, i.e. edit the array
;; Also edit the dates you want to loop over

; code names for satellites, used for reading diagnostic files
sat_array = (/"cris_npp","atms_npp","amsua_aqua","airs_aqua","ssmis_f16","amsua_n15","amsub_n17","amsua_n18","amsua_n19","hirs3_n17","hirs4_n19","mhs_n18","mhs_n19","amsua_metop-a","sndrd4_g13","sndrd3_g13","sndrd2_g13","sndrd1_g13","seviri_m09","mhs_metop-a","iasi_metop-a","hirs4_metop-a"/)

; user friendly names of satellites for plotting routines
sat_names = (/"Cris NPP","Atms NPP","Amsua Aqua","Airs Aqua","DMSP F16","NOAA-15","NOAA-17","NOAA-18","NOAA-19","NOAA-17-hirs3","NOAA-19-hirs4","NOAA-18-mhs","NOAA-19-mhs","Amsua_metop-a","sndrd4_g13","sndrd3_g13","sndrd2_g13","sndrd1_g13","seviri_m09","mhs_metop-a","iasi_metop-a","hirs4_metop-a"/)

; dates to loop over for 3 OSSE cases
dates_array = (/"2006013000","2006013006","2006013012","2006013018","2006013100","2006013106","2006013112","2006022400","2006022406","2006022412","2006022418","2006022500"/)

; regional domain grid box to count number of obs per 6hr cycle
lon_w = 130 ; west end of grid box to include
lon_e = 235 ; east end of grid box
lat_s = 38  ; south end of grid box
lat_n = 80  ; north end of grid box

; flight domain grid box to count number of obs per 6hr cycle
lon_w_f = 180 ; west end of grid box to include
lon_e_f = 211 ; east end of grid box
lat_s_f = 37  ; south end	of grid	box
lat_n_f = 60  ; north end	of grid	box

; satellite dot size
sat_dot_size = 0.005

; path for diagnostic netcdf files...
path_files = "/scratch4/BMC/shout/Andrew.Kren/rad_sat/diag_radiance/"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; end of user's settings
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

sum_all_reg = new((/dimsizes(dates_array)/),"integer")
sum_all_fl = new((/dimsizes(dates_array)/),"integer")
sum_all_reg@_FillValue = 0
sum_all_fl@_FillValue = 0

sum_sat_reg = new((/dimsizes(dates_array),dimsizes(sat_names)/),"integer")
sum_sat_fl = new((/dimsizes(dates_array),dimsizes(sat_names)/),"integer")
sum_sat_reg@_FillValue = 0
sum_sat_fl@_FillValue = 0


; loop over dates first
do i=0,dimsizes(dates_array)-1  ; begin loop over dates
print("Current date: " + dates_array(i))

count_lat_reg = 0 ; initialize to zero at beginning of each time step
count_lon_reg = 0
count_total_lon_reg = 0 ; sum
count_total_lat_reg = 0

count_lat_fl = 0 ; initialize to zero at beginning of each time step
count_lon_fl = 0
count_total_lon_fl = 0 ; sum
count_total_lat_fl = 0

; loop over satellites next
do j=0,dimsizes(sat_array)-1 ; begin loop over satellites

; read file
filename = path_files + "diags_" + sat_array(j) + "_" + dates_array(i) + ".nc"

; check if file exists since not some satellites do not exist for every 6h timestep
if (fileexists(filename)) then
print("Reading file " + filename + "...")
w = addfile(filename,"r")
lat = w->lat
lon = w->lon

; sum up this individual satellite for this time period...
lat_x =	lat
lon_x =	lon
;qsort(lat_x)
;qsort(lon_x)
lat_x@_FillValue = -999
lon_x@_FillValue = -999

do y=0,dimsizes(lat_x)-1

 if (lat_x(y).ge.lat_s .and. lat_x(y).le.lat_n) then
  if (lon_x(y).ge.lon_w .and. lon_x(y).le.lon_e) then
   lat_x(y) = lat_x(y)
   lon_x(y) = lon_x(y)
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if

end do

count_lat_reg = num(.not.ismissing(lat_x)) ; number of satellite data for this particular satellite at this time step
count_lon_reg = num(.not.ismissing(lon_x))

count_total_lat_reg = count_lat_reg + count_total_lat_reg
count_total_lon_reg = count_lon_reg + count_total_lon_reg

lat_x = lat
lon_x = lon

do y=0,dimsizes(lat_x)-1

 if (lat_x(y).ge.lat_s_f .and. lat_x(y).le.lat_n_f) then
  if (lon_x(y).ge.lon_w_f .and. lon_x(y).le.lon_e_f) then
   lat_x(y) = lat_x(y)
   lon_x(y) = lon_x(y)
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if

end do

count_lat_fl = num(.not.ismissing(lat_x)) ; number of satellite data for this particular satellite at this time step
count_lon_fl = num(.not.ismissing(lon_x))

count_total_lat_fl = count_lat_fl + count_total_lat_fl
count_total_lon_fl = count_lon_fl + count_total_lon_fl

; delete variables
delete([/lat,lon,w,filename,lat_x,lon_x/])

else
print("No data for this date: " + dates_array(i) + " and satellite: " + sat_array(j)) 
end if

end do ; end loop over satellites

; sum up the total number of data points from all satellites at this time period...

sum_lat_reg = count_total_lat_reg ; total at this time period over the specified box region of ETS, specified in beginning of program
sum_lon_reg = count_total_lon_reg
sum_lat_fl = count_total_lat_fl ; total at this time period over the specified box region of ETS, specified in beginning of program
sum_lon_fl = count_total_lon_fl

print("Total regional satellite observations at: " + dates_array(i) + ": " + sum_lat_reg)
print("Total flight satellite observations at: " + dates_array(i) + ": " + sum_lat_fl)

; storing all satellites in regional and flight box domains
sum_all_reg(i) = sum_lat_reg
sum_all_fl(i) = sum_lat_fl


;;;;;;;;;; make plots of individual satellites ;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; loop over satellites next
do j=0,dimsizes(sat_array)-1 ; begin loop over satellites

count_lat_reg_sub = 0 ; initialize to zero at beginning of each time step
count_lon_reg_sub = 0
count_lat_fl_sub = 0 ; initialize to zero at beginning of each time step
count_lon_fl_sub = 0

; read file
filename = path_files + "diags_" + sat_array(j) + "_" + dates_array(i) + ".nc"

; check if file exists, if so make plot
if (fileexists(filename)) then
;print("Reading file " + filename + "...")
w = addfile(filename,"r")
lat = w->lat
lon = w->lon

; plot the map with satellite data
print("Making individual satellite pass plot " + sat_array(j) + " for current date:" + dates_array(i))
outfile = "satellite_radiance_" + sat_array(j) + "_" + dates_array(i)     ; output file

wks = gsn_open_wks("pdf",outfile)

;
; Set up some map resources.
;
  mpres              = True
  mpres@gsnMaximize  = True             ; Maximize plot in frame.
  mpres@gsnFrame     = False            ; Don't advance the frame.
  mpres@gsnDraw      = False
  mpres@mpMinLatF    = 0
  mpres@mpMaxLatF    = 90
  mpres@mpMinLonF    = 110
  mpres@mpMaxLonF    = 290
  mpres@mpCenterLonF = 220
  mpres@mpGridAndLimbOn = True
  mpres@mpGridLatSpacingF = 10.0
  mpres@mpGridLonSpacingF = 10.0
  mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates"
  mpres@mpGeophysicalLineThicknessF = 1.5 ; thickness of boundaries
  mpres@mpUSStateLineThicknessF = 1.5
  mpres@mpFillColor = (/"gray"/)
  mpres@tiMainString = "Satellite Data for: " + sat_names(j) + " " + dates_array(i)
  
  map = gsn_csm_map(wks,mpres)

  gsres               = True
  gsres@gsMarkerIndex = 16          ; Use filled dots for markers.

  gres = True
 gres@gsLineColor = "Blue"
 gres@gsLineThicknessF = 5.0
 lat2x = (/ lat_s_f, lat_s_f, lat_n_f, lat_n_f, lat_s_f/)
 lon2x = (/lon_w_f, lon_e_f, lon_e_f, lon_w_f, lon_w_f/)

 dum2 = gsn_add_polyline(wks,map,lon2x,lat2x,gres)


 gres = True
 gres@gsLineColor = "Red"
 gres@gsLineThicknessF = 5.0
 lat1x = (/ lat_s, lat_s, lat_n, lat_n, lat_s/)
 lon1x = (/lon_w, lon_e, lon_e, lon_w, lon_w/)

 dum = gsn_add_polyline(wks,map,lon1x,lat1x,gres)


draw(map)

; plot lat/lon of satellite
gsres@gsMarkerColor	 = "tomato3"
gsres@gsMarkerOpacityF = 0.5
gsres@gsMarkerThicknessF = sat_dot_size
gsn_polymarker(wks,map,lon,lat,gsres)


; sum up this individual satellite for this time period...
lat_x = lat
lon_x = lon
;qsort(lat_x) ; sorted in ascending order
;qsort(lon_x)

; sum up this individual satellite for this time period...
lat_x =	lat
lon_x =	lon
;qsort(lat_x)
;qsort(lon_x)
lat_x@_FillValue = -999
lon_x@_FillValue = -999

do y=0,dimsizes(lat_x)-1

 if (lat_x(y).ge.lat_s .and. lat_x(y).le.lat_n) then
  if (lon_x(y).ge.lon_w .and. lon_x(y).le.lon_e) then
   lat_x(y) = lat_x(y)
   lon_x(y) = lon_x(y)
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if

end do

count_lat_reg_sub = num(.not.ismissing(lat_x)) ; number of satellite data for this particular satellite at this time step
count_lon_reg_sub = num(.not.ismissing(lon_x))

lat_x = lat
lon_x = lon

do y=0,dimsizes(lat_x)-1

 if (lat_x(y).ge.lat_s_f .and. lat_x(y).le.lat_n_f) then
  if (lon_x(y).ge.lon_w_f .and. lon_x(y).le.lon_e_f) then
   lat_x(y) = lat_x(y)
   lon_x(y) = lon_x(y)
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if
  else
   lat_x(y) = -999
   lon_x(y) = -999
  end if

end do

count_lat_fl_sub = num(.not.ismissing(lat_x)) ; number of satellite data for this particular satellite at this time step
count_lon_fl_sub = num(.not.ismissing(lon_x))

print("Total regional " + sat_array(j) + " for date: " + dates_array(i) + " " + count_lat_reg_sub)
print("Total flight " + sat_array(j) + " for date: " + dates_array(i) + " " + count_lat_fl_sub)

sum_sat_reg(i,j) = count_lat_reg_sub
sum_sat_fl(i,j) = count_lat_fl_sub

; delete variables
delete([/lat,lat_x,lon_x,lon,w,filename/])

frame(wks)

else
print("No data for this date for individual satellite image: " + dates_array(i) + " and satellite: " + sat_array(j)) 
end if

end do ; end loop over satellites

end do ; end loop over dates

print("Writing satellite statistics to output text files...")

; write satellite statistics to output file
output_text_file = "all_satellites_ets_domain.txt"
system("rm -rf " + output_text_file)

do i=0,dimsizes(dates_array)-1

   if(i.eq.0) then
     header = [/"Date", "All-Sat","Cris NPP","Atms NPP","Amsua Aqua","Airs Aqua","DMSP F16","NOAA-15","NOAA-17","NOAA-18","NOAA-19","NOAA-17-hirs3","NOAA-19-hirs4","NOAA-18-mhs","NOAA-19-mhs","Amsua_metop-a","sndrd4_g13","sndrd3_g13","sndrd2_g13","sndrd1_g13","seviri_m09","mhs_metop-a","iasi_metop-a","hirs4_metop-a"/]
     write_table(output_text_file, "w", header, "%13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s")
   end if

   alist   = [/dates_array(i),sum_all_reg(i),sum_sat_reg(i,0),sum_sat_reg(i,1),sum_sat_reg(i,2),sum_sat_reg(i,3),sum_sat_reg(i,4),sum_sat_reg(i,5),sum_sat_reg(i,6),sum_sat_reg(i,7),sum_sat_reg(i,8),sum_sat_reg(i,9),sum_sat_reg(i,10),sum_sat_reg(i,11),sum_sat_reg(i,12),sum_sat_reg(i,13),sum_sat_reg(i,14),sum_sat_reg(i,15),sum_sat_reg(i,16),sum_sat_reg(i,17),sum_sat_reg(i,18),sum_sat_reg(i,19),sum_sat_reg(i,20),sum_sat_reg(i,21)/]
   write_table(output_text_file, "a", alist,"%13s %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13s %13i %13i %13i %13i %13i %13i %13i %13i")
exit
end do

output_text_file = "all_satellites_flight_domain.txt"
system("rm -rf " + output_text_file)

do i=0,dimsizes(dates_array)-1

   if(i.eq.0) then
     header = [/"Date", "All-Sat","Cris NPP","Atms NPP","Amsua Aqua","Airs Aqua","DMSP F16","NOAA-15","NOAA-17","NOAA-18","NOAA-19","NOAA-17-hirs3","NOAA-19-hirs4","NOAA-18-mhs","NOAA-19-mhs","Amsua_metop-a","sndrd4_g13","sndrd3_g13","sndrd2_g13","sndrd1_g13","seviri_m09","mhs_metop-a","iasi_metop-a","hirs4_metop-a"/]
     write_table(output_text_file, "w", header, "%13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s %13s")
   end if

   alist   = [/dates_array(i), sum_all_fl(i),sum_sat_fl(i,0),sum_sat_fl(i,1),sum_sat_fl(i,2),sum_sat_fl(i,3),sum_sat_fl(i,4),sum_sat_fl(i,5),sum_sat_fl(i,6),sum_sat_fl(i,7),sum_sat_fl(i,8),sum_sat_fl(i,9),sum_sat_fl(i,10),sum_sat_fl(i,11),sum_sat_fl(i,12),sum_sat_fl(i,13),sum_sat_fl(i,14),sum_sat_fl(i,15),sum_sat_fl(i,16),sum_sat_fl(i,17),sum_sat_fl(i,18),sum_sat_fl(i,19),sum_sat_fl(i,20),sum_sat_fl(i,21)/]
   write_table(output_text_file, "a", alist,"%13s %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13i %13s %13i %13i %13i %13i %13i %13i %13i %13i")
end do

print("Program ended successfully...")

end