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

Last change on this file since 3629 was 2837, checked in by gronemeier, 7 years ago

palmplot: config file of ncl scripts can be chosen in shell command

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