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

Last change on this file since 1914 was 1626, checked in by heinze, 9 years ago

Bugfix to plot xy-cross sections in case variables contain only one level (like lwp*_xy, shf*_xy, ...)

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