FFEWS Daily Rainfall Plot using NCL

FFEWSRainDaily.ncl

; wrf_gsmcsn_raindaily.ncl
; wrf_gsmcsn_raindaily plots rainfall map on a 24 hour time period. 
; The script implicitly expect that the result is from 00Z, and the 
; netxt 24 hour interval is also on 00Z. 
; Author: Jamal Uddin Khan
; Email: [email protected]

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
	fid = "/home/ffews/FFEWS/WRFV3/run/wrfout_d01.nc";

	a = addfile(fid, "r");
	
	wks = gsn_open_wks("png", "/home/ffews/FFEWS/ffews.github.io/forecasts/rain/rain");
	
	firsttime = True;
	times = wrf_user_getvar(a, "times", -1);
	ntimes = dimsizes(times);
	
	do it = 0, ntimes-1, 8
		print("Processing: " + times(it));
		
		xlat = wrf_user_getvar(a,"XLAT", it);
		xlong = wrf_user_getvar(a,"XLONG", it);
		
		rainc = wrf_user_getvar(a, "RAINC", it);
		rainc@lat2d = xlat;
		rainc@lon2d = xlong;
		rainnc = wrf_user_getvar(a, "RAINNC", it);
		rainnc@lat2d = xlat;
		rainnc@lon2d = xlat;
		
		rain = rainc + rainnc;
		rain@lat2d = xlat;
		rain@lon2d = xlong;
		rain@units = "mm";
		rain@description = "24 Hour Rainfall";
		
		if(firsttime) then
			inittime = times(it);
			inittimestr = str_split(inittime, "_");
			firsttime = False;
			rainprev = rain;
		else
			rainplot = rain - rainprev;
			rainplot@lat2d = xlat;
			rainplot@lon2d = xlong;
			wrf_smooth_2d(rainplot, 1);
			rainprev = rain;
			validtime = times(it);
			validtimestr = str_split(validtime, "_");
			
			
			res := True;
			res@tiMainString = "WRF Rainfall (mm) Forecast";
			res@gsnCenterString = "Valid for " + inittimestr(0) + " " + "00Z" + " to " + validtimestr(0) + " " + "00Z";
			
			inittime = times(it);
			inittimestr = str_split(inittime, "_");
			
			res@gsnMaximize = True;
			
			res@cnFillOn = True;
			res@cnLevelSelectionMode = "ExplicitLevels";
			res@cnLevels = (/ 0, 0.1, 1, 2, 4, 8, 16, 32, 64, 128, 256/);
			res@cnFillColors = (/"White", "White", "lightskyblue1", "DarkOliveGreen1", "Green", "ForestGreen", "Yellow", "Orange", "Red", "Red4", "purple4"/)
			res@cnLinesOn = False;
			
			res@mpProjection = "Mercator";
			res@mpDataBaseVersion = "MediumRes";
			res@mpOutlineOn = True;
			res@mpOutlineBoundarySets = "AllBoundaries";
			res@mpPerimLineColor = "Black";
			res@mpGeophysicalLineThicknessF = 2.0;
			res@mpNationalLineThicknessF = 4.0;
			
			res@mpLimitMode = "LatLon";
			res@mpMinLatF = 20;min(rain@lat2d);22.875;
			res@mpMaxLatF = 28;max(rain@lat2d);26.125;
			res@mpMinLonF = 86;min(rain@lon2d);89.875;
			res@mpMaxLonF = 95;max(rain@lon2d);94.625;
			
			res@tmXBMinorOn = True
			res@tmXTMinorOn = True
			res@tmYLMinorOn = True
			res@tmYRMinorOn = True
			res@tmXTLabelsOn = True
			res@tmYRLabelsOn = True
			
			res@gsnAddCyclic = False
			
			res@pmTickMarkDisplayMode = "Always"
			
			res@lbOrientation = "vertical"
			res@lbLabelJust = "CenterLeft"
			
			contour = gsn_csm_contour_map(wks, rainplot, res); 
		end if
	end do
end