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

Last change on this file since 1262 was 1248, checked in by heinze, 11 years ago

NCL function getfilevarnames changed with NCL version 6.1.1 and higher -> adaptation required

File size: 78.1 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;***************************************************
[418]8; load .ncl.config or .ncl.config.default
[190]9;***************************************************
[194]10
11function which_script()
12local script
13begin
14   script="cross_section"
15   return(script)
16end
[418]17
18if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
19   loadscript("$PALM_BIN/../../.ncl.config")
[190]20else
[418]21  if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
22     loadscript( "$PALM_BIN/NCL/.ncl.config.default")
[190]23  else
[418]24      palm_bin_path = getenv("PALM_BIN")
[190]25      print(" ")
[418]26      print("Neither the personal configuration file '.ncl.config' exists in")
27      print("~/palm/current_version")
[566]28      print("nor the default configuration file '.ncl.config.default' "+\
29            "exists in")
[418]30      print(palm_bin_path + "/NCL")
[190]31      print(" ")
32      exit
33   end if
[418]34end if   
[194]35
[154]36begin
37
[532]38   ;***************************************************
39   ; Retrieving the NCL version used
40   ;***************************************************
41   
42   ncl_version_ch = systemfunc("ncl -V")
43   ncl_version    = stringtofloat(ncl_version_ch)
44
[154]45   ; ***************************************************
[526]46   ; Retrieving the double quote character
47   ; ***************************************************
48   
49   dq=str_get_dq()
50
51   ; ***************************************************
[566]52   ; set default parameter values and strings if not assigned
53   ; in prompt or parameter list
[154]54   ; ***************************************************
[566]55
56   ; Selection of type of cross section is done in palmplot
57   if (.not. isvar("xyc"))then       
58      xyc = 0   
59   end if
60   if (.not. isvar("xzc"))then 
61      xzc = 0   
62   end if
63   if (.not. isvar("yzc"))then     
64      yzc = 0
65   end if
[154]66   
[190]67   if (file_1 .EQ. "File in") then
68      print(" ")
[418]69      print("Declare input file 'file_1=' in '.ncl.config' or prompt")
[190]70      print(" ")
71      exit
[175]72   else
73      file_in = file_1   
[154]74   end if
[175]75
[566]76   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND.   \
77       format_out .NE. "eps" .AND. format_out .NE. "ps" .AND.    \
78       format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
79       format_out .NE. "png") then
[190]80      print(" ")
81      print("'format_out = "+format_out+"' is invalid and set to'x11'")
82      print(" ")
83      format_out="x11"
[154]84   end if
[175]85
[532]86   if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
87      print(" ")
88      print("Output of png files not available")
89      print("png output is avaiable with NCL version 5.2.0 and higher ")
90      print("NCL version used: " + ncl_version_ch)
91      print(" ")
92      exit
93   end if
94
[190]95   if (sort .NE. "layer" .AND. sort .NE. "time") then 
96      print(" ")
97      print("'sort'= "+sort+" is invalid and set to 'layer'")
98      print(" ")
99      sort = "layer" 
[154]100   end if
[190]101   
102   if (mode .NE. "Fill" .AND. mode .NE. "Line" .AND. mode .NE. "Both") then
103      print(" ")
104      print("'mode'= "+mode+" is invalid and set to 'Fill'")
105      print(" ")
[154]106      mode = "Fill"
107   end if
[175]108
[566]109   if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND.\
110       fill_mode .NE. "CellFill") then
[190]111      print(" ")
[218]112      print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
[190]113      print(" ")
[154]114      fill_mode = "AreaFill"
115   end if
[190]116   
117   if (shape .NE. 0 .AND. shape .NE. 1) then
118      print(" ")
119      print("'shape'= "+shape+" is invalid and set to 1")
120      print(" ")
[154]121      shape = 1
[175]122   end if
[513]123
[190]124   if (xyc .NE. 0 .AND. xyc .NE. 1)then
125      print(" ")
126      print("'xyc'= "+xyc+" is invalid and set to 0")
127      print(" ")
128      xyc = 0 
[154]129   end if
[190]130   
131   if (xzc .NE. 0 .AND. xzc .NE. 1)then
132      print(" ")
133      print("'xzc'= "+xzc+" is invalid and set to 0")
134      print(" ")
135      xyc = 0 
[154]136   end if
[190]137   
138   if (yzc .NE. 0 .AND. yzc .NE. 1)then
139      print(" ")
140      print("'yzc'= "+yzc+" is invalid and set to 0")
141      print(" ")
142      xyc = 0 
[154]143   end if
[190]144   
[154]145   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
146      print(" ")
[218]147      print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]148      print(" ")
149      exit
150   end if
151   if (xyc .EQ. 1 ) then
152      if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
153         print(" ")
[218]154         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]155         print(" ")
156         exit
157      end if
158   end if
159   if (xzc .EQ. 1) then
160      if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
161         print(" ")
[218]162         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]163         print(" ")
164         exit
165      end if
166   end if
167   if (yzc .EQ. 1) then
168      if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
169         print(" ")
[218]170         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
[154]171         print(" ")
172         exit
173      end if
174   end if
175
[190]176   if (vector .NE. 0 .AND. vector .NE. 1) then
177      print(" ")
[218]178      print("Set 'vector' to 0 or 1")
[190]179      print(" ")
180      exit
[218]181   end if 
182   
183   if (norm_x .EQ. 0) then
184      print(" ")
185      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
186      print(" ")
187      norm_x = 1.0
188   end if
189   if (norm_y .EQ. 0) then
190      print(" ")
191      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
192      print(" ")
193      norm_y = 1.0
194   end if
195   if (norm_z .EQ. 0) then
196      print(" ")
197      print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
198      print(" ")
199      norm_z = 1.0
200   end if
[190]201 
[154]202   ; ***************************************************
203   ; open input file
204   ; ***************************************************
[218]205   
206   file_in_1 = False
207   if (isStrSubset(file_in, ".nc"))then
208      start_f = -2
209      end_f = -2
210      file_in_1 = True     
211   end if
[154]212
[218]213   if (start_f .EQ. -1)then
214      print(" ")
[566]215      print("'start_f' must be one of the cyclic numbers (at least 0) "+\
216            "of your input file(s)")
[218]217      print(" ") 
218      exit
219   end if
220   if (end_f .EQ. -1)then
221      print(" ")
[566]222      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
223            "your input file(s)")
[218]224      print(" ") 
225      exit
226   end if
[154]227
[218]228   files=new(end_f-start_f+1,string)
229
230   if (file_in_1)then
231      if (isfilepresent(file_in))then
232         files(0)=file_in
233      else
234         print(" ")
235         print("1st input file: '"+file_in+"' does not exist")
236         print(" ")
237         exit
238      end if
239   else   
240      if (start_f .EQ. 0)then
241         if (isfilepresent(file_in+".nc"))then
242            files(0)=file_in+".nc"
243            do i=1,end_f
244               if (isfilepresent(file_in+"."+i+".nc"))then   
245                  files(i)=file_in+"."+i+".nc"
246               else
247                  print(" ")
248                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
249                  print(" ")
250                  exit 
251               end if       
252            end do         
253         else
254            print(" ")
255            print("Input file: '"+file_in+".nc' does not exist")
256            print(" ")
257            exit
258         end if
259      else
260         do i=start_f,end_f
261            if (isfilepresent(file_in+"."+i+".nc"))then   
262               files(i-start_f)=file_in+"."+i+".nc"
263            else
264               print(" ")
265               print("Input file: '"+file_in+"."+i+".nc' does not exist")
266               print(" ")
267               exit 
268            end if
269         end do
270      end if
271   end if
272   
273   f=addfiles(files,"r")
274   f_att=addfile(files(0),"r")
275   ListSetType(f,"cat")
276
277   vNam  = getfilevarnames(f_att)
[1248]278   if( ncl_version .GE. 6.1 ) then
279      vNam = vNam(::-1)
280   end if
[585]281   vType = getfilevartypes(f_att,vNam)
282
283   if ((all(vType .eq. "double"))) then ;distinction if data is double or float
284      check_vType = True
285   else
286      check_vType = False
287   end if
288
[154]289   print(" ")
[190]290   print("Variables in input file:")
291   print("- "+ vNam)
[154]292   print(" ")
293   dim   = dimsizes(vNam)
[162]294
295   ; ***************************************************
296   ; check for kind of input file
297   ; ***************************************************
298
299   check_3d = True
300   do varn=0,dim-1
301      checkxy = isStrSubset( vNam(varn),"xy")
302      checkxz = isStrSubset( vNam(varn),"xz")
303      checkyz = isStrSubset( vNam(varn),"yz")
304      if (checkxy .OR. checkxz .OR. checkyz) then
305         check_3d = False
306         break
307      end if
308   end do
309   if (.not. check_3d) then
[195]310      if (xyc .EQ. 1)then
311         if (.not. checkxy) then
312            print(" ")
313            print("Your input file doesn't have values for xy cross-sections")
314            if (checkxz)then
315               print("Select another input file or xz cross-sections")
316            else
317               print("Select another input file or yz cross-sections")
318            end if
319            print(" ")
320            exit
[162]321         else
[195]322            print(" ")
323            print("Your input file contains xy data")
324            print(" ")
325         end if
[162]326      end if
[195]327      if (xzc .EQ. 1)then
328         if (.not. checkxz) then
329            print(" ")
330            print("Your input file doesn't have values for xz cross-sections")
331            if (checkxy)then
332               print("Select another input file or xy cross-sections")
333            else
334               print("Select another input file or yz cross-sections")
335            end if
336            print(" ")
337            exit
[162]338         else
[195]339            print(" ")
340            print("Your input file contains xz data")
341            print(" ")
342         end if
[162]343      end if
[195]344      if (yzc .EQ. 1)then
345         if (.not. checkyz) then
346            print(" ")
347            print("Your input file doesn't have values for yz cross-sections")
348            if (checkxy)then
349               print("Select another input file or xy cross-sections")
350            else
351               print("Select another input file or xz cross-sections")
352            end if
353            print(" ")
354            exit
[162]355         else
[195]356            print(" ")
357            print("Your input file contains yz data")
358            print(" ")
359         end if
[162]360      end if
361   else
362      print(" ")
[218]363      print("Your input file contains 3d or other data")
[162]364      print(" ")
365   end if
366   
[161]367   if (dim .EQ. 0) then
368      print(" ")
[162]369      print("There is no data on file")
[161]370      print(" ")
[162]371      exit
[161]372   end if
373
[154]374   ; ***************************************************
375   ; set up recourses
376   ; ***************************************************
377
378   cs_res                          = True
[218]379   vecres                          = True
[154]380   cs_res@gsnYAxisIrregular2Linear = True 
[566]381   vecres@gsnYAxisIrregular2Linear = True
[154]382 
383   cs_res@gsnDraw                 = False
384   cs_res@gsnFrame                = False
[566]385   cs_res@gsnMaximize                = True
[218]386   
387   cs_res@tmXBLabelFontHeightF   = font_size
388   cs_res@tmYLLabelFontHeightF   = font_size
389   cs_res@tiXAxisFontHeightF     = font_size
390   cs_res@tiYAxisFontHeightF     = font_size
[154]391   cs_res@tmXBMode                ="Automatic"
392   cs_res@tmYLMode                ="Automatic"
[218]393 
[162]394   cs_res@cnLevelSelectionMode    = "ManualLevels"
[218]395   cs_res@lbLabelFontHeightF = font_size_legend
396   cs_res@lbLabelStride = legend_label_stride
397   
[162]398
[154]399   cs_resP = True
[983]400   cs_resP@gsnMaximize                = True
401   cs_resP@gsnPanelXWhiteSpacePercent = 4.0
402   cs_resP@gsnPanelYWhiteSpacePercent = 4.0
403   cs_resP@txFont                     = "helvetica"
404   cs_resP@txString                   = f_att@title
405   cs_resP@txFuncCode                 = "~"             
406   cs_resP@txFontHeightF              = 0.0105       
[154]407 
408   if ( mode .eq. "Fill" ) then
409      cs_res@cnFillOn             = True
410      cs_res@gsnSpreadColors      = True
[218]411      cs_res@cnFillMode           = fill_mode   
[154]412      cs_res@cnLinesOn            = False       
413      cs_res@cnLineLabelsOn       = False
414   end if
415
416   if ( mode .eq. "Both" ) then
417      cs_res@cnFillOn             = True
418      cs_res@gsnSpreadColors      = True
[218]419      cs_res@cnFillMode           = fill_mode
[154]420      cs_res@cnLinesOn            = True
421      cs_res@cnLineLabelsOn       = True
422   end if
423
[157]424   ; *********************************************
425   ; set up of start_time_step and end_time_step
426   ; *********************************************
[154]427
[827]428   t_all         = f[:]->time
429   nt            = dimsizes(t_all)
[162]430
[827]431   if (nt .EQ. 1)then
432      delta_t=t_all(nt-1)/nt
433   else
434      delta_t_array = new(nt-1,double)
435
436      do i=0,nt-2
437         delta_t_array(i) = t_all(i+1)-t_all(i)
438      end do
439
440      delta_t = min(delta_t_array)
441      delete(delta_t_array)
442   end if
443
444
[162]445   ; ****************************************************       
[157]446   ; start of time step and different types of mistakes that could be done
[162]447   ; ****************************************************
[154]448
[190]449   if (start_time_step .EQ. -1.) then           
450      start_time_step=t_all(0)/3600     
[157]451   else
[162]452      if (start_time_step .GT. t_all(nt-1)/3600)
[157]453         print(" ")
[566]454         print("'start_time_step' = "+ start_time_step +"h is greater than "+\
455               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[157]456         print(" ")
[218]457         print("Select another 'start_time_step'")
[162]458         print(" ")
[157]459         exit
460      end if
[162]461      if (start_time_step .LT. t_all(0)/3600)
[157]462         print(" ")
[566]463         print("'start_time_step' = "+ start_time_step +"h is lower than "+\
464               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]465         print(" ")
[218]466         print("Select another 'start_time_step'")
[162]467         print(" ")
[157]468         exit
469      end if
470   end if
[218]471   do i=0,nt-1   
[566]472      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
473          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[162]474         st=i
475         break
[769]476      else
477         st=0
[162]478      end if
479   end do
[769]480
[162]481       
[157]482   ; ****************************************************
483   ; end of time step and different types of mistakes that could be done
484   ; ****************************************************
485
[190]486   if (end_time_step .EQ. -1.) then             
487      end_time_step = t_all(nt-1)/3600 
488   else 
[162]489      if (end_time_step .LT. t_all(0)/3600)
[157]490         print(" ")
[566]491         print("'end_time_step' = "+end_time_step+ "h is lower than "+\
492               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
[157]493         print(" ")
[218]494         print("Select another 'end_time_step'")
[162]495         print(" ")
[157]496         exit
497      end if
[162]498      if (end_time_step .GT. t_all(nt-1)/3600)
[157]499         print(" ")
[566]500         print("'end_time_step' = "+ end_time_step +"h is greater than "+\
501               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
[162]502         print(" ")
[218]503         print("Select another 'end_time_step'") 
[162]504         print(" ")
[157]505         exit
506      end if
[162]507      if (end_time_step .LT. start_time_step/3600)
[157]508         print(" ")
[566]509         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
510               "'start_time_step' = "+start_time_step+"h")
[157]511         print(" ")
[218]512         print("Select another 'start_time_step' or 'end_time_step'")
[162]513         print(" ")
[157]514         exit
515      end if
516   end if
[162]517   do i=0,nt-1     
[566]518      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
519          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
[162]520         et=i
521         break
[769]522      else
523         et=0
[162]524      end if
525   end do
526 
527   delete(start_time_step)
528   start_time_step=round(st,3)
529   delete(end_time_step)
530   end_time_step=round(et,3)
[157]531
[162]532   print(" ")
[566]533   print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
534                t_all(start_time_step)+" s => index = "+start_time_step)
535   print("                     till "+t_all(end_time_step)/3600+" h = "  +\
536                t_all(end_time_step)+" s => index = "+end_time_step)
[162]537   print(" ")
538 
[161]539   no_time=(end_time_step-start_time_step)+1
540
[157]541   ; ****************************************************
542   ; set up variables for vector plot if required
543   ; ****************************************************   
544
545   if (vector .EQ. 1) then
[190]546      if (vec1 .EQ. "component1") then
547         print(" ")
[566]548         print("Indicate Vector 1 ('vec1') for Vector-Plot or "+\
549               "set 'vector' to 0")
[190]550         print(" ")
551         exit
[157]552      end if
[190]553      if (vec2 .EQ. "component2") then
554         print(" ")
[566]555         print("Indicate Vector 2 ('vec2') for Vector-Plot or "+\
556               "set 'vector' to 0")
[190]557         print(" ")
558         exit
559      end if           
[157]560   end if
561
562   check_vec1 = False
563   check_vec2 = False
564   check_vecp = False
565
566   ; ****************************************************
567   ; get data for plots
568   ; ****************************************************
569
[162]570   if (xyc .EQ. 1) then
571      do varn=0,dim-1
[190]572         if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then   
[218]573            x_d     = f_att->$vNam(varn)$
[162]574            xdim    = dimsizes(x_d)-1
575            delta_x = x_d(1)-x_d(0)
576            break
577         end if
578      end do
579      do varn=0,dim-1
580         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
[218]581            y_d     = f_att->$vNam(varn)$
[162]582            ydim    = dimsizes(y_d)-1
583            delta_y = y_d(1)-y_d(0)
584            break
585         end if
586      end do
587      do varn=0,dim-1
588         if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
[218]589            z_d     = f_att->$vNam(varn)$
[162]590            zdim    = dimsizes(z_d)-1
591            delta_z = 0
592            break
593         else
594            if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
[218]595               z_d     = f_att->$vNam(varn)$
[162]596               zdim    = dimsizes(z_d)-1
597               delta_z = -1.d
598               break
599            else
600               zdim    = 0
601               delta_z = -1.d
602            end if
603         end if
[357]604         if (vNam(varn) .eq. "zu1_xy") then
605            zu1    = f_att->$vNam(varn)$
606         end if
[162]607      end do
608   end if
609   if (xzc .EQ. 1) then
610      do varn=0,dim-1
611         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then
[218]612            x_d     = f_att->$vNam(varn)$
[162]613            xdim    = dimsizes(x_d)-1
614            delta_x = x_d(1)-x_d(0)
615            break
616         end if
617      end do
618      do varn=0,dim-1   
[566]619         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
620             vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
[218]621            z_d     = f_att->$vNam(varn)$
[162]622            zdim    = dimsizes(z_d)-1
623            delta_z = z_d(1)-z_d(0)
624            break
625         end if
626      end do
627      do varn=0,dim-1
628         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
[218]629            y_d     = f_att->$vNam(varn)$
[162]630            ydim    = dimsizes(y_d)-1
631            delta_y = y_d(1)-y_d(0)
632            break
633         else
634            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
[218]635               y_d     = f_att->$vNam(varn)$
[162]636               ydim    = dimsizes(y_d)-1
637               delta_y = -1.d
638               break
639            else
640               ydim    = 0
641               delta_y = -1.d
642            end if
643         end if
644      end do
645   end if
646   if (yzc .EQ. 1) then
647      do varn=0,dim-1 
648         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
[218]649            y_d     = f_att->$vNam(varn)$
[162]650            ydim    = dimsizes(y_d)-1
651            delta_y = y_d(1)-y_d(0)
652            break
653         end if
654      end do
655      do varn=0,dim-1
[566]656         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
657             vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
[218]658            z_d     = f_att->$vNam(varn)$
[162]659            zdim    = dimsizes(z_d)-1
660            delta_z = z_d(1)-z_d(0)
661            break
662         end if
663      end do
664      do varn=0,dim-1
665         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
[218]666            x_d     = f_att->$vNam(varn)$
[162]667            xdim    = dimsizes(x_d)-1
668            delta_x = x_d(1)-x_d(0)
669            break
670         else
671            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
[218]672               x_d     = f_att->$vNam(varn)$
[162]673               xdim    = dimsizes(x_d)-1
674               delta_x = -1.d
675               break
676            else
677               xdim    = 0
678               delta_x = -1.d
679            end if
680         end if
681      end do
682   end if
[195]683 
[162]684   ; ****************************************************
685   ; set up ranges of x-, y- and z-coordinates
686   ; ****************************************************         
687                   
[190]688   if (xs .EQ. -1.d) then               
689      xs = x_d(0)
[195]690      if (delta_x .EQ. -1) then
691         print(" ")
[566]692         print("You cannot choose a start value for x, "+\
693               "there are preset layers for x")
[195]694         print(" ")
695         xstart=0
696      end if
[190]697   else
[162]698      if (delta_x .EQ. -1) then
699         print(" ")
[566]700         print("You cannot choose a start value for x, "+\
701               "there are preset layers for x")
[162]702         print(" ")
703         xstart=0
[161]704      else
[162]705         if (xs .LT. 0-delta_x/2) then
706            print(" ")
[566]707            print("range start for x-coordinate = "+\
708                   xs+"m is lower than first value = "+(0-delta_x/2)+"m")
[162]709            print(" ")
710            exit
711         end if
712         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
713            if (xs .GE. xdim*delta_x) then
714               print(" ")
[566]715               print("range start for x-coordinate = "+\
716                     xs+"m is equal or greater than last value = "+\
717                     xdim*delta_x+"m")
[162]718               print(" ")
719               exit
720            end if
721         else
722            if (xs .GT. xdim*delta_x) then
723               print(" ")
[566]724               print("range start for x-coordinate = "+\
725                     xs+"m is greater than last value = "+\
726                     xdim*delta_x+"m")
[162]727               print(" ")
728               exit
729            end if
730         end if
[154]731      end if
[162]732   end if
[154]733
[190]734   if (ys .EQ. -1.d) then   
735      ys = y_d(0)
[195]736      if (delta_y .EQ. -1) then
737         print(" ")
[566]738         print("You cannot choose a start value for y, "+\
739               "there are preset layers for y")
[195]740         print(" ")
741         ystart=0
742      end if
[162]743   else
744      if (delta_y .EQ. -1) then
745         print(" ")
[566]746         print("You cannot choose a start value for y, "+\
747                "there are preset layers for y")
[162]748         print(" ")
749         ystart=0
750      else
751         if (ys .LT. 0-delta_y/2) then
752            print(" ")
[566]753            print("range start for y-coordinate = "+\
754                  ys+"m is lower than first value = "+0-delta_y/2+"m")
[162]755            print(" ")
756            exit
757         end if
758         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
759            if (ys .GE. ydim*delta_y) then
760               print(" ")
[566]761               print("range start for y-coordinate = "+\
762                     ys+"m is equal or greater than last value = "+\
763                     ydim*delta_y+"m")
[162]764               print(" ")
765               exit
766            end if
[154]767         else
[162]768            if (ys .GT. ydim*delta_y) then
[157]769               print(" ")
[566]770               print("range start for y-coordinate = "+\
771                     ys+"m is greater than last value = "+ydim*delta_y+"m")
[157]772               print(" ")
773               exit
774            end if
[162]775         end if
776      end if
777   end if
778 
[190]779   if (zs .EQ. -1) then                         
[162]780      zs = 0
[195]781      if (delta_z .EQ. -1) then
782         print(" ")
[566]783         print("You cannot choose a start value for z, "+\
784               "there are preset layers for z")
[195]785         print(" ")
786      end if
[162]787   else
788      if (delta_z .EQ. -1) then
789         print(" ")
[566]790         print("You cannot choose a start value for z, "+\
791               "there are preset layers for z")
[162]792         print(" ")
[358]793         zs = 0
[162]794      else
795         if (zs .LT. 0) then
796            print(" ")
[566]797            print("range start for z-coordinate = "+\
798                  zs+" is lower than first gridpoint = 0")
[162]799            print(" ")
800            exit
801         end if
802         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
803            if (zs .GE. zdim) then
804               print(" ")
[566]805               print("range start for z-coordinate = "+\
806                     zs+" is equal or greater than last gridpoint = "+zdim)
[162]807               print(" ")
808               exit
809            end if
810         else
811            if (zs .GT. zdim) then
812               print(" ")
[566]813               print("range start for z-coordinate = "+\
814                     zs+" is greater than last gridpoint = "+zdim)
[162]815               print(" ")
816               exit
817            end if
818         end if
819      end if
820   end if 
[157]821
[190]822   if (xe .EQ. -1) then   
823      xe = x_d(xdim)
[195]824      if (delta_x .EQ. -1) then
825         print(" ")
[566]826         print("You cannot choose an end value for x, "+\
827               "there are preset layers for x")
[195]828         print(" ")
829         xend=xdim
830      end if
[162]831   else
832      if (delta_x .EQ. -1) then
833         print(" ")
[566]834         print("You cannot choose an end value for x, "+\
835               "there are preset layers for x")
[162]836         print(" ")
837         xend=xdim
838      else
839         if (xe .GT. xdim*delta_x) then
840            print(" ")
[566]841            print("range end for x-coordinate = "+\
842                  xe+"m is greater than last value = "+xdim*delta_x+"m")
[162]843            print(" ")
844            exit
845         end if
846         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
847            if (xe .LE. 0-delta_x/2)
848               print(" ")
[566]849               print("range end for x-coordinate = "+\
850                     xe+"m is equal or lower than first value = "+\
851                     (0-delta_x/2)+"m")
[162]852               print(" ")
853               exit
854            end if
855            if (xe .LE. xs) then
856               print(" ")
[566]857               print("range end for x-coordinate = "+\
858                     xe+"m is equal or lower than start range = "+xs+"m")
[162]859               print(" ")
860               exit
861            end if
[190]862            if (xe .EQ. xs+1) then
[162]863               print(" ")
[566]864               print("range end for x-coordinate = "+\
865                     xe+"m must be at least two more gridpoints "+\
866                     "greater than start range = "+xs+"m")
[162]867               print(" ")
868               exit
869            end if
[154]870         else
[162]871            if (xe .LT. 0-delta_x/2)
[154]872               print(" ")
[566]873               print("range end for x-coordinate = "+\
874                     xe+"m is lower than first value = "+(0-delta_x/2)+"m")
[154]875               print(" ")
[162]876               exit
[154]877            end if
[162]878            if (xe .LT. xs) then
[154]879               print(" ")
[566]880               print("range end for x-coordinate = "+\
881                     xe+"m is lower than start range = "+xs+"m")
[154]882               print(" ")
883               exit
884            end if
[162]885         end if
886      end if               
887   end if
[190]888   
889   if (ye .EQ. -1) then 
890      ye = y_d(ydim)
[195]891      if (delta_y .EQ. -1) then
892         print(" ")
[566]893         print("You cannot choose an end value for y, "+\
894               "there are preset layers for y")
[195]895         print(" ")
896         yend=ydim
897      end if
[162]898   else
899      if (delta_y .EQ. -1) then
900         print(" ")
[566]901         print("You cannot choose an end value for y, "+\
902               "there are preset layers for y")
[162]903         print(" ")
904         yend=ydim
905      else
906         if (ye .GT. ydim*delta_y) then
907            print(" ")
[566]908            print("range end for y-coordinate = "+\
909                  ye+"m is greater than last value = "+ydim*delta_y+"m")
[162]910            print(" ")
911            exit
912         end if
913         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
914            if (ye .LE. 0-delta_y/2)
915               print(" ")
[566]916               print("range end for y-coordinate = "+\
917                     ye+"m is equal or lower than first value = "+\
918                     (0-delta_y/2)+"m")
[162]919               print(" ")
920               exit
921            end if
922            if (ye .LE. ys) then
923               print(" ")
[566]924               print("range end for y-coordinate = "+\
925                    ye+"m is equal or lower than start range = "+ys+"m")
[162]926               print(" ")
927               exit
928            end if
929            if (ye .EQ. ys+1) then
930               print(" ")
[566]931               print("range end for y-coordinate = "+\
932                     ye+"m must be at least two more gridpoints greater "+\
933                     "than start range = "+ys+"m")
[162]934               print(" ")
935               exit
936            end if
[154]937         else
[162]938            if (ye .LT. 0-delta_y/2)
[154]939               print(" ")
[566]940               print("range end for y-coordinate = "+\
941                     ye+"m is lower than first value = "+(0-delta_y/2)+"m")
[162]942               print(" ")
943               exit
944            end if
945            if (ye .LT. ys) then
946               print(" ")
[566]947               print("range end for y-coordinate = "+\
948                     ye+"m is lower than start range = "+ys+"m")
[162]949               print(" ")
950               exit
951            end if
952         end if
953      end if
954   end if
955 
[190]956   if (ze .EQ. -1) then 
[162]957      ze = zdim
[195]958      if (delta_z .EQ. -1) then
959         print(" ")
[566]960         print("You cannot choose an end value for z, "+\
961               "there are preset layers for z")
[195]962         print(" ")
963         ze = zdim
964      end if
[162]965   else
966      if (delta_z .EQ. -1) then
967         print(" ")
[566]968         print("You cannot choose an end value for z, "+\
969               "there are preset layers for z")
[162]970         print(" ")
971         ze = zdim
972      else
973         if (ze .GT. zdim) then
974            print(" ")
[566]975            print("range end for z-coordinate = "+\
976                  ze+" is greater than last gridpoint = "+zdim)
[162]977            print(" ")
978            exit
979         end if
980         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
981            if (ze .LE. 0)
982               print(" ")
[566]983               print("range end for z-coordinate = "+\
984                     ze+" is equal or lower than first gridpoint = 0")
[162]985               print(" ")
986               exit
987            end if
988            if (ze .LE. zs) then
989               print(" ")
[566]990               print("range end for z-coordinate = "+\
991                     ze+" is equal or lower than start range = "+zs)
[162]992               print(" ")
993               exit
994            end if
995            if (ze .EQ. zs+1) then
996               print(" ")
[566]997               print("range end for z-coordinate = "+\
998                     ze+" must be at least two more gridpoints "+\
999                     "greater than start range = "+zs)
[162]1000               print(" ")
1001               exit
1002            end if
1003         else
1004            if (ze .LT. 0)
1005               print(" ")
[566]1006               print("range end for z-coordinate = "+\
1007                     ze+" is lower than first gridpoint = 0")
[162]1008               print(" ")
1009               exit
1010            end if
1011            if (ze .LT. zs) then
1012               print(" ")
[566]1013               print("range end for z-coordinate = "+\
1014                     ze+" is lower than start range = "+zs)
[162]1015               print(" ")
1016               exit
1017            end if
1018         end if
1019      end if
1020   end if
[157]1021
[162]1022   if (delta_x .NE. -1) then
1023      do i=0,xdim   
1024         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
1025            xstart=i
1026            break
1027         end if
1028      end do
1029      do i=0,xdim   
1030         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
1031            xend=i
1032            break
1033         end if
1034      end do
1035   end if
1036   if (delta_y .NE. -1) then
[190]1037      do i=0,ydim   
[162]1038         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
1039            ystart=i
1040            break
1041         end if
1042      end do
[190]1043      do i=0,ydim   
[162]1044         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
1045            yend=i
1046            break
1047         end if
1048      end do
1049   end if
[286]1050
[190]1051   if( shape .eq. 1 ) then
[195]1052      if (xyc .EQ. 1)then
1053         cs_res@vpWidthF    = (xe-xs)/(ye-ys)       
1054         cs_res@vpHeightF   = 1
[218]1055         vecres@vpWidthF    = (xe-xs)/(ye-ys)       
1056         vecres@vpHeightF   = 1
1057         if (xe-xs .GT. ye-ys)then
[286]1058            if (no_rows .LE. no_columns)then
1059               cs_res@gsnPaperOrientation     = "landscape"
1060               vecres@gsnPaperOrientation     = "landscape"
1061               cs_res@lbOrientation = "Horizontal"
1062            else
1063               cs_res@gsnPaperOrientation     = "portrait"
1064               vecres@gsnPaperOrientation     = "portrait"
1065               cs_res@lbOrientation = "Vertical"
1066            end if 
[218]1067         else
[286]1068            if (no_rows .GE. no_columns)then
1069               cs_res@gsnPaperOrientation     = "portrait"
1070               vecres@gsnPaperOrientation     = "portrait"
1071               cs_res@lbOrientation = "Vertical"
1072            else
1073               cs_res@gsnPaperOrientation     = "landscape"
1074               vecres@gsnPaperOrientation     = "landscape"
1075               cs_res@lbOrientation = "Horizontal"
1076            end if
[218]1077         end if
[195]1078      end if
1079      if (xzc .EQ. 1)then
1080         cs_res@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
1081         cs_res@vpHeightF   = 1
[218]1082         vecres@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))       
1083         vecres@vpHeightF   = 1
1084         if (xe-xs .GT. (delta_x*(ze-zs)))then
[286]1085            if (no_rows .LE. no_columns)then
1086               cs_res@gsnPaperOrientation     = "landscape"
1087               vecres@gsnPaperOrientation     = "landscape"
1088               cs_res@lbOrientation = "Horizontal"
1089            else
1090               cs_res@gsnPaperOrientation     = "portrait"
1091               vecres@gsnPaperOrientation     = "portrait"
1092               cs_res@lbOrientation = "Vertical"
1093            end if 
[218]1094         else
[286]1095            if (no_rows .GE. no_columns)then
1096               cs_res@gsnPaperOrientation     = "portrait"
1097               vecres@gsnPaperOrientation     = "portrait"
1098               cs_res@lbOrientation = "Vertical"
1099            else
1100               cs_res@gsnPaperOrientation     = "landscape"
1101               vecres@gsnPaperOrientation     = "landscape"
1102               cs_res@lbOrientation = "Horizontal"
1103            end if
[218]1104         end if
[195]1105      end if
1106      if (yzc .EQ. 1)then
1107         cs_res@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
1108         cs_res@vpHeightF   = 1
[218]1109         vecres@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))       
1110         vecres@vpHeightF   = 1
1111         if (ye-ys .GT. (delta_y*(ze-zs)))then
[286]1112            if (no_rows .LE. no_columns)then
1113               cs_res@gsnPaperOrientation     = "landscape"
1114               vecres@gsnPaperOrientation     = "landscape"
1115               cs_res@lbOrientation = "Horizontal"
1116            else
1117               cs_res@gsnPaperOrientation     = "portrait"
1118               vecres@gsnPaperOrientation     = "portrait"
1119               cs_res@lbOrientation = "Vertical"
1120            end if 
[218]1121         else
[286]1122            if (no_rows .GE. no_columns)then
1123               cs_res@gsnPaperOrientation     = "portrait"
1124               vecres@gsnPaperOrientation     = "portrait"
1125               cs_res@lbOrientation = "Vertical"
1126            else
1127               cs_res@gsnPaperOrientation     = "landscape"
1128               vecres@gsnPaperOrientation     = "landscape"
1129               cs_res@lbOrientation = "Horizontal"
1130            end if
[218]1131         end if
[195]1132      end if
[190]1133   end if
[218]1134
[162]1135   delete(xs)
1136   delete(xe)   
1137   delete(ys)
1138   delete(ye)
1139
1140   xs=xstart
1141   xe=xend
1142   ys=ystart
1143   ye=yend
1144
[218]1145   if (xyc .EQ. 1) then
[329]1146      d= (ye-ys+1)/(major_ticks_y-1)
[218]1147      e=(xe-xs+1)/(major_ticks_x-1)
1148      array_yl       =new(major_ticks_y,integer)
1149      array_yl_labels=new(major_ticks_y,double)
1150      array_minor_yl =new((major_ticks_y-1)*4,double)
1151      array_xb       =new(major_ticks_x,integer)
1152      array_xb_labels=new(major_ticks_x,double)
1153      array_minor_xb =new((major_ticks_x-1)*4,double)
1154      array_yl(0)=ys
1155      array_xb(0)=xs
1156      array_yl_labels(0)=y_d(ys)
1157      array_xb_labels(0)=x_d(xs)
1158      do ar=1,major_ticks_y-1
1159         array_yl(ar)=d*(ar-1)+d-1
1160         array_yl_labels(ar) = y_d(array_yl(ar))
1161         if (ar .GT. 0)
1162            do min_ar=0,3
[566]1163                array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+\
1164                         (y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1)
[218]1165            end do
1166         end if 
1167      end do
1168      do br=1,major_ticks_x-1
1169         array_xb(br)=e*(br-1)+e-1
1170         array_xb_labels(br) = x_d(array_xb(br))
1171         if (br .GT. 0)
1172            do min_br=0,3
[566]1173               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1174                         (x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
[218]1175            end do
1176         end if
1177      end do
1178   end if
1179 
1180   if (xzc .EQ. 1) then
1181      d=(ze-zs+1)/(major_ticks_y-1)
1182      e=(xe-xs+1)/(major_ticks_x-1)
1183      array_yl       =new(major_ticks_y,integer)
1184      array_yl_labels=new(major_ticks_y,double)
1185      array_minor_yl =new((major_ticks_y-1)*4,double)
1186      array_xb       =new(major_ticks_x,integer)
1187      array_xb_labels=new(major_ticks_x,double)
1188      array_minor_xb =new((major_ticks_x-1)*4,double)
1189      array_yl(0)=zs
1190      array_xb(0)=xs
1191      array_yl_labels(0)=z_d(zs)
1192      array_xb_labels(0)=x_d(xs)
1193      do ar=1,major_ticks_y-1
1194         array_yl(ar)=d*(ar-1)+d-1
1195         array_yl_labels(ar) = z_d(array_yl(ar))
1196         if (ar .GT. 0)
1197            do min_ar=0,3
[566]1198               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1199                          (z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
[218]1200            end do
1201         end if 
1202      end do
1203      do br=1,major_ticks_x-1
1204         array_xb(br)=e*(br-1)+e-1
1205         array_xb_labels(br) = x_d(array_xb(br))
1206         if (br .GT. 0)
1207            do min_br=0,3
[566]1208               array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1209                         (x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
[218]1210            end do
1211         end if
1212      end do
1213   end if
1214
1215   if (yzc .EQ. 1) then
1216      d=(ze-zs+1)/(major_ticks_y-1)
1217      e=(ye-ys+1)/(major_ticks_x-1)
1218      array_yl       =new(major_ticks_y,integer)
1219      array_yl_labels=new(major_ticks_y,double)
1220      array_minor_yl =new((major_ticks_y-1)*4,double)
1221      array_xb       =new(major_ticks_x,integer)
1222      array_xb_labels=new(major_ticks_x,double)
1223      array_minor_xb =new((major_ticks_x-1)*4,double)
1224      array_yl(0)=zs
1225      array_xb(0)=ys
1226      array_yl_labels(0)=z_d(zs)
1227      array_xb_labels(0)=y_d(ys)
1228      do ar=1,major_ticks_y-1
1229         array_yl(ar)=d*(ar-1)+d-1
[329]1230         array_yl_labels(ar) = z_d(array_yl(ar))
[218]1231         if (ar .GT. 0)
1232            do min_ar=0,3
[566]1233               array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1234                         (z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
[218]1235            end do
1236         end if 
1237      end do
1238      do br=1,major_ticks_x-1
1239         array_xb(br)=e*(br-1)+e-1
[329]1240         array_xb_labels(br) = y_d(array_xb(br))
[218]1241         if (br .GT. 0)
1242            do min_br=0,3
[566]1243               array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+\
1244                         (y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1)
[218]1245            end do
1246         end if
1247      end do
1248   end if
[329]1249 
[218]1250   if (axes_explicit .EQ. 1)then
1251      cs_res@tmYLMode = "Explicit"
1252      cs_res@tmXBMode = "Explicit"
[329]1253      if (xyc .EQ. 1)then
1254         cs_res@tmYLValues = y_d(array_yl)
1255         cs_res@tmXBValues = x_d(array_xb)     
[218]1256         cs_res@tmYLLabels = array_yl_labels/norm_y
1257         cs_res@tmXBLabels = array_xb_labels/norm_x
1258         if (norm_x .NE. 1.)then
[534]1259            cs_res@tiXAxisString = "x ("+unit_x+")"
[218]1260         else
[534]1261            cs_res@tiXAxisString = "x (m)"
[218]1262         end if
1263         if (norm_y .NE. 1.)then
[534]1264            cs_res@tiYAxisString = "y ("+unit_y+")"
[218]1265         else
[534]1266            cs_res@tiYAxisString = "y (m)"
[218]1267         end if   
1268      end if
1269      if (xzc .EQ. 1)then
[329]1270         cs_res@tmYLValues = z_d(array_yl)
1271         cs_res@tmXBValues = x_d(array_xb) 
[218]1272         cs_res@tmYLLabels = array_yl_labels/norm_z
1273         cs_res@tmXBLabels = array_xb_labels/norm_x
1274         if (norm_x .NE. 1.)then
[534]1275            cs_res@tiXAxisString = "x ("+unit_x+")"
[218]1276         else
[534]1277            cs_res@tiXAxisString = "x (m)"
[218]1278         end if
1279         if (norm_z .NE. 1.)then
[534]1280            cs_res@tiYAxisString = "z ("+unit_z+")"
[218]1281         else
[534]1282            cs_res@tiYAxisString = "z (m)"
[218]1283         end if
1284      end if
1285      if (yzc .EQ. 1)then
[329]1286         cs_res@tmYLValues = z_d(array_yl)
1287         cs_res@tmXBValues = y_d(array_xb) 
[218]1288         cs_res@tmYLLabels = array_yl_labels/norm_z
1289         cs_res@tmXBLabels = array_xb_labels/norm_y
1290         if (norm_y .NE. 1.)then
[534]1291            cs_res@tiXAxisString = "y ("+unit_y+")"
[218]1292         else
[534]1293            cs_res@tiXAxisString = "y (m)"
[218]1294         end if
1295         if (norm_z .NE. 1.)then
[534]1296            cs_res@tiYAxisString = "z ("+unit_z+")"
[218]1297         else
[534]1298            cs_res@tiYAxisString = "z (m)"
[218]1299         end if
1300      end if
1301      cs_res@tmYLMinorValues = array_minor_yl
1302      cs_res@tmXBMinorValues = array_minor_xb
1303   else
1304      if (axes_explicit .EQ. 0)then 
1305         cs_res@tmYLMinorPerMajor = 4
1306         cs_res@tmXBMinorPerMajor = 4
1307      else
1308         print(" ")
1309         print("'axes_explicit' is invalid and set to 0")
1310         print(" ")
1311         axes_explicit = 0
1312         cs_res@tmYLMinorPerMajor = 4
1313         cs_res@tmXBMinorPerMajor = 4
1314      end if
1315      if (xyc .EQ. 1)then
[534]1316         cs_res@tiXAxisString = "x (m)"
1317         cs_res@tiYAxisString = "y (m)"
[218]1318      end if
1319      if (xzc .EQ. 1)then
[534]1320         cs_res@tiXAxisString = "x (m)"
1321         cs_res@tiYAxisString = "z (m)"
[218]1322      end if
1323      if (yzc .EQ. 1)then
[534]1324         cs_res@tiXAxisString = "y (m)"
1325         cs_res@tiYAxisString = "z (m)"
[218]1326      end if
1327   end if
1328
[162]1329   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1330      if (xe .EQ. xs+1) then
1331         print(" ")
[566]1332         print("range end for x-coordinate="+xe*delta_x+\
1333               "m must be at least two")
1334         print("more gridpoints(="+2*delta_x+"m) greater than start range="+\
1335               xs*delta_x+"m)")
[162]1336         print(" ")
1337         exit
1338      end if
1339   end if
1340   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1341      if (ye .EQ. ys+1) then
1342         print(" ")
[566]1343         print("range end for y-coordinate="+ye*delta_y+\
1344               "m must be at least two")
1345         print("more gridpoints(="+2*delta_y+"m greater than start range="+\
1346               ys*delta_y+"m)")
[162]1347         print(" ")
1348         exit
1349      end if
1350   end if
1351   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1352      if (ze .EQ. zs+1) then
1353         print(" ")
[566]1354         print("range end for z-coordinate="+ze+" must be at least two")
[190]1355         print("more gridpoints greater than start range="+zs+")")
[162]1356         print(" ")
1357         exit
1358      end if
1359   end if
1360 
[161]1361   if (xyc .EQ. 1) then
1362      no_layer = (ze-zs)+1
[585]1363      if (check_vType) then
1364         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1365      else
1366         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1367      end if
[161]1368   end if
1369   if (xzc .EQ. 1) then
1370      no_layer = (ye-ys)+1
[585]1371      if (check_vType) then
1372         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1373      else
1374         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1375      end if
[161]1376   end if
1377   if (yzc .EQ. 1) then
1378      no_layer = (xe-xs)+1
[585]1379      if (check_vType) then
1380         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1381      else
1382         data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1383      end if
[161]1384   end if
1385
[585]1386   if (check_vType) then
1387      MinVal = new(dim,double)
1388      MaxVal = new(dim,double)
1389   else
1390      MinVal = new(dim,float)
1391      MaxVal = new(dim,float)
1392   end if
[218]1393   unit   = new(dim,string)
[157]1394
1395   ; ****************************************************
1396   ; define inner and outer loops depending on "sort"
1397   ; ****************************************************       
1398   
1399   if ( xyc .eq. 1 ) then
1400      lays = zs
1401      laye = ze
1402   end if
1403   if ( xzc .eq. 1 ) then
1404      lays = ys
1405      laye = ye
1406   end if
1407   if ( yzc .eq. 1) then
1408      lays = xs
1409      laye = xe
1410   end if
1411
1412   if ( sort .eq. "time" ) then
1413      lis = start_time_step
1414      lie = end_time_step
[358]1415      los = lays
1416      loe = laye
[157]1417   else
[358]1418      lis = lays
1419      lie = laye
[157]1420      los = start_time_step
1421      loe = end_time_step
1422   end if
1423
[758]1424   check  = True
1425   v1     = 0
1426   v2     = 0
1427   no_var = 0
1428   n      = 0
1429   no_zu1 = 0
[769]1430   no_topo= 0
[157]1431
[566]1432   ;****************************************************
1433   ; Preparation of vector plots
1434   ; ****************************************************       
1435
1436   do varn=dim-1,0,1
1437
1438     if (vector .EQ. 1) then   
1439        check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1440        if (check_vec1) then
1441            temp = f[:]->$vNam(varn)$
1442            data_att = f_att->$vNam(varn)$
1443            vect1=temp(:,zs:ze,ys:ye,xs:xe)
1444            delete(temp)
1445         end if
1446
1447         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1448         if (check_vec2) then
1449            temp = f[:]->$vNam(varn)$
1450            data_att = f_att->$vNam(varn)$
1451            vect2=temp(:,zs:ze,ys:ye,xs:xe)
1452            delete(temp)
1453         end if
[162]1454   
[566]1455         if (","+vNam(varn)+"," .EQ. vec1) then
1456            v1=v1+1
1457         end if
1458         if (","+vNam(varn)+"," .EQ. vec2) then
1459             v2=v2+1
1460          end if
1461     end if
1462
1463   end do
1464
1465
1466
1467   do varn=dim-1,0,1     
[157]1468           
[566]1469      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1470           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1471           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1472           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1473           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1474           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1475           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1476           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1477           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1478           vNam(varn) .eq. "ind_x_yz") then
[157]1479         check = False
[161]1480      else
1481         check = True
1482      end if 
[157]1483
[357]1484      if (var .NE. "all") then
1485         check = isStrSubset( var,","+vNam(varn)+"," )
1486      end if
1487
[162]1488      if(check) then
[194]1489     
[161]1490         no_var=no_var+1
1491
[162]1492         if (xyc .EQ. 1) then
[218]1493            temp = f[:]->$vNam(varn)$
[769]1494            if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1495               dummy=0
1496            else
1497               data_att = f_att->$vNam(varn)$
1498            end if
[514]1499            if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy"     \
1500              .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1501              .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy"   \
[526]1502              .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy"    \
1503              .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \
1504              .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\
1505              .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy"   \
1506              .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then
[566]1507              ;these variables depend on zu1_xy and that's why they have
1508              ;only one z-layer
[329]1509              data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
[357]1510              no_zu1=no_zu1+1
[329]1511            else
[769]1512              if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1513                  ;these variables depend on x and y
1514                  data(varn,0,0,:,:)=doubletofloat(temp(ys:ye,xs:xe))
1515                  no_topo=no_topo+1
1516               else
1517                  data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1518               end if
[329]1519            end if
1520            delete(temp)               
[162]1521         end if
1522         if ( xzc .eq. 1 ) then
[218]1523            temp = f[:]->$vNam(varn)$
1524            data_att = f_att->$vNam(varn)$
1525            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[357]1526            delete(temp)
[162]1527         end if
1528         if ( yzc .eq. 1) then
[218]1529            temp = f[:]->$vNam(varn)$
1530            data_att = f_att->$vNam(varn)$
1531            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
[357]1532            delete(temp)
[162]1533         end if
[154]1534
[157]1535         data!0 = "var"
1536         data!1 = "t"
1537         data!2 = "z"
1538         data!3 = "y"
[162]1539         data!4 = "x" 
[154]1540
[566]1541         MinVal(varn) = min(data(varn,start_time_step:end_time_step,\
1542                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1543         MaxVal(varn) = max(data(varn,start_time_step:end_time_step,\
1544                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
[218]1545         
[769]1546         if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1547            unit(varn) = "meters"
1548          else
1549            unit(varn) = data_att@units
1550            delete(data_att)
1551          end if
[162]1552
[157]1553      end if
[154]1554
[190]1555   end do
[154]1556
[357]1557   if (no_var .EQ. 0) then
1558      print(" ")
1559      print("The variables var='"+var+"' do not exist on your input file;")
1560      print("be sure to have one comma before and after each variable")
1561      print(" ")
1562      exit
1563   end if
1564
[190]1565   var_input=new(no_var,string)
1566   numb=0
[194]1567   no_var=0
1568   
[190]1569   do varn=dim-1,0,1   
1570           
[566]1571      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1572           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1573           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1574           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1575           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1576           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1577           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1578           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1579           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1580           vNam(varn) .eq. "ind_x_yz") then
[190]1581         check = False
1582      else
1583         check = True
1584      end if 
1585
[357]1586      if (var .NE. "all") then
1587         check = isStrSubset( var,","+vNam(varn)+"," )
1588      end if
1589
[190]1590      if(check) then     
1591         var_input(numb)=vNam(varn)
1592         numb=numb+1     
1593      end if
1594     
[194]1595      if (check) then
1596         no_var = no_var+1
1597      end if
1598     
[161]1599   end do
[154]1600
[190]1601   if (no_var .EQ. 0) then
[162]1602      print(" ")
[190]1603      print("The variables var='"+var+"' do not exist on your input file;")
[357]1604      print("be sure to have one comma before and after each variable")
[162]1605      print(" ")
1606      exit
1607   end if
1608
[157]1609   if (vector .EQ. 1) then
1610      if (v1 .EQ. 0)then
1611         print(" ")
[566]1612         print("Component 1 for the vector-plot ('vec1') must be one "+\
1613               "of the variables on the input file:")
[190]1614         print("- "+var_input)
[357]1615         print("be sure to have one comma before and after the variable")
[157]1616         print(" ")
1617         exit
1618      end if
[154]1619
[157]1620      if (v2 .EQ. 0)then
1621         print(" ")
[566]1622         print("Component 2 for the vector-plot ('vec2') must be one "+\
1623               "of the variables on the input file:")
[190]1624         print("- "+var_input)
[357]1625         print("be sure to have one comma before and after the variable")
[157]1626         print(" ")
1627         exit
1628      end if
1629   end if
[154]1630
[157]1631   ; ***************************************************
[161]1632   ; open workstation(s)
1633   ; ***************************************************
1634
[532]1635   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1636      format_out@wkPaperSize = "A4"
1637   end if
1638   if ( format_out .EQ. "png" ) then
1639      format_out@wkWidth  = 1000
1640      format_out@wkHeight = 1000
1641   end if
1642
[161]1643   wks_ps  = gsn_open_wks(format_out,file_out)
1644   gsn_define_colormap(wks_ps,"rainbow+white")
1645
[218]1646   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[758]1647      plot=new((/no_time*no_layer/),graphic)
[218]1648   else
[769]1649      plot=new((/no_time*no_layer*(no_var-no_topo) + no_topo - \
1650                       no_time*(no_layer-1)*no_zu1/),graphic)
[218]1651   end if
[357]1652   dim_plot=dimsizes(plot)
[218]1653
[161]1654   page = 0
1655   n=0
1656   print(" ")
[162]1657   print("Plots sorted by '"+sort+"'")
1658   print(" ")
[161]1659
1660   ; ***************************************************
[157]1661   ; create plots
1662   ; ***************************************************
1663
[357]1664
[190]1665   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[162]1666      do lo = los, loe                                         
1667         do li = lis, lie
[218]1668            if (xyc .EQ. 1)then
1669               if (sort .EQ. "time")then
[358]1670                  if(z_d(lo) .eq. -1.d) then
1671                    level = "z-average"
1672                  else
1673                    level = "z=" + z_d(lo) + "m"
1674                  end if
[218]1675               else
[358]1676                  if(z_d(li) .eq. -1.d) then
1677                    level = "z-average"
1678                  else
1679                    level = "z=" + z_d(li) + "m"
1680                  end if
[218]1681               end if
1682            end if
1683            if (xzc .EQ. 1)then
1684               if (sort .EQ. "time")then
[358]1685                 if(y_d(lo) .eq. -1.d) then
1686                   level = "y-average"
1687                 else
1688                   level = "y=" + y_d(lo) + "m"
1689                 end if
[218]1690               else
[358]1691                  if(y_d(li) .eq. -1.d) then
1692                    level = "y-average"
1693                  else
1694                    level = "y=" + y_d(li) + "m"
1695                  end if
[218]1696               end if
1697            end if
1698            if (yzc .EQ. 1)then
1699               if (sort .EQ. "time")then
[358]1700                 if(x_d(lo) .eq. -1.d) then
1701                    level = "x-average"
1702                 else
1703                    level = "x=" + x_d(lo) + "m"
1704                 end if
[218]1705               else
[358]1706                 if(x_d(li) .eq. -1.d) then
1707                    level = "x-average"
1708                 else
1709                    level = "x=" + x_d(li) + "m"
1710                 end if
[218]1711               end if
1712            end if               
[162]1713            vecres                  = True            ; vector only resources
1714            vecres@gsnDraw          = False           ; don't draw
1715            vecres@gsnFrame         = False           ; don't advance frame
1716            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1717            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
[566]1718            vecres@vcRefLengthF     = 0.05            ; define length of
1719                                                      ; vec ref
[162]1720            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
[218]1721            vecres@tmXBLabelFontHeightF   = font_size
1722            vecres@tmYLLabelFontHeightF   = font_size
1723            vecres@tiXAxisFontHeightF     = font_size
1724            vecres@tiYAxisFontHeightF     = font_size
1725            if (sort .EQ. "time")then
[566]1726               vecres@gsnLeftString = "t=" + \
1727                             decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
[218]1728            else
[566]1729               vecres@gsnLeftString = "t=" + \
1730                             decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
[218]1731            end if
[529]1732            if (xyc .EQ. 1) then 
[534]1733               vecres@tiXAxisString = "x (m)"
1734               vecres@tiYAxisString = "y (m)"   
[566]1735               if (sort .EQ. "time")then
1736                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
1737                                               vect2(li,lo-los,:,:),vecres)
[218]1738               else
[566]1739                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1740                                                vect2(lo,li-lis,:,:),vecres)
[218]1741               end if
1742            end if
1743            if (xzc .EQ. 1) then
[534]1744               vecres@tiXAxisString = "x (m)"
1745               vecres@tiYAxisString = "z (m)"
[218]1746               if (sort .EQ. "time")then
[566]1747                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1748                                               vect2(li,:,lo-los,:),vecres)
[218]1749               else
[566]1750                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
1751                                                vect2(lo,:,li-lis,:),vecres)
[218]1752               end if
1753            end if
1754            if (yzc .EQ. 1) then
[534]1755               vecres@tiXAxisString = "y (m)"
1756               vecres@tiYAxisString = "z (m)"
[218]1757               if (sort .EQ. "time")then
[566]1758                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
1759                                                vect2(li,:,:,lo-los),vecres)
[218]1760               else
[566]1761                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
1762                                                vect2(lo,:,:,li-lis),vecres)
[218]1763               end if
1764            end if
[162]1765            n=n+1
1766         end do
1767      end do
1768   end if
1769 
[161]1770   do varn=dim-1,0,1
[157]1771
[529]1772      if (vector .EQ. 1 ) then   
[566]1773         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
[157]1774      end if
[358]1775     
[566]1776      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1777           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1778           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1779           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1780           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1781           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1782           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1783           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1784           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1785           vNam(varn) .eq. "ind_x_yz") then
[161]1786         check = False
1787      else
1788         check = True
1789      end if   
[190]1790 
1791      if (var .NE. "all") then
[161]1792         check = isStrSubset( var,","+vNam(varn)+"," )
[157]1793      end if
1794   
1795      if(check) then
[218]1796
[162]1797         space=(MaxVal(varn)-MinVal(varn))/24
1798 
1799         cs_res@cnMinLevelValF = MinVal(varn)
1800         cs_res@cnMaxLevelValF = MaxVal(varn)
1801
1802         cs_res@cnLevelSpacingF = space
1803     
[157]1804         ; ****************************************************
1805         ; loops over time and layer
1806         ; ****************************************************
1807
[162]1808         
[161]1809         do lo = los, loe                                       
[162]1810            do li = lis, lie
1811
[157]1812               ; ****************************************************
[154]1813               ; xy cross section
[157]1814               ; ****************************************************
[154]1815
[157]1816               if ( xyc .eq. 1 ) then
[218]1817         
[157]1818                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[218]1819                  cs_res@gsnRightString = unit(varn)
[162]1820                 
[154]1821                  if ( sort .eq. "time" ) then
[358]1822                     if ( z_d(lo) .eq. -1)then
[194]1823                        if (delta_z .EQ. -1) then
1824                           level = "-average"
1825                        else
[358]1826                           level = "=" + z_d(lo) + "m"
[194]1827                        end if
[154]1828                     else
[358]1829                        level = "=" + z_d(lo) + "m"
[154]1830                     end if
[357]1831
[566]1832                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1833                        vNam(varn) .eq. "pras_xy"  .or. \
1834                        vNam(varn) .eq. "prrs_xy"  .or. \
1835                        vNam(varn) .eq. "qsws_xy"  .or. \
1836                        vNam(varn) .eq. "shfs_xy"  .or. \
1837                        vNam(varn) .eq. "ts_xy"    .or. \
1838                        vNam(varn) .eq. "us_xy"    .or. \
1839                        vNam(varn) .eq. "z0s_xy"   .or. \
1840                        vNam(varn) .eq. "lwp*_xy"  .or. \
1841                        vNam(varn) .eq. "pra*_xy"  .or. \
1842                        vNam(varn) .eq. "prr*_xy"  .or. \
1843                        vNam(varn) .eq. "qsws*_xy" .or. \
1844                        vNam(varn) .eq. "shf*_xy"  .or. \
1845                        vNam(varn) .eq. "t*_xy"    .or. \
1846                        vNam(varn) .eq. "u*_xy"    .or. \
1847                        vNam(varn) .eq. "z0*_xy") then
[357]1848                         loe = 0
1849                         level = "=" + zu1(0) + "m"
[769]1850                      else
1851                         loe = laye
[357]1852                     end if
[769]1853
1854                     if(vNam(varn) .eq. "zwwi"  .or. \
1855                        vNam(varn) .eq. "zusi") then
1856                         loe = 0
1857                         los = 0
1858                         lie  = 0
1859                         lis = 0
1860                         level = ""
1861                      else
1862                         lis = start_time_step
1863                         lie = end_time_step
1864                         los = lays
1865                         loe = laye
[770]1866                     end if                 
1867
[769]1868                     if(vNam(varn) .eq. "zwwi"  .or. \
1869                        vNam(varn) .eq. "zusi") then
1870                         cs_res@gsnCenterString = ""
1871                     else
1872                         cs_res@gsnCenterString = "t=" + \
1873                          decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1874                     end if
[758]1875
1876                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1877                        ;nothing to be done
1878                     else
1879                         plot(n) = gsn_csm_contour(wks_ps,\
[566]1880                                           data(varn,li,lo-los,:,:),cs_res)
[758]1881                     end if
1882
1883
[162]1884                     if (vector .EQ. 1 .AND. check_vecp) then
[566]1885                        vecres                 = True           ; vector only
1886                                                                ; resources
1887                        vecres@gsnDraw         = False          ; don't draw
1888                        vecres@gsnFrame        = False          ; don't
1889                                                                ; advance frame
1890                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly
1891                                                                ; vectors
1892                        vecres@vcRefMagnitudeF = ref_mag        ; define
1893                                                                ; vector ref
1894                                                                ; mag
1895                        vecres@vcRefLengthF    = 0.05           ; define
1896                                                                ; length of
1897                                                                ; vec ref 
1898                        vecres@gsnRightString  = " "
1899                        vecres@gsnLeftString   = " "
1900                        vecres@tiXAxisString   = " "   
1901                        plot_vec=gsn_csm_vector(wks_ps,\
1902                              vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
[162]1903                        overlay(plot(n), plot_vec)
[161]1904                     end if                         
[154]1905                  end if
[162]1906                 
[154]1907                  if ( sort .eq. "layer" ) then
[513]1908                     if (z_d(li) .eq. -1) then
[194]1909                        if (delta_z .EQ. -1) then
1910                           level = "-average"
1911                        else
[358]1912                           level = "=" + z_d(li) + "m"
[194]1913                        end if
[161]1914                     else
[358]1915                        level = "=" + z_d(li) + "m"
[190]1916                     end if
[357]1917       
[566]1918                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1919                        vNam(varn) .eq. "pras_xy"  .or. \
1920                        vNam(varn) .eq. "prrs_xy"  .or. \
1921                        vNam(varn) .eq. "qsws_xy"  .or. \
1922                        vNam(varn) .eq. "shfs_xy"  .or. \
1923                        vNam(varn) .eq. "ts_xy"    .or. \
1924                        vNam(varn) .eq. "us_xy"    .or. \
1925                        vNam(varn) .eq. "z0s_xy"   .or. \
1926                        vNam(varn) .eq. "lwp*_xy"  .or. \
1927                        vNam(varn) .eq. "pra*_xy"  .or. \
1928                        vNam(varn) .eq. "prr*_xy"  .or. \
1929                        vNam(varn) .eq. "qsws*_xy" .or. \
1930                        vNam(varn) .eq. "shf*_xy"  .or. \
1931                        vNam(varn) .eq. "t*_xy"    .or. \
1932                        vNam(varn) .eq. "u*_xy"    .or. \
1933                        vNam(varn) .eq. "z0*_xy") then
[357]1934                         lie = 0
1935                         level = "=" + zu1(0) + "m"
[769]1936                     else
1937                         lie = laye
[357]1938                     end if
[758]1939
[769]1940                     if(vNam(varn) .eq. "zwwi"  .or. \
1941                        vNam(varn) .eq. "zusi") then
1942                         lie = 0
1943                         lis = 0
1944                         loe = 0
1945                         los = 0
1946                         level = ""
1947                     else
1948                         lis = lays
1949                         lie = laye
1950                         los = start_time_step
1951                         loe = end_time_step
1952                     end if
1953
1954                     if(vNam(varn) .eq. "zwwi"  .or. \
1955                        vNam(varn) .eq. "zusi") then
1956                         cs_res@gsnCenterString = ""
1957                     else
1958                         cs_res@gsnCenterString = "t=" + \
[770]1959                          decimalPlaces(t_all(lo)/3600,2,True) +"h  z"+level
[769]1960                     end if
1961
[758]1962                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1963                        ;nothing to be done
1964                     else
1965                        plot(n) = gsn_csm_contour(wks_ps,\
[566]1966                                            data(varn,lo,li-lis,:,:),cs_res)
[758]1967                     end if
1968
[162]1969                     if (vector .EQ. 1 .AND. check_vecp) then
[566]1970                        vecres                 = True           ; vector only
1971                                                                ; resources
1972                        vecres@gsnDraw         = False          ; don't draw
1973                        vecres@gsnFrame        = False          ; don't
1974                                                                ; advance frame
1975                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
1976                        vecres@vcRefMagnitudeF = ref_mag        ; define
1977                                                                ; vector ref
1978                                                                ; mag
1979                        vecres@vcRefLengthF    = 0.05           ; define
1980                                                                ; length of
1981                                                                ; vec ref
1982                        vecres@gsnRightString  = " "            ; turn off
1983                                                                ; right string
1984                        vecres@gsnLeftString   = " "            ; turn off
1985                                                                ; left string
1986                        vecres@tiXAxisString   = " "   
1987                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1988                                                   vect2(lo,li-lis,:,:),vecres)
[162]1989                        overlay(plot(n), plot_vec)
[161]1990                     end if
[154]1991                  end if
1992               end if
1993
[157]1994               ; ****************************************************
[154]1995               ; xz cross section
[157]1996               ; ****************************************************
[154]1997
1998               if ( xzc .eq. 1 ) then
[218]1999                 
[157]2000                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
[162]2001               
[154]2002                  if ( sort .eq. "time" ) then
[358]2003                     if ( y_d(lo) .eq. -1 ) then
[154]2004                        level = "-average"
2005                     else
[358]2006                        level = "=" + y_d(lo) + "m"
[154]2007                     end if
[566]2008                     cs_res@gsnCenterString = "t=" + \
2009                           decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
[758]2010
2011                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2012                        ;nothing to be done
2013                     else
2014                         plot(n) = gsn_csm_contour(wks_ps,\
[566]2015                                              data(varn,li,:,lo-los,:),cs_res)
[758]2016                     end if
2017
[162]2018                     if (vector .EQ. 1 .AND. check_vecp) then
[566]2019                        vecres                 = True           ; vector only
2020                                                                ; resources
2021                        vecres@gsnDraw         = False          ; don't draw
2022                        vecres@gsnFrame        = False          ; don't
2023                                                                ; advance frame
2024                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2025                        vecres@vcRefMagnitudeF = ref_mag        ; define
2026                                                                ; vector ref
2027                                                                ; mag
2028                        vecres@vcRefLengthF    = 0.05           ; define
2029                                                                ; length of
2030                                                                ; vec ref
2031                        vecres@gsnRightString  = " "            ; turn off
2032                                                                ; right string
2033                        vecres@gsnLeftString   = " "            ; turn off
2034                                                                ; left string
2035                        vecres@tiXAxisString   = " "   
2036                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2037                                                   vect2(li,:,lo-los,:),vecres)
[162]2038                        overlay(plot(n), plot_vec)
[161]2039                     end if
[154]2040                  end if
2041
2042                  if ( sort .eq. "layer" ) then
[358]2043                     if ( y_d(li) .eq. -1 ) then
[161]2044                        level = "-average"
2045                     else
[358]2046                        level = "=" + y_d(li) + "m"
[161]2047                     end if
[566]2048                     cs_res@gsnCenterString = "t=" + \
2049                           decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
[758]2050
2051                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2052                        ;nothing to be done
2053                     else
2054                         plot(n) = gsn_csm_contour(wks_ps,\
2055                                              data(varn,lo,:,li-lis,:),cs_res) 
2056                     end if
2057
[162]2058                     if (vector .EQ. 1 .AND. check_vecp) then
[566]2059                        vecres                 = True           ; vector only
2060                                                                ; resources
2061                        vecres@gsnDraw         = False          ; don't draw
2062                        vecres@gsnFrame        = False          ; don't
2063                                                                ; advance frame
2064                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2065                        vecres@vcRefMagnitudeF = ref_mag        ; define
2066                                                                ; vector ref
2067                                                                ; mag
2068                        vecres@vcRefLengthF    = 0.05           ; define
2069                                                                ; length of
2070                                                                ; vec ref
2071                        vecres@gsnRightString  = " "            ; turn off
2072                                                                ; right string
2073                        vecres@gsnLeftString   = " "            ; turn off
2074                                                                ; left string
2075                        vecres@tiXAxisString   = " "   
2076                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2077                                                   vect2(lo,:,li-lis,:),vecres)
[162]2078                        overlay(plot(n), plot_vec)
[161]2079                     end if
[157]2080                  end if                 
[154]2081               end if
2082
[157]2083               ; ****************************************************
[154]2084               ; yz cross section
[157]2085               ; ****************************************************
[154]2086
2087               if ( yzc .eq. 1 ) then
[218]2088                                 
[157]2089                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2090
[154]2091                  if ( sort .eq. "time" ) then
[513]2092                     if ( x_d(lo) .eq. -1 ) then
[154]2093                        level = "-average"
2094                     else
[358]2095                        level = "=" + x_d(lo) + "m"
[154]2096                     end if
[566]2097                     cs_res@gsnCenterString = "t=" + \
2098                          decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
[758]2099
2100                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2101                        ;nothing to be done
2102                     else
2103                        plot(n) = gsn_csm_contour(wks_ps,\
[566]2104                                          data(varn,li,:,:,lo-los),cs_res)
[758]2105                     end if
2106
[162]2107                     if (vector .EQ. 1 .AND. check_vecp) then
[566]2108                        vecres                 = True           ; vector only
2109                                                                ; resources
2110                        vecres@gsnDraw         = False          ; don't draw
2111                        vecres@gsnFrame        = False          ; don't
2112                                                                ; advance frame
2113                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2114                        vecres@vcRefMagnitudeF = ref_mag        ; define
2115                                                                ; vector ref
2116                                                                ; mag
2117                        vecres@vcRefLengthF    = 0.05           ; define
2118                                                                ; length of
2119                                                                ; vec ref
2120                        vecres@gsnRightString  = " "            ; turn off
2121                                                                ; right string
2122                        vecres@gsnLeftString   = " "            ; turn off
2123                                                                ; left string
2124                        vecres@tiXAxisString   = " "   
2125                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2126                                                   vect2(li,:,:,lo-los),vecres)
[162]2127                        overlay(plot(n), plot_vec)
[161]2128                     end if
[154]2129                  end if
2130
2131                  if ( sort .eq. "layer" ) then
[513]2132                     if ( x_d(li) .eq. -1 ) then
[161]2133                        level = "-average"
2134                     else
[358]2135                        level = "=" + x_d(li) + "m"
[161]2136                     end if
[566]2137                     cs_res@gsnCenterString = "t=" + \
2138                           decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
[758]2139
2140                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2141                        ;nothing to be done
2142                     else
2143                        plot(n) = gsn_csm_contour(wks_ps,\
[566]2144                                          data(varn,lo,:,:,li-lis),cs_res)
[758]2145                     end if
2146
[162]2147                     if (vector .EQ. 1 .AND. check_vecp)then
[566]2148                        vecres                 = True           ; vector only
2149                                                                ;resources
2150                        vecres@gsnDraw         = False          ; don't draw
2151                        vecres@gsnFrame        = False          ; don't
2152                                                                ; advance frame
2153                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2154                        vecres@vcRefMagnitudeF = ref_mag        ; define
2155                                                                ; vector ref
2156                                                                ; mag
2157                        vecres@vcRefLengthF    = 0.05           ; define
2158                                                                ; length of
2159                                                                ; vec ref
2160                        vecres@gsnRightString  = " "            ; turn off
2161                                                                ; right string
2162                        vecres@gsnLeftString   = " "            ; turn off
2163                                                                ; left string
2164                        vecres@tiXAxisString   = " "   
2165                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2166                                                   vect2(lo,:,:,li-lis),vecres)
[162]2167                        overlay(plot(n), plot_vec)
[161]2168                     end if
[154]2169                  end if
[161]2170               end if         
[218]2171               n=n+1 
[154]2172            end do     
2173         end do 
[162]2174      end if
[161]2175   end do
[154]2176
[161]2177   ; ***************************************************
2178   ; merge plots onto one page
[357]2179   ; ***************************************************   
[532]2180
2181   no_frames = 0
2182
[529]2183   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
[566]2184      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2185           no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) then
2186         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),\
2187                                       (/no_var+1,no_layer*no_time/),cs_resP)
[162]2188         print(" ")
2189         print("Outputs to .eps or .epsi have only one frame")
2190         print(" ")
2191      else
[357]2192         do np = 0,dim_plot-1,no_rows*no_columns   
2193            if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]2194               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2195                                                                      cs_resP)
[532]2196               no_frames = no_frames + 1
[162]2197            else
[566]2198               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2199                                              (/no_rows,no_columns/),cs_resP)
[532]2200               no_frames = no_frames + 1
[162]2201            end if
2202         end do
2203      end if
[267]2204   else       
[566]2205      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2206                               .AND. dim_plot .gt. no_rows*no_columns) then
[357]2207         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
[162]2208         print(" ")
2209         print("Outputs to .eps or .epsi have only one frame")
2210         print(" ")
2211      else
[357]2212         do np = 0,dim_plot-1,no_rows*no_columns 
2213            if ( np + no_rows*no_columns .gt. dim_plot-1) then
[566]2214               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2215                                                                     cs_resP)
[532]2216               no_frames = no_frames + 1
[162]2217            else
[566]2218               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2219                                              (/no_rows,no_columns/),cs_resP)
[532]2220               no_frames = no_frames + 1
[162]2221            end if
2222         end do
2223      end if
[161]2224   end if
2225
[532]2226   if (format_out .EQ. "png" ) then
2227     png_output = new((/no_frames/), string)
2228     j = 0
2229
[957]2230     if (no_frames .eq. 1) then
2231        if (ncl_version .GE. 6 ) then
2232           png_output(0) = file_out+".png"
2233        else
2234           png_output(0) = file_out+".00000"+1+".png"
2235        end if
2236        ;using imagemagick's convert for reducing the white
2237        ;space around the plot
2238        cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2239              png_output(0) + " " + png_output(0)
[532]2240       system(cmd)
[957]2241     else
[532]2242
[957]2243       do i=0, no_frames-1
2244         j = i + 1
2245         if (j .LE. 9) then
2246           png_output(i) = file_out+".00000"+j+".png"
2247         end if
2248         if (j .GT. 9 .AND. j .LE. 99) then
2249           png_output(i) = file_out+".0000"+j+".png"
2250         end if
2251         if (j .GT. 99 .AND. j .LE. 999) then
2252           png_output(i) = file_out+".000"+j+".png"
2253         end if
2254         if (j .GT. 999) then
2255           png_output(i) = file_out+".00"+j+".png"
2256         end if
2257
2258         ;using imagemagick's convert for reducing the white
2259         ;space around the plot
2260         cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2261                png_output(i) + " " + png_output(i)
2262         system(cmd)
2263       end do
2264     end if
2265
[532]2266     print(" ")
2267     print("Output to: "+ png_output)
2268     print(" ")
2269   else
2270     print(" ")
2271     print("Output to: " + file_out +"."+ format_out)
2272     print(" ")
2273   end if
[957]2274
[190]2275end
Note: See TracBrowser for help on using the repository browser.