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

Last change on this file since 1858 was 1626, checked in by heinze, 9 years ago

Bugfix to plot xy-cross sections in case variables contain only one level (like lwp*_xy, shf*_xy, ...)

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