Home

Context Navigation

spectra.ncl @ 4013

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

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

Download in other formats:

Original Format