Home

Context Navigation

cross_sections.ncl @ 1621

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

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

Download in other formats:

Original Format