Home

Context Navigation

cross_sections.ncl @ 1317

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

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

Download in other formats:

Original Format