Home

Context Navigation

cross_sections.ncl @ 2480

Last change on this file since 2480 was 2147, checked in by scharf, 8 years ago
improved palmplot, added two imuk hosts, changed allocation of resources for lcbullhh and corrected land surface parameters
File size: 79.5 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
[2030]	32	if( version .GE. "6.2.1" ) then
[1559]	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
[2147]	418	cs_res@cnLevelSelectionMode = "AutomaticLevels"
[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
[2147]	1821	space=(decimalPlaces(MaxVal(varn),3,True)-decimalPlaces(MinVal(varn),3,True))/24
[162]	1822
[2147]	1823	;cs_res@cnMinLevelValF = decimalPlaces(MinVal(varn),3,True)
	1824	;cs_res@cnMaxLevelValF = decimalPlaces(MaxVal(varn),3,True)
[162]	1825
[2147]	1826	;cs_res@cnLevelSpacingF = space
	1827	cs_res@cnMaxLevelCount = 26
[162]	1828
[157]	1829	; ****************************************************
	1830	; loops over time and layer
	1831	; ****************************************************
	1832
[162]	1833
[161]	1834	do lo = los, loe
[162]	1835	do li = lis, lie
	1836
[157]	1837	; ****************************************************
[154]	1838	; xy cross section
[157]	1839	; ****************************************************
[154]	1840
[157]	1841	if ( xyc .eq. 1 ) then
[218]	1842
[157]	1843	cs_res@gsnLeftString = "Plot of "+vNam(varn)
[218]	1844	cs_res@gsnRightString = unit(varn)
[162]	1845
[154]	1846	if ( sort .eq. "time" ) then
[358]	1847	if ( z_d(lo) .eq. -1)then
[194]	1848	if (delta_z .EQ. -1) then
	1849	level = "-average"
	1850	else
[358]	1851	level = "=" + z_d(lo) + "m"
[194]	1852	end if
[154]	1853	else
[358]	1854	level = "=" + z_d(lo) + "m"
[154]	1855	end if
[357]	1856
[1626]	1857	if(vNam(varn) .eq. "zwwi" .or. \
	1858	vNam(varn) .eq. "zusi") then
	1859	loe = 0
	1860	los = 0
	1861	lie = 0
	1862	lis = 0
	1863	level = ""
	1864	end if
	1865
[566]	1866	if(vNam(varn) .eq. "lwps_xy" .or. \
	1867	vNam(varn) .eq. "pras_xy" .or. \
	1868	vNam(varn) .eq. "prrs_xy" .or. \
	1869	vNam(varn) .eq. "qsws_xy" .or. \
	1870	vNam(varn) .eq. "shfs_xy" .or. \
	1871	vNam(varn) .eq. "ts_xy" .or. \
	1872	vNam(varn) .eq. "us_xy" .or. \
	1873	vNam(varn) .eq. "z0s_xy" .or. \
	1874	vNam(varn) .eq. "lwp*_xy" .or. \
	1875	vNam(varn) .eq. "pra*_xy" .or. \
	1876	vNam(varn) .eq. "prr*_xy" .or. \
	1877	vNam(varn) .eq. "qsws*_xy" .or. \
	1878	vNam(varn) .eq. "shf*_xy" .or. \
	1879	vNam(varn) .eq. "t*_xy" .or. \
	1880	vNam(varn) .eq. "u*_xy" .or. \
	1881	vNam(varn) .eq. "z0*_xy") then
[357]	1882	loe = 0
	1883	level = "=" + zu1(0) + "m"
[769]	1884	else
	1885	lis = start_time_step
	1886	lie = end_time_step
	1887	los = lays
	1888	loe = laye
[770]	1889	end if
	1890
[769]	1891	if(vNam(varn) .eq. "zwwi" .or. \
	1892	vNam(varn) .eq. "zusi") then
	1893	cs_res@gsnCenterString = ""
	1894	else
	1895	cs_res@gsnCenterString = "t=" + \
	1896	decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
	1897	end if
[758]	1898
	1899	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	1900	;nothing to be done
	1901	else
	1902	plot(n) = gsn_csm_contour(wks_ps,\
[566]	1903	data(varn,li,lo-los,:,:),cs_res)
[758]	1904	end if
	1905
	1906
[162]	1907	if (vector .EQ. 1 .AND. check_vecp) then
[566]	1908	vecres = True ; vector only
	1909	; resources
	1910	vecres@gsnDraw = False ; don't draw
	1911	vecres@gsnFrame = False ; don't
	1912	; advance frame
	1913	vecres@vcGlyphStyle = "CurlyVector" ; curly
	1914	; vectors
	1915	vecres@vcRefMagnitudeF = ref_mag ; define
	1916	; vector ref
	1917	; mag
	1918	vecres@vcRefLengthF = 0.05 ; define
	1919	; length of
	1920	; vec ref
	1921	vecres@gsnRightString = " "
	1922	vecres@gsnLeftString = " "
	1923	vecres@tiXAxisString = " "
	1924	plot_vec=gsn_csm_vector(wks_ps,\
	1925	vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
[162]	1926	overlay(plot(n), plot_vec)
[161]	1927	end if
[154]	1928	end if
[162]	1929
[154]	1930	if ( sort .eq. "layer" ) then
[513]	1931	if (z_d(li) .eq. -1) then
[194]	1932	if (delta_z .EQ. -1) then
	1933	level = "-average"
	1934	else
[358]	1935	level = "=" + z_d(li) + "m"
[194]	1936	end if
[161]	1937	else
[358]	1938	level = "=" + z_d(li) + "m"
[190]	1939	end if
[1626]	1940
	1941	if(vNam(varn) .eq. "zwwi" .or. \
	1942	vNam(varn) .eq. "zusi") then
	1943	lie = 0
	1944	lis = 0
	1945	loe = 0
	1946	los = 0
	1947	level = ""
	1948	end if
	1949
[566]	1950	if(vNam(varn) .eq. "lwps_xy" .or. \
	1951	vNam(varn) .eq. "pras_xy" .or. \
	1952	vNam(varn) .eq. "prrs_xy" .or. \
	1953	vNam(varn) .eq. "qsws_xy" .or. \
	1954	vNam(varn) .eq. "shfs_xy" .or. \
	1955	vNam(varn) .eq. "ts_xy" .or. \
	1956	vNam(varn) .eq. "us_xy" .or. \
	1957	vNam(varn) .eq. "z0s_xy" .or. \
	1958	vNam(varn) .eq. "lwp*_xy" .or. \
	1959	vNam(varn) .eq. "pra*_xy" .or. \
	1960	vNam(varn) .eq. "prr*_xy" .or. \
	1961	vNam(varn) .eq. "qsws*_xy" .or. \
	1962	vNam(varn) .eq. "shf*_xy" .or. \
	1963	vNam(varn) .eq. "t*_xy" .or. \
	1964	vNam(varn) .eq. "u*_xy" .or. \
	1965	vNam(varn) .eq. "z0*_xy") then
[357]	1966	lie = 0
	1967	level = "=" + zu1(0) + "m"
[769]	1968	else
	1969	lis = lays
	1970	lie = laye
	1971	los = start_time_step
	1972	loe = end_time_step
	1973	end if
	1974
	1975	if(vNam(varn) .eq. "zwwi" .or. \
	1976	vNam(varn) .eq. "zusi") then
	1977	cs_res@gsnCenterString = ""
	1978	else
	1979	cs_res@gsnCenterString = "t=" + \
[770]	1980	decimalPlaces(t_all(lo)/3600,2,True) +"h z"+level
[769]	1981	end if
	1982
[2147]	1983	;if (topo .EQ. 0) then
	1984	;nothing to be done
	1985	;else
	1986	; print("'topo' is set to " + topo)
	1987
	1988	;plot(n) = gsn_csm_contour(wks_ps,\
	1989	;data(varn,lo,li-lis,:,:),cs_res)
	1990	;if(vNam(varn) .EQ. "zwwi" .OR. \
	1991	; vNam(varn) .EQ. "zusi") then
	1992	; to_res@cnFillPalette = "GMT_gray"
	1993	;plot_topo = gsn_csm_contour(wks_ps,\
	1994	; data(varn,lo,li-lis,:,:),to_res)
	1995	;end if
	1996	;overlay(plot(n), plot_topo)
	1997	;end if
	1998
[758]	1999	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2000	;nothing to be done
	2001	else
	2002	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2003	data(varn,lo,li-lis,:,:),cs_res)
[758]	2004	end if
	2005
[162]	2006	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2007	vecres = True ; vector only
	2008	; resources
	2009	vecres@gsnDraw = False ; don't draw
	2010	vecres@gsnFrame = False ; don't
	2011	; advance frame
	2012	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2013	vecres@vcRefMagnitudeF = ref_mag ; define
	2014	; vector ref
	2015	; mag
	2016	vecres@vcRefLengthF = 0.05 ; define
	2017	; length of
	2018	; vec ref
	2019	vecres@gsnRightString = " " ; turn off
	2020	; right string
	2021	vecres@gsnLeftString = " " ; turn off
	2022	; left string
	2023	vecres@tiXAxisString = " "
	2024	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
	2025	vect2(lo,li-lis,:,:),vecres)
[162]	2026	overlay(plot(n), plot_vec)
[161]	2027	end if
[154]	2028	end if
	2029	end if
	2030
[157]	2031	; ****************************************************
[154]	2032	; xz cross section
[157]	2033	; ****************************************************
[154]	2034
	2035	if ( xzc .eq. 1 ) then
[218]	2036
[157]	2037	cs_res@gsnLeftString = "Plot of "+vNam(varn)
[162]	2038
[154]	2039	if ( sort .eq. "time" ) then
[358]	2040	if ( y_d(lo) .eq. -1 ) then
[154]	2041	level = "-average"
	2042	else
[358]	2043	level = "=" + y_d(lo) + "m"
[154]	2044	end if
[566]	2045	cs_res@gsnCenterString = "t=" + \
	2046	decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
[758]	2047
	2048	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2049	;nothing to be done
	2050	else
	2051	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2052	data(varn,li,:,lo-los,:),cs_res)
[758]	2053	end if
	2054
[162]	2055	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2056	vecres = True ; vector only
	2057	; resources
	2058	vecres@gsnDraw = False ; don't draw
	2059	vecres@gsnFrame = False ; don't
	2060	; advance frame
	2061	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2062	vecres@vcRefMagnitudeF = ref_mag ; define
	2063	; vector ref
	2064	; mag
	2065	vecres@vcRefLengthF = 0.05 ; define
	2066	; length of
	2067	; vec ref
	2068	vecres@gsnRightString = " " ; turn off
	2069	; right string
	2070	vecres@gsnLeftString = " " ; turn off
	2071	; left string
	2072	vecres@tiXAxisString = " "
	2073	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
	2074	vect2(li,:,lo-los,:),vecres)
[162]	2075	overlay(plot(n), plot_vec)
[161]	2076	end if
[154]	2077	end if
	2078
	2079	if ( sort .eq. "layer" ) then
[358]	2080	if ( y_d(li) .eq. -1 ) then
[161]	2081	level = "-average"
	2082	else
[358]	2083	level = "=" + y_d(li) + "m"
[161]	2084	end if
[566]	2085	cs_res@gsnCenterString = "t=" + \
	2086	decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
[758]	2087
	2088	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2089	;nothing to be done
	2090	else
	2091	plot(n) = gsn_csm_contour(wks_ps,\
	2092	data(varn,lo,:,li-lis,:),cs_res)
	2093	end if
	2094
[162]	2095	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2096	vecres = True ; vector only
	2097	; resources
	2098	vecres@gsnDraw = False ; don't draw
	2099	vecres@gsnFrame = False ; don't
	2100	; advance frame
	2101	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2102	vecres@vcRefMagnitudeF = ref_mag ; define
	2103	; vector ref
	2104	; mag
	2105	vecres@vcRefLengthF = 0.05 ; define
	2106	; length of
	2107	; vec ref
	2108	vecres@gsnRightString = " " ; turn off
	2109	; right string
	2110	vecres@gsnLeftString = " " ; turn off
	2111	; left string
	2112	vecres@tiXAxisString = " "
	2113	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
	2114	vect2(lo,:,li-lis,:),vecres)
[162]	2115	overlay(plot(n), plot_vec)
[161]	2116	end if
[157]	2117	end if
[154]	2118	end if
	2119
[157]	2120	; ****************************************************
[154]	2121	; yz cross section
[157]	2122	; ****************************************************
[154]	2123
	2124	if ( yzc .eq. 1 ) then
[218]	2125
[157]	2126	cs_res@gsnLeftString = "Plot of "+vNam(varn)
	2127
[154]	2128	if ( sort .eq. "time" ) then
[513]	2129	if ( x_d(lo) .eq. -1 ) then
[154]	2130	level = "-average"
	2131	else
[358]	2132	level = "=" + x_d(lo) + "m"
[154]	2133	end if
[566]	2134	cs_res@gsnCenterString = "t=" + \
	2135	decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
[758]	2136
	2137	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2138	;nothing to be done
	2139	else
	2140	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2141	data(varn,li,:,:,lo-los),cs_res)
[758]	2142	end if
	2143
[162]	2144	if (vector .EQ. 1 .AND. check_vecp) then
[566]	2145	vecres = True ; vector only
	2146	; resources
	2147	vecres@gsnDraw = False ; don't draw
	2148	vecres@gsnFrame = False ; don't
	2149	; advance frame
	2150	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2151	vecres@vcRefMagnitudeF = ref_mag ; define
	2152	; vector ref
	2153	; mag
	2154	vecres@vcRefLengthF = 0.05 ; define
	2155	; length of
	2156	; vec ref
	2157	vecres@gsnRightString = " " ; turn off
	2158	; right string
	2159	vecres@gsnLeftString = " " ; turn off
	2160	; left string
	2161	vecres@tiXAxisString = " "
	2162	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
	2163	vect2(li,:,:,lo-los),vecres)
[162]	2164	overlay(plot(n), plot_vec)
[161]	2165	end if
[154]	2166	end if
	2167
	2168	if ( sort .eq. "layer" ) then
[513]	2169	if ( x_d(li) .eq. -1 ) then
[161]	2170	level = "-average"
	2171	else
[358]	2172	level = "=" + x_d(li) + "m"
[161]	2173	end if
[566]	2174	cs_res@gsnCenterString = "t=" + \
	2175	decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
[758]	2176
	2177	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
	2178	;nothing to be done
	2179	else
	2180	plot(n) = gsn_csm_contour(wks_ps,\
[566]	2181	data(varn,lo,:,:,li-lis),cs_res)
[758]	2182	end if
	2183
[162]	2184	if (vector .EQ. 1 .AND. check_vecp)then
[566]	2185	vecres = True ; vector only
	2186	;resources
	2187	vecres@gsnDraw = False ; don't draw
	2188	vecres@gsnFrame = False ; don't
	2189	; advance frame
	2190	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	2191	vecres@vcRefMagnitudeF = ref_mag ; define
	2192	; vector ref
	2193	; mag
	2194	vecres@vcRefLengthF = 0.05 ; define
	2195	; length of
	2196	; vec ref
	2197	vecres@gsnRightString = " " ; turn off
	2198	; right string
	2199	vecres@gsnLeftString = " " ; turn off
	2200	; left string
	2201	vecres@tiXAxisString = " "
	2202	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
	2203	vect2(lo,:,:,li-lis),vecres)
[162]	2204	overlay(plot(n), plot_vec)
[161]	2205	end if
[154]	2206	end if
[161]	2207	end if
[218]	2208	n=n+1
[154]	2209	end do
	2210	end do
[162]	2211	end if
[161]	2212	end do
[154]	2213
[161]	2214	; ***************************************************
	2215	; merge plots onto one page
[357]	2216	; ***************************************************
[532]	2217
	2218	no_frames = 0
	2219
[529]	2220	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[566]	2221	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
	2222	no_timeno_layer(no_var+1) .gt. no_rows*no_columns) then
	2223	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),\
	2224	(/no_var+1,no_layer*no_time/),cs_resP)
[162]	2225	print(" ")
	2226	print("Outputs to .eps or .epsi have only one frame")
	2227	print(" ")
	2228	else
[357]	2229	do np = 0,dim_plot-1,no_rows*no_columns
	2230	if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]	2231	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
	2232	cs_resP)
[532]	2233	no_frames = no_frames + 1
[162]	2234	else
[566]	2235	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
	2236	(/no_rows,no_columns/),cs_resP)
[532]	2237	no_frames = no_frames + 1
[162]	2238	end if
	2239	end do
	2240	end if
[267]	2241	else
[566]	2242	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
	2243	.AND. dim_plot .gt. no_rows*no_columns) then
[357]	2244	gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
[162]	2245	print(" ")
	2246	print("Outputs to .eps or .epsi have only one frame")
	2247	print(" ")
	2248	else
[357]	2249	do np = 0,dim_plot-1,no_rows*no_columns
	2250	if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]	2251	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
	2252	cs_resP)
[532]	2253	no_frames = no_frames + 1
[162]	2254	else
[566]	2255	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
	2256	(/no_rows,no_columns/),cs_resP)
[532]	2257	no_frames = no_frames + 1
[162]	2258	end if
	2259	end do
	2260	end if
[161]	2261	end if
	2262
[532]	2263	if (format_out .EQ. "png" ) then
	2264	png_output = new((/no_frames/), string)
	2265	j = 0
	2266
[957]	2267	if (no_frames .eq. 1) then
	2268	if (ncl_version .GE. 6 ) then
	2269	png_output(0) = file_out+".png"
	2270	else
	2271	png_output(0) = file_out+".00000"+1+".png"
	2272	end if
	2273	;using imagemagick's convert for reducing the white
	2274	;space around the plot
	2275	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
	2276	png_output(0) + " " + png_output(0)
[532]	2277	system(cmd)
[957]	2278	else
[532]	2279
[957]	2280	do i=0, no_frames-1
	2281	j = i + 1
	2282	if (j .LE. 9) then
	2283	png_output(i) = file_out+".00000"+j+".png"
	2284	end if
	2285	if (j .GT. 9 .AND. j .LE. 99) then
	2286	png_output(i) = file_out+".0000"+j+".png"
	2287	end if
	2288	if (j .GT. 99 .AND. j .LE. 999) then
	2289	png_output(i) = file_out+".000"+j+".png"
	2290	end if
	2291	if (j .GT. 999) then
	2292	png_output(i) = file_out+".00"+j+".png"
	2293	end if
	2294
	2295	;using imagemagick's convert for reducing the white
	2296	;space around the plot
	2297	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
	2298	png_output(i) + " " + png_output(i)
	2299	system(cmd)
	2300	end do
	2301	end if
	2302
[532]	2303	print(" ")
	2304	print("Output to: "+ png_output)
	2305	print(" ")
	2306	else
	2307	print(" ")
	2308	print("Output to: " + file_out +"."+ format_out)
	2309	print(" ")
	2310	end if
[957]	2311
[190]	2312	end

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

Download in other formats:

Original Format