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

Last change on this file since 196 was 195, checked in by steinfeld, 16 years ago

Bug fix in the NCL scripts

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