Home

Context Navigation

cross_sections.ncl @ 161

Last change on this file since 161 was 161, checked in by letzel, 16 years ago
NCL scripts in trunk/SCRIPTS/NCL updated for vector plots
File size: 38.7 KB

Rev	Line
[154]	1	load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
[157]	2	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
[154]	3	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
	4	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
	5	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
	6
	7	begin
	8
	9	; ***************************************************
	10	; read parameter_list
	11	; ***************************************************
	12
	13	if (isfilepresent("~/.ncl_preferences")) then
[161]	14	parameter = asciiread("~/.ncl_preferences",75,"string")
[154]	15	delete(parameter@_FillValue)
	16	else
	17	print(" ")
	18	print("Please copy '.ncl_preferences' into your $home dircetory")
	19	print(" ")
	20	exit
	21	end if
	22
	23	; ***************************************************
	24	; set default parameter values and strings if not assigned in prompt or parameter list
	25	; ***************************************************
	26
	27	if ( .not. isvar("file_in") ) then ; path+name of input file
	28	if (parameter(7) .EQ. "input file") then
	29	print(" ")
	30	print("Please provide input file 'file_in = ' either in prompt or parameter_list")
	31	print(" ")
	32	exit
	33	else
	34	file_in = parameter(7)
	35	end if
	36	end if
	37	if ( .not. isvar("format_out") ) then ; format of output file
	38	format_out = "x11"
	39	if (parameter(9) .NE. "x11") then
	40	format_out = parameter(9)
	41	end if
	42	end if
	43	if ( .not. isvar("file_out") ) then ; path+name of output file
	44	file_out = "test"
	45	if (parameter(11) .NE. "test") then
	46	file_out = parameter(11)
	47	end if
	48	end if
	49	if ( .not. isvar("no_columns") ) then ; number of plots in one row
	50	no_columns = 1
	51	if (parameter(17) .NE. "1") then
	52	no_columns = stringtointeger(parameter(17))
	53	end if
	54	end if
	55	if ( .not. isvar("no_lines") ) then ; number of plot-lines on one sheet
	56	no_lines = 2
	57	if (parameter(19) .NE. "2") then
	58	no_lines = stringtointeger(parameter(19))
	59	end if
	60	end if
	61	if ( .not. isvar("sort") ) then ; sort of plots
	62	sort = "time"
	63	if (parameter(37) .NE. "time") then
	64	sort = parameter(37)
	65	end if
	66	end if
	67	if ( .not. isvar("mode") ) then ; pattern of contour plots
	68	mode = "Fill"
	69	if (parameter(39) .NE. "Fill") then
	70	mode = parameter(39)
	71	end if
	72	end if
	73	if ( .not. isvar("fill_mode") ) then ; pattern of filling
	74	fill_mode = "AreaFill"
	75	if (parameter(41) .NE. "AreaFill") then
	76	mode = parameter(41)
	77	end if
	78	end if
	79	if ( .not. isvar("shape") ) then ; keeping of aspect ratio
	80	shape = 1
	81	if (parameter(43) .NE. "1") then
	82	shape = stringtointeger(parameter(43))
	83	if (stringtointeger(parameter(43)) .NE. 0) then
	84	print(" ")
	85	print("Please set 'shape' to 0 or 1")
	86	print(" ")
	87	exit
	88	end if
	89	end if
	90	end if
	91	if ( .not. isvar("var") ) then ; set output of all variables
	92	check = True
	93	end if
	94
	95	if ( .not. isvar("xyc") ) then ; turn xy cross-section on or off
	96	xyc = 0
	97	if (parameter(45) .NE. "0") then
	98	xyc = stringtointeger(parameter(45))
	99	if (stringtointeger(parameter(45)) .NE. 1) then
	100	print(" ")
	101	print("Please set 'xyc' to 0 or 1")
	102	print(" ")
	103	exit
	104	end if
	105	end if
	106	end if
	107	if ( .not. isvar("xzc") ) then ; turn xz cross-section on or off
	108	xzc = 0
	109	if (parameter(47) .NE. "0") then
	110	xzc = stringtointeger(parameter(47))
	111	if (stringtointeger(parameter(47)) .NE. 1) then
	112	print(" ")
	113	print("Please set 'xzc' to 0 or 1")
	114	print(" ")
	115	exit
	116	end if
	117	end if
	118	end if
	119	if ( .not. isvar("yzc") ) then ; turn yz cross-section on or off
	120	yzc = 0
	121	if (parameter(49) .NE. "0") then
	122	yzc = stringtointeger(parameter(49))
	123	if (stringtointeger(parameter(49)) .NE. 1) then
	124	print(" ")
	125	print("Please set 'yzc' to 0 or 1")
	126	print(" ")
	127	exit
	128	end if
	129	end if
	130	end if
	131	if ( .not. isvar("xs") ) then ; start of x-coordinate range
	132	xs = 0
	133	if (parameter(51) .NE. "0") then
	134	xs = stringtointeger(parameter(51))
	135	end if
	136	end if
	137	if ( .not. isvar("ys") ) then ; start of y-coordinate range
	138	ys = 0
	139	if (parameter(55) .NE. "0") then
	140	ys = stringtointeger(parameter(55))
	141	end if
	142	end if
	143	if ( .not. isvar("zs") ) then ; start of z-coordinate range
	144	zs = 0
	145	if (parameter(59) .NE. "0") then
	146	zs = stringtointeger(parameter(59))
	147	end if
	148	end if
	149
	150	if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
	151	print(" ")
	152	print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)")
	153	print(" ")
	154	exit
	155	end if
	156	if (xyc .EQ. 1 ) then
	157	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
	158	print(" ")
	159	print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
	160	print(" ")
	161	exit
	162	end if
	163	end if
	164	if (xzc .EQ. 1) then
	165	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
	166	print(" ")
	167	print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
	168	print(" ")
	169	exit
	170	end if
	171	end if
	172	if (yzc .EQ. 1) then
	173	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
	174	print(" ")
	175	print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
	176	print(" ")
	177	exit
	178	end if
	179	end if
[157]	180	if ( .not. isvar("vector") ) then ; sort of plots
	181	vector = 0
	182	if (parameter(63) .NE. "0") then
	183	vector = stringtointeger(parameter(63))
	184	if (stringtointeger(parameter(63)) .NE. 1) then
	185	print(" ")
	186	print("Please set 'vector' to 0 or 1")
	187	print(" ")
	188	exit
	189	end if
	190	end if
	191	end if
	192	if ( .not. isvar("ref_mag") ) then ; sort of plots
	193	ref_mag = 0.05
	194	if (parameter(71) .NE. "0.05") then
	195	ref_mag = stringtofloat(parameter(71))
	196	end if
	197	end if
[154]	198
	199	; ***************************************************
	200	; open input file
	201	; ***************************************************
	202
	203	f = addfile( file_in, "r" )
	204
	205	vNam = getfilevarnames(f)
	206	print(" ")
	207	print("Variable on netCDF file: " + vNam)
	208	print(" ")
	209	dim = dimsizes(vNam)
[161]	210	if (dim .EQ. 0) then
	211	print(" ")
	212	print("There are no data on file")
	213	print(" ")
	214	end if
	215
[154]	216	; ***************************************************
	217	; set up recourses
	218	; ***************************************************
	219
	220	cs_res = True
	221	cs_res@gsnYAxisIrregular2Linear = True
	222
	223	if( shape .eq. 1 ) then ; keep the shape
	224	cs_res@gsnShape = True
	225	end if
	226
	227	cs_res@gsnDraw = False
	228	cs_res@gsnFrame = False
	229	cs_res@gsnMaximize = True
	230	cs_res@gsnPaperOrientation = "portrait"
	231	cs_res@gsnPaperWidth = 8.27
	232	cs_res@gsnPaperHeight = 11.69
	233	cs_res@gsnPaperMargin = 0.79
	234	cs_res@tmXBLabelFontHeightF = .02
	235	cs_res@tmYLLabelFontHeightF = .02
	236	cs_res@tiXAxisFontHeightF = .02
	237	cs_res@tiYAxisFontHeightF = .02
	238	cs_res@tmXBMode ="Automatic"
	239	cs_res@tmYLMode ="Automatic"
	240	cs_res@lgTitleFontHeightF = .02
	241	cs_res@lgLabelFontHeightF = .02
	242	cs_res@txFontHeightF = .02
	243	cs_resP = True
	244	cs_resP@txString = f@title
	245	cs_resP@txFuncCode = "~" ; necessary for title
	246	cs_resP@txFontHeightF = .02
	247
	248	if ( mode .eq. "Fill" ) then
	249	cs_res@cnFillOn = True
	250	cs_res@gsnSpreadColors = True
	251	cs_res@cnFillMode = fill_mode
	252	cs_res@lbOrientation = "Vertical" ; vertical label bar
	253	cs_res@cnLinesOn = False
	254	cs_res@cnLineLabelsOn = False
	255	end if
	256
	257	if ( mode .eq. "Both" ) then
	258	cs_res@cnFillOn = True
	259	cs_res@gsnSpreadColors = True
	260	cs_res@cnFillMode = fill_mode
	261	cs_res@lbOrientation = "Vertical" ; vertical label bar
	262	cs_res@cnLinesOn = True
	263	cs_res@cnLineLabelsOn = True
	264	end if
	265
[157]	266	; *********************************************
	267	; set up of start_time_step and end_time_step
	268	; *********************************************
[154]	269
[161]	270	t = f->time
[157]	271	nt = dimsizes(t)
	272
	273	; *********************************************
	274	; start of time step and different types of mistakes that could be done
	275	; *********************************************
[154]	276
[157]	277	if ( .not. isvar("start_time_step") ) then
[161]	278	start_time_step = 0
[157]	279	if (parameter(13) .NE. "1") then
[161]	280	if (parameter(13) .LT. "1")
[157]	281	print(" ")
	282	print("Begin with time step 1")
	283	print(" ")
	284	exit
	285	end if
	286	if (stringtointeger(parameter(13)) .GT. nt)
	287	print(" ")
	288	print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt))
	289	print(" ")
	290	exit
	291	end if
[161]	292	start_time_step = stringtointeger(parameter(13))-1
[157]	293	end if
	294	else
	295	if (start_time_step .LE. 0)
	296	print(" ")
	297	print("Begin with time step 1")
	298	print(" ")
	299	exit
	300	end if
	301	if (start_time_step .GT. nt)
	302	print(" ")
	303	print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt))
	304	print(" ")
	305	exit
	306	end if
[161]	307	start_time_step = start_time_step - 1
[157]	308	end if
	309
	310	; ****************************************************
	311	; end of time step and different types of mistakes that could be done
	312	; ****************************************************
	313
	314	if ( .not. isvar("end_time_step") ) then
[161]	315	end_time_step = nt-1
[157]	316	if (parameter(15) .NE. "nt") then
	317	if (parameter(15) .LE. "0")
	318	print(" ")
	319	print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
	320	print(" ")
	321	exit
	322	end if
	323	if (stringtointeger(parameter(15)) .GT. nt)
	324	print(" ")
	325	print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt))
	326	print(" ")
	327	exit
	328	end if
	329	if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
	330	print(" ")
	331	print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(15))
	332	print(" ")
	333	exit
	334	end if
[161]	335	end_time_step = stringtointeger(parameter(15))-1
[157]	336	end if
	337	else
	338	if (end_time_step .LE. 0)
	339	print(" ")
	340	print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
	341	print(" ")
	342	exit
	343	end if
	344	if (end_time_step .GT. nt)
	345	print(" ")
	346	print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt))
	347	print(" ")
	348	exit
	349	end if
	350	if (end_time_step .LT. start_time_step)
	351	print(" ")
	352	print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
	353	print(" ")
	354	exit
	355	end if
[161]	356	end_time_step = end_time_step-1
[157]	357	end if
	358
[161]	359	no_time=(end_time_step-start_time_step)+1
	360
[157]	361	; ****************************************************
	362	; set up variables for vector plot if required
	363	; ****************************************************
	364
	365	if (vector .EQ. 1) then
	366	if ( .not. isvar("plotvec") ) then
	367	plotvec = parameter(69)
	368	end if
	369	if ( .not. isvar("vec1") ) then
	370	vec1 = parameter(65)
	371	if (parameter(65) .EQ. "vec1") then
	372	print(" ")
	373	print("Please indicate Vector 1 ('vec1') for Vector-Plot")
	374	print(" ")
	375	exit
	376	end if
	377	end if
	378	if ( .not. isvar("vec2") ) then
	379	vec2 = parameter(67)
	380	if (parameter(67) .EQ. "vec2") then
	381	print(" ")
	382	print("Please indicate Vector 2 ('vec2') for Vector-Plot")
	383	print(" ")
	384	exit
	385	end if
	386	end if
	387	end if
	388
	389	check_vec1 = False
	390	check_vec2 = False
	391	check_vecp = False
	392
	393	; ****************************************************
	394	; get data for plots
	395	; ****************************************************
	396
[161]	397	do varn=dim-1,0,1
[157]	398
[161]	399	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
[154]	400	check = False
[161]	401	else
	402	check = True
	403	end if
[154]	404	if ( isvar("var") ) then
[161]	405	check = isStrSubset( var,","+vNam(varn)+",")
[154]	406	end if
	407	if (parameter(21) .NE. "variables") then
	408	var = parameter(21)
[161]	409	check = isStrSubset( var,","+vNam(varn)+"," )
[157]	410	end if
[154]	411
	412	if(check) then
[161]	413	print(vNam(varn))
[157]	414	data_all = f->$vNam(varn)$
[154]	415
	416	; ****************************************************
	417	; set up ranges of x-, y- and z-coordinates
	418	; ****************************************************
	419
	420	xdim = dimsizes(data_all(0,0,0,:)) - 1
	421	ydim = dimsizes(data_all(0,0,:,0)) - 1
[157]	422	zdim = dimsizes(data_all(0,:,0,0)) - 1
	423
[154]	424	if ( .not. isvar("xe")) then ; output x-coordinate range end (in index)
	425	xe = xdim
	426	if (parameter(53) .NE. "xdim") then
	427	if (stringtointeger(parameter(53)) .GT. xdim) then
	428	print(" ")
	429	print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
	430	print(" ")
	431	exit
	432	end if
	433	if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
	434	print(" ")
	435	print("range end for x-coordinate = "+parameter(53)+" is too small")
	436	print(" ")
	437	exit
	438	end if
	439	xe = stringtointeger(parameter(53))
	440	end if
	441	else
[157]	442	if (xe .GT. xdim) then
	443	print(" ")
	444	print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
	445	print(" ")
	446	exit
	447	end if
	448	if (xe .LT. 0 .OR. xe .LT. xs) then
	449	print(" ")
	450	print("range end for x-coordinate = "+xe+" is too small")
	451	print(" ")
	452	exit
	453	end if
[154]	454	end if
[157]	455
[154]	456	if ( .not. isvar("ye")) then ; output y-coordinate range end (in index)
	457	ye = ydim
	458	if (parameter(57) .NE. "ydim") then
	459	if (stringtointeger(parameter(57)) .GT. ydim) then
	460	print(" ")
	461	print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
	462	print(" ")
	463	exit
	464	end if
	465	if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
	466	print(" ")
	467	print("range end for y-coordinate = "+parameter(57)+" is too small")
	468	print(" ")
	469	exit
	470	end if
	471	ye = stringtointeger(parameter(57))
	472	end if
	473	else
	474	if (ye .GT. ydim) then
	475	print(" ")
	476	print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
	477	print(" ")
[157]	478	exit
[154]	479	end if
	480	if (ye .LT. 0 .OR. ye .LT. ys) then
	481	print(" ")
	482	print("range end for y-coordinate = "+ye+" is too small")
	483	print(" ")
	484	exit
	485	end if
	486	end if
[157]	487
[154]	488	if ( .not. isvar("ze")) then ; output z-coordinate range end (in index)
	489	ze = zdim
	490	if (parameter(61) .NE. "zdim") then
	491	if (stringtointeger(parameter(61)) .GT. zdim) then
	492	print(" ")
	493	print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
	494	print(" ")
	495	exit
	496	end if
	497	if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
	498	print(" ")
	499	print("range end for z-coordinate = "+parameter(61)+" is too small")
	500	print(" ")
	501	exit
	502	end if
	503	ze = stringtointeger(parameter(61))
	504	end if
	505	else
	506	if (ze .GT. zdim) then
	507	print(" ")
	508	print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
	509	print(" ")
	510	exit
	511	end if
	512	if (ze .LT. 0 .OR. ze .LT. zs) then
	513	print(" ")
	514	print("range end for z-coordinate = "+ze+" is too small")
	515	print(" ")
	516	exit
	517	end if
[157]	518	end if
	519	delete(data_all)
	520	end if
	521	end do
	522
[161]	523	if (xyc .EQ. 1) then
	524	no_layer = (ze-zs)+1
	525	end if
	526	if (xzc .EQ. 1) then
	527	no_layer = (ye-ys)+1
	528	end if
	529	if (yzc .EQ. 1) then
	530	no_layer = (xe-xs)+1
	531	end if
	532
[157]	533	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
	534
	535	; ****************************************************
	536	; define inner and outer loops depending on "sort"
	537	; ****************************************************
	538
	539	if ( xyc .eq. 1 ) then
	540	lays = zs
	541	laye = ze
	542	end if
	543	if ( xzc .eq. 1 ) then
	544	lays = ys
	545	laye = ye
	546	end if
	547	if ( yzc .eq. 1) then
	548	lays = xs
	549	laye = xe
	550	end if
	551
	552	if ( sort .eq. "time" ) then
	553	lis = start_time_step
	554	lie = end_time_step
	555	los = lays
	556	loe = laye
	557	else
	558	lis = lays
	559	lie = laye
	560	los = start_time_step
	561	loe = end_time_step
	562	end if
	563
	564	check = True
	565	v1=0
	566	v2=0
[161]	567	no_var=0
	568	n=0
[157]	569
[161]	570	do varn=dim-1,0,1
[157]	571
	572	data_all = f->$vNam(varn)$
	573
	574	if (vector .EQ. 1) then
[161]	575	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
	576	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
[157]	577	end if
	578
[161]	579	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]	580	check = False
[161]	581	else
	582	check = True
	583	end if
[157]	584	if ( isvar("var") ) then
[161]	585	check = isStrSubset( var,","+vNam(varn)+",")
	586	end if
[157]	587	if (parameter(21) .NE. "variables") then
	588	var = parameter(21)
[161]	589	check = isStrSubset( var,","+vNam(varn)+"," )
[157]	590	end if
	591
	592	if(check) then
	593
[161]	594	no_var=no_var+1
	595
[157]	596	if (vector .EQ. 1) then
[161]	597	if (","+vNam(varn)+"," .EQ. vec1) then
[157]	598	v1=v1+1
	599	end if
[161]	600	if (","+vNam(varn)+"," .EQ. vec2) then
[157]	601	v2=v2+1
	602	end if
	603	end if
	604
	605	data(varn,:,:,:,:)=data_all(0:nt-1,zs:ze,ys:ye,xs:xe)
	606	if (check_vec1) then
	607	vect1=data_all
[154]	608	end if
[157]	609	if (check_vec2) then
	610	vect2=data_all
[154]	611	end if
	612
[157]	613	data!0 = "var"
	614	data!1 = "t"
	615	data!2 = "z"
	616	data!3 = "y"
	617	data!4 = "x"
[154]	618
[157]	619	end if
[154]	620
[157]	621	delete(data_all)
[154]	622
[161]	623	end do
	624
	625	plot_panel=new((/no_timeno_layerno_var/),graphic)
[154]	626
[157]	627	if (vector .EQ. 1) then
	628	if (v1 .EQ. 0)then
	629	print(" ")
	630	print("Varible for vector-plot ('vec1') must be one of the varibles for plot ('var')")
	631	print(" ")
	632	exit
	633	end if
[154]	634
[157]	635	if (v2 .EQ. 0)then
	636	print(" ")
	637	print("Varible for vector-plot ('vec2') must be one of the varibles for plot ('var')")
	638	print(" ")
	639	exit
	640	end if
	641	end if
[154]	642
[157]	643	; ***************************************************
[161]	644	; open workstation(s)
	645	; ***************************************************
	646
	647	wks_ps = gsn_open_wks(format_out,file_out)
	648	gsn_define_colormap(wks_ps,"rainbow+white")
	649
	650	plot=new((/no_timeno_layerno_var/),graphic)
	651
	652	page = 0
	653	n=0
	654
	655	print("plots sorted by '"+sort+"'")
	656	print(" ")
	657
	658	; ***************************************************
[157]	659	; create plots
	660	; ***************************************************
	661
	662	check = True
	663
[161]	664	do varn=dim-1,0,1
[157]	665
	666	if (vector .EQ. 1) then
[161]	667	check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
[157]	668	end if
	669
[161]	670	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
	671	check = False
	672	else
	673	check = True
	674	end if
[157]	675	if ( isvar("var") ) then
[161]	676	check = isStrSubset( var,","+vNam(varn)+",")
[157]	677	end if
	678	if (parameter(21) .NE. "variables") then
	679	var = parameter(21)
[161]	680	check = isStrSubset( var,","+vNam(varn)+"," )
[157]	681	end if
	682
	683	if(check) then
[161]	684
[157]	685	; ****************************************************
	686	; loops over time and layer
	687	; ****************************************************
	688
[161]	689	do lo = los, loe
	690	do li = lis, lie
	691
[157]	692	; ****************************************************
[154]	693	; xy cross section
[157]	694	; ****************************************************
[154]	695
[157]	696	if ( xyc .eq. 1 ) then
	697
[161]	698	cs_res@tiXAxisString = "x [m]"
	699	cs_res@tiYAxisString = "y [m]"
[157]	700	cs_res@gsnLeftString = "Plot of "+vNam(varn)
[161]	701
[154]	702	if ( sort .eq. "time" ) then
[157]	703	if ( data&z(lo) .eq. -1 ) then
[154]	704	level = "-average"
	705	else
[157]	706	level = "=" + data&z(lo) + "m"
[154]	707	end if
[161]	708	cs_res@gsnCenterString = "time=" + data&t(li) +"s z"+level
	709	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
	710	if (vector .EQ. 1) then
[157]	711	if (check_vecp)
	712	vecres = True ; vector only resources
	713	vecres@gsnDraw = False ; don't draw
[161]	714	vecres@gsnFrame = False ; don't advance frame
	715	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	716	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	717	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	718	vecres@gsnRightString = " " ; turn off right string
	719	vecres@gsnLeftString = " " ; turn off left string
	720	vecres@tiXAxisString = " "
	721	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	722	overlay(plot(n), plot_vec)
	723	else
	724	vecres = True ; vector only resources
	725	vecres@gsnDraw = False ; don't draw
[157]	726	vecres@gsnFrame = False ; don't advance frame
	727	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	728	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	729	vecres@vcRefLengthF = 0.05 ; define length of vec ref
[161]	730	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
[157]	731	vecres@gsnLeftString = " " ; turn off left string
	732	vecres@tiXAxisString = " "
[161]	733	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	734	plot(n) = plot_vec
	735	end if
	736	end if
[154]	737	end if
	738
	739	if ( sort .eq. "layer" ) then
[161]	740	if ( data&z(li) .eq. -1 ) then
	741	level = "-average"
	742	else
	743	level = "=" + data&z(li) + "m"
	744	end if
	745	cs_res@gsnCenterString = "t=" + data&t(lo) + "s z"+ level
	746	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
	747	if (vector .EQ. 1) then
[157]	748	if (check_vecp)
	749	vecres = True ; vector only resources
	750	vecres@gsnDraw = False ; don't draw
	751	vecres@gsnFrame = False ; don't advance frame
	752	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	753	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	754	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	755	vecres@gsnRightString = " " ; turn off right string
	756	vecres@gsnLeftString = " " ; turn off left string
	757	vecres@tiXAxisString = " "
[161]	758	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	759	overlay(plot(n), plot_vec)
	760	else
	761	vecres = True ; vector only resources
	762	vecres@gsnDraw = False ; don't draw
	763	vecres@gsnFrame = False ; don't advance frame
	764	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	765	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	766	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	767	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
	768	vecres@gsnLeftString = " " ; turn off left string
	769	vecres@tiXAxisString = " "
	770	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	771	plot(n) = plot_vec
	772	end if
	773	end if
[154]	774	end if
	775	end if
	776
[157]	777	; ****************************************************
[154]	778	; xz cross section
[157]	779	; ****************************************************
[154]	780
	781	if ( xzc .eq. 1 ) then
[157]	782
[154]	783	cs_res@tiXAxisString = "x [m]"
	784	cs_res@tiYAxisString = "z [m]"
[157]	785	cs_res@gsnLeftString = "Plot of "+vNam(varn)
	786
[154]	787	if ( sort .eq. "time" ) then
[157]	788	if ( data&z(lo) .eq. -1 ) then
[154]	789	level = "-average"
	790	else
[157]	791	level = "=" + data&y(lo) + "m"
[154]	792	end if
[161]	793	cs_res@gsnCenterString = "t=" + data&t(li) + "s y"+ level
	794	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
	795	if (vector .EQ. 1) then
[157]	796	if (check_vecp)
	797	vecres = True ; vector only resources
	798	vecres@gsnDraw = False ; don't draw
	799	vecres@gsnFrame = False ; don't advance frame
	800	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	801	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	802	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	803	vecres@gsnRightString = " " ; turn off right string
	804	vecres@gsnLeftString = " " ; turn off left string
	805	vecres@tiXAxisString = " "
[161]	806	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	807	overlay(plot(n), plot_vec)
	808	else
	809	vecres = True ; vector only resources
	810	vecres@gsnDraw = False ; don't draw
	811	vecres@gsnFrame = False ; don't advance frame
	812	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	813	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	814	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	815	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
	816	vecres@gsnLeftString = " " ; turn off left string
	817	vecres@tiXAxisString = " "
	818	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	819	plot(n) = plot_vec
[157]	820	end if
[161]	821	end if
[154]	822	end if
	823
	824	if ( sort .eq. "layer" ) then
[161]	825	if ( data&z(li) .eq. -1 ) then
	826	level = "-average"
	827	else
	828	level = "=" + data&y(li) + "m"
	829	end if
	830	cs_res@gsnCenterString = "t=" + data&t(lo) + "s y"+ level
	831	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
	832	if (vector .EQ. 1) then
[157]	833	if (check_vecp)
	834	vecres = True ; vector only resources
	835	vecres@gsnDraw = False ; don't draw
	836	vecres@gsnFrame = False ; don't advance frame
	837	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	838	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	839	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	840	vecres@gsnRightString = " " ; turn off right string
	841	vecres@gsnLeftString = " " ; turn off left string
	842	vecres@tiXAxisString = " "
[161]	843	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	844	overlay(plot(n), plot_vec)
	845	else
	846	vecres = True ; vector only resources
	847	vecres@gsnDraw = False ; don't draw
	848	vecres@gsnFrame = False ; don't advance frame
	849	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	850	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	851	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	852	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
	853	vecres@gsnLeftString = " " ; turn off left string
	854	vecres@tiXAxisString = " "
	855	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	856	plot(n) = plot_vec
[157]	857	end if
[161]	858	end if
[157]	859	end if
[154]	860	end if
	861
[157]	862	; ****************************************************
[154]	863	; yz cross section
[157]	864	; ****************************************************
[154]	865
	866	if ( yzc .eq. 1 ) then
[157]	867
[154]	868	cs_res@tiXAxisString = "y [m]"
	869	cs_res@tiYAxisString = "z [m]"
[157]	870	cs_res@gsnLeftString = "Plot of "+vNam(varn)
	871
[154]	872	if ( sort .eq. "time" ) then
[157]	873	if ( data&z(lo) .eq. -1 ) then
[154]	874	level = "-average"
	875	else
[157]	876	level = "=" + data&z(lo) + "m"
[154]	877	end if
[161]	878	cs_res@gsnCenterString = "t=" + data&t(li) + "s x"+ level
	879	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
	880	if (vector .EQ. 1) then
[157]	881	if (check_vecp)
	882	vecres = True ; vector only resources
	883	vecres@gsnDraw = False ; don't draw
	884	vecres@gsnFrame = False ; don't advance frame
	885	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	886	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	887	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	888	vecres@gsnRightString = " " ; turn off right string
	889	vecres@gsnLeftString = " " ; turn off left string
	890	vecres@tiXAxisString = " "
[161]	891	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	892	overlay(plot(n), plot_vec)
	893	else
	894	vecres = True ; vector only resources
	895	vecres@gsnDraw = False ; don't draw
	896	vecres@gsnFrame = False ; don't advance frame
	897	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	898	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	899	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	900	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
	901	vecres@gsnLeftString = " " ; turn off left string
	902	vecres@tiXAxisString = " "
	903	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
	904	plot(n) = plot_vec
[157]	905	end if
[161]	906	end if
[154]	907	end if
	908
	909	if ( sort .eq. "layer" ) then
[161]	910	if ( data&z(li) .eq. -1 ) then
	911	level = "-average"
	912	else
	913	level = "=" + data&z(li) + "m"
	914	end if
	915	cs_res@gsnCenterString = "t=" + data&t(lo) + "s x"+ level
	916	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
	917	if (vector .EQ. 1) then
[157]	918	if (check_vecp)
	919	vecres = True ; vector only resources
	920	vecres@gsnDraw = False ; don't draw
	921	vecres@gsnFrame = False ; don't advance frame
	922	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	923	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	924	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	925	vecres@gsnRightString = " " ; turn off right string
	926	vecres@gsnLeftString = " " ; turn off left string
	927	vecres@tiXAxisString = " "
[161]	928	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	929	overlay(plot(n), plot_vec)
	930	else
	931	vecres = True ; vector only resources
	932	vecres@gsnDraw = False ; don't draw
	933	vecres@gsnFrame = False ; don't advance frame
	934	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
	935	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
	936	vecres@vcRefLengthF = 0.05 ; define length of vec ref
	937	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
	938	vecres@gsnLeftString = " " ; turn off left string
	939	vecres@tiXAxisString = " "
	940	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
	941	plot(n) = plot_vec
[157]	942	end if
[161]	943	end if
[154]	944	end if
[161]	945	end if
	946	n=n+1
[154]	947	end do
	948	end do
[157]	949	end if
[161]	950	end do
[154]	951
[161]	952	; ***************************************************
	953	; merge plots onto one page
	954	; ***************************************************
	955
	956	if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
	957	gsn_panel(wks_ps,plot(0:(no_timeno_layerno_var)-1),(/no_var,no_layer*no_time/),cs_resP)
	958	print(" ")
	959	print("Outputs to .eps or .epsi have only one frame")
	960	print(" ")
	961	else
	962	do np = 0,no_layerno_timeno_var-1,no_lines*no_columns
	963	if ( np + no_linesno_columns .gt. (no_layerno_time*no_var)-1) then
	964	gsn_panel(wks_ps, plot(np:(no_layerno_timeno_var)-1),(/no_lines,no_columns/),cs_resP)
	965	else
	966	gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
	967	end if
	968	end do
	969	end if
	970
[154]	971	print(" ")
	972	print("Output to: " + file_out +"."+ format_out)
	973	print(" ")
	974
	975	end

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

Download in other formats:

Original Format