Home

Context Navigation

cross_sections.ncl @ 566

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

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

Download in other formats:

Original Format