Home

Context Navigation

spectra.ncl @ 1873

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

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

Download in other formats:

Original Format