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

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