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

Last change on this file since 769 was 769, checked in by heinze, 11 years ago

Bugfixes in case of plot of t=0h and plot of topography zusi/zwwi possible

File size: 77.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
7;***************************************************
8; load .ncl.config or .ncl.config.default
9;***************************************************
10
11function which_script()
12local script
13begin
14   script="cross_section"
15   return(script)
16end
17
18if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
19   loadscript("$PALM_BIN/../../.ncl.config")
20else
21  if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
22     loadscript( "$PALM_BIN/NCL/.ncl.config.default")
23  else
24      palm_bin_path = getenv("PALM_BIN")
25      print(" ")
26      print("Neither the personal configuration file '.ncl.config' exists in")
27      print("~/palm/current_version")
28      print("nor the default configuration file '.ncl.config.default' "+\
29            "exists in")
30      print(palm_bin_path + "/NCL")
31      print(" ")
32      exit
33   end if
34end if   
35
36begin
37
38   ;***************************************************
39   ; Retrieving the NCL version used
40   ;***************************************************
41   
42   ncl_version_ch = systemfunc("ncl -V")
43   ncl_version    = stringtofloat(ncl_version_ch)
44
45   ; ***************************************************
46   ; Retrieving the double quote character
47   ; ***************************************************
48   
49   dq=str_get_dq()
50
51   ; ***************************************************
52   ; set default parameter values and strings if not assigned
53   ; in prompt or parameter list
54   ; ***************************************************
55
56   ; Selection of type of cross section is done in palmplot
57   if (.not. isvar("xyc"))then       
58      xyc = 0   
59   end if
60   if (.not. isvar("xzc"))then 
61      xzc = 0   
62   end if
63   if (.not. isvar("yzc"))then     
64      yzc = 0
65   end if
66   
67   if (file_1 .EQ. "File in") then
68      print(" ")
69      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
70      print(" ")
71      exit
72   else
73      file_in = file_1   
74   end if
75
76   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND.   \
77       format_out .NE. "eps" .AND. format_out .NE. "ps" .AND.    \
78       format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
79       format_out .NE. "png") then
80      print(" ")
81      print("'format_out = "+format_out+"' is invalid and set to'x11'")
82      print(" ")
83      format_out="x11"
84   end if
85
86   if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
87      print(" ")
88      print("Output of png files not available")
89      print("png output is avaiable with NCL version 5.2.0 and higher ")
90      print("NCL version used: " + ncl_version_ch)
91      print(" ")
92      exit
93   end if
94
95   if (sort .NE. "layer" .AND. sort .NE. "time") then 
96      print(" ")
97      print("'sort'= "+sort+" is invalid and set to 'layer'")
98      print(" ")
99      sort = "layer" 
100   end if
101   
102   if (mode .NE. "Fill" .AND. mode .NE. "Line" .AND. mode .NE. "Both") then
103      print(" ")
104      print("'mode'= "+mode+" is invalid and set to 'Fill'")
105      print(" ")
106      mode = "Fill"
107   end if
108
109   if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND.\
110       fill_mode .NE. "CellFill") then
111      print(" ")
112      print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
113      print(" ")
114      fill_mode = "AreaFill"
115   end if
116   
117   if (shape .NE. 0 .AND. shape .NE. 1) then
118      print(" ")
119      print("'shape'= "+shape+" is invalid and set to 1")
120      print(" ")
121      shape = 1
122   end if
123
124   if (xyc .NE. 0 .AND. xyc .NE. 1)then
125      print(" ")
126      print("'xyc'= "+xyc+" is invalid and set to 0")
127      print(" ")
128      xyc = 0 
129   end if
130   
131   if (xzc .NE. 0 .AND. xzc .NE. 1)then
132      print(" ")
133      print("'xzc'= "+xzc+" is invalid and set to 0")
134      print(" ")
135      xyc = 0 
136   end if
137   
138   if (yzc .NE. 0 .AND. yzc .NE. 1)then
139      print(" ")
140      print("'yzc'= "+yzc+" is invalid and set to 0")
141      print(" ")
142      xyc = 0 
143   end if
144   
145   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
146      print(" ")
147      print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
148      print(" ")
149      exit
150   end if
151   if (xyc .EQ. 1 ) then
152      if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
153         print(" ")
154         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
155         print(" ")
156         exit
157      end if
158   end if
159   if (xzc .EQ. 1) then
160      if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
161         print(" ")
162         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
163         print(" ")
164         exit
165      end if
166   end if
167   if (yzc .EQ. 1) then
168      if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
169         print(" ")
170         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
171         print(" ")
172         exit
173      end if
174   end if
175
176   if (vector .NE. 0 .AND. vector .NE. 1) then
177      print(" ")
178      print("Set 'vector' to 0 or 1")
179      print(" ")
180      exit
181   end if 
182   
183   if (norm_x .EQ. 0) then
184      print(" ")
185      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
186      print(" ")
187      norm_x = 1.0
188   end if
189   if (norm_y .EQ. 0) then
190      print(" ")
191      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
192      print(" ")
193      norm_y = 1.0
194   end if
195   if (norm_z .EQ. 0) then
196      print(" ")
197      print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
198      print(" ")
199      norm_z = 1.0
200   end if
201 
202   ; ***************************************************
203   ; open input file
204   ; ***************************************************
205   
206   file_in_1 = False
207   if (isStrSubset(file_in, ".nc"))then
208      start_f = -2
209      end_f = -2
210      file_in_1 = True     
211   end if
212
213   if (start_f .EQ. -1)then
214      print(" ")
215      print("'start_f' must be one of the cyclic numbers (at least 0) "+\
216            "of your input file(s)")
217      print(" ") 
218      exit
219   end if
220   if (end_f .EQ. -1)then
221      print(" ")
222      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
223            "your input file(s)")
224      print(" ") 
225      exit
226   end if
227
228   files=new(end_f-start_f+1,string)
229
230   if (file_in_1)then
231      if (isfilepresent(file_in))then
232         files(0)=file_in
233      else
234         print(" ")
235         print("1st input file: '"+file_in+"' does not exist")
236         print(" ")
237         exit
238      end if
239   else   
240      if (start_f .EQ. 0)then
241         if (isfilepresent(file_in+".nc"))then
242            files(0)=file_in+".nc"
243            do i=1,end_f
244               if (isfilepresent(file_in+"."+i+".nc"))then   
245                  files(i)=file_in+"."+i+".nc"
246               else
247                  print(" ")
248                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
249                  print(" ")
250                  exit 
251               end if       
252            end do         
253         else
254            print(" ")
255            print("Input file: '"+file_in+".nc' does not exist")
256            print(" ")
257            exit
258         end if
259      else
260         do i=start_f,end_f
261            if (isfilepresent(file_in+"."+i+".nc"))then   
262               files(i-start_f)=file_in+"."+i+".nc"
263            else
264               print(" ")
265               print("Input file: '"+file_in+"."+i+".nc' does not exist")
266               print(" ")
267               exit 
268            end if
269         end do
270      end if
271   end if
272   
273   f=addfiles(files,"r")
274   f_att=addfile(files(0),"r")
275   ListSetType(f,"cat")
276
277   vNam  = getfilevarnames(f_att)
278   vType = getfilevartypes(f_att,vNam)
279
280   if ((all(vType .eq. "double"))) then ;distinction if data is double or float
281      check_vType = True
282   else
283      check_vType = False
284   end if
285
286   print(" ")
287   print("Variables in input file:")
288   print("- "+ vNam)
289   print(" ")
290   dim   = dimsizes(vNam)
291
292   ; ***************************************************
293   ; check for kind of input file
294   ; ***************************************************
295
296   check_3d = True
297   do varn=0,dim-1
298      checkxy = isStrSubset( vNam(varn),"xy")
299      checkxz = isStrSubset( vNam(varn),"xz")
300      checkyz = isStrSubset( vNam(varn),"yz")
301      if (checkxy .OR. checkxz .OR. checkyz) then
302         check_3d = False
303         break
304      end if
305   end do
306   if (.not. check_3d) then
307      if (xyc .EQ. 1)then
308         if (.not. checkxy) then
309            print(" ")
310            print("Your input file doesn't have values for xy cross-sections")
311            if (checkxz)then
312               print("Select another input file or xz cross-sections")
313            else
314               print("Select another input file or yz cross-sections")
315            end if
316            print(" ")
317            exit
318         else
319            print(" ")
320            print("Your input file contains xy data")
321            print(" ")
322         end if
323      end if
324      if (xzc .EQ. 1)then
325         if (.not. checkxz) then
326            print(" ")
327            print("Your input file doesn't have values for xz cross-sections")
328            if (checkxy)then
329               print("Select another input file or xy cross-sections")
330            else
331               print("Select another input file or yz cross-sections")
332            end if
333            print(" ")
334            exit
335         else
336            print(" ")
337            print("Your input file contains xz data")
338            print(" ")
339         end if
340      end if
341      if (yzc .EQ. 1)then
342         if (.not. checkyz) then
343            print(" ")
344            print("Your input file doesn't have values for yz cross-sections")
345            if (checkxy)then
346               print("Select another input file or xy cross-sections")
347            else
348               print("Select another input file or xz cross-sections")
349            end if
350            print(" ")
351            exit
352         else
353            print(" ")
354            print("Your input file contains yz data")
355            print(" ")
356         end if
357      end if
358   else
359      print(" ")
360      print("Your input file contains 3d or other data")
361      print(" ")
362   end if
363   
364   if (dim .EQ. 0) then
365      print(" ")
366      print("There is no data on file")
367      print(" ")
368      exit
369   end if
370
371   ; ***************************************************
372   ; set up recourses
373   ; ***************************************************
374
375   cs_res                          = True
376   vecres                          = True
377   cs_res@gsnYAxisIrregular2Linear = True 
378   vecres@gsnYAxisIrregular2Linear = True
379 
380   cs_res@gsnDraw                 = False
381   cs_res@gsnFrame                = False
382   cs_res@gsnMaximize                = True
383   
384   cs_res@tmXBLabelFontHeightF   = font_size
385   cs_res@tmYLLabelFontHeightF   = font_size
386   cs_res@tiXAxisFontHeightF     = font_size
387   cs_res@tiYAxisFontHeightF     = font_size
388   cs_res@tmXBMode                ="Automatic"
389   cs_res@tmYLMode                ="Automatic"
390 
391   cs_res@cnLevelSelectionMode    = "ManualLevels"
392   cs_res@lbLabelFontHeightF = font_size_legend
393   cs_res@lbLabelStride = legend_label_stride
394   
395
396   cs_resP = True
397   cs_resP@txString               = f_att@title
398   cs_resP@txFuncCode             = "~"                 
399   cs_resP@txFontHeightF          = 0.015       
400 
401   if ( mode .eq. "Fill" ) then
402      cs_res@cnFillOn             = True
403      cs_res@gsnSpreadColors      = True
404      cs_res@cnFillMode           = fill_mode   
405      cs_res@cnLinesOn            = False       
406      cs_res@cnLineLabelsOn       = False
407   end if
408
409   if ( mode .eq. "Both" ) then
410      cs_res@cnFillOn             = True
411      cs_res@gsnSpreadColors      = True
412      cs_res@cnFillMode           = fill_mode
413      cs_res@cnLinesOn            = True
414      cs_res@cnLineLabelsOn       = True
415   end if
416
417   ; *********************************************
418   ; set up of start_time_step and end_time_step
419   ; *********************************************
420
421   t_all = f[:]->time
422   nt = dimsizes(t_all)
423   delta_t = t_all(nt-1)/nt
424
425   ; ****************************************************       
426   ; start of time step and different types of mistakes that could be done
427   ; ****************************************************
428
429   if (start_time_step .EQ. -1.) then           
430      start_time_step=t_all(0)/3600     
431   else
432      if (start_time_step .GT. t_all(nt-1)/3600)
433         print(" ")
434         print("'start_time_step' = "+ start_time_step +"h is greater than "+\
435               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
436         print(" ")
437         print("Select another 'start_time_step'")
438         print(" ")
439         exit
440      end if
441      if (start_time_step .LT. t_all(0)/3600)
442         print(" ")
443         print("'start_time_step' = "+ start_time_step +"h is lower than "+\
444               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
445         print(" ")
446         print("Select another 'start_time_step'")
447         print(" ")
448         exit
449      end if
450   end if
451   do i=0,nt-1   
452      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
453          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
454         st=i
455         break
456      else
457         st=0
458      end if
459   end do
460
461       
462   ; ****************************************************
463   ; end of time step and different types of mistakes that could be done
464   ; ****************************************************
465
466   if (end_time_step .EQ. -1.) then             
467      end_time_step = t_all(nt-1)/3600 
468   else 
469      if (end_time_step .LT. t_all(0)/3600)
470         print(" ")
471         print("'end_time_step' = "+end_time_step+ "h is lower than "+\
472               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
473         print(" ")
474         print("Select another 'end_time_step'")
475         print(" ")
476         exit
477      end if
478      if (end_time_step .GT. t_all(nt-1)/3600)
479         print(" ")
480         print("'end_time_step' = "+ end_time_step +"h is greater than "+\
481               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
482         print(" ")
483         print("Select another 'end_time_step'") 
484         print(" ")
485         exit
486      end if
487      if (end_time_step .LT. start_time_step/3600)
488         print(" ")
489         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
490               "'start_time_step' = "+start_time_step+"h")
491         print(" ")
492         print("Select another 'start_time_step' or 'end_time_step'")
493         print(" ")
494         exit
495      end if
496   end if
497   do i=0,nt-1     
498      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
499          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
500         et=i
501         break
502      else
503         et=0
504      end if
505   end do
506 
507   delete(start_time_step)
508   start_time_step=round(st,3)
509   delete(end_time_step)
510   end_time_step=round(et,3)
511
512   print(" ")
513   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
514                t_all(start_time_step)+" s => index = "+start_time_step)
515   print("                     till "+t_all(end_time_step)/3600+" h = "  +\
516                t_all(end_time_step)+" s => index = "+end_time_step)
517   print(" ")
518 
519   no_time=(end_time_step-start_time_step)+1
520
521   ; ****************************************************
522   ; set up variables for vector plot if required
523   ; ****************************************************   
524
525   if (vector .EQ. 1) then
526      if (vec1 .EQ. "component1") then
527         print(" ")
528         print("Indicate Vector 1 ('vec1') for Vector-Plot or "+\
529               "set 'vector' to 0")
530         print(" ")
531         exit
532      end if
533      if (vec2 .EQ. "component2") then
534         print(" ")
535         print("Indicate Vector 2 ('vec2') for Vector-Plot or "+\
536               "set 'vector' to 0")
537         print(" ")
538         exit
539      end if           
540   end if
541
542   check_vec1 = False
543   check_vec2 = False
544   check_vecp = False
545
546   ; ****************************************************
547   ; get data for plots
548   ; ****************************************************
549
550   if (xyc .EQ. 1) then
551      do varn=0,dim-1
552         if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then   
553            x_d     = f_att->$vNam(varn)$
554            xdim    = dimsizes(x_d)-1
555            delta_x = x_d(1)-x_d(0)
556            break
557         end if
558      end do
559      do varn=0,dim-1
560         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
561            y_d     = f_att->$vNam(varn)$
562            ydim    = dimsizes(y_d)-1
563            delta_y = y_d(1)-y_d(0)
564            break
565         end if
566      end do
567      do varn=0,dim-1
568         if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
569            z_d     = f_att->$vNam(varn)$
570            zdim    = dimsizes(z_d)-1
571            delta_z = 0
572            break
573         else
574            if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
575               z_d     = f_att->$vNam(varn)$
576               zdim    = dimsizes(z_d)-1
577               delta_z = -1.d
578               break
579            else
580               zdim    = 0
581               delta_z = -1.d
582            end if
583         end if
584         if (vNam(varn) .eq. "zu1_xy") then
585            zu1    = f_att->$vNam(varn)$
586         end if
587      end do
588   end if
589   if (xzc .EQ. 1) then
590      do varn=0,dim-1
591         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then
592            x_d     = f_att->$vNam(varn)$
593            xdim    = dimsizes(x_d)-1
594            delta_x = x_d(1)-x_d(0)
595            break
596         end if
597      end do
598      do varn=0,dim-1   
599         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
600             vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
601            z_d     = f_att->$vNam(varn)$
602            zdim    = dimsizes(z_d)-1
603            delta_z = z_d(1)-z_d(0)
604            break
605         end if
606      end do
607      do varn=0,dim-1
608         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
609            y_d     = f_att->$vNam(varn)$
610            ydim    = dimsizes(y_d)-1
611            delta_y = y_d(1)-y_d(0)
612            break
613         else
614            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
615               y_d     = f_att->$vNam(varn)$
616               ydim    = dimsizes(y_d)-1
617               delta_y = -1.d
618               break
619            else
620               ydim    = 0
621               delta_y = -1.d
622            end if
623         end if
624      end do
625   end if
626   if (yzc .EQ. 1) then
627      do varn=0,dim-1 
628         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
629            y_d     = f_att->$vNam(varn)$
630            ydim    = dimsizes(y_d)-1
631            delta_y = y_d(1)-y_d(0)
632            break
633         end if
634      end do
635      do varn=0,dim-1
636         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
637             vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
638            z_d     = f_att->$vNam(varn)$
639            zdim    = dimsizes(z_d)-1
640            delta_z = z_d(1)-z_d(0)
641            break
642         end if
643      end do
644      do varn=0,dim-1
645         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
646            x_d     = f_att->$vNam(varn)$
647            xdim    = dimsizes(x_d)-1
648            delta_x = x_d(1)-x_d(0)
649            break
650         else
651            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
652               x_d     = f_att->$vNam(varn)$
653               xdim    = dimsizes(x_d)-1
654               delta_x = -1.d
655               break
656            else
657               xdim    = 0
658               delta_x = -1.d
659            end if
660         end if
661      end do
662   end if
663 
664   ; ****************************************************
665   ; set up ranges of x-, y- and z-coordinates
666   ; ****************************************************         
667                   
668   if (xs .EQ. -1.d) then               
669      xs = x_d(0)
670      if (delta_x .EQ. -1) then
671         print(" ")
672         print("You cannot choose a start value for x, "+\
673               "there are preset layers for x")
674         print(" ")
675         xstart=0
676      end if
677   else
678      if (delta_x .EQ. -1) then
679         print(" ")
680         print("You cannot choose a start value for x, "+\
681               "there are preset layers for x")
682         print(" ")
683         xstart=0
684      else
685         if (xs .LT. 0-delta_x/2) then
686            print(" ")
687            print("range start for x-coordinate = "+\
688                   xs+"m is lower than first value = "+(0-delta_x/2)+"m")
689            print(" ")
690            exit
691         end if
692         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
693            if (xs .GE. xdim*delta_x) then
694               print(" ")
695               print("range start for x-coordinate = "+\
696                     xs+"m is equal or greater than last value = "+\
697                     xdim*delta_x+"m")
698               print(" ")
699               exit
700            end if
701         else
702            if (xs .GT. xdim*delta_x) then
703               print(" ")
704               print("range start for x-coordinate = "+\
705                     xs+"m is greater than last value = "+\
706                     xdim*delta_x+"m")
707               print(" ")
708               exit
709            end if
710         end if
711      end if
712   end if
713
714   if (ys .EQ. -1.d) then   
715      ys = y_d(0)
716      if (delta_y .EQ. -1) then
717         print(" ")
718         print("You cannot choose a start value for y, "+\
719               "there are preset layers for y")
720         print(" ")
721         ystart=0
722      end if
723   else
724      if (delta_y .EQ. -1) then
725         print(" ")
726         print("You cannot choose a start value for y, "+\
727                "there are preset layers for y")
728         print(" ")
729         ystart=0
730      else
731         if (ys .LT. 0-delta_y/2) then
732            print(" ")
733            print("range start for y-coordinate = "+\
734                  ys+"m is lower than first value = "+0-delta_y/2+"m")
735            print(" ")
736            exit
737         end if
738         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
739            if (ys .GE. ydim*delta_y) then
740               print(" ")
741               print("range start for y-coordinate = "+\
742                     ys+"m is equal or greater than last value = "+\
743                     ydim*delta_y+"m")
744               print(" ")
745               exit
746            end if
747         else
748            if (ys .GT. ydim*delta_y) then
749               print(" ")
750               print("range start for y-coordinate = "+\
751                     ys+"m is greater than last value = "+ydim*delta_y+"m")
752               print(" ")
753               exit
754            end if
755         end if
756      end if
757   end if
758 
759   if (zs .EQ. -1) then                         
760      zs = 0
761      if (delta_z .EQ. -1) then
762         print(" ")
763         print("You cannot choose a start value for z, "+\
764               "there are preset layers for z")
765         print(" ")
766      end if
767   else
768      if (delta_z .EQ. -1) then
769         print(" ")
770         print("You cannot choose a start value for z, "+\
771               "there are preset layers for z")
772         print(" ")
773         zs = 0
774      else
775         if (zs .LT. 0) then
776            print(" ")
777            print("range start for z-coordinate = "+\
778                  zs+" is lower than first gridpoint = 0")
779            print(" ")
780            exit
781         end if
782         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
783            if (zs .GE. zdim) then
784               print(" ")
785               print("range start for z-coordinate = "+\
786                     zs+" is equal or greater than last gridpoint = "+zdim)
787               print(" ")
788               exit
789            end if
790         else
791            if (zs .GT. zdim) then
792               print(" ")
793               print("range start for z-coordinate = "+\
794                     zs+" is greater than last gridpoint = "+zdim)
795               print(" ")
796               exit
797            end if
798         end if
799      end if
800   end if 
801
802   if (xe .EQ. -1) then   
803      xe = x_d(xdim)
804      if (delta_x .EQ. -1) then
805         print(" ")
806         print("You cannot choose an end value for x, "+\
807               "there are preset layers for x")
808         print(" ")
809         xend=xdim
810      end if
811   else
812      if (delta_x .EQ. -1) then
813         print(" ")
814         print("You cannot choose an end value for x, "+\
815               "there are preset layers for x")
816         print(" ")
817         xend=xdim
818      else
819         if (xe .GT. xdim*delta_x) then
820            print(" ")
821            print("range end for x-coordinate = "+\
822                  xe+"m is greater than last value = "+xdim*delta_x+"m")
823            print(" ")
824            exit
825         end if
826         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
827            if (xe .LE. 0-delta_x/2)
828               print(" ")
829               print("range end for x-coordinate = "+\
830                     xe+"m is equal or lower than first value = "+\
831                     (0-delta_x/2)+"m")
832               print(" ")
833               exit
834            end if
835            if (xe .LE. xs) then
836               print(" ")
837               print("range end for x-coordinate = "+\
838                     xe+"m is equal or lower than start range = "+xs+"m")
839               print(" ")
840               exit
841            end if
842            if (xe .EQ. xs+1) then
843               print(" ")
844               print("range end for x-coordinate = "+\
845                     xe+"m must be at least two more gridpoints "+\
846                     "greater than start range = "+xs+"m")
847               print(" ")
848               exit
849            end if
850         else
851            if (xe .LT. 0-delta_x/2)
852               print(" ")
853               print("range end for x-coordinate = "+\
854                     xe+"m is lower than first value = "+(0-delta_x/2)+"m")
855               print(" ")
856               exit
857            end if
858            if (xe .LT. xs) then
859               print(" ")
860               print("range end for x-coordinate = "+\
861                     xe+"m is lower than start range = "+xs+"m")
862               print(" ")
863               exit
864            end if
865         end if
866      end if               
867   end if
868   
869   if (ye .EQ. -1) then 
870      ye = y_d(ydim)
871      if (delta_y .EQ. -1) then
872         print(" ")
873         print("You cannot choose an end value for y, "+\
874               "there are preset layers for y")
875         print(" ")
876         yend=ydim
877      end if
878   else
879      if (delta_y .EQ. -1) then
880         print(" ")
881         print("You cannot choose an end value for y, "+\
882               "there are preset layers for y")
883         print(" ")
884         yend=ydim
885      else
886         if (ye .GT. ydim*delta_y) then
887            print(" ")
888            print("range end for y-coordinate = "+\
889                  ye+"m is greater than last value = "+ydim*delta_y+"m")
890            print(" ")
891            exit
892         end if
893         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
894            if (ye .LE. 0-delta_y/2)
895               print(" ")
896               print("range end for y-coordinate = "+\
897                     ye+"m is equal or lower than first value = "+\
898                     (0-delta_y/2)+"m")
899               print(" ")
900               exit
901            end if
902            if (ye .LE. ys) then
903               print(" ")
904               print("range end for y-coordinate = "+\
905                    ye+"m is equal or lower than start range = "+ys+"m")
906               print(" ")
907               exit
908            end if
909            if (ye .EQ. ys+1) then
910               print(" ")
911               print("range end for y-coordinate = "+\
912                     ye+"m must be at least two more gridpoints greater "+\
913                     "than start range = "+ys+"m")
914               print(" ")
915               exit
916            end if
917         else
918            if (ye .LT. 0-delta_y/2)
919               print(" ")
920               print("range end for y-coordinate = "+\
921                     ye+"m is lower than first value = "+(0-delta_y/2)+"m")
922               print(" ")
923               exit
924            end if
925            if (ye .LT. ys) then
926               print(" ")
927               print("range end for y-coordinate = "+\
928                     ye+"m is lower than start range = "+ys+"m")
929               print(" ")
930               exit
931            end if
932         end if
933      end if
934   end if
935 
936   if (ze .EQ. -1) then 
937      ze = zdim
938      if (delta_z .EQ. -1) then
939         print(" ")
940         print("You cannot choose an end value for z, "+\
941               "there are preset layers for z")
942         print(" ")
943         ze = zdim
944      end if
945   else
946      if (delta_z .EQ. -1) then
947         print(" ")
948         print("You cannot choose an end value for z, "+\
949               "there are preset layers for z")
950         print(" ")
951         ze = zdim
952      else
953         if (ze .GT. zdim) then
954            print(" ")
955            print("range end for z-coordinate = "+\
956                  ze+" is greater than last gridpoint = "+zdim)
957            print(" ")
958            exit
959         end if
960         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
961            if (ze .LE. 0)
962               print(" ")
963               print("range end for z-coordinate = "+\
964                     ze+" is equal or lower than first gridpoint = 0")
965               print(" ")
966               exit
967            end if
968            if (ze .LE. zs) then
969               print(" ")
970               print("range end for z-coordinate = "+\
971                     ze+" is equal or lower than start range = "+zs)
972               print(" ")
973               exit
974            end if
975            if (ze .EQ. zs+1) then
976               print(" ")
977               print("range end for z-coordinate = "+\
978                     ze+" must be at least two more gridpoints "+\
979                     "greater than start range = "+zs)
980               print(" ")
981               exit
982            end if
983         else
984            if (ze .LT. 0)
985               print(" ")
986               print("range end for z-coordinate = "+\
987                     ze+" is lower than first gridpoint = 0")
988               print(" ")
989               exit
990            end if
991            if (ze .LT. zs) then
992               print(" ")
993               print("range end for z-coordinate = "+\
994                     ze+" is lower than start range = "+zs)
995               print(" ")
996               exit
997            end if
998         end if
999      end if
1000   end if
1001
1002   if (delta_x .NE. -1) then
1003      do i=0,xdim   
1004         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
1005            xstart=i
1006            break
1007         end if
1008      end do
1009      do i=0,xdim   
1010         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
1011            xend=i
1012            break
1013         end if
1014      end do
1015   end if
1016   if (delta_y .NE. -1) then
1017      do i=0,ydim   
1018         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
1019            ystart=i
1020            break
1021         end if
1022      end do
1023      do i=0,ydim   
1024         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
1025            yend=i
1026            break
1027         end if
1028      end do
1029   end if
1030
1031   if( shape .eq. 1 ) then
1032      if (xyc .EQ. 1)then
1033         cs_res@vpWidthF    = (xe-xs)/(ye-ys)       
1034         cs_res@vpHeightF   = 1
1035         vecres@vpWidthF    = (xe-xs)/(ye-ys)       
1036         vecres@vpHeightF   = 1
1037         if (xe-xs .GT. ye-ys)then
1038            if (no_rows .LE. no_columns)then
1039               cs_res@gsnPaperOrientation     = "landscape"
1040               vecres@gsnPaperOrientation     = "landscape"
1041               cs_res@lbOrientation = "Horizontal"
1042            else
1043               cs_res@gsnPaperOrientation     = "portrait"
1044               vecres@gsnPaperOrientation     = "portrait"
1045               cs_res@lbOrientation = "Vertical"
1046            end if 
1047         else
1048            if (no_rows .GE. no_columns)then
1049               cs_res@gsnPaperOrientation     = "portrait"
1050               vecres@gsnPaperOrientation     = "portrait"
1051               cs_res@lbOrientation = "Vertical"
1052            else
1053               cs_res@gsnPaperOrientation     = "landscape"
1054               vecres@gsnPaperOrientation     = "landscape"
1055               cs_res@lbOrientation = "Horizontal"
1056            end if
1057         end if
1058      end if
1059      if (xzc .EQ. 1)then
1060         cs_res@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
1061         cs_res@vpHeightF   = 1
1062         vecres@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
1063         vecres@vpHeightF   = 1
1064         if (xe-xs .GT. (delta_x*(ze-zs)))then
1065            if (no_rows .LE. no_columns)then
1066               cs_res@gsnPaperOrientation     = "landscape"
1067               vecres@gsnPaperOrientation     = "landscape"
1068               cs_res@lbOrientation = "Horizontal"
1069            else
1070               cs_res@gsnPaperOrientation     = "portrait"
1071               vecres@gsnPaperOrientation     = "portrait"
1072               cs_res@lbOrientation = "Vertical"
1073            end if 
1074         else
1075            if (no_rows .GE. no_columns)then
1076               cs_res@gsnPaperOrientation     = "portrait"
1077               vecres@gsnPaperOrientation     = "portrait"
1078               cs_res@lbOrientation = "Vertical"
1079            else
1080               cs_res@gsnPaperOrientation     = "landscape"
1081               vecres@gsnPaperOrientation     = "landscape"
1082               cs_res@lbOrientation = "Horizontal"
1083            end if
1084         end if
1085      end if
1086      if (yzc .EQ. 1)then
1087         cs_res@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
1088         cs_res@vpHeightF   = 1
1089         vecres@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
1090         vecres@vpHeightF   = 1
1091         if (ye-ys .GT. (delta_y*(ze-zs)))then
1092            if (no_rows .LE. no_columns)then
1093               cs_res@gsnPaperOrientation     = "landscape"
1094               vecres@gsnPaperOrientation     = "landscape"
1095               cs_res@lbOrientation = "Horizontal"
1096            else
1097               cs_res@gsnPaperOrientation     = "portrait"
1098               vecres@gsnPaperOrientation     = "portrait"
1099               cs_res@lbOrientation = "Vertical"
1100            end if 
1101         else
1102            if (no_rows .GE. no_columns)then
1103               cs_res@gsnPaperOrientation     = "portrait"
1104               vecres@gsnPaperOrientation     = "portrait"
1105               cs_res@lbOrientation = "Vertical"
1106            else
1107               cs_res@gsnPaperOrientation     = "landscape"
1108               vecres@gsnPaperOrientation     = "landscape"
1109               cs_res@lbOrientation = "Horizontal"
1110            end if
1111         end if
1112      end if
1113   end if
1114
1115   delete(xs)
1116   delete(xe)   
1117   delete(ys)
1118   delete(ye)
1119
1120   xs=xstart
1121   xe=xend
1122   ys=ystart
1123   ye=yend
1124
1125   if (xyc .EQ. 1) then
1126      d= (ye-ys+1)/(major_ticks_y-1)
1127      e=(xe-xs+1)/(major_ticks_x-1)
1128      array_yl       =new(major_ticks_y,integer)
1129      array_yl_labels=new(major_ticks_y,double)
1130      array_minor_yl =new((major_ticks_y-1)*4,double)
1131      array_xb       =new(major_ticks_x,integer)
1132      array_xb_labels=new(major_ticks_x,double)
1133      array_minor_xb =new((major_ticks_x-1)*4,double)
1134      array_yl(0)=ys
1135      array_xb(0)=xs
1136      array_yl_labels(0)=y_d(ys)
1137      array_xb_labels(0)=x_d(xs)
1138      do ar=1,major_ticks_y-1
1139         array_yl(ar)=d*(ar-1)+d-1
1140         array_yl_labels(ar) = y_d(array_yl(ar))
1141         if (ar .GT. 0)
1142            do min_ar=0,3
1143                array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+\
1144                         (y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1)
1145            end do
1146         end if 
1147      end do
1148      do br=1,major_ticks_x-1
1149         array_xb(br)=e*(br-1)+e-1
1150         array_xb_labels(br) = x_d(array_xb(br))
1151         if (br .GT. 0)
1152            do min_br=0,3
1153               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1154                         (x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1155            end do
1156         end if
1157      end do
1158   end if
1159 
1160   if (xzc .EQ. 1) then
1161      d=(ze-zs+1)/(major_ticks_y-1)
1162      e=(xe-xs+1)/(major_ticks_x-1)
1163      array_yl       =new(major_ticks_y,integer)
1164      array_yl_labels=new(major_ticks_y,double)
1165      array_minor_yl =new((major_ticks_y-1)*4,double)
1166      array_xb       =new(major_ticks_x,integer)
1167      array_xb_labels=new(major_ticks_x,double)
1168      array_minor_xb =new((major_ticks_x-1)*4,double)
1169      array_yl(0)=zs
1170      array_xb(0)=xs
1171      array_yl_labels(0)=z_d(zs)
1172      array_xb_labels(0)=x_d(xs)
1173      do ar=1,major_ticks_y-1
1174         array_yl(ar)=d*(ar-1)+d-1
1175         array_yl_labels(ar) = z_d(array_yl(ar))
1176         if (ar .GT. 0)
1177            do min_ar=0,3
1178               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1179                          (z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1180            end do
1181         end if 
1182      end do
1183      do br=1,major_ticks_x-1
1184         array_xb(br)=e*(br-1)+e-1
1185         array_xb_labels(br) = x_d(array_xb(br))
1186         if (br .GT. 0)
1187            do min_br=0,3
1188               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1189                         (x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1190            end do
1191         end if
1192      end do
1193   end if
1194
1195   if (yzc .EQ. 1) then
1196      d=(ze-zs+1)/(major_ticks_y-1)
1197      e=(ye-ys+1)/(major_ticks_x-1)
1198      array_yl       =new(major_ticks_y,integer)
1199      array_yl_labels=new(major_ticks_y,double)
1200      array_minor_yl =new((major_ticks_y-1)*4,double)
1201      array_xb       =new(major_ticks_x,integer)
1202      array_xb_labels=new(major_ticks_x,double)
1203      array_minor_xb =new((major_ticks_x-1)*4,double)
1204      array_yl(0)=zs
1205      array_xb(0)=ys
1206      array_yl_labels(0)=z_d(zs)
1207      array_xb_labels(0)=y_d(ys)
1208      do ar=1,major_ticks_y-1
1209         array_yl(ar)=d*(ar-1)+d-1
1210         array_yl_labels(ar) = z_d(array_yl(ar))
1211         if (ar .GT. 0)
1212            do min_ar=0,3
1213               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1214                         (z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1215            end do
1216         end if 
1217      end do
1218      do br=1,major_ticks_x-1
1219         array_xb(br)=e*(br-1)+e-1
1220         array_xb_labels(br) = y_d(array_xb(br))
1221         if (br .GT. 0)
1222            do min_br=0,3
1223               array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+\
1224                         (y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1)
1225            end do
1226         end if
1227      end do
1228   end if
1229 
1230   if (axes_explicit .EQ. 1)then
1231      cs_res@tmYLMode = "Explicit"
1232      cs_res@tmXBMode = "Explicit"
1233      if (xyc .EQ. 1)then
1234         cs_res@tmYLValues = y_d(array_yl)
1235         cs_res@tmXBValues = x_d(array_xb)     
1236         cs_res@tmYLLabels = array_yl_labels/norm_y
1237         cs_res@tmXBLabels = array_xb_labels/norm_x
1238         if (norm_x .NE. 1.)then
1239            cs_res@tiXAxisString = "x ("+unit_x+")"
1240         else
1241            cs_res@tiXAxisString = "x (m)"
1242         end if
1243         if (norm_y .NE. 1.)then
1244            cs_res@tiYAxisString = "y ("+unit_y+")"
1245         else
1246            cs_res@tiYAxisString = "y (m)"
1247         end if   
1248      end if
1249      if (xzc .EQ. 1)then
1250         cs_res@tmYLValues = z_d(array_yl)
1251         cs_res@tmXBValues = x_d(array_xb) 
1252         cs_res@tmYLLabels = array_yl_labels/norm_z
1253         cs_res@tmXBLabels = array_xb_labels/norm_x
1254         if (norm_x .NE. 1.)then
1255            cs_res@tiXAxisString = "x ("+unit_x+")"
1256         else
1257            cs_res@tiXAxisString = "x (m)"
1258         end if
1259         if (norm_z .NE. 1.)then
1260            cs_res@tiYAxisString = "z ("+unit_z+")"
1261         else
1262            cs_res@tiYAxisString = "z (m)"
1263         end if
1264      end if
1265      if (yzc .EQ. 1)then
1266         cs_res@tmYLValues = z_d(array_yl)
1267         cs_res@tmXBValues = y_d(array_xb) 
1268         cs_res@tmYLLabels = array_yl_labels/norm_z
1269         cs_res@tmXBLabels = array_xb_labels/norm_y
1270         if (norm_y .NE. 1.)then
1271            cs_res@tiXAxisString = "y ("+unit_y+")"
1272         else
1273            cs_res@tiXAxisString = "y (m)"
1274         end if
1275         if (norm_z .NE. 1.)then
1276            cs_res@tiYAxisString = "z ("+unit_z+")"
1277         else
1278            cs_res@tiYAxisString = "z (m)"
1279         end if
1280      end if
1281      cs_res@tmYLMinorValues = array_minor_yl
1282      cs_res@tmXBMinorValues = array_minor_xb
1283   else
1284      if (axes_explicit .EQ. 0)then 
1285         cs_res@tmYLMinorPerMajor = 4
1286         cs_res@tmXBMinorPerMajor = 4
1287      else
1288         print(" ")
1289         print("'axes_explicit' is invalid and set to 0")
1290         print(" ")
1291         axes_explicit = 0
1292         cs_res@tmYLMinorPerMajor = 4
1293         cs_res@tmXBMinorPerMajor = 4
1294      end if
1295      if (xyc .EQ. 1)then
1296         cs_res@tiXAxisString = "x (m)"
1297         cs_res@tiYAxisString = "y (m)"
1298      end if
1299      if (xzc .EQ. 1)then
1300         cs_res@tiXAxisString = "x (m)"
1301         cs_res@tiYAxisString = "z (m)"
1302      end if
1303      if (yzc .EQ. 1)then
1304         cs_res@tiXAxisString = "y (m)"
1305         cs_res@tiYAxisString = "z (m)"
1306      end if
1307   end if
1308
1309   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1310      if (xe .EQ. xs+1) then
1311         print(" ")
1312         print("range end for x-coordinate="+xe*delta_x+\
1313               "m must be at least two")
1314         print("more gridpoints(="+2*delta_x+"m) greater than start range="+\
1315               xs*delta_x+"m)")
1316         print(" ")
1317         exit
1318      end if
1319   end if
1320   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1321      if (ye .EQ. ys+1) then
1322         print(" ")
1323         print("range end for y-coordinate="+ye*delta_y+\
1324               "m must be at least two")
1325         print("more gridpoints(="+2*delta_y+"m greater than start range="+\
1326               ys*delta_y+"m)")
1327         print(" ")
1328         exit
1329      end if
1330   end if
1331   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1332      if (ze .EQ. zs+1) then
1333         print(" ")
1334         print("range end for z-coordinate="+ze+" must be at least two")
1335         print("more gridpoints greater than start range="+zs+")")
1336         print(" ")
1337         exit
1338      end if
1339   end if
1340 
1341   if (xyc .EQ. 1) then
1342      no_layer = (ze-zs)+1
1343      if (check_vType) then
1344         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1345      else
1346         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1347      end if
1348   end if
1349   if (xzc .EQ. 1) then
1350      no_layer = (ye-ys)+1
1351      if (check_vType) then
1352         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1353      else
1354         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1355      end if
1356   end if
1357   if (yzc .EQ. 1) then
1358      no_layer = (xe-xs)+1
1359      if (check_vType) then
1360         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1361      else
1362         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1363      end if
1364   end if
1365
1366   if (check_vType) then
1367      MinVal = new(dim,double)
1368      MaxVal = new(dim,double)
1369   else
1370      MinVal = new(dim,float)
1371      MaxVal = new(dim,float)
1372   end if
1373   unit   = new(dim,string)
1374
1375   ; ****************************************************
1376   ; define inner and outer loops depending on "sort"
1377   ; ****************************************************       
1378   
1379   if ( xyc .eq. 1 ) then
1380      lays = zs
1381      laye = ze
1382   end if
1383   if ( xzc .eq. 1 ) then
1384      lays = ys
1385      laye = ye
1386   end if
1387   if ( yzc .eq. 1) then
1388      lays = xs
1389      laye = xe
1390   end if
1391
1392   if ( sort .eq. "time" ) then
1393      lis = start_time_step
1394      lie = end_time_step
1395      los = lays
1396      loe = laye
1397   else
1398      lis = lays
1399      lie = laye
1400      los = start_time_step
1401      loe = end_time_step
1402   end if
1403
1404   check  = True
1405   v1     = 0
1406   v2     = 0
1407   no_var = 0
1408   n      = 0
1409   no_zu1 = 0
1410   no_topo= 0
1411
1412   ;****************************************************
1413   ; Preparation of vector plots
1414   ; ****************************************************       
1415
1416   do varn=dim-1,0,1
1417
1418     if (vector .EQ. 1) then   
1419        check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1420        if (check_vec1) then
1421            temp = f[:]->$vNam(varn)$
1422            data_att = f_att->$vNam(varn)$
1423            vect1=temp(:,zs:ze,ys:ye,xs:xe)
1424            delete(temp)
1425         end if
1426
1427         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1428         if (check_vec2) then
1429            temp = f[:]->$vNam(varn)$
1430            data_att = f_att->$vNam(varn)$
1431            vect2=temp(:,zs:ze,ys:ye,xs:xe)
1432            delete(temp)
1433         end if
1434   
1435         if (","+vNam(varn)+"," .EQ. vec1) then
1436            v1=v1+1
1437         end if
1438         if (","+vNam(varn)+"," .EQ. vec2) then
1439             v2=v2+1
1440          end if
1441     end if
1442
1443   end do
1444
1445
1446
1447   do varn=dim-1,0,1     
1448           
1449      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1450           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1451           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1452           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1453           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1454           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1455           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1456           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1457           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1458           vNam(varn) .eq. "ind_x_yz") then
1459         check = False
1460      else
1461         check = True
1462      end if 
1463
1464      if (var .NE. "all") then
1465         check = isStrSubset( var,","+vNam(varn)+"," )
1466      end if
1467
1468      if(check) then
1469     
1470         no_var=no_var+1
1471
1472         if (xyc .EQ. 1) then
1473            temp = f[:]->$vNam(varn)$
1474            if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1475               dummy=0
1476            else
1477               data_att = f_att->$vNam(varn)$
1478            end if
1479            if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy"     \
1480              .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1481              .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy"   \
1482              .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy"    \
1483              .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \
1484              .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\
1485              .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy"   \
1486              .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then
1487              ;these variables depend on zu1_xy and that's why they have
1488              ;only one z-layer
1489              data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1490              no_zu1=no_zu1+1
1491            else
1492              if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1493                  ;these variables depend on x and y
1494                  data(varn,0,0,:,:)=doubletofloat(temp(ys:ye,xs:xe))
1495                  no_topo=no_topo+1
1496               else
1497                  data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1498               end if
1499            end if
1500            delete(temp)               
1501         end if
1502         if ( xzc .eq. 1 ) then
1503            temp = f[:]->$vNam(varn)$
1504            data_att = f_att->$vNam(varn)$
1505            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1506            delete(temp)
1507         end if
1508         if ( yzc .eq. 1) then
1509            temp = f[:]->$vNam(varn)$
1510            data_att = f_att->$vNam(varn)$
1511            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1512            delete(temp)
1513         end if
1514
1515         data!0 = "var"
1516         data!1 = "t"
1517         data!2 = "z"
1518         data!3 = "y"
1519         data!4 = "x" 
1520
1521         MinVal(varn) = min(data(varn,start_time_step:end_time_step,\
1522                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1523         MaxVal(varn) = max(data(varn,start_time_step:end_time_step,\
1524                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1525         
1526         if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1527            unit(varn) = "meters"
1528          else
1529            unit(varn) = data_att@units
1530            delete(data_att)
1531          end if
1532
1533      end if
1534
1535   end do
1536
1537   if (no_var .EQ. 0) then
1538      print(" ")
1539      print("The variables var='"+var+"' do not exist on your input file;")
1540      print("be sure to have one comma before and after each variable")
1541      print(" ")
1542      exit
1543   end if
1544
1545   var_input=new(no_var,string)
1546   numb=0
1547   no_var=0
1548   
1549   do varn=dim-1,0,1   
1550           
1551      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1552           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1553           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1554           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1555           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1556           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1557           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1558           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1559           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1560           vNam(varn) .eq. "ind_x_yz") then
1561         check = False
1562      else
1563         check = True
1564      end if 
1565
1566      if (var .NE. "all") then
1567         check = isStrSubset( var,","+vNam(varn)+"," )
1568      end if
1569
1570      if(check) then     
1571         var_input(numb)=vNam(varn)
1572         numb=numb+1     
1573      end if
1574     
1575      if (check) then
1576         no_var = no_var+1
1577      end if
1578     
1579   end do
1580
1581   if (no_var .EQ. 0) then
1582      print(" ")
1583      print("The variables var='"+var+"' do not exist on your input file;")
1584      print("be sure to have one comma before and after each variable")
1585      print(" ")
1586      exit
1587   end if
1588
1589   if (vector .EQ. 1) then
1590      if (v1 .EQ. 0)then
1591         print(" ")
1592         print("Component 1 for the vector-plot ('vec1') must be one "+\
1593               "of the variables on the input file:")
1594         print("- "+var_input)
1595         print("be sure to have one comma before and after the variable")
1596         print(" ")
1597         exit
1598      end if
1599
1600      if (v2 .EQ. 0)then
1601         print(" ")
1602         print("Component 2 for the vector-plot ('vec2') must be one "+\
1603               "of the variables on the input file:")
1604         print("- "+var_input)
1605         print("be sure to have one comma before and after the variable")
1606         print(" ")
1607         exit
1608      end if
1609   end if
1610
1611   ; ***************************************************
1612   ; open workstation(s)
1613   ; ***************************************************
1614
1615   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1616      format_out@wkPaperSize = "A4"
1617   end if
1618   if ( format_out .EQ. "png" ) then
1619      format_out@wkWidth  = 1000
1620      format_out@wkHeight = 1000
1621   end if
1622
1623   wks_ps  = gsn_open_wks(format_out,file_out)
1624   gsn_define_colormap(wks_ps,"rainbow+white")
1625
1626   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1627      plot=new((/no_time*no_layer/),graphic)
1628   else
1629      plot=new((/no_time*no_layer*(no_var-no_topo) + no_topo - \
1630                       no_time*(no_layer-1)*no_zu1/),graphic)
1631   end if
1632   dim_plot=dimsizes(plot)
1633
1634   page = 0
1635   n=0
1636   print(" ")
1637   print("Plots sorted by '"+sort+"'")
1638   print(" ")
1639
1640   ; ***************************************************
1641   ; create plots
1642   ; ***************************************************
1643
1644
1645   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1646      do lo = los, loe                                         
1647         do li = lis, lie
1648            if (xyc .EQ. 1)then
1649               if (sort .EQ. "time")then
1650                  if(z_d(lo) .eq. -1.d) then
1651                    level = "z-average"
1652                  else
1653                    level = "z=" + z_d(lo) + "m"
1654                  end if
1655               else
1656                  if(z_d(li) .eq. -1.d) then
1657                    level = "z-average"
1658                  else
1659                    level = "z=" + z_d(li) + "m"
1660                  end if
1661               end if
1662            end if
1663            if (xzc .EQ. 1)then
1664               if (sort .EQ. "time")then
1665                 if(y_d(lo) .eq. -1.d) then
1666                   level = "y-average"
1667                 else
1668                   level = "y=" + y_d(lo) + "m"
1669                 end if
1670               else
1671                  if(y_d(li) .eq. -1.d) then
1672                    level = "y-average"
1673                  else
1674                    level = "y=" + y_d(li) + "m"
1675                  end if
1676               end if
1677            end if
1678            if (yzc .EQ. 1)then
1679               if (sort .EQ. "time")then
1680                 if(x_d(lo) .eq. -1.d) then
1681                    level = "x-average"
1682                 else
1683                    level = "x=" + x_d(lo) + "m"
1684                 end if
1685               else
1686                 if(x_d(li) .eq. -1.d) then
1687                    level = "x-average"
1688                 else
1689                    level = "x=" + x_d(li) + "m"
1690                 end if
1691               end if
1692            end if               
1693            vecres                  = True            ; vector only resources
1694            vecres@gsnDraw          = False           ; don't draw
1695            vecres@gsnFrame         = False           ; don't advance frame
1696            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1697            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1698            vecres@vcRefLengthF     = 0.05            ; define length of
1699                                                      ; vec ref
1700            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1701            vecres@tmXBLabelFontHeightF   = font_size
1702            vecres@tmYLLabelFontHeightF   = font_size
1703            vecres@tiXAxisFontHeightF     = font_size
1704            vecres@tiYAxisFontHeightF     = font_size
1705            if (sort .EQ. "time")then
1706               vecres@gsnLeftString = "t=" + \
1707                             decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1708            else
1709               vecres@gsnLeftString = "t=" + \
1710                             decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1711            end if
1712            if (xyc .EQ. 1) then 
1713               vecres@tiXAxisString = "x (m)"
1714               vecres@tiYAxisString = "y (m)"   
1715               if (sort .EQ. "time")then
1716                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
1717                                               vect2(li,lo-los,:,:),vecres)
1718               else
1719                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1720                                                vect2(lo,li-lis,:,:),vecres)
1721               end if
1722            end if
1723            if (xzc .EQ. 1) then
1724               vecres@tiXAxisString = "x (m)"
1725               vecres@tiYAxisString = "z (m)"
1726               if (sort .EQ. "time")then
1727                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1728                                               vect2(li,:,lo-los,:),vecres)
1729               else
1730                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
1731                                                vect2(lo,:,li-lis,:),vecres)
1732               end if
1733            end if
1734            if (yzc .EQ. 1) then
1735               vecres@tiXAxisString = "y (m)"
1736               vecres@tiYAxisString = "z (m)"
1737               if (sort .EQ. "time")then
1738                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
1739                                                vect2(li,:,:,lo-los),vecres)
1740               else
1741                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
1742                                                vect2(lo,:,:,li-lis),vecres)
1743               end if
1744            end if
1745            n=n+1
1746         end do
1747      end do
1748   end if
1749 
1750   do varn=dim-1,0,1
1751
1752      if (vector .EQ. 1 ) then   
1753         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1754      end if
1755     
1756      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1757           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1758           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1759           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1760           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1761           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1762           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1763           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1764           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1765           vNam(varn) .eq. "ind_x_yz") then
1766         check = False
1767      else
1768         check = True
1769      end if   
1770 
1771      if (var .NE. "all") then
1772         check = isStrSubset( var,","+vNam(varn)+"," )
1773      end if
1774   
1775      if(check) then
1776
1777         space=(MaxVal(varn)-MinVal(varn))/24
1778 
1779         cs_res@cnMinLevelValF = MinVal(varn)
1780         cs_res@cnMaxLevelValF = MaxVal(varn)
1781
1782         cs_res@cnLevelSpacingF = space
1783     
1784         ; ****************************************************
1785         ; loops over time and layer
1786         ; ****************************************************
1787
1788         
1789         do lo = los, loe                                       
1790            do li = lis, lie
1791
1792               ; ****************************************************
1793               ; xy cross section
1794               ; ****************************************************
1795
1796               if ( xyc .eq. 1 ) then
1797         
1798                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1799                  cs_res@gsnRightString = unit(varn)
1800                 
1801                  if ( sort .eq. "time" ) then
1802                     if ( z_d(lo) .eq. -1)then
1803                        if (delta_z .EQ. -1) then
1804                           level = "-average"
1805                        else
1806                           level = "=" + z_d(lo) + "m"
1807                        end if
1808                     else
1809                        level = "=" + z_d(lo) + "m"
1810                     end if
1811
1812                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1813                        vNam(varn) .eq. "pras_xy"  .or. \
1814                        vNam(varn) .eq. "prrs_xy"  .or. \
1815                        vNam(varn) .eq. "qsws_xy"  .or. \
1816                        vNam(varn) .eq. "shfs_xy"  .or. \
1817                        vNam(varn) .eq. "ts_xy"    .or. \
1818                        vNam(varn) .eq. "us_xy"    .or. \
1819                        vNam(varn) .eq. "z0s_xy"   .or. \
1820                        vNam(varn) .eq. "lwp*_xy"  .or. \
1821                        vNam(varn) .eq. "pra*_xy"  .or. \
1822                        vNam(varn) .eq. "prr*_xy"  .or. \
1823                        vNam(varn) .eq. "qsws*_xy" .or. \
1824                        vNam(varn) .eq. "shf*_xy"  .or. \
1825                        vNam(varn) .eq. "t*_xy"    .or. \
1826                        vNam(varn) .eq. "u*_xy"    .or. \
1827                        vNam(varn) .eq. "z0*_xy") then
1828                         loe = 0
1829                         level = "=" + zu1(0) + "m"
1830                      else
1831                         loe = laye
1832                     end if
1833
1834                     if(vNam(varn) .eq. "zwwi"  .or. \
1835                        vNam(varn) .eq. "zusi") then
1836                         loe = 0
1837                         los = 0
1838                         lie  = 0
1839                         lis = 0
1840                         level = ""
1841                      else
1842                         lis = start_time_step
1843                         lie = end_time_step
1844                         los = lays
1845                         loe = laye
1846                     end if
1847                   
1848                     if(vNam(varn) .eq. "zwwi"  .or. \
1849                        vNam(varn) .eq. "zusi") then
1850                         cs_res@gsnCenterString = ""
1851                     else
1852                         cs_res@gsnCenterString = "t=" + \
1853                          decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1854                     end if
1855
1856                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1857                        ;nothing to be done
1858                     else
1859                         plot(n) = gsn_csm_contour(wks_ps,\
1860                                           data(varn,li,lo-los,:,:),cs_res)
1861                     end if
1862
1863
1864                     if (vector .EQ. 1 .AND. check_vecp) then
1865                        vecres                 = True           ; vector only
1866                                                                ; resources
1867                        vecres@gsnDraw         = False          ; don't draw
1868                        vecres@gsnFrame        = False          ; don't
1869                                                                ; advance frame
1870                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly
1871                                                                ; vectors
1872                        vecres@vcRefMagnitudeF = ref_mag        ; define
1873                                                                ; vector ref
1874                                                                ; mag
1875                        vecres@vcRefLengthF    = 0.05           ; define
1876                                                                ; length of
1877                                                                ; vec ref 
1878                        vecres@gsnRightString  = " "
1879                        vecres@gsnLeftString   = " "
1880                        vecres@tiXAxisString   = " "   
1881                        plot_vec=gsn_csm_vector(wks_ps,\
1882                              vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1883                        overlay(plot(n), plot_vec)
1884                     end if                         
1885                  end if
1886                 
1887                  if ( sort .eq. "layer" ) then
1888                     if (z_d(li) .eq. -1) then
1889                        if (delta_z .EQ. -1) then
1890                           level = "-average"
1891                        else
1892                           level = "=" + z_d(li) + "m"
1893                        end if
1894                     else
1895                        level = "=" + z_d(li) + "m"
1896                     end if
1897       
1898                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1899                        vNam(varn) .eq. "pras_xy"  .or. \
1900                        vNam(varn) .eq. "prrs_xy"  .or. \
1901                        vNam(varn) .eq. "qsws_xy"  .or. \
1902                        vNam(varn) .eq. "shfs_xy"  .or. \
1903                        vNam(varn) .eq. "ts_xy"    .or. \
1904                        vNam(varn) .eq. "us_xy"    .or. \
1905                        vNam(varn) .eq. "z0s_xy"   .or. \
1906                        vNam(varn) .eq. "lwp*_xy"  .or. \
1907                        vNam(varn) .eq. "pra*_xy"  .or. \
1908                        vNam(varn) .eq. "prr*_xy"  .or. \
1909                        vNam(varn) .eq. "qsws*_xy" .or. \
1910                        vNam(varn) .eq. "shf*_xy"  .or. \
1911                        vNam(varn) .eq. "t*_xy"    .or. \
1912                        vNam(varn) .eq. "u*_xy"    .or. \
1913                        vNam(varn) .eq. "z0*_xy") then
1914                         lie = 0
1915                         level = "=" + zu1(0) + "m"
1916                     else
1917                         lie = laye
1918                     end if
1919
1920                     if(vNam(varn) .eq. "zwwi"  .or. \
1921                        vNam(varn) .eq. "zusi") then
1922                         lie = 0
1923                         lis = 0
1924                         loe = 0
1925                         los = 0
1926                         level = ""
1927                     else
1928                         lis = lays
1929                         lie = laye
1930                         los = start_time_step
1931                         loe = end_time_step
1932                     end if
1933
1934                     if(vNam(varn) .eq. "zwwi"  .or. \
1935                        vNam(varn) .eq. "zusi") then
1936                         cs_res@gsnCenterString = ""
1937                     else
1938                         cs_res@gsnCenterString = "t=" + \
1939                          decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1940                     end if
1941
1942                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1943                        ;nothing to be done
1944                     else
1945                        plot(n) = gsn_csm_contour(wks_ps,\
1946                                            data(varn,lo,li-lis,:,:),cs_res)
1947                     end if
1948
1949                     if (vector .EQ. 1 .AND. check_vecp) then
1950                        vecres                 = True           ; vector only
1951                                                                ; resources
1952                        vecres@gsnDraw         = False          ; don't draw
1953                        vecres@gsnFrame        = False          ; don't
1954                                                                ; advance frame
1955                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
1956                        vecres@vcRefMagnitudeF = ref_mag        ; define
1957                                                                ; vector ref
1958                                                                ; mag
1959                        vecres@vcRefLengthF    = 0.05           ; define
1960                                                                ; length of
1961                                                                ; vec ref
1962                        vecres@gsnRightString  = " "            ; turn off
1963                                                                ; right string
1964                        vecres@gsnLeftString   = " "            ; turn off
1965                                                                ; left string
1966                        vecres@tiXAxisString   = " "   
1967                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1968                                                   vect2(lo,li-lis,:,:),vecres)
1969                        overlay(plot(n), plot_vec)
1970                     end if
1971                  end if
1972               end if
1973
1974               ; ****************************************************
1975               ; xz cross section
1976               ; ****************************************************
1977
1978               if ( xzc .eq. 1 ) then
1979                 
1980                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1981               
1982                  if ( sort .eq. "time" ) then
1983                     if ( y_d(lo) .eq. -1 ) then
1984                        level = "-average"
1985                     else
1986                        level = "=" + y_d(lo) + "m"
1987                     end if
1988                     cs_res@gsnCenterString = "t=" + \
1989                           decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
1990
1991                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1992                        ;nothing to be done
1993                     else
1994                         plot(n) = gsn_csm_contour(wks_ps,\
1995                                              data(varn,li,:,lo-los,:),cs_res)
1996                     end if
1997
1998                     if (vector .EQ. 1 .AND. check_vecp) then
1999                        vecres                 = True           ; vector only
2000                                                                ; resources
2001                        vecres@gsnDraw         = False          ; don't draw
2002                        vecres@gsnFrame        = False          ; don't
2003                                                                ; advance frame
2004                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2005                        vecres@vcRefMagnitudeF = ref_mag        ; define
2006                                                                ; vector ref
2007                                                                ; mag
2008                        vecres@vcRefLengthF    = 0.05           ; define
2009                                                                ; length of
2010                                                                ; vec ref
2011                        vecres@gsnRightString  = " "            ; turn off
2012                                                                ; right string
2013                        vecres@gsnLeftString   = " "            ; turn off
2014                                                                ; left string
2015                        vecres@tiXAxisString   = " "   
2016                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2017                                                   vect2(li,:,lo-los,:),vecres)
2018                        overlay(plot(n), plot_vec)
2019                     end if
2020                  end if
2021
2022                  if ( sort .eq. "layer" ) then
2023                     if ( y_d(li) .eq. -1 ) then
2024                        level = "-average"
2025                     else
2026                        level = "=" + y_d(li) + "m"
2027                     end if
2028                     cs_res@gsnCenterString = "t=" + \
2029                           decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
2030
2031                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2032                        ;nothing to be done
2033                     else
2034                         plot(n) = gsn_csm_contour(wks_ps,\
2035                                              data(varn,lo,:,li-lis,:),cs_res) 
2036                     end if
2037
2038                     if (vector .EQ. 1 .AND. check_vecp) then
2039                        vecres                 = True           ; vector only
2040                                                                ; resources
2041                        vecres@gsnDraw         = False          ; don't draw
2042                        vecres@gsnFrame        = False          ; don't
2043                                                                ; advance frame
2044                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2045                        vecres@vcRefMagnitudeF = ref_mag        ; define
2046                                                                ; vector ref
2047                                                                ; mag
2048                        vecres@vcRefLengthF    = 0.05           ; define
2049                                                                ; length of
2050                                                                ; vec ref
2051                        vecres@gsnRightString  = " "            ; turn off
2052                                                                ; right string
2053                        vecres@gsnLeftString   = " "            ; turn off
2054                                                                ; left string
2055                        vecres@tiXAxisString   = " "   
2056                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2057                                                   vect2(lo,:,li-lis,:),vecres)
2058                        overlay(plot(n), plot_vec)
2059                     end if
2060                  end if                 
2061               end if
2062
2063               ; ****************************************************
2064               ; yz cross section
2065               ; ****************************************************
2066
2067               if ( yzc .eq. 1 ) then
2068                                 
2069                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2070
2071                  if ( sort .eq. "time" ) then
2072                     if ( x_d(lo) .eq. -1 ) then
2073                        level = "-average"
2074                     else
2075                        level = "=" + x_d(lo) + "m"
2076                     end if
2077                     cs_res@gsnCenterString = "t=" + \
2078                          decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
2079
2080                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2081                        ;nothing to be done
2082                     else
2083                        plot(n) = gsn_csm_contour(wks_ps,\
2084                                          data(varn,li,:,:,lo-los),cs_res)
2085                     end if
2086
2087                     if (vector .EQ. 1 .AND. check_vecp) then
2088                        vecres                 = True           ; vector only
2089                                                                ; resources
2090                        vecres@gsnDraw         = False          ; don't draw
2091                        vecres@gsnFrame        = False          ; don't
2092                                                                ; advance frame
2093                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2094                        vecres@vcRefMagnitudeF = ref_mag        ; define
2095                                                                ; vector ref
2096                                                                ; mag
2097                        vecres@vcRefLengthF    = 0.05           ; define
2098                                                                ; length of
2099                                                                ; vec ref
2100                        vecres@gsnRightString  = " "            ; turn off
2101                                                                ; right string
2102                        vecres@gsnLeftString   = " "            ; turn off
2103                                                                ; left string
2104                        vecres@tiXAxisString   = " "   
2105                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2106                                                   vect2(li,:,:,lo-los),vecres)
2107                        overlay(plot(n), plot_vec)
2108                     end if
2109                  end if
2110
2111                  if ( sort .eq. "layer" ) then
2112                     if ( x_d(li) .eq. -1 ) then
2113                        level = "-average"
2114                     else
2115                        level = "=" + x_d(li) + "m"
2116                     end if
2117                     cs_res@gsnCenterString = "t=" + \
2118                           decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
2119
2120                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2121                        ;nothing to be done
2122                     else
2123                        plot(n) = gsn_csm_contour(wks_ps,\
2124                                          data(varn,lo,:,:,li-lis),cs_res)
2125                     end if
2126
2127                     if (vector .EQ. 1 .AND. check_vecp)then
2128                        vecres                 = True           ; vector only
2129                                                                ;resources
2130                        vecres@gsnDraw         = False          ; don't draw
2131                        vecres@gsnFrame        = False          ; don't
2132                                                                ; advance frame
2133                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2134                        vecres@vcRefMagnitudeF = ref_mag        ; define
2135                                                                ; vector ref
2136                                                                ; mag
2137                        vecres@vcRefLengthF    = 0.05           ; define
2138                                                                ; length of
2139                                                                ; vec ref
2140                        vecres@gsnRightString  = " "            ; turn off
2141                                                                ; right string
2142                        vecres@gsnLeftString   = " "            ; turn off
2143                                                                ; left string
2144                        vecres@tiXAxisString   = " "   
2145                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2146                                                   vect2(lo,:,:,li-lis),vecres)
2147                        overlay(plot(n), plot_vec)
2148                     end if
2149                  end if
2150               end if         
2151               n=n+1 
2152            end do     
2153         end do 
2154      end if
2155   end do
2156
2157   ; ***************************************************
2158   ; merge plots onto one page
2159   ; ***************************************************   
2160
2161   no_frames = 0
2162
2163   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2164      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2165           no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) then
2166         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),\
2167                                       (/no_var+1,no_layer*no_time/),cs_resP)
2168         print(" ")
2169         print("Outputs to .eps or .epsi have only one frame")
2170         print(" ")
2171      else
2172         do np = 0,dim_plot-1,no_rows*no_columns   
2173            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2174               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2175                                                                      cs_resP)
2176               no_frames = no_frames + 1
2177            else
2178               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2179                                              (/no_rows,no_columns/),cs_resP)
2180               no_frames = no_frames + 1
2181            end if
2182         end do
2183      end if
2184   else       
2185      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2186                               .AND. dim_plot .gt. no_rows*no_columns) then
2187         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2188         print(" ")
2189         print("Outputs to .eps or .epsi have only one frame")
2190         print(" ")
2191      else
2192         do np = 0,dim_plot-1,no_rows*no_columns 
2193            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2194               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2195                                                                     cs_resP)
2196               no_frames = no_frames + 1
2197            else
2198               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2199                                              (/no_rows,no_columns/),cs_resP)
2200               no_frames = no_frames + 1
2201            end if
2202         end do
2203      end if
2204   end if
2205
2206   if (format_out .EQ. "png" ) then
2207     png_output = new((/no_frames/), string)
2208     j = 0
2209     do i=0, no_frames-1
2210       j = i + 1
2211       if (j .LE. 9) then
2212         png_output(i) = file_out+".00000"+j+".png"
2213       end if
2214       if (j .GT. 9 .AND. j .LE. 99) then
2215         png_output(i) = file_out+".0000"+j+".png"
2216       end if
2217       if (j .GT. 99 .AND. j .LE. 999) then
2218         png_output(i) = file_out+".000"+j+".png"
2219       end if
2220       if (j .GT. 999) then
2221         png_output(i) = file_out+".00"+j+".png"
2222       end if
2223
2224       ;using imagemagick's convert for reducing the white
2225       ;space around the plot
2226       cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2227              png_output(i) + " " + png_output(i)
2228       system(cmd)
2229     end do
2230
2231     print(" ")
2232     print("Output to: "+ png_output)
2233     print(" ")
2234   else
2235     print(" ")
2236     print("Output to: " + file_out +"."+ format_out)
2237     print(" ")
2238   end if
2239 
2240end
Note: See TracBrowser for help on using the repository browser.