Home

Context Navigation

cross_sections.ncl @ 858

Last change on this file since 858 was 827, checked in by heinze, 13 years ago
allow plotting of data with very small time increments
File size: 77.4 KB

Rev	Line
[154]	1	load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
[157]	2	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
[154]	3	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
	4	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
	5	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
	6
[190]	7	;***************************************************
[418]	8	; load .ncl.config or .ncl.config.default
[190]	9	;***************************************************
[194]	10
	11	function which_script()
	12	local script
	13	begin
	14	script="cross_section"
	15	return(script)
	16	end
[418]	17
	18	if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
	19	loadscript("$PALM_BIN/../../.ncl.config")
[190]	20	else
[418]	21	if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
	22	loadscript( "$PALM_BIN/NCL/.ncl.config.default")
[190]	23	else
[418]	24	palm_bin_path = getenv("PALM_BIN")
[190]	25	print(" ")
[418]	26	print("Neither the personal configuration file '.ncl.config' exists in")
	27	print("~/palm/current_version")
[566]	28	print("nor the default configuration file '.ncl.config.default' "+\
	29	"exists in")
[418]	30	print(palm_bin_path + "/NCL")
[190]	31	print(" ")
	32	exit
	33	end if
[418]	34	end if
[194]	35
[154]	36	begin
	37
[532]	38	;***************************************************
	39	; Retrieving the NCL version used
	40	;***************************************************
	41
	42	ncl_version_ch = systemfunc("ncl -V")
	43	ncl_version = stringtofloat(ncl_version_ch)
	44
[154]	45	; ***************************************************
[526]	46	; Retrieving the double quote character
	47	; ***************************************************
	48
	49	dq=str_get_dq()
	50
	51	; ***************************************************
[566]	52	; set default parameter values and strings if not assigned
	53	; in prompt or parameter list
[154]	54	; ***************************************************
[566]	55
	56	; Selection of type of cross section is done in palmplot
	57	if (.not. isvar("xyc"))then
	58	xyc = 0
	59	end if
	60	if (.not. isvar("xzc"))then
	61	xzc = 0
	62	end if
	63	if (.not. isvar("yzc"))then
	64	yzc = 0
	65	end if
[154]	66
[190]	67	if (file_1 .EQ. "File in") then
	68	print(" ")
[418]	69	print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]	70	print(" ")
	71	exit
[175]	72	else
	73	file_in = file_1
[154]	74	end if
[175]	75
[566]	76	if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. \
	77	format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. \
	78	format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
	79	format_out .NE. "png") then
[190]	80	print(" ")
	81	print("'format_out = "+format_out+"' is invalid and set to'x11'")
	82	print(" ")
	83	format_out="x11"
[154]	84	end if
[175]	85
[532]	86	if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
	87	print(" ")
	88	print("Output of png files not available")
	89	print("png output is avaiable with NCL version 5.2.0 and higher ")
	90	print("NCL version used: " + ncl_version_ch)
	91	print(" ")
	92	exit
	93	end if
	94
[190]	95	if (sort .NE. "layer" .AND. sort .NE. "time") then
	96	print(" ")
	97	print("'sort'= "+sort+" is invalid and set to 'layer'")
	98	print(" ")
	99	sort = "layer"
[154]	100	end if
[190]	101
	102	if (mode .NE. "Fill" .AND. mode .NE. "Line" .AND. mode .NE. "Both") then
	103	print(" ")
	104	print("'mode'= "+mode+" is invalid and set to 'Fill'")
	105	print(" ")
[154]	106	mode = "Fill"
	107	end if
[175]	108
[566]	109	if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND.\
	110	fill_mode .NE. "CellFill") then
[190]	111	print(" ")
[218]	112	print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
[190]	113	print(" ")
[154]	114	fill_mode = "AreaFill"
	115	end if
[190]	116
	117	if (shape .NE. 0 .AND. shape .NE. 1) then
	118	print(" ")
	119	print("'shape'= "+shape+" is invalid and set to 1")
	120	print(" ")
[154]	121	shape = 1
[175]	122	end if
[513]	123
[190]	124	if (xyc .NE. 0 .AND. xyc .NE. 1)then
	125	print(" ")
	126	print("'xyc'= "+xyc+" is invalid and set to 0")
	127	print(" ")
	128	xyc = 0
[154]	129	end if
[190]	130
	131	if (xzc .NE. 0 .AND. xzc .NE. 1)then
	132	print(" ")
	133	print("'xzc'= "+xzc+" is invalid and set to 0")
	134	print(" ")
	135	xyc = 0
[154]	136	end if
[190]	137
	138	if (yzc .NE. 0 .AND. yzc .NE. 1)then
	139	print(" ")
	140	print("'yzc'= "+yzc+" is invalid and set to 0")
	141	print(" ")
	142	xyc = 0
[154]	143	end if
[190]	144
[154]	145	if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
	146	print(" ")
[218]	147	print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]	148	print(" ")
	149	exit
	150	end if
	151	if (xyc .EQ. 1 ) then
	152	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
	153	print(" ")
[218]	154	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]	155	print(" ")
	156	exit
	157	end if
	158	end if
	159	if (xzc .EQ. 1) then
	160	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
	161	print(" ")
[218]	162	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]	163	print(" ")
	164	exit
	165	end if
	166	end if
	167	if (yzc .EQ. 1) then
	168	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
	169	print(" ")
[218]	170	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]	171	print(" ")
	172	exit
	173	end if
	174	end if
	175
[190]	176	if (vector .NE. 0 .AND. vector .NE. 1) then
	177	print(" ")
[218]	178	print("Set 'vector' to 0 or 1")
[190]	179	print(" ")
	180	exit
[218]	181	end if
	182
	183	if (norm_x .EQ. 0) then
	184	print(" ")
	185	print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
	186	print(" ")
	187	norm_x = 1.0
	188	end if
	189	if (norm_y .EQ. 0) then
	190	print(" ")
	191	print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
	192	print(" ")
	193	norm_y = 1.0
	194	end if
	195	if (norm_z .EQ. 0) then
	196	print(" ")
	197	print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
	198	print(" ")
	199	norm_z = 1.0
	200	end if
[190]	201
[154]	202	; ***************************************************
	203	; open input file
	204	; ***************************************************
[218]	205
	206	file_in_1 = False
	207	if (isStrSubset(file_in, ".nc"))then
	208	start_f = -2
	209	end_f = -2
	210	file_in_1 = True
	211	end if
[154]	212
[218]	213	if (start_f .EQ. -1)then
	214	print(" ")
[566]	215	print("'start_f' must be one of the cyclic numbers (at least 0) "+\
	216	"of your input file(s)")
[218]	217	print(" ")
	218	exit
	219	end if
	220	if (end_f .EQ. -1)then
	221	print(" ")
[566]	222	print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
	223	"your input file(s)")
[218]	224	print(" ")
	225	exit
	226	end if
[154]	227
[218]	228	files=new(end_f-start_f+1,string)
	229
	230	if (file_in_1)then
	231	if (isfilepresent(file_in))then
	232	files(0)=file_in
	233	else
	234	print(" ")
	235	print("1st input file: '"+file_in+"' does not exist")
	236	print(" ")
	237	exit
	238	end if
	239	else
	240	if (start_f .EQ. 0)then
	241	if (isfilepresent(file_in+".nc"))then
	242	files(0)=file_in+".nc"
	243	do i=1,end_f
	244	if (isfilepresent(file_in+"."+i+".nc"))then
	245	files(i)=file_in+"."+i+".nc"
	246	else
	247	print(" ")
	248	print("Input file: '"+file_in+"."+i+".nc' does not exist")
	249	print(" ")
	250	exit
	251	end if
	252	end do
	253	else
	254	print(" ")
	255	print("Input file: '"+file_in+".nc' does not exist")
	256	print(" ")
	257	exit
	258	end if
	259	else
	260	do i=start_f,end_f
	261	if (isfilepresent(file_in+"."+i+".nc"))then
	262	files(i-start_f)=file_in+"."+i+".nc"
	263	else
	264	print(" ")
	265	print("Input file: '"+file_in+"."+i+".nc' does not exist")
	266	print(" ")
	267	exit
	268	end if
	269	end do
	270	end if
	271	end if
	272
	273	f=addfiles(files,"r")
	274	f_att=addfile(files(0),"r")
	275	ListSetType(f,"cat")
	276
	277	vNam = getfilevarnames(f_att)
[585]	278	vType = getfilevartypes(f_att,vNam)
	279
	280	if ((all(vType .eq. "double"))) then ;distinction if data is double or float
	281	check_vType = True
	282	else
	283	check_vType = False
	284	end if
	285
[154]	286	print(" ")
[190]	287	print("Variables in input file:")
	288	print("- "+ vNam)
[154]	289	print(" ")
	290	dim = dimsizes(vNam)
[162]	291
	292	; ***************************************************
	293	; check for kind of input file
	294	; ***************************************************
	295
	296	check_3d = True
	297	do varn=0,dim-1
	298	checkxy = isStrSubset( vNam(varn),"xy")
	299	checkxz = isStrSubset( vNam(varn),"xz")
	300	checkyz = isStrSubset( vNam(varn),"yz")
	301	if (checkxy .OR. checkxz .OR. checkyz) then
	302	check_3d = False
	303	break
	304	end if
	305	end do
	306	if (.not. check_3d) then
[195]	307	if (xyc .EQ. 1)then
	308	if (.not. checkxy) then
	309	print(" ")
	310	print("Your input file doesn't have values for xy cross-sections")
	311	if (checkxz)then
	312	print("Select another input file or xz cross-sections")
	313	else
	314	print("Select another input file or yz cross-sections")
	315	end if
	316	print(" ")
	317	exit
[162]	318	else
[195]	319	print(" ")
	320	print("Your input file contains xy data")
	321	print(" ")
	322	end if
[162]	323	end if
[195]	324	if (xzc .EQ. 1)then
	325	if (.not. checkxz) then
	326	print(" ")
	327	print("Your input file doesn't have values for xz cross-sections")
	328	if (checkxy)then
	329	print("Select another input file or xy cross-sections")
	330	else
	331	print("Select another input file or yz cross-sections")
	332	end if
	333	print(" ")
	334	exit
[162]	335	else
[195]	336	print(" ")
	337	print("Your input file contains xz data")
	338	print(" ")
	339	end if
[162]	340	end if
[195]	341	if (yzc .EQ. 1)then
	342	if (.not. checkyz) then
	343	print(" ")
	344	print("Your input file doesn't have values for yz cross-sections")
	345	if (checkxy)then
	346	print("Select another input file or xy cross-sections")
	347	else
	348	print("Select another input file or xz cross-sections")
	349	end if
	350	print(" ")
	351	exit
[162]	352	else
[195]	353	print(" ")
	354	print("Your input file contains yz data")
	355	print(" ")
	356	end if
[162]	357	end if
	358	else
	359	print(" ")
[218]	360	print("Your input file contains 3d or other data")
[162]	361	print(" ")
	362	end if
	363
[161]	364	if (dim .EQ. 0) then
	365	print(" ")
[162]	366	print("There is no data on file")
[161]	367	print(" ")
[162]	368	exit
[161]	369	end if
	370
[154]	371	; ***************************************************
	372	; set up recourses
	373	; ***************************************************
	374
	375	cs_res = True
[218]	376	vecres = True
[154]	377	cs_res@gsnYAxisIrregular2Linear = True
[566]	378	vecres@gsnYAxisIrregular2Linear = True
[154]	379
	380	cs_res@gsnDraw = False
	381	cs_res@gsnFrame = False
[566]	382	cs_res@gsnMaximize = True
[218]	383
	384	cs_res@tmXBLabelFontHeightF = font_size
	385	cs_res@tmYLLabelFontHeightF = font_size
	386	cs_res@tiXAxisFontHeightF = font_size
	387	cs_res@tiYAxisFontHeightF = font_size
[154]	388	cs_res@tmXBMode ="Automatic"
	389	cs_res@tmYLMode ="Automatic"
[218]	390
[162]	391	cs_res@cnLevelSelectionMode = "ManualLevels"
[218]	392	cs_res@lbLabelFontHeightF = font_size_legend
	393	cs_res@lbLabelStride = legend_label_stride
	394
[162]	395
[154]	396	cs_resP = True
[218]	397	cs_resP@txString = f_att@title
[175]	398	cs_resP@txFuncCode = "~"
[218]	399	cs_resP@txFontHeightF = 0.015
[154]	400
	401	if ( mode .eq. "Fill" ) then
	402	cs_res@cnFillOn = True
	403	cs_res@gsnSpreadColors = True
[218]	404	cs_res@cnFillMode = fill_mode
[154]	405	cs_res@cnLinesOn = False
	406	cs_res@cnLineLabelsOn = False
	407	end if
	408
	409	if ( mode .eq. "Both" ) then
	410	cs_res@cnFillOn = True
	411	cs_res@gsnSpreadColors = True
[218]	412	cs_res@cnFillMode = fill_mode
[154]	413	cs_res@cnLinesOn = True
	414	cs_res@cnLineLabelsOn = True
	415	end if
	416
[157]	417	; *********************************************
	418	; set up of start_time_step and end_time_step
	419	; *********************************************
[154]	420
[827]	421	t_all = f[:]->time
	422	nt = dimsizes(t_all)
[162]	423
[827]	424	if (nt .EQ. 1)then
	425	delta_t=t_all(nt-1)/nt
	426	else
	427	delta_t_array = new(nt-1,double)
	428
	429	do i=0,nt-2
	430	delta_t_array(i) = t_all(i+1)-t_all(i)
	431	end do
	432
	433	delta_t = min(delta_t_array)
	434	delete(delta_t_array)
	435	end if
	436
	437
[162]	438	; ****************************************************
[157]	439	; start of time step and different types of mistakes that could be done
[162]	440	; ****************************************************
[154]	441
[190]	442	if (start_time_step .EQ. -1.) then
	443	start_time_step=t_all(0)/3600
[157]	444	else
[162]	445	if (start_time_step .GT. t_all(nt-1)/3600)
[157]	446	print(" ")
[566]	447	print("'start_time_step' = "+ start_time_step +"h is greater than "+\
	448	"last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[157]	449	print(" ")
[218]	450	print("Select another 'start_time_step'")
[162]	451	print(" ")
[157]	452	exit
	453	end if
[162]	454	if (start_time_step .LT. t_all(0)/3600)
[157]	455	print(" ")
[566]	456	print("'start_time_step' = "+ start_time_step +"h is lower than "+\
	457	"first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]	458	print(" ")
[218]	459	print("Select another 'start_time_step'")
[162]	460	print(" ")
[157]	461	exit
	462	end if
	463	end if
[218]	464	do i=0,nt-1
[566]	465	if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
	466	start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[162]	467	st=i
	468	break
[769]	469	else
	470	st=0
[162]	471	end if
	472	end do
[769]	473
[162]	474
[157]	475	; ****************************************************
	476	; end of time step and different types of mistakes that could be done
	477	; ****************************************************
	478
[190]	479	if (end_time_step .EQ. -1.) then
	480	end_time_step = t_all(nt-1)/3600
	481	else
[162]	482	if (end_time_step .LT. t_all(0)/3600)
[157]	483	print(" ")
[566]	484	print("'end_time_step' = "+end_time_step+ "h is lower than "+\
	485	"first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]	486	print(" ")
[218]	487	print("Select another 'end_time_step'")
[162]	488	print(" ")
[157]	489	exit
	490	end if
[162]	491	if (end_time_step .GT. t_all(nt-1)/3600)
[157]	492	print(" ")
[566]	493	print("'end_time_step' = "+ end_time_step +"h is greater than "+\
	494	"last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[162]	495	print(" ")
[218]	496	print("Select another 'end_time_step'")
[162]	497	print(" ")
[157]	498	exit
	499	end if
[162]	500	if (end_time_step .LT. start_time_step/3600)
[157]	501	print(" ")
[566]	502	print("'end_time_step' = "+ end_time_step +"h is lower than "+\
	503	"'start_time_step' = "+start_time_step+"h")
[157]	504	print(" ")
[218]	505	print("Select another 'start_time_step' or 'end_time_step'")
[162]	506	print(" ")
[157]	507	exit
	508	end if
	509	end if
[162]	510	do i=0,nt-1
[566]	511	if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
	512	end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[162]	513	et=i
	514	break
[769]	515	else
	516	et=0
[162]	517	end if
	518	end do
	519
	520	delete(start_time_step)
	521	start_time_step=round(st,3)
	522	delete(end_time_step)
	523	end_time_step=round(et,3)
[157]	524
[162]	525	print(" ")
[566]	526	print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
	527	t_all(start_time_step)+" s => index = "+start_time_step)
	528	print(" till "+t_all(end_time_step)/3600+" h = " +\
	529	t_all(end_time_step)+" s => index = "+end_time_step)
[162]	530	print(" ")
	531
[161]	532	no_time=(end_time_step-start_time_step)+1
	533
[157]	534	; ****************************************************
	535	; set up variables for vector plot if required
	536	; ****************************************************
	537
	538	if (vector .EQ. 1) then
[190]	539	if (vec1 .EQ. "component1") then
	540	print(" ")
[566]	541	print("Indicate Vector 1 ('vec1') for Vector-Plot or "+\
	542	"set 'vector' to 0")
[190]	543	print(" ")
	544	exit
[157]	545	end if
[190]	546	if (vec2 .EQ. "component2") then
	547	print(" ")
[566]	548	print("Indicate Vector 2 ('vec2') for Vector-Plot or "+\
	549	"set 'vector' to 0")
[190]	550	print(" ")
	551	exit
	552	end if
[157]	553	end if
	554
	555	check_vec1 = False
	556	check_vec2 = False
	557	check_vecp = False
	558
	559	; ****************************************************
	560	; get data for plots
	561	; ****************************************************
	562
[162]	563	if (xyc .EQ. 1) then
	564	do varn=0,dim-1
[190]	565	if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then
[218]	566	x_d = f_att->$vNam(varn)$
[162]	567	xdim = dimsizes(x_d)-1
	568	delta_x = x_d(1)-x_d(0)
	569	break
	570	end if
	571	end do
	572	do varn=0,dim-1
	573	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then
[218]	574	y_d = f_att->$vNam(varn)$
[162]	575	ydim = dimsizes(y_d)-1
	576	delta_y = y_d(1)-y_d(0)
	577	break
	578	end if
	579	end do
	580	do varn=0,dim-1
	581	if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
[218]	582	z_d = f_att->$vNam(varn)$
[162]	583	zdim = dimsizes(z_d)-1
	584	delta_z = 0
	585	break
	586	else
	587	if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
[218]	588	z_d = f_att->$vNam(varn)$
[162]	589	zdim = dimsizes(z_d)-1
	590	delta_z = -1.d
	591	break
	592	else
	593	zdim = 0
	594	delta_z = -1.d
	595	end if
	596	end if
[357]	597	if (vNam(varn) .eq. "zu1_xy") then
	598	zu1 = f_att->$vNam(varn)$
	599	end if
[162]	600	end do
	601	end if
	602	if (xzc .EQ. 1) then
	603	do varn=0,dim-1
	604	if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x") then
[218]	605	x_d = f_att->$vNam(varn)$
[162]	606	xdim = dimsizes(x_d)-1
	607	delta_x = x_d(1)-x_d(0)
	608	break
	609	end if
	610	end do
	611	do varn=0,dim-1
[566]	612	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
	613	vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
[218]	614	z_d = f_att->$vNam(varn)$
[162]	615	zdim = dimsizes(z_d)-1
	616	delta_z = z_d(1)-z_d(0)
	617	break
	618	end if
	619	end do
	620	do varn=0,dim-1
	621	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
[218]	622	y_d = f_att->$vNam(varn)$
[162]	623	ydim = dimsizes(y_d)-1
	624	delta_y = y_d(1)-y_d(0)
	625	break
	626	else
	627	if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
[218]	628	y_d = f_att->$vNam(varn)$
[162]	629	ydim = dimsizes(y_d)-1
	630	delta_y = -1.d
	631	break
	632	else
	633	ydim = 0
	634	delta_y = -1.d
	635	end if
	636	end if
	637	end do
	638	end if
	639	if (yzc .EQ. 1) then
	640	do varn=0,dim-1
	641	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
[218]	642	y_d = f_att->$vNam(varn)$
[162]	643	ydim = dimsizes(y_d)-1
	644	delta_y = y_d(1)-y_d(0)
	645	break
	646	end if
	647	end do
	648	do varn=0,dim-1
[566]	649	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
	650	vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
[218]	651	z_d = f_att->$vNam(varn)$
[162]	652	zdim = dimsizes(z_d)-1
	653	delta_z = z_d(1)-z_d(0)
	654	break
	655	end if
	656	end do
	657	do varn=0,dim-1
	658	if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
[218]	659	x_d = f_att->$vNam(varn)$
[162]	660	xdim = dimsizes(x_d)-1
	661	delta_x = x_d(1)-x_d(0)
	662	break
	663	else
	664	if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
[218]	665	x_d = f_att->$vNam(varn)$
[162]	666	xdim = dimsizes(x_d)-1
	667	delta_x = -1.d
	668	break
	669	else
	670	xdim = 0
	671	delta_x = -1.d
	672	end if
	673	end if
	674	end do
	675	end if
[195]	676
[162]	677	; ****************************************************
	678	; set up ranges of x-, y- and z-coordinates
	679	; ****************************************************
	680
[190]	681	if (xs .EQ. -1.d) then
	682	xs = x_d(0)
[195]	683	if (delta_x .EQ. -1) then
	684	print(" ")
[566]	685	print("You cannot choose a start value for x, "+\
	686	"there are preset layers for x")
[195]	687	print(" ")
	688	xstart=0
	689	end if
[190]	690	else
[162]	691	if (delta_x .EQ. -1) then
	692	print(" ")
[566]	693	print("You cannot choose a start value for x, "+\
	694	"there are preset layers for x")
[162]	695	print(" ")
	696	xstart=0
[161]	697	else
[162]	698	if (xs .LT. 0-delta_x/2) then
	699	print(" ")
[566]	700	print("range start for x-coordinate = "+\
	701	xs+"m is lower than first value = "+(0-delta_x/2)+"m")
[162]	702	print(" ")
	703	exit
	704	end if
	705	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
	706	if (xs .GE. xdim*delta_x) then
	707	print(" ")
[566]	708	print("range start for x-coordinate = "+\
	709	xs+"m is equal or greater than last value = "+\
	710	xdim*delta_x+"m")
[162]	711	print(" ")
	712	exit
	713	end if
	714	else
	715	if (xs .GT. xdim*delta_x) then
	716	print(" ")
[566]	717	print("range start for x-coordinate = "+\
	718	xs+"m is greater than last value = "+\
	719	xdim*delta_x+"m")
[162]	720	print(" ")
	721	exit
	722	end if
	723	end if
[154]	724	end if
[162]	725	end if
[154]	726
[190]	727	if (ys .EQ. -1.d) then
	728	ys = y_d(0)
[195]	729	if (delta_y .EQ. -1) then
	730	print(" ")
[566]	731	print("You cannot choose a start value for y, "+\
	732	"there are preset layers for y")
[195]	733	print(" ")
	734	ystart=0
	735	end if
[162]	736	else
	737	if (delta_y .EQ. -1) then
	738	print(" ")
[566]	739	print("You cannot choose a start value for y, "+\
	740	"there are preset layers for y")
[162]	741	print(" ")
	742	ystart=0
	743	else
	744	if (ys .LT. 0-delta_y/2) then
	745	print(" ")
[566]	746	print("range start for y-coordinate = "+\
	747	ys+"m is lower than first value = "+0-delta_y/2+"m")
[162]	748	print(" ")
	749	exit
	750	end if
	751	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
	752	if (ys .GE. ydim*delta_y) then
	753	print(" ")
[566]	754	print("range start for y-coordinate = "+\
	755	ys+"m is equal or greater than last value = "+\
	756	ydim*delta_y+"m")
[162]	757	print(" ")
	758	exit
	759	end if
[154]	760	else
[162]	761	if (ys .GT. ydim*delta_y) then
[157]	762	print(" ")
[566]	763	print("range start for y-coordinate = "+\
	764	ys+"m is greater than last value = "+ydim*delta_y+"m")
[157]	765	print(" ")
	766	exit
	767	end if
[162]	768	end if
	769	end if
	770	end if
	771
[190]	772	if (zs .EQ. -1) then
[162]	773	zs = 0
[195]	774	if (delta_z .EQ. -1) then
	775	print(" ")
[566]	776	print("You cannot choose a start value for z, "+\
	777	"there are preset layers for z")
[195]	778	print(" ")
	779	end if
[162]	780	else
	781	if (delta_z .EQ. -1) then
	782	print(" ")
[566]	783	print("You cannot choose a start value for z, "+\
	784	"there are preset layers for z")
[162]	785	print(" ")
[358]	786	zs = 0
[162]	787	else
	788	if (zs .LT. 0) then
	789	print(" ")
[566]	790	print("range start for z-coordinate = "+\
	791	zs+" is lower than first gridpoint = 0")
[162]	792	print(" ")
	793	exit
	794	end if
	795	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
	796	if (zs .GE. zdim) then
	797	print(" ")
[566]	798	print("range start for z-coordinate = "+\
	799	zs+" is equal or greater than last gridpoint = "+zdim)
[162]	800	print(" ")
	801	exit
	802	end if
	803	else
	804	if (zs .GT. zdim) then
	805	print(" ")
[566]	806	print("range start for z-coordinate = "+\
	807	zs+" is greater than last gridpoint = "+zdim)
[162]	808	print(" ")
	809	exit
	810	end if
	811	end if
	812	end if
	813	end if
[157]	814
[190]	815	if (xe .EQ. -1) then
	816	xe = x_d(xdim)
[195]	817	if (delta_x .EQ. -1) then
	818	print(" ")
[566]	819	print("You cannot choose an end value for x, "+\
	820	"there are preset layers for x")
[195]	821	print(" ")
	822	xend=xdim
	823	end if
[162]	824	else
	825	if (delta_x .EQ. -1) then
	826	print(" ")
[566]	827	print("You cannot choose an end value for x, "+\
	828	"there are preset layers for x")
[162]	829	print(" ")
	830	xend=xdim
	831	else
	832	if (xe .GT. xdim*delta_x) then
	833	print(" ")
[566]	834	print("range end for x-coordinate = "+\
	835	xe+"m is greater than last value = "+xdim*delta_x+"m")
[162]	836	print(" ")
	837	exit
	838	end if
	839	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
	840	if (xe .LE. 0-delta_x/2)
	841	print(" ")
[566]	842	print("range end for x-coordinate = "+\
	843	xe+"m is equal or lower than first value = "+\
	844	(0-delta_x/2)+"m")
[162]	845	print(" ")
	846	exit
	847	end if
	848	if (xe .LE. xs) then
	849	print(" ")
[566]	850	print("range end for x-coordinate = "+\
	851	xe+"m is equal or lower than start range = "+xs+"m")
[162]	852	print(" ")
	853	exit
	854	end if
[190]	855	if (xe .EQ. xs+1) then
[162]	856	print(" ")
[566]	857	print("range end for x-coordinate = "+\
	858	xe+"m must be at least two more gridpoints "+\
	859	"greater than start range = "+xs+"m")
[162]	860	print(" ")
	861	exit
	862	end if
[154]	863	else
[162]	864	if (xe .LT. 0-delta_x/2)
[154]	865	print(" ")
[566]	866	print("range end for x-coordinate = "+\
	867	xe+"m is lower than first value = "+(0-delta_x/2)+"m")
[154]	868	print(" ")
[162]	869	exit
[154]	870	end if
[162]	871	if (xe .LT. xs) then
[154]	872	print(" ")
[566]	873	print("range end for x-coordinate = "+\
	874	xe+"m is lower than start range = "+xs+"m")
[154]	875	print(" ")
	876	exit
	877	end if
[162]	878	end if
	879	end if
	880	end if
[190]	881
	882	if (ye .EQ. -1) then
	883	ye = y_d(ydim)
[195]	884	if (delta_y .EQ. -1) then
	885	print(" ")
[566]	886	print("You cannot choose an end value for y, "+\
	887	"there are preset layers for y")
[195]	888	print(" ")
	889	yend=ydim
	890	end if
[162]	891	else
	892	if (delta_y .EQ. -1) then
	893	print(" ")
[566]	894	print("You cannot choose an end value for y, "+\
	895	"there are preset layers for y")
[162]	896	print(" ")
	897	yend=ydim
	898	else
	899	if (ye .GT. ydim*delta_y) then
	900	print(" ")
[566]	901	print("range end for y-coordinate = "+\
	902	ye+"m is greater than last value = "+ydim*delta_y+"m")
[162]	903	print(" ")
	904	exit
	905	end if
	906	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
	907	if (ye .LE. 0-delta_y/2)
	908	print(" ")
[566]	909	print("range end for y-coordinate = "+\
	910	ye+"m is equal or lower than first value = "+\
	911	(0-delta_y/2)+"m")
[162]	912	print(" ")
	913	exit
	914	end if
	915	if (ye .LE. ys) then
	916	print(" ")
[566]	917	print("range end for y-coordinate = "+\
	918	ye+"m is equal or lower than start range = "+ys+"m")
[162]	919	print(" ")
	920	exit
	921	end if
	922	if (ye .EQ. ys+1) then
	923	print(" ")
[566]	924	print("range end for y-coordinate = "+\
	925	ye+"m must be at least two more gridpoints greater "+\
	926	"than start range = "+ys+"m")
[162]	927	print(" ")
	928	exit
	929	end if
[154]	930	else
[162]	931	if (ye .LT. 0-delta_y/2)
[154]	932	print(" ")
[566]	933	print("range end for y-coordinate = "+\
	934	ye+"m is lower than first value = "+(0-delta_y/2)+"m")
[162]	935	print(" ")
	936	exit
	937	end if
	938	if (ye .LT. ys) then
	939	print(" ")
[566]	940	print("range end for y-coordinate = "+\
	941	ye+"m is lower than start range = "+ys+"m")
[162]	942	print(" ")
	943	exit
	944	end if
	945	end if
	946	end if
	947	end if
	948
[190]	949	if (ze .EQ. -1) then
[162]	950	ze = zdim
[195]	951	if (delta_z .EQ. -1) then
	952	print(" ")
[566]	953	print("You cannot choose an end value for z, "+\
	954	"there are preset layers for z")
[195]	955	print(" ")
	956	ze = zdim
	957	end if
[162]	958	else
	959	if (delta_z .EQ. -1) then
	960	print(" ")
[566]	961	print("You cannot choose an end value for z, "+\
	962	"there are preset layers for z")
[162]	963	print(" ")
	964	ze = zdim
	965	else
	966	if (ze .GT. zdim) then
	967	print(" ")
[566]	968	print("range end for z-coordinate = "+\
	969	ze+" is greater than last gridpoint = "+zdim)
[162]	970	print(" ")
	971	exit
	972	end if
	973	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
	974	if (ze .LE. 0)
	975	print(" ")
[566]	976	print("range end for z-coordinate = "+\
	977	ze+" is equal or lower than first gridpoint = 0")
[162]	978	print(" ")
	979	exit
	980	end if
	981	if (ze .LE. zs) then
	982	print(" ")
[566]	983	print("range end for z-coordinate = "+\
	984	ze+" is equal or lower than start range = "+zs)
[162]	985	print(" ")
	986	exit
	987	end if
	988	if (ze .EQ. zs+1) then
	989	print(" ")
[566]	990	print("range end for z-coordinate = "+\
	991	ze+" must be at least two more gridpoints "+\
	992	"greater than start range = "+zs)
[162]	993	print(" ")
	994	exit
	995	end if
	996	else
	997	if (ze .LT. 0)
	998	print(" ")
[566]	999	print("range end for z-coordinate = "+\
	1000	ze+" is lower than first gridpoint = 0")
[162]	1001	print(" ")
	1002	exit
	1003	end if
	1004	if (ze .LT. zs) then
	1005	print(" ")
[566]	1006	print("range end for z-coordinate = "+\
	1007	ze+" is lower than start range = "+zs)
[162]	1008	print(" ")
	1009	exit
	1010	end if
	1011	end if
	1012	end if
	1013	end if
[157]	1014
[162]	1015	if (delta_x .NE. -1) then
	1016	do i=0,xdim
	1017	if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
	1018	xstart=i
	1019	break
	1020	end if
	1021	end do
	1022	do i=0,xdim
	1023	if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
	1024	xend=i
	1025	break
	1026	end if
	1027	end do
	1028	end if
	1029	if (delta_y .NE. -1) then
[190]	1030	do i=0,ydim
[162]	1031	if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
	1032	ystart=i
	1033	break
	1034	end if
	1035	end do
[190]	1036	do i=0,ydim
[162]	1037	if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
	1038	yend=i
	1039	break
	1040	end if
	1041	end do
	1042	end if
[286]	1043
[190]	1044	if( shape .eq. 1 ) then
[195]	1045	if (xyc .EQ. 1)then
	1046	cs_res@vpWidthF = (xe-xs)/(ye-ys)
	1047	cs_res@vpHeightF = 1
[218]	1048	vecres@vpWidthF = (xe-xs)/(ye-ys)
	1049	vecres@vpHeightF = 1
	1050	if (xe-xs .GT. ye-ys)then
[286]	1051	if (no_rows .LE. no_columns)then
	1052	cs_res@gsnPaperOrientation = "landscape"
	1053	vecres@gsnPaperOrientation = "landscape"
	1054	cs_res@lbOrientation = "Horizontal"
	1055	else
	1056	cs_res@gsnPaperOrientation = "portrait"
	1057	vecres@gsnPaperOrientation = "portrait"
	1058	cs_res@lbOrientation = "Vertical"
	1059	end if
[218]	1060	else
[286]	1061	if (no_rows .GE. no_columns)then
	1062	cs_res@gsnPaperOrientation = "portrait"
	1063	vecres@gsnPaperOrientation = "portrait"
	1064	cs_res@lbOrientation = "Vertical"
	1065	else
	1066	cs_res@gsnPaperOrientation = "landscape"
	1067	vecres@gsnPaperOrientation = "landscape"
	1068	cs_res@lbOrientation = "Horizontal"
	1069	end if
[218]	1070	end if
[195]	1071	end if
	1072	if (xzc .EQ. 1)then
	1073	cs_res@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
	1074	cs_res@vpHeightF = 1
[218]	1075	vecres@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
	1076	vecres@vpHeightF = 1
	1077	if (xe-xs .GT. (delta_x*(ze-zs)))then
[286]	1078	if (no_rows .LE. no_columns)then
	1079	cs_res@gsnPaperOrientation = "landscape"
	1080	vecres@gsnPaperOrientation = "landscape"
	1081	cs_res@lbOrientation = "Horizontal"
	1082	else
	1083	cs_res@gsnPaperOrientation = "portrait"
	1084	vecres@gsnPaperOrientation = "portrait"
	1085	cs_res@lbOrientation = "Vertical"
	1086	end if
[218]	1087	else
[286]	1088	if (no_rows .GE. no_columns)then
	1089	cs_res@gsnPaperOrientation = "portrait"
	1090	vecres@gsnPaperOrientation = "portrait"
	1091	cs_res@lbOrientation = "Vertical"
	1092	else
	1093	cs_res@gsnPaperOrientation = "landscape"
	1094	vecres@gsnPaperOrientation = "landscape"
	1095	cs_res@lbOrientation = "Horizontal"
	1096	end if
[218]	1097	end if
[195]	1098	end if
	1099	if (yzc .EQ. 1)then
	1100	cs_res@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
	1101	cs_res@vpHeightF = 1
[218]	1102	vecres@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
	1103	vecres@vpHeightF = 1
	1104	if (ye-ys .GT. (delta_y*(ze-zs)))then
[286]	1105	if (no_rows .LE. no_columns)then
	1106	cs_res@gsnPaperOrientation = "landscape"
	1107	vecres@gsnPaperOrientation = "landscape"
	1108	cs_res@lbOrientation = "Horizontal"
	1109	else
	1110	cs_res@gsnPaperOrientation = "portrait"
	1111	vecres@gsnPaperOrientation = "portrait"
	1112	cs_res@lbOrientation = "Vertical"
	1113	end if
[218]	1114	else
[286]	1115	if (no_rows .GE. no_columns)then
	1116	cs_res@gsnPaperOrientation = "portrait"
	1117	vecres@gsnPaperOrientation = "portrait"
	1118	cs_res@lbOrientation = "Vertical"
	1119	else
	1120	cs_res@gsnPaperOrientation = "landscape"
	1121	vecres@gsnPaperOrientation = "landscape"
	1122	cs_res@lbOrientation = "Horizontal"
	1123	end if
[218]	1124	end if
[195]	1125	end if
[190]	1126	end if
[218]	1127
[162]	1128	delete(xs)
	1129	delete(xe)
	1130	delete(ys)
	1131	delete(ye)
	1132
	1133	xs=xstart
	1134	xe=xend
	1135	ys=ystart
	1136	ye=yend
	1137
[218]	1138	if (xyc .EQ. 1) then
[329]	1139	d= (ye-ys+1)/(major_ticks_y-1)
[218]	1140	e=(xe-xs+1)/(major_ticks_x-1)
	1141	array_yl =new(major_ticks_y,integer)
	1142	array_yl_labels=new(major_ticks_y,double)
	1143	array_minor_yl =new((major_ticks_y-1)*4,double)
	1144	array_xb =new(major_ticks_x,integer)
	1145	array_xb_labels=new(major_ticks_x,double)
	1146	array_minor_xb =new((major_ticks_x-1)*4,double)
	1147	array_yl(0)=ys
	1148	array_xb(0)=xs
	1149	array_yl_labels(0)=y_d(ys)
	1150	array_xb_labels(0)=x_d(xs)
	1151	do ar=1,major_ticks_y-1
	1152	array_yl(ar)=d*(ar-1)+d-1
	1153	array_yl_labels(ar) = y_d(array_yl(ar))
	1154	if (ar .GT. 0)
	1155	do min_ar=0,3
[566]	1156	array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+\
	1157	(y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1)
[218]	1158	end do
	1159	end if
	1160	end do
	1161	do br=1,major_ticks_x-1
	1162	array_xb(br)=e*(br-1)+e-1
	1163	array_xb_labels(br) = x_d(array_xb(br))
	1164	if (br .GT. 0)
	1165	do min_br=0,3
[566]	1166	array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
	1167	(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
[218]	1168	end do
	1169	end if
	1170	end do
	1171	end if
	1172
	1173	if (xzc .EQ. 1) then
	1174	d=(ze-zs+1)/(major_ticks_y-1)
	1175	e=(xe-xs+1)/(major_ticks_x-1)
	1176	array_yl =new(major_ticks_y,integer)
	1177	array_yl_labels=new(major_ticks_y,double)
	1178	array_minor_yl =new((major_ticks_y-1)*4,double)
	1179	array_xb =new(major_ticks_x,integer)
	1180	array_xb_labels=new(major_ticks_x,double)
	1181	array_minor_xb =new((major_ticks_x-1)*4,double)
	1182	array_yl(0)=zs
	1183	array_xb(0)=xs
	1184	array_yl_labels(0)=z_d(zs)
	1185	array_xb_labels(0)=x_d(xs)
	1186	do ar=1,major_ticks_y-1
	1187	array_yl(ar)=d*(ar-1)+d-1
	1188	array_yl_labels(ar) = z_d(array_yl(ar))
	1189	if (ar .GT. 0)
	1190	do min_ar=0,3
[566]	1191	array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
	1192	(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
[218]	1193	end do
	1194	end if
	1195	end do
	1196	do br=1,major_ticks_x-1
	1197	array_xb(br)=e*(br-1)+e-1
	1198	array_xb_labels(br) = x_d(array_xb(br))
	1199	if (br .GT. 0)
	1200	do min_br=0,3
[566]	1201	array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
	1202	(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
[218]	1203	end do
	1204	end if
	1205	end do
	1206	end if
	1207
	1208	if (yzc .EQ. 1) then
	1209	d=(ze-zs+1)/(major_ticks_y-1)
	1210	e=(ye-ys+1)/(major_ticks_x-1)
	1211	array_yl =new(major_ticks_y,integer)
	1212	array_yl_labels=new(major_ticks_y,double)
	1213	array_minor_yl =new((major_ticks_y-1)*4,double)
	1214	array_xb =new(major_ticks_x,integer)
	1215	array_xb_labels=new(major_ticks_x,double)
	1216	array_minor_xb =new((major_ticks_x-1)*4,double)
	1217	array_yl(0)=zs
	1218	array_xb(0)=ys
	1219	array_yl_labels(0)=z_d(zs)
	1220	array_xb_labels(0)=y_d(ys)
	1221	do ar=1,major_ticks_y-1
	1222	array_yl(ar)=d*(ar-1)+d-1
[329]	1223	array_yl_labels(ar) = z_d(array_yl(ar))
[218]	1224	if (ar .GT. 0)
	1225	do min_ar=0,3
[566]	1226	array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
	1227	(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
[218]	1228	end do
	1229	end if
	1230	end do
	1231	do br=1,major_ticks_x-1
	1232	array_xb(br)=e*(br-1)+e-1
[329]	1233	array_xb_labels(br) = y_d(array_xb(br))
[218]	1234	if (br .GT. 0)
	1235	do min_br=0,3
[566]	1236	array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+\
	1237	(y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1)
[218]	1238	end do
	1239	end if
	1240	end do
	1241	end if
[329]	1242
[218]	1243	if (axes_explicit .EQ. 1)then
	1244	cs_res@tmYLMode = "Explicit"
	1245	cs_res@tmXBMode = "Explicit"
[329]	1246	if (xyc .EQ. 1)then
	1247	cs_res@tmYLValues = y_d(array_yl)
	1248	cs_res@tmXBValues = x_d(array_xb)
[218]	1249	cs_res@tmYLLabels = array_yl_labels/norm_y
	1250	cs_res@tmXBLabels = array_xb_labels/norm_x
	1251	if (norm_x .NE. 1.)then
[534]	1252	cs_res@tiXAxisString = "x ("+unit_x+")"
[218]	1253	else
[534]	1254	cs_res@tiXAxisString = "x (m)"
[218]	1255	end if
	1256	if (norm_y .NE. 1.)then
[534]	1257	cs_res@tiYAxisString = "y ("+unit_y+")"
[218]	1258	else
[534]	1259	cs_res@tiYAxisString = "y (m)"
[218]	1260	end if
	1261	end if
	1262	if (xzc .EQ. 1)then
[329]	1263	cs_res@tmYLValues = z_d(array_yl)
	1264	cs_res@tmXBValues = x_d(array_xb)
[218]	1265	cs_res@tmYLLabels = array_yl_labels/norm_z
	1266	cs_res@tmXBLabels = array_xb_labels/norm_x
	1267	if (norm_x .NE. 1.)then
[534]	1268	cs_res@tiXAxisString = "x ("+unit_x+")"
[218]	1269	else
[534]	1270	cs_res@tiXAxisString = "x (m)"
[218]	1271	end if
	1272	if (norm_z .NE. 1.)then
[534]	1273	cs_res@tiYAxisString = "z ("+unit_z+")"
[218]	1274	else
[534]	1275	cs_res@tiYAxisString = "z (m)"
[218]	1276	end if
	1277	end if
	1278	if (yzc .EQ. 1)then
[329]	1279	cs_res@tmYLValues = z_d(array_yl)
	1280	cs_res@tmXBValues = y_d(array_xb)
[218]	1281	cs_res@tmYLLabels = array_yl_labels/norm_z
	1282	cs_res@tmXBLabels = array_xb_labels/norm_y
	1283	if (norm_y .NE. 1.)then
[534]	1284	cs_res@tiXAxisString = "y ("+unit_y+")"
[218]	1285	else
[534]	1286	cs_res@tiXAxisString = "y (m)"
[218]	1287	end if
	1288	if (norm_z .NE. 1.)then
[534]	1289	cs_res@tiYAxisString = "z ("+unit_z+")"
[218]	1290	else
[534]	1291	cs_res@tiYAxisString = "z (m)"
[218]	1292	end if
	1293	end if
	1294	cs_res@tmYLMinorValues = array_minor_yl
	1295	cs_res@tmXBMinorValues = array_minor_xb
	1296	else
	1297	if (axes_explicit .EQ. 0)then
	1298	cs_res@tmYLMinorPerMajor = 4
	1299	cs_res@tmXBMinorPerMajor = 4
	1300	else
	1301	print(" ")
	1302	print("'axes_explicit' is invalid and set to 0")
	1303	print(" ")
	1304	axes_explicit = 0
	1305	cs_res@tmYLMinorPerMajor = 4
	1306	cs_res@tmXBMinorPerMajor = 4
	1307	end if
	1308	if (xyc .EQ. 1)then
[534]	1309	cs_res@tiXAxisString = "x (m)"
	1310	cs_res@tiYAxisString = "y (m)"
[218]	1311	end if
	1312	if (xzc .EQ. 1)then
[534]	1313	cs_res@tiXAxisString = "x (m)"
	1314	cs_res@tiYAxisString = "z (m)"
[218]	1315	end if
	1316	if (yzc .EQ. 1)then
[534]	1317	cs_res@tiXAxisString = "y (m)"
	1318	cs_res@tiYAxisString = "z (m)"
[218]	1319	end if
	1320	end if
	1321
[162]	1322	if (xyc .EQ. 1 .OR. xzc .EQ.1)then
	1323	if (xe .EQ. xs+1) then
	1324	print(" ")
[566]	1325	print("range end for x-coordinate="+xe*delta_x+\
	1326	"m must be at least two")
	1327	print("more gridpoints(="+2*delta_x+"m) greater than start range="+\
	1328	xs*delta_x+"m)")
[162]	1329	print(" ")
	1330	exit
	1331	end if
	1332	end if
	1333	if (xyc .EQ. 1 .OR. yzc .EQ.1)then
	1334	if (ye .EQ. ys+1) then
	1335	print(" ")
[566]	1336	print("range end for y-coordinate="+ye*delta_y+\
	1337	"m must be at least two")
	1338	print("more gridpoints(="+2*delta_y+"m greater than start range="+\
	1339	ys*delta_y+"m)")
[162]	1340	print(" ")
	1341	exit
	1342	end if
	1343	end if
	1344	if (xzc .EQ. 1 .OR. yzc .EQ.1)then
	1345	if (ze .EQ. zs+1) then
	1346	print(" ")
[566]	1347	print("range end for z-coordinate="+ze+" must be at least two")
[190]	1348	print("more gridpoints greater than start range="+zs+")")
[162]	1349	print(" ")
	1350	exit
	1351	end if
	1352	end if
	1353
[161]	1354	if (xyc .EQ. 1) then
	1355	no_layer = (ze-zs)+1
[585]	1356	if (check_vType) then
	1357	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
	1358	else
	1359	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
	1360	end if
[161]	1361	end if
	1362	if (xzc .EQ. 1) then
	1363	no_layer = (ye-ys)+1
[585]	1364	if (check_vType) then
	1365	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
	1366	else
	1367	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
	1368	end if
[161]	1369	end if
	1370	if (yzc .EQ. 1) then
	1371	no_layer = (xe-xs)+1
[585]	1372	if (check_vType) then
	1373	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
	1374	else
	1375	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
	1376	end if
[161]	1377	end if
	1378
[585]	1379	if (check_vType) then
	1380	MinVal = new(dim,double)
	1381	MaxVal = new(dim,double)
	1382	else
	1383	MinVal = new(dim,float)
	1384	MaxVal = new(dim,float)
	1385	end if
[218]	1386	unit = new(dim,string)
[157]	1387
	1388	; ****************************************************
	1389	; define inner and outer loops depending on "sort"
	1390	; ****************************************************
	1391
	1392	if ( xyc .eq. 1 ) then
	1393	lays = zs
	1394	laye = ze
	1395	end if
	1396	if ( xzc .eq. 1 ) then
	1397	lays = ys
	1398	laye = ye
	1399	end if
	1400	if ( yzc .eq. 1) then
	1401	lays = xs
	1402	laye = xe
	1403	end if
	1404
	1405	if ( sort .eq. "time" ) then
	1406	lis = start_time_step
	1407	lie = end_time_step
[358]	1408	los = lays
	1409	loe = laye
[157]	1410	else
[358]	1411	lis = lays
	1412	lie = laye
[157]	1413	los = start_time_step
	1414	loe = end_time_step
	1415	end if
	1416
[758]	1417	check = True
	1418	v1 = 0
	1419	v2 = 0
	1420	no_var = 0
	1421	n = 0
	1422	no_zu1 = 0
[769]	1423	no_topo= 0
[157]	1424
[566]	1425	;****************************************************
	1426	; Preparation of vector plots
	1427	; ****************************************************
	1428
	1429	do varn=dim-1,0,1
	1430
	1431	if (vector .EQ. 1) then
	1432	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
	1433	if (check_vec1) then
	1434	temp = f[:]->$vNam(varn)$
	1435	data_att = f_att->$vNam(varn)$
	1436	vect1=temp(:,zs:ze,ys:ye,xs:xe)
	1437	delete(temp)
	1438	end if
	1439
	1440	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
	1441	if (check_vec2) then
	1442	temp = f[:]->$vNam(varn)$
	1443	data_att = f_att->$vNam(varn)$
	1444	vect2=temp(:,zs:ze,ys:ye,xs:xe)
	1445	delete(temp)
	1446	end if
[162]	1447
[566]	1448	if (","+vNam(varn)+"," .EQ. vec1) then
	1449	v1=v1+1
	1450	end if
	1451	if (","+vNam(varn)+"," .EQ. vec2) then
	1452	v2=v2+1
	1453	end if
	1454	end if
	1455
	1456	end do
	1457
	1458
	1459
	1460	do varn=dim-1,0,1
[157]	1461
[566]	1462	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
	1463	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
	1464	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
	1465	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
	1466	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
	1467	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
	1468	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
	1469	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
	1470	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
	1471	vNam(varn) .eq. "ind_x_yz") then
[157]	1472	check = False
[161]	1473	else
	1474	check = True
	1475	end if
[157]	1476
[357]	1477	if (var .NE. "all") then
	1478	check = isStrSubset( var,","+vNam(varn)+"," )
	1479	end if
	1480
[162]	1481	if(check) then
[194]	1482
[161]	1483	no_var=no_var+1
	1484
[162]	1485	if (xyc .EQ. 1) then
[218]	1486	temp = f[:]->$vNam(varn)$
[769]	1487	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
	1488	dummy=0
	1489	else
	1490	data_att = f_att->$vNam(varn)$
	1491	end if
[514]	1492	if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \
	1493	.or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
	1494	.or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \
[526]	1495	.or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy" \
	1496	.or. vNam(varn) .eq. "lwp_xy" .or. vNam(varn) .eq. "pra_xy" \
	1497	.or. vNam(varn) .eq. "prr_xy" .or. vNam(varn) .eq. "qsws_xy"\
	1498	.or. vNam(varn) .eq. "shf_xy" .or. vNam(varn) .eq. "t_xy" \
	1499	.or. vNam(varn) .eq. "u_xy" .or. vNam(varn) .eq. "z0_xy") then
[566]	1500	;these variables depend on zu1_xy and that's why they have
	1501	;only one z-layer
[329]	1502	data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
[357]	1503	no_zu1=no_zu1+1
[329]	1504	else
[769]	1505	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
	1506	;these variables depend on x and y
	1507	data(varn,0,0,:,:)=doubletofloat(temp(ys:ye,xs:xe))
	1508	no_topo=no_topo+1
	1509	else
	1510	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
	1511	end if
[329]	1512	end if
	1513	delete(temp)
[162]	1514	end if
	1515	if ( xzc .eq. 1 ) then
[218]	1516	temp = f[:]->$vNam(varn)$
	1517	data_att = f_att->$vNam(varn)$
	1518	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[357]	1519	delete(temp)
[162]	1520	end if
	1521	if ( yzc .eq. 1) then
[218]	1522	temp = f[:]->$vNam(varn)$
	1523	data_att = f_att->$vNam(varn)$
	1524	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[357]	1525	delete(temp)
[162]	1526	end if
[154]	1527
[157]	1528	data!0 = "var"
	1529	data!1 = "t"
	1530	data!2 = "z"
	1531	data!3 = "y"
[162]	1532	data!4 = "x"
[154]	1533
[566]	1534	MinVal(varn) = min(data(varn,start_time_step:end_time_step,\
	1535	0:(ze-zs),0:(ye-ys),0:(xe-xs)))
	1536	MaxVal(varn) = max(data(varn,start_time_step:end_time_step,\
	1537	0:(ze-zs),0:(ye-ys),0:(xe-xs)))
[218]	1538
[769]	1539	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
	1540	unit(varn) = "meters"
	1541	else
	1542	unit(varn) = data_att@units
	1543	delete(data_att)
	1544	end if
[162]	1545
[157]	1546	end if
[154]	1547
[190]	1548	end do
[154]	1549
[357]	1550	if (no_var .EQ. 0) then
	1551	print(" ")
	1552	print("The variables var='"+var+"' do not exist on your input file;")
	1553	print("be sure to have one comma before and after each variable")
	1554	print(" ")
	1555	exit
	1556	end if
	1557
[190]	1558	var_input=new(no_var,string)
	1559	numb=0
[194]	1560	no_var=0
	1561
[190]	1562	do varn=dim-1,0,1
	1563
[566]	1564	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
	1565	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
	1566	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
	1567	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
	1568	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
	1569	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
	1570	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
	1571	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
	1572	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
	1573	vNam(varn) .eq. "ind_x_yz") then
[190]	1574	check = False
	1575	else
	1576	check = True
	1577	end if
	1578
[357]	1579	if (var .NE. "all") then
	1580	check = isStrSubset( var,","+vNam(varn)+"," )
	1581	end if
	1582
[190]	1583	if(check) then
	1584	var_input(numb)=vNam(varn)
	1585	numb=numb+1
	1586	end if
	1587
[194]	1588	if (check) then
	1589	no_var = no_var+1
	1590	end if
	1591
[161]	1592	end do
[154]	1593
[190]	1594	if (no_var .EQ. 0) then
[162]	1595	print(" ")
[190]	1596	print("The variables var='"+var+"' do not exist on your input file;")
[357]	1597	print("be sure to have one comma before and after each variable")
[162]	1598	print(" ")
	1599	exit
	1600	end if
	1601
[157]	1602	if (vector .EQ. 1) then
	1603	if (v1 .EQ. 0)then
	1604	print(" ")
[566]	1605	print("Component 1 for the vector-plot ('vec1') must be one "+\
	1606	"of the variables on the input file:")
[190]	1607	print("- "+var_input)
[357]	1608	print("be sure to have one comma before and after the variable")
[157]	1609	print(" ")
	1610	exit
	1611	end if
[154]	1612
[157]	1613	if (v2 .EQ. 0)then
	1614	print(" ")
[566]	1615	print("Component 2 for the vector-plot ('vec2') must be one "+\
	1616	"of the variables on the input file:")
[190]	1617	print("- "+var_input)
[357]	1618	print("be sure to have one comma before and after the variable")
[157]	1619	print(" ")
	1620	exit
	1621	end if
	1622	end if
[154]	1623
[157]	1624	; ***************************************************
[161]	1625	; open workstation(s)
	1626	; ***************************************************
	1627
[532]	1628	if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
	1629	format_out@wkPaperSize = "A4"
	1630	end if
	1631	if ( format_out .EQ. "png" ) then
	1632	format_out@wkWidth = 1000
	1633	format_out@wkHeight = 1000
	1634	end if
	1635
[161]	1636	wks_ps = gsn_open_wks(format_out,file_out)
	1637	gsn_define_colormap(wks_ps,"rainbow+white")
	1638
[218]	1639	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[758]	1640	plot=new((/no_time*no_layer/),graphic)
[218]	1641	else
[769]	1642	plot=new((/no_timeno_layer(no_var-no_topo) + no_topo - \
	1643	no_time(no_layer-1)no_zu1/),graphic)
[218]	1644	end if
[357]	1645	dim_plot=dimsizes(plot)
[218]	1646
[161]	1647	page = 0
	1648	n=0
	1649	print(" ")
[162]	1650	print("Plots sorted by '"+sort+"'")
	1651	print(" ")
[161]	1652
	1653	; ***************************************************
[157]	1654	; create plots
	1655	; ***************************************************
	1656
[357]	1657
[190]	1658	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[162]	1659	do lo = los, loe
	1660	do li = lis, lie
[218]	1661	if (xyc .EQ. 1)then
	1662	if (sort .EQ. "time")then
[358]	1663	if(z_d(lo) .eq. -1.d) then
	1664	level = "z-average"
	1665	else
	1666	level = "z=" + z_d(lo) + "m"
	1667	end if
[218]	1668	else
[358]	1669	if(z_d(li) .eq. -1.d) then
	1670	level = "z-average"
	1671	else
	1672	level = "z=" + z_d(li) + "m"
	1673	end if
[218]	1674	end if
	1675	end if
	1676	if (xzc .EQ. 1)then
	1677	if (sort .EQ. "time")then
[358]	1678	if(y_d(lo) .eq. -1.d) then
	1679	level = "y-average"
	1680	else
	1681	level = "y=" + y_d(lo) + "m"
	1682	end if
[218]	1683	else
[358]	1684	if(y_d(li) .eq. -1.d) then
	1685	level = "y-average"
	1686	else
	1687	level = "y=" + y_d(li) + "m"
	1688	end if
[218]	1689	end if
	1690	end if
	1691	if (yzc .EQ. 1)then
	1692	if (sort .EQ. "time")then
[358]	1693	if(x_d(lo) .eq. -1.d) then
	1694	level = "x-average"
	1695	else
	1696	level = "x=" + x_d(lo) + "m"
	1697	end if
[218]	1698	else
[358]	1699	if(x_d(li) .eq. -1.d) then
	1700	level = "x-average"
	1701	else
	1702	level = "x=" + x_d(li) + "m"
	1703	end if
[218]	1704	end if
	1705	end if
[162]	1706	vecres = True ; vector only resources
	1707	vecres@gsnDraw = False ; don't draw
	1708	vecres@gsnFrame = False ; don't advance frame
	1709	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	1710	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
[566]	1711	vecres@vcRefLengthF = 0.05 ; define length of
	1712	; vec ref
[162]	1713	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
[218]	1714	vecres@tmXBLabelFontHeightF = font_size
	1715	vecres@tmYLLabelFontHeightF = font_size
	1716	vecres@tiXAxisFontHeightF = font_size
	1717	vecres@tiYAxisFontHeightF = font_size
	1718	if (sort .EQ. "time")then
[566]	1719	vecres@gsnLeftString = "t=" + \
	1720	decimalPlaces(t_all(li)/3600,2,True) +"h "+level
[218]	1721	else
[566]	1722	vecres@gsnLeftString = "t=" + \
	1723	decimalPlaces(t_all(lo)/3600,2,True) +"h "+level
[218]	1724	end if
[529]	1725	if (xyc .EQ. 1) then
[534]	1726	vecres@tiXAxisString = "x (m)"
	1727	vecres@tiYAxisString = "y (m)"
[566]	1728	if (sort .EQ. "time")then
	1729	plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
	1730	vect2(li,lo-los,:,:),vecres)
[218]	1731	else
[566]	1732	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
	1733	vect2(lo,li-lis,:,:),vecres)
[218]	1734	end if
	1735	end if
	1736	if (xzc .EQ. 1) then
[534]	1737	vecres@tiXAxisString = "x (m)"
	1738	vecres@tiYAxisString = "z (m)"
[218]	1739	if (sort .EQ. "time")then
[566]	1740	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
	1741	vect2(li,:,lo-los,:),vecres)
[218]	1742	else
[566]	1743	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
	1744	vect2(lo,:,li-lis,:),vecres)
[218]	1745	end if
	1746	end if
	1747	if (yzc .EQ. 1) then
[534]	1748	vecres@tiXAxisString = "y (m)"
	1749	vecres@tiYAxisString = "z (m)"
[218]	1750	if (sort .EQ. "time")then
[566]	1751	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
	1752	vect2(li,:,:,lo-los),vecres)
[218]	1753	else
[566]	1754	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
	1755	vect2(lo,:,:,li-lis),vecres)
[218]	1756	end if
	1757	end if
[162]	1758	n=n+1
	1759	end do
	1760	end do
	1761	end if
	1762
[161]	1763	do varn=dim-1,0,1
[157]	1764
[529]	1765	if (vector .EQ. 1 ) then
[566]	1766	check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
[157]	1767	end if
[358]	1768
[566]	1769	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
	1770	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
	1771	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
	1772	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
	1773	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
	1774	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
	1775	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
	1776	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
	1777	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
	1778	vNam(varn) .eq. "ind_x_yz") then
[161]	1779	check = False
	1780	else
	1781	check = True
	1782	end if
[190]	1783
	1784	if (var .NE. "all") then
[161]	1785	check = isStrSubset( var,","+vNam(varn)+"," )
[157]	1786	end if
	1787
	1788	if(check) then
[218]	1789
[162]	1790	space=(MaxVal(varn)-MinVal(varn))/24
	1791
	1792	cs_res@cnMinLevelValF = MinVal(varn)
	1793	cs_res@cnMaxLevelValF = MaxVal(varn)
	1794
	1795	cs_res@cnLevelSpacingF = space
	1796
[157]	1797	; ****************************************************
	1798	; loops over time and layer
	1799	; ****************************************************
	1800
[162]	1801
[161]	1802	do lo = los, loe
[162]	1803	do li = lis, lie
	1804
[157]	1805	; ****************************************************
[154]	1806	; xy cross section
[157]	1807	; ****************************************************
[154]	1808
[157]	1809	if ( xyc .eq. 1 ) then
[218]	1810
[157]	1811	cs_res@gsnLeftString = "Plot of "+vNam(varn)
[218]	1812	cs_res@gsnRightString = unit(varn)
[162]	1813
[154]	1814	if ( sort .eq. "time" ) then
[358]	1815	if ( z_d(lo) .eq. -1)then
[194]	1816	if (delta_z .EQ. -1) then
	1817	level = "-average"
	1818	else
[358]	1819	level = "=" + z_d(lo) + "m"
[194]	1820	end if
[154]	1821	else
[358]	1822	level = "=" + z_d(lo) + "m"
[154]	1823	end if
[357]	1824
[566]	1825	if(vNam(varn) .eq. "lwps_xy" .or. \
	1826	vNam(varn) .eq. "pras_xy" .or. \
	1827	vNam(varn) .eq. "prrs_xy" .or. \
	1828	vNam(varn) .eq. "qsws_xy" .or. \
	1829	vNam(varn) .eq. "shfs_xy" .or. \
	1830	vNam(varn) .eq. "ts_xy" .or. \
	1831	vNam(varn) .eq. "us_xy" .or. \
	1832	vNam(varn) .eq. "z0s_xy" .or. \
	1833	vNam(varn) .eq. "lwp*_xy" .or. \
	1834	vNam(varn) .eq. "pra*_xy" .or. \
	1835	vNam(varn) .eq. "prr*_xy" .or. \
	1836	vNam(varn) .eq. "qsws*_xy" .or. \
	1837	vNam(varn) .eq. "shf*_xy" .or. \
	1838	vNam(varn) .eq. "t*_xy" .or. \
	1839	vNam(varn) .eq. "u*_xy" .or. \
	1840	vNam(varn) .eq. "z0*_xy") then
[357]	1841	loe = 0
	1842	level = "=" + zu1(0) + "m"
[769]	1843	else
	1844	loe = laye
[357]	1845	end if
[769]	1846
	1847	if(vNam(varn) .eq. "zwwi" .or. \
	1848	vNam(varn) .eq. "zusi") then
	1849	loe = 0
	1850	los = 0
	1851	lie = 0
	1852	lis = 0
	1853	level = ""
	1854	else
	1855	lis = start_time_step
	1856	lie = end_time_step
	1857	los = lays
	1858	loe = laye
[770]	1859	end if
	1860
[769]	1861	if(vNam(varn) .eq. "zwwi" .or. \
	1862	vNam(varn) .eq. "zusi") then
	1863	cs_res@gsnCenterString = ""
	1864	else
	1865	cs_res@gsnCenterString = "t=" + \
	1866	decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
	1867	end if
[758]	1868
	1869	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	1870	;nothing to be done
	1871	else
	1872	plot(n) = gsn_csm_contour(wks_ps,\
[566]	1873	data(varn,li,lo-los,:,:),cs_res)
[758]	1874	end if
	1875
	1876
[162]	1877	if (vector .EQ. 1 .AND. check_vecp) then
[566]	1878	vecres = True ; vector only
	1879	; resources
	1880	vecres@gsnDraw = False ; don't draw
	1881	vecres@gsnFrame = False ; don't
	1882	; advance frame
	1883	vecres@vcGlyphStyle = "CurlyVector" ; curly
	1884	; vectors
	1885	vecres@vcRefMagnitudeF = ref_mag ; define
	1886	; vector ref
	1887	; mag
	1888	vecres@vcRefLengthF = 0.05 ; define
	1889	; length of
	1890	; vec ref
	1891	vecres@gsnRightString = " "
	1892	vecres@gsnLeftString = " "
	1893	vecres@tiXAxisString = " "
	1894	plot_vec=gsn_csm_vector(wks_ps,\
	1895	vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
[162]	1896	overlay(plot(n), plot_vec)
[161]	1897	end if
[154]	1898	end if
[162]	1899
[154]	1900	if ( sort .eq. "layer" ) then
[513]	1901	if (z_d(li) .eq. -1) then
[194]	1902	if (delta_z .EQ. -1) then
	1903	level = "-average"
	1904	else
[358]	1905	level = "=" + z_d(li) + "m"
[194]	1906	end if
[161]	1907	else
[358]	1908	level = "=" + z_d(li) + "m"
[190]	1909	end if
[357]	1910
[566]	1911	if(vNam(varn) .eq. "lwps_xy" .or. \
	1912	vNam(varn) .eq. "pras_xy" .or. \
	1913	vNam(varn) .eq. "prrs_xy" .or. \
	1914	vNam(varn) .eq. "qsws_xy" .or. \
	1915	vNam(varn) .eq. "shfs_xy" .or. \
	1916	vNam(varn) .eq. "ts_xy" .or. \
	1917	vNam(varn) .eq. "us_xy" .or. \
	1918	vNam(varn) .eq. "z0s_xy" .or. \
	1919	vNam(varn) .eq. "lwp*_xy" .or. \
	1920	vNam(varn) .eq. "pra*_xy" .or. \
	1921	vNam(varn) .eq. "prr*_xy" .or. \
	1922	vNam(varn) .eq. "qsws*_xy" .or. \
	1923	vNam(varn) .eq. "shf*_xy" .or. \
	1924	vNam(varn) .eq. "t*_xy" .or. \
	1925	vNam(varn) .eq. "u*_xy" .or. \
	1926	vNam(varn) .eq. "z0*_xy") then
[357]	1927	lie = 0
	1928	level = "=" + zu1(0) + "m"
[769]	1929	else
	1930	lie = laye
[357]	1931	end if
[758]	1932
[769]	1933	if(vNam(varn) .eq. "zwwi" .or. \
	1934	vNam(varn) .eq. "zusi") then
	1935	lie = 0
	1936	lis = 0
	1937	loe = 0
	1938	los = 0
	1939	level = ""
	1940	else
	1941	lis = lays
	1942	lie = laye
	1943	los = start_time_step
	1944	loe = end_time_step
	1945	end if
	1946
	1947	if(vNam(varn) .eq. "zwwi" .or. \
	1948	vNam(varn) .eq. "zusi") then
	1949	cs_res@gsnCenterString = ""
	1950	else
	1951	cs_res@gsnCenterString = "t=" + \
[770]	1952	decimalPlaces(t_all(lo)/3600,2,True) +"h z"+level
[769]	1953	end if
	1954
[758]	1955	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	1956	;nothing to be done
	1957	else
	1958	plot(n) = gsn_csm_contour(wks_ps,\
[566]	1959	data(varn,lo,li-lis,:,:),cs_res)
[758]	1960	end if
	1961
[162]	1962	if (vector .EQ. 1 .AND. check_vecp) then
[566]	1963	vecres = True ; vector only
	1964	; resources
	1965	vecres@gsnDraw = False ; don't draw
	1966	vecres@gsnFrame = False ; don't
	1967	; advance frame
	1968	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	1969	vecres@vcRefMagnitudeF = ref_mag ; define
	1970	; vector ref
	1971	; mag
	1972	vecres@vcRefLengthF = 0.05 ; define
	1973	; length of
	1974	; vec ref
	1975	vecres@gsnRightString = " " ; turn off
	1976	; right string
	1977	vecres@gsnLeftString = " " ; turn off
	1978	; left string
	1979	vecres@tiXAxisString = " "
	1980	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
	1981	vect2(lo,li-lis,:,:),vecres)
[162]	1982	overlay(plot(n), plot_vec)
[161]	1983	end if
[154]	1984	end if
	1985	end if
	1986
[157]	1987	; ****************************************************
[154]	1988	; xz cross section
[157]	1989	; ****************************************************
[154]	1990
	1991	if ( xzc .eq. 1 ) then
[218]	1992
[157]	1993	cs_res@gsnLeftString = "Plot of "+vNam(varn)
[162]	1994
[154]	1995	if ( sort .eq. "time" ) then
[358]	1996	if ( y_d(lo) .eq. -1 ) then
[154]	1997	level = "-average"
	1998	else
[358]	1999	level = "=" + y_d(lo) + "m"
[154]	2000	end if
[566]	2001	cs_res@gsnCenterString = "t=" + \
	2002	decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
[758]	2003
	2004	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2005	;nothing to be done
	2006	else
	2007	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2008	data(varn,li,:,lo-los,:),cs_res)
[758]	2009	end if
	2010
[162]	2011	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2012	vecres = True ; vector only
	2013	; resources
	2014	vecres@gsnDraw = False ; don't draw
	2015	vecres@gsnFrame = False ; don't
	2016	; advance frame
	2017	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2018	vecres@vcRefMagnitudeF = ref_mag ; define
	2019	; vector ref
	2020	; mag
	2021	vecres@vcRefLengthF = 0.05 ; define
	2022	; length of
	2023	; vec ref
	2024	vecres@gsnRightString = " " ; turn off
	2025	; right string
	2026	vecres@gsnLeftString = " " ; turn off
	2027	; left string
	2028	vecres@tiXAxisString = " "
	2029	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
	2030	vect2(li,:,lo-los,:),vecres)
[162]	2031	overlay(plot(n), plot_vec)
[161]	2032	end if
[154]	2033	end if
	2034
	2035	if ( sort .eq. "layer" ) then
[358]	2036	if ( y_d(li) .eq. -1 ) then
[161]	2037	level = "-average"
	2038	else
[358]	2039	level = "=" + y_d(li) + "m"
[161]	2040	end if
[566]	2041	cs_res@gsnCenterString = "t=" + \
	2042	decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
[758]	2043
	2044	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2045	;nothing to be done
	2046	else
	2047	plot(n) = gsn_csm_contour(wks_ps,\
	2048	data(varn,lo,:,li-lis,:),cs_res)
	2049	end if
	2050
[162]	2051	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2052	vecres = True ; vector only
	2053	; resources
	2054	vecres@gsnDraw = False ; don't draw
	2055	vecres@gsnFrame = False ; don't
	2056	; advance frame
	2057	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2058	vecres@vcRefMagnitudeF = ref_mag ; define
	2059	; vector ref
	2060	; mag
	2061	vecres@vcRefLengthF = 0.05 ; define
	2062	; length of
	2063	; vec ref
	2064	vecres@gsnRightString = " " ; turn off
	2065	; right string
	2066	vecres@gsnLeftString = " " ; turn off
	2067	; left string
	2068	vecres@tiXAxisString = " "
	2069	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
	2070	vect2(lo,:,li-lis,:),vecres)
[162]	2071	overlay(plot(n), plot_vec)
[161]	2072	end if
[157]	2073	end if
[154]	2074	end if
	2075
[157]	2076	; ****************************************************
[154]	2077	; yz cross section
[157]	2078	; ****************************************************
[154]	2079
	2080	if ( yzc .eq. 1 ) then
[218]	2081
[157]	2082	cs_res@gsnLeftString = "Plot of "+vNam(varn)
	2083
[154]	2084	if ( sort .eq. "time" ) then
[513]	2085	if ( x_d(lo) .eq. -1 ) then
[154]	2086	level = "-average"
	2087	else
[358]	2088	level = "=" + x_d(lo) + "m"
[154]	2089	end if
[566]	2090	cs_res@gsnCenterString = "t=" + \
	2091	decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
[758]	2092
	2093	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2094	;nothing to be done
	2095	else
	2096	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2097	data(varn,li,:,:,lo-los),cs_res)
[758]	2098	end if
	2099
[162]	2100	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2101	vecres = True ; vector only
	2102	; resources
	2103	vecres@gsnDraw = False ; don't draw
	2104	vecres@gsnFrame = False ; don't
	2105	; advance frame
	2106	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2107	vecres@vcRefMagnitudeF = ref_mag ; define
	2108	; vector ref
	2109	; mag
	2110	vecres@vcRefLengthF = 0.05 ; define
	2111	; length of
	2112	; vec ref
	2113	vecres@gsnRightString = " " ; turn off
	2114	; right string
	2115	vecres@gsnLeftString = " " ; turn off
	2116	; left string
	2117	vecres@tiXAxisString = " "
	2118	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
	2119	vect2(li,:,:,lo-los),vecres)
[162]	2120	overlay(plot(n), plot_vec)
[161]	2121	end if
[154]	2122	end if
	2123
	2124	if ( sort .eq. "layer" ) then
[513]	2125	if ( x_d(li) .eq. -1 ) then
[161]	2126	level = "-average"
	2127	else
[358]	2128	level = "=" + x_d(li) + "m"
[161]	2129	end if
[566]	2130	cs_res@gsnCenterString = "t=" + \
	2131	decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
[758]	2132
	2133	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2134	;nothing to be done
	2135	else
	2136	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2137	data(varn,lo,:,:,li-lis),cs_res)
[758]	2138	end if
	2139
[162]	2140	if (vector .EQ. 1 .AND. check_vecp)then
[566]	2141	vecres = True ; vector only
	2142	;resources
	2143	vecres@gsnDraw = False ; don't draw
	2144	vecres@gsnFrame = False ; don't
	2145	; advance frame
	2146	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2147	vecres@vcRefMagnitudeF = ref_mag ; define
	2148	; vector ref
	2149	; mag
	2150	vecres@vcRefLengthF = 0.05 ; define
	2151	; length of
	2152	; vec ref
	2153	vecres@gsnRightString = " " ; turn off
	2154	; right string
	2155	vecres@gsnLeftString = " " ; turn off
	2156	; left string
	2157	vecres@tiXAxisString = " "
	2158	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
	2159	vect2(lo,:,:,li-lis),vecres)
[162]	2160	overlay(plot(n), plot_vec)
[161]	2161	end if
[154]	2162	end if
[161]	2163	end if
[218]	2164	n=n+1
[154]	2165	end do
	2166	end do
[162]	2167	end if
[161]	2168	end do
[154]	2169
[161]	2170	; ***************************************************
	2171	; merge plots onto one page
[357]	2172	; ***************************************************
[532]	2173
	2174	no_frames = 0
	2175
[529]	2176	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[566]	2177	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
	2178	no_timeno_layer(no_var+1) .gt. no_rows*no_columns) then
	2179	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),\
	2180	(/no_var+1,no_layer*no_time/),cs_resP)
[162]	2181	print(" ")
	2182	print("Outputs to .eps or .epsi have only one frame")
	2183	print(" ")
	2184	else
[357]	2185	do np = 0,dim_plot-1,no_rows*no_columns
	2186	if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]	2187	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
	2188	cs_resP)
[532]	2189	no_frames = no_frames + 1
[162]	2190	else
[566]	2191	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
	2192	(/no_rows,no_columns/),cs_resP)
[532]	2193	no_frames = no_frames + 1
[162]	2194	end if
	2195	end do
	2196	end if
[267]	2197	else
[566]	2198	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
	2199	.AND. dim_plot .gt. no_rows*no_columns) then
[357]	2200	gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
[162]	2201	print(" ")
	2202	print("Outputs to .eps or .epsi have only one frame")
	2203	print(" ")
	2204	else
[357]	2205	do np = 0,dim_plot-1,no_rows*no_columns
	2206	if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]	2207	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
	2208	cs_resP)
[532]	2209	no_frames = no_frames + 1
[162]	2210	else
[566]	2211	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
	2212	(/no_rows,no_columns/),cs_resP)
[532]	2213	no_frames = no_frames + 1
[162]	2214	end if
	2215	end do
	2216	end if
[161]	2217	end if
	2218
[532]	2219	if (format_out .EQ. "png" ) then
	2220	png_output = new((/no_frames/), string)
	2221	j = 0
	2222	do i=0, no_frames-1
	2223	j = i + 1
	2224	if (j .LE. 9) then
	2225	png_output(i) = file_out+".00000"+j+".png"
	2226	end if
	2227	if (j .GT. 9 .AND. j .LE. 99) then
	2228	png_output(i) = file_out+".0000"+j+".png"
	2229	end if
	2230	if (j .GT. 99 .AND. j .LE. 999) then
	2231	png_output(i) = file_out+".000"+j+".png"
	2232	end if
	2233	if (j .GT. 999) then
	2234	png_output(i) = file_out+".00"+j+".png"
	2235	end if
	2236
	2237	;using imagemagick's convert for reducing the white
	2238	;space around the plot
	2239	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
	2240	png_output(i) + " " + png_output(i)
	2241	system(cmd)
	2242	end do
	2243
	2244	print(" ")
	2245	print("Output to: "+ png_output)
	2246	print(" ")
	2247	else
	2248	print(" ")
	2249	print("Output to: " + file_out +"."+ format_out)
	2250	print(" ")
	2251	end if
[154]	2252
[190]	2253	end

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

Download in other formats:

Original Format