Home

Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

cross_sections.ncl

Last change on this file was 2837, checked in by gronemeier, 6 years ago
palmplot: config file of ncl scripts can be chosen in shell command
File size: 79.2 KB

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

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

Download in other formats:

Original Format