Home

Context Navigation

cross_sections.ncl @ 182

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

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

Download in other formats:

Original Format