Home

Context Navigation

cross_sections.ncl @ 710

Last change on this file since 710 was 585, checked in by heinze, 14 years ago
Bugfix: enable plot of data if it is of kind double instead of kind float
File size: 74.8 KB

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

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

Download in other formats:

Original Format