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

Last change on this file since 1018 was 983, checked in by hoffmann, 12 years ago

Bugfix in netcdf.f90 and improvements in palmplot.

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