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

Last change on this file since 289 was 286, checked in by weinreis, 16 years ago

ncl-script cross_sections.ncl modified

File size: 59.5 KB
RevLine 
[154]1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
[157]2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
[154]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
[190]7;***************************************************
8; load ncl_preferences.ncl
9;***************************************************
[194]10
11function which_script()
12local script
13begin
14   script="cross_section"
15   return(script)
16end
17       
[190]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
[194]30
[154]31begin
32
33   ; ***************************************************
34   ; set default parameter values and strings if not assigned in prompt or parameter list
35   ; ***************************************************
36   
[190]37   if (file_1 .EQ. "File in") then
38      print(" ")
[218]39      print("Declare input file 'file_1=' in 'ncl_preferences.ncl' or prompt")
[190]40      print(" ")
41      exit
[175]42   else
43      file_in = file_1   
[154]44   end if
[175]45
[190]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"
[154]51   end if
[175]52
[190]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" 
[154]58   end if
[190]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(" ")
[154]64      mode = "Fill"
65   end if
[175]66
[190]67   if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND. fill_mode .NE. "CellFill") then
68      print(" ")
[218]69      print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
[190]70      print(" ")
[154]71      fill_mode = "AreaFill"
72   end if
[190]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(" ")
[154]78      shape = 1
[175]79   end if
[190]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 
[154]86   end if
[190]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 
[154]93   end if
[190]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 
[154]100   end if
[190]101   
[154]102   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
103      print(" ")
[218]104      print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]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(" ")
[218]111         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]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(" ")
[218]119         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]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(" ")
[218]127         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]128         print(" ")
129         exit
130      end if
131   end if
132
[190]133   if (vector .NE. 0 .AND. vector .NE. 1) then
134      print(" ")
[218]135      print("Set 'vector' to 0 or 1")
[190]136      print(" ")
137      exit
[218]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
[190]158 
[154]159   ; ***************************************************
160   ; open input file
161   ; ***************************************************
[218]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
[154]169
[218]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
[154]182
[218]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)
[154]233   print(" ")
[190]234   print("Variables in input file:")
235   print("- "+ vNam)
[154]236   print(" ")
237   dim   = dimsizes(vNam)
[162]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
[195]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
[162]265         else
[195]266            print(" ")
267            print("Your input file contains xy data")
268            print(" ")
269         end if
[162]270      end if
[195]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
[162]282         else
[195]283            print(" ")
284            print("Your input file contains xz data")
285            print(" ")
286         end if
[162]287      end if
[195]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
[162]299         else
[195]300            print(" ")
301            print("Your input file contains yz data")
302            print(" ")
303         end if
[162]304      end if
305   else
306      print(" ")
[218]307      print("Your input file contains 3d or other data")
[162]308      print(" ")
309   end if
310   
[161]311   if (dim .EQ. 0) then
312      print(" ")
[162]313      print("There is no data on file")
[161]314      print(" ")
[162]315      exit
[161]316   end if
317
[154]318   ; ***************************************************
319   ; set up recourses
320   ; ***************************************************
321
322   cs_res                          = True
[218]323   vecres                          = True
[154]324   cs_res@gsnYAxisIrregular2Linear = True 
325 
326   cs_res@gsnDraw                 = False
327   cs_res@gsnFrame                = False
328   cs_res@gsnMaximize             = True
[218]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
[154]334   cs_res@tmXBMode                ="Automatic"
335   cs_res@tmYLMode                ="Automatic"
[218]336 
[162]337   cs_res@cnLevelSelectionMode    = "ManualLevels"
[218]338   cs_res@lbLabelFontHeightF = font_size_legend
339   cs_res@lbLabelStride = legend_label_stride
340   
[162]341
[154]342   cs_resP = True
[218]343   cs_resP@txString               = f_att@title
[175]344   cs_resP@txFuncCode             = "~"                 
[218]345   cs_resP@txFontHeightF          = 0.015       
[154]346 
347   if ( mode .eq. "Fill" ) then
348      cs_res@cnFillOn             = True
349      cs_res@gsnSpreadColors      = True
[218]350      cs_res@cnFillMode           = fill_mode   
[154]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
[218]358      cs_res@cnFillMode           = fill_mode
[154]359      cs_res@cnLinesOn            = True
360      cs_res@cnLineLabelsOn       = True
361   end if
362
[157]363   ; *********************************************
364   ; set up of start_time_step and end_time_step
365   ; *********************************************
[154]366
[218]367   t_all = f[:]->time
[162]368   nt = dimsizes(t_all)
369   delta_t = t_all(nt-1)/nt
370
371   ; ****************************************************       
[157]372   ; start of time step and different types of mistakes that could be done
[162]373   ; ****************************************************
[154]374
[190]375   if (start_time_step .EQ. -1.) then           
376      start_time_step=t_all(0)/3600     
[157]377   else
[162]378      if (start_time_step .GT. t_all(nt-1)/3600)
[157]379         print(" ")
[162]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")
[157]381         print(" ")
[218]382         print("Select another 'start_time_step'")
[162]383         print(" ")
[157]384         exit
385      end if
[162]386      if (start_time_step .LT. t_all(0)/3600)
[157]387         print(" ")
[162]388         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]389         print(" ")
[218]390         print("Select another 'start_time_step'")
[162]391         print(" ")
[157]392         exit
393      end if
394   end if
[218]395   do i=0,nt-1   
[190]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
[162]397         st=i
398         break
399      end if
400   end do
[190]401   
402   if (.not. isvar("st"))then
403      print(" ")
404      print("'start_time_step' = "+ start_time_step +"h is invalid")
405      print(" ")
[218]406      print("Select another 'start_time_step'")
[190]407      print(" ")
408      exit
409   end if
[162]410       
[157]411   ; ****************************************************
412   ; end of time step and different types of mistakes that could be done
413   ; ****************************************************
414
[190]415   if (end_time_step .EQ. -1.) then             
416      end_time_step = t_all(nt-1)/3600 
417   else 
[162]418      if (end_time_step .LT. t_all(0)/3600)
[157]419         print(" ")
[162]420         print("'end_time_step' = "+end_time_step+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]421         print(" ")
[218]422         print("Select another 'end_time_step'")
[162]423         print(" ")
[157]424         exit
425      end if
[162]426      if (end_time_step .GT. t_all(nt-1)/3600)
[157]427         print(" ")
[162]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(" ")
[218]430         print("Select another 'end_time_step'") 
[162]431         print(" ")
[157]432         exit
433      end if
[162]434      if (end_time_step .LT. start_time_step/3600)
[157]435         print(" ")
[162]436         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
[157]437         print(" ")
[218]438         print("Select another 'start_time_step' or 'end_time_step'")
[162]439         print(" ")
[157]440         exit
441      end if
442   end if
[162]443   do i=0,nt-1     
[190]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
[162]445         et=i
446         break
447      end if
448   end do
[190]449   
450   if (.not. isvar("et"))then
451      print(" ")
452      print("'end_time_step' = "+ end_time_step +"h is invalid")
453      print(" ")
[218]454      print("Select another 'end_time_step'")
[190]455      print(" ")
456      exit
457   end if
[162]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)
[157]463
[162]464   print(" ")
[175]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)
[162]467   print(" ")
468 
[161]469   no_time=(end_time_step-start_time_step)+1
470
[157]471   ; ****************************************************
472   ; set up variables for vector plot if required
473   ; ****************************************************   
474
475   if (vector .EQ. 1) then
[190]476      if (vec1 .EQ. "component1") then
477         print(" ")
[218]478         print("Indicate Vector 1 ('vec1') for Vector-Plot or set 'vector' to 0")
[190]479         print(" ")
480         exit
[157]481      end if
[190]482      if (vec2 .EQ. "component2") then
483         print(" ")
[218]484         print("Indicate Vector 2 ('vec2') for Vector-Plot or set 'vector' to 0")
[190]485         print(" ")
486         exit
487      end if           
[157]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
[162]498   if (xyc .EQ. 1) then
499      do varn=0,dim-1
[190]500         if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then   
[218]501            x_d     = f_att->$vNam(varn)$
[162]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   
[218]509            y_d     = f_att->$vNam(varn)$
[162]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
[218]517            z_d     = f_att->$vNam(varn)$
[162]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
[218]523               z_d     = f_att->$vNam(varn)$
[162]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
[218]537            x_d     = f_att->$vNam(varn)$
[162]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
[218]545            z_d     = f_att->$vNam(varn)$
[162]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
[218]553            y_d     = f_att->$vNam(varn)$
[162]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
[218]559               y_d     = f_att->$vNam(varn)$
[162]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
[218]573            y_d     = f_att->$vNam(varn)$
[162]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
[218]581            z_d     = f_att->$vNam(varn)$
[162]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")
[218]589            x_d     = f_att->$vNam(varn)$
[162]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
[218]595               x_d     = f_att->$vNam(varn)$
[162]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
[195]606 
[162]607   ; ****************************************************
608   ; set up ranges of x-, y- and z-coordinates
609   ; ****************************************************         
610                   
[190]611   if (xs .EQ. -1.d) then               
612      xs = x_d(0)
[195]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
[190]619   else
[162]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
[161]625      else
[162]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
[154]647      end if
[162]648   end if
[154]649
[190]650   if (ys .EQ. -1.d) then   
651      ys = y_d(0)
[195]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
[162]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
[154]678         else
[162]679            if (ys .GT. ydim*delta_y) then
[157]680               print(" ")
[162]681               print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
[157]682               print(" ")
683               exit
684            end if
[162]685         end if
686      end if
687   end if
688 
[190]689   if (zs .EQ. -1) then                         
[162]690      zs = 0
[195]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
[162]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 
[157]725
[190]726   if (xe .EQ. -1) then   
727      xe = x_d(xdim)
[195]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
[162]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
[190]760            if (xe .EQ. xs+1) then
[162]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
[154]766         else
[162]767            if (xe .LT. 0-delta_x/2)
[154]768               print(" ")
[162]769               print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
[154]770               print(" ")
[162]771               exit
[154]772            end if
[162]773            if (xe .LT. xs) then
[154]774               print(" ")
[162]775               print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
[154]776               print(" ")
777               exit
778            end if
[162]779         end if
780      end if               
781   end if
[190]782   
783   if (ye .EQ. -1) then 
784      ye = y_d(ydim)
[195]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
[162]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
[154]823         else
[162]824            if (ye .LT. 0-delta_y/2)
[154]825               print(" ")
[162]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 
[190]840   if (ze .EQ. -1) then 
[162]841      ze = zdim
[195]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
[162]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
[157]896
[162]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
[190]912      do i=0,ydim   
[162]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
[190]918      do i=0,ydim   
[162]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
[286]925
[190]926   if( shape .eq. 1 ) then
[195]927      if (xyc .EQ. 1)then
928         cs_res@vpWidthF    = (xe-xs)/(ye-ys)       
929         cs_res@vpHeightF   = 1
[218]930         vecres@vpWidthF    = (xe-xs)/(ye-ys)       
931         vecres@vpHeightF   = 1
932         if (xe-xs .GT. ye-ys)then
[286]933            if (no_rows .LE. no_columns)then
934               cs_res@gsnPaperOrientation     = "landscape"
935               vecres@gsnPaperOrientation     = "landscape"
936               cs_res@lbOrientation = "Horizontal"
937            else
938               cs_res@gsnPaperOrientation     = "portrait"
939               vecres@gsnPaperOrientation     = "portrait"
940               cs_res@lbOrientation = "Vertical"
941            end if 
[218]942         else
[286]943            if (no_rows .GE. no_columns)then
944               cs_res@gsnPaperOrientation     = "portrait"
945               vecres@gsnPaperOrientation     = "portrait"
946               cs_res@lbOrientation = "Vertical"
947            else
948               cs_res@gsnPaperOrientation     = "landscape"
949               vecres@gsnPaperOrientation     = "landscape"
950               cs_res@lbOrientation = "Horizontal"
951            end if
[218]952         end if
[195]953      end if
954      if (xzc .EQ. 1)then
955         cs_res@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
956         cs_res@vpHeightF   = 1
[218]957         vecres@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
958         vecres@vpHeightF   = 1
959         if (xe-xs .GT. (delta_x*(ze-zs)))then
[286]960            if (no_rows .LE. no_columns)then
961               cs_res@gsnPaperOrientation     = "landscape"
962               vecres@gsnPaperOrientation     = "landscape"
963               cs_res@lbOrientation = "Horizontal"
964            else
965               cs_res@gsnPaperOrientation     = "portrait"
966               vecres@gsnPaperOrientation     = "portrait"
967               cs_res@lbOrientation = "Vertical"
968            end if 
[218]969         else
[286]970            if (no_rows .GE. no_columns)then
971               cs_res@gsnPaperOrientation     = "portrait"
972               vecres@gsnPaperOrientation     = "portrait"
973               cs_res@lbOrientation = "Vertical"
974            else
975               cs_res@gsnPaperOrientation     = "landscape"
976               vecres@gsnPaperOrientation     = "landscape"
977               cs_res@lbOrientation = "Horizontal"
978            end if
[218]979         end if
[195]980      end if
981      if (yzc .EQ. 1)then
982         cs_res@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
983         cs_res@vpHeightF   = 1
[218]984         vecres@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
985         vecres@vpHeightF   = 1
986         if (ye-ys .GT. (delta_y*(ze-zs)))then
[286]987            if (no_rows .LE. no_columns)then
988               cs_res@gsnPaperOrientation     = "landscape"
989               vecres@gsnPaperOrientation     = "landscape"
990               cs_res@lbOrientation = "Horizontal"
991            else
992               cs_res@gsnPaperOrientation     = "portrait"
993               vecres@gsnPaperOrientation     = "portrait"
994               cs_res@lbOrientation = "Vertical"
995            end if 
[218]996         else
[286]997            if (no_rows .GE. no_columns)then
998               cs_res@gsnPaperOrientation     = "portrait"
999               vecres@gsnPaperOrientation     = "portrait"
1000               cs_res@lbOrientation = "Vertical"
1001            else
1002               cs_res@gsnPaperOrientation     = "landscape"
1003               vecres@gsnPaperOrientation     = "landscape"
1004               cs_res@lbOrientation = "Horizontal"
1005            end if
[218]1006         end if
[195]1007      end if
[190]1008   end if
[218]1009
[162]1010   delete(xs)
1011   delete(xe)   
1012   delete(ys)
1013   delete(ye)
1014
1015   xs=xstart
1016   xe=xend
1017   ys=ystart
1018   ye=yend
1019
[218]1020   if (xyc .EQ. 1) then
1021      d=(ye-ys+1)/(major_ticks_y-1)
1022      e=(xe-xs+1)/(major_ticks_x-1)
1023      array_yl       =new(major_ticks_y,integer)
1024      array_yl_labels=new(major_ticks_y,double)
1025      array_minor_yl =new((major_ticks_y-1)*4,double)
1026      array_xb       =new(major_ticks_x,integer)
1027      array_xb_labels=new(major_ticks_x,double)
1028      array_minor_xb =new((major_ticks_x-1)*4,double)
1029      array_yl(0)=ys
1030      array_xb(0)=xs
1031      array_yl_labels(0)=y_d(ys)
1032      array_xb_labels(0)=x_d(xs)
1033      do ar=1,major_ticks_y-1
1034         array_yl(ar)=d*(ar-1)+d-1
1035         array_yl_labels(ar) = y_d(array_yl(ar))
1036         if (ar .GT. 0)
1037            do min_ar=0,3
1038               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)
1039            end do
1040         end if 
1041      end do
1042      do br=1,major_ticks_x-1
1043         array_xb(br)=e*(br-1)+e-1
1044         array_xb_labels(br) = x_d(array_xb(br))
1045         if (br .GT. 0)
1046            do min_br=0,3
1047               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)
1048            end do
1049         end if
1050      end do
1051   end if
1052 
1053   if (xzc .EQ. 1) then
1054      d=(ze-zs+1)/(major_ticks_y-1)
1055      e=(xe-xs+1)/(major_ticks_x-1)
1056      array_yl       =new(major_ticks_y,integer)
1057      array_yl_labels=new(major_ticks_y,double)
1058      array_minor_yl =new((major_ticks_y-1)*4,double)
1059      array_xb       =new(major_ticks_x,integer)
1060      array_xb_labels=new(major_ticks_x,double)
1061      array_minor_xb =new((major_ticks_x-1)*4,double)
1062      array_yl(0)=zs
1063      array_xb(0)=xs
1064      array_yl_labels(0)=z_d(zs)
1065      array_xb_labels(0)=x_d(xs)
1066      do ar=1,major_ticks_y-1
1067         array_yl(ar)=d*(ar-1)+d-1
1068         array_yl_labels(ar) = z_d(array_yl(ar))
1069         if (ar .GT. 0)
1070            do min_ar=0,3
1071               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)
1072            end do
1073         end if 
1074      end do
1075      do br=1,major_ticks_x-1
1076         array_xb(br)=e*(br-1)+e-1
1077         array_xb_labels(br) = x_d(array_xb(br))
1078         if (br .GT. 0)
1079            do min_br=0,3
1080               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)
1081            end do
1082         end if
1083      end do
1084   end if
1085
1086   if (yzc .EQ. 1) then
1087      d=(ze-zs+1)/(major_ticks_y-1)
1088      e=(ye-ys+1)/(major_ticks_x-1)
1089      array_yl       =new(major_ticks_y,integer)
1090      array_yl_labels=new(major_ticks_y,double)
1091      array_minor_yl =new((major_ticks_y-1)*4,double)
1092      array_xb       =new(major_ticks_x,integer)
1093      array_xb_labels=new(major_ticks_x,double)
1094      array_minor_xb =new((major_ticks_x-1)*4,double)
1095      array_yl(0)=zs
1096      array_xb(0)=ys
1097      array_yl_labels(0)=z_d(zs)
1098      array_xb_labels(0)=y_d(ys)
1099      do ar=1,major_ticks_y-1
1100         array_yl(ar)=d*(ar-1)+d-1
1101         array_yl_labels(ar) = y_d(array_yl(ar))
1102         if (ar .GT. 0)
1103            do min_ar=0,3
1104               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)
1105            end do
1106         end if 
1107      end do
1108      do br=1,major_ticks_x-1
1109         array_xb(br)=e*(br-1)+e-1
1110         array_xb_labels(br) = x_d(array_xb(br))
1111         if (br .GT. 0)
1112            do min_br=0,3
1113               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)
1114            end do
1115         end if
1116      end do
1117   end if
1118
1119   if (axes_explicit .EQ. 1)then
1120      cs_res@tmYLMode = "Explicit"
1121      cs_res@tmXBMode = "Explicit"
1122      cs_res@tmYLValues = array_yl
1123      if (xyc .EQ. 1)then
1124         cs_res@tmYLLabels = array_yl_labels/norm_y
1125         cs_res@tmXBLabels = array_xb_labels/norm_x
1126         if (norm_x .NE. 1.)then
1127            cs_res@tiXAxisString = "x ["+unit_x+"]"
1128         else
1129            cs_res@tiXAxisString = "x [m]"
1130         end if
1131         if (norm_y .NE. 1.)then
1132            cs_res@tiYAxisString = "y ["+unit_y+"]"
1133         else
1134            cs_res@tiYAxisString = "y [m]"
1135         end if   
1136      end if
1137      if (xzc .EQ. 1)then
1138         cs_res@tmYLLabels = array_yl_labels/norm_z
1139         cs_res@tmXBLabels = array_xb_labels/norm_x
1140         if (norm_x .NE. 1.)then
1141            cs_res@tiXAxisString = "x ["+unit_x+"]"
1142         else
1143            cs_res@tiXAxisString = "x [m]"
1144         end if
1145         if (norm_z .NE. 1.)then
1146            cs_res@tiYAxisString = "z ["+unit_z+"]"
1147         else
1148            cs_res@tiYAxisString = "z [m]"
1149         end if
1150      end if
1151      if (yzc .EQ. 1)then
1152         cs_res@tmYLLabels = array_yl_labels/norm_z
1153         cs_res@tmXBLabels = array_xb_labels/norm_y
1154         if (norm_y .NE. 1.)then
1155            cs_res@tiXAxisString = "y ["+unit_y+"]"
1156         else
1157            cs_res@tiXAxisString = "y [m]"
1158         end if
1159         if (norm_z .NE. 1.)then
1160            cs_res@tiYAxisString = "z ["+unit_z+"]"
1161         else
1162            cs_res@tiYAxisString = "z [m]"
1163         end if
1164      end if
1165      cs_res@tmXBValues = array_xb     
1166      cs_res@tmYLMinorValues = array_minor_yl
1167      cs_res@tmXBMinorValues = array_minor_xb
1168   else
1169      if (axes_explicit .EQ. 0)then 
1170         cs_res@tmYLMinorPerMajor = 4
1171         cs_res@tmXBMinorPerMajor = 4
1172      else
1173         print(" ")
1174         print("'axes_explicit' is invalid and set to 0")
1175         print(" ")
1176         axes_explicit = 0
1177         cs_res@tmYLMinorPerMajor = 4
1178         cs_res@tmXBMinorPerMajor = 4
1179      end if
1180      if (xyc .EQ. 1)then
[267]1181         cs_res@tiXAxisString = "x [m]"
1182         cs_res@tiYAxisString = "y [m]"
[218]1183      end if
1184      if (xzc .EQ. 1)then
[267]1185         cs_res@tiXAxisString = "x [m]"
1186         cs_res@tiYAxisString = "z [m]"
[218]1187      end if
1188      if (yzc .EQ. 1)then
[267]1189         cs_res@tiXAxisString = "y [m]"
1190         cs_res@tiYAxisString = "z [m]"
[218]1191      end if
1192   end if
1193
[162]1194   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1195      if (xe .EQ. xs+1) then
1196         print(" ")
[190]1197         print("range end for x-coordinate="+xe*delta_x+"m) must be at least two")
1198         print("more gridpoints(="+2*delta_x+"m) greater than start range="+xs*delta_x+"m)")
[162]1199         print(" ")
1200         exit
1201      end if
1202   end if
1203   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1204      if (ye .EQ. ys+1) then
1205         print(" ")
[190]1206         print("range end for y-coordinate="+ye*delta_y+"m) must be at least two")
1207         print("more gridpoints(="+2*delta_y+"m) greater than start range="+ys*delta_y+"m)")
[162]1208         print(" ")
1209         exit
1210      end if
1211   end if
1212   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1213      if (ze .EQ. zs+1) then
1214         print(" ")
[190]1215         print("range end for x-coordinate="+ze+") must be at least two")
1216         print("more gridpoints greater than start range="+zs+")")
[162]1217         print(" ")
1218         exit
1219      end if
1220   end if
1221 
[161]1222   if (xyc .EQ. 1) then
1223      no_layer = (ze-zs)+1
[190]1224      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1225   end if
1226   if (xzc .EQ. 1) then
1227      no_layer = (ye-ys)+1
[190]1228      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1229   end if
1230   if (yzc .EQ. 1) then
1231      no_layer = (xe-xs)+1
[190]1232      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1233   end if
1234
[162]1235   MinVal = new(dim,float)
[218]1236   MaxVal = new(dim,float)
1237   unit   = new(dim,string)
[157]1238
1239   ; ****************************************************
1240   ; define inner and outer loops depending on "sort"
1241   ; ****************************************************       
1242   
1243   if ( xyc .eq. 1 ) then
1244      lays = zs
1245      laye = ze
1246   end if
1247   if ( xzc .eq. 1 ) then
1248      lays = ys
1249      laye = ye
1250   end if
1251   if ( yzc .eq. 1) then
1252      lays = xs
1253      laye = xe
1254   end if
1255
1256   if ( sort .eq. "time" ) then
1257      lis = start_time_step
1258      lie = end_time_step
[190]1259      los = 0
1260      loe = laye-lays
[157]1261   else
[190]1262      lis = 0
1263      lie = laye-lays
[157]1264      los = start_time_step
1265      loe = end_time_step
1266   end if
1267
1268   check = True
1269   v1=0
1270   v2=0
[161]1271   no_var=0
1272   n=0
[157]1273
[190]1274   do varn=dim-1,0,1       
[162]1275   
[157]1276      if (vector .EQ. 1) then   
[161]1277         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1278         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
[157]1279      end if
1280           
[161]1281      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
[157]1282         check = False
[161]1283      else
1284         check = True
1285      end if 
[157]1286
[162]1287      if(check) then
[194]1288     
[161]1289         no_var=no_var+1
1290
[157]1291         if (vector .EQ. 1) then
[161]1292            if (","+vNam(varn)+"," .EQ. vec1) then
[157]1293               v1=v1+1
1294            end if
[161]1295            if (","+vNam(varn)+"," .EQ. vec2) then
[157]1296               v2=v2+1
1297            end if
1298         end if
1299
[162]1300         if (xyc .EQ. 1) then
[218]1301            temp = f[:]->$vNam(varn)$
1302            data_att = f_att->$vNam(varn)$
1303            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)               
[162]1304         end if
1305         if ( xzc .eq. 1 ) then
[218]1306            temp = f[:]->$vNam(varn)$
1307            data_att = f_att->$vNam(varn)$
1308            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[162]1309         end if
1310         if ( yzc .eq. 1) then
[218]1311            temp = f[:]->$vNam(varn)$
1312            data_att = f_att->$vNam(varn)$
1313            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[162]1314         end if
[157]1315         if (check_vec1) then
[190]1316            vect1=data(varn,:,:,:,:)
[154]1317         end if
[157]1318         if (check_vec2) then
[190]1319            vect2=data(varn,:,:,:,:)
[154]1320         end if
1321
[157]1322         data!0 = "var"
1323         data!1 = "t"
1324         data!2 = "z"
1325         data!3 = "y"
[162]1326         data!4 = "x" 
[154]1327
[162]1328         MinVal(varn) = min(data(varn,:,:,:,:))
1329         MaxVal(varn) = max(data(varn,:,:,:,:))
[218]1330         
1331         unit(varn) = data_att@units
[162]1332
[157]1333      end if
[154]1334
[190]1335   end do
[154]1336
[190]1337   var_input=new(no_var,string)
1338   numb=0
[194]1339   no_var=0
1340   
[190]1341   do varn=dim-1,0,1   
1342   
1343      if (vector .EQ. 1) then   
1344         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1345         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1346      end if
1347           
1348      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
1349         check = False
1350      else
1351         check = True
1352      end if 
1353
1354      if(check) then     
1355         var_input(numb)=vNam(varn)
1356         numb=numb+1     
1357      end if
[194]1358                 
1359      if (var .NE. "all") then
1360         check = isStrSubset( var,","+vNam(varn)+"," )
1361      end if
[190]1362     
[194]1363      if (check) then
1364         no_var = no_var+1
1365      end if
1366     
[161]1367   end do
[154]1368
[190]1369   if (no_var .EQ. 0) then
[162]1370      print(" ")
[190]1371      print("The variables var='"+var+"' do not exist on your input file;")
1372      print("be sure to have one comma berfore and after each variable")
[162]1373      print(" ")
1374      exit
1375   end if
1376
[157]1377   if (vector .EQ. 1) then
1378      if (v1 .EQ. 0)then
1379         print(" ")
[194]1380         print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
[190]1381         print("- "+var_input)
1382         print("be sure to have one comma berfore and after the variable")
[157]1383         print(" ")
1384         exit
1385      end if
[154]1386
[157]1387      if (v2 .EQ. 0)then
1388         print(" ")
[194]1389         print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
[190]1390         print("- "+var_input)
1391         print("be sure to have one comma berfore and after the variable")
[157]1392         print(" ")
1393         exit
1394      end if
1395   end if
[154]1396
[157]1397   ; ***************************************************
[161]1398   ; open workstation(s)
1399   ; ***************************************************
1400
1401   wks_ps  = gsn_open_wks(format_out,file_out)
1402   gsn_define_colormap(wks_ps,"rainbow+white")
1403
[218]1404   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1405      plot=new((/no_time*no_layer*no_var+no_time*no_layer/),graphic)
1406   else
1407      plot=new((/no_time*no_layer*no_var/),graphic)
1408   end if
1409
[161]1410   page = 0
1411   n=0
1412   print(" ")
[162]1413   print("Plots sorted by '"+sort+"'")
1414   print(" ")
[161]1415
1416   ; ***************************************************
[157]1417   ; create plots
1418   ; ***************************************************
1419
[190]1420   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[162]1421      do lo = los, loe                                         
1422         do li = lis, lie
[218]1423            if (xyc .EQ. 1)then
1424               if (sort .EQ. "time")then
1425                  level = "z=" + z_d(lo) + "m"
1426               else
1427                  level = "z=" + z_d(li) + "m"
1428               end if
1429            end if
1430            if (xzc .EQ. 1)then
1431               if (sort .EQ. "time")then
1432                  level = "y=" + y_d(lo) + "m"
1433               else
1434                  level = "y=" + y_d(li) + "m"
1435               end if
1436            end if
1437            if (yzc .EQ. 1)then
1438               if (sort .EQ. "time")then
1439                  level = "x=" + x_d(lo) + "m"
1440               else
1441                  level = "x=" + x_d(li) + "m"
1442               end if
1443            end if               
[162]1444            vecres                  = True            ; vector only resources
1445            vecres@gsnDraw          = False           ; don't draw
1446            vecres@gsnFrame         = False           ; don't advance frame
1447            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1448            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1449            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1450            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
[218]1451            vecres@tmXBLabelFontHeightF   = font_size
1452            vecres@tmYLLabelFontHeightF   = font_size
1453            vecres@tiXAxisFontHeightF     = font_size
1454            vecres@tiYAxisFontHeightF     = font_size
1455            if (sort .EQ. "time")then
1456               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1457            else
1458               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1459            end if
1460            vecres@tiXAxisString    = " "
1461            if (xyc .EQ. 1)then 
1462               if (sort .EQ. "time")then                                         
1463                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1464               else
1465                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1466               end if
1467            end if
1468            if (xzc .EQ. 1) then
1469               if (sort .EQ. "time")then
1470                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1471               else
1472                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1473               end if
1474            end if
1475            if (yzc .EQ. 1) then
1476               if (sort .EQ. "time")then
1477                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1478               else
1479                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1480               end if
1481            end if
[162]1482            n=n+1
1483         end do
1484      end do
1485   end if
1486 
[161]1487   do varn=dim-1,0,1
[157]1488
1489      if (vector .EQ. 1) then   
[161]1490         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
[157]1491      end if
1492
[161]1493      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
1494         check = False
1495      else
1496         check = True
1497      end if   
[190]1498 
1499      if (var .NE. "all") then
[161]1500         check = isStrSubset( var,","+vNam(varn)+"," )
[157]1501      end if
1502   
1503      if(check) then
[218]1504
[162]1505         space=(MaxVal(varn)-MinVal(varn))/24
1506 
1507         cs_res@cnMinLevelValF = MinVal(varn)
1508         cs_res@cnMaxLevelValF = MaxVal(varn)
1509
1510         cs_res@cnLevelSpacingF = space
1511     
[157]1512         ; ****************************************************
1513         ; loops over time and layer
1514         ; ****************************************************
1515
[162]1516         
[161]1517         do lo = los, loe                                       
[162]1518            do li = lis, lie
1519
[157]1520               ; ****************************************************
[154]1521               ; xy cross section
[157]1522               ; ****************************************************
[154]1523
[157]1524               if ( xyc .eq. 1 ) then
[218]1525         
[157]1526                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[218]1527                  cs_res@gsnRightString = unit(varn)
[162]1528                 
[154]1529                  if ( sort .eq. "time" ) then
[267]1530                     if ( z_d(zs+lo) .eq. -1)then
[194]1531                        if (delta_z .EQ. -1) then
1532                           level = "-average"
1533                        else
[267]1534                           level = "=" + z_d(zs+lo) + "m"
[194]1535                        end if
[154]1536                     else
[267]1537                        level = "=" + z_d(zs+lo) + "m"
[154]1538                     end if
[190]1539                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level             
[161]1540                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
[162]1541                     if (vector .EQ. 1 .AND. check_vecp) then
1542                        vecres                  = True            ; vector only resources
1543                        vecres@gsnDraw          = False           ; don't draw
1544                        vecres@gsnFrame         = False           ; don't advance frame
1545                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1546                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1547                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
1548                        vecres@gsnRightString   = " "
1549                        vecres@gsnLeftString    = " "
1550                        vecres@tiXAxisString    = " "   
1551                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1552                        overlay(plot(n), plot_vec)
[161]1553                     end if                         
[154]1554                  end if
[162]1555                 
[154]1556                  if ( sort .eq. "layer" ) then
[267]1557                     if (z_d(zs+li) .eq. -1) then
[194]1558                        if (delta_z .EQ. -1) then
1559                           level = "-average"
1560                        else
[267]1561                           level = "=" + z_d(zs+li) + "m"
[194]1562                        end if
[161]1563                     else
[267]1564                        level = "=" + z_d(zs+li) + "m"
[190]1565                     end if
1566                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level               
[161]1567                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
[162]1568                     if (vector .EQ. 1 .AND. check_vecp) then
1569                        vecres                  = True            ; vector only resources
1570                        vecres@gsnDraw          = False           ; don't draw
1571                        vecres@gsnFrame         = False           ; don't advance frame
1572                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1573                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1574                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1575                        vecres@gsnRightString   = " "             ; turn off right string
1576                        vecres@gsnLeftString    = " "             ; turn off left string
1577                        vecres@tiXAxisString    = " "   
1578                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1579                        overlay(plot(n), plot_vec)
[161]1580                     end if
[154]1581                  end if
1582               end if
1583
[157]1584               ; ****************************************************
[154]1585               ; xz cross section
[157]1586               ; ****************************************************
[154]1587
1588               if ( xzc .eq. 1 ) then
[218]1589                 
[157]1590                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[162]1591               
[154]1592                  if ( sort .eq. "time" ) then
[267]1593                     if ( y_d(ys+lo) .eq. -1 ) then
[154]1594                        level = "-average"
1595                     else
[267]1596                        level = "=" + y_d(ys+lo) + "m"
[154]1597                     end if
[218]1598                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
[161]1599                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
[162]1600                     if (vector .EQ. 1 .AND. check_vecp) then
1601                        vecres                  = True            ; vector only resources
1602                        vecres@gsnDraw          = False           ; don't draw
1603                        vecres@gsnFrame         = False           ; don't advance frame
1604                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1605                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1606                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1607                        vecres@gsnRightString   = " "             ; turn off right string
1608                        vecres@gsnLeftString    = " "             ; turn off left string
1609                        vecres@tiXAxisString    = " "   
[218]1610                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
[162]1611                        overlay(plot(n), plot_vec)
[161]1612                     end if
[154]1613                  end if
1614
1615                  if ( sort .eq. "layer" ) then
[267]1616                     if ( y_d(ys+li) .eq. -1 ) then
[161]1617                        level = "-average"
1618                     else
[267]1619                        level = "=" + y_d(ys+li) + "m"
[161]1620                     end if
[218]1621                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
[161]1622                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
[162]1623                     if (vector .EQ. 1 .AND. check_vecp) then
1624                        vecres                  = True            ; vector only resources
1625                        vecres@gsnDraw          = False           ; don't draw
1626                        vecres@gsnFrame         = False           ; don't advance frame
1627                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1628                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1629                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1630                        vecres@gsnRightString   = " "             ; turn off right string
1631                        vecres@gsnLeftString    = " "             ; turn off left string
1632                        vecres@tiXAxisString    = " "   
[218]1633                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
[162]1634                        overlay(plot(n), plot_vec)
[161]1635                     end if
[157]1636                  end if                 
[154]1637               end if
1638
[157]1639               ; ****************************************************
[154]1640               ; yz cross section
[157]1641               ; ****************************************************
[154]1642
1643               if ( yzc .eq. 1 ) then
[218]1644                                 
[157]1645                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1646
[154]1647                  if ( sort .eq. "time" ) then
[267]1648                     if ( x_d(xs+lo) .eq. -1 ) then
[154]1649                        level = "-average"
1650                     else
[267]1651                        level = "=" + x_d(xs+lo) + "m"
[154]1652                     end if
[218]1653                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
[161]1654                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
[162]1655                     if (vector .EQ. 1 .AND. check_vecp) then
1656                        vecres                  = True            ; vector only resources
1657                        vecres@gsnDraw          = False           ; don't draw
1658                        vecres@gsnFrame         = False           ; don't advance frame
1659                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1660                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1661                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1662                        vecres@gsnRightString   = " "             ; turn off right string
1663                        vecres@gsnLeftString    = " "             ; turn off left string
1664                        vecres@tiXAxisString    = " "   
[218]1665                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
[162]1666                        overlay(plot(n), plot_vec)
[161]1667                     end if
[154]1668                  end if
1669
1670                  if ( sort .eq. "layer" ) then
[267]1671                     if ( x_d(xs+li) .eq. -1 ) then
[161]1672                        level = "-average"
1673                     else
[267]1674                        level = "=" + x_d(xs+li) + "m"
[161]1675                     end if
[218]1676                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
[161]1677                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
[162]1678                     if (vector .EQ. 1 .AND. check_vecp)then
1679                        vecres                  = True            ; vector only resources
1680                        vecres@gsnDraw          = False           ; don't draw
1681                        vecres@gsnFrame         = False           ; don't advance frame
1682                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1683                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1684                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1685                        vecres@gsnRightString   = " "             ; turn off right string
1686                        vecres@gsnLeftString    = " "             ; turn off left string
1687                        vecres@tiXAxisString    = " "   
[218]1688                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
[162]1689                        overlay(plot(n), plot_vec)
[161]1690                     end if
[154]1691                  end if
[161]1692               end if         
[218]1693               n=n+1 
[154]1694            end do     
1695         end do 
[162]1696      end if
[161]1697   end do
[154]1698
[161]1699   ; ***************************************************
1700   ; merge plots onto one page
1701   ; ***************************************************
[194]1702   
[190]1703   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
[162]1704      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1705         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1706         print(" ")
1707         print("Outputs to .eps or .epsi have only one frame")
1708         print(" ")
1709      else
[218]1710         do np = 0,no_layer*no_time*(no_var+1)-1,no_rows*no_columns   
1711            if ( np + no_rows*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then
1712               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_rows,no_columns/),cs_resP)
[162]1713            else
[218]1714               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
[162]1715            end if
1716         end do
1717      end if
[267]1718   else       
[162]1719      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1720         gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
1721         print(" ")
1722         print("Outputs to .eps or .epsi have only one frame")
1723         print(" ")
1724      else
[267]1725         do np = 0,no_layer*no_time*no_var-1,no_rows*no_columns 
[218]1726            if ( np + no_rows*no_columns .gt. (no_layer*no_time*no_var)-1) then
1727               gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_rows,no_columns/),cs_resP)
[162]1728            else
[218]1729               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
[162]1730            end if
1731         end do
1732      end if
[161]1733   end if
1734
[154]1735   print(" ")
1736   print("Output to: " + file_out +"."+ format_out)
1737   print(" ")   
1738 
[190]1739end
Note: See TracBrowser for help on using the repository browser.