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

Last change on this file since 154 was 154, checked in by letzel, 14 years ago
  • new NCL scripts added to trunk/SCRIPTS/NCL
  • previous NCL scripts moved to trunk/SCRIPTS/NCL/Archive
File size: 24.2 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
5
6begin
7
8   ; ***************************************************
9   ; read parameter_list
10   ; ***************************************************
11
12   if (isfilepresent("~/.ncl_preferences")) then
13      parameter = asciiread("~/.ncl_preferences",63,"string")
14      delete(parameter@_FillValue)
15   else
16      print(" ")
17      print("Please copy '.ncl_preferences' into your $home dircetory")
18      print(" ")
19      exit
20   end if
21
22   ; ***************************************************
23   ; set default parameter values and strings if not assigned in prompt or parameter list
24   ; ***************************************************
25   
26   if ( .not. isvar("file_in") ) then                   ; path+name of input file     
27      if (parameter(7) .EQ. "input file") then
28         print(" ")
29         print("Please provide input file 'file_in = ' either in prompt or parameter_list")
30         print(" ")
31         exit
32      else
33         file_in = parameter(7)
34      end if     
35   end if
36   if ( .not. isvar("format_out") ) then                ; format of output file
37      format_out = "x11"
38      if (parameter(9) .NE. "x11") then
39         format_out = parameter(9) 
40      end if
41   end if
42   if ( .not. isvar("file_out") ) then                  ; path+name of output file
43      file_out = "test"
44      if (parameter(11) .NE. "test") then
45         file_out = parameter(11) 
46     end if     
47   end if
48   if ( .not. isvar("no_columns") ) then                ; number of plots in one row
49      no_columns = 1
50      if (parameter(17) .NE. "1") then
51         no_columns = stringtointeger(parameter(17)) 
52      end if
53   end if
54   if ( .not. isvar("no_lines") ) then                  ; number of plot-lines on one sheet
55      no_lines = 2
56      if (parameter(19) .NE. "2") then
57         no_lines = stringtointeger(parameter(19)) 
58      end if
59   end if
60   if ( .not. isvar("sort") ) then                      ; sort of plots
61      sort = "time"
62      if (parameter(37) .NE. "time") then
63         sort = parameter(37) 
64      end if
65   end if
66   if ( .not. isvar("mode") ) then                      ; pattern of contour plots
67      mode = "Fill"
68      if (parameter(39) .NE. "Fill") then
69         mode = parameter(39) 
70      end if
71   end if
72   if ( .not. isvar("fill_mode") ) then                 ; pattern of filling
73      fill_mode = "AreaFill"
74      if (parameter(41) .NE. "AreaFill") then
75         mode = parameter(41) 
76      end if
77   end if
78   if ( .not. isvar("shape") ) then                     ; keeping of aspect ratio
79      shape = 1
80      if (parameter(43) .NE. "1") then
81         shape = stringtointeger(parameter(43))
82         if (stringtointeger(parameter(43)) .NE. 0) then
83            print(" ")
84            print("Please set 'shape' to 0 or 1")
85            print(" ")
86            exit
87         end if
88      end if 
89   end if 
90   if ( .not. isvar("var") ) then                       ; set output of all variables
91      check = True
92   end if
93
94   if ( .not. isvar("xyc") ) then                       ; turn xy cross-section on or off
95      xyc = 0
96      if (parameter(45) .NE. "0") then
97         xyc = stringtointeger(parameter(45))
98         if (stringtointeger(parameter(45)) .NE. 1) then
99            print(" ")
100            print("Please set 'xyc' to 0 or 1")
101            print(" ")
102            exit
103         end if
104      end if
105   end if
106   if ( .not. isvar("xzc") ) then                       ; turn xz cross-section on or off
107      xzc = 0
108      if (parameter(47) .NE. "0") then
109         xzc = stringtointeger(parameter(47))
110         if (stringtointeger(parameter(47)) .NE. 1) then
111            print(" ")
112            print("Please set 'xzc' to 0 or 1")
113            print(" ")
114            exit
115         end if
116      end if
117   end if
118   if ( .not. isvar("yzc") ) then                       ; turn yz cross-section on or off
119      yzc = 0
120      if (parameter(49) .NE. "0") then
121         yzc = stringtointeger(parameter(49))
122         if (stringtointeger(parameter(49)) .NE. 1) then
123            print(" ")
124            print("Please set 'yzc' to 0 or 1")
125            print(" ")
126            exit
127         end if
128      end if
129   end if
130   if ( .not. isvar("xs") ) then                        ; start of x-coordinate range
131      xs = 0
132      if (parameter(51) .NE. "0") then
133         xs = stringtointeger(parameter(51))
134      end if
135   end if
136   if ( .not. isvar("ys") ) then                        ; start of y-coordinate range
137      ys = 0
138      if (parameter(55) .NE. "0") then
139         ys = stringtointeger(parameter(55))
140      end if
141   end if
142   if ( .not. isvar("zs") ) then                        ; start of z-coordinate range
143      zs = 0
144      if (parameter(59) .NE. "0") then
145         zs = stringtointeger(parameter(59))
146      end if
147   end if 
148
149   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
150      print(" ")
151      print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)")
152      print(" ")
153      exit
154   end if
155   if (xyc .EQ. 1 ) then
156      if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
157         print(" ")
158         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
159         print(" ")
160         exit
161      end if
162   end if
163   if (xzc .EQ. 1) then
164      if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
165         print(" ")
166         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
167         print(" ")
168         exit
169      end if
170   end if
171   if (yzc .EQ. 1) then
172      if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
173         print(" ")
174         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
175         print(" ")
176         exit
177      end if
178   end if
179
180   ; ***************************************************
181   ; open input file
182   ; ***************************************************
183
184   f  = addfile( file_in, "r" )
185
186   vNam  = getfilevarnames(f)
187   print(" ")
188   print("Variable on netCDF file: " + vNam)
189   print(" ")
190   dim   = dimsizes(vNam)
191 
192   ; ***************************************************
193   ; set up recourses
194   ; ***************************************************
195
196   cs_res                          = True
197   cs_res@gsnYAxisIrregular2Linear = True 
198
199   if( shape .eq. 1 ) then                              ; keep the shape
200      cs_res@gsnShape = True
201   end if
202 
203   cs_res@gsnDraw                 = False
204   cs_res@gsnFrame                = False
205   cs_res@gsnMaximize             = True
206   cs_res@gsnPaperOrientation     = "portrait"
207   cs_res@gsnPaperWidth           = 8.27
208   cs_res@gsnPaperHeight          = 11.69
209   cs_res@gsnPaperMargin          = 0.79
210   cs_res@tmXBLabelFontHeightF   = .02
211   cs_res@tmYLLabelFontHeightF   = .02
212   cs_res@tiXAxisFontHeightF     = .02
213   cs_res@tiYAxisFontHeightF     = .02
214   cs_res@tmXBMode                ="Automatic"
215   cs_res@tmYLMode                ="Automatic"
216   cs_res@lgTitleFontHeightF     = .02
217   cs_res@lgLabelFontHeightF     = .02
218   cs_res@txFontHeightF          = .02
219   cs_resP = True
220   cs_resP@txString               = f@title
221   cs_resP@txFuncCode             = "~"                 ; necessary for title
222   cs_resP@txFontHeightF          = .02
223 
224   if ( mode .eq. "Fill" ) then
225      cs_res@cnFillOn             = True
226      cs_res@gsnSpreadColors      = True
227      cs_res@cnFillMode           = fill_mode
228      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
229      cs_res@cnLinesOn            = False       
230      cs_res@cnLineLabelsOn       = False
231   end if
232
233   if ( mode .eq. "Both" ) then
234      cs_res@cnFillOn             = True
235      cs_res@gsnSpreadColors      = True
236      cs_res@cnFillMode           = fill_mode
237      cs_res@lbOrientation        = "Vertical"          ; vertical label bar
238      cs_res@cnLinesOn            = True
239      cs_res@cnLineLabelsOn       = True
240   end if
241
242   ; ***************************************************
243   ; open workstation(s)
244   ; ***************************************************
245
246   wks_ps  = gsn_open_wks(format_out,file_out)
247   gsn_define_colormap(wks_ps,"rainbow+white")
248
249   plot=new((/no_columns*no_lines/),graphic)
250
251   page = 0
252
253   nc1 = no_columns
254   nl1 = no_lines
255
256   print("plots sorted by '"+sort+"'")
257   print(" ")
258
259   do varn=0,dim-1   
260
261      data_all = f->$vNam(varn)$
262
263      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
264         check = False
265      end if
266      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
267         check = False
268      end if
269     
270      if (  isvar("var") ) then 
271         check = isStrSubset( var,vNam(varn)+",")
272      end if
273      if (parameter(21) .NE. "variables") then
274         var = parameter(21)
275         check = isStrSubset( var,vNam(varn)+"," )
276      end if
277
278      if(check) then
279
280         print("plot of "+vNam(varn))
281         print(" ")
282
283         ; *********************************************
284         ; set up of start_time_step and end_time_step
285         ; *********************************************
286
287         t = f->time
288         nt = dimsizes(t)
289         
290         ; *********************************************       
291         ; start of time step and different types of mistakes that could be done
292         ; *********************************************
293
294         if ( .not. isvar("start_time_step") ) then     
295            start_time_step = 1
296            if (parameter(13) .NE. "1") then
297               if (parameter(13) .LE. "0")
298                  print(" ")
299                  print("Begin with time step 1")
300                  print(" ")
301                  exit
302               end if
303               if (stringtointeger(parameter(13)) .GE. nt)
304                  print(" ")
305                  print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
306                  print(" ")
307                  exit
308               end if
309               start_time_step = stringtointeger(parameter(13)) 
310            end if
311         else
312            if (start_time_step .LE. 0)
313               print(" ")
314               print("Begin with time step 1")
315               print(" ")
316               exit
317            end if
318            if (start_time_step .GE. nt)
319               print(" ")
320               print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
321               print(" ")
322               exit
323            end if
324         end if
325
326         ; ****************************************************
327         ; end of time step and different types of mistakes that could be done
328         ; ****************************************************
329
330         if ( .not. isvar("end_time_step") ) then       
331            end_time_step = nt-1
332            if (parameter(15) .NE. "nt-1") then
333               if (parameter(15) .LE. "0")
334                  print(" ")
335                  print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
336                  print(" ")
337                  exit
338               end if
339               if (stringtointeger(parameter(15)) .GE. nt)
340                  print(" ")
341                  print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1))
342                  print(" ")
343                  exit
344               end if
345               if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
346                  print(" ")
347                  print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(15))
348                  print(" ")
349                  exit
350               end if
351               end_time_step = stringtointeger(parameter(15)) 
352            end if   
353         else
354            if (end_time_step .LE. 0)
355               print(" ")
356               print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
357               print(" ")
358               exit
359            end if
360            if (end_time_step .GE. nt)
361               print(" ")
362               print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
363               print(" ")
364               exit
365            end if
366            if (end_time_step .LT. start_time_step)
367               print(" ")
368               print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
369               print(" ")
370               exit
371            end if
372         end if
373
374         ; ****************************************************
375         ; set up ranges of x-, y- and z-coordinates
376         ; ****************************************************         
377                   
378         xdim = dimsizes(data_all(0,0,0,:)) - 1   
379         ydim = dimsizes(data_all(0,0,:,0)) - 1       
380         zdim = dimsizes(data_all(0,:,0,0)) - 1   
381   
382         if ( .not. isvar("xe")) then     ; output x-coordinate range end (in index)
383            xe = xdim
384            if (parameter(53) .NE. "xdim") then
385               if (stringtointeger(parameter(53)) .GT. xdim) then
386                  print(" ")
387                  print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
388                  print(" ")
389                  exit
390               end if
391               if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
392                  print(" ")
393                  print("range end for x-coordinate = "+parameter(53)+" is too small")
394                  print(" ")
395                  exit
396               end if
397               xe = stringtointeger(parameter(53))
398            end if
399         else
400           if (xe .GT. xdim) then
401              print(" ")
402              print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
403              print(" ")
404              exit
405           end if
406           if (xe .LT. 0 .OR. xe .LT. xs) then
407              print(" ")
408              print("range end for x-coordinate = "+xe+" is too small")
409              print(" ")
410              exit
411           end if                 
412         end if
413         if ( .not. isvar("ye")) then     ; output y-coordinate range end (in index)
414            ye = ydim
415            if (parameter(57) .NE. "ydim") then
416               if (stringtointeger(parameter(57)) .GT. ydim) then
417                  print(" ")
418                  print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
419                  print(" ")
420                  exit
421               end if
422               if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
423                  print(" ")
424                  print("range end for y-coordinate = "+parameter(57)+" is too small")
425                  print(" ")
426                  exit
427               end if
428               ye = stringtointeger(parameter(57))
429            end if
430         else
431            if (ye .GT. ydim) then
432               print(" ")
433               print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
434               print(" ")
435               exit
436            end if
437            if (ye .LT. 0 .OR. ye .LT. ys) then
438               print(" ")
439               print("range end for y-coordinate = "+ye+" is too small")
440               print(" ")
441               exit
442            end if
443         end if
444         if ( .not. isvar("ze")) then     ; output z-coordinate range end (in index)
445            ze = zdim
446            if (parameter(61) .NE. "zdim") then
447               if (stringtointeger(parameter(61)) .GT. zdim) then
448                  print(" ")
449                  print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
450                  print(" ")
451                  exit
452               end if
453               if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
454                  print(" ")
455                  print("range end for z-coordinate = "+parameter(61)+" is too small")
456                  print(" ")
457                  exit
458               end if
459               ze = stringtointeger(parameter(61))
460            end if
461         else
462            if (ze .GT. zdim) then
463               print(" ")
464               print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
465               print(" ")
466               exit
467            end if
468            if (ze .LT. 0 .OR. ze .LT. zs) then
469               print(" ")
470               print("range end for z-coordinate = "+ze+" is too small")
471               print(" ")
472               exit
473            end if
474         end if
475 
476         ; Defining inner and outer loop depending on variable
477         ; "sort" (values: time or layers)
478   
479         if ( xyc .eq. 1 ) then
480            lays = zs
481            laye = ze
482         end if
483         if ( xzc .eq. 1 ) then
484            lays = ys
485            laye = ye
486         end if
487         if ( yzc .eq. 1) then
488            lays = xs
489            laye = xe
490         end if
491
492         if ( sort .eq. "time" ) then
493            lis = start_time_step
494            lie = end_time_step
495            los = lays
496            loe = laye
497         else
498            lis = lays
499            lie = laye
500            los = start_time_step
501            loe = end_time_step
502         end if
503
504
505         ; loops over time and layer (order depends on "sort")
506
507         do lo = los, loe                                       ; lo: loop outer
508            do li = lis, lie, no_lines*no_columns               ; li: loop inner
509
510               print("page "+page)
511               page = page + 1
512
513               ; Defining indices for reading data needed to plot one page
514               ; The indices depend on the parameter "sort"
515
516               if ( sort .eq. "time" ) then
517                  ta = li
518                  tb = min( (/li + no_lines*no_columns - 1, end_time_step/) )
519                  if ( xyc .eq. 1 ) then
520                     xa = xs
521                     xb = xe
522                     ya = ys
523                     yb = ye
524                     za = lo
525                     zb = lo
526                  end if
527                  if ( xzc .eq. 1 ) then
528                     xa = xs
529                     xb = xe
530                     ya = lo
531                     yb = lo
532                     za = zs
533                     zb = ze
534                  end if
535                  if ( yzc .eq. 1 ) then
536                     xa = lo
537                     xb = lo
538                     ya = ys
539                     yb = ye
540                     za = zs
541                     zb = ze
542                  end if
543               end if
544
545               if ( sort .eq. "layer" ) then
546                  ta = lo
547                  tb = lo
548                  if ( xyc .eq. 1 ) then
549                     xa = xs
550                     xb = xe
551                     ya = ys
552                     yb = ye
553                     za = li
554                     zb = min( (/li + no_lines*no_columns - 1, ze/) ) 
555                  end if
556                  if ( xzc .eq. 1 ) then
557                     xa = xs
558                     xb = xe
559                     ya = li
560                     yb = min( (/li + no_lines*no_columns - 1, ye/) )
561                     za = zs
562                     zb = ze
563                  end if
564                  if ( yzc .eq. 1 ) then
565                     xa = li
566                     xb = min( (/li + no_lines*no_columns - 1, xe/) )
567                     ya = ys
568                     yb = ye
569                     za = zs
570                     zb = ze
571                  end if
572               end if
573                           
574               data = data_all(ta:tb,za:zb,ya:yb,xa:xb)
575       
576               data!0 = "t"
577               data!1 = "z"
578               data!2 = "y"
579               data!3 = "x"
580
581               ; set up labels and plot
582
583               ; xy cross section
584
585               if ( xyc .eq. 1 ) then
586                  cs_res@tiXAxisString   = "x [m]"
587                  cs_res@tiYAxisString   = "y [m]"
588
589                  if ( sort .eq. "time" ) then
590                     if ( data&z(0) .eq. -1 ) then
591                        level = "-average"
592                     else
593                        level = "=" + data&z(0) + "m"
594                     end if
595                     do n = 0, tb-ta
596                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  z"+ level
597                        plot(n) = gsn_csm_contour(wks_ps,data(n,0,:,:),cs_res) 
598                     end do
599                  end if
600
601                  if ( sort .eq. "layer" ) then
602                     do n = 0, zb-za
603                        if ( data&z(n) .eq. -1 ) then
604                           level = "-average"
605                        else
606                           level = "=" + data&z(n) + "m"
607                        end if
608                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  z"+ level
609                        plot(n) = gsn_csm_contour(wks_ps,data(0,n,:,:),cs_res)
610                     end do
611                  end if
612               end if
613
614               ; xz cross section
615
616               if ( xzc .eq. 1 ) then
617                  cs_res@tiXAxisString   = "x [m]"
618                  cs_res@tiYAxisString   = "z [m]"
619 
620                  if ( sort .eq. "time" ) then
621                     if ( data&y(0) .eq. -1 ) then
622                        level = "-average"
623                     else
624                        level = "=" + data&y(0) + "m"
625                     end if
626                     do n = 0, tb-ta
627                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  z"+ level
628                        plot(n) = gsn_csm_contour(wks_ps,data(n,:,0,:),cs_res) 
629                     end do
630                  end if
631
632                  if ( sort .eq. "layer" ) then
633                     do n = 0, yb-ya
634                        if ( data&y(n) .eq. -1 ) then
635                           level = "-average"
636                        else
637                           level = "=" + data&y(n) + "m"
638                        end if
639                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  z"+ level
640                        plot(n) = gsn_csm_contour(wks_ps,data(0,:,n,:),cs_res)
641                     end do
642                  end if
643               end if
644
645               ; yz cross section
646
647               if ( yzc .eq. 1 ) then
648                  cs_res@tiXAxisString   = "y [m]"
649                  cs_res@tiYAxisString   = "z [m]"
650 
651                  if ( sort .eq. "time" ) then
652                     if ( data&x(0) .eq. -1 ) then
653                        level = "-average"
654                     else
655                        level = "=" + data&x(0) + "m"
656                     end if
657                     do n = 0, tb-ta
658                        cs_res@gsnCenterString = "t=" + data&t(n) + "s  x"+ level
659                        plot(n) = gsn_csm_contour(wks_ps,data(n,:,:,0),cs_res) 
660                     end do
661                  end if
662
663                  if ( sort .eq. "layer" ) then
664                     do n = 0, xb-xa
665                        if ( data&x(n) .eq. -1 ) then
666                           level = "-average"
667                        else
668                           level = "=" + data&x(n) + "m"
669                        end if
670                        cs_res@gsnCenterString = "t=" + data&t(0) + "s  x"+ level
671                        plot(n) = gsn_csm_contour(wks_ps,data(0,:,:,n),cs_res)
672                     end do
673                  end if
674               end if
675
676               delete(data)
677
678               ; ***************************************************
679               ; merge plots onto one page
680               ; ***************************************************
681         
682               gsn_panel(wks_ps, plot(0:n-1),(/no_lines,no_columns/),cs_resP)
683
684               ; To hold no_lines and no_columns constant, it must be defined again:
685
686               no_lines   = nl1
687               no_columns = nc1
688
689            end do     
690         end do 
691      end if 
692 
693      delete(data_all)
694
695   end do   
696
697   print(" ")
698   print("Output to: " + file_out +"."+ format_out)
699   print(" ")   
700 
701end
Note: See TracBrowser for help on using the repository browser.