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

Last change on this file since 377 was 358, checked in by heinze, 15 years ago

Bugfix in cross_sections.ncl concerning plot layers and vector plots fixed.

File size: 62.0 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         if (vNam(varn) .eq. "zu1_xy") then
533            zu1    = f_att->$vNam(varn)$
534         end if
535      end do
536   end if
537   if (xzc .EQ. 1) then
538      do varn=0,dim-1
539         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then
540            x_d     = f_att->$vNam(varn)$
541            xdim    = dimsizes(x_d)-1
542            delta_x = x_d(1)-x_d(0)
543            break
544         end if
545      end do
546      do varn=0,dim-1   
547         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
548            z_d     = f_att->$vNam(varn)$
549            zdim    = dimsizes(z_d)-1
550            delta_z = z_d(1)-z_d(0)
551            break
552         end if
553      end do
554      do varn=0,dim-1
555         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
556            y_d     = f_att->$vNam(varn)$
557            ydim    = dimsizes(y_d)-1
558            delta_y = y_d(1)-y_d(0)
559            break
560         else
561            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
562               y_d     = f_att->$vNam(varn)$
563               ydim    = dimsizes(y_d)-1
564               delta_y = -1.d
565               break
566            else
567               ydim    = 0
568               delta_y = -1.d
569            end if
570         end if
571      end do
572   end if
573   if (yzc .EQ. 1) then
574      do varn=0,dim-1 
575         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
576            y_d     = f_att->$vNam(varn)$
577            ydim    = dimsizes(y_d)-1
578            delta_y = y_d(1)-y_d(0)
579            break
580         end if
581      end do
582      do varn=0,dim-1
583         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
584            z_d     = f_att->$vNam(varn)$
585            zdim    = dimsizes(z_d)-1
586            delta_z = z_d(1)-z_d(0)
587            break
588         end if
589      end do
590      do varn=0,dim-1
591         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
592            x_d     = f_att->$vNam(varn)$
593            xdim    = dimsizes(x_d)-1
594            delta_x = x_d(1)-x_d(0)
595            break
596         else
597            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
598               x_d     = f_att->$vNam(varn)$
599               xdim    = dimsizes(x_d)-1
600               delta_x = -1.d
601               break
602            else
603               xdim    = 0
604               delta_x = -1.d
605            end if
606         end if
607      end do
608   end if
609 
610   ; ****************************************************
611   ; set up ranges of x-, y- and z-coordinates
612   ; ****************************************************         
613                   
614   if (xs .EQ. -1.d) then               
615      xs = x_d(0)
616      if (delta_x .EQ. -1) then
617         print(" ")
618         print("You cannot choose a start value for x, there are preseted layers for x")
619         print(" ")
620         xstart=0
621      end if
622   else
623      if (delta_x .EQ. -1) then
624         print(" ")
625         print("You cannot choose a start value for x, there are preseted layers for x")
626         print(" ")
627         xstart=0
628      else
629         if (xs .LT. 0-delta_x/2) then
630            print(" ")
631            print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m")
632            print(" ")
633            exit
634         end if
635         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
636            if (xs .GE. xdim*delta_x) then
637               print(" ")
638               print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m")
639               print(" ")
640               exit
641            end if
642         else
643            if (xs .GT. xdim*delta_x) then
644               print(" ")
645               print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m")
646               print(" ")
647               exit
648            end if
649         end if
650      end if
651   end if
652
653   if (ys .EQ. -1.d) then   
654      ys = y_d(0)
655      if (delta_y .EQ. -1) then
656         print(" ")
657         print("You cannot choose a start value for y, there are preseted layers for y")
658         print(" ")
659         ystart=0
660      end if
661   else
662      if (delta_y .EQ. -1) then
663         print(" ")
664         print("You cannot choose a start value for y, there are preseted layers for y")
665         print(" ")
666         ystart=0
667      else
668         if (ys .LT. 0-delta_y/2) then
669            print(" ")
670            print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m")
671            print(" ")
672            exit
673         end if
674         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
675            if (ys .GE. ydim*delta_y) then
676               print(" ")
677               print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m")
678               print(" ")
679               exit
680            end if
681         else
682            if (ys .GT. ydim*delta_y) then
683               print(" ")
684               print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
685               print(" ")
686               exit
687            end if
688         end if
689      end if
690   end if
691 
692   if (zs .EQ. -1) then                         
693      zs = 0
694      if (delta_z .EQ. -1) then
695         print(" ")
696         print("You cannot choose a start value for z, there are preseted layers for z")
697         print(" ")
698      end if
699   else
700      if (delta_z .EQ. -1) then
701         print(" ")
702         print("You cannot choose a start value for z, there are preseted layers for z")
703         print(" ")
704         zs = 0
705      else
706         if (zs .LT. 0) then
707            print(" ")
708            print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0")
709            print(" ")
710            exit
711         end if
712         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
713            if (zs .GE. zdim) then
714               print(" ")
715               print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim)
716               print(" ")
717               exit
718            end if
719         else
720            if (zs .GT. zdim) then
721               print(" ")
722               print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim)
723               print(" ")
724               exit
725            end if
726         end if
727      end if
728   end if 
729
730   if (xe .EQ. -1) then   
731      xe = x_d(xdim)
732      if (delta_x .EQ. -1) then
733         print(" ")
734         print("You cannot choose an end value for x, there are preseted layers for x")
735         print(" ")
736         xend=xdim
737      end if
738   else
739      if (delta_x .EQ. -1) then
740         print(" ")
741         print("You cannot choose an end value for x, there are preseted layers for x")
742         print(" ")
743         xend=xdim
744      else
745         if (xe .GT. xdim*delta_x) then
746            print(" ")
747            print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m")
748            print(" ")
749            exit
750         end if
751         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
752            if (xe .LE. 0-delta_x/2)
753               print(" ")
754               print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
755               print(" ")
756               exit
757            end if
758            if (xe .LE. xs) then
759               print(" ")
760               print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m")
761               print(" ")
762               exit
763            end if
764            if (xe .EQ. xs+1) then
765               print(" ")
766               print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m")
767               print(" ")
768               exit
769            end if
770         else
771            if (xe .LT. 0-delta_x/2)
772               print(" ")
773               print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
774               print(" ")
775               exit
776            end if
777            if (xe .LT. xs) then
778               print(" ")
779               print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
780               print(" ")
781               exit
782            end if
783         end if
784      end if               
785   end if
786   
787   if (ye .EQ. -1) then 
788      ye = y_d(ydim)
789      if (delta_y .EQ. -1) then
790         print(" ")
791         print("You cannot choose an end value for y, there are preseted layers for y")
792         print(" ")
793         yend=ydim
794      end if
795   else
796      if (delta_y .EQ. -1) then
797         print(" ")
798         print("You cannot choose an end value for y, there are preseted layers for y")
799         print(" ")
800         yend=ydim
801      else
802         if (ye .GT. ydim*delta_y) then
803            print(" ")
804            print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m")
805            print(" ")
806            exit
807         end if
808         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
809            if (ye .LE. 0-delta_y/2)
810               print(" ")
811               print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
812               print(" ")
813               exit
814            end if
815            if (ye .LE. ys) then
816               print(" ")
817               print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m")
818               print(" ")
819               exit
820            end if
821            if (ye .EQ. ys+1) then
822               print(" ")
823               print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m")
824               print(" ")
825               exit
826            end if
827         else
828            if (ye .LT. 0-delta_y/2)
829               print(" ")
830               print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m")
831               print(" ")
832               exit
833            end if
834            if (ye .LT. ys) then
835               print(" ")
836               print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m")
837               print(" ")
838               exit
839            end if
840         end if
841      end if
842   end if
843 
844   if (ze .EQ. -1) then 
845      ze = zdim
846      if (delta_z .EQ. -1) then
847         print(" ")
848         print("You cannot choose an end value for z, there are preseted layers for z")
849         print(" ")
850         ze = zdim
851      end if
852   else
853      if (delta_z .EQ. -1) then
854         print(" ")
855         print("You cannot choose an end value for z, there are preseted layers for z")
856         print(" ")
857         ze = zdim
858      else
859         if (ze .GT. zdim) then
860            print(" ")
861            print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim)
862            print(" ")
863            exit
864         end if
865         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
866            if (ze .LE. 0)
867               print(" ")
868               print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0")
869               print(" ")
870               exit
871            end if
872            if (ze .LE. zs) then
873               print(" ")
874               print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs)
875               print(" ")
876               exit
877            end if
878            if (ze .EQ. zs+1) then
879               print(" ")
880               print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs)
881               print(" ")
882               exit
883            end if
884         else
885            if (ze .LT. 0)
886               print(" ")
887               print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0")
888               print(" ")
889               exit
890            end if
891            if (ze .LT. zs) then
892               print(" ")
893               print("range end for z-coordinate = "+ze+" is lower than start range = "+zs)
894               print(" ")
895               exit
896            end if
897         end if
898      end if
899   end if
900
901   if (delta_x .NE. -1) then
902      do i=0,xdim   
903         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
904            xstart=i
905            break
906         end if
907      end do
908      do i=0,xdim   
909         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
910            xend=i
911            break
912         end if
913      end do
914   end if
915   if (delta_y .NE. -1) then
916      do i=0,ydim   
917         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
918            ystart=i
919            break
920         end if
921      end do
922      do i=0,ydim   
923         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
924            yend=i
925            break
926         end if
927      end do
928   end if
929
930   if( shape .eq. 1 ) then
931      if (xyc .EQ. 1)then
932         cs_res@vpWidthF    = (xe-xs)/(ye-ys)       
933         cs_res@vpHeightF   = 1
934         vecres@vpWidthF    = (xe-xs)/(ye-ys)       
935         vecres@vpHeightF   = 1
936         if (xe-xs .GT. ye-ys)then
937            if (no_rows .LE. no_columns)then
938               cs_res@gsnPaperOrientation     = "landscape"
939               vecres@gsnPaperOrientation     = "landscape"
940               cs_res@lbOrientation = "Horizontal"
941            else
942               cs_res@gsnPaperOrientation     = "portrait"
943               vecres@gsnPaperOrientation     = "portrait"
944               cs_res@lbOrientation = "Vertical"
945            end if 
946         else
947            if (no_rows .GE. no_columns)then
948               cs_res@gsnPaperOrientation     = "portrait"
949               vecres@gsnPaperOrientation     = "portrait"
950               cs_res@lbOrientation = "Vertical"
951            else
952               cs_res@gsnPaperOrientation     = "landscape"
953               vecres@gsnPaperOrientation     = "landscape"
954               cs_res@lbOrientation = "Horizontal"
955            end if
956         end if
957      end if
958      if (xzc .EQ. 1)then
959         cs_res@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
960         cs_res@vpHeightF   = 1
961         vecres@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
962         vecres@vpHeightF   = 1
963         if (xe-xs .GT. (delta_x*(ze-zs)))then
964            if (no_rows .LE. no_columns)then
965               cs_res@gsnPaperOrientation     = "landscape"
966               vecres@gsnPaperOrientation     = "landscape"
967               cs_res@lbOrientation = "Horizontal"
968            else
969               cs_res@gsnPaperOrientation     = "portrait"
970               vecres@gsnPaperOrientation     = "portrait"
971               cs_res@lbOrientation = "Vertical"
972            end if 
973         else
974            if (no_rows .GE. no_columns)then
975               cs_res@gsnPaperOrientation     = "portrait"
976               vecres@gsnPaperOrientation     = "portrait"
977               cs_res@lbOrientation = "Vertical"
978            else
979               cs_res@gsnPaperOrientation     = "landscape"
980               vecres@gsnPaperOrientation     = "landscape"
981               cs_res@lbOrientation = "Horizontal"
982            end if
983         end if
984      end if
985      if (yzc .EQ. 1)then
986         cs_res@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
987         cs_res@vpHeightF   = 1
988         vecres@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
989         vecres@vpHeightF   = 1
990         if (ye-ys .GT. (delta_y*(ze-zs)))then
991            if (no_rows .LE. no_columns)then
992               cs_res@gsnPaperOrientation     = "landscape"
993               vecres@gsnPaperOrientation     = "landscape"
994               cs_res@lbOrientation = "Horizontal"
995            else
996               cs_res@gsnPaperOrientation     = "portrait"
997               vecres@gsnPaperOrientation     = "portrait"
998               cs_res@lbOrientation = "Vertical"
999            end if 
1000         else
1001            if (no_rows .GE. no_columns)then
1002               cs_res@gsnPaperOrientation     = "portrait"
1003               vecres@gsnPaperOrientation     = "portrait"
1004               cs_res@lbOrientation = "Vertical"
1005            else
1006               cs_res@gsnPaperOrientation     = "landscape"
1007               vecres@gsnPaperOrientation     = "landscape"
1008               cs_res@lbOrientation = "Horizontal"
1009            end if
1010         end if
1011      end if
1012   end if
1013
1014   delete(xs)
1015   delete(xe)   
1016   delete(ys)
1017   delete(ye)
1018
1019   xs=xstart
1020   xe=xend
1021   ys=ystart
1022   ye=yend
1023
1024   if (xyc .EQ. 1) then
1025      d= (ye-ys+1)/(major_ticks_y-1)
1026      e=(xe-xs+1)/(major_ticks_x-1)
1027      array_yl       =new(major_ticks_y,integer)
1028      array_yl_labels=new(major_ticks_y,double)
1029      array_minor_yl =new((major_ticks_y-1)*4,double)
1030      array_xb       =new(major_ticks_x,integer)
1031      array_xb_labels=new(major_ticks_x,double)
1032      array_minor_xb =new((major_ticks_x-1)*4,double)
1033      array_yl(0)=ys
1034      array_xb(0)=xs
1035      array_yl_labels(0)=y_d(ys)
1036      array_xb_labels(0)=x_d(xs)
1037      do ar=1,major_ticks_y-1
1038         array_yl(ar)=d*(ar-1)+d-1
1039         array_yl_labels(ar) = y_d(array_yl(ar))
1040         if (ar .GT. 0)
1041            do min_ar=0,3
1042                array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+(y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1)
1043            end do
1044         end if 
1045      end do
1046      do br=1,major_ticks_x-1
1047         array_xb(br)=e*(br-1)+e-1
1048         array_xb_labels(br) = x_d(array_xb(br))
1049         if (br .GT. 0)
1050            do min_br=0,3
1051               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1052            end do
1053         end if
1054      end do
1055   end if
1056 
1057   if (xzc .EQ. 1) then
1058      d=(ze-zs+1)/(major_ticks_y-1)
1059      e=(xe-xs+1)/(major_ticks_x-1)
1060      array_yl       =new(major_ticks_y,integer)
1061      array_yl_labels=new(major_ticks_y,double)
1062      array_minor_yl =new((major_ticks_y-1)*4,double)
1063      array_xb       =new(major_ticks_x,integer)
1064      array_xb_labels=new(major_ticks_x,double)
1065      array_minor_xb =new((major_ticks_x-1)*4,double)
1066      array_yl(0)=zs
1067      array_xb(0)=xs
1068      array_yl_labels(0)=z_d(zs)
1069      array_xb_labels(0)=x_d(xs)
1070      do ar=1,major_ticks_y-1
1071         array_yl(ar)=d*(ar-1)+d-1
1072         array_yl_labels(ar) = z_d(array_yl(ar))
1073         if (ar .GT. 0)
1074            do min_ar=0,3
1075               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1076            end do
1077         end if 
1078      end do
1079      do br=1,major_ticks_x-1
1080         array_xb(br)=e*(br-1)+e-1
1081         array_xb_labels(br) = x_d(array_xb(br))
1082         if (br .GT. 0)
1083            do min_br=0,3
1084               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1085            end do
1086         end if
1087      end do
1088   end if
1089
1090   if (yzc .EQ. 1) then
1091      d=(ze-zs+1)/(major_ticks_y-1)
1092      e=(ye-ys+1)/(major_ticks_x-1)
1093      array_yl       =new(major_ticks_y,integer)
1094      array_yl_labels=new(major_ticks_y,double)
1095      array_minor_yl =new((major_ticks_y-1)*4,double)
1096      array_xb       =new(major_ticks_x,integer)
1097      array_xb_labels=new(major_ticks_x,double)
1098      array_minor_xb =new((major_ticks_x-1)*4,double)
1099      array_yl(0)=zs
1100      array_xb(0)=ys
1101      array_yl_labels(0)=z_d(zs)
1102      array_xb_labels(0)=y_d(ys)
1103      do ar=1,major_ticks_y-1
1104         array_yl(ar)=d*(ar-1)+d-1
1105         array_yl_labels(ar) = z_d(array_yl(ar))
1106         if (ar .GT. 0)
1107            do min_ar=0,3
1108               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1109            end do
1110         end if 
1111      end do
1112      do br=1,major_ticks_x-1
1113         array_xb(br)=e*(br-1)+e-1
1114         array_xb_labels(br) = y_d(array_xb(br))
1115         if (br .GT. 0)
1116            do min_br=0,3
1117               array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+(y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1)
1118            end do
1119         end if
1120      end do
1121   end if
1122 
1123   if (axes_explicit .EQ. 1)then
1124      cs_res@tmYLMode = "Explicit"
1125      cs_res@tmXBMode = "Explicit"
1126      if (xyc .EQ. 1)then
1127         cs_res@tmYLValues = y_d(array_yl)
1128         cs_res@tmXBValues = x_d(array_xb)     
1129         cs_res@tmYLLabels = array_yl_labels/norm_y
1130         cs_res@tmXBLabels = array_xb_labels/norm_x
1131         if (norm_x .NE. 1.)then
1132            cs_res@tiXAxisString = "x ["+unit_x+"]"
1133         else
1134            cs_res@tiXAxisString = "x [m]"
1135         end if
1136         if (norm_y .NE. 1.)then
1137            cs_res@tiYAxisString = "y ["+unit_y+"]"
1138         else
1139            cs_res@tiYAxisString = "y [m]"
1140         end if   
1141      end if
1142      if (xzc .EQ. 1)then
1143         cs_res@tmYLValues = z_d(array_yl)
1144         cs_res@tmXBValues = x_d(array_xb) 
1145         cs_res@tmYLLabels = array_yl_labels/norm_z
1146         cs_res@tmXBLabels = array_xb_labels/norm_x
1147         if (norm_x .NE. 1.)then
1148            cs_res@tiXAxisString = "x ["+unit_x+"]"
1149         else
1150            cs_res@tiXAxisString = "x [m]"
1151         end if
1152         if (norm_z .NE. 1.)then
1153            cs_res@tiYAxisString = "z ["+unit_z+"]"
1154         else
1155            cs_res@tiYAxisString = "z [m]"
1156         end if
1157      end if
1158      if (yzc .EQ. 1)then
1159         cs_res@tmYLValues = z_d(array_yl)
1160         cs_res@tmXBValues = y_d(array_xb) 
1161         cs_res@tmYLLabels = array_yl_labels/norm_z
1162         cs_res@tmXBLabels = array_xb_labels/norm_y
1163         if (norm_y .NE. 1.)then
1164            cs_res@tiXAxisString = "y ["+unit_y+"]"
1165         else
1166            cs_res@tiXAxisString = "y [m]"
1167         end if
1168         if (norm_z .NE. 1.)then
1169            cs_res@tiYAxisString = "z ["+unit_z+"]"
1170         else
1171            cs_res@tiYAxisString = "z [m]"
1172         end if
1173      end if
1174      cs_res@tmYLMinorValues = array_minor_yl
1175      cs_res@tmXBMinorValues = array_minor_xb
1176   else
1177      if (axes_explicit .EQ. 0)then 
1178         cs_res@tmYLMinorPerMajor = 4
1179         cs_res@tmXBMinorPerMajor = 4
1180      else
1181         print(" ")
1182         print("'axes_explicit' is invalid and set to 0")
1183         print(" ")
1184         axes_explicit = 0
1185         cs_res@tmYLMinorPerMajor = 4
1186         cs_res@tmXBMinorPerMajor = 4
1187      end if
1188      if (xyc .EQ. 1)then
1189         cs_res@tiXAxisString = "x [m]"
1190         cs_res@tiYAxisString = "y [m]"
1191      end if
1192      if (xzc .EQ. 1)then
1193         cs_res@tiXAxisString = "x [m]"
1194         cs_res@tiYAxisString = "z [m]"
1195      end if
1196      if (yzc .EQ. 1)then
1197         cs_res@tiXAxisString = "y [m]"
1198         cs_res@tiYAxisString = "z [m]"
1199      end if
1200   end if
1201
1202   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1203      if (xe .EQ. xs+1) then
1204         print(" ")
1205         print("range end for x-coordinate="+xe*delta_x+"m) must be at least two")
1206         print("more gridpoints(="+2*delta_x+"m) greater than start range="+xs*delta_x+"m)")
1207         print(" ")
1208         exit
1209      end if
1210   end if
1211   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1212      if (ye .EQ. ys+1) then
1213         print(" ")
1214         print("range end for y-coordinate="+ye*delta_y+"m) must be at least two")
1215         print("more gridpoints(="+2*delta_y+"m) greater than start range="+ys*delta_y+"m)")
1216         print(" ")
1217         exit
1218      end if
1219   end if
1220   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1221      if (ze .EQ. zs+1) then
1222         print(" ")
1223         print("range end for x-coordinate="+ze+") must be at least two")
1224         print("more gridpoints greater than start range="+zs+")")
1225         print(" ")
1226         exit
1227      end if
1228   end if
1229 
1230   if (xyc .EQ. 1) then
1231      no_layer = (ze-zs)+1
1232      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1233   end if
1234   if (xzc .EQ. 1) then
1235      no_layer = (ye-ys)+1
1236      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1237   end if
1238   if (yzc .EQ. 1) then
1239      no_layer = (xe-xs)+1
1240      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1241   end if
1242
1243   MinVal = new(dim,float)
1244   MaxVal = new(dim,float)
1245   unit   = new(dim,string)
1246
1247   ; ****************************************************
1248   ; define inner and outer loops depending on "sort"
1249   ; ****************************************************       
1250   
1251   if ( xyc .eq. 1 ) then
1252      lays = zs
1253      laye = ze
1254   end if
1255   if ( xzc .eq. 1 ) then
1256      lays = ys
1257      laye = ye
1258   end if
1259   if ( yzc .eq. 1) then
1260      lays = xs
1261      laye = xe
1262   end if
1263
1264   if ( sort .eq. "time" ) then
1265      lis = start_time_step
1266      lie = end_time_step
1267      los = lays
1268      loe = laye
1269   else
1270      lis = lays
1271      lie = laye
1272      los = start_time_step
1273      loe = end_time_step
1274   end if
1275
1276   check = True
1277   v1=0
1278   v2=0
1279   no_var=0
1280   n=0
1281   no_zu1=0
1282
1283   do varn=dim-1,0,1       
1284   
1285      if (vector .EQ. 1) then   
1286         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1287         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1288      end if
1289           
1290      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
1291         check = False
1292      else
1293         check = True
1294      end if 
1295
1296      if (var .NE. "all") then
1297         check = isStrSubset( var,","+vNam(varn)+"," )
1298      end if
1299
1300      if(check) then
1301     
1302         no_var=no_var+1
1303
1304         if (vector .EQ. 1) then
1305            if (","+vNam(varn)+"," .EQ. vec1) then
1306               v1=v1+1
1307            end if
1308            if (","+vNam(varn)+"," .EQ. vec2) then
1309               v2=v2+1
1310            end if
1311         end if
1312
1313         if (xyc .EQ. 1) then
1314            temp = f[:]->$vNam(varn)$
1315            data_att = f_att->$vNam(varn)$
1316            if(vNam(varn) .eq. "ts_xy" .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy")
1317              ;these variables depend zu1_xy and that's why they have only one z-layer
1318              data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1319              no_zu1=no_zu1+1
1320            else
1321              data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1322            end if
1323            delete(temp)               
1324         end if
1325         if ( xzc .eq. 1 ) then
1326            temp = f[:]->$vNam(varn)$
1327            data_att = f_att->$vNam(varn)$
1328            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1329            delete(temp)
1330         end if
1331         if ( yzc .eq. 1) then
1332            temp = f[:]->$vNam(varn)$
1333            data_att = f_att->$vNam(varn)$
1334            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1335            delete(temp)
1336         end if
1337         if (check_vec1) then
1338            vect1=data(varn,:,:,:,:)
1339         end if
1340         if (check_vec2) then
1341            vect2=data(varn,:,:,:,:)
1342         end if
1343
1344         data!0 = "var"
1345         data!1 = "t"
1346         data!2 = "z"
1347         data!3 = "y"
1348         data!4 = "x" 
1349
1350         MinVal(varn) = min(data(varn,:,:,:,:))
1351         MaxVal(varn) = max(data(varn,:,:,:,:))
1352         
1353         unit(varn) = data_att@units
1354         delete(data_att)
1355
1356      end if
1357
1358   end do
1359
1360   if (no_var .EQ. 0) then
1361      print(" ")
1362      print("The variables var='"+var+"' do not exist on your input file;")
1363      print("be sure to have one comma before and after each variable")
1364      print(" ")
1365      exit
1366   end if
1367
1368   var_input=new(no_var,string)
1369   numb=0
1370   no_var=0
1371   
1372   do varn=dim-1,0,1   
1373   
1374      if (vector .EQ. 1) then   
1375         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1376         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1377      end if
1378           
1379      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
1380         check = False
1381      else
1382         check = True
1383      end if 
1384
1385      if (var .NE. "all") then
1386         check = isStrSubset( var,","+vNam(varn)+"," )
1387      end if
1388
1389      if(check) then     
1390         var_input(numb)=vNam(varn)
1391         numb=numb+1     
1392      end if
1393     
1394      if (check) then
1395         no_var = no_var+1
1396      end if
1397     
1398   end do
1399
1400   if (no_var .EQ. 0) then
1401      print(" ")
1402      print("The variables var='"+var+"' do not exist on your input file;")
1403      print("be sure to have one comma before and after each variable")
1404      print(" ")
1405      exit
1406   end if
1407
1408   if (vector .EQ. 1) then
1409      if (v1 .EQ. 0)then
1410         print(" ")
1411         print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
1412         print("- "+var_input)
1413         print("be sure to have one comma before and after the variable")
1414         print(" ")
1415         exit
1416      end if
1417
1418      if (v2 .EQ. 0)then
1419         print(" ")
1420         print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
1421         print("- "+var_input)
1422         print("be sure to have one comma before and after the variable")
1423         print(" ")
1424         exit
1425      end if
1426   end if
1427
1428   ; ***************************************************
1429   ; open workstation(s)
1430   ; ***************************************************
1431
1432   wks_ps  = gsn_open_wks(format_out,file_out)
1433   gsn_define_colormap(wks_ps,"rainbow+white")
1434
1435   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1436      plot=new((/no_time*no_layer*no_var - no_time*(no_layer-1)*no_zu1 \
1437                                         + no_time*no_layer/),graphic)
1438   else
1439      plot=new((/no_time*no_layer*no_var - no_time*(no_layer-1)*no_zu1/),graphic)
1440   end if
1441   dim_plot=dimsizes(plot)
1442
1443   page = 0
1444   n=0
1445   print(" ")
1446   print("Plots sorted by '"+sort+"'")
1447   print(" ")
1448
1449   ; ***************************************************
1450   ; create plots
1451   ; ***************************************************
1452
1453
1454   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1455      do lo = los, loe                                         
1456         do li = lis, lie
1457            if (xyc .EQ. 1)then
1458               if (sort .EQ. "time")then
1459                  if(z_d(lo) .eq. -1.d) then
1460                    level = "z-average"
1461                  else
1462                    level = "z=" + z_d(lo) + "m"
1463                  end if
1464               else
1465                  if(z_d(li) .eq. -1.d) then
1466                    level = "z-average"
1467                  else
1468                    level = "z=" + z_d(li) + "m"
1469                  end if
1470               end if
1471            end if
1472            if (xzc .EQ. 1)then
1473               if (sort .EQ. "time")then
1474                 if(y_d(lo) .eq. -1.d) then
1475                   level = "y-average"
1476                 else
1477                   level = "y=" + y_d(lo) + "m"
1478                 end if
1479               else
1480                  if(y_d(li) .eq. -1.d) then
1481                    level = "y-average"
1482                  else
1483                    level = "y=" + y_d(li) + "m"
1484                  end if
1485               end if
1486            end if
1487            if (yzc .EQ. 1)then
1488               if (sort .EQ. "time")then
1489                 if(x_d(lo) .eq. -1.d) then
1490                    level = "x-average"
1491                 else
1492                    level = "x=" + x_d(lo) + "m"
1493                 end if
1494               else
1495                 if(x_d(li) .eq. -1.d) then
1496                    level = "x-average"
1497                 else
1498                    level = "x=" + x_d(li) + "m"
1499                 end if
1500               end if
1501            end if               
1502            vecres                  = True            ; vector only resources
1503            vecres@gsnDraw          = False           ; don't draw
1504            vecres@gsnFrame         = False           ; don't advance frame
1505            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1506            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1507            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1508            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1509            vecres@tmXBLabelFontHeightF   = font_size
1510            vecres@tmYLLabelFontHeightF   = font_size
1511            vecres@tiXAxisFontHeightF     = font_size
1512            vecres@tiYAxisFontHeightF     = font_size
1513            if (sort .EQ. "time")then
1514               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1515            else
1516               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1517            end if
1518            vecres@tiXAxisString    = " "
1519            if (xyc .EQ. 1)then 
1520               if (sort .EQ. "time")then                                         
1521                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1522               else
1523                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),vect2(lo,li-lis,:,:),vecres)
1524               end if
1525            end if
1526            if (xzc .EQ. 1) then
1527               if (sort .EQ. "time")then
1528                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),vect2(li,:,lo-los,:),vecres)
1529               else
1530                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),vect2(lo,:,li-lis,:),vecres)
1531               end if
1532            end if
1533            if (yzc .EQ. 1) then
1534               if (sort .EQ. "time")then
1535                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),vect2(li,:,:,lo-los),vecres)
1536               else
1537                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),vect2(lo,:,:,li-lis),vecres)
1538               end if
1539            end if
1540            n=n+1
1541         end do
1542      end do
1543   end if
1544 
1545   do varn=dim-1,0,1
1546
1547      if (vector .EQ. 1) then   
1548         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1549      end if
1550     
1551      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
1552         check = False
1553      else
1554         check = True
1555      end if   
1556 
1557      if (var .NE. "all") then
1558         check = isStrSubset( var,","+vNam(varn)+"," )
1559      end if
1560   
1561      if(check) then
1562
1563         space=(MaxVal(varn)-MinVal(varn))/24
1564 
1565         cs_res@cnMinLevelValF = MinVal(varn)
1566         cs_res@cnMaxLevelValF = MaxVal(varn)
1567
1568         cs_res@cnLevelSpacingF = space
1569     
1570         ; ****************************************************
1571         ; loops over time and layer
1572         ; ****************************************************
1573
1574         
1575         do lo = los, loe                                       
1576            do li = lis, lie
1577
1578               ; ****************************************************
1579               ; xy cross section
1580               ; ****************************************************
1581
1582               if ( xyc .eq. 1 ) then
1583         
1584                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1585                  cs_res@gsnRightString = unit(varn)
1586                 
1587                  if ( sort .eq. "time" ) then
1588                     if ( z_d(lo) .eq. -1)then
1589                        if (delta_z .EQ. -1) then
1590                           level = "-average"
1591                        else
1592                           level = "=" + z_d(lo) + "m"
1593                        end if
1594                     else
1595                        level = "=" + z_d(lo) + "m"
1596                     end if
1597
1598                     if(vNam(varn) .eq. "ts_xy" .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy")then
1599                         loe = 0
1600                         level = "=" + zu1(0) + "m"
1601                     end if
1602                   
1603                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level                               
1604                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo-los,:,:),cs_res)
1605                     if (vector .EQ. 1 .AND. check_vecp) then
1606                        vecres                  = True            ; vector only resources
1607                        vecres@gsnDraw          = False           ; don't draw
1608                        vecres@gsnFrame         = False           ; don't advance frame
1609                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1610                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1611                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
1612                        vecres@gsnRightString   = " "
1613                        vecres@gsnLeftString    = " "
1614                        vecres@tiXAxisString    = " "   
1615                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1616                        overlay(plot(n), plot_vec)
1617                     end if                         
1618                  end if
1619                 
1620                  if ( sort .eq. "layer" ) then
1621                     if (z_d(zs+li) .eq. -1) then
1622                        if (delta_z .EQ. -1) then
1623                           level = "-average"
1624                        else
1625                           level = "=" + z_d(li) + "m"
1626                        end if
1627                     else
1628                        level = "=" + z_d(li) + "m"
1629                     end if
1630       
1631                     if(vNam(varn) .eq. "ts_xy" .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy")then
1632                         lie = 0
1633                         level = "=" + zu1(0) + "m"
1634                     end if
1635               
1636                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level
1637                     
1638                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li-lis,:,:),cs_res)
1639                     if (vector .EQ. 1 .AND. check_vecp) then
1640                        vecres                  = True            ; vector only resources
1641                        vecres@gsnDraw          = False           ; don't draw
1642                        vecres@gsnFrame         = False           ; don't advance frame
1643                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1644                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1645                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1646                        vecres@gsnRightString   = " "             ; turn off right string
1647                        vecres@gsnLeftString    = " "             ; turn off left string
1648                        vecres@tiXAxisString    = " "   
1649                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1650                        overlay(plot(n), plot_vec)
1651                     end if
1652                  end if
1653               end if
1654
1655               ; ****************************************************
1656               ; xz cross section
1657               ; ****************************************************
1658
1659               if ( xzc .eq. 1 ) then
1660                 
1661                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1662               
1663                  if ( sort .eq. "time" ) then
1664                     if ( y_d(lo) .eq. -1 ) then
1665                        level = "-average"
1666                     else
1667                        level = "=" + y_d(lo) + "m"
1668                     end if
1669                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
1670                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo-los,:),cs_res)
1671                     if (vector .EQ. 1 .AND. check_vecp) then
1672                        vecres                  = True            ; vector only resources
1673                        vecres@gsnDraw          = False           ; don't draw
1674                        vecres@gsnFrame         = False           ; don't advance frame
1675                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1676                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1677                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1678                        vecres@gsnRightString   = " "             ; turn off right string
1679                        vecres@gsnLeftString    = " "             ; turn off left string
1680                        vecres@tiXAxisString    = " "   
1681                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),vect2(li,:,lo-los,:),vecres)
1682                        overlay(plot(n), plot_vec)
1683                     end if
1684                  end if
1685
1686                  if ( sort .eq. "layer" ) then
1687                     if ( y_d(li) .eq. -1 ) then
1688                        level = "-average"
1689                     else
1690                        level = "=" + y_d(li) + "m"
1691                     end if
1692                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
1693                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li-lis,:),cs_res)
1694                     if (vector .EQ. 1 .AND. check_vecp) then
1695                        vecres                  = True            ; vector only resources
1696                        vecres@gsnDraw          = False           ; don't draw
1697                        vecres@gsnFrame         = False           ; don't advance frame
1698                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1699                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1700                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1701                        vecres@gsnRightString   = " "             ; turn off right string
1702                        vecres@gsnLeftString    = " "             ; turn off left string
1703                        vecres@tiXAxisString    = " "   
1704                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),vect2(lo,:,li-lis,:),vecres)
1705                        overlay(plot(n), plot_vec)
1706                     end if
1707                  end if                 
1708               end if
1709
1710               ; ****************************************************
1711               ; yz cross section
1712               ; ****************************************************
1713
1714               if ( yzc .eq. 1 ) then
1715                                 
1716                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1717
1718                  if ( sort .eq. "time" ) then
1719                     if ( x_d(xs+lo) .eq. -1 ) then
1720                        level = "-average"
1721                     else
1722                        level = "=" + x_d(lo) + "m"
1723                     end if
1724                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
1725                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo-los),cs_res)
1726                     if (vector .EQ. 1 .AND. check_vecp) then
1727                        vecres                  = True            ; vector only resources
1728                        vecres@gsnDraw          = False           ; don't draw
1729                        vecres@gsnFrame         = False           ; don't advance frame
1730                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1731                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1732                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1733                        vecres@gsnRightString   = " "             ; turn off right string
1734                        vecres@gsnLeftString    = " "             ; turn off left string
1735                        vecres@tiXAxisString    = " "   
1736                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),vect2(li,:,:,lo-los),vecres)
1737                        overlay(plot(n), plot_vec)
1738                     end if
1739                  end if
1740
1741                  if ( sort .eq. "layer" ) then
1742                     if ( x_d(xs+li) .eq. -1 ) then
1743                        level = "-average"
1744                     else
1745                        level = "=" + x_d(li) + "m"
1746                     end if
1747                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
1748                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li-lis),cs_res)
1749                     if (vector .EQ. 1 .AND. check_vecp)then
1750                        vecres                  = True            ; vector only resources
1751                        vecres@gsnDraw          = False           ; don't draw
1752                        vecres@gsnFrame         = False           ; don't advance frame
1753                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1754                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1755                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1756                        vecres@gsnRightString   = " "             ; turn off right string
1757                        vecres@gsnLeftString    = " "             ; turn off left string
1758                        vecres@tiXAxisString    = " "   
1759                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),vect2(lo,:,:,li-lis),vecres)
1760                        overlay(plot(n), plot_vec)
1761                     end if
1762                  end if
1763               end if         
1764               n=n+1 
1765            end do     
1766         end do 
1767      end if
1768   end do
1769
1770   ; ***************************************************
1771   ; merge plots onto one page
1772   ; ***************************************************   
1773   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
1774      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1775         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1776         print(" ")
1777         print("Outputs to .eps or .epsi have only one frame")
1778         print(" ")
1779      else
1780         do np = 0,dim_plot-1,no_rows*no_columns   
1781            if ( np + no_rows*no_columns .gt. dim_plot-1) then
1782               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP)
1783            else
1784               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1785            end if
1786         end do
1787      end if
1788   else       
1789      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1790         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
1791         print(" ")
1792         print("Outputs to .eps or .epsi have only one frame")
1793         print(" ")
1794      else
1795         do np = 0,dim_plot-1,no_rows*no_columns 
1796            if ( np + no_rows*no_columns .gt. dim_plot-1) then
1797               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP)
1798            else
1799               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1800            end if
1801         end do
1802      end if
1803   end if
1804
1805   print(" ")
1806   print("Output to: " + file_out +"."+ format_out)
1807   print(" ")   
1808 
1809end
Note: See TracBrowser for help on using the repository browser.