Home

Context Navigation

cross_sections.ncl @ 537

Last change on this file since 537 was 534, checked in by heinze, 15 years ago
Bugfix in spectra.ncl concerning output in png and small changes in layout
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
[534]	1158	cs_res@tiXAxisString = "x ("+unit_x+")"
[218]	1159	else
[534]	1160	cs_res@tiXAxisString = "x (m)"
[218]	1161	end if
	1162	if (norm_y .NE. 1.)then
[534]	1163	cs_res@tiYAxisString = "y ("+unit_y+")"
[218]	1164	else
[534]	1165	cs_res@tiYAxisString = "y (m)"
[218]	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
[534]	1174	cs_res@tiXAxisString = "x ("+unit_x+")"
[218]	1175	else
[534]	1176	cs_res@tiXAxisString = "x (m)"
[218]	1177	end if
	1178	if (norm_z .NE. 1.)then
[534]	1179	cs_res@tiYAxisString = "z ("+unit_z+")"
[218]	1180	else
[534]	1181	cs_res@tiYAxisString = "z (m)"
[218]	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
[534]	1190	cs_res@tiXAxisString = "y ("+unit_y+")"
[218]	1191	else
[534]	1192	cs_res@tiXAxisString = "y (m)"
[218]	1193	end if
	1194	if (norm_z .NE. 1.)then
[534]	1195	cs_res@tiYAxisString = "z ("+unit_z+")"
[218]	1196	else
[534]	1197	cs_res@tiYAxisString = "z (m)"
[218]	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
[534]	1215	cs_res@tiXAxisString = "x (m)"
	1216	cs_res@tiYAxisString = "y (m)"
[218]	1217	end if
	1218	if (xzc .EQ. 1)then
[534]	1219	cs_res@tiXAxisString = "x (m)"
	1220	cs_res@tiYAxisString = "z (m)"
[218]	1221	end if
	1222	if (yzc .EQ. 1)then
[534]	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
[534]	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
[534]	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
[534]	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