Home

Context Navigation

cross_sections.ncl @ 195

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

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

Download in other formats:

Original Format