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

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