Home

Context Navigation

cross_sections.ncl @ 1091

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

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

Download in other formats:

Original Format