Home

Context Navigation

cross_sections.ncl @ 416

Last change on this file since 416 was 358, checked in by heinze, 15 years ago
Bugfix in cross_sections.ncl concerning plot layers and vector plots fixed.
File size: 62.0 KB

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

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

Download in other formats:

Original Format