; Plot the observed precipitation for the Feb 21 storm of the third OSE case in 2016
; Uses radar and station data to plot the 24 hour and accumulated precipitation plots over Alaska and CONUS

; Andrew Kren
; May 20, 2016
;;;;;;;;;;;;;;;;;;;;;;;;;;;;


begin

  obs_dir = "/scratch4/BMC/shout/Andrew.Kren/ENRR_obs_precip/"
  
 ; list files in directory to determine how many dropsondes to read
 region = "conus"  ; region to plot precipitation, ak for alaska, conus for conus, this is set in the filename of netcdf
 
 
 fils = systemfunc("csh -c 'cd " + obs_dir + " ; ls *" + region +"*.nc'")
 
 
 if (region.eq."conus") then
  nlat = 813
  nlon = 1051
 else
  nlat = 530
  nlon = 460
 end if
 
 if (region.eq."conus") then
  lat_low = 32    ; coordinates used to plot the verification region for CA
  lat_high = 46
  lon_low = 231
  lon_high = 244
 else
  lat_low = 55    ; coordinates used to plot the verification region for AK
  lat_high = 65
  lon_low = 193
  lon_high = 220
 end if
 
; --------------------------- end of user's settings -----------------------------------------


hourly_pcp = new((/dimsizes(fils),nlat,nlon/),float) ; 24 hr accumulated precip for the current day
       hourly_pcp!0="hour"
       hourly_pcp!1="lat"
       hourly_pcp!2="lon"
       hourly_pcp@_FillValue = -999
accum_pcp = new((/dimsizes(fils),nlat,nlon/),float)  ; accumulated precipitation over each day
       accum_pcp!0="hour"
       accum_pcp!1="lat"
       accum_pcp!2="lon"
       accum_pcp@_FillValue = -999
       
if (region.eq."conus") then
year = tointeger(str_get_cols(fils,17,20))    ; year, month, day, hour, minute, sec in UTC, same for hamsr
month = tointeger(str_get_cols(fils,21,22))
day = tointeger(str_get_cols(fils,23,24))
else
year = tointeger(str_get_cols(fils,14,17))    ; year, month, day, hour, minute, sec in UTC, same for hamsr
month = tointeger(str_get_cols(fils,18,19))
day = tointeger(str_get_cols(fils,20,21))
end if

do i=0,dimsizes(fils)-1                                 ; loop over all daily files

print("Reading YMD " + tostring(year(i))+tostring(month(i))+tostring(day(i)))


print("Reading file: " + fils(i))

        w = addfile(obs_dir+fils(i),"r")
        lat_corners = w->lat
        lon_corners = w->lon
	true_lat = w->true_lat
	true_lon = w->true_lon
	factor = 0.0003937008
        pcp = factor*tofloat(w->amountofprecip) ; to inches
      
hourly_pcp(i,:,:) = (/pcp/)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      ; beginning of plot section
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

outfile = "pcp_" + region + "_" + tostring(year(i)) + tostring(month(i)) + tostring(day(i))

wks = gsn_open_wks("eps",outfile)
gsn_define_colormap(wks,"WhBlGrYeRe")


res = True
res@mpProjection = "Stereographic"
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat_corners(0)
res@mpLeftCornerLonF = lon_corners(0)
res@mpRightCornerLatF = lat_corners(2)
res@mpRightCornerLonF = lon_corners(2)
res@mpEllipticalBoundary = True

res@mpRelativeCenterLon         = True           ; set a center lon
res@mpCenterLonF                = true_lon       ; center lon
res@mpRelativeCenterLat         = True           ; set a center lat
res@mpCenterLatF                = true_lat       ; center lat
res@tfDoNDCOverlay              = True 

res@cnLinesOn            = False
res@gsnMaximize          = True
res@cnFillOn             = True               ; color plot desired
res@cnLineLabelsOn	 = False              ; turn off contour lines
res@tiMainString         = "24-hr Accumulated Precipitation: " + tostring(month(i)) + "/" + tostring(day(i)) + "/" + tostring(year(i))
;res@mpMinLatF            = 15                 ; specify min lat
;res@mpMaxLatF            = 70
;res@mpMinLonF            = 150
;res@mpMaxLonF            = 270
;res@mpCenterLonF         = 220
res@gsnRightString = ""
res@gsnLeftString = ""
res@mpGridAndLimbOn = True
res@mpGridLatSpacingF = 10.0
res@mpGridLonSpacingF = 10.0
res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
;res@mpOutlineBoundarySets       = "National"   ; turn on country boundaries
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
res@cnFillOpacityF = 0.9
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/ 0.01,0.05,0.10,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4/)
res@lbLabelStrings = (/"0.01","0.05","0.10", "0.25", "0.5", "0.75", "1", "1.5", "2","2.5","3","3.5","4"/)
res@cnExplicitLabelBarLabelsOn = True
res@lbTitleOn = True
res@lbTitleString    = "Precipitation (in)"                ; title string
res@lbTitlePosition  = "Top"              ; title position
res@lbLabelFontHeightF =.022                 ; make labels larger
res@lbTitleFontHeightF= .015                 ; make title smaller
res@pmLabelBarOrthogonalPosF = .10           ; move whole thing down

plot = gsn_csm_contour_map(wks,hourly_pcp(i,:,:),res)

 gres = True
 gres@gsLineColor = "Red"
 gres@gsLineThicknessF = 5.0
 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/)
 lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/)

 dum = gsn_add_polyline(wks,plot,lon2,lat2,gres)

      draw(plot)
      frame(wks) ; advance frame

delete(res)


delete([/lat_corners,lon_corners,true_lat,true_lon,pcp/])

end do

print("Calculating accumulated precipitation")
ntstrt = 0 ; counter variables for accumulated precipitation
ntlast = 0
    do t=0,dimsizes(fils)-1
    accum_pcp(t,:,:) = dim_sum_n(hourly_pcp(ntstrt:ntlast,:,:),0)
    ntlast = ntlast + 1   ; increment for accumulated precipitation
    end do
    
    exit

print("Making accumulated pcp plots...")

;;;;;;;;;;;;;;;;;;;
; Plot accumulated precipitation
;;;;;;;;;;;;;;;;;;;
;do i=0,fc_files-1 




;end do


print("Program ended successfully")
end