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

Last change on this file since 174 was 174, checked in by letzel, 16 years ago
  • NCL scripts in trunk/SCRIPTS/NCL updated
File size: 58.7 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
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
7begin
8
9   ; ***************************************************
10   ; read parameter_list
11   ; ***************************************************
12
13   if (isfilepresent(".ncl_preferences")) then
14      parameter = asciiread(".ncl_preferences",77,"string")
15      delete(parameter@_FillValue)
16   else
17      print(" ")
18      print("'.ncl_preferences' is not existent")
19      print(" ")
20      exit
21   end if
22
23   ; ***************************************************
24   ; set default parameter values and strings if not assigned in prompt or parameter list
25   ; ***************************************************
26   
27   if ( .not. isvar("file_in") ) then                   ; path+name of input file     
28      if (parameter(7) .EQ. "input file") then
29         print(" ")
30         print("Please provide input file 'file_in = ' either in prompt or parameter_list")
31         print(" ")
32         exit
33      else
34         file_in = parameter(7)
35      end if     
36   end if
37   if ( .not. isvar("format_out") ) then                ; format of output file
38      format_out = "x11"
39      if (parameter(9) .NE. "x11") then
40         format_out = parameter(9) 
41      end if
42   end if
43   if ( .not. isvar("file_out") ) then                  ; path+name of output file
44      file_out = "test"
45      if (parameter(11) .NE. "test") then
46         file_out = parameter(11) 
47     end if     
48   end if
49   if ( .not. isvar("no_columns") ) then                ; number of plots in one row
50      no_columns = 1
51      if (parameter(17) .NE. "1") then
52         no_columns = stringtointeger(parameter(17)) 
53      end if
54   end if
55   if ( .not. isvar("no_lines") ) then                  ; number of plot-lines on one sheet
56      no_lines = 2
57      if (parameter(19) .NE. "2") then
58         no_lines = stringtointeger(parameter(19)) 
59      end if
60   end if
61   if ( .not. isvar("sort") ) then                      ; sort of plots
62      sort = "time"
63      if (parameter(51) .NE. "time") then
64         sort = parameter(51) 
65      end if
66   end if
67   if ( .not. isvar("mode") ) then                      ; pattern of contour plots
68      mode = "Fill"
69      if (parameter(49) .NE. "Fill") then
70         mode = parameter(49) 
71      end if
72   end if
73   if ( .not. isvar("fill_mode") ) then                 ; pattern of filling
74      fill_mode = "AreaFill"
75      if (parameter(53) .NE. "AreaFill") then
76         mode = parameter(53) 
77      end if
78   end if
79   if ( .not. isvar("shape") ) then                     ; keeping of aspect ratio
80      shape = 1
81      if (parameter(55) .NE. "1") then
82         shape = stringtointeger(parameter(55))
83         if (stringtointeger(parameter(55)) .NE. 0) then
84            print(" ")
85            print("Please set 'shape' to 0 or 1")
86            print(" ")
87            exit
88         end if
89      end if 
90   end if 
91   if ( .not. isvar("var") ) then                       ; set output of all variables
92      check = True
93   end if
94
95   if ( .not. isvar("xyc") ) then                       ; turn xy cross-section on or off
96      xyc = 0
97      if (parameter(57) .NE. "0") then
98         xyc = stringtointeger(parameter(57))
99         if (stringtointeger(parameter(57)) .NE. 1) then
100            print(" ")
101            print("Please set 'xyc' to 0 or 1")
102            print(" ")
103            exit
104         end if
105      end if
106   end if
107   if ( .not. isvar("xzc") ) then                       ; turn xz cross-section on or off
108      xzc = 0
109      if (parameter(59) .NE. "0") then
110         xzc = stringtointeger(parameter(59))
111         if (stringtointeger(parameter(59)) .NE. 1) then
112            print(" ")
113            print("Please set 'xzc' to 0 or 1")
114            print(" ")
115            exit
116         end if
117      end if
118   end if
119   if ( .not. isvar("yzc") ) then                       ; turn yz cross-section on or off
120      yzc = 0
121      if (parameter(61) .NE. "0") then
122         yzc = stringtointeger(parameter(61))
123         if (stringtointeger(parameter(61)) .NE. 1) then
124            print(" ")
125            print("Please set 'yzc' to 0 or 1")
126            print(" ")
127            exit
128         end if
129      end if
130   end if
131
132   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
133      print(" ")
134      print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)")
135      print(" ")
136      exit
137   end if
138   if (xyc .EQ. 1 ) then
139      if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
140         print(" ")
141         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
142         print(" ")
143         exit
144      end if
145   end if
146   if (xzc .EQ. 1) then
147      if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
148         print(" ")
149         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
150         print(" ")
151         exit
152      end if
153   end if
154   if (yzc .EQ. 1) then
155      if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
156         print(" ")
157         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
158         print(" ")
159         exit
160      end if
161   end if
162   if ( .not. isvar("vector") ) then                    ; sort of plots
163      vector = 0
164      if (parameter(39) .NE. "0") then
165         vector = stringtointeger(parameter(39))
166         if (stringtointeger(parameter(39)) .NE. 1) then
167            print(" ")
168            print("Please set 'vector' to 0 or 1")
169            print(" ")
170            exit
171         end if   
172      end if
173   end if
174   if ( .not. isvar("ref_mag") ) then                           ; sort of plots
175      ref_mag = 0.05
176      if (parameter(47) .NE. "0.05") then
177         ref_mag = stringtofloat(parameter(47))   
178      end if
179   end if
180
181   ; ***************************************************
182   ; open input file
183   ; ***************************************************
184
185   f  = addfile( file_in, "r" )
186
187   vNam  = getfilevarnames(f)
188   print(" ")
189   print("Variable on netCDF file: " + vNam)
190   print(" ")
191   dim   = dimsizes(vNam)
192
193   ; ***************************************************
194   ; check for kind of input file
195   ; ***************************************************
196
197   check_3d = True
198   do varn=0,dim-1
199      checkxy = isStrSubset( vNam(varn),"xy")
200      checkxz = isStrSubset( vNam(varn),"xz")
201      checkyz = isStrSubset( vNam(varn),"yz")
202      if (checkxy .OR. checkxz .OR. checkyz) then
203         check_3d = False
204         break
205      end if
206   end do
207   if (.not. check_3d) then
208      if (xyc .EQ. 1 .AND. .not. checkxy) then
209         print(" ")
210         print("Your input file doesn't have values for xy cross-sections")
211         if (checkxz)then
212            print("Select another input file or xz cross-sections")
213         else
214            print("Select another input file or yz cross-sections")
215         end if
216         print(" ")
217         exit
218      else
219         print(" ")
220         print("Your input file contains xy data")
221         print(" ")
222      end if
223      if (xzc .EQ. 1 .AND. .not. checkxz) then
224         print(" ")
225         print("Your input file doesn't have values for xz cross-sections")
226         if (checkxy)then
227            print("Select another input file or xy cross-sections")
228         else
229            print("Select another input file or yz cross-sections")
230         end if
231         print(" ")
232         exit
233      else
234         print(" ")
235         print("Your input file contains xz data")
236         print(" ")
237      end if
238      if (yzc .EQ. 1 .AND. .not. checkyz) then
239         print(" ")
240         print("Your input file doesn't have values for yz cross-sections")
241         if (checkxy)then
242            print("Select another input file or xy cross-sections")
243         else
244            print("Select another input file or xz cross-sections")
245         end if
246         print(" ")
247         exit
248      else
249         print(" ")
250         print("Your input file contains yz data")
251         print(" ")
252      end if
253   else
254      print(" ")
255      print("Your input file: '"+file_in+"'")
256      print("contains 3d or other data")
257      print(" ")
258   end if
259   
260   if (dim .EQ. 0) then
261      print(" ")
262      print("There is no data on file")
263      print(" ")
264      exit
265   end if
266
267   ; ***************************************************
268   ; set up recourses
269   ; ***************************************************
270
271   cs_res                          = True
272   cs_res@gsnYAxisIrregular2Linear = True 
273
274   if( shape .eq. 1 ) then                              ; keep the shape
275      cs_res@gsnShape = True
276   end if
277 
278   cs_res@gsnDraw                 = False
279   cs_res@gsnFrame                = False
280   cs_res@gsnMaximize             = True
281   cs_res@gsnPaperOrientation     = "portrait"
282   cs_res@gsnPaperWidth           = 8.27
283   cs_res@gsnPaperHeight          = 11.69
284   cs_res@gsnPaperMargin          = 0.79
285   cs_res@tmXBLabelFontHeightF   = .02
286   cs_res@tmYLLabelFontHeightF   = .02
287   cs_res@tiXAxisFontHeightF     = .02
288   cs_res@tiYAxisFontHeightF     = .02
289   cs_res@tmXBMode                ="Automatic"
290   cs_res@tmYLMode                ="Automatic"
291   cs_res@lgTitleFontHeightF     = .02
292   cs_res@lgLabelFontHeightF     = .02
293   cs_res@txFontHeightF          = .02
294   cs_res@cnLevelSelectionMode    = "ManualLevels"
295
296   cs_resP = True
297   cs_resP@txString               = f@title
298   cs_resP@txFuncCode             = "~"                 ; necessary for title
299   cs_resP@txFontHeightF          = .02
300   
301 
302   if ( mode .eq. "Fill" ) then
303      cs_res@cnFillOn             = True
304      cs_res@gsnSpreadColors      = True
305      cs_res@cnFillMode           = fill_mode
306      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
307      cs_res@cnLinesOn            = False       
308      cs_res@cnLineLabelsOn       = False
309   end if
310
311   if ( mode .eq. "Both" ) then
312      cs_res@cnFillOn             = True
313      cs_res@gsnSpreadColors      = True
314      cs_res@cnFillMode           = fill_mode
315      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
316      cs_res@cnLinesOn            = True
317      cs_res@cnLineLabelsOn       = True
318   end if
319
320   ; *********************************************
321   ; set up of start_time_step and end_time_step
322   ; *********************************************
323
324   t_all = f->time
325   nt = dimsizes(t_all)
326   delta_t = t_all(nt-1)/nt
327
328   ; ****************************************************       
329   ; start of time step and different types of mistakes that could be done
330   ; ****************************************************
331
332   if ( .not. isvar("start_time_step") ) then           
333      start_time_step=t_all(0)/3600
334      if (parameter(13) .NE. "t(0)") then
335         if (stringtodouble(parameter(13)) .GT. t_all(nt-1)/3600)
336            print(" ")
337            print("'start_time_step' = "+ parameter(13) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
338            print(" ")
339            print("Please select another 'start_time_step'")
340            print(" ")
341            exit
342         end if
343         if (stringtofloat(parameter(13)) .LT. t_all(0)/3600)
344            print(" ")
345            print("'start_time_step' = "+ parameter(13) +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
346            print(" ")
347            print("Please select another 'start_time_step'")
348            print(" ")
349            exit
350         end if
351         start_time_step=stringtodouble(parameter(13))
352      end if
353   else
354      if (start_time_step .GT. t_all(nt-1)/3600)
355         print(" ")
356         print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
357         print(" ")
358         print("Please select another 'start_time_step'")
359         print(" ")
360         exit
361      end if
362      if (start_time_step .LT. t_all(0)/3600)
363         print(" ")
364         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
365         print(" ")
366         print("Please select another 'start_time_step'")
367         print(" ")
368         exit
369      end if
370   end if
371   start_time_step = start_time_step*3600
372   do i=0,nt-1     
373      if (start_time_step .GE. t_all(i)-delta_t/2 .AND. start_time_step .LT. t_all(i)+delta_t/2)then
374         st=i
375         break
376      end if
377   end do
378       
379   ; ****************************************************
380   ; end of time step and different types of mistakes that could be done
381   ; ****************************************************
382
383   if ( .not. isvar("end_time_step") ) then             
384      end_time_step = t_all(nt-1)/3600
385      if (parameter(15) .NE. "t(end)") then
386         if (stringtodouble(parameter(15)) .LT. t_all(0)/3600)
387            print(" ")
388            print("'end_time_step' = "+parameter(15)+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
389            print(" ")
390            print("Please select another 'end_time_step'")
391            print(" ")
392            exit
393         end if
394         if (stringtodouble(parameter(15)) .GT. t_all(nt-1)/3600)
395            print(" ")
396            print("'end_time_step' = "+ parameter(15) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
397            print(" ")
398            print("Please select another 'end_time_step'") 
399            print(" ")
400            exit
401         end if
402         if (stringtodouble(parameter(15)) .LT. start_time_step/3600)
403            print(" ")
404            print("'end_time_step' = "+ parameter(15) +"h is lower than 'start_time_step' = "+parameter(13)+"h")
405            print(" ")
406            print("Please select another 'start_time_step' or 'end_time_step'")
407            print(" ")
408            exit
409         end if
410         end_time_step = stringtodouble(parameter(15))
411      end if   
412   else
413      if (end_time_step .LT. t_all(0)/3600)
414         print(" ")
415         print("'end_time_step' = "+end_time_step+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
416         print(" ")
417         print("Please select another 'end_time_step'")
418         print(" ")
419         exit
420      end if
421      if (end_time_step .GT. t_all(nt-1)/3600)
422         print(" ")
423         print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
424         print(" ")
425         print("Please select another 'end_time_step'") 
426         print(" ")
427         exit
428      end if
429      if (end_time_step .LT. start_time_step/3600)
430         print(" ")
431         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
432         print(" ")
433         print("Please select another 'start_time_step' or 'end_time_step'")
434         print(" ")
435         exit
436      end if
437   end if
438   end_time_step = end_time_step*3600
439   do i=0,nt-1     
440      if (end_time_step .GE. t_all(i)-delta_t/2 .AND. end_time_step .LT. t_all(i)+delta_t/2)then
441         et=i
442         break
443      end if
444   end do
445 
446   delete(start_time_step)
447   start_time_step=round(st,3)
448   delete(end_time_step)
449   end_time_step=round(et,3)
450
451   print(" ")
452   print("Output from "+t_all(start_time_step)/3600+"h = "+t_all(start_time_step)+"s => index = "+start_time_step)
453   print("       till "+t_all(end_time_step)/3600+"h = "+t_all(end_time_step)+"s => index = "+end_time_step)
454   print(" ")
455 
456   no_time=(end_time_step-start_time_step)+1
457
458   ; ****************************************************
459   ; set up variables for vector plot if required
460   ; ****************************************************   
461
462   if (vector .EQ. 1) then
463      if ( .not. isvar("plotvec") ) then
464         plotvec = parameter(45)
465      end if
466      if ( .not. isvar("vec1") ) then
467         vec1 = parameter(41)
468         if (parameter(41) .EQ. "vec1") then
469            print(" ")
470            print("Please indicate Vector 1 ('vec1') for Vector-Plot")
471            print(" ")
472            exit
473         end if
474      end if
475      if ( .not. isvar("vec2") ) then
476         vec2 = parameter(43)
477         if (parameter(43) .EQ. "vec2") then
478            print(" ")
479            print("Please indicate Vector 2 ('vec2') for Vector-Plot")
480            print(" ")
481            exit
482         end if
483      end if             
484   end if
485
486   check_vec1 = False
487   check_vec2 = False
488   check_vecp = False
489
490   ; ****************************************************
491   ; get data for plots
492   ; ****************************************************
493
494   if (xyc .EQ. 1) then
495      do varn=0,dim-1
496         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x")then   
497            x_d     = f->$vNam(varn)$
498            xdim    = dimsizes(x_d)-1
499            delta_x = x_d(1)-x_d(0)
500            break
501         end if
502      end do
503      do varn=0,dim-1
504         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
505            y_d     = f->$vNam(varn)$
506            ydim    = dimsizes(y_d)-1
507            delta_y = y_d(1)-y_d(0)
508            break
509         end if
510      end do
511      do varn=0,dim-1
512         if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
513            z_d     = f->$vNam(varn)$
514            zdim    = dimsizes(z_d)-1
515            delta_z = 0
516            break
517         else
518            if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
519               z_d     = f->$vNam(varn)$
520               zdim    = dimsizes(z_d)-1
521               delta_z = -1.d
522               break
523            else
524               zdim    = 0
525               delta_z = -1.d
526            end if
527         end if
528      end do
529   end if
530   if (xzc .EQ. 1) then
531      do varn=0,dim-1
532         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then
533            x_d     = f->$vNam(varn)$
534            xdim    = dimsizes(x_d)-1
535            delta_x = x_d(1)-x_d(0)
536            break
537         end if
538      end do
539      do varn=0,dim-1   
540         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
541            z_d     = f->$vNam(varn)$
542            zdim    = dimsizes(z_d)-1
543            delta_z = z_d(1)-z_d(0)
544            break
545         end if
546      end do
547      do varn=0,dim-1
548         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
549            y_d     = f->$vNam(varn)$
550            ydim    = dimsizes(y_d)-1
551            delta_y = y_d(1)-y_d(0)
552            break
553         else
554            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
555               y_d     = f->$vNam(varn)$
556               ydim    = dimsizes(y_d)-1
557               delta_y = -1.d
558               break
559            else
560               ydim    = 0
561               delta_y = -1.d
562            end if
563         end if
564      end do
565   end if
566   if (yzc .EQ. 1) then
567      do varn=0,dim-1 
568         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
569            y_d     = f->$vNam(varn)$
570            ydim    = dimsizes(y_d)-1
571            delta_y = y_d(1)-y_d(0)
572            break
573         end if
574      end do
575      do varn=0,dim-1
576         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
577            z_d     = f->$vNam(varn)$
578            zdim    = dimsizes(z_d)-1
579            delta_z = z_d(1)-z_d(0)
580            break
581         end if
582      end do
583      do varn=0,dim-1
584         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
585            x_d     = f->$vNam(varn)$
586            xdim    = dimsizes(x_d)-1
587            delta_x = x_d(1)-x_d(0)
588            break
589         else
590            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
591               x_d     = f->$vNam(varn)$
592               xdim    = dimsizes(x_d)-1
593               delta_x = -1.d
594               break
595            else
596               xdim    = 0
597               delta_x = -1.d
598            end if
599         end if
600      end do
601   end if
602   
603   ; ****************************************************
604   ; set up ranges of x-, y- and z-coordinates
605   ; ****************************************************         
606                   
607   if ( .not. isvar("xs") ) then               
608      xs     = 0.0d
609      xstart = 0
610      if (parameter(63) .NE. "x0") then
611         if (delta_x .EQ. -1) then
612            print(" ")
613            print("You cannot choose a start value for x, there are preseted layers for x")
614            print(" ")
615            xstart=0
616         else
617            if (stringtodouble(parameter(63)) .LT. 0-delta_x/2) then
618               print(" ")
619               print("range start for x-coordinate = "+parameter(63)+"m is lower than first value x = "+0+"m or xu = "+(0-delta_x/2)+"m")
620               print(" ")
621               exit
622            end if
623            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
624               if (stringtodouble(parameter(63)) .GE. xdim*delta_x) then
625                  print(" ")
626                  print("range start for x-coordinate = "+parameter(63)+"m is equal or greater than last value = "+xdim*delta_x+"m")
627                  print(" ")
628                  exit
629               end if
630            else
631               if (stringtodouble(parameter(63)) .GT. xdim*delta_x) then
632                  print(" ")
633                  print("range start for x-coordinate = "+parameter(63)+"m is greater than last value = "+xdim*delta_x+"m")
634                  print(" ")
635                  exit
636               end if
637            end if
638            xs = stringtodouble(parameter(63))
639         end if
640      end if
641   else
642      if (delta_x .EQ. -1) then
643         print(" ")
644         print("You cannot choose a start value for x, there are preseted layers for x")
645         print(" ")
646         xstart=0
647      else
648         if (xs .LT. 0-delta_x/2) then
649            print(" ")
650            print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m")
651            print(" ")
652            exit
653         end if
654         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
655            if (xs .GE. xdim*delta_x) then
656               print(" ")
657               print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m")
658               print(" ")
659               exit
660            end if
661         else
662            if (xs .GT. xdim*delta_x) then
663               print(" ")
664               print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m")
665               print(" ")
666               exit
667            end if
668         end if
669      end if
670   end if
671
672   if ( .not. isvar("ys") ) then   
673      ys = 0.0d
674      ystart = 0
675      if (parameter(67) .NE. "y0") then
676         if (delta_y .EQ. -1) then
677            print(" ")
678            print("You cannot choose a start value for y, there are preseted layers for y")
679            print(" ")
680            ystart=0
681         else
682            if (stringtodouble(parameter(67)) .LT. 0-delta_y/2) then
683               print(" ")
684               print("range start for y-coordinate = "+parameter(67)+"m is lower than first value = "+0-delta_y/2+"m")
685               print(" ")
686               exit
687            end if
688            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
689               if (stringtodouble(parameter(67)) .GE. ydim*delta_y) then
690                  print(" ")
691                  print("range start for y-coordinate = "+parameter(67)+"m is equal or greater than last value = "+ydim*delta_y+"m")
692                  print(" ")
693                  exit
694               end if
695            else
696               if (stringtodouble(parameter(67)) .GT. ydim*delta_y) then
697                  print(" ")
698                  print("range start for y-coordinate = "+parameter(67)+"m is greater than last value = "+ydim*delta_y+"m")
699                  print(" ")
700                  exit
701               end if
702            end if
703            ys = stringtodouble(parameter(67))
704         end if
705      end if
706   else
707      if (delta_y .EQ. -1) then
708         print(" ")
709         print("You cannot choose a start value for y, there are preseted layers for y")
710         print(" ")
711         ystart=0
712      else
713         if (ys .LT. 0-delta_y/2) then
714            print(" ")
715            print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m")
716            print(" ")
717            exit
718         end if
719         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
720            if (ys .GE. ydim*delta_y) then
721               print(" ")
722               print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m")
723               print(" ")
724               exit
725            end if
726         else
727            if (ys .GT. ydim*delta_y) then
728               print(" ")
729               print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
730               print(" ")
731               exit
732            end if
733         end if
734      end if
735   end if
736 
737   if ( .not. isvar("zs") ) then                        ; start of z-coordinate range
738      zs = 0
739      if (parameter(71) .NE. "z0") then
740         if (delta_z .EQ. -1) then
741            print(" ")
742            print("You cannot choose a start value for z, there are preseted layers for z")
743            print(" ")
744         else
745            print(" ")
746            print("Please mind to indicate start and end ranges for the z-coordinate in")
747            print("indices not in 'meters'. Corresponding index and meter:")
748            print(" ")
749            print(" = "+z_d+" m")
750            print(" ")
751            if (stringtointeger(parameter(71)) .LT. 0) then
752               print(" ")
753               print("range start for z-coordinate = "+parameter(71)+" is lower than first gridpoint = 0")
754               print(" ")
755               exit
756            end if
757            if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
758               if (stringtointeger(parameter(71)) .GE. zdim) then
759                  print(" ")
760                  print("range start for z-coordinate = "+parameter(71)+" is equal or greater than last gridpoint = "+zdim)
761                  print(" ")
762                  exit
763               end if
764            else
765               if (stringtodouble(parameter(71)) .GT. zdim) then
766                  print(" ")
767                  print("range start for z-coordinate = "+parameter(71)+" is greater than last gridpoint = "+zdim)
768                  print(" ")
769                  exit
770               end if
771            end if
772            zs = stringtointeger(parameter(71))
773         end if
774      end if
775   else
776      if (delta_z .EQ. -1) then
777         print(" ")
778         print("You cannot choose a start value for z, there are preseted layers for z")
779         print(" ")
780      else
781         if (zs .LT. 0) then
782            print(" ")
783            print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0")
784            print(" ")
785            exit
786         end if
787         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
788            if (zs .GE. zdim) then
789               print(" ")
790               print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim)
791               print(" ")
792               exit
793            end if
794         else
795            if (zs .GT. zdim) then
796               print(" ")
797               print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim)
798               print(" ")
799               exit
800            end if
801         end if
802      end if
803   end if 
804
805   if ( .not. isvar("xe")) then   
806      xe   = xdim*delta_x
807      xend = xdim
808      if (parameter(65) .NE. "xdim") then
809         if (delta_x .EQ. -1) then
810            print(" ")
811            print("You cannot choose an end value for x, there are preseted layers for x")
812            print(" ")
813            xend=xdim           
814         else
815            if (stringtodouble(parameter(65)) .GT. xdim*delta_x) then
816               print(" ")
817               print("range end for x-coordinate = "+parameter(65)+"m is greater than last value = "+xdim*delta_x+"m")
818               print(" ")
819               exit
820            end if
821            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
822               if (stringtodouble(parameter(65)) .LE. 0-delta_x/2)
823                  print(" ")
824                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
825                  print(" ")
826                  exit
827               end if
828               if (stringtodouble(parameter(65)) .LE. xs) then
829                  print(" ")
830                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than start range = "+xs+"m")
831                  print(" ")
832                  exit
833               end if
834            else
835               if (stringtodouble(parameter(65)) .LT. 0-delta_x/2)
836                  print(" ")
837                  print("range end for x-coordinate = "+parameter(65)+"m is lower than first value = "+(0-delta_x/2)+"m")
838                  print(" ")
839                  exit
840               end if
841               if (stringtodouble(parameter(65)) .LT. xs) then
842                  print(" ")
843                  print("range end for x-coordinate = "+parameter(65)+"m is lower than start range = "+xs+"m")
844                  print(" ")
845                  exit
846               end if
847            end if         
848            xe = stringtodouble(parameter(65))
849         end if
850      end if
851   else
852      if (delta_x .EQ. -1) then
853         print(" ")
854         print("You cannot choose an end value for x, there are preseted layers for x")
855         print(" ")
856         xend=xdim
857      else
858         if (xe .GT. xdim*delta_x) then
859            print(" ")
860            print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m")
861            print(" ")
862            exit
863         end if
864         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
865            if (xe .LE. 0-delta_x/2)
866               print(" ")
867               print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
868               print(" ")
869               exit
870            end if
871            if (xe .LE. xs) then
872               print(" ")
873               print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m")
874               print(" ")
875               exit
876            end if
877            if (stringtodouble(xe .EQ. xs+1)) then
878               print(" ")
879               print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m")
880               print(" ")
881               exit
882            end if
883         else
884            if (xe .LT. 0-delta_x/2)
885               print(" ")
886               print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
887               print(" ")
888               exit
889            end if
890            if (xe .LT. xs) then
891               print(" ")
892               print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
893               print(" ")
894               exit
895            end if
896         end if
897      end if               
898   end if
899
900   if ( .not. isvar("ye")) then 
901      ye   = ydim*delta_y
902      yend = ydim
903      if (parameter(69) .NE. "ydim") then
904         if (delta_y .EQ. -1) then
905            print(" ")
906            print("You cannot choose an end value for y, there are preseted layers for y")
907            print(" ")
908            yend=ydim
909         else
910            if (stringtodouble(parameter(69)) .GT. ydim*delta_y) then
911               print(" ")
912               print("range end for y-coordinate = "+parameter(69)+"m is greater than last value = "+ydim*delta_y+"m")
913               print(" ")
914               exit
915            end if
916            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
917               if (stringtodouble(parameter(69)) .LE. 0-delta_y/2)
918                  print(" ")
919                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
920                  print(" ")
921                  exit
922               end if
923               if (stringtodouble(parameter(69)) .LE. ys) then
924                  print(" ")
925                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than start range = "+ys+"m")
926                  print(" ")
927                  exit
928               end if
929               if (stringtodouble(parameter(69)) .EQ. ys+1) then
930                  print(" ")
931                  print("range end for y-coordinate = "+parameter(69)+"m must be at least two more gridpoints greater than start range = "+ys+"m")
932                  print(" ")
933                  exit
934               end if
935            else
936               if (stringtodouble(parameter(69)) .LT. 0-delta_y/2)
937                  print(" ")
938                  print("range end for y-coordinate = "+parameter(69)+"m is lower than first value = "+(0-delta_y/2)+"m")
939                  print(" ")
940                  exit
941               end if
942               if (stringtodouble(parameter(69)) .LT. ys) then
943                  print(" ")
944                  print("range end for y-coordinate = "+parameter(69)+"m is lower than start range = "+ys+"m")
945                  print(" ")
946                  exit
947               end if
948            end if         
949            ye = stringtodouble(parameter(69))
950         end if
951      end if
952   else
953      if (delta_y .EQ. -1) then
954         print(" ")
955         print("You cannot choose an end value for y, there are preseted layers for y")
956         print(" ")
957         yend=ydim
958      else
959         if (ye .GT. ydim*delta_y) then
960            print(" ")
961            print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m")
962            print(" ")
963            exit
964         end if
965         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
966            if (ye .LE. 0-delta_y/2)
967               print(" ")
968               print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
969               print(" ")
970               exit
971            end if
972            if (ye .LE. ys) then
973               print(" ")
974               print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m")
975               print(" ")
976               exit
977            end if
978            if (ye .EQ. ys+1) then
979               print(" ")
980               print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m")
981               print(" ")
982               exit
983            end if
984         else
985            if (ye .LT. 0-delta_y/2)
986               print(" ")
987               print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m")
988               print(" ")
989               exit
990            end if
991            if (ye .LT. ys) then
992               print(" ")
993               print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m")
994               print(" ")
995               exit
996            end if
997         end if
998      end if
999   end if
1000 
1001   if ( .not. isvar("ze")) then 
1002      ze = zdim
1003      if (parameter(73) .NE. "zdim") then
1004         if (delta_z .EQ. -1) then
1005            print(" ")
1006            print("You cannot choose an end value for z, there are preseted layers for z")
1007            print(" ")           
1008         else
1009            if (stringtointeger(parameter(73)) .GT. zdim) then
1010               print(" ")
1011               print("range end for z-coordinate = "+parameter(73)+" is greater than last gridpoint = "+zdim)
1012               print(" ")
1013               exit
1014            end if
1015         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
1016            if (stringtointeger(parameter(73)) .LE. 0)
1017               print(" ")
1018               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than first gridpoint = 0")
1019               print(" ")
1020               exit
1021            end if
1022            if (stringtointeger(parameter(73)) .LE. zs) then
1023               print(" ")
1024               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than start range = "+zs)
1025               print(" ")
1026               exit
1027            end if
1028            if (stringtodouble(parameter(73)) .EQ. zs+1) then
1029               print(" ")
1030               print("range end for z-coordinate = "+parameter(73)+" must be at least two more gridpoints greater than start range = "+zs)
1031               print(" ")
1032               exit
1033            end if
1034         else
1035            if (stringtointeger(parameter(73)) .LT. 0)
1036               print(" ")
1037               print("range end for z-coordinate = "+parameter(73)+" is lower than first gridpoint = 0")
1038               print(" ")
1039               exit
1040            end if
1041            if (stringtointeger(parameter(73)) .LT. zs) then
1042               print(" ")
1043               print("range end for z-coordinate = "+parameter(73)+" is lower than start range = "+zs)
1044               print(" ")
1045               exit
1046            end if
1047         end if         
1048         ze = stringtointeger(parameter(73))
1049      end if
1050      end if
1051   else
1052      if (delta_z .EQ. -1) then
1053         print(" ")
1054         print("You cannot choose an end value for z, there are preseted layers for z")
1055         print(" ")
1056         ze = zdim
1057      else
1058         if (ze .GT. zdim) then
1059            print(" ")
1060            print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim)
1061            print(" ")
1062            exit
1063         end if
1064         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
1065            if (ze .LE. 0)
1066               print(" ")
1067               print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0")
1068               print(" ")
1069               exit
1070            end if
1071            if (ze .LE. zs) then
1072               print(" ")
1073               print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs)
1074               print(" ")
1075               exit
1076            end if
1077            if (ze .EQ. zs+1) then
1078               print(" ")
1079               print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs)
1080               print(" ")
1081               exit
1082            end if
1083         else
1084            if (ze .LT. 0)
1085               print(" ")
1086               print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0")
1087               print(" ")
1088               exit
1089            end if
1090            if (ze .LT. zs) then
1091               print(" ")
1092               print("range end for z-coordinate = "+ze+" is lower than start range = "+zs)
1093               print(" ")
1094               exit
1095            end if
1096         end if
1097      end if
1098   end if
1099
1100   if (delta_x .NE. -1) then
1101      do i=0,xdim   
1102         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
1103            xstart=i
1104            break
1105         end if
1106      end do
1107      do i=0,xdim   
1108         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
1109            xend=i
1110            break
1111         end if
1112      end do
1113   end if
1114   if (delta_y .NE. -1) then
1115      do i=0,xdim   
1116         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
1117            ystart=i
1118            break
1119         end if
1120      end do
1121      do i=0,xdim   
1122         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
1123            yend=i
1124            break
1125         end if
1126      end do
1127   end if
1128   
1129   delete(xs)
1130   delete(xe)   
1131   delete(ys)
1132   delete(ye)
1133
1134   xs=xstart
1135   xe=xend
1136   ys=ystart
1137   ye=yend
1138
1139   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1140      if (xe .EQ. xs+1) then
1141         print(" ")
1142         print("range end for x-coordinate="+parameter(65)+"m(=>value on file="+xe*delta_x+"m) must be at least two")
1143         print("more gridpoints(="+2*delta_x+"m) greater than start range="+parameter(63)+"m(=>value on file="+xs*delta_x+"m)")
1144         print(" ")
1145         exit
1146      end if
1147   end if
1148   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1149      if (ye .EQ. ys+1) then
1150         print(" ")
1151         print("range end for y-coordinate="+parameter(69)+"m(=>value on file="+ye*delta_y+"m) must be at least two")
1152         print("more gridpoints(="+2*delta_y+"m) greater than start range="+parameter(67)+"m(=>value on file="+ys*delta_y+"m)")
1153         print(" ")
1154         exit
1155      end if
1156   end if
1157   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1158      if (ze .EQ. zs+1) then
1159         print(" ")
1160         print("range end for x-coordinate="+parameter(73)+"(=> value on file="+ze+") must be at least two")
1161         print("more gridpoints greater than start range="+parameter(71)+"(=>value on file="+zs+")")
1162         print(" ")
1163         exit
1164      end if
1165   end if
1166 
1167   if (xyc .EQ. 1) then
1168      no_layer = (ze-zs)+1
1169      data = new((/dim,nt,zdim+1,(ye-ys)+1,(xe-xs)+1/),float)
1170   end if
1171   if (xzc .EQ. 1) then
1172      no_layer = (ye-ys)+1
1173      data = new((/dim,nt,(ze-zs)+1,ydim+1,(xe-xs)+1/),float)
1174   end if
1175   if (yzc .EQ. 1) then
1176      no_layer = (xe-xs)+1
1177      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,xdim+1/),float)
1178   end if
1179
1180   MinVal = new(dim,float)
1181   MaxVal = new(dim,float)
1182
1183   ; ****************************************************
1184   ; define inner and outer loops depending on "sort"
1185   ; ****************************************************       
1186   
1187   if ( xyc .eq. 1 ) then
1188      lays = zs
1189      laye = ze
1190   end if
1191   if ( xzc .eq. 1 ) then
1192      lays = ys
1193      laye = ye
1194   end if
1195   if ( yzc .eq. 1) then
1196      lays = xs
1197      laye = xe
1198   end if
1199
1200   if ( sort .eq. "time" ) then
1201      lis = start_time_step
1202      lie = end_time_step
1203      los = lays
1204      loe = laye
1205   else
1206      lis = lays
1207      lie = laye
1208      los = start_time_step
1209      loe = end_time_step
1210   end if
1211
1212   check = True
1213   v1=0
1214   v2=0
1215   no_var=0
1216   n=0
1217
1218   do varn=dim-1,0,1 
1219     
1220      data_all = f->$vNam(varn)$
1221   
1222      if (vector .EQ. 1) then   
1223         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1224         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1225      end if
1226           
1227      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
1228         check = False
1229      else
1230         check = True
1231      end if 
1232      if (  isvar("var") ) then 
1233         check = isStrSubset( var,","+vNam(varn)+",")
1234      end if 
1235      if (parameter(21) .NE. "variables") then
1236         var = parameter(21)
1237         check = isStrSubset( var,","+vNam(varn)+"," )
1238      end if 
1239
1240      if(check) then
1241
1242         no_var=no_var+1
1243
1244         if (vector .EQ. 1) then
1245            if (","+vNam(varn)+"," .EQ. vec1) then
1246               v1=v1+1
1247            end if
1248            if (","+vNam(varn)+"," .EQ. vec2) then
1249               v2=v2+1
1250            end if
1251         end if
1252
1253         if (xyc .EQ. 1) then
1254            data(varn,:,:,:,:)=data_all(:,:,ys:ye,xs:xe)
1255         end if
1256         if ( xzc .eq. 1 ) then
1257            data(varn,:,:,:,:)=data_all(:,zs:ze,:,xs:xe)
1258         end if
1259         if ( yzc .eq. 1) then
1260            data(varn,:,:,:,:)=data_all(:,zs:ze,ys:ye,:)
1261         end if
1262
1263         if (check_vec1) then
1264            vect1=data_all
1265         end if
1266         if (check_vec2) then
1267            vect2=data_all
1268         end if
1269
1270         data!0 = "var"
1271         data!1 = "t"
1272         data!2 = "z"
1273         data!3 = "y"
1274         data!4 = "x" 
1275
1276         MinVal(varn) = min(data(varn,:,:,:,:))
1277         MaxVal(varn) = max(data(varn,:,:,:,:))
1278
1279      end if
1280
1281      delete(data_all)
1282
1283   end do
1284
1285   if (no_var .EQ. 0) then
1286      print(" ")
1287      print("Please select a variable 'var=' or use the default value")
1288      print(" ")
1289      print("Your selection '"+var+"' does not exist on the input file")
1290      print(" ")
1291      exit
1292   end if
1293
1294   if (vector .EQ. 1) then
1295      if (v1 .EQ. 0)then
1296         print(" ")
1297         print("Varible for vector-plot ('vec1') must be one of the varibles for plot ('var')")
1298         print(" ")
1299         exit
1300      end if
1301
1302      if (v2 .EQ. 0)then
1303         print(" ")
1304         print("Varible for vector-plot ('vec2') must be one of the varibles for plot ('var')")
1305         print(" ")
1306         exit
1307      end if
1308   end if
1309
1310   ; ***************************************************
1311   ; open workstation(s)
1312   ; ***************************************************
1313
1314   wks_ps  = gsn_open_wks(format_out,file_out)
1315   gsn_define_colormap(wks_ps,"rainbow+white")
1316
1317   plot=new((/no_time*no_layer*no_var*2/),graphic)
1318 
1319   page = 0
1320   n=0
1321   print(" ")
1322   print("Plots sorted by '"+sort+"'")
1323   print(" ")
1324
1325   ; ***************************************************
1326   ; create plots
1327   ; ***************************************************
1328
1329   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec") then
1330      do lo = los, loe                                         
1331         do li = lis, lie
1332            level = "=" + data&z(lo) + "m"       
1333            vecres                  = True            ; vector only resources
1334            vecres@gsnDraw          = False           ; don't draw
1335            vecres@gsnFrame         = False           ; don't advance frame
1336            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1337            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1338            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1339            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1340            vecres@gsnLeftString = "t=" + data&t(li) +"s  z"+level
1341            vecres@tiXAxisString    = " "                                               
1342            plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1343            n=n+1
1344         end do
1345      end do
1346   end if
1347
1348   check = True
1349 
1350   do varn=dim-1,0,1
1351
1352      if (vector .EQ. 1) then   
1353         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1354      end if
1355
1356      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
1357         check = False
1358      else
1359         check = True
1360      end if   
1361      if (  isvar("var") ) then 
1362         check = isStrSubset( var,","+vNam(varn)+",")
1363      end if
1364      if (parameter(21) .NE. "variables") then
1365         var = parameter(21)
1366         check = isStrSubset( var,","+vNam(varn)+"," )
1367      end if
1368   
1369      if(check) then
1370         print(vNam(varn))
1371         space=(MaxVal(varn)-MinVal(varn))/24
1372 
1373         cs_res@cnMinLevelValF = MinVal(varn)
1374         cs_res@cnMaxLevelValF = MaxVal(varn)
1375
1376         cs_res@cnLevelSpacingF = space
1377     
1378         ; ****************************************************
1379         ; loops over time and layer
1380         ; ****************************************************
1381
1382         
1383         do lo = los, loe                                       
1384            do li = lis, lie
1385
1386               ; ****************************************************
1387               ; xy cross section
1388               ; ****************************************************
1389
1390               if ( xyc .eq. 1 ) then
1391               
1392                  cs_res@tiXAxisString = "x [m]"
1393                  cs_res@tiYAxisString = "y [m]"
1394                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1395                 
1396                  if ( sort .eq. "time" ) then
1397                     if ( data&z(lo) .eq. -1 ) then
1398                        level = "-average"
1399                     else
1400                        level = "=" + data&z(lo) + "m"
1401                     end if
1402                     cs_res@gsnCenterString = "t=" + data&t(li) +"s  z"+level
1403                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
1404                     if (vector .EQ. 1 .AND. check_vecp) then
1405                        vecres                  = True            ; vector only resources
1406                        vecres@gsnDraw          = False           ; don't draw
1407                        vecres@gsnFrame         = False           ; don't advance frame
1408                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1409                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1410                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
1411                        vecres@gsnRightString   = " "
1412                        vecres@gsnLeftString    = " "
1413                        vecres@tiXAxisString    = " "   
1414                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1415                        overlay(plot(n), plot_vec)
1416                     end if                         
1417                  end if
1418                 
1419                  if ( sort .eq. "layer" ) then
1420                     if ( data&z(li) .eq. -1 ) then
1421                        level = "-average"
1422                     else
1423                        level = "=" + data&z(li) + "m"
1424                     end if
1425                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  z"+ level
1426                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
1427                     if (vector .EQ. 1 .AND. check_vecp) then
1428                        vecres                  = True            ; vector only resources
1429                        vecres@gsnDraw          = False           ; don't draw
1430                        vecres@gsnFrame         = False           ; don't advance frame
1431                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1432                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1433                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1434                        vecres@gsnRightString   = " "             ; turn off right string
1435                        vecres@gsnLeftString    = " "             ; turn off left string
1436                        vecres@tiXAxisString    = " "   
1437                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1438                        overlay(plot(n), plot_vec)
1439                     end if
1440                  end if
1441               end if
1442
1443               ; ****************************************************
1444               ; xz cross section
1445               ; ****************************************************
1446
1447               if ( xzc .eq. 1 ) then
1448       
1449                  cs_res@tiXAxisString   = "x [m]"
1450                  cs_res@tiYAxisString   = "z [m]"
1451                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1452               
1453                  if ( sort .eq. "time" ) then
1454                     if ( data&z(lo) .eq. -1 ) then
1455                        level = "-average"
1456                     else
1457                        level = "=" + data&y(lo) + "m"
1458                     end if
1459                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  y"+ level
1460                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
1461                     if (vector .EQ. 1 .AND. check_vecp) then
1462                        vecres                  = True            ; vector only resources
1463                        vecres@gsnDraw          = False           ; don't draw
1464                        vecres@gsnFrame         = False           ; don't advance frame
1465                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1466                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1467                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1468                        vecres@gsnRightString   = " "             ; turn off right string
1469                        vecres@gsnLeftString    = " "             ; turn off left string
1470                        vecres@tiXAxisString    = " "   
1471                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1472                        overlay(plot(n), plot_vec)
1473                     end if
1474                  end if
1475
1476                  if ( sort .eq. "layer" ) then
1477                     if ( data&z(li) .eq. -1 ) then
1478                        level = "-average"
1479                     else
1480                        level = "=" + data&y(li) + "m"
1481                     end if
1482                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  y"+ level
1483                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
1484                     if (vector .EQ. 1 .AND. check_vecp) then
1485                        vecres                  = True            ; vector only resources
1486                        vecres@gsnDraw          = False           ; don't draw
1487                        vecres@gsnFrame         = False           ; don't advance frame
1488                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1489                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1490                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1491                        vecres@gsnRightString   = " "             ; turn off right string
1492                        vecres@gsnLeftString    = " "             ; turn off left string
1493                        vecres@tiXAxisString    = " "   
1494                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1495                        overlay(plot(n), plot_vec)
1496                     end if
1497                  end if                 
1498               end if
1499
1500               ; ****************************************************
1501               ; yz cross section
1502               ; ****************************************************
1503
1504               if ( yzc .eq. 1 ) then
1505               
1506                  cs_res@tiXAxisString   = "y [m]"
1507                  cs_res@tiYAxisString   = "z [m]"
1508                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1509
1510                  if ( sort .eq. "time" ) then
1511                     if ( data&z(lo) .eq. -1 ) then
1512                        level = "-average"
1513                     else
1514                        level = "=" + data&z(lo) + "m"
1515                     end if
1516                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  x"+ level
1517                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
1518                     if (vector .EQ. 1 .AND. check_vecp) then
1519                        vecres                  = True            ; vector only resources
1520                        vecres@gsnDraw          = False           ; don't draw
1521                        vecres@gsnFrame         = False           ; don't advance frame
1522                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1523                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1524                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1525                        vecres@gsnRightString   = " "             ; turn off right string
1526                        vecres@gsnLeftString    = " "             ; turn off left string
1527                        vecres@tiXAxisString    = " "   
1528                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1529                        overlay(plot(n), plot_vec)
1530                     end if
1531                  end if
1532
1533                  if ( sort .eq. "layer" ) then
1534                     if ( data&z(li) .eq. -1 ) then
1535                        level = "-average"
1536                     else
1537                        level = "=" + data&z(li) + "m"
1538                     end if
1539                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  x"+ level
1540                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
1541                     if (vector .EQ. 1 .AND. check_vecp)then
1542                        vecres                  = True            ; vector only resources
1543                        vecres@gsnDraw          = False           ; don't draw
1544                        vecres@gsnFrame         = False           ; don't advance frame
1545                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1546                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1547                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1548                        vecres@gsnRightString   = " "             ; turn off right string
1549                        vecres@gsnLeftString    = " "             ; turn off left string
1550                        vecres@tiXAxisString    = " "   
1551                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1552                        overlay(plot(n), plot_vec)
1553                     end if
1554                  end if
1555               end if         
1556               n=n+1   
1557            end do     
1558         end do 
1559      end if
1560   end do
1561
1562   ; ***************************************************
1563   ; merge plots onto one page
1564   ; ***************************************************
1565
1566   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec")then
1567      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1568         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1569         print(" ")
1570         print("Outputs to .eps or .epsi have only one frame")
1571         print(" ")
1572      else
1573         do np = 0,no_layer*no_time*(no_var+1)-1,no_lines*no_columns   
1574            if ( np + no_lines*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then
1575               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_lines,no_columns/),cs_resP)
1576            else
1577               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
1578            end if
1579         end do
1580      end if
1581   else         
1582      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1583         gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
1584         print(" ")
1585         print("Outputs to .eps or .epsi have only one frame")
1586         print(" ")
1587      else
1588         do np = 0,no_layer*no_time*no_var-1,no_lines*no_columns   
1589            if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then
1590               gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_lines,no_columns/),cs_resP)
1591            else
1592               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
1593            end if
1594         end do
1595      end if
1596   end if
1597
1598   print(" ")
1599   print("Output to: " + file_out +"."+ format_out)
1600   print(" ")   
1601 
1602end
Note: See TracBrowser for help on using the repository browser.