load "/usr/local/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "/usr/local/lib/ncarg/nclex/gsun/gsn_csm.ncl"
load "/usr/local/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  setvalues NhlGetWorkspaceObjectId()
	"wsMaximumSize": 134217728
  end setvalues

  ;diri   = "/home/Jacob05/grib/test/ens/gfs/"					; input directory
  ;fils   = systemfunc ("dir /b "+diri+"geavg.t00z.pgrbaf*")	      ; file paths

  data_dir = "/home/Jacob05/grib/test/gfs_20070213_06z/"
  gfs_HH   = "gfs_3_20070213_0600_"							;NCDC NOMADS archive of ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.2007021306/gfs.t06z.pgrbfHHH
  	;gfs0    = data_dir+gfs_HH+"000"
  	gfs1    = data_dir+gfs_HH+"003"
  	gfs2    = data_dir+gfs_HH+"006"
  	gfs3    = data_dir+gfs_HH+"009"
  	gfs4    = data_dir+gfs_HH+"012"
  	gfs5    = data_dir+gfs_HH+"015"
  	gfs6    = data_dir+gfs_HH+"018"
  	gfs7    = data_dir+gfs_HH+"021"
  	gfs8    = data_dir+gfs_HH+"024"
  	gfs9    = data_dir+gfs_HH+"027"
  	gfs10   = data_dir+gfs_HH+"030"
  	gfs11   = data_dir+gfs_HH+"033"
  	gfs12   = data_dir+gfs_HH+"036"
  	gfs13   = data_dir+gfs_HH+"039"
  	gfs14   = data_dir+gfs_HH+"042"
  	gfs15   = data_dir+gfs_HH+"045"
  	gfs16   = data_dir+gfs_HH+"048"
  	gfs17   = data_dir+gfs_HH+"051"
  	gfs18   = data_dir+gfs_HH+"054"
  	gfs19   = data_dir+gfs_HH+"057"
  	gfs20   = data_dir+gfs_HH+"060"
  	gfs21   = data_dir+gfs_HH+"063"
  	gfs22   = data_dir+gfs_HH+"066"
  	gfs23   = data_dir+gfs_HH+"069"
  	gfs24   = data_dir+gfs_HH+"072"
  	gfs25   = data_dir+gfs_HH+"075"
  	gfs26   = data_dir+gfs_HH+"078"
  	gfs27   = data_dir+gfs_HH+"081"
  	gfs28   = data_dir+gfs_HH+"084"
  	gfs29   = data_dir+gfs_HH+"087"
  	gfs30   = data_dir+gfs_HH+"090"
  	gfs31   = data_dir+gfs_HH+"093"
  	gfs32   = data_dir+gfs_HH+"096"
  	gfs33   = data_dir+gfs_HH+"099"
  	gfs34   = data_dir+gfs_HH+"102"
  	gfs35   = data_dir+gfs_HH+"105"
  	gfs36   = data_dir+gfs_HH+"108"
  	gfs37   = data_dir+gfs_HH+"111"
  	gfs38   = data_dir+gfs_HH+"114"
  	gfs39   = data_dir+gfs_HH+"117"
  	gfs40   = data_dir+gfs_HH+"120"
  	gfs41   = data_dir+gfs_HH+"123"
  	gfs42   = data_dir+gfs_HH+"126"
  	gfs43   = data_dir+gfs_HH+"129"
  	gfs44   = data_dir+gfs_HH+"132"
  	gfs45   = data_dir+gfs_HH+"135"
  	gfs46   = data_dir+gfs_HH+"138"
  	gfs47   = data_dir+gfs_HH+"141"
  	gfs48   = data_dir+gfs_HH+"144"
  	gfs49   = data_dir+gfs_HH+"147"
  	gfs50   = data_dir+gfs_HH+"150"
  	gfs51   = data_dir+gfs_HH+"153"
  	gfs52   = data_dir+gfs_HH+"156"
  	gfs53   = data_dir+gfs_HH+"159"
  	gfs54   = data_dir+gfs_HH+"162"
  	gfs55   = data_dir+gfs_HH+"165"
  	gfs56   = data_dir+gfs_HH+"168"
  	gfs57   = data_dir+gfs_HH+"171"
  	gfs58   = data_dir+gfs_HH+"174"
  	gfs59   = data_dir+gfs_HH+"177"
  	gfs60   = data_dir+gfs_HH+"180"

    ;fils = (/gfs0, gfs1, gfs2, gfs3, gfs4, gfs5, gfs6, gfs7, gfs8, gfs9, gfs10, gfs11, gfs12, gfs13, gfs14, gfs15, gfs16, gfs17, gfs18, gfs19, gfs20, gfs21, gfs22, gfs23, gfs24, gfs25, gfs26, gfs27, gfs28, gfs29, gfs30, gfs31, gfs32, gfs33, gfs34, gfs35, gfs36, gfs37, gfs38, gfs39, gfs40, gfs41, gfs42, gfs43, gfs44, gfs45, gfs46, gfs47, gfs48, gfs49, gfs50, gfs51, gfs52, gfs53, gfs54, gfs55, gfs56, gfs57, gfs58, gfs59, gfs60/) 
    fils = (/gfs1, gfs2, gfs3, gfs4, gfs5, gfs6, gfs7, gfs8, gfs9, gfs10, gfs11, gfs12, gfs13, gfs14, gfs15, gfs16, gfs17, gfs18, gfs19, gfs20, gfs21, gfs22, gfs23, gfs24, gfs25, gfs26, gfs27, gfs28, gfs29, gfs30, gfs31, gfs32, gfs33, gfs34, gfs35, gfs36, gfs37, gfs38, gfs39, gfs40, gfs41, gfs42, gfs43, gfs44, gfs45, gfs46, gfs47, gfs48, gfs49, gfs50, gfs51, gfs52, gfs53, gfs54, gfs55, gfs56, gfs57, gfs58, gfs59, gfs60/)

  print(fils)
  
  f      = addfiles(fils+".grb","r")  
  ListSetType (f, "join") 

  ; "Dimensions and sizes: [case | 60] x [lat_3 | 181] x [lon_3 | 360]"

  u	    = addfiles_GetVar(f,fils+".grb","U_GRD_3_HTGL_10")
  v 	    = addfiles_GetVar(f,fils+".grb","V_GRD_3_HTGL_10")
  tem     = addfiles_GetVar(f,fils+".grb","TMP_3_HTGL_10")

print(f)

  if (isfilevar(f, "CICEPL_3_SFC_ave6h"))
   	cpl = addfiles_GetVar(f,fils+".grb","CICEPL_3_SFC_ave6h")
  else 
   	cpl = addfiles_GetVar(f,fils+".grb","CICEPL_3_SFC_ave3h")
  end if 

  if (isfilevar(f, "CFRZR_3_SFC_ave6h"))
   	cfzra = addfiles_GetVar(f,fils+".grb","CFRZR_3_SFC_ave6h")
  else 
   	cfzra = addfiles_GetVar(f,fils+".grb","CFRZR_3_SFC_ave3h")
  end if 

  if (isfilevar(f, "A_PCP_3_SFC_acc6h"))
   	qpf = addfiles_GetVar(f,fils+".grb","A_PCP_3_SFC_acc6h")
  else 
   	qpf = addfiles_GetVar(f,fils+".grb","A_PCP_3_SFC_acc3h")
  end if 

  wnd     = (u^2+v^2)^(.5) 
  itime   = (/u@initial_time/)

print ("*************START ICE CALCS***************")

LuIce = qpf(:,:,:) ; preallocate space for new daily ppt array (metadata information kept)
JonesIce = qpf(:,:,:) ; preallocate space for new daily ppt array (metadata information kept)
JonesPIce = qpf(:,:,:) ; preallocate space for new daily ppt array (metadata information kept)

do t = 0,58,2
	do x = 0,180
		do y = 0,359
		;print(t+","+x+","+y+" : First Do Loop 't,x,y'")
			if (qpf(t,x,y) .gt. 0.0)
				if (tem(t,x,y) .lt. 273.3)
					if (cfzra(t,x,y) .gt. 0)
						if(cpl(t,x,y) .lt. 1)
							LuIce(t,x,y) = ((0.9999/0.92)*(qpf(t,x,y)/3.14159)*sqrt(1+(((wnd(t,x,y)*0.4)/4)^2)))
							JonesIce(t,x,y) = (3/(0.92*3.14159))*sqrt((((qpf(t,x,y)/3)*0.9999)^2)+((wnd(t,x,y)*0.2592*((qpf(t,x,y)/3)^0.88))^2))
							JonesPIce(t,x,y) = (3/(0.92*3.14159))*sqrt((((qpf(t,x,y)/3)*0.9999)^2)+((0*0.2592*((qpf(t,x,y)/3)^0.88))^2))
							;print(t+","+x+","+y+","+cra(t,x,y)+","+csn(t,x,y)+","+qpf(t,x,y)+","+cpl(t,x,y))
						else
							LuIce(t,x,y) = ((0.9999/0.92)*((qpf(t,x,y)/2)/3.14159)*sqrt(1+(((wnd(t,x,y)*0.4)/4)^2)))
							JonesIce(t,x,y) = (3/(0.92*3.14159))*sqrt((((qpf(t,x,y)/6)*0.9999)^2)+((wnd(t,x,y)*0.2592*((qpf(t,x,y)/6)^0.88))^2))
                  	                  JonesPIce(t,x,y) = (3/(0.92*3.14159))*sqrt((((qpf(t,x,y)/6)*0.9999)^2)+((0*0.2592*((qpf(t,x,y)/6)^0.88))^2))
							;print(t+","+x+","+y+","+cra(t,x,y)+","+csn(t,x,y)+","+qpf(t,x,y)+","+cpl(t,x,y))	
				end if
					else
						LuIce(t,x,y)     = 0.0
						JonesIce(t,x,y)  = 0.0
						JonesPIce(t,x,y) = 0.0
					end if 
				else
					LuIce(t,x,y)     = 0.0
					JonesIce(t,x,y)  = 0.0
					JonesPIce(t,x,y) = 0.0
				end if 
			else
				LuIce(t,x,y)     = 0.0
				JonesIce(t,x,y)  = 0.0
				JonesPIce(t,x,y) = 0.0
			end if
		end do
	end do
end do
delete(t)
delete(x)
delete(y)

do t = 1,59,2
	do x = 0,180
		do y = 0,359
		;print(t+","+x+","+y+" : First Do Loop 't,x,y'")
			if ((qpf(t,x,y)-qpf(t-1,x,y)) .gt. 0.0)
				if (tem(t,x,y) .lt. 273.3)
					if (cfzra(t,x,y) .gt. 0)
						if(cpl(t,x,y) .lt. 1)
							LuIce(t,x,y) = ((0.9999/0.92)*((qpf(t,x,y)-qpf(t-1,x,y))/3.14159)*sqrt(1+(((wnd(t,x,y)*0.4)/4)^2)))
							JonesIce(t,x,y) = (3/(0.92*3.14159))*sqrt(((((qpf(t,x,y)-qpf(t-1,x,y))/3)*0.9999)^2)+((wnd(t,x,y)*0.2592*(((qpf(t,x,y)-qpf(t-1,x,y))/3)^0.88))^2))
							JonesPIce(t,x,y) = (3/(0.92*3.14159))*sqrt(((((qpf(t,x,y)-qpf(t-1,x,y))/3)*0.9999)^2)+((0*0.2592*(((qpf(t,x,y)-qpf(t-1,x,y))/3)^0.88))^2))
							;print(t+","+x+","+y+","+cra(t,x,y)+","+csn(t,x,y)+","+(qpf(t,x,y)-qpf(t-1,x,y))+","+cpl(t,x,y))
						else
							LuIce(t,x,y) = ((0.9999/0.92)*(((qpf(t,x,y)-qpf(t-1,x,y))/2)/3.14159)*sqrt(1+(((wnd(t,x,y)*0.4)/4)^2)))
							JonesIce(t,x,y) = (3/(0.92*3.14159))*sqrt(((((qpf(t,x,y)-qpf(t-1,x,y))/6)*0.9999)^2)+((wnd(t,x,y)*0.2592*(((qpf(t,x,y)-qpf(t-1,x,y))/6)^0.88))^2))
                  	                  JonesPIce(t,x,y) = (3/(0.92*3.14159))*sqrt(((((qpf(t,x,y)-qpf(t-1,x,y))/6)*0.9999)^2)+((0*0.2592*(((qpf(t,x,y)-qpf(t-1,x,y))/6)^0.88))^2))
							;print(t+","+x+","+y+","+cra(t,x,y)+","+csn(t,x,y)+","+(qpf(t,x,y)-qpf(t-1,x,y))+","+cpl(t,x,y))	
				end if
					else
						LuIce(t,x,y)     = 0.0
						JonesIce(t,x,y)  = 0.0
						JonesPIce(t,x,y) = 0.0
					end if 
				else
					LuIce(t,x,y)     = 0.0
					JonesIce(t,x,y)  = 0.0
					JonesPIce(t,x,y) = 0.0
				end if 
			else
				LuIce(t,x,y)     = 0.0
				JonesIce(t,x,y)  = 0.0
				JonesPIce(t,x,y) = 0.0
			end if
		end do
	end do
end do
delete(t)
delete(x)
delete(y)

delme = LuIce(lat_3|:,lon_3|:,case|:) ; reorder LuIce
LuIceRunAccum = qpf(:,:,:) ; preallocate space for new array (metadata information kept)

do t = 0,33
	LuIceRunAccum(t,:,:) = (/dim_sum(delme(:,:,0:t))/)
end do

delete(delme) 

;-------

delme = JonesIce(lat_3|:,lon_3|:,case|:) ; reorder JonesIce
JonesIceRunAccum = qpf(:,:,:) ; preallocate space for new array (metadata information kept)

do t = 0,33
	JonesIceRunAccum(t,:,:) = (/dim_sum(delme(:,:,0:t))/)
end do

delete(delme) 

;-------

delme = JonesPIce(lat_3|:,lon_3|:,case|:) ; reorder JonesPIce
JonesPIceRunAccum = qpf(:,:,:) ; preallocate space for new array (metadata information kept)

do t = 0,33
	JonesPIceRunAccum(t,:,:) = (/dim_sum(delme(:,:,0:t))/)
end do

delete(delme)


print ("*************END ICE CALCS******************")

print ("*************create plots*******************")

  wks  = gsn_open_wks ("pdf","/home/Jacob05/results/Accum_Ice_Plot_GFS")          ; open workstation  
  gsn_define_colormap (wks,"gui_default")     						    ; choose color map
  ;gsn_define_colormap (wks,"WhViBlGrYeOrRe")     					    ; choose color map
  
  res                        = True           						    ; plot mods desired for original grid
  res@gsnMaximize            = True
  res@cnFillOn               = True            		     				    ; color fill  
  res@cnLinesOn              = False          						    
  res@gsnSpreadColors        = True               					    ; use total colormap (if True)
  res@gsnSpreadColorStart    = 0 
  res@gsnSpreadColorEnd      = 23   
  res@cnLevelSpacingF	     = 0.125  						     	    ; contour levels 
  res@mpGridAndLimbOn        = True
  res@pmTickMarkDisplayMode  = "Always"          					    ; turn on tickmarks
  res@tmXTOn                 = False  
  res@gsnAddCyclic           = False                                              ; regional data
  res@mpFillOn               = False            

; These next three resources set up the drawing of county boundaries.
;   res@mpOutlineBoundarySets  = "AllBoundaries"
;   res@mpDataBaseVersion      = "Ncarg4_1"
;   res@mpDataSetName          = "Earth..2" 

  res@mpMinLatF               = -90		       					    ; specify region to be plotted
  res@mpMaxLatF               = 90
  res@mpMinLonF               = 0
  res@mpMaxLonF               = 359
  
  res@mpGridLonSpacingF       = 10 								    ; Lon spacing (degrees)
  res@mpGridLatSpacingF       = 10 								    ; Lon spacing (degrees)
  res@mpGridLineDashPattern   = 2								    ; 0 (solid) default (2=dots) 

  res@gsnStringFontHeightF    = 0.0125
  res@lbLabelFontHeightF      = 0.015

  x = 33												    ; Accumulation through x period (x = x_hr/3)

; Lu Prime Plot
  res@gsnRightString         = "Lu' Req(in) Hr 00 to "+((cpl&case(x)+1)*6)			 		; draw center subtitle
  res@gsnLeftString           = itime+"Z GEFS"				       					; draw left subtitle
  plot = gsn_csm_contour_map(wks,LuIceRunAccum(x,:,:)/25.4,res)
  print("ploting Lu Prime")

; Jones (w/ Marshal-Palmer LWC) Plot
  res@gsnRightString         = "Jones Req(in) Hr 00 to "+((cpl&case(x)+1)*6)					; draw center subtitle
  res@gsnLeftString           = itime+"Z GEFS"				       					; draw left subtitle
  plot = gsn_csm_contour_map(wks,JonesIceRunAccum(x,:,:)/25.4,res)
  print("ploting Jones")

; Jones Prime (No Wind) (w/ Marshal-Palmer LWC) Plot
  res@gsnRightString         = "Jones' Req(in) Hr 00 to "+((cpl&case(x)+1)*6) 				; draw center subtitle
  res@gsnLeftString           = itime+"Z GEFS"				       					; draw left subtitle
  plot = gsn_csm_contour_map(wks,JonesPIceRunAccum(x,:,:)/25.4,res)
  print("ploting Jones Prime")

; Wtd Mean Plot
  res@gsnStringFontHeightF    = 0.0075
  res@gsnRightString         = "(J*0.10)+(J'*0.60)+(L'*0.30) Req(in) Hr 00 to "+((cpl&case(x)+1)*6)	; draw center subtitle
  res@gsnLeftString           = itime+"Z GEFS"				       					; draw left subtitle
  plot = gsn_csm_contour_map(wks,((LuIceRunAccum(x,:,:)*0.30)+(JonesPIceRunAccum(x,:,:)*0.60)+(JonesIceRunAccum(x,:,:)*0.10))/25.4,res)
  print("ploting Wtd Mean")

;do x = 0,33	
	;res@cnLevelSelectionMode    = "AutomaticLevels"
	;res@cnInfoLabelOrthogonalPosF = -0.07          		 							; move label inside plot
	;res@gsnRightString          = "qpf at "+((cpl&case(x)+1)*6)				 		; draw center subtitle
	;res@gsnLeftString           = itime+"Z GEFS"	
	;	plot = gsn_csm_contour_map(wks,qpf(x,:,:)/25.4,res)
	;res@gsnRightString          = "Temperature at "+((cpl&case(x)+1)*6)				 		; draw center subtitle
	;res@gsnLeftString           = itime+"Z GEFS"
	;	plot = gsn_csm_contour_map(wks,tem(x,:,:),res)
	;res@gsnRightString          = "CPOFZRA at "+((cpl&case(x)+1)*6)				 		; draw center subtitle
	;res@gsnLeftString           = itime+"Z GEFS"
	;	plot = gsn_csm_contour_map(wks,cfzra(x,:,:),res)
	;res@gsnRightString          = "CPPL at "+((cpl&case(x)+1)*6)				 		; draw center subtitle
	;res@gsnLeftString           = itime+"Z GEFS"
	;	plot = gsn_csm_contour_map(wks,cpl(x,:,:),res)
;end do	

end