Home

Context Navigation

cross_sections.ncl @ 518

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

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

Download in other formats:

Original Format