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

Last change on this file since 4640 was 2837, checked in by gronemeier, 7 years ago

palmplot: config file of ncl scripts can be chosen in shell command

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