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

Last change on this file since 160 was 157, checked in by letzel, 17 years ago
  • NCL scripts in trunk/SCRIPTS/NCL updated for vector plots
File size: 34.1 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",73,"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 
211   ; ***************************************************
212   ; set up recourses
213   ; ***************************************************
214
215   cs_res                          = True
216   cs_res@gsnYAxisIrregular2Linear = True 
217
218   if( shape .eq. 1 ) then                              ; keep the shape
219      cs_res@gsnShape = True
220   end if
221 
222   cs_res@gsnDraw                 = False
223   cs_res@gsnFrame                = False
224   cs_res@gsnMaximize             = True
225   cs_res@gsnPaperOrientation     = "portrait"
226   cs_res@gsnPaperWidth           = 8.27
227   cs_res@gsnPaperHeight          = 11.69
228   cs_res@gsnPaperMargin          = 0.79
229   cs_res@tmXBLabelFontHeightF   = .02
230   cs_res@tmYLLabelFontHeightF   = .02
231   cs_res@tiXAxisFontHeightF     = .02
232   cs_res@tiYAxisFontHeightF     = .02
233   cs_res@tmXBMode                ="Automatic"
234   cs_res@tmYLMode                ="Automatic"
235   cs_res@lgTitleFontHeightF     = .02
236   cs_res@lgLabelFontHeightF     = .02
237   cs_res@txFontHeightF          = .02
238   cs_resP = True
239   cs_resP@txString               = f@title
240   cs_resP@txFuncCode             = "~"                 ; necessary for title
241   cs_resP@txFontHeightF          = .02
242 
243   if ( mode .eq. "Fill" ) then
244      cs_res@cnFillOn             = True
245      cs_res@gsnSpreadColors      = True
246      cs_res@cnFillMode           = fill_mode
247      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
248      cs_res@cnLinesOn            = False       
249      cs_res@cnLineLabelsOn       = False
250   end if
251
252   if ( mode .eq. "Both" ) then
253      cs_res@cnFillOn             = True
254      cs_res@gsnSpreadColors      = True
255      cs_res@cnFillMode           = fill_mode
256      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
257      cs_res@cnLinesOn            = True
258      cs_res@cnLineLabelsOn       = True
259   end if
260
261   ; ***************************************************
262   ; open workstation(s)
263   ; ***************************************************
264
265   wks_ps  = gsn_open_wks(format_out,file_out)
266   gsn_define_colormap(wks_ps,"rainbow+white")
267
268   plot=new((/no_columns*no_lines/),graphic)
269
270   page = 0
271
272   nc1 = no_columns
273   nl1 = no_lines
274
275   print("plots sorted by '"+sort+"'")
276   print(" ")
277
278   ; *********************************************
279   ; set up of start_time_step and end_time_step
280   ; *********************************************
281
282   t = f->time
283   nt = dimsizes(t)
284     
285   ; *********************************************     
286   ; start of time step and different types of mistakes that could be done
287   ; *********************************************
288
289   if ( .not. isvar("start_time_step") ) then           
290      start_time_step = 1
291      if (parameter(13) .NE. "1") then
292         if (parameter(13) .LE. "0")
293            print(" ")
294            print("Begin with time step 1")
295            print(" ")
296            exit
297         end if
298         if (stringtointeger(parameter(13)) .GT. nt)
299            print(" ")
300            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt))
301            print(" ")
302            exit
303         end if
304         start_time_step = stringtointeger(parameter(13)) 
305      end if
306   else
307      if (start_time_step .LE. 0)
308         print(" ")
309         print("Begin with time step 1")
310         print(" ")
311         exit
312      end if
313      if (start_time_step .GT. nt)
314         print(" ")
315         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt))
316         print(" ")
317         exit
318      end if
319   end if
320
321   ; ****************************************************
322   ; end of time step and different types of mistakes that could be done
323   ; ****************************************************
324
325   if ( .not. isvar("end_time_step") ) then             
326      end_time_step = nt
327      if (parameter(15) .NE. "nt") then
328         if (parameter(15) .LE. "0")
329            print(" ")
330            print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
331            print(" ")
332            exit
333         end if
334         if (stringtointeger(parameter(15)) .GT. nt)
335            print(" ")
336            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt))
337            print(" ")
338            exit
339         end if
340         if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
341            print(" ")
342            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(15))
343            print(" ")
344            exit
345         end if
346         end_time_step = stringtointeger(parameter(15)) 
347      end if   
348   else
349      if (end_time_step .LE. 0)
350         print(" ")
351         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
352         print(" ")
353         exit
354      end if
355      if (end_time_step .GT. nt)
356         print(" ")
357         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt))
358         print(" ")
359         exit
360      end if
361      if (end_time_step .LT. start_time_step)
362         print(" ")
363         print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
364         print(" ")
365         exit
366      end if
367   end if
368
369   ; ****************************************************
370   ; set up variables for vector plot if required
371   ; ****************************************************   
372
373   if (vector .EQ. 1) then
374      if ( .not. isvar("plotvec") ) then
375         plotvec = parameter(69)
376         if (parameter(69) .EQ. "plotvec") then
377            print(" ")
378            print("Please indicate the variable where the vector plot shall overlay ('plotvec')")
379            print(" ")
380            exit
381         end if
382      end if
383      if ( .not. isvar("vec1") ) then
384         vec1 = parameter(65)
385         if (parameter(65) .EQ. "vec1") then
386            print(" ")
387            print("Please indicate Vector 1 ('vec1') for Vector-Plot")
388            print(" ")
389            exit
390         end if
391      end if
392      if ( .not. isvar("vec2") ) then
393         vec2 = parameter(67)
394         if (parameter(67) .EQ. "vec2") then
395            print(" ")
396            print("Please indicate Vector 2 ('vec2') for Vector-Plot")
397            print(" ")
398            exit
399         end if
400      end if             
401   end if
402
403   check_vec1 = False
404   check_vec2 = False
405   check_vecp = False
406
407   ; ****************************************************
408   ; get data for plots
409   ; ****************************************************
410
411   do varn=0,dim-1 
412
413      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
414         check = False
415      end if
416      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
417         check = False
418      end if     
419      if (  isvar("var") ) then 
420         check = isStrSubset( var,vNam(varn)+",")
421      end if
422      if (parameter(21) .NE. "variables") then
423         var = parameter(21)
424         check = isStrSubset( var,vNam(varn)+"," )
425      end if   
426
427      if(check) then
428       
429         data_all = f->$vNam(varn)$
430
431         ; ****************************************************
432         ; set up ranges of x-, y- and z-coordinates
433         ; ****************************************************         
434                   
435         xdim = dimsizes(data_all(0,0,0,:)) - 1   
436         ydim = dimsizes(data_all(0,0,:,0)) - 1       
437         zdim = dimsizes(data_all(0,:,0,0)) - 1 
438       
439         if ( .not. isvar("xe")) then     ; output x-coordinate range end (in index)
440            xe = xdim
441            if (parameter(53) .NE. "xdim") then
442               if (stringtointeger(parameter(53)) .GT. xdim) then
443                  print(" ")
444                  print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
445                  print(" ")
446                  exit
447               end if
448               if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
449                  print(" ")
450                  print("range end for x-coordinate = "+parameter(53)+" is too small")
451                  print(" ")
452                  exit
453               end if
454               xe = stringtointeger(parameter(53))
455            end if
456         else
457            if (xe .GT. xdim) then
458               print(" ")
459               print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
460               print(" ")
461               exit
462            end if
463            if (xe .LT. 0 .OR. xe .LT. xs) then
464               print(" ")
465               print("range end for x-coordinate = "+xe+" is too small")
466               print(" ")
467               exit
468            end if                 
469         end if
470
471         if ( .not. isvar("ye")) then     ; output y-coordinate range end (in index)
472            ye = ydim
473            if (parameter(57) .NE. "ydim") then
474               if (stringtointeger(parameter(57)) .GT. ydim) then
475                  print(" ")
476                  print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
477                  print(" ")
478                  exit
479               end if
480               if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
481                  print(" ")
482                  print("range end for y-coordinate = "+parameter(57)+" is too small")
483                  print(" ")
484                  exit
485               end if
486               ye = stringtointeger(parameter(57))
487            end if
488         else
489            if (ye .GT. ydim) then
490               print(" ")
491               print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
492               print(" ")
493               exit
494            end if
495            if (ye .LT. 0 .OR. ye .LT. ys) then
496               print(" ")
497               print("range end for y-coordinate = "+ye+" is too small")
498               print(" ")
499               exit
500            end if
501         end if
502
503         if ( .not. isvar("ze")) then     ; output z-coordinate range end (in index)
504            ze = zdim
505            if (parameter(61) .NE. "zdim") then
506               if (stringtointeger(parameter(61)) .GT. zdim) then
507                  print(" ")
508                  print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
509                  print(" ")
510                  exit
511               end if
512               if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
513                  print(" ")
514                  print("range end for z-coordinate = "+parameter(61)+" is too small")
515                  print(" ")
516                  exit
517               end if
518               ze = stringtointeger(parameter(61))
519            end if
520         else
521            if (ze .GT. zdim) then
522               print(" ")
523               print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
524               print(" ")
525               exit
526            end if
527            if (ze .LT. 0 .OR. ze .LT. zs) then
528               print(" ")
529               print("range end for z-coordinate = "+ze+" is too small")
530               print(" ")
531               exit
532            end if
533         end if   
534         delete(data_all)
535      end if 
536   end do
537
538   data  = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
539
540   ; ****************************************************
541   ; define inner and outer loops depending on "sort"
542   ; ****************************************************       
543   
544   if ( xyc .eq. 1 ) then
545      lays = zs
546      laye = ze
547   end if
548   if ( xzc .eq. 1 ) then
549      lays = ys
550      laye = ye
551   end if
552   if ( yzc .eq. 1) then
553      lays = xs
554      laye = xe
555   end if
556
557   if ( sort .eq. "time" ) then
558      lis = start_time_step
559      lie = end_time_step
560      los = lays
561      loe = laye
562   else
563      lis = lays
564      lie = laye
565      los = start_time_step
566      loe = end_time_step
567   end if
568
569   check = True
570   v1=0
571   v2=0
572
573   do varn=0,dim-1 
574     
575      data_all = f->$vNam(varn)$
576
577      if (vector .EQ. 1) then   
578         check_vec1 = isStrSubset( vec1,vNam(varn)+",")
579         check_vec2 = isStrSubset( vec2,vNam(varn)+",")
580      end if
581           
582      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
583         check = False
584      end if
585      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
586         check = False
587      end if     
588      if (  isvar("var") ) then 
589         check = isStrSubset( var,vNam(varn)+",")
590      end if
591 
592      if (parameter(21) .NE. "variables") then
593         var = parameter(21)
594         check = isStrSubset( var,vNam(varn)+"," )
595      end if 
596
597      if(check) then
598
599         if (vector .EQ. 1) then
600            if (vNam(varn)+"," .EQ. vec1) then
601               v1=v1+1
602            end if
603            if (vNam(varn)+"," .EQ. vec2) then
604               v2=v2+1
605            end if
606         end if
607
608         data(varn,:,:,:,:)=data_all(0:nt-1,zs:ze,ys:ye,xs:xe)
609         if (check_vec1) then
610            print("in vec1")
611            vect1=data_all
612         end if
613         if (check_vec2) then
614             print("in vec2")
615            vect2=data_all
616         end if
617
618         data!0 = "var"
619         data!1 = "t"
620         data!2 = "z"
621         data!3 = "y"
622         data!4 = "x"   
623
624      end if
625
626      delete(data_all)
627
628   end do 
629
630   if (vector .EQ. 1) then
631      if (v1 .EQ. 0)then
632         print(" ")
633         print("Varible for vector-plot ('vec1') must be one of the varibles for plot ('var')")
634         print(" ")
635         exit
636      end if
637
638      if (v2 .EQ. 0)then
639         print(" ")
640         print("Varible for vector-plot ('vec2') must be one of the varibles for plot ('var')")
641         print(" ")
642         exit
643      end if
644   end if
645
646   ; ***************************************************
647   ; create plots
648   ; ***************************************************
649
650   check = True
651
652   do varn=0,dim-1
653
654      if (vector .EQ. 1) then   
655         check_vecp = isStrSubset( plotvec,vNam(varn)+",")
656      end if
657
658      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
659         check = False
660      end if
661      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
662         check = False
663      end if     
664      if (  isvar("var") ) then 
665         check = isStrSubset( var,vNam(varn)+",")
666      end if
667      if (parameter(21) .NE. "variables") then
668         var = parameter(21)
669         check = isStrSubset( var,vNam(varn)+"," )
670      end if
671   
672      if(check) then
673     
674         ; ****************************************************
675         ; loops over time and layer
676         ; ****************************************************
677
678         do lo = los, loe                                       ; lo: loop outer
679            do li = lis, lie, no_lines*no_columns               ; li: loop inner
680               
681               if ( sort .eq. "time" ) then
682                  ta = li
683                  tb = min( (/li + no_lines*no_columns - 1, end_time_step/) )
684                  if ( xyc .eq. 1 ) then
685                     xa = xs
686                     xb = xe
687                     ya = ys
688                     yb = ye
689                     za = lo
690                     zb = lo
691                  end if
692                  if ( xzc .eq. 1 ) then
693                     xa = xs
694                     xb = xe
695                     ya = lo
696                     yb = lo
697                     za = zs
698                     zb = ze
699                  end if
700                  if ( yzc .eq. 1 ) then
701                     xa = lo
702                     xb = lo
703                     ya = ys
704                     yb = ye
705                     za = zs
706                     zb = ze
707                  end if
708               end if
709
710               if ( sort .eq. "layer" ) then
711                  ta = lo
712                  tb = lo
713                  if ( xyc .eq. 1 ) then
714                     xa = xs
715                     xb = xe
716                     ya = ys
717                     yb = ye
718                     za = li
719                     zb = min( (/li + no_lines*no_columns - 1, ze/) ) 
720                  end if
721                  if ( xzc .eq. 1 ) then
722                     xa = xs
723                     xb = xe
724                     ya = li
725                     yb = min( (/li + no_lines*no_columns - 1, ye/) )
726                     za = zs
727                     zb = ze
728                  end if
729                  if ( yzc .eq. 1 ) then
730                     xa = li
731                     xb = min( (/li + no_lines*no_columns - 1, xe/) )
732                     ya = ys
733                     yb = ye
734                     za = zs
735                     zb = ze
736                  end if
737               end if   
738
739               ; ****************************************************
740               ; xy cross section
741               ; ****************************************************
742
743               if ( xyc .eq. 1 ) then
744               
745                  cs_res@tiXAxisString   = "x [m]"
746                  cs_res@tiYAxisString   = "y [m]"
747                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
748                 
749                  if ( sort .eq. "time" ) then
750                     if ( data&z(lo) .eq. -1 ) then
751                        level = "-average"
752                     else
753                        level = "=" + data&z(lo) + "m"
754                     end if
755                     do n = 0,tb-ta
756                        cs_res@gsnCenterString = "time=" + data&t(li+n-1) +"s  z"+level
757                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,lo,:,:),cs_res) 
758                        if (check_vecp)
759                           vecres                  = True            ; vector only resources
760                           vecres@gsnDraw          = False           ; don't draw
761                           vecres@gsnFrame         = False           ; don't advance frame
762                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
763                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
764                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
765                           vecres@gsnRightString   = " "             ; turn off right string
766                           vecres@gsnLeftString    = " "             ; turn off left string
767                           vecres@tiXAxisString    = " "   
768                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
769                           overlay(plot1, plot_vec)
770                        end if                     
771                        plot(n) = plot1                         
772                     end do
773                  end if
774
775                  if ( sort .eq. "layer" ) then
776                     do n = 0, zb-za
777                        if ( data&z(li+n) .eq. -1 ) then
778                           level = "-average"
779                        else
780                           level = "=" + data&z(li+n) + "m"
781                        end if
782                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  z"+ level
783                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,li+n,:,:),cs_res)
784                        if (check_vecp)
785                           vecres                  = True            ; vector only resources
786                           vecres@gsnDraw          = False           ; don't draw
787                           vecres@gsnFrame         = False           ; don't advance frame
788                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
789                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
790                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
791                           vecres@gsnRightString   = " "             ; turn off right string
792                           vecres@gsnLeftString    = " "             ; turn off left string
793                           vecres@tiXAxisString    = " "   
794                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo-1,li+n,:,:),vect2(lo-1,li+n,:,:),vecres)
795                           overlay(plot1, plot_vec)
796                        end if
797                        plot(n) = plot1
798                     end do
799                  end if
800
801               end if
802
803               ; ****************************************************
804               ; xz cross section
805               ; ****************************************************
806
807               if ( xzc .eq. 1 ) then
808       
809                  cs_res@tiXAxisString   = "x [m]"
810                  cs_res@tiYAxisString   = "z [m]"
811                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
812
813                  if ( sort .eq. "time" ) then
814                     if ( data&z(lo) .eq. -1 ) then
815                        level = "-average"
816                     else
817                        level = "=" + data&y(lo) + "m"
818                     end if
819                     do n = 0, tb-ta
820                        cs_res@gsnCenterString = "t=" + data&t(li+n-1) + "s  y"+ level
821                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,:,lo,:),cs_res)
822                        if (check_vecp)
823                           vecres                  = True            ; vector only resources
824                           vecres@gsnDraw          = False           ; don't draw
825                           vecres@gsnFrame         = False           ; don't advance frame
826                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
827                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
828                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
829                           vecres@gsnRightString   = " "             ; turn off right string
830                           vecres@gsnLeftString    = " "             ; turn off left string
831                           vecres@tiXAxisString    = " "   
832                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
833                           overlay(plot1, plot_vec)
834                        end if
835                        plot(n) = plot1 
836                     end do
837                  end if
838
839                  if ( sort .eq. "layer" ) then
840                     do n = 0, yb-ya
841                        if ( data&z(li+n) .eq. -1 ) then
842                           level = "-average"
843                        else
844                           level = "=" + data&y(li+n) + "m"
845                        end if
846                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  y"+ level
847                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,:,li+n,:),cs_res)
848                        if (check_vecp)
849                           vecres                  = True            ; vector only resources
850                           vecres@gsnDraw          = False           ; don't draw
851                           vecres@gsnFrame         = False           ; don't advance frame
852                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
853                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
854                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
855                           vecres@gsnRightString   = " "             ; turn off right string
856                           vecres@gsnLeftString    = " "             ; turn off left string
857                           vecres@tiXAxisString    = " "   
858                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
859                           overlay(plot1, plot_vec)
860                        end if
861                        plot(n) = plot1
862                     end do
863                  end if                 
864                 
865               end if
866
867               ; ****************************************************
868               ; yz cross section
869               ; ****************************************************
870
871               if ( yzc .eq. 1 ) then
872               
873                  cs_res@tiXAxisString   = "y [m]"
874                  cs_res@tiYAxisString   = "z [m]"
875                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
876
877                  if ( sort .eq. "time" ) then
878                     if ( data&z(lo) .eq. -1 ) then
879                        level = "-average"
880                     else
881                        level = "=" + data&z(lo) + "m"
882                     end if
883                     do n = 0, tb-ta
884                        cs_res@gsnCenterString = "t=" + data&t(li+n-1) + "s  x"+ level
885                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,:,:,lo),cs_res)
886                        if (check_vecp)
887                           vecres                  = True            ; vector only resources
888                           vecres@gsnDraw          = False           ; don't draw
889                           vecres@gsnFrame         = False           ; don't advance frame
890                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
891                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
892                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
893                           vecres@gsnRightString   = " "             ; turn off right string
894                           vecres@gsnLeftString    = " "             ; turn off left string
895                           vecres@tiXAxisString    = " "   
896                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
897                           overlay(plot1, plot_vec)
898                        end if
899                        plot(n) = plot1 
900                     end do
901                  end if
902
903                  if ( sort .eq. "layer" ) then
904                     do n = 0, xb-xa
905                        if ( data&z(li+n) .eq. -1 ) then
906                           level = "-average"
907                        else
908                           level = "=" + data&z(li+n) + "m"
909                        end if
910                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  x"+ level
911                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,:,:,li+n),cs_res)
912                        if (check_vecp)
913                           vecres                  = True            ; vector only resources
914                           vecres@gsnDraw          = False           ; don't draw
915                           vecres@gsnFrame         = False           ; don't advance frame
916                           vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
917                           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
918                           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
919                           vecres@gsnRightString   = " "             ; turn off right string
920                           vecres@gsnLeftString    = " "             ; turn off left string
921                           vecres@tiXAxisString    = " "   
922                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
923                           overlay(plot1, plot_vec)
924                        end if
925                        plot(n) = plot1
926                     end do
927                  end if
928
929               end if
930
931               ; ***************************************************
932               ; merge plots onto one page
933               ; ***************************************************
934         
935               gsn_panel(wks_ps, plot(0:n-1),(/no_lines,no_columns/),cs_resP)
936
937               ; To hold no_lines and no_columns constant, it must be defined again:
938
939               no_lines   = nl1
940               no_columns = nc1
941
942            end do     
943         end do 
944      end if
945   end do   
946
947   print(" ")
948   print("Output to: " + file_out +"."+ format_out)
949   print(" ")   
950 
951end
Note: See TracBrowser for help on using the repository browser.