Home

Context Navigation

cross_sections.ncl @ 170

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

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

Download in other formats:

Original Format