Home

Context Navigation

cross_sections.ncl @ 794

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

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

Download in other formats:

Original Format