source: palm/trunk/SCRIPTS/NCL/cross_sections.ncl @ 161

Last change on this file since 161 was 161, checked in by letzel, 14 years ago
  • NCL scripts in trunk/SCRIPTS/NCL updated for vector plots
File size: 38.7 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
5load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
6
7begin
8
9   ; ***************************************************
10   ; read parameter_list
11   ; ***************************************************
12
13   if (isfilepresent("~/.ncl_preferences")) then
14      parameter = asciiread("~/.ncl_preferences",75,"string")
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
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
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)
210   if (dim .EQ. 0) then
211      print(" ")
212      print("There are no data on file")
213      print(" ")
214   end if
215
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
266   ; *********************************************
267   ; set up of start_time_step and end_time_step
268   ; *********************************************
269
270   t    = f->time
271   nt = dimsizes(t)
272     
273   ; *********************************************     
274   ; start of time step and different types of mistakes that could be done
275   ; *********************************************
276
277   if ( .not. isvar("start_time_step") ) then           
278      start_time_step = 0
279      if (parameter(13) .NE. "1") then
280         if (parameter(13) .LT. "1")
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
292         start_time_step = stringtointeger(parameter(13))-1 
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
307      start_time_step = start_time_step - 1
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             
315      end_time_step = nt-1
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
335         end_time_step = stringtointeger(parameter(15))-1 
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
356      end_time_step = end_time_step-1
357   end if
358
359   no_time=(end_time_step-start_time_step)+1
360
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
397   do varn=dim-1,0,1 
398
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
400         check = False
401      else
402         check = True
403      end if   
404      if (  isvar("var") ) then 
405         check = isStrSubset( var,","+vNam(varn)+",")
406      end if
407      if (parameter(21) .NE. "variables") then
408         var = parameter(21)
409         check = isStrSubset( var,","+vNam(varn)+"," )
410      end if   
411
412      if(check) then
413         print(vNam(varn))
414         data_all = f->$vNam(varn)$
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       
422         zdim = dimsizes(data_all(0,:,0,0)) - 1 
423       
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
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                 
454         end if
455
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(" ")
478               exit
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
487
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
518         end if   
519         delete(data_all)
520      end if 
521   end do
522
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
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
567   no_var=0
568   n=0
569
570   do varn=dim-1,0,1 
571     
572      data_all = f->$vNam(varn)$
573
574      if (vector .EQ. 1) then   
575         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
576         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
577      end if
578           
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
580         check = False
581      else
582         check = True
583      end if 
584      if (  isvar("var") ) then 
585         check = isStrSubset( var,","+vNam(varn)+",")
586      end if 
587      if (parameter(21) .NE. "variables") then
588         var = parameter(21)
589         check = isStrSubset( var,","+vNam(varn)+"," )
590      end if 
591
592      if(check) then
593
594         no_var=no_var+1
595
596         if (vector .EQ. 1) then
597            if (","+vNam(varn)+"," .EQ. vec1) then
598               v1=v1+1
599            end if
600            if (","+vNam(varn)+"," .EQ. vec2) then
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
608         end if
609         if (check_vec2) then
610            vect2=data_all
611         end if
612
613         data!0 = "var"
614         data!1 = "t"
615         data!2 = "z"
616         data!3 = "y"
617         data!4 = "x"   
618
619      end if
620
621      delete(data_all)
622
623   end do
624 
625   plot_panel=new((/no_time*no_layer*no_var/),graphic)
626
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
634
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
642
643   ; ***************************************************
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_time*no_layer*no_var/),graphic)
651
652   page = 0
653   n=0
654
655   print("plots sorted by '"+sort+"'")
656   print(" ")
657
658   ; ***************************************************
659   ; create plots
660   ; ***************************************************
661
662   check = True
663
664   do varn=dim-1,0,1
665
666      if (vector .EQ. 1) then   
667         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
668      end if
669
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   
675      if (  isvar("var") ) then 
676         check = isStrSubset( var,","+vNam(varn)+",")
677      end if
678      if (parameter(21) .NE. "variables") then
679         var = parameter(21)
680         check = isStrSubset( var,","+vNam(varn)+"," )
681      end if
682   
683      if(check) then
684   
685         ; ****************************************************
686         ; loops over time and layer
687         ; ****************************************************
688
689         do lo = los, loe                                       
690            do li = lis, lie                           
691                 
692               ; ****************************************************
693               ; xy cross section
694               ; ****************************************************
695
696               if ( xyc .eq. 1 ) then
697               
698                  cs_res@tiXAxisString = "x [m]"
699                  cs_res@tiYAxisString = "y [m]"
700                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
701     
702                  if ( sort .eq. "time" ) then
703                     if ( data&z(lo) .eq. -1 ) then
704                        level = "-average"
705                     else
706                        level = "=" + data&z(lo) + "m"
707                     end if
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
711                        if (check_vecp)
712                           vecres                  = True            ; vector only resources
713                           vecres@gsnDraw          = False           ; don't draw
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
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
730                           vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
731                           vecres@gsnLeftString    = " "             ; turn off left string
732                           vecres@tiXAxisString    = " "   
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                         
737                  end if
738
739                  if ( sort .eq. "layer" ) then
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
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    = " "   
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
774                  end if
775               end if
776
777               ; ****************************************************
778               ; xz cross section
779               ; ****************************************************
780
781               if ( xzc .eq. 1 ) then
782       
783                  cs_res@tiXAxisString   = "x [m]"
784                  cs_res@tiYAxisString   = "z [m]"
785                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
786
787                  if ( sort .eq. "time" ) then
788                     if ( data&z(lo) .eq. -1 ) then
789                        level = "-average"
790                     else
791                        level = "=" + data&y(lo) + "m"
792                     end if
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
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    = " "   
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
820                        end if
821                     end if
822                  end if
823
824                  if ( sort .eq. "layer" ) then
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
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    = " "   
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
857                        end if
858                     end if
859                  end if                 
860               end if
861
862               ; ****************************************************
863               ; yz cross section
864               ; ****************************************************
865
866               if ( yzc .eq. 1 ) then
867               
868                  cs_res@tiXAxisString   = "y [m]"
869                  cs_res@tiYAxisString   = "z [m]"
870                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
871
872                  if ( sort .eq. "time" ) then
873                     if ( data&z(lo) .eq. -1 ) then
874                        level = "-average"
875                     else
876                        level = "=" + data&z(lo) + "m"
877                     end if
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
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    = " "   
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
905                        end if
906                     end if
907                  end if
908
909                  if ( sort .eq. "layer" ) then
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
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    = " "   
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
942                        end if
943                     end if
944                  end if
945               end if         
946               n=n+1   
947            end do     
948         end do 
949      end if
950   end do
951
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_time*no_layer*no_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_layer*no_time*no_var-1,no_lines*no_columns   
963         if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then
964            gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_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
971   print(" ")
972   print("Output to: " + file_out +"."+ format_out)
973   print(" ")   
974 
975end
Note: See TracBrowser for help on using the repository browser.