Index: configs/dcmip2016/PLOTS/copy_graphs.sh
===================================================================
--- configs/dcmip2016/PLOTS/copy_graphs.sh	(revision 437)
+++ configs/dcmip2016/PLOTS/copy_graphs.sh	(revision 437)
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+# collects all graphs from subdirectories
+# run first "make_all_graphs" to generate graphs
+
+ROOT=$(pwd)
+
+for CASE in test* ; do
+    cd $ROOT/$CASE/data
+    EXPS=$(ls)
+    for EXP in $EXPS ; do
+	echo $EXP
+	cd $ROOT/$CASE/data/$EXP
+	for IMG in $(ls *.png *.gif) ; do
+#	    echo $IMG ../../graphs/${CASE}_${EXP}_${IMG}
+	    cp $IMG $ROOT/graphs/${CASE}_${EXP}_${IMG}
+	done
+    done
+done
Index: configs/dcmip2016/PLOTS/make_graphs.sh
===================================================================
--- configs/dcmip2016/PLOTS/make_graphs.sh	(revision 437)
+++ configs/dcmip2016/PLOTS/make_graphs.sh	(revision 437)
@@ -0,0 +1,28 @@
+#!/bin/bash
+
+# this script should be run from a test case directory (test1, test2, test3)
+# not from the main directory (PLOTS), where it is stored
+# from PLOTS, run "make_all_graphs" to make all graphs
+
+# execute "bash -l" before running this script
+
+module load ncl
+
+ROOT=$(pwd)
+CASE=$(basename $ROOT)
+
+cd data
+EXPS=$(ls)
+
+cd $ROOT
+for EXP in $EXPS ; do
+    echo $EXP
+    cd $ROOT/data/$EXP
+    rm -f *.png *.gif
+    ln -s ../../include/* .
+    for NCL in ../../*.ncl ; do
+	echo $NCL
+	ncl $NCL &
+    done
+    wait
+done
Index: configs/dcmip2016/PLOTS/make_links.sh
===================================================================
--- configs/dcmip2016/PLOTS/make_links.sh	(revision 437)
+++ configs/dcmip2016/PLOTS/make_links.sh	(revision 437)
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+# run this script first 
+
+# adapt the line below as needed
+RUN_DIR=/glade/scratch/yannm/DCMIP2016/RUNS/
+
+# the runs must be organized as : $RUN_DIR/TESTCASE/VARIANT
+# with TESTCASE = SUPERCELL ...
+# and VARIANT = dry ...
+# this structure will be reflected in the form tescase/TESTCASE/variant/VARIANT
+# data/TESTCASE/ contains NCL scripts for each test case
+
+ROOT=$(pwd)
+
+cd $RUN_DIR
+TESTS=$(ls)
+echo "Test cases : " $TESTS
+for CASE in $TESTS ; do
+    cd $RUN_DIR/$CASE
+    VARIANTS=$(ls)
+    for VARIANT in $VARIANTS ; do
+	DIR=$ROOT/testcase/$CASE/variant/$VARIANT
+	mkdir -p $DIR
+	cd $DIR
+	for x in $RUN_DIR/$CASE/$VARIANT/*.nc ; do
+	    echo $x
+	    ln -sf $x .
+	done
+    done
+done
Index: configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS.ncl	(revision 437)
@@ -0,0 +1,265 @@
+;================================================;
+;	Example ncl script to produce the set of
+;	lat-lon plots for DCMIP-2016 test case 1
+; Created by James Kent, Christiane Jablonowski
+;	and Paul Ullrich (University of Michigan) for DCMIP-2012 
+;
+; Modified by Colin Zarzycki for DCMIP-2016
+;
+; v1.01 - 6/7/16 
+; v1.02 - 6/8/16 - Typo in var_choice options
+;================================================;
+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/csm/contributed.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
+; ================================================;
+
+; NOTE, if your model does not output T850 OR P at level
+; midpoints OR you cannot calculate P from hybrid levels, you
+; will need to make model specific additions denoted by
+; CALCPHERE below
+
+; PRECL needs to be in m/s for contours to plot correctly
+
+begin
+
+;=================================================;
+; open file and read in data 
+; GIVEN AS (time,lev,lat,lon) from 0 to n-1
+;=================================================;
+
+  ; NOTE, that this file needs to be a catted file including all times at daily interval
+  ; if you are outputting 1 file per time, you can run "ncrcat dcmip1.files.*.nc cat.nc"
+  ; for one file
+
+; TD	f = addfile("/home/zarzycki/dcmip/results/dcmip1/newcat.nc","r")
+   f = addfile("output_dcmip2016_regular.nc","r")
+		; Input useful parameters
+
+	lat  = f->lat
+  lon  = f->lon
+	lev  = f->lev
+	nlat   = getfilevardimsizes(f, "lat" )
+	nlon   = getfilevardimsizes(f, "lon" )
+	nlev   = getfilevardimsizes(f, "lev" )
+
+  ; We want plots of temperature, vertical velocity and relative vorticity 
+  ; at the 850 hPa level. We also want surface pressure.
+
+  ; Select var_choice 	1 - Temperature at 850 hPa
+  ; 		        2 - Large-scale precipitation
+  ;			3 - Surface Pressure
+
+	var_choice = 3
+	
+	if (var_choice .eq. 1) then				; Select T
+
+		if (isfilevar(f,"T850")) then
+
+			varload = f->T850
+
+		else	
+
+			; we interp T to the 850 hPa level
+			varload0 = f->T				; T is a 4D var
+			varload  = f->T(:,0,:,:)
+			; Calculate P using hybrid coefficients
+			; If you have pressure as a variable just do
+			; P = f->P
+      ; CALCPHERE
+			hyam = f->hyam
+  		hybm = f->hybm
+  		P0   = f->P0
+			PS   = f->PS
+			P    = 0.0*varload0
+			P = pres_hybrid_ccm(PS,P0,hyam,hybm)
+			P!0 = "time"
+			P!1 = "lev"
+			P!2 = "lat"
+			P!3 = "lon"
+			; interp to 850 hPa level
+			plevel = 85000.0
+			varload = int2p_n(P,varload0,plevel,2,1)
+
+		end if
+
+	else if (var_choice .eq. 2) then			; Select Large-scale precipitation
+
+		varload = f->PRECL
+    varload=varload*8.64e7		
+
+	else if (var_choice .eq. 3) then			; Select PS
+
+		varload = f->PS
+
+	end if
+	end if	
+	end if
+
+		; We want the output at days 6, 9, 12 and 15. If your model 
+		; output is in 6 hourly intervals this becomes indices 24, 36, 48 and 60
+
+	output_per_day=8
+	day6  = 6*output_per_day
+	day9  = 9*output_per_day
+	day12 = 12*output_per_day
+	day15 = 15*output_per_day
+
+	var1 = varload(day6,:,:)
+	var2 = varload(day9,:,:)
+	var3 = varload(day12,:,:)
+	var4 = varload(day15,:,:)
+
+		; We don't want long-name in the plots
+
+	var1@long_name = " "
+	var2@long_name = " "
+	var3@long_name = " "
+	var4@long_name = " "
+
+		; Delete loaded data
+
+	delete(varload)
+
+		; We now produce the plot
+
+	plot  = new (4, graphic)               			; define plot - need 4 panels
+
+	res1                      = True
+	res1@gsnDraw              = False          		; panel plot 
+	res1@gsnFrame             = False			; don't draw yet
+	res1@cnFillOn             = True
+	res1@cnLinesOn            = True
+	res1@gsnSpreadColors      = True
+	res1@lbLabelAutoStride    = True
+	res1@gsnCenterString      = ""
+	res1@tiMainString         = ""
+	res1@vpWidthF             = 0.38
+	res1@vpHeightF            = 0.19
+	res1@cnLevelSelectionMode = "ManualLevels"
+	res1@cnInfoLabelOn        = False			; don't give each
+	res1@cnLineLabelsOn       = False			; panel its own
+	res1@lbLabelBarOn         = False			; label bar
+
+	pres                          = True
+	pres@gsnMaximize              = True  
+	pres@gsnPanelLabelBar         = True			; Communal label bar
+	pres@gsnPanelLeft             = 0.1
+	pres@gsnPanelRight            = 0.9
+	pres@pmLabelBarOrthogonalPosF = -0.03
+	pres@gsnFrame                 = False
+	pres@lbLabelStride            = 1
+
+	res1@sfXArray        = lon            			; uses lon as plot x-axis
+	res1@sfYArray        = lat            			; uses lat for y axis
+	res1@trYReverse      = False          			; reverses y-axis, false
+	res1@tiYAxisString   = ""   	       			; y-axis title
+	res1@tiXAxisString   = ""             			; x-axis title
+
+	if (var_choice .eq. 1) then 	;================ Temperature plot ====================
+ 
+	pltTitle="DYNAMICO Test DCMIP-2016-1, T850 hPa" 			; Plot title if required
+	pres@txString = pltTitle
+
+		; Change the output type and name
+
+	;wks  = gsn_open_wks("X11","cam-fv-test42latlonT850")   	; output using X11
+	;wks  = gsn_open_wks("eps","cam-fv-test42latlonT850")   	; output using eps
+	wks  = gsn_open_wks("png","DYNAMICO-test42latlonT850")   	; output using png
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 300.0				; max contour color label
+	res1@cnMinLevelValF  = 230.0				; min contour color label
+	res1@cnLevelSpacingF = 10.0				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 2) then 	;================ Large-scale precipitation plot ====================
+ 
+	pltTitle="DYNAMIXO Test 42, Large-scale precipitation" 	; Plot title if required
+	pres@txString = pltTitle
+
+	var1@units = " "
+	var2@units = " "
+	var3@units = " "
+	var4@units = " "
+
+	res1@cnLinesOn            = False				; Turn lines off for clearer plot
+
+		; Change the output type and name
+
+	;wks  = gsn_open_wks("X11","cam-fv-test42latlonPRECL")   	; output using X11
+	;wks  = gsn_open_wks("eps","cam-fv-test42latlonPRECL")   	; output using eps
+	wks  = gsn_open_wks("png","test42latlonPRECL")   	; output using png
+	gsn_define_colormap(wks,"gui_default")  
+
+		; Note that the maximum might be larger than these contour spacing values
+
+	res1@cnMaxLevelValF  = 50.0				; max contour color label
+	res1@cnMinLevelValF  = 0.0				; min contour color label
+	res1@cnLevelSpacingF = 4.0				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 3) then 	;================ Surface Pressure plot ====================
+ 
+	pltTitle="DYNAMICO Test 42, PS" 			; Plot title if required
+	pres@txString = pltTitle
+
+		; Convert to hPa 
+
+	var1=var1/100.0
+	var2=var2/100.0
+	var3=var3/100.0
+	var4=var4/100.0
+
+	var1@units = "hPa "
+	var2@units = "hPa "
+	var3@units = "hPa "
+	var4@units = "hPa "
+
+		; Change the output type and name
+
+	;wks  = gsn_open_wks("X11","cam-fv-test42latlonPS")   	; output using X11
+	;wks  = gsn_open_wks("eps","cam-fv-test42latlonPS")   	; output using eps
+	wks  = gsn_open_wks("png","test42latlonPS")   	; output using png
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 1020.0				; max contour color label
+	res1@cnMinLevelValF  = 920.0				; min contour color label
+	res1@cnLevelSpacingF = 10.0				; contour color spacing			; choose a colormap
+
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+
+	res1@gsnCenterString = "t = 6"
+	plot(0) = gsn_csm_contour(wks,var1(:,:),res1)		; plot var1
+
+	res1@gsnCenterString = "t = 9"
+     	plot(1) = gsn_csm_contour(wks,var2(:,:),res1)		; plot var2
+
+	res1@gsnCenterString = "t = 12"
+	plot(2) = gsn_csm_contour(wks,var3(:,:),res1)		; plot var3
+
+	res1@gsnCenterString = "t = 15"
+     	plot(3) = gsn_csm_contour(wks,var4(:,:),res1)		; plot var4
+
+	gsn_panel(wks,plot,(/2,2/),pres)			; 2x2 plot
+
+		; Add latitude and longitude labels
+
+	txres3                = True
+	txres3@txAngleF       = 90.
+	txres3@txFontHeightF  = 0.02
+	gsn_text_ndc(wks,"Latitude",0.08,0.49,txres3)
+
+	txres3@txAngleF       = 0.
+	gsn_text_ndc(wks,"Longitude",0.5,0.22,txres3)
+
+	frame(wks)
+
+end
+
+
+
+
+
+
Index: configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS_LOOP.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS_LOOP.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS_LOOP.ncl	(revision 437)
@@ -0,0 +1,280 @@
+;================================================;
+;	Example ncl script to produce the set of
+;	lat-lon plots for DCMIP-2016 test case 1
+; Created by James Kent, Christiane Jablonowski
+;	and Paul Ullrich (University of Michigan) for DCMIP-2012 
+;
+; Modified by Colin Zarzycki for DCMIP-2016
+; Further modified by Maximo Menchaca for creating gifs
+;
+; v1.01 - 6/7/16 
+; v1.02 - 6/8/16 - Typo in var_choice options
+;================================================;
+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/csm/contributed.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
+; ================================================;
+
+; NOTE, if your model does not output T850 OR P at level
+; midpoints OR you cannot calculate P from hybrid levels, you
+; will need to make model specific additions denoted by
+; CALCPHERE below
+
+; PRECL needs to be in m/s for contours to plot correctly
+
+begin
+
+;=================================================;
+; open file and read in data 
+; GIVEN AS (time,lev,lat,lon) from 0 to n-1
+;=================================================;
+
+  ; NOTE, that this file needs to be a catted file including all times at daily interval
+  ; if you are outputting 1 file per time, you can run "ncrcat dcmip1.files.*.nc cat.nc"
+  ; for one file
+
+	f = addfile("output_dcmip2016_regular.nc","r")
+
+use_pdf = False
+
+		; Input useful parameters
+
+	lat  = f->lat
+        lon  = f->lon
+	lev  = f->lev
+	nlat   = getfilevardimsizes(f, "lat" )
+	nlon   = getfilevardimsizes(f, "lon" )
+	nlev   = getfilevardimsizes(f, "lev" )
+
+  ; We want plots of temperature, vertical velocity and relative vorticity 
+  ; at the 850 hPa level. We also want surface pressure.
+
+  ; Select var_choice 	1 - Temperature at 850 hPa
+  ; 		        2 - Large-scale precipitation
+  ;			3 - Surface Pressure
+
+	do var_choice = 1,3
+;	var_choice = 1 
+	
+	if (var_choice .eq. 1) then				; Select T
+
+		if (isfilevar(f,"T850")) then
+
+			varload = f->T850
+
+		else	
+
+			; we interp T to the 850 hPa level
+			varload0 = f->T				; T is a 4D var
+			varload  = f->T(:,0,:,:)
+			; Calculate P using hybrid coefficients
+			; If you have pressure as a variable just do
+			; P = f->P
+      ; CALCPHERE
+			hyam = f->hyam
+  		hybm = f->hybm
+  		P0   = f->P0
+			PS   = f->PS
+			P    = 0.0*varload0
+			P = pres_hybrid_ccm(PS,P0,hyam,hybm)
+			P!0 = "time"
+			P!1 = "lev"
+			P!2 = "lat"
+			P!3 = "lon"
+			; interp to 850 hPa level
+			plevel = 85000.0
+			varload = int2p_n(P,varload0,plevel,2,1)
+
+		end if
+
+	filename_string = "DYNAMICO-case161latlonT850"
+
+	else if (var_choice .eq. 2) then			; Select Large-scale precipitation
+
+		varload = f->PRECT
+    varload=varload*8.64e7		
+	filename_string = "DYNAMICO-case161latlonPRECT"
+
+	else if (var_choice .eq. 3) then			; Select PS
+	filename_string = "DYNAMICO-case161latlonPS"
+
+		varload = f->PS
+
+	end if
+        end if
+        end if
+
+if(use_pdf) then
+	wks  = gsn_open_wks("pdf",filename_string)   	; output using X11
+else
+	wks  = gsn_open_wks("png",filename_string)   	; output using png
+end if
+
+		; We want the output at days 6, 9, 12 and 15. If your model 
+		; output is in 6 hourly intervals this becomes indices 24, 36, 48 and 60
+
+	do d = 1,240,2
+	  day_curr = d/8
+
+;	day6  = 24 
+;	day9  = 36
+;	day12 = 48
+;	day15 = 60
+
+	var1 = varload(d,:,:)
+;	var2 = varload(day9,:,:)
+;	var3 = varload(day12,:,:)
+;	var4 = varload(day15,:,:)
+
+		; We don't want long-name in the plots
+
+	var1@long_name = " "
+;	var2@long_name = " "
+;	var3@long_name = " "
+;	var4@long_name = " "
+
+		; Delete loaded data
+
+;	delete(varload)
+
+		; We now produce the plot
+
+;	plot  = new (4, graphic)               			; define plot - need 4 panels
+
+	res1                      = True
+	res1@gsnDraw              = False          		; panel plot 
+	res1@gsnFrame             = False			; don't draw yet
+	res1@cnFillOn             = True
+	res1@cnLinesOn            = True
+	res1@gsnSpreadColors      = True
+	res1@lbLabelAutoStride    = True
+	res1@gsnCenterString      = ""
+	res1@tiMainString         = ""
+	res1@vpWidthF             = 0.38
+	res1@vpHeightF            = 0.19
+	res1@cnLevelSelectionMode = "ManualLevels"
+	res1@cnInfoLabelOn        = False			; don't give each
+	res1@cnLineLabelsOn       = False			; panel its own
+	res1@lbLabelBarOn         = False			; label bar
+
+	pres                          = True
+	pres@gsnMaximize              = True  
+	pres@gsnPanelLabelBar         = True			; Communal label bar
+	pres@gsnPanelLeft             = 0.1
+	pres@gsnPanelRight            = 0.9
+	pres@pmLabelBarOrthogonalPosF = -0.03
+	pres@gsnFrame                 = False
+	pres@lbLabelStride            = 1
+
+	res1@sfXArray        = lon            			; uses lon as plot x-axis
+	res1@sfYArray        = lat            			; uses lat for y axis
+	res1@trYReverse      = False          			; reverses y-axis, false
+	res1@tiYAxisString   = ""   	       			; y-axis title
+	res1@tiXAxisString   = ""             			; x-axis title
+
+	if (var_choice .eq. 1) then 	;================ Temperature plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, T850 hPa" 			; Plot title if required
+	pres@txString = pltTitle
+
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 300.0				; max contour color label
+	res1@cnMinLevelValF  = 230.0				; min contour color label
+	res1@cnLevelSpacingF = 10.0				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 2) then 	;================ Large-scale precipitation plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, Large-scale precipitation" 	; Plot title if required
+	pres@txString = pltTitle
+
+	var1@units = " "
+;	var2@units = " "
+;	var3@units = " "
+;	var4@units = " "
+
+	res1@cnLinesOn            = False				; Turn lines off for clearer plot
+
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+		; Note that the maximum might be larger than these contour spacing values
+
+	res1@cnMaxLevelValF  = 50.0				; max contour color label
+	res1@cnMinLevelValF  = 0.0				; min contour color label
+	res1@cnLevelSpacingF = 4.0				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 3) then 	;================ Surface Pressure plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, PS" 			; Plot title if required
+	pres@txString = pltTitle
+
+		; Convert to hPa 
+
+	var1=var1/100.0
+;	var2=var2/100.0
+;	var3=var3/100.0
+;	var4=var4/100.0
+
+	var1@units = "hPa "
+;	var2@units = "hPa "
+;	var3@units = "hPa "
+;	var4@units = "hPa "
+
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 1020.0				; max contour color label
+	res1@cnMinLevelValF  = 920.0				; min contour color label
+	res1@cnLevelSpacingF = 10.0				; contour color spacing			; choose a colormap
+
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+
+	res1@gsnCenterString = "t = " + day_curr + " days"
+	plot = gsn_csm_contour(wks,var1(:,:),res1)		; plot var1
+
+;	res1@gsnCenterString = "t = 9"
+;     	plot(1) = gsn_csm_contour(wks,var2(:,:),res1)		; plot var2
+;
+;	res1@gsnCenterString = "t = 12"
+;	plot(2) = gsn_csm_contour(wks,var3(:,:),res1)		; plot var3
+;
+;	res1@gsnCenterString = "t = 15"
+;     	plot(3) = gsn_csm_contour(wks,var4(:,:),res1)		; plot var4
+
+	gsn_panel(wks,plot,(/1,1/),pres)			; 2x2 plot
+
+		; Add latitude and longitude labels
+
+	txres3                = True
+	txres3@txAngleF       = 90.
+	txres3@txFontHeightF  = 0.02
+	gsn_text_ndc(wks,"Latitude",0.08,0.49,txres3)
+
+	txres3@txAngleF       = 0.
+	gsn_text_ndc(wks,"Longitude",0.5,0.22,txres3)
+
+	frame(wks)
+
+	end do
+	delete(varload)
+ cmd = "convert -delay 25 "+ filename_string+"*.png "+filename_string+".gif"
+ system(cmd)
+  cmd = "rm "+filename_string+"*.png"
+ system(cmd)
+	end do
+
+end
+
+
+
+
+
+
Index: configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-terminator_LOOP.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-terminator_LOOP.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-terminator_LOOP.ncl	(revision 437)
@@ -0,0 +1,231 @@
+;================================================;
+;	Example ncl script to produce the set of
+;	lat-lon plots for DCMIP-2016 test case 1
+; Created by James Kent, Christiane Jablonowski
+;	and Paul Ullrich (University of Michigan) for DCMIP-2012 
+;
+; Modified by Colin Zarzycki for DCMIP-2016
+;
+; Modified by Kevin Viner for terminator ouput
+; Further modified by Maximo Menchaca for creating gifs
+;
+; v1.01 - 6/7/16 
+; v1.02 - 6/8/16 - Typo in var_choice options
+;================================================;
+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/csm/contributed.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
+; ================================================;
+
+begin
+
+;=================================================;
+; open file and read in data 
+; GIVEN AS (time,lev,lat,lon) from 0 to n-1
+;=================================================;
+
+  ; NOTE, that this file needs to be a catted file including all times at daily interval
+  ; if you are outputting 1 file per time, you can run "ncrcat dcmip1.files.*.nc cat.nc"
+  ; for one file
+
+	f = addfile("output_dcmip2016_regular.nc","r")
+
+use_pdf = False
+
+		; Input useful parameters
+
+	lat  = f->lat
+        lon  = f->lon
+	lev  = f->lev
+	nlat   = getfilevardimsizes(f, "lat" )
+	nlon   = getfilevardimsizes(f, "lon" )
+	nlev   = getfilevardimsizes(f, "lev" )
+
+  ; We want plots of temperature, vertical velocity and relative vorticity 
+  ; at the 850 hPa level. We also want surface pressure.
+
+  ; Select var_choice 	1 - average integrated atomic Cl
+  ; 		        2 - average integrated diatomic Cl
+  ;			3 - average integrated Cly
+	do var_choice = 1,3
+
+	if (var_choice .eq. 1) then				; Select Cl
+  		filename_string = "DYNAMICO-case161latlonQ1int"   	; output using X11
+	else if (var_choice .eq. 2) then			; Select Cl2
+		filename_string = "DYNAMICO-case161latlonQ2int"   	; output using png
+	else if (var_choice .eq. 3) then			; Select Cly
+		filename_string = "DYNAMICO-case161latlonQtot_int"
+	end if
+        end if
+        end if
+	
+if(use_pdf) then
+	wks  = gsn_open_wks("pdf",filename_string)   	; output using X11
+else
+	wks  = gsn_open_wks("png",filename_string)   	; output using png
+end if
+
+		; We want the output at days 6, 9, 12 and 15. If your model 
+		; output is in 6 hourly intervals this becomes indices 24, 36, 48 and 60
+
+	output_per_day=8 ; 3-hourly outputs
+	do d = 1,30*output_per_day, (output_per_day/4)
+	day_curr = d/output_per_day
+
+	varload1 = f->Q1_col_int
+	varload2 = f->Q2_col_int
+	if (var_choice .eq. 1) then				; Select Cl
+		var1 = varload1(d,:,:)/4e-6
+  		filename_string = "DYNAMICO-case161latlonQ1int"   	; output using X11
+	else if (var_choice .eq. 2) then			; Select Cl2
+		var1 = varload2(d,:,:)/4e-6
+		filename_string = "DYNAMICO-case161latlonQ2int"   	; output using png
+	else if (var_choice .eq. 3) then				; Select Cly
+		var1 = (varload1(d,:,:)+2*varload2(d,:,:))/4e-6 ; should equal 1.
+		print("Max Cly error :")
+		print(max(abs(var1-1.)))
+		var1 = 1e5*abs(var1-1.)
+		filename_string = "DYNAMICO-case161latlonQtot_int"
+	end if
+        end if
+        end if
+
+
+	var1@long_name = " "
+;	var2@long_name = " "
+;	var3@long_name = " "
+;	var4@long_name = " "
+
+		; Delete loaded data
+
+;	delete(varload)
+
+		; We now produce the plot
+
+;	plot  = new (4, graphic)               			; define plot - need 4 panels
+
+	res1                      = True
+	res1@gsnDraw              = False          		; panel plot 
+	res1@gsnFrame             = False			; don't draw yet
+	res1@cnFillOn             = True
+	res1@cnLinesOn            = True
+	res1@gsnSpreadColors      = True
+	res1@lbLabelAutoStride    = True
+	res1@gsnCenterString      = ""
+	res1@tiMainString         = ""
+	res1@vpWidthF             = 0.38
+	res1@vpHeightF            = 0.19
+	res1@cnLevelSelectionMode = "ManualLevels"
+	res1@cnInfoLabelOn        = False			; don't give each
+	res1@cnLineLabelsOn       = False			; panel its own
+	res1@lbLabelBarOn         = False			; label bar
+
+	pres                          = True
+	pres@gsnMaximize              = True  
+	pres@gsnPanelLabelBar         = True			; Communal label bar
+	pres@gsnPanelLeft             = 0.1
+	pres@gsnPanelRight            = 0.9
+	pres@pmLabelBarOrthogonalPosF = -0.03
+	pres@gsnFrame                 = False
+	pres@lbLabelStride            = 1
+
+	res1@sfXArray        = lon            			; uses lon as plot x-axis
+	res1@sfYArray        = lat            			; uses lat for y axis
+	res1@trYReverse      = False          			; reverses y-axis, false
+	res1@tiYAxisString   = ""   	       			; y-axis title
+	res1@tiXAxisString   = ""             			; x-axis title
+
+	if (var_choice .eq. 1) then 	;================ Temperature plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, Integrated Average Cl" 			; Plot title if required
+	pres@txString = pltTitle
+
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 1.2				; max contour color label
+	res1@cnMinLevelValF  = -0.1				; min contour color label
+	res1@cnLevelSpacingF = 0.1				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 2) then 	;================ Large-scale precipitation plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, Integrated Average Cl2" 	; Plot title if required
+	pres@txString = pltTitle
+
+	var1@units = " "
+;	var2@units = " "
+;	var3@units = " "
+;	var4@units = " "
+
+	res1@cnLinesOn            = False				; Turn lines off for clearer plot
+
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+		; Note that the maximum might be larger than these contour spacing values
+
+	res1@cnMaxLevelValF  = 0.48			; max contour color label
+	res1@cnMinLevelValF  = -0.04				; min contour color label
+	res1@cnLevelSpacingF = 0.04				; contour color spacing			; choose a colormap
+ 
+	else if (var_choice .eq. 3) then 	;================ Surface Pressure plot ====================
+ 
+	pltTitle="DYNAMICO Case 161, Integrated Average Cly" 			; Plot title if required
+	pres@txString = pltTitle
+
+
+	var1@units = " "
+;	var2@units = " "
+;	var3@units = " "
+;	var4@units = " "
+		; Change the output type and name
+
+	gsn_define_colormap(wks,"gui_default")  
+
+	res1@cnMaxLevelValF  = 1.24				; max contour color label
+	res1@cnMinLevelValF  = 0.78				; min contour color label
+	res1@cnLevelSpacingF = 0.04				; contour color spacing			; choose a colormap
+
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+	end if				;============== end plot choice ===============
+
+	res1@gsnCenterString = "t = " + day_curr + " days"
+	plot = gsn_csm_contour(wks,var1(:,:),res1)		; plot var1
+;
+;	res1@gsnCenterString = "t = 9"
+;     	plot(1) = gsn_csm_contour(wks,var2(:,:),res1)		; plot var2
+;
+;	res1@gsnCenterString = "t = 12"
+;	plot(2) = gsn_csm_contour(wks,var3(:,:),res1)		; plot var3
+;
+;	res1@gsnCenterString = "t = 15"
+;     	plot(3) = gsn_csm_contour(wks,var4(:,:),res1)		; plot var4
+
+	gsn_panel(wks,plot,(/1,1/),pres)			; 2x2 plot
+
+		; Add latitude and longitude labels
+
+	txres3                = True
+	txres3@txAngleF       = 90.
+	txres3@txFontHeightF  = 0.02
+	gsn_text_ndc(wks,"Latitude",0.08,0.52,txres3)
+
+	txres3@txAngleF       = 0.
+	gsn_text_ndc(wks,"Longitude",0.54,0.28,txres3)
+
+;	draw(plot)
+	frame(wks)
+
+	end do
+;	delete(varload)
+ cmd = "convert -delay 25 "+ filename_string+"*.png "+filename_string+".gif"
+ system(cmd)
+  cmd = "rm "+filename_string+"*.png"
+ system(cmd)
+	end do
+
+end
Index: configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-horiz-crossx.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-horiz-crossx.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-horiz-crossx.ncl	(revision 437)
@@ -0,0 +1,303 @@
+;=======================================================================================
+; This NCL code calculates radially-averaged tangential and radial wind components
+; as well as T anomaly for DCMIP test case #2 (cyclone)
+; this code requires the accompanying function set "radialAvg.ncl"
+;
+; Usage: User should modify "user options" for their particular data set. Currently,
+; U, V, T, PS are required as variables.
+; If variables are on constant Z surfaces, life is easy.
+;
+; Grepping for "MODELSPEC" will point to possible areas of the code in need of modification
+; for model specific output
+;
+; Written by Colin Zarzycki (zarzycki@ucar.edu)
+; Version 0.1 (6/5/2016) - DCMIP-2016 release
+;=======================================================================================
+
+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/csm/contributed.ncl"  
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"  
+
+begin
+
+;=======================================================================================
+; User options
+;=======================================================================================
+;filename = "/glade/scratch/zarzycki/cam5_work/dcmip2/kessler_bryan/dcmip2.cam.h0.2000-01-05-00000.nc"
+
+filename="output_dcmip2016_regular.nc"
+data_on_constant_Z=False     ; is data already on CONSTANT Z surfaces?
+convert_hybridP_to_Z=False   ; is data on hybrid pressure levels?
+hasTimeIx=True               ; does file have time index?
+timeStep=79                   ; If yes, what time index do you want to plot?
+model="DYNAMICO"              ; used for mainStr, but also for model specific if statements
+mainStr=model+" Day 8"       ; main string for plot titles
+out_type="png"                ; output format, popular options are x11, png, pdf, eps
+
+Uname="U"                    ; Variable name of zonal wind (typically "U" or "ua")
+Vname="V"                    ; Variable name of meridional wind (typically "V" or "ua")
+PSname="PS"                  ; Variable name of surface pressure (typically "PS","PSL","ps",or "slp")
+Tname="T"                    ; Variable name of air temperature (typically "T" or "ta")
+Zname="Z"                    ; Variable name of height (typically "Z" or "za")
+PRECTname="PRECL"            ; Variable name of height (typically "Z" or "za")
+
+;=======================================================================================
+
+;=======================================================================================
+; Other settings, required to remain constant for intercomparison
+; Generally, don't touch unless experimenting/debugging
+;=======================================================================================
+stride=1                       ; subset in horizontal, leave @ 1 unless debugging
+minLat=0.                      ; max lat to read in (deg)
+maxLat=80.                     ; min lat to read in (deg)
+thetaWindMax=60.0              ; max theta wind to plot (m/s)
+radMinMax=18.0                 ; min/max for radial plot (m/s)
+tAnomMax=15.0                  ; max for tAnom plots (K)
+offAnomDist=5.0                ; offset distance for anom calcs (deg)
+zintmin=20.0                   ; min height for z interpolation (m)
+zintmax=18000.0                ; max height for z interpolation (m)
+nzint=100                      ; num points for z interpolation
+Zmidpt=2500.0                  ; half of Z points BELOW this height, half ABOVE
+;=======================================================================================
+
+;=======================================================================================
+; Input data checks
+;=======================================================================================
+if (data_on_constant_Z .and. convert_hybridP_to_Z) then
+  print("Both data_on_constant_Z and convert_hybridP_to_Z cannot be true, exiting...")
+  exit
+end if
+
+;=======================================================================================
+; Get file, coordinate variables, and U, V, PSL
+;=======================================================================================
+print("Loading data from file...")
+
+f = addfile(filename,"r")
+
+lat = f->lat({minLat:maxLat:stride})
+lon = f->lon(::stride)
+lev = f->lev(:)
+nlat = dimsizes(lat)
+nlon = dimsizes(lon)
+nlev = dimsizes(lev)
+
+if (hasTimeIx) then
+  U =  f->$Uname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  V =  f->$Vname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  T =  f->$Tname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  PS = f->$PSname$(timeStep,{minLat:maxLat:stride},::stride)
+  PRECT = f->$PRECTname$(timeStep,{minLat:maxLat:stride},::stride)
+else
+  U =  f->$Uname$(:,{minLat:maxLat:stride},::stride)
+  V = f->$Vname$(:,{minLat:maxLat:stride},::stride)
+  T = f->$Tname$(:,{minLat:maxLat:stride},::stride)
+  PS = f->$PSname$({minLat:maxLat:stride},::stride)
+  PRECT = f->$PRECTname$({minLat:maxLat:stride},::stride)
+end if
+; If U and V are not m/s, please convert here
+U@units="m/s"
+V@units="m/s"
+T@units="K"
+T@long_name="Temperature"
+;=======================================================================================
+
+;=======================================================================================
+; Find center of storm by minimizing PS
+;=======================================================================================
+print("Finding minimum PS location...")
+a = new((/nlat,nlon/),typeof(PS))
+a(:,:) = PS(:,:)
+a1D      = ndtooned(a)
+dsizes_a = dimsizes(a)
+a_indices  = ind_resolve(minind(a1D),dsizes_a) ; Resolve 1D indices to original array
+psminlat = lat(a_indices(0,0))
+psminlon = lon(a_indices(0,1))
+delete([/a,a1D,dsizes_a,a_indices/])
+print("... PS minimum found at lat: "+psminlat+" lon: "+psminlon)
+;=======================================================================================
+
+
+;=======================================================================================
+; Calculate temperature anomaly
+;=======================================================================================
+Tanom = T
+Tref = T(:,{psminlat+offAnomDist},{psminlon+offAnomDist})
+Tanom = T - conform(T,Tref,0)
+Tanom@long_name="Temperature anomaly"
+;=======================================================================================
+
+
+;=======================================================================================
+; Do Z calculations/interpolations if necessary
+;=======================================================================================
+if (.not. data_on_constant_Z) then
+  ;=======================================================================================
+  ; Convert from hybrid levels to Z (CAM)
+  ;=======================================================================================
+  if (convert_hybridP_to_Z)
+    print("Converting from hybrid P to Z...")
+    ; If need_to_convert_P_to_Z is true, variables T,Q,hyam,hybm,hyai,hybm,P0 needed
+    ; all hybrid coefs and TV need to be TOP TO BOTTOM!
+
+    hyam=f->hyam
+    hybm=f->hybm
+    hyai=f->hyai
+    hybi=f->hybi
+    P0=f->P0
+    Tconv = f->T(timeStep,:,{minLat:maxLat:stride},::stride)
+    Qconv = f->Q(timeStep,:,{minLat:maxLat:stride},::stride)
+
+    ; Calculate TV from T (K) and Q (kg/kg)
+    TV=Tconv*(1.+0.61*Qconv)
+
+    ; PHIS is nlatxnlon array = 0.0 for DCMIP Test 2
+    PHIS=PS
+    PHIS=0.0
+
+    Z = cz2ccm(PS,PHIS,TV,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1))
+    Z@units="m"
+
+    delete([/Tconv,Qconv,TV,PHIS/])
+  end if
+
+  ;=======================================================================================
+  ; Example of calculating Z at each level
+  ; MODELSPEC: (if you do not put output directly on either Z levels OR hybrid levels, you will
+  ; need to modify the code here (Tempest shown as example).
+  ;=======================================================================================
+  if (model .eq. "tempest") then
+    Zs=f->Zs
+    Ztop=30000.0
+    Z = U
+    do kk = 0,nlev-1
+      do ii = 0,nlat-1
+        do jj = 0,nlon-1
+          Z(kk,ii,jj) = Zs(ii,jj) + lev(kk)*(Ztop - Zs(ii,jj))
+        end do
+      end do
+    end do
+    Z@units="m"
+  end if
+
+  if (model .eq. "DYNAMICO") then
+    print("DYNAMICO : computing Z from geopotential ...")
+    ; NB : geopotential is defined at level interfaces, average from interfaces to full levels
+    Phi=f->PHI
+    Z=U ; create Z with correct shape
+    do kk = 0,nlev-1
+      do ii = 0,nlat-1
+        do jj = 0,nlon-1
+          Z(kk,ii,jj) = (.5/9.81)*(Phi(timeStep,kk,ii,jj)+Phi(timeStep,kk+1,ii,jj))
+        end do
+      end do
+    end do
+    Z@units="m"
+    print("...done")
+  end if
+
+
+  ;=======================================================================================
+  ; If all else fails, try to load Z directly from file
+  ;=======================================================================================
+  if(.not. isdefined("Z")) then
+    if (isfilevar(f, Zname)) then
+      print("Found Z on file...")
+      Z = f->Z
+    else
+      print("FATAL: Z needs to be either loaded from a file or otherwise defined before continuing...")
+      exit
+    end if
+  end if
+  ;=======================================================================================
+
+end if
+
+;=======================================================================================
+; Generate Z levels of interest
+;=======================================================================================
+print("Generating constant Z levels to interpolate to")
+ZlevPBL = fspan(zintmin,Zmidpt,nzint/2)
+ZlevABL = fspan(Zmidpt,zintmax,nzint/2)
+
+Zlev=(/100.,1000.,1500.,2500.,5000.,10000.,15000./)
+Zlev@units = "m"
+Zlev!0     = "lev"
+Zlev&lev = Zlev
+
+;=======================================================================================
+; Interpolate lat/lon variables to constant Z levels
+;=======================================================================================
+print("Interpolating to Z surfaces")
+U_Z     = int2p_n_Wrap(Z,U,Zlev,2,0)
+V_Z     = int2p_n_Wrap(Z,V,Zlev,2,0)
+Tanom_Z = int2p_n_Wrap(Z,Tanom,Zlev,2,0)
+;=======================================================================================
+
+WIND_Z = U_Z
+WIND_Z = sqrt(U_Z^2+V_Z^2)
+WIND_Z@long_name="Horizontal wind"
+
+PRECT = PRECT*8.64e7   ; convert m/s to mm/day
+PRECT@long_name="Precipitation rate"
+PRECT@units="mm/d"
+
+print("Plotting...")
+
+wks   = gsn_open_wks (out_type,"horiz_crossx")               ; send graphics to PNG file
+contour = new(4,"graphic")
+
+gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
+
+res  = True
+res@gsnDraw = False
+res@gsnFrame = False
+res@gsnSpreadColors  = True        ; Span full color map
+res@cnFillOn         = True        ; Turn on contour fill
+res@cnLinesOn        = False
+res@cnLineLabelsOn   = False
+res@cnInfoLabelOn    = False
+res@gsnAddCyclic = True
+
+latWidth = 10.0
+lonWidth = 10.0
+res@mpOutlineOn = False
+res@mpMaxLatF = psminlat + latWidth
+res@mpMinLatF = psminlat - latWidth
+res@mpMinLonF = psminlon - lonWidth
+res@mpMaxLonF = psminlon + lonWidth
+
+res@cnLevelSelectionMode = "ManualLevels"
+res@cnLevelSpacingF      =  5.0
+res@cnMinLevelValF       =  10.0
+res@cnMaxLevelValF       =  80.0
+res@tiMainString="1500m Horiz. Wind"
+contour(0) = gsn_csm_contour_map(wks,WIND_Z({1500.0},:,:),res)  ; create the plot
+
+res@cnLevelSelectionMode = "ManualLevels"
+res@cnLevelSpacingF      =  10.0
+res@cnMinLevelValF       =  -60.0
+res@cnMaxLevelValF       =  60.0
+res@tiMainString="100m U-Wind"
+contour(1) = gsn_csm_contour_map(wks,U_Z({100.0},:,:),res)  ; create the plot
+
+res@cnLevelSelectionMode = "ManualLevels"
+res@cnLevelSpacingF      =  1.0
+res@cnMinLevelValF       =  0.0
+res@cnMaxLevelValF       =  10.0
+res@tiMainString="5000m T anom."
+contour(2) = gsn_csm_contour_map(wks,Tanom_Z({5000.0},:,:),res)  ; create the plot
+
+res@cnLevelSelectionMode = "ManualLevels"
+res@cnLevelSpacingF      =  200.0
+res@cnMinLevelValF       =  200.0
+res@cnMaxLevelValF       =  4000.0
+res@tiMainString="Precip. rate"
+contour(3) = gsn_csm_contour_map(wks,PRECT(:,:),res)  ; create the plot
+
+resP                     = True                ; modify the panel plot
+
+gsn_panel(wks,contour,(/2,2/),resP)             ; now draw as one plot
+
+end
Index: configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-intensity-trace.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-intensity-trace.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-intensity-trace.ncl	(revision 437)
@@ -0,0 +1,82 @@
+;=======================================================================================
+; This NCL 
+;
+; Written by Colin Zarzycki (zarzycki@ucar.edu)
+; Version 0.1 (6/5/2016) - DCMIP-2016 release
+;=======================================================================================
+
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
+begin
+
+;=======================================================================================
+
+
+;whichVar="PS"          ; PS or WIND
+top_to_bottom=True     ; is data top to bottom? (false for bottom to top)
+model="DYNAMICO"
+out_type="png"
+
+; needs to be a catted file ordered time, lev, lat, lon (or time, lat, lon for PS)
+filename="output_dcmip2016_regular.nc"
+;=======================================================================================
+
+maxLat=90.0
+minLat=0.0
+
+;=======================================================================================
+
+a = addfile(filename,"r")
+
+; figure out what is surface level for wind speed...
+if (top_to_bottom) then
+  sfcLev=dimsizes(a->lev)-1                  ; dimsizes(a->lev)-1 subsets bottom level
+else
+  sfcLev=0
+end if
+
+do varid=1,2
+    
+if( varid .eq. 1) then
+  whichVar = "PS"
+  PS=a->PS(:,{minLat:maxLat},:)              ; just take out NH since we know storm is there
+  PS=PS/100.                                 ; convert from Pa to hPa/mb
+  PS@long_name="Surface pressure"
+  PS@units="hPa"
+  var_vtime = dim_min_n_Wrap(PS,(/1,2/))     ; for now, just find min over lat/lon domain at each time
+else if (varid .eq. 2) then
+  whichVar = "WIND"
+  UBOT=a->U(:,sfcLev,{minLat:maxLat},:)
+  VBOT=a->V(:,sfcLev,{minLat:maxLat},:)
+  WIND=UBOT
+  WIND=sqrt(UBOT^2+VBOT^2)
+  WIND@long_name="Lowest level wind"
+  WIND@units="m/s"
+  var_vtime = dim_max_n_Wrap(WIND,(/1,2/))
+else
+  print("unsupported var choice... exiting...")
+  exit
+end if
+end if
+
+wks   = gsn_open_wks (out_type,whichVar)               ; send graphics to PNG file
+
+res                  = True                     ; plot mods desired
+res@gsnDraw=False
+res@gsnFrame=False
+res@tiMainString     = "DCMIP2 "+model+" "+whichVar+" v. time"          ; add title
+
+res@tiXAxisString        = "Time (days)"  
+res@tiYAxisString        = var_vtime@long_name+" ("+var_vtime@units+")"
+
+res@xyLineThicknessF=5.0
+res@xyLineColor="Blue"
+
+time = var_vtime&time_counter/86400
+plot  = gsn_csm_xy (wks,time,var_vtime,res) ; create plot
+draw(plot)
+frame(wks)
+
+end do
+
+end
Index: configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-radial-avgs.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-radial-avgs.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/CYCLONE/dcmip2-radial-avgs.ncl	(revision 437)
@@ -0,0 +1,317 @@
+;=======================================================================================
+; This NCL code calculates radially-averaged tangential and radial wind components
+; as well as T anomaly for DCMIP test case #2 (cyclone)
+; this code requires the accompanying function set "radialAvg.ncl"
+;
+; Usage: User should modify "user options" for their particular data set. Currently,
+; U, V, T, PS are required as variables.
+; If variables are on constant Z surfaces, life is easy.
+;
+; Grepping for "MODELSPEC" will point to possible areas of the code in need of modification
+; for model specific output
+;
+; Written by Colin Zarzycki (zarzycki@ucar.edu)
+; Version 0.1 (6/5/2016) - DCMIP-2016 release
+;=======================================================================================
+
+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/csm/contributed.ncl"  
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"  
+load "radialAvg.ncl"
+
+begin
+
+;=======================================================================================
+; User options
+;=======================================================================================
+filename = "output_dcmip2016_regular.nc"
+data_on_constant_Z=False     ; is data already on CONSTANT Z surfaces?
+convert_hybridP_to_Z=False   ; is data on hybrid pressure levels?
+hasTimeIx=True               ; does file have time index?
+timeStep=64                   ; If yes, what time index do you want to plot?
+model="DYNAMICO"              ; used for mainStr, but also for model specific if statements
+mainStr=model+" Day 8"       ; main string for plot titles
+outType="png"                ; output format, popular options are x11, png, pdf, eps
+
+Uname="U"                    ; Variable name of zonal wind (typically "U" or "ua")
+Vname="V"                    ; Variable name of meridional wind (typically "V" or "ua")
+Rhoname="Rho"                  ; Variable name of Rho
+PSname="PS"                  ; Variable name of surface pressure (typically "PS","PSL","ps",or "slp")
+Tname="T"                    ; Variable name of air temperature (typically "T" or "ta")
+Zname="Z"                    ; Variable name of height (typically "Z" or "za")
+;=======================================================================================
+
+;=======================================================================================
+; Other settings, required to remain constant for intercomparison
+; Generally, don't touch unless experimenting/debugging
+;=======================================================================================
+stride=1                       ; subset in horizontal, leave @ 1 unless debugging
+minLat=0.                      ; max lat to read in (deg)
+maxLat=80.                     ; min lat to read in (deg)
+thetaWindMax=60.0              ; max theta wind to plot (m/s)
+radMinMax=18.0                 ; min/max for radial plot (m/s)
+tAnomMax=15.0                  ; max for tAnom plots (K)
+offAnomDist=5.0                ; offset distance for anom calcs (deg)
+zintmin=20.0                   ; min height for z interpolation (m)
+zintmax=18000.0                ; max height for z interpolation (m)
+nzint=100                      ; num points for z interpolation
+Zmidpt=2500.0                  ; half of Z points BELOW this height, half ABOVE
+;=======================================================================================
+
+;=======================================================================================
+; Input data checks
+;=======================================================================================
+if (data_on_constant_Z .and. convert_hybridP_to_Z) then
+  print("Both data_on_constant_Z and convert_hybridP_to_Z cannot be true, exiting...")
+  exit
+end if
+
+;=======================================================================================
+; Get file, coordinate variables, and U, V, PSL
+;=======================================================================================
+print("Loading data from file...")
+
+f = addfile(filename,"r")
+
+lat = f->lat({minLat:maxLat:stride})
+lon = f->lon(::stride)
+lev = f->lev(:)
+nlat = dimsizes(lat)
+nlon = dimsizes(lon)
+nlev = dimsizes(lev)
+
+if (hasTimeIx) then
+  U =  f->$Uname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  V =  f->$Vname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  T =  f->$Tname$(timeStep,:,{minLat:maxLat:stride},::stride)
+  if (model .eq. "DYNAMICO") then
+    P0 = f->P(timeStep,0,{minLat:maxLat:stride},::stride)
+    T0 =  f->$Tname$(timeStep,0,{minLat:maxLat:stride},::stride)
+    Rho = P0/(287*T0) ;  not quite exact (FIXME)
+    Rho@units="kg/m3"
+  else
+    Rho = f->$Rhoname$(timeStep,0,{minLat:maxLat:stride},::stride)
+  end if
+  PS = f->$PSname$(timeStep,{minLat:maxLat:stride},::stride)
+else
+  U =  f->$Uname$(:,{minLat:maxLat:stride},::stride)
+  V = f->$Vname$(:,{minLat:maxLat:stride},::stride)
+  T = f->$Tname$(:,{minLat:maxLat:stride},::stride)
+  Rho = f->$Rhoname$(0,{minLat:maxLat:stride},::stride)
+  ;PS = f->$PSname$({minLat:maxLat:stride},::stride)
+end if
+; If U and V are not m/s, please convert here
+U@units="m/s"
+V@units="m/s"
+T@units="K"
+T@long_name="Temperature"
+;=======================================================================================
+
+;=======================================================================================
+; Find center of storm by minimizing PS
+;=======================================================================================
+print("Finding minimum PS location...")
+a = new((/nlat,nlon/),typeof(PS))
+a(:,:) = PS(:,:)
+a1D      = ndtooned(a)
+dsizes_a = dimsizes(a)
+a_indices  = ind_resolve(minind(a1D),dsizes_a) ; Resolve 1D indices to original array
+psminlat = lat(a_indices(0,0))
+psminlon = lon(a_indices(0,1))
+delete([/a,a1D,dsizes_a,a_indices/])
+print("... minimum found at lat: "+psminlat+" lon: "+psminlon)
+;=======================================================================================
+
+
+;=======================================================================================
+; Calculate temperature anomaly
+;=======================================================================================
+Tanom = T
+Tref = T(:,{psminlat+offAnomDist},{psminlon+offAnomDist})
+Tanom = T - conform(T,Tref,0)
+Tanom@long_name="Temperature anomaly"
+;=======================================================================================
+
+
+;=======================================================================================
+; Do Z calculations/interpolations if necessary
+;=======================================================================================
+if (.not. data_on_constant_Z) then
+  ;=======================================================================================
+  ; Convert from hybrid levels to Z (CAM)
+  ;=======================================================================================
+  if (convert_hybridP_to_Z)
+    print("Converting from hybrid P to Z...")
+    ; If need_to_convert_P_to_Z is true, variables T,Q,hyam,hybm,hyai,hybm,P0 needed
+    ; all hybrid coefs and TV need to be TOP TO BOTTOM!
+
+    hyam=f->hyam
+    hybm=f->hybm
+    hyai=f->hyai
+    hybi=f->hybi
+    P0=f->P0
+    Tconv = f->T(timeStep,:,{minLat:maxLat:stride},::stride)
+    Qconv = f->Q(timeStep,:,{minLat:maxLat:stride},::stride)
+
+    ; Calculate TV from T (K) and Q (kg/kg)
+    TV=Tconv*(1.+0.61*Qconv)
+
+    ; PHIS is nlatxnlon array = 0.0 for DCMIP Test 2
+    PHIS=Rho(0,:,:)
+    PHIS=0.0
+
+    Z = cz2ccm(Rho,PHIS,TV,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1))
+    Z@units="m"
+
+    delete([/Tconv,Qconv,TV,PHIS/])
+  end if
+
+  ;=======================================================================================
+  ; Example of calculating Z at each level
+  ; MODELSPEC: (if you do not put output directly on either Z levels OR hybrid levels, you will
+  ; need to modify the code here (Tempest shown as example).
+  ;=======================================================================================
+  if (model .eq. "DYNAMICO") then
+    print("DYNAMICO : computing Z from geopotential ...")
+    ; NB : geopotential is defined at level interfaces, average from interfaces to full levels
+    Phi=f->PHI
+    Z=U ; create Z with correct shape
+    do kk = 0,nlev-1
+      do ii = 0,nlat-1
+        do jj = 0,nlon-1
+          Z(kk,ii,jj) = (.5/9.81)*(Phi(timeStep,kk,ii,jj)+Phi(timeStep,kk+1,ii,jj))
+        end do
+      end do
+    end do
+    Z@units="m"
+    print("...done")
+  end if
+
+  ;=======================================================================================
+  ; If all else fails, try to load Z directly from file
+  ;=======================================================================================
+  if(.not. isdefined("Z")) then
+    if (isfilevar(f, Zname)) then
+      print("Found Z on file...")
+      Z = f->Z
+    else
+      print("FATAL: Z needs to be either loaded from a file or otherwise defined before continuing...")
+      exit
+    end if
+  end if
+
+  ;=======================================================================================
+  ; Generate uniform Z levs to interpolate to
+  ; Extra density (half of pts) is specified in PBL, controlled by Zmidpt
+  ;=======================================================================================
+  print("Generating constant Z levels to interpolate to")
+  ZlevPBL = fspan(zintmin,Zmidpt,nzint/2)
+  ZlevABL = fspan(Zmidpt,zintmax,nzint/2)
+
+  Zlev=new(nzint-1,"float")
+  Zlev(0:(nzint/2)-1)=ZlevPBL
+  Zlev((nzint/2)-1:nzint-2)=ZlevABL
+  Zlev@units = "m"
+  Zlev!0     = "lev"
+  Zlev&lev = Zlev
+  delete([/ZlevPBL,ZlevABL/])
+  ;=======================================================================================
+
+  ;=======================================================================================
+  ; Interpolate lat/lon variables to constant Z levels
+  ;=======================================================================================
+  print("Interpolating from Z model levels to constant Z surfaces")
+  U_Z     = int2p_n_Wrap(Z,U,Zlev,2,0)
+  V_Z     = int2p_n_Wrap(Z,V,Zlev,2,0)
+  Tanom_Z = int2p_n_Wrap(Z,Tanom,Zlev,2,0)
+  ;=======================================================================================
+else
+  print("Data already on Z surfaces, copying to new vars")
+  U_Z     = U
+  V_Z     = V
+  Tanom_Z = Tanom
+end if
+
+;=======================================================================================
+; Calculate radial and tangential wind components on native grid 
+;=======================================================================================
+print("Separating winds into radial and tangential components given Rho min loc...")
+vComps = calcWindComponents(U_Z,V_Z,lat,lon,psminlat,psminlon)
+v_rad = vComps[0]
+v_theta = vComps[1]
+delete(vComps)
+;=======================================================================================
+
+
+;=======================================================================================
+; Perform radial integrations
+;=======================================================================================
+print("Do radial integrations...")
+rad_v_theta  = radialAvg3D(v_theta,lat,lon,v_theta&lev,psminlat,psminlon,800.,False)
+rad_v_rad    = radialAvg3D(v_rad  ,lat,lon,v_rad&lev  ,psminlat,psminlon,800.,False)
+rad_t_anom   = radialAvg3D(Tanom_Z,lat,lon,Tanom_Z&lev  ,psminlat,psminlon,800.,False)
+;=======================================================================================
+
+
+print(v_theta&lev)
+print(rad_v_theta&radius)
+
+print("Plotting...")
+
+res  = True
+res@gsnDraw = False
+res@gsnFrame = False
+res@gsnSpreadColors  = True        ; Span full color map
+res@cnFillOn         = True        ; Turn on contour fill
+res@cnLinesOn        = False
+res@cnLineLabelsOn   = False
+res@cnInfoLabelOn    = False
+res@tiYAxisString    = "Height (m)"
+res@gsnYAxisIrregular2Linear = True
+res@gsnXAxisIrregular2Linear = True
+
+wks   = gsn_open_wks (outType,"t_anom_"+model)
+gsn_define_colormap(wks,"matlab_jet")
+res_rad = res
+res_rad@cnLevelSelectionMode = "ExplicitLevels"
+res_rad@cnLevels = ispan(0,toint(tAnomMax),1)
+plot = gsn_csm_contour(wks,rad_t_anom,res_rad)
+draw(plot)
+frame(wks)
+
+delete(plot)
+delete(wks)
+delete(res_rad)
+
+wks   = gsn_open_wks (outType,"v_theta_"+model)
+gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
+res_rad = res
+res_rad@tiMainString     = mainStr+" tangential wind"
+res_rad@cnLevelSelectionMode = "ExplicitLevels"
+res_rad@cnLevels = fspan(0.0,thetaWindMax,21)
+plot = gsn_csm_contour(wks,rad_v_theta,res_rad)
+draw(plot)
+frame(wks)
+
+delete(plot)
+delete(wks)
+delete(res_rad)
+
+wks   = gsn_open_wks (outType,"v_rad_"+model)
+gsn_define_colormap(wks,"hotcolr_19lev")
+res_rad = res
+res_rad@tiMainString     = mainStr+" radial wind"
+res_rad@cnLevelSelectionMode = "ExplicitLevels"
+res_rad@cnLevels = ispan(toint(-radMinMax),toint(radMinMax),1)
+plot = gsn_csm_contour(wks,rad_v_rad,res_rad)
+draw(plot)
+frame(wks)
+
+delete(plot)
+delete(wks)
+delete(res_rad)
+
+
+
+end
Index: configs/dcmip2016/PLOTS/testcase/CYCLONE/include/radialAvg.ncl
===================================================================
--- configs/dcmip2016/PLOTS/testcase/CYCLONE/include/radialAvg.ncl	(revision 437)
+++ configs/dcmip2016/PLOTS/testcase/CYCLONE/include/radialAvg.ncl	(revision 437)
@@ -0,0 +1,466 @@
+; This file contains various functions relating to TC plotting in radial space
+; Colin Zarzycki (zarzycki@ucar.edu)
+; History:
+; 6/5/2016 - Released for DCMIP-2016
+;
+; contains: radialAvg3D, radialAvg2D, calcWindComponents, radialInt3D, calc_vertmax_from_rad
+;
+;
+; radialAvg3D
+; Usage: radialAvg3D(data,lat,lon,lev,psminlat,psminlon,outerRad,mergeInnerBins)
+; data (numeric): 3-D (lev,lat,lon) of variable to transform to radial space
+; lat (numeric): 1-D array of lats (deg)
+; lon (numeric): 1-D array of lons (deg)
+; lev (numeric): 1-D array of levs (Pa or m)
+; psminlat (numeric): lat of TC center (deg)
+; psminlon (numeric): lon of TC center (deg)
+; outerRad (numeric) : outermost radius to calculate (km) (typically 500-800km)
+; mergeInnerBins (logical) :merge bins at storm center for sampling purpose? (generally True)
+;
+; ex: radialAvg3D(v_theta,lat,lon,v_theta&lev,psminlat,psminlon,500.,True)
+;
+;
+; radialAvg2D
+; Usage: radialAvg2D(data,lat,lon,psminlat,psminlon,outerRad,mergeInnerBins)
+; data (numeric): 2-D (lat,lon) of variable to transform to radial space
+; lat (numeric): 1-D array of lats (deg)
+; lon (numeric): 1-D array of lons (deg)
+; psminlat (numeric): lat of TC center (deg)
+; psminlon (numeric): lon of TC center (deg)
+; outerRad (numeric) : outermost radius to calculate (km) (typically 500-800km)
+; mergeInnerBins (logical) :merge bins at storm center for sampling purpose? (generally True)
+;
+; ex: 
+;
+;
+; calcWindComponents
+; Usage: radialAvg3D(U,V,lat,lon,psminlat,psminlon)
+; U (numeric): 2-D (lat,lon) or 3-D (lev, lat, lon) or 4-D (time, lev, lat, lon) zonal wind
+; V (numeric): 2-D (lat,lon) or 3-D (lev, lat, lon) or 4-D (time, lev, lat, lon) zonal wind
+; lat (numeric): 1-D array of lats (deg)
+; lon (numeric): 1-D array of lons (deg)
+; psminlat (numeric): lat of TC center (deg)
+; psminlon (numeric): lon of TC center (deg)
+;
+; ex: calcWindComponents(U_Z,V_Z,lat,lon,psminlat,psminlon)
+;
+;
+; radialInt3D
+; Usage: radialInt3D(data,lat,lon,psminlat,psminlon,outerRad)
+; data (numeric): 3-D (lev,lat,lon)
+; lat (numeric): 1-D array of lats (deg)
+; lon (numeric): 1-D array of lons (deg)
+; lev (numeric): 1-D array of levs (Pa or m)
+; psminlat (numeric): lat of TC center (deg)
+; psminlon (numeric): lon of TC center (deg)
+; outerRad (numeric) : outermost radius to integrate (km) (typically 500-800km)
+;
+; ex: radialInt3D(tofloat(DTCOND),tofloat(lat),tofloat(lon),tofloat(lev),psminlat,psminlon,radius_3D)
+
+undef ( "radialAvg3D" )
+function radialAvg3D ( \
+  data[*][*][*]    : numeric, \
+  lat[*]           : numeric, \
+  lon[*]           : numeric, \
+  lev[*]           : numeric, \
+  psminlat         : numeric, \
+  psminlon         : numeric, \
+  outerRad         : numeric, \
+  mergeInnerBins   : logical \
+)
+local None
+
+; lat (deg)
+; lon (deg)
+; center_lat(deg)
+; center_lon(deg)
+; outerRad (km)
+
+begin
+
+deltaLat = lat(1)-lat(0)
+deltaLon = lon(1)-lon(0)
+deltaMax = max((/deltaLon, deltaLat/))
+deltaMin = min((/deltaLon, deltaLat/))
+kmInDeg = 111.32
+kmGrid = kmInDeg * deltaMax
+print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
+
+nlat = dimsizes(lat)
+nlon = dimsizes(lon)
+nlev = dimsizes(lev)
+
+pi = 3.14159
+d2r = pi/180.
+r2d = 180./pi
+lonr = lon*d2r
+latr = lat*d2r
+copy_VarCoords(lon,lonr)
+copy_VarCoords(lat,latr)
+
+psminLocY = ind_nearest_coord(psminlat,lat,0)
+psminLocX = ind_nearest_coord(psminlon,lon,0)
+
+;;;;; DO RADIAL AVERAGES
+
+loopWidth = toint(ceil((outerRad*1.10) / (kmInDeg * deltaMin)))
+; for standard lat/lon grid
+;1.25 gives n ~= 5 in inner grid if numMerge = 2
+;1.4 gives n ~= 9 in inner grid if numMerge = 2
+timesGrid = 1.1 ; we want each radius bin to be timeGrid times kmGrid
+
+nx = toint(outerRad/(timesGrid*kmGrid))
+print("Number of bins is equal to "+nx)
+  
+if (mergeInnerBins) then
+;  origRadiusArr = fspan(0,outerRad,nx)
+;  radiusArr = new(dimsizes(origRadiusArr)-1,double)
+;  radiusArr(0) = origRadiusArr(0)
+;  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(2:dimsizes(origRadiusArr)-1)
+  numMerge = 2 ; numMerge is the number of innermost radial bins we want to merge
+  print("Merging innermost "+numMerge+" bins because radial average is so small.")
+  numMergeMinusOne = numMerge-1
+  origRadiusArr = fspan(0,outerRad,nx+numMergeMinusOne)
+  radiusArr = new((dimsizes(origRadiusArr)-numMergeMinusOne),double)
+  radiusArr(0) = origRadiusArr(0)
+  print("Going from "+dimsizes(origRadiusArr)+" to "+dimsizes(radiusArr)+" bins...")
+  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(numMerge:dimsizes(origRadiusArr)-1)
+else
+  print("Not merging any innermost bins -- be careful that your inner bins have > 1 pt.")
+  radiusArr = fspan(0,outerRad,nx)
+end if
+
+numRadBins = dimsizes(radiusArr)
+
+rad_thevar_hit  = new((/nlev,numRadBins/),integer)
+rad_thevar_cum  = new((/nlev,numRadBins/),double)
+
+; NEEDS to be 1 or 0
+rad_thevar_hit!1 = "radius"
+rad_thevar_cum!1 = "radius"
+rad_thevar_hit&radius = radiusArr
+rad_thevar_cum&radius = radiusArr
+
+rad_thevar_hit!0 = "lev"
+rad_thevar_cum!0 = "lev"
+rad_thevar_hit&lev = lev
+rad_thevar_cum&lev = lev
+
+rad_thevar_hit = 0
+rad_thevar_cum = 0
+
+; NEED TO CORRECT NLEV
+print("Starting loop")
+do i = psminLocY-loopWidth,psminLocY+loopWidth
+  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
+  do j = psminLocX-loopWidth,psminLocX+loopWidth
+    do k = 0,dimsizes(lev)-1
+      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
+      if (gcdist .le. outerRad)
+        rad_thevar_hit(k,{gcdist}) = rad_thevar_hit(k,{gcdist}) + 1 
+        rad_thevar_cum(k,{gcdist}) = rad_thevar_cum(k,{gcdist}) + data(k,i,j)
+      end if
+    end do
+  end do
+end do
+
+print("Minimum number of hits per gridbox: "+min(rad_thevar_hit))
+print("Maximum number of hits per gridbox: "+max(rad_thevar_hit))
+rad_thevar = rad_thevar_cum/rad_thevar_hit
+
+copy_VarCoords(rad_thevar_hit,rad_thevar)
+
+rad_thevar@units = data@units
+rad_thevar@long_name=data@long_name
+
+return(rad_thevar)
+
+end
+
+; ====================================================================================================
+
+undef ( "radialAvg2D" )
+function radialAvg2D ( \
+  data[*][*]       : numeric, \
+  lat[*]           : numeric, \
+  lon[*]           : numeric, \
+  psminlat         : numeric, \
+  psminlon         : numeric, \
+  outerRad         : numeric, \
+  mergeInnerBins   : logical \
+)
+local None
+
+; lat (deg)
+; lon (deg)
+; center_lat(deg)
+; center_lon(deg)
+; outerRad (km)
+
+begin
+
+deltaLat = lat(1)-lat(0)
+deltaLon = lon(1)-lon(0)
+deltaMax = max((/deltaLon, deltaLat/))
+deltaMin = min((/deltaLon, deltaLat/))
+kmInDeg = 111.32
+kmGrid = kmInDeg * deltaMax
+print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
+
+nlat = dimsizes(lat)
+nlon = dimsizes(lon)
+
+pi = 3.14159
+d2r = pi/180.
+r2d = 180./pi
+lonr = lon*d2r
+latr = lat*d2r
+copy_VarCoords(lon,lonr)
+copy_VarCoords(lat,latr)
+
+psminLocY = ind_nearest_coord(psminlat,lat,0)
+psminLocX = ind_nearest_coord(psminlon,lon,0)
+
+;;;;; DO RADIAL AVERAGES
+
+loopWidth = toint(ceil((outerRad*1.10) / (kmInDeg * deltaMin)))
+; for standard lat/lon grid
+;1.25 gives n ~= 5 in inner grid if numMerge = 2
+;1.4 gives n ~= 9 in inner grid if numMerge = 2
+timesGrid = 1.1 ; we want each radius bin to be timeGrid times kmGrid
+
+nx = toint(outerRad/(timesGrid*kmGrid))
+print("Number of bins is equal to "+nx)
+  
+if (mergeInnerBins) then
+;  origRadiusArr = fspan(0,outerRad,nx)
+;  radiusArr = new(dimsizes(origRadiusArr)-1,double)
+;  radiusArr(0) = origRadiusArr(0)
+;  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(2:dimsizes(origRadiusArr)-1)
+  numMerge = 2 ; numMerge is the number of innermost radial bins we want to merge
+  print("Merging innermost "+numMerge+" bins because radial average is so small.")
+  numMergeMinusOne = numMerge-1
+  origRadiusArr = fspan(0,outerRad,nx+numMergeMinusOne)
+  radiusArr = new((dimsizes(origRadiusArr)-numMergeMinusOne),double)
+  radiusArr(0) = origRadiusArr(0)
+  print("Going from "+dimsizes(origRadiusArr)+" to "+dimsizes(radiusArr)+" bins...")
+  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(numMerge:dimsizes(origRadiusArr)-1)
+else
+  print("Not merging any innermost bins -- be careful that your inner bins have > 1 pt.")
+  radiusArr = fspan(0,outerRad,nx)
+end if
+
+numRadBins = dimsizes(radiusArr)
+
+rad_thevar_hit  = new((/numRadBins/),integer)
+rad_thevar_cum  = new((/numRadBins/),double)
+
+; NEEDS to be 1 or 0
+rad_thevar_hit!0 = "radius"
+rad_thevar_cum!0 = "radius"
+rad_thevar_hit&radius = radiusArr
+rad_thevar_cum&radius = radiusArr
+
+rad_thevar_hit = 0
+rad_thevar_cum = 0
+
+; NEED TO CORRECT NLEV
+print("Starting loop")
+do i = psminLocY-loopWidth,psminLocY+loopWidth
+  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
+  do j = psminLocX-loopWidth,psminLocX+loopWidth
+      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
+      if (gcdist .le. outerRad)
+        rad_thevar_hit({gcdist}) = rad_thevar_hit({gcdist}) + 1 
+        rad_thevar_cum({gcdist}) = rad_thevar_cum({gcdist}) + data(i,j)
+      end if
+  end do
+end do
+
+print("Minimum number of hits per gridbox: "+min(rad_thevar_hit))
+print("Maximum number of hits per gridbox: "+max(rad_thevar_hit))
+rad_thevar = rad_thevar_cum/rad_thevar_hit
+
+copy_VarCoords(rad_thevar_hit,rad_thevar)
+
+rad_thevar@units = data@units
+rad_thevar@long_name=data@long_name
+
+return(rad_thevar)
+
+end
+
+; ====================================================================================================
+
+undef ( "calcWindComponents" )
+function calcWindComponents ( \
+  U[*][*][*]       : numeric, \
+  V[*][*][*]       : numeric, \
+  lat[*]           : numeric, \
+  lon[*]           : numeric, \
+  psminlat         : numeric, \
+  psminlon         : numeric \
+)
+local None
+begin
+
+pi = 3.14159
+d2r = pi/180.
+r2d = 180./pi
+lonr = lon*d2r
+latr = lat*d2r
+copy_VarCoords(lon,lonr)
+copy_VarCoords(lat,latr)
+
+deltalong = lonr(:)-(psminlon*d2r)
+
+arr1 = conform(U,sin(deltalong),2) * conform(U,cos(latr),1)
+arr2 = cos(psminlat*d2r)*conform(U,sin(latr),1) - sin(psminlat*d2r)*conform(U,cos(latr),1)*conform(U,cos(deltalong),2)
+dir_angle_r = atan2(arr1,arr2)
+; need to conform here
+dir_angle_d = dir_angle_r*r2d
+
+copy_VarCoords(U,dir_angle_r)
+copy_VarCoords(U,dir_angle_d)
+
+eps = 10.e-6
+phi_r = atan2(U,V)
+phi_d = phi_r*r2d
+
+copy_VarCoords(U,phi_r)
+copy_VarCoords(U,phi_d)
+
+WIND = U
+WIND = sqrt(U^2+V^2)
+
+v_theta = WIND*sin(dir_angle_r - phi_r)
+v_rad = -WIND*cos(dir_angle_r - phi_r)
+
+copy_VarCoords(U,v_theta)
+copy_VarCoords(U,v_rad)
+
+v_theta@units = U@units
+v_theta@long_name="Tangential wind"
+v_rad@units = U@units
+v_rad@long_name="Radial wind"
+
+return( [/ v_rad, v_theta /] )
+
+end
+
+
+undef ( "radialInt3D" )
+function radialInt3D ( \
+  data[*][*][*]    : numeric, \
+  lat[*]           : numeric, \
+  lon[*]           : numeric, \
+  lev[*]           : numeric, \
+  psminlat         : numeric, \
+  psminlon         : numeric, \
+  outerRad         : numeric \
+)
+local None
+
+; lat (deg)
+; lon (deg)
+; center_lat(deg)
+; center_lon(deg)
+; outerRad (km)
+
+begin
+
+deltaLat = lat(1)-lat(0)
+deltaLon = lon(1)-lon(0)
+deltaMax = max((/deltaLon, deltaLat/))
+deltaMin = min((/deltaLon, deltaLat/))
+kmInDeg = 111.32
+kmGrid = kmInDeg * deltaMax
+print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
+
+nlat = dimsizes(lat)
+nlon = dimsizes(lon)
+nlev = dimsizes(lev)
+
+pi = 3.14159
+d2r = pi/180.
+r2d = 180./pi
+lonr = lon*d2r
+latr = lat*d2r
+copy_VarCoords(lon,lonr)
+copy_VarCoords(lat,latr)
+
+psminLocY = ind_nearest_coord(psminlat,lat,0)
+psminLocX = ind_nearest_coord(psminlon,lon,0)
+
+; Get cosine weights
+rad    = 4.0*atan(1.0)/180.0
+clat   = cos(lat*rad)
+
+loopWidth = toint(ceil((outerRad*1.25) / (kmInDeg * deltaMin)))
+; make sure we are more than 1 to do a few rows of buffer cells
+
+rad_thevar_hit  = new((/nlev/),float)
+rad_thevar_cum  = new((/nlev/),float)
+
+rad_thevar_hit!0 = "lev"
+rad_thevar_cum!0 = "lev"
+rad_thevar_hit&lev = lev
+rad_thevar_cum&lev = lev
+
+rad_thevar_hit = 0
+rad_thevar_cum = 0
+
+; NEED TO CORRECT NLEV
+print("Starting loop")
+do i = psminLocY-loopWidth,psminLocY+loopWidth
+  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
+  do j = psminLocX-loopWidth,psminLocX+loopWidth
+    do k = 0,dimsizes(lev)-1
+      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
+      if (gcdist .le. outerRad)
+        rad_thevar_hit(k) = rad_thevar_hit(k) + data(k,i,j)*clat(i)
+        rad_thevar_cum(k) = rad_thevar_cum(k) + clat(i)
+      end if
+    end do
+  end do
+end do
+
+rad_thevar = rad_thevar_hit/rad_thevar_cum
+
+copy_VarCoords(rad_thevar_hit,rad_thevar)
+
+rad_thevar@units = data@units
+rad_thevar@long_name=data@long_name
+
+return(rad_thevar)
+
+end
+
+
+
+undef ( "calc_vertmax_from_rad" )
+function calc_vertmax_from_rad ( \
+  vrad[*][*]       : numeric \
+)
+local None
+
+begin
+
+lev = vrad&lev
+nlev = dimsizes(lev)
+rad = vrad&radius
+nrad = dimsizes(rad)
+
+rmw = new((/nlev/),"float")
+
+do k = 0,dimsizes(lev)-1
+  maxix = maxind(vrad(k,:))
+  rmw(k) = tofloat(rad(maxix))
+end do
+
+rmw!0 = "lev"
+rmw&lev = lev
+
+return(rmw)
+
+end
