Home

Context Navigation

spectra.ncl @ 1433

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

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

Download in other formats:

Original Format