Home

Context Navigation

cross_sections.ncl @ 526

Last change on this file since 526 was 526, checked in by heinze, 14 years ago
Adjustment of the NCL scripts and palmplot to allow for the use of special characters in NetCDF variable names
File size: 64.2 KB

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

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

Download in other formats:

Original Format