;================================================;
;	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