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

Last change on this file since 183 was 175, checked in by steinfeld, 16 years ago

NCL scripts for spectra added. Bug fix concerning spectra in netcdf.f90

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