Home

Context Navigation

spectra.ncl @ 560

Last change on this file since 560 was 534, checked in by heinze, 15 years ago
Bugfix in spectra.ncl concerning output in png and small changes in layout
File size: 21.5 KB

Rev	Line
[176]	1	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
	2	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
	3	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
	4
[190]	5	;***************************************************
[418]	6	; load .ncl.config or .ncl.config.default
[190]	7	;***************************************************
[194]	8
	9	function which_script()
	10	local script
	11	begin
	12	script="spectra"
	13	return(script)
	14	end
[190]	15
[418]	16	if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
	17	loadscript("$PALM_BIN/../../.ncl.config")
[190]	18	else
[418]	19	if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
	20	loadscript( "$PALM_BIN/NCL/.ncl.config.default")
[190]	21	else
[418]	22	palm_bin_path = getenv("PALM_BIN")
[190]	23	print(" ")
[418]	24	print("Neither the personal configuration file '.ncl.config' exists in")
	25	print("~/palm/current_version")
	26	print("nor the default configuration file '.ncl.config.default' exists in")
	27	print(palm_bin_path + "/NCL")
[190]	28	print(" ")
	29	exit
[176]	30	end if
[418]	31	end if
[218]	32
[190]	33	begin
[176]	34
[190]	35	;***************************************************
[532]	36	; Retrieving the NCL version used
	37	;***************************************************
	38
	39	ncl_version_ch = systemfunc("ncl -V")
	40	ncl_version = stringtofloat(ncl_version_ch)
	41
	42	;***************************************************
[526]	43	; Retrieving the double quote character
	44	;***************************************************
	45
	46	dq=str_get_dq()
	47
	48	;***************************************************
[190]	49	; set up default parameter values and strings
	50	;***************************************************
	51
	52	if (file_1 .EQ. "File in") then
	53	print(" ")
[418]	54	print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]	55	print(" ")
	56	exit
[176]	57	else
	58	file_in = file_1
	59	end if
	60
[532]	61	if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. format_out .NE. "png")then
[190]	62	print(" ")
	63	print("'format_out = "+format_out+"' is invalid and set to'x11'")
	64	print(" ")
	65	format_out="x11"
	66	end if
[532]	67
	68	if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
	69	print(" ")
	70	print("Output of png files not available")
	71	print("png output is avaiable with NCL version 5.2.0 and higher ")
	72	print("NCL version used: " + ncl_version_ch)
	73	print(" ")
	74	exit
	75	end if
[190]	76
[218]	77	if (log_x .NE. 0 .AND. log_x .NE. 1)then
[190]	78	print(" ")
[218]	79	print("'log_x'= "+log_x+" is invalid and set to 1")
[190]	80	print(" ")
[218]	81	log_x = 1
[190]	82	end if
	83
[218]	84	if (log_y .NE. 0 .AND. log_y .NE. 1)then
[190]	85	print(" ")
[218]	86	print("'log_y'= "+log_y+" is invalid and set to 1")
[190]	87	print(" ")
[218]	88	log_y = 1
[190]	89	end if
	90
[218]	91	if (norm_y .EQ. 0.) then
[190]	92	print(" ")
[218]	93	print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
[190]	94	print(" ")
[218]	95	norm_y = 1.0
[176]	96	end if
[218]	97	if (norm_x .EQ. 0.) then
[190]	98	print(" ")
[218]	99	print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
[190]	100	print(" ")
[218]	101	norm_x= 1.0
[176]	102	end if
[190]	103
	104	if (sort .NE. "height" .AND. sort .NE. "time") then
	105	print(" ")
	106	print("'sort'= "+sort+" is invalid and set to 'height'")
	107	print(" ")
	108	sort = "height"
[176]	109	end if
[190]	110
	111	if (black .NE. 0 .AND. black .NE. 1)then
	112	print(" ")
	113	print("'black'= "+black+" is invalid and set to 0")
	114	print(" ")
[176]	115	black = 0
	116	end if
[190]	117
	118	if (dash .NE. 0 .AND. dash .NE. 1)then
	119	print(" ")
	120	print("'dash'= "+dash+" is invalid and set to 0")
	121	print(" ")
[176]	122	dash = 0
	123	end if
	124
[190]	125	;***************************************************
	126	; open input file
	127	;***************************************************
[218]	128
	129	file_in_1 = False
	130	if (isStrSubset(file_in, ".nc"))then
	131	start_f = -2
	132	end_f = -2
	133	file_in_1 = True
	134	end if
[176]	135
[218]	136	if (start_f .EQ. -1)then
	137	print(" ")
	138	print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
	139	print(" ")
	140	exit
	141	end if
	142	if (end_f .EQ. -1)then
	143	print(" ")
	144	print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
	145	print(" ")
	146	exit
	147	end if
	148
	149	files=new(end_f-start_f+1,string)
[176]	150
[218]	151	if (file_in_1)then
	152	if (isfilepresent(file_in))then
	153	files(0)=file_in
	154	else
	155	print(" ")
	156	print("1st input file: '"+file_in+"' does not exist")
	157	print(" ")
	158	exit
	159	end if
	160	else
	161	if (start_f .EQ. 0)then
	162	if (isfilepresent(file_in+".nc"))then
	163	files(0)=file_in+".nc"
	164	do i=1,end_f
	165	if (isfilepresent(file_in+"."+i+".nc"))then
	166	files(i)=file_in+"."+i+".nc"
	167	else
	168	print(" ")
	169	print("Input file: '"+file_in+"."+i+".nc' does not exist")
	170	print(" ")
	171	exit
	172	end if
	173	end do
	174	else
	175	print(" ")
	176	print("Input file: '"+file_in+".nc' does not exist")
	177	print(" ")
	178	exit
	179	end if
	180	else
	181	do i=start_f,end_f
	182	if (isfilepresent(file_in+"."+i+".nc"))then
	183	files(i-start_f)=file_in+"."+i+".nc"
	184	else
	185	print(" ")
	186	print("Input file: '"+file_in+"."+i+".nc' does not exist")
	187	print(" ")
	188	exit
	189	end if
	190	end do
	191	end if
	192	end if
	193
	194	f=addfiles(files,"r")
	195	f_att=addfile(files(0),"r")
	196	ListSetType(f,"cat")
	197
	198	vNam = getfilevarnames(f_att)
[176]	199	print(" ")
[190]	200	print("Variables in input file:")
	201	print("- "+ vNam)
[176]	202	print(" ")
	203	dim = dimsizes(vNam)
[218]	204	vDim = getfiledimsizes(f_att)
[176]	205
[218]	206	t_all = f[:]->time
[176]	207	nt = dimsizes(t_all)
	208	delta_t=t_all(nt-1)/nt
	209
[218]	210	k_x=f_att->k_x
[190]	211	dimx=dimsizes(k_x)
[218]	212	k_y=f_att->k_y
[190]	213	dimy=dimsizes(k_y)
	214
	215
	216	dim_level=dimsizes(height_level)
	217
[176]	218	do i=0,dim-1
	219	if (vNam(i) .EQ. "zu_sp")then
[218]	220	zu=f_att->zu_sp
[190]	221	if (height_level(0) .EQ. -1)then
	222	dimz=dimsizes(zu)
	223	else
	224	if (dim_level .GT. dimsizes(zu))then
	225	print(" ")
	226	print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
	227	print(" ")
	228	exit
	229	else
	230	zuh=new(dim_level,double)
	231	do le=0,dim_level-1
	232	zuh(le)=zu(height_level(le))
	233	end do
	234	dimz=dim_level
	235	end if
	236	end if
[176]	237	else
	238	if (vNam(i) .EQ. "zw_sp")then
[218]	239	zw=f_att->zw_sp
[190]	240	if (height_level(0) .EQ. -1)then
	241	dimz=dimsizes(zw)
	242	else
	243	if (dim_level .GT. dimsizes(zw))then
	244	print(" ")
	245	print("'height_level' has more elements than available height levels in input file (= "+dimz+")")
	246	print(" ")
	247	exit
	248	else
	249	zwh=new(dim_level,double)
	250	do le=0,dim_level-1
	251	zwh(le)=zw(height_level(le))
	252	end do
	253	dimz=dim_level
	254	end if
	255	end if
[176]	256	end if
	257	end if
	258	end do
[190]	259
	260	;****************************************************
[176]	261	; start of time step and different types of mistakes that could be done
[190]	262	;****************************************************
[176]	263
[190]	264	if (start_time_step .EQ. -1.d) then
[176]	265	start_time_step=t_all(0)/3600
[190]	266	else
[176]	267	if (start_time_step .GT. t_all(nt-1)/3600)then
	268	print(" ")
	269	print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
	270	print(" ")
[218]	271	print("Select another 'start_time_step'")
[176]	272	print(" ")
	273	exit
	274	end if
	275	if (start_time_step .LT. t_all(0)/3600)then
	276	print(" ")
	277	print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
	278	exit
	279	end if
	280	end if
	281
	282	do i=0,nt-1
[190]	283	if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]	284	st=i
	285	break
	286	end if
	287	end do
	288
[190]	289	if (.not. isvar("st"))then
	290	print(" ")
	291	print("'start_time_step' = "+ start_time_step +"h is invalid")
	292	print(" ")
[218]	293	print("Select another 'start_time_step'")
[190]	294	print(" ")
	295	exit
	296	end if
	297
	298	;****************************************************
[176]	299	; end of time step and different types of mistakes that could be done
[190]	300	;****************************************************
[176]	301
[190]	302	if (end_time_step .EQ. -1.d) then
[176]	303	end_time_step = t_all(nt-1)/3600
	304	else
	305	if (end_time_step .GT. t_all(nt-1)/3600)then
	306	print(" ")
	307	print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
	308	print(" ")
[218]	309	print("Select another 'end_time_step'")
[176]	310	print(" ")
	311	exit
	312	end if
	313	if (end_time_step .LT. start_time_step/3600)then
	314	print(" ")
	315	print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
	316	print(" ")
[218]	317	print("Select another 'start_time_step' or 'end_time_step'")
[176]	318	print(" ")
	319	exit
	320	end if
	321	end if
	322
	323	do i=0,nt-1
[190]	324	if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[176]	325	et=i
	326	break
	327	end if
	328	end do
[190]	329
	330	if (.not. isvar("et"))then
	331	print(" ")
	332	print("'end_time_step' = "+ end_time_step +"h is invalid")
	333	print(" ")
[218]	334	print("Select another 'end_time_step'")
[190]	335	print(" ")
	336	exit
	337	end if
[176]	338
	339	delete(start_time_step)
	340	start_time_step=round(st,3)
	341	delete(end_time_step)
	342	end_time_step=round(et,3)
	343
	344	print(" ")
	345	print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step)
	346	print(" till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step)
	347	print(" ")
	348
	349	dimt = end_time_step-start_time_step+1
	350
[190]	351	;***************************************************
[176]	352	; set up recourses
[190]	353	;***************************************************
[176]	354
	355	res = True
	356	res@gsnDraw = False
	357	res@gsnFrame = False
	358	res@txFont = "helvetica"
	359	res@tiMainFont = "helvetica"
	360	res@tiXAxisFont = "helvetica"
	361	res@tiYAxisFont = "helvetica"
	362	res@tmXBLabelFont = "helvetica"
	363	res@tmYLLabelFont = "helvetica"
	364	res@lgLabelFont = "helvetica"
	365	res@tmLabelAutoStride = True
	366
[218]	367	res@lgLabelFontHeightF = font_size_legend
	368	res@lgTitleFontHeightF = font_size
	369	res@txFontHeightF = font_size
	370	res@tiXAxisFontHeightF = font_size
	371	res@tiYAxisFontHeightF = font_size
	372	res@tmXBLabelFontHeightF = font_size
	373	res@tmYLLabelFontHeightF = font_size
	374
	375	res@tmXBMinorPerMajor = 4
	376	res@tmYLMinorPerMajor = 4
	377
	378	if (log_x .EQ. 1) then
[176]	379	res@trXLog = True
[190]	380	else
	381	res@trXLog = False
	382	end if
[218]	383	if (log_y .EQ. 1)then
[176]	384	res@trYLog = True
[190]	385	else
[178]	386	res@trYLog = False
[176]	387	end if
	388
	389	legend_label=new(dimt,double)
	390	legend_label_zu=new(dimz,double)
	391	legend_label_zw=new(dimz,double)
	392	legend_label_z=new(dimz,double)
	393	do p=start_time_step,end_time_step
	394	if (t_all(p)/3600 .LT. 1) then
	395	legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True)
	396	else
	397	legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True)
	398	end if
	399	end do
	400	if (sort .EQ. "time")
	401	plot = new(dim*dimz,graphic)
	402	np=dimt
[534]	403	res@lgTitleString = "Time (h)"
[176]	404	else
	405	plot = new(dim*dimt,graphic)
	406	np=dimz
[218]	407
[176]	408	do p=0,dimz-1
[190]	409	if (height_level(0) .EQ. -1)then
	410	legend_label_zu(p)=round(zu(p),3)
	411	legend_label_zw(p)=round(zw(p),3)
	412	else
	413	legend_label_zu(p)=round(zuh(p),3)
	414	legend_label_zw(p)=round(zwh(p),3)
	415	end if
[176]	416	end do
	417	end if
[194]	418
	419	if (black .eq. 0 ) then
	420	if (np .EQ. 1)then
[218]	421	color = 237
[194]	422	else
	423	step=round(235/(np-1),3)
[218]	424	color = ispan(2,237,step)
[194]	425	end if
[324]	426	else
	427	color = 2
[176]	428	end if
	429	if ( dash .eq. 0 ) then
[190]	430	res@xyMonoDashPattern = True
[176]	431	end if
	432
[532]	433	if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
	434	format_out@wkPaperSize = "A4"
	435	end if
	436	if ( format_out .EQ. "png" ) then
	437	format_out@wkWidth = 1000
	438	format_out@wkHeight = 1000
	439	end if
	440
[176]	441	wks=gsn_open_wks(format_out,file_out)
	442	gsn_define_colormap(wks,"rainbow+white")
	443
	444	n=0
	445	do varn =dim-1,0,1
[178]	446
[176]	447	check = True
	448
	449	if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "zu_sp") .OR. isStrSubset( vNam(varn), "zw_sp") .OR. isStrSubset( vNam(varn), "k_x") .OR. isStrSubset( vNam(varn), "k_y")) then
	450	check = False
	451	end if
	452
[190]	453	if (var .NE. "all") then
[176]	454	check = isStrSubset( var,","+vNam(varn)+"," )
	455	end if
	456
[190]	457	if(check) then
[218]	458
	459	temp = f[:]->$vNam(varn)$
	460	data = temp(start_time_step:end_time_step,0:dimz-1,:)
	461
	462	temp_att = f_att->$vNam(varn)$
	463	a=getvardims(temp_att)
[190]	464	b=dimsizes(a)
[534]	465	delete(temp_att)
[176]	466
[190]	467	if (height_level(0) .NE. -1)then
	468	do te=0,dimz-1
[218]	469	print(height_level(te))
	470	data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)
[190]	471	end do
	472	end if
	473
[218]	474	data=data/(norm_y*norm_x) ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
[190]	475
[176]	476	do i=0,b-1
	477	if (isStrSubset( a(i),"zu_sp" ))then
	478	legend_label_z=legend_label_zu
[190]	479	if (height_level(0) .NE. -1)then
	480	z=zuh
	481	else
	482	z=zu
	483	end if
[176]	484	else
	485	if (isStrSubset( a(i),"zw_sp" ))then
[190]	486	legend_label_z=legend_label_zw
	487	if (height_level(0) .NE. -1)then
	488	z=zwh
	489	else
	490	z=zw
	491	end if
[176]	492	end if
	493	end if
	494	end do
[218]	495
	496	min_x=new(dimz,double)
	497	max_x=new(dimz,double)
	498	min_y=new(dimz,float)
	499	max_y=new(dimz,float)
	500	plot_h = new(dimz,graphic)
	501
[176]	502	if (isStrSubset(vNam(varn),"x"))then
[218]	503	x_axis = new((/dimz,dimx/),double)
	504	do q=0,dimz-1
	505	x_axis(q,:) = f_att->k_x
	506	x_axis = x_axis/norm_x
	507	end do
	508	if (norm_x .NE. 1.)then
[534]	509	res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]	510	else
[218]	511	if (norm_height .EQ. 1)then
[534]	512	res@tiXAxisString = "k~B~x~N~ x z (1"
[218]	513	else
[534]	514	res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]	515	end if
[190]	516	end if
[218]	517	dim_r=dimx
[176]	518	else
[218]	519	x_axis=new((/dimz,dimy/),double)
	520	do q=0,dimz-1
	521	x_axis(q,:) = f_att->k_y
	522	x_axis = x_axis/norm_x
	523	end do
	524	if (norm_x .NE. 1.)then
[534]	525	res@tiXAxisString = "k~B~x~N~ ("+unit_x+")"
[190]	526	else
[218]	527	if (norm_height .EQ. 1)then
[534]	528	res@tiXAxisString = "k~B~x~N~ x z (1)"
[218]	529	else
[534]	530	res@tiXAxisString = "k~B~x~N~ (1/m)"
[218]	531	end if
[190]	532	end if
[218]	533	dim_r=dimy
[176]	534	end if
	535
	536	if (sort .EQ. "time")
[218]	537	res@xyLineColors = color
	538	res@pmLegendDisplayMode = "Always"
	539	res@pmLegendSide = "Top"
	540	res@pmLegendParallelPosF = 1.2
	541	res@pmLegendOrthogonalPosF = -1.0
	542	res@pmLegendWidthF = 0.12
	543	res@pmLegendHeightF = 0.04*(end_time_step-start_time_step+1)
[190]	544	do p=dimz-1,0,1
[218]	545	if (log_y .EQ. 1)then
[190]	546	do q=0,dimt-1
[218]	547	do r=0,dim_r-1
	548	if (data(q,p,r) .EQ. 0)then
[190]	549	st=p+start_time_step
	550	print(" ")
	551	print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used")
	552	print(" ")
	553	res@trYLog = False
	554	end if
	555	end do
[176]	556	end do
[190]	557	end if
[176]	558	res@gsnLeftString = vNam(varn)
[218]	559	res@gsnRightString = "Height = "+z(p)+"m"
[534]	560	res@tiYAxisString = "("+unit_y+")"
[218]	561	res@xyExplicitLegendLabels = legend_label
	562	if (norm_height .EQ. 1)then
	563	data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
	564	x_axis(p,:) = x_axis(p,:)*z(p)
	565	end if
	566	res@trXMinF = min(x_axis(p,:))
	567	res@trXMaxF = max(x_axis(p,:))
	568	plot(n) = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
[176]	569	n=n+1
	570	end do
	571	else
[190]	572	if (sort .EQ. "height")
	573	do p=0,dimt-1
[176]	574	do q=0,dimz-1
[218]	575	do r=0,dim_r-1
	576	if (data(p,q,r) .EQ. 0)then
[176]	577	st=p+start_time_step
	578	print(" ")
[190]	579	print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used")
[176]	580	print(" ")
	581	res@trYLog = False
	582	end if
	583	end do
[218]	584	if (norm_height .EQ. 1 .AND. p .EQ. 0)then
	585	data(p,q,:) = data(p,q,:)*doubletofloat(legend_label_z(q))
	586	x_axis(q,:) = x_axis(q,:)*doubletofloat(legend_label_z(q))
	587	end if
	588	max_y(q)=max(data(p,q,:))
	589	min_y(q)=min(data(p,q,:))
	590	min_x(q)=min(x_axis(q,:))
	591	max_x(q)=max(x_axis(q,:))
	592	end do
	593	do q=0,dimz-1
	594	res@xyLineColor = color(q)
	595	if (dash .EQ. 1)then
	596	res@xyDashPattern = q
	597	end if
	598	if (q .EQ. 0)then
[534]	599	res@tiYAxisString = "("+unit_y+")"
[218]	600	res@gsnLeftString = vNam(varn)
	601	res@gsnRightString = "Time = "+legend_label(p)+"h"
	602	res@trXMinF = min(min_x)
	603	res@trXMaxF = max(max_x)
	604	res@trYMinF = min(min_y)
	605	res@trYMaxF = max(max_y)
	606
	607	plot_h(q) = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
	608
	609	lgres = True
	610	if (dash .EQ. 0)then
	611	lgres@lgMonoDashIndex = True
	612	else
	613	lgres@lgDashIndexes = ispan(0,dimz-1,1)
	614	end if
	615	if (black .EQ. 1)then
	616	lgres@lgMonoLineColors = True
	617	else
	618	lgres@lgLineColors = color
	619	end if
[534]	620	lgres@lgTitleString = "Height (m)"
[218]	621	lgres@lgLabelFont = "helvetica"
[534]	622	lgres@lgLabelFontHeightF = font_size_legend*6
[218]	623	lgres@lgTitleFontHeightF = font_size
	624	lgres@vpWidthF = 0.12
[534]	625	lgres@vpHeightF = font_size_legend*(dimz+3)
[218]	626
	627	lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)
	628
	629	amres = True
	630	amres@amParallelPosF = 0.75
	631	amres@amOrthogonalPosF = 0.15
	632	annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
	633	else
	634	plot_h(q) = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
	635	overlay(plot_h(0),plot_h(q))
	636	end if
	637	end do
	638	plot(n)=plot_h(0)
	639	n=n+1
[176]	640	end do
	641	end if
[190]	642	end if
[218]	643	delete(data)
	644	delete(temp)
[178]	645	delete(x_axis)
[218]	646	delete(min_x)
	647	delete(max_x)
	648	delete(min_y)
	649	delete(max_y)
	650	delete(plot_h)
[176]	651	end if
	652	end do
	653
	654	if (n .EQ. 0) then
	655	print(" ")
[190]	656	print("The variables 'var="+var+"' do not exist on your input file;")
	657	print("be sure to have one comma berfore and after each variable")
[176]	658	print(" ")
	659	exit
	660	end if
	661
	662	; ***************************************************
	663	; merge plots onto one page
	664	; ***************************************************
	665
	666	resP = True
	667	resP@txFont = "helvetica"
[218]	668	resP@txString = f_att@title
[176]	669	resP@txFuncCode = "~"
[218]	670	resP@txFontHeightF = 0.015
[176]	671
[534]	672	no_frames = 0
	673
[418]	674	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. n .gt. no_rows*no_columns) then
[183]	675	gsn_panel(wks,plot,(/n,1/),resP)
[250]	676	print(" ")
	677	print("Outputs to .eps or .epsi have only one frame")
	678	print(" ")
[176]	679	else
[218]	680	do i = 0,n-1, no_rows*no_columns
	681	if( (i+no_rows*no_columns) .gt. (n-1)) then
	682	gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
[534]	683	no_frames = no_frames + 1
[176]	684	else
[534]	685	gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
	686	no_frames = no_frames + 1
[176]	687	end if
	688	end do
	689	end if
	690
[532]	691	if (format_out .EQ. "png" ) then
	692	png_output = new((/no_frames/), string)
	693	j = 0
	694	do i=0, no_frames-1
	695	j = i + 1
	696	if (j .LE. 9) then
	697	png_output(i) = file_out+".00000"+j+".png"
	698	end if
	699	if (j .GT. 9 .AND. j .LE. 99) then
	700	png_output(i) = file_out+".0000"+j+".png"
	701	end if
	702	if (j .GT. 99 .AND. j .LE. 999) then
	703	png_output(i) = file_out+".000"+j+".png"
	704	end if
	705	if (j .GT. 999) then
	706	png_output(i) = file_out+".00"+j+".png"
	707	end if
	708
	709	;using imagemagick's convert for reducing the white
	710	;space around the plot
	711	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
	712	png_output(i) + " " + png_output(i)
	713	system(cmd)
	714	end do
	715
	716	print(" ")
	717	print("Output to: "+ png_output)
	718	print(" ")
	719	else
	720	print(" ")
	721	print("Output to: " + file_out +"."+ format_out)
	722	print(" ")
	723	end if
[176]	724
[190]	725	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format