Home

Context Navigation

spectra.ncl @ 2699

Last change on this file since 2699 was 2563, checked in by Giersch, 7 years ago
Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data
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
[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
	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)
[2563]	687	res@trYMaxF = max(max_y*10)
[218]	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