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

Last change on this file since 268 was 267, checked in by weinreis, 16 years ago

NCL scripts profiles cross_sections ncl_preference and UserGuide? modified

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