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

Last change on this file since 346 was 329, checked in by heinze, 16 years ago

Small change in formatting in message.f90. Bugfix in cross_sections.ncl in case of normalizing.

File size: 60.0 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
[329]1021      d= (ye-ys+1)/(major_ticks_y-1)
[218]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
[329]1038                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)
[218]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
[329]1047               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)
[218]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
[329]1071               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)
[218]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
[329]1080               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)
[218]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
[329]1101         array_yl_labels(ar) = z_d(array_yl(ar))
[218]1102         if (ar .GT. 0)
1103            do min_ar=0,3
[329]1104               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)
[218]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
[329]1110         array_xb_labels(br) = y_d(array_xb(br))
[218]1111         if (br .GT. 0)
1112            do min_br=0,3
[329]1113               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)
[218]1114            end do
1115         end if
1116      end do
1117   end if
[329]1118 
[218]1119   if (axes_explicit .EQ. 1)then
1120      cs_res@tmYLMode = "Explicit"
1121      cs_res@tmXBMode = "Explicit"
[329]1122      if (xyc .EQ. 1)then
1123         cs_res@tmYLValues = y_d(array_yl)
1124         cs_res@tmXBValues = x_d(array_xb)     
[218]1125         cs_res@tmYLLabels = array_yl_labels/norm_y
1126         cs_res@tmXBLabels = array_xb_labels/norm_x
1127         if (norm_x .NE. 1.)then
1128            cs_res@tiXAxisString = "x ["+unit_x+"]"
1129         else
1130            cs_res@tiXAxisString = "x [m]"
1131         end if
1132         if (norm_y .NE. 1.)then
1133            cs_res@tiYAxisString = "y ["+unit_y+"]"
1134         else
1135            cs_res@tiYAxisString = "y [m]"
1136         end if   
1137      end if
1138      if (xzc .EQ. 1)then
[329]1139         cs_res@tmYLValues = z_d(array_yl)
1140         cs_res@tmXBValues = x_d(array_xb) 
[218]1141         cs_res@tmYLLabels = array_yl_labels/norm_z
1142         cs_res@tmXBLabels = array_xb_labels/norm_x
1143         if (norm_x .NE. 1.)then
1144            cs_res@tiXAxisString = "x ["+unit_x+"]"
1145         else
1146            cs_res@tiXAxisString = "x [m]"
1147         end if
1148         if (norm_z .NE. 1.)then
1149            cs_res@tiYAxisString = "z ["+unit_z+"]"
1150         else
1151            cs_res@tiYAxisString = "z [m]"
1152         end if
1153      end if
1154      if (yzc .EQ. 1)then
[329]1155         cs_res@tmYLValues = z_d(array_yl)
1156         cs_res@tmXBValues = y_d(array_xb) 
[218]1157         cs_res@tmYLLabels = array_yl_labels/norm_z
1158         cs_res@tmXBLabels = array_xb_labels/norm_y
1159         if (norm_y .NE. 1.)then
1160            cs_res@tiXAxisString = "y ["+unit_y+"]"
1161         else
1162            cs_res@tiXAxisString = "y [m]"
1163         end if
1164         if (norm_z .NE. 1.)then
1165            cs_res@tiYAxisString = "z ["+unit_z+"]"
1166         else
1167            cs_res@tiYAxisString = "z [m]"
1168         end if
1169      end if
1170      cs_res@tmYLMinorValues = array_minor_yl
1171      cs_res@tmXBMinorValues = array_minor_xb
1172   else
1173      if (axes_explicit .EQ. 0)then 
1174         cs_res@tmYLMinorPerMajor = 4
1175         cs_res@tmXBMinorPerMajor = 4
1176      else
1177         print(" ")
1178         print("'axes_explicit' is invalid and set to 0")
1179         print(" ")
1180         axes_explicit = 0
1181         cs_res@tmYLMinorPerMajor = 4
1182         cs_res@tmXBMinorPerMajor = 4
1183      end if
1184      if (xyc .EQ. 1)then
[267]1185         cs_res@tiXAxisString = "x [m]"
1186         cs_res@tiYAxisString = "y [m]"
[218]1187      end if
1188      if (xzc .EQ. 1)then
[267]1189         cs_res@tiXAxisString = "x [m]"
1190         cs_res@tiYAxisString = "z [m]"
[218]1191      end if
1192      if (yzc .EQ. 1)then
[267]1193         cs_res@tiXAxisString = "y [m]"
1194         cs_res@tiYAxisString = "z [m]"
[218]1195      end if
1196   end if
1197
[162]1198   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1199      if (xe .EQ. xs+1) then
1200         print(" ")
[190]1201         print("range end for x-coordinate="+xe*delta_x+"m) must be at least two")
1202         print("more gridpoints(="+2*delta_x+"m) greater than start range="+xs*delta_x+"m)")
[162]1203         print(" ")
1204         exit
1205      end if
1206   end if
1207   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1208      if (ye .EQ. ys+1) then
1209         print(" ")
[190]1210         print("range end for y-coordinate="+ye*delta_y+"m) must be at least two")
1211         print("more gridpoints(="+2*delta_y+"m) greater than start range="+ys*delta_y+"m)")
[162]1212         print(" ")
1213         exit
1214      end if
1215   end if
1216   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1217      if (ze .EQ. zs+1) then
1218         print(" ")
[190]1219         print("range end for x-coordinate="+ze+") must be at least two")
1220         print("more gridpoints greater than start range="+zs+")")
[162]1221         print(" ")
1222         exit
1223      end if
1224   end if
1225 
[161]1226   if (xyc .EQ. 1) then
1227      no_layer = (ze-zs)+1
[190]1228      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1229   end if
1230   if (xzc .EQ. 1) then
1231      no_layer = (ye-ys)+1
[190]1232      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1233   end if
1234   if (yzc .EQ. 1) then
1235      no_layer = (xe-xs)+1
[190]1236      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
[161]1237   end if
1238
[162]1239   MinVal = new(dim,float)
[218]1240   MaxVal = new(dim,float)
1241   unit   = new(dim,string)
[157]1242
1243   ; ****************************************************
1244   ; define inner and outer loops depending on "sort"
1245   ; ****************************************************       
1246   
1247   if ( xyc .eq. 1 ) then
1248      lays = zs
1249      laye = ze
1250   end if
1251   if ( xzc .eq. 1 ) then
1252      lays = ys
1253      laye = ye
1254   end if
1255   if ( yzc .eq. 1) then
1256      lays = xs
1257      laye = xe
1258   end if
1259
1260   if ( sort .eq. "time" ) then
1261      lis = start_time_step
1262      lie = end_time_step
[190]1263      los = 0
1264      loe = laye-lays
[157]1265   else
[190]1266      lis = 0
1267      lie = laye-lays
[157]1268      los = start_time_step
1269      loe = end_time_step
1270   end if
1271
1272   check = True
1273   v1=0
1274   v2=0
[161]1275   no_var=0
1276   n=0
[157]1277
[190]1278   do varn=dim-1,0,1       
[162]1279   
[157]1280      if (vector .EQ. 1) then   
[161]1281         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1282         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
[157]1283      end if
1284           
[161]1285      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]1286         check = False
[161]1287      else
1288         check = True
1289      end if 
[157]1290
[162]1291      if(check) then
[194]1292     
[161]1293         no_var=no_var+1
1294
[157]1295         if (vector .EQ. 1) then
[161]1296            if (","+vNam(varn)+"," .EQ. vec1) then
[157]1297               v1=v1+1
1298            end if
[161]1299            if (","+vNam(varn)+"," .EQ. vec2) then
[157]1300               v2=v2+1
1301            end if
1302         end if
1303
[162]1304         if (xyc .EQ. 1) then
[218]1305            temp = f[:]->$vNam(varn)$
1306            data_att = f_att->$vNam(varn)$
[329]1307            if(vNam(varn) .eq. "ts_xy" .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy")
1308              ;these variables depend von zu1_xy and that's why they have only one z-layer
1309              data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1310            else
1311              data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1312            end if
1313            delete(temp)               
[162]1314         end if
1315         if ( xzc .eq. 1 ) then
[218]1316            temp = f[:]->$vNam(varn)$
1317            data_att = f_att->$vNam(varn)$
1318            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[162]1319         end if
1320         if ( yzc .eq. 1) then
[218]1321            temp = f[:]->$vNam(varn)$
1322            data_att = f_att->$vNam(varn)$
1323            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[162]1324         end if
[157]1325         if (check_vec1) then
[190]1326            vect1=data(varn,:,:,:,:)
[154]1327         end if
[157]1328         if (check_vec2) then
[190]1329            vect2=data(varn,:,:,:,:)
[154]1330         end if
1331
[157]1332         data!0 = "var"
1333         data!1 = "t"
1334         data!2 = "z"
1335         data!3 = "y"
[162]1336         data!4 = "x" 
[154]1337
[162]1338         MinVal(varn) = min(data(varn,:,:,:,:))
1339         MaxVal(varn) = max(data(varn,:,:,:,:))
[218]1340         
1341         unit(varn) = data_att@units
[329]1342         delete(data_att)
[162]1343
[157]1344      end if
[154]1345
[190]1346   end do
[154]1347
[190]1348   var_input=new(no_var,string)
1349   numb=0
[194]1350   no_var=0
1351   
[190]1352   do varn=dim-1,0,1   
1353   
1354      if (vector .EQ. 1) then   
1355         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1356         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1357      end if
1358           
1359      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
1360         check = False
1361      else
1362         check = True
1363      end if 
1364
1365      if(check) then     
1366         var_input(numb)=vNam(varn)
1367         numb=numb+1     
1368      end if
[194]1369                 
1370      if (var .NE. "all") then
1371         check = isStrSubset( var,","+vNam(varn)+"," )
1372      end if
[190]1373     
[194]1374      if (check) then
1375         no_var = no_var+1
1376      end if
1377     
[161]1378   end do
[154]1379
[190]1380   if (no_var .EQ. 0) then
[162]1381      print(" ")
[190]1382      print("The variables var='"+var+"' do not exist on your input file;")
1383      print("be sure to have one comma berfore and after each variable")
[162]1384      print(" ")
1385      exit
1386   end if
1387
[157]1388   if (vector .EQ. 1) then
1389      if (v1 .EQ. 0)then
1390         print(" ")
[194]1391         print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
[190]1392         print("- "+var_input)
1393         print("be sure to have one comma berfore and after the variable")
[157]1394         print(" ")
1395         exit
1396      end if
[154]1397
[157]1398      if (v2 .EQ. 0)then
1399         print(" ")
[194]1400         print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
[190]1401         print("- "+var_input)
1402         print("be sure to have one comma berfore and after the variable")
[157]1403         print(" ")
1404         exit
1405      end if
1406   end if
[154]1407
[157]1408   ; ***************************************************
[161]1409   ; open workstation(s)
1410   ; ***************************************************
1411
1412   wks_ps  = gsn_open_wks(format_out,file_out)
1413   gsn_define_colormap(wks_ps,"rainbow+white")
1414
[218]1415   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1416      plot=new((/no_time*no_layer*no_var+no_time*no_layer/),graphic)
1417   else
1418      plot=new((/no_time*no_layer*no_var/),graphic)
1419   end if
1420
[161]1421   page = 0
1422   n=0
1423   print(" ")
[162]1424   print("Plots sorted by '"+sort+"'")
1425   print(" ")
[161]1426
1427   ; ***************************************************
[157]1428   ; create plots
1429   ; ***************************************************
1430
[190]1431   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[162]1432      do lo = los, loe                                         
1433         do li = lis, lie
[218]1434            if (xyc .EQ. 1)then
1435               if (sort .EQ. "time")then
1436                  level = "z=" + z_d(lo) + "m"
1437               else
1438                  level = "z=" + z_d(li) + "m"
1439               end if
1440            end if
1441            if (xzc .EQ. 1)then
1442               if (sort .EQ. "time")then
1443                  level = "y=" + y_d(lo) + "m"
1444               else
1445                  level = "y=" + y_d(li) + "m"
1446               end if
1447            end if
1448            if (yzc .EQ. 1)then
1449               if (sort .EQ. "time")then
1450                  level = "x=" + x_d(lo) + "m"
1451               else
1452                  level = "x=" + x_d(li) + "m"
1453               end if
1454            end if               
[162]1455            vecres                  = True            ; vector only resources
1456            vecres@gsnDraw          = False           ; don't draw
1457            vecres@gsnFrame         = False           ; don't advance frame
1458            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1459            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1460            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1461            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
[218]1462            vecres@tmXBLabelFontHeightF   = font_size
1463            vecres@tmYLLabelFontHeightF   = font_size
1464            vecres@tiXAxisFontHeightF     = font_size
1465            vecres@tiYAxisFontHeightF     = font_size
1466            if (sort .EQ. "time")then
1467               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1468            else
1469               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1470            end if
1471            vecres@tiXAxisString    = " "
1472            if (xyc .EQ. 1)then 
1473               if (sort .EQ. "time")then                                         
1474                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1475               else
1476                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1477               end if
1478            end if
1479            if (xzc .EQ. 1) then
1480               if (sort .EQ. "time")then
1481                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1482               else
1483                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1484               end if
1485            end if
1486            if (yzc .EQ. 1) then
1487               if (sort .EQ. "time")then
1488                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1489               else
1490                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1491               end if
1492            end if
[162]1493            n=n+1
1494         end do
1495      end do
1496   end if
1497 
[161]1498   do varn=dim-1,0,1
[157]1499
1500      if (vector .EQ. 1) then   
[161]1501         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
[157]1502      end if
1503
[161]1504      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
1505         check = False
1506      else
1507         check = True
1508      end if   
[190]1509 
1510      if (var .NE. "all") then
[161]1511         check = isStrSubset( var,","+vNam(varn)+"," )
[157]1512      end if
1513   
1514      if(check) then
[218]1515
[162]1516         space=(MaxVal(varn)-MinVal(varn))/24
1517 
1518         cs_res@cnMinLevelValF = MinVal(varn)
1519         cs_res@cnMaxLevelValF = MaxVal(varn)
1520
1521         cs_res@cnLevelSpacingF = space
1522     
[157]1523         ; ****************************************************
1524         ; loops over time and layer
1525         ; ****************************************************
1526
[162]1527         
[161]1528         do lo = los, loe                                       
[162]1529            do li = lis, lie
1530
[157]1531               ; ****************************************************
[154]1532               ; xy cross section
[157]1533               ; ****************************************************
[154]1534
[157]1535               if ( xyc .eq. 1 ) then
[218]1536         
[157]1537                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[218]1538                  cs_res@gsnRightString = unit(varn)
[162]1539                 
[154]1540                  if ( sort .eq. "time" ) then
[267]1541                     if ( z_d(zs+lo) .eq. -1)then
[194]1542                        if (delta_z .EQ. -1) then
1543                           level = "-average"
1544                        else
[267]1545                           level = "=" + z_d(zs+lo) + "m"
[194]1546                        end if
[154]1547                     else
[267]1548                        level = "=" + z_d(zs+lo) + "m"
[154]1549                     end if
[190]1550                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level             
[161]1551                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
[162]1552                     if (vector .EQ. 1 .AND. check_vecp) then
1553                        vecres                  = True            ; vector only resources
1554                        vecres@gsnDraw          = False           ; don't draw
1555                        vecres@gsnFrame         = False           ; don't advance frame
1556                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1557                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1558                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
1559                        vecres@gsnRightString   = " "
1560                        vecres@gsnLeftString    = " "
1561                        vecres@tiXAxisString    = " "   
1562                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1563                        overlay(plot(n), plot_vec)
[161]1564                     end if                         
[154]1565                  end if
[162]1566                 
[154]1567                  if ( sort .eq. "layer" ) then
[267]1568                     if (z_d(zs+li) .eq. -1) then
[194]1569                        if (delta_z .EQ. -1) then
1570                           level = "-average"
1571                        else
[267]1572                           level = "=" + z_d(zs+li) + "m"
[194]1573                        end if
[161]1574                     else
[267]1575                        level = "=" + z_d(zs+li) + "m"
[190]1576                     end if
1577                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level               
[161]1578                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
[162]1579                     if (vector .EQ. 1 .AND. check_vecp) then
1580                        vecres                  = True            ; vector only resources
1581                        vecres@gsnDraw          = False           ; don't draw
1582                        vecres@gsnFrame         = False           ; don't advance frame
1583                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1584                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1585                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1586                        vecres@gsnRightString   = " "             ; turn off right string
1587                        vecres@gsnLeftString    = " "             ; turn off left string
1588                        vecres@tiXAxisString    = " "   
1589                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1590                        overlay(plot(n), plot_vec)
[161]1591                     end if
[154]1592                  end if
1593               end if
1594
[157]1595               ; ****************************************************
[154]1596               ; xz cross section
[157]1597               ; ****************************************************
[154]1598
1599               if ( xzc .eq. 1 ) then
[218]1600                 
[157]1601                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[162]1602               
[154]1603                  if ( sort .eq. "time" ) then
[267]1604                     if ( y_d(ys+lo) .eq. -1 ) then
[154]1605                        level = "-average"
1606                     else
[267]1607                        level = "=" + y_d(ys+lo) + "m"
[154]1608                     end if
[218]1609                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
[161]1610                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
[162]1611                     if (vector .EQ. 1 .AND. check_vecp) then
1612                        vecres                  = True            ; vector only resources
1613                        vecres@gsnDraw          = False           ; don't draw
1614                        vecres@gsnFrame         = False           ; don't advance frame
1615                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1616                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1617                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1618                        vecres@gsnRightString   = " "             ; turn off right string
1619                        vecres@gsnLeftString    = " "             ; turn off left string
1620                        vecres@tiXAxisString    = " "   
[218]1621                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
[162]1622                        overlay(plot(n), plot_vec)
[161]1623                     end if
[154]1624                  end if
1625
1626                  if ( sort .eq. "layer" ) then
[267]1627                     if ( y_d(ys+li) .eq. -1 ) then
[161]1628                        level = "-average"
1629                     else
[267]1630                        level = "=" + y_d(ys+li) + "m"
[161]1631                     end if
[218]1632                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
[161]1633                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
[162]1634                     if (vector .EQ. 1 .AND. check_vecp) then
1635                        vecres                  = True            ; vector only resources
1636                        vecres@gsnDraw          = False           ; don't draw
1637                        vecres@gsnFrame         = False           ; don't advance frame
1638                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1639                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1640                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1641                        vecres@gsnRightString   = " "             ; turn off right string
1642                        vecres@gsnLeftString    = " "             ; turn off left string
1643                        vecres@tiXAxisString    = " "   
[218]1644                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
[162]1645                        overlay(plot(n), plot_vec)
[161]1646                     end if
[157]1647                  end if                 
[154]1648               end if
1649
[157]1650               ; ****************************************************
[154]1651               ; yz cross section
[157]1652               ; ****************************************************
[154]1653
1654               if ( yzc .eq. 1 ) then
[218]1655                                 
[157]1656                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1657
[154]1658                  if ( sort .eq. "time" ) then
[267]1659                     if ( x_d(xs+lo) .eq. -1 ) then
[154]1660                        level = "-average"
1661                     else
[267]1662                        level = "=" + x_d(xs+lo) + "m"
[154]1663                     end if
[218]1664                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
[161]1665                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
[162]1666                     if (vector .EQ. 1 .AND. check_vecp) then
1667                        vecres                  = True            ; vector only resources
1668                        vecres@gsnDraw          = False           ; don't draw
1669                        vecres@gsnFrame         = False           ; don't advance frame
1670                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1671                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1672                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1673                        vecres@gsnRightString   = " "             ; turn off right string
1674                        vecres@gsnLeftString    = " "             ; turn off left string
1675                        vecres@tiXAxisString    = " "   
[218]1676                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
[162]1677                        overlay(plot(n), plot_vec)
[161]1678                     end if
[154]1679                  end if
1680
1681                  if ( sort .eq. "layer" ) then
[267]1682                     if ( x_d(xs+li) .eq. -1 ) then
[161]1683                        level = "-average"
1684                     else
[267]1685                        level = "=" + x_d(xs+li) + "m"
[161]1686                     end if
[218]1687                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
[161]1688                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
[162]1689                     if (vector .EQ. 1 .AND. check_vecp)then
1690                        vecres                  = True            ; vector only resources
1691                        vecres@gsnDraw          = False           ; don't draw
1692                        vecres@gsnFrame         = False           ; don't advance frame
1693                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1694                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1695                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1696                        vecres@gsnRightString   = " "             ; turn off right string
1697                        vecres@gsnLeftString    = " "             ; turn off left string
1698                        vecres@tiXAxisString    = " "   
[218]1699                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
[162]1700                        overlay(plot(n), plot_vec)
[161]1701                     end if
[154]1702                  end if
[161]1703               end if         
[218]1704               n=n+1 
[154]1705            end do     
1706         end do 
[162]1707      end if
[161]1708   end do
[154]1709
[161]1710   ; ***************************************************
1711   ; merge plots onto one page
1712   ; ***************************************************
[194]1713   
[190]1714   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
[162]1715      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1716         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1717         print(" ")
1718         print("Outputs to .eps or .epsi have only one frame")
1719         print(" ")
1720      else
[218]1721         do np = 0,no_layer*no_time*(no_var+1)-1,no_rows*no_columns   
1722            if ( np + no_rows*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then
1723               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_rows,no_columns/),cs_resP)
[162]1724            else
[218]1725               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
[162]1726            end if
1727         end do
1728      end if
[267]1729   else       
[162]1730      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1731         gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
1732         print(" ")
1733         print("Outputs to .eps or .epsi have only one frame")
1734         print(" ")
1735      else
[267]1736         do np = 0,no_layer*no_time*no_var-1,no_rows*no_columns 
[218]1737            if ( np + no_rows*no_columns .gt. (no_layer*no_time*no_var)-1) then
1738               gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_rows,no_columns/),cs_resP)
[162]1739            else
[218]1740               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
[162]1741            end if
1742         end do
1743      end if
[161]1744   end if
1745
[154]1746   print(" ")
1747   print("Output to: " + file_out +"."+ format_out)
1748   print(" ")   
1749 
[190]1750end
Note: See TracBrowser for help on using the repository browser.