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

Last change on this file since 2429 was 2147, checked in by scharf, 8 years ago

improved palmplot, added two imuk hosts, changed allocation of resources for lcbullhh and corrected land surface parameters

File size: 79.5 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   
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    = "AutomaticLevels"
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=(decimalPlaces(MaxVal(varn),3,True)-decimalPlaces(MinVal(varn),3,True))/24
1822 
1823         ;cs_res@cnMinLevelValF = decimalPlaces(MinVal(varn),3,True)
1824         ;cs_res@cnMaxLevelValF = decimalPlaces(MaxVal(varn),3,True)
1825
1826         ;cs_res@cnLevelSpacingF = space
1827         cs_res@cnMaxLevelCount = 26
1828     
1829         ; ****************************************************
1830         ; loops over time and layer
1831         ; ****************************************************
1832
1833         
1834         do lo = los, loe                                       
1835            do li = lis, lie
1836
1837               ; ****************************************************
1838               ; xy cross section
1839               ; ****************************************************
1840
1841               if ( xyc .eq. 1 ) then
1842         
1843                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1844                  cs_res@gsnRightString = unit(varn)
1845                 
1846                  if ( sort .eq. "time" ) then
1847                     if ( z_d(lo) .eq. -1)then
1848                        if (delta_z .EQ. -1) then
1849                           level = "-average"
1850                        else
1851                           level = "=" + z_d(lo) + "m"
1852                        end if
1853                     else
1854                        level = "=" + z_d(lo) + "m"
1855                     end if
1856
1857                     if(vNam(varn) .eq. "zwwi"  .or. \
1858                        vNam(varn) .eq. "zusi") then
1859                         loe = 0
1860                         los = 0
1861                         lie  = 0
1862                         lis = 0
1863                         level = ""
1864                     end if
1865
1866                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1867                        vNam(varn) .eq. "pras_xy"  .or. \
1868                        vNam(varn) .eq. "prrs_xy"  .or. \
1869                        vNam(varn) .eq. "qsws_xy"  .or. \
1870                        vNam(varn) .eq. "shfs_xy"  .or. \
1871                        vNam(varn) .eq. "ts_xy"    .or. \
1872                        vNam(varn) .eq. "us_xy"    .or. \
1873                        vNam(varn) .eq. "z0s_xy"   .or. \
1874                        vNam(varn) .eq. "lwp*_xy"  .or. \
1875                        vNam(varn) .eq. "pra*_xy"  .or. \
1876                        vNam(varn) .eq. "prr*_xy"  .or. \
1877                        vNam(varn) .eq. "qsws*_xy" .or. \
1878                        vNam(varn) .eq. "shf*_xy"  .or. \
1879                        vNam(varn) .eq. "t*_xy"    .or. \
1880                        vNam(varn) .eq. "u*_xy"    .or. \
1881                        vNam(varn) .eq. "z0*_xy") then
1882                         loe = 0
1883                         level = "=" + zu1(0) + "m"
1884                      else
1885                         lis = start_time_step
1886                         lie = end_time_step
1887                         los = lays
1888                         loe = laye
1889                     end if                 
1890
1891                     if(vNam(varn) .eq. "zwwi"  .or. \
1892                        vNam(varn) .eq. "zusi") then
1893                         cs_res@gsnCenterString = ""
1894                     else
1895                         cs_res@gsnCenterString = "t=" + \
1896                          decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1897                     end if
1898
1899                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1900                        ;nothing to be done
1901                     else
1902                         plot(n) = gsn_csm_contour(wks_ps,\
1903                                           data(varn,li,lo-los,:,:),cs_res)
1904                     end if
1905
1906
1907                     if (vector .EQ. 1 .AND. check_vecp) then
1908                        vecres                 = True           ; vector only
1909                                                                ; resources
1910                        vecres@gsnDraw         = False          ; don't draw
1911                        vecres@gsnFrame        = False          ; don't
1912                                                                ; advance frame
1913                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly
1914                                                                ; vectors
1915                        vecres@vcRefMagnitudeF = ref_mag        ; define
1916                                                                ; vector ref
1917                                                                ; mag
1918                        vecres@vcRefLengthF    = 0.05           ; define
1919                                                                ; length of
1920                                                                ; vec ref 
1921                        vecres@gsnRightString  = " "
1922                        vecres@gsnLeftString   = " "
1923                        vecres@tiXAxisString   = " "   
1924                        plot_vec=gsn_csm_vector(wks_ps,\
1925                              vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1926                        overlay(plot(n), plot_vec)
1927                     end if                         
1928                  end if
1929                 
1930                  if ( sort .eq. "layer" ) then
1931                     if (z_d(li) .eq. -1) then
1932                        if (delta_z .EQ. -1) then
1933                           level = "-average"
1934                        else
1935                           level = "=" + z_d(li) + "m"
1936                        end if
1937                     else
1938                        level = "=" + z_d(li) + "m"
1939                     end if
1940
1941                     if(vNam(varn) .eq. "zwwi"  .or. \
1942                        vNam(varn) .eq. "zusi") then
1943                         lie = 0
1944                         lis = 0
1945                         loe = 0
1946                         los = 0
1947                         level = ""
1948                      end if
1949
1950                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1951                        vNam(varn) .eq. "pras_xy"  .or. \
1952                        vNam(varn) .eq. "prrs_xy"  .or. \
1953                        vNam(varn) .eq. "qsws_xy"  .or. \
1954                        vNam(varn) .eq. "shfs_xy"  .or. \
1955                        vNam(varn) .eq. "ts_xy"    .or. \
1956                        vNam(varn) .eq. "us_xy"    .or. \
1957                        vNam(varn) .eq. "z0s_xy"   .or. \
1958                        vNam(varn) .eq. "lwp*_xy"  .or. \
1959                        vNam(varn) .eq. "pra*_xy"  .or. \
1960                        vNam(varn) .eq. "prr*_xy"  .or. \
1961                        vNam(varn) .eq. "qsws*_xy" .or. \
1962                        vNam(varn) .eq. "shf*_xy"  .or. \
1963                        vNam(varn) .eq. "t*_xy"    .or. \
1964                        vNam(varn) .eq. "u*_xy"    .or. \
1965                        vNam(varn) .eq. "z0*_xy") then
1966                         lie = 0
1967                         level = "=" + zu1(0) + "m"
1968                     else
1969                         lis = lays
1970                         lie = laye
1971                         los = start_time_step
1972                         loe = end_time_step
1973                     end if
1974
1975                     if(vNam(varn) .eq. "zwwi"  .or. \
1976                        vNam(varn) .eq. "zusi") then
1977                         cs_res@gsnCenterString = ""
1978                     else
1979                         cs_res@gsnCenterString = "t=" + \
1980                          decimalPlaces(t_all(lo)/3600,2,True) +"h  z"+level
1981                     end if
1982
1983                     ;if (topo .EQ. 0) then
1984                        ;nothing to be done
1985                     ;else
1986                     ;   print("'topo' is set to " + topo)
1987
1988                        ;plot(n) = gsn_csm_contour(wks_ps,\
1989                                            ;data(varn,lo,li-lis,:,:),cs_res)
1990                        ;if(vNam(varn) .EQ. "zwwi" .OR. \
1991                         ;  vNam(varn) .EQ. "zusi") then
1992                          ; to_res@cnFillPalette = "GMT_gray"
1993                           ;plot_topo = gsn_csm_contour(wks_ps,\
1994                            ;  data(varn,lo,li-lis,:,:),to_res)
1995                         ;end if
1996                        ;overlay(plot(n), plot_topo)
1997                     ;end if
1998
1999                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2000                        ;nothing to be done
2001                     else
2002                        plot(n) = gsn_csm_contour(wks_ps,\
2003                                            data(varn,lo,li-lis,:,:),cs_res)
2004                     end if
2005
2006                     if (vector .EQ. 1 .AND. check_vecp) then
2007                        vecres                 = True           ; vector only
2008                                                                ; resources
2009                        vecres@gsnDraw         = False          ; don't draw
2010                        vecres@gsnFrame        = False          ; don't
2011                                                                ; advance frame
2012                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2013                        vecres@vcRefMagnitudeF = ref_mag        ; define
2014                                                                ; vector ref
2015                                                                ; mag
2016                        vecres@vcRefLengthF    = 0.05           ; define
2017                                                                ; length of
2018                                                                ; vec ref
2019                        vecres@gsnRightString  = " "            ; turn off
2020                                                                ; right string
2021                        vecres@gsnLeftString   = " "            ; turn off
2022                                                                ; left string
2023                        vecres@tiXAxisString   = " "   
2024                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
2025                                                   vect2(lo,li-lis,:,:),vecres)
2026                        overlay(plot(n), plot_vec)
2027                     end if
2028                  end if
2029               end if
2030
2031               ; ****************************************************
2032               ; xz cross section
2033               ; ****************************************************
2034
2035               if ( xzc .eq. 1 ) then
2036                 
2037                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2038               
2039                  if ( sort .eq. "time" ) then
2040                     if ( y_d(lo) .eq. -1 ) then
2041                        level = "-average"
2042                     else
2043                        level = "=" + y_d(lo) + "m"
2044                     end if
2045                     cs_res@gsnCenterString = "t=" + \
2046                           decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
2047
2048                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2049                        ;nothing to be done
2050                     else
2051                         plot(n) = gsn_csm_contour(wks_ps,\
2052                                              data(varn,li,:,lo-los,:),cs_res)
2053                     end if
2054
2055                     if (vector .EQ. 1 .AND. check_vecp) then
2056                        vecres                 = True           ; vector only
2057                                                                ; resources
2058                        vecres@gsnDraw         = False          ; don't draw
2059                        vecres@gsnFrame        = False          ; don't
2060                                                                ; advance frame
2061                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2062                        vecres@vcRefMagnitudeF = ref_mag        ; define
2063                                                                ; vector ref
2064                                                                ; mag
2065                        vecres@vcRefLengthF    = 0.05           ; define
2066                                                                ; length of
2067                                                                ; vec ref
2068                        vecres@gsnRightString  = " "            ; turn off
2069                                                                ; right string
2070                        vecres@gsnLeftString   = " "            ; turn off
2071                                                                ; left string
2072                        vecres@tiXAxisString   = " "   
2073                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2074                                                   vect2(li,:,lo-los,:),vecres)
2075                        overlay(plot(n), plot_vec)
2076                     end if
2077                  end if
2078
2079                  if ( sort .eq. "layer" ) then
2080                     if ( y_d(li) .eq. -1 ) then
2081                        level = "-average"
2082                     else
2083                        level = "=" + y_d(li) + "m"
2084                     end if
2085                     cs_res@gsnCenterString = "t=" + \
2086                           decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
2087
2088                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2089                        ;nothing to be done
2090                     else
2091                         plot(n) = gsn_csm_contour(wks_ps,\
2092                                              data(varn,lo,:,li-lis,:),cs_res) 
2093                     end if
2094
2095                     if (vector .EQ. 1 .AND. check_vecp) then
2096                        vecres                 = True           ; vector only
2097                                                                ; resources
2098                        vecres@gsnDraw         = False          ; don't draw
2099                        vecres@gsnFrame        = False          ; don't
2100                                                                ; advance frame
2101                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2102                        vecres@vcRefMagnitudeF = ref_mag        ; define
2103                                                                ; vector ref
2104                                                                ; mag
2105                        vecres@vcRefLengthF    = 0.05           ; define
2106                                                                ; length of
2107                                                                ; vec ref
2108                        vecres@gsnRightString  = " "            ; turn off
2109                                                                ; right string
2110                        vecres@gsnLeftString   = " "            ; turn off
2111                                                                ; left string
2112                        vecres@tiXAxisString   = " "   
2113                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2114                                                   vect2(lo,:,li-lis,:),vecres)
2115                        overlay(plot(n), plot_vec)
2116                     end if
2117                  end if                 
2118               end if
2119
2120               ; ****************************************************
2121               ; yz cross section
2122               ; ****************************************************
2123
2124               if ( yzc .eq. 1 ) then
2125                                 
2126                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2127
2128                  if ( sort .eq. "time" ) then
2129                     if ( x_d(lo) .eq. -1 ) then
2130                        level = "-average"
2131                     else
2132                        level = "=" + x_d(lo) + "m"
2133                     end if
2134                     cs_res@gsnCenterString = "t=" + \
2135                          decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
2136
2137                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2138                        ;nothing to be done
2139                     else
2140                        plot(n) = gsn_csm_contour(wks_ps,\
2141                                          data(varn,li,:,:,lo-los),cs_res)
2142                     end if
2143
2144                     if (vector .EQ. 1 .AND. check_vecp) then
2145                        vecres                 = True           ; vector only
2146                                                                ; resources
2147                        vecres@gsnDraw         = False          ; don't draw
2148                        vecres@gsnFrame        = False          ; don't
2149                                                                ; advance frame
2150                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2151                        vecres@vcRefMagnitudeF = ref_mag        ; define
2152                                                                ; vector ref
2153                                                                ; mag
2154                        vecres@vcRefLengthF    = 0.05           ; define
2155                                                                ; length of
2156                                                                ; vec ref
2157                        vecres@gsnRightString  = " "            ; turn off
2158                                                                ; right string
2159                        vecres@gsnLeftString   = " "            ; turn off
2160                                                                ; left string
2161                        vecres@tiXAxisString   = " "   
2162                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2163                                                   vect2(li,:,:,lo-los),vecres)
2164                        overlay(plot(n), plot_vec)
2165                     end if
2166                  end if
2167
2168                  if ( sort .eq. "layer" ) then
2169                     if ( x_d(li) .eq. -1 ) then
2170                        level = "-average"
2171                     else
2172                        level = "=" + x_d(li) + "m"
2173                     end if
2174                     cs_res@gsnCenterString = "t=" + \
2175                           decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
2176
2177                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2178                        ;nothing to be done
2179                     else
2180                        plot(n) = gsn_csm_contour(wks_ps,\
2181                                          data(varn,lo,:,:,li-lis),cs_res)
2182                     end if
2183
2184                     if (vector .EQ. 1 .AND. check_vecp)then
2185                        vecres                 = True           ; vector only
2186                                                                ;resources
2187                        vecres@gsnDraw         = False          ; don't draw
2188                        vecres@gsnFrame        = False          ; don't
2189                                                                ; advance frame
2190                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2191                        vecres@vcRefMagnitudeF = ref_mag        ; define
2192                                                                ; vector ref
2193                                                                ; mag
2194                        vecres@vcRefLengthF    = 0.05           ; define
2195                                                                ; length of
2196                                                                ; vec ref
2197                        vecres@gsnRightString  = " "            ; turn off
2198                                                                ; right string
2199                        vecres@gsnLeftString   = " "            ; turn off
2200                                                                ; left string
2201                        vecres@tiXAxisString   = " "   
2202                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2203                                                   vect2(lo,:,:,li-lis),vecres)
2204                        overlay(plot(n), plot_vec)
2205                     end if
2206                  end if
2207               end if         
2208               n=n+1 
2209            end do     
2210         end do 
2211      end if
2212   end do
2213
2214   ; ***************************************************
2215   ; merge plots onto one page
2216   ; ***************************************************   
2217
2218   no_frames = 0
2219
2220   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2221      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2222           no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) then
2223         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),\
2224                                       (/no_var+1,no_layer*no_time/),cs_resP)
2225         print(" ")
2226         print("Outputs to .eps or .epsi have only one frame")
2227         print(" ")
2228      else
2229         do np = 0,dim_plot-1,no_rows*no_columns   
2230            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2231               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2232                                                                      cs_resP)
2233               no_frames = no_frames + 1
2234            else
2235               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2236                                              (/no_rows,no_columns/),cs_resP)
2237               no_frames = no_frames + 1
2238            end if
2239         end do
2240      end if
2241   else       
2242      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2243                               .AND. dim_plot .gt. no_rows*no_columns) then
2244         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2245         print(" ")
2246         print("Outputs to .eps or .epsi have only one frame")
2247         print(" ")
2248      else
2249         do np = 0,dim_plot-1,no_rows*no_columns 
2250            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2251               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2252                                                                     cs_resP)
2253               no_frames = no_frames + 1
2254            else
2255               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2256                                              (/no_rows,no_columns/),cs_resP)
2257               no_frames = no_frames + 1
2258            end if
2259         end do
2260      end if
2261   end if
2262
2263   if (format_out .EQ. "png" ) then
2264     png_output = new((/no_frames/), string)
2265     j = 0
2266
2267     if (no_frames .eq. 1) then
2268        if (ncl_version .GE. 6 ) then
2269           png_output(0) = file_out+".png"
2270        else
2271           png_output(0) = file_out+".00000"+1+".png"
2272        end if
2273        ;using imagemagick's convert for reducing the white
2274        ;space around the plot
2275        cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2276              png_output(0) + " " + png_output(0)
2277       system(cmd)
2278     else
2279
2280       do i=0, no_frames-1
2281         j = i + 1
2282         if (j .LE. 9) then
2283           png_output(i) = file_out+".00000"+j+".png"
2284         end if
2285         if (j .GT. 9 .AND. j .LE. 99) then
2286           png_output(i) = file_out+".0000"+j+".png"
2287         end if
2288         if (j .GT. 99 .AND. j .LE. 999) then
2289           png_output(i) = file_out+".000"+j+".png"
2290         end if
2291         if (j .GT. 999) then
2292           png_output(i) = file_out+".00"+j+".png"
2293         end if
2294
2295         ;using imagemagick's convert for reducing the white
2296         ;space around the plot
2297         cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2298                png_output(i) + " " + png_output(i)
2299         system(cmd)
2300       end do
2301     end if
2302
2303     print(" ")
2304     print("Output to: "+ png_output)
2305     print(" ")
2306   else
2307     print(" ")
2308     print("Output to: " + file_out +"."+ format_out)
2309     print(" ")
2310   end if
2311
2312end
Note: See TracBrowser for help on using the repository browser.