Home

Context Navigation

cross_sections.ncl @ 494

Last change on this file since 494 was 418, checked in by heinze, 15 years ago
ncl_preferences is renamed .ncl.config.default
File size: 62.2 KB

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

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

Download in other formats:

Original Format