Home

Context Navigation

cross_sections.ncl @ 532

Last change on this file since 532 was 532, checked in by heinze, 14 years ago
NCL scripts allow the output of png files
File size: 66.3 KB

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

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

Download in other formats:

Original Format