Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Dec 10, 2008 9:14:34 AM (15 years ago)
Author:
letzel
Message:
  • User can check user parameters and deduce further quantities in user_check_parameters
  • Topography definition according to new user parameter topography_grid_convention (init_grid, modules, user_header, user_init_grid, user_parin)
  • New subdirectory trunk/EXAMPLES with two examples highlighting the topography_grid_convention
File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SCRIPTS/NCL/spectra.ncl

    r194 r218  
    2626   end if
    2727end if
    28    
     28
    2929begin
    3030
    … …  
    3535   if (file_1 .EQ. "File in") then
    3636      print(" ")
    37       print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'")
     37      print("Declare input file 'file_1=' in 'ncl_preferences.ncl' or prompt")
    3838      print(" ")
    3939      exit
    … …  
    4141      file_in = file_1   
    4242   end if
    43    if (.not. isfilepresent(file_in)) then
    44       print(" ")
    45       print("1st input file: '"+file_in+"' does not exist")
    46       print(" ")
    47       exit
    48    end if
    4943
    5044   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
    … …  
    5549   end if
    5650   
    57    if (logx .NE. 0 .AND. logx .NE. 1)then
    58       print(" ")
    59       print("'logx'= "+logx+" is invalid and set to 1")
    60       print(" ")
    61       logx = 1
     51   if (log_x .NE. 0 .AND. log_x .NE. 1)then
     52      print(" ")
     53      print("'log_x'= "+log_x+" is invalid and set to 1")
     54      print(" ")
     55      log_x = 1
    6256   end if 
    6357   
    64    if (logy .NE. 0 .AND. logy .NE. 1)then
    65       print(" ")
    66       print("'logy'= "+logy+" is invalid and set to 1")
    67       print(" ")
    68       logy = 1
     58   if (log_y .NE. 0 .AND. log_y .NE. 1)then
     59      print(" ")
     60      print("'log_y'= "+log_y+" is invalid and set to 1")
     61      print(" ")
     62      log_y = 1
    6963   end if   
    7064 
    71    if (normy .EQ. 0.) then
    72       print(" ")
    73       print("You cannot normalise the y-axis with 0, 'normy' is set to 1.0")
    74       print(" ")
    75       normy = 1.0
    76    end if
    77    if (normx .EQ. 0.) then
    78       print(" ")
    79       print("You cannot normalise the x-axis with 0, 'normx' is set to 1.0")
    80       print(" ")
    81       normx= 1.0
     65   if (norm_y .EQ. 0.) then
     66      print(" ")
     67      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
     68      print(" ")
     69      norm_y = 1.0
     70   end if
     71   if (norm_x .EQ. 0.) then
     72      print(" ")
     73      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
     74      print(" ")
     75      norm_x= 1.0
    8276   end if
    8377   
    … …  
    106100   ; open input file
    107101   ;***************************************************
    108 
    109    f=addfile(file_in,"r")
    110    
    111    vNam = getfilevarnames(f)
     102   
     103   file_in_1 = False
     104   if (isStrSubset(file_in, ".nc"))then
     105      start_f = -2
     106      end_f = -2
     107      file_in_1 = True     
     108   end if
     109
     110   if (start_f .EQ. -1)then
     111      print(" ")
     112      print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
     113      print(" ") 
     114      exit
     115   end if
     116   if (end_f .EQ. -1)then
     117      print(" ")
     118      print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
     119      print(" ") 
     120      exit
     121   end if
     122
     123   files=new(end_f-start_f+1,string)
     124   
     125   if (file_in_1)then
     126      if (isfilepresent(file_in))then
     127         files(0)=file_in
     128      else
     129         print(" ")
     130         print("1st input file: '"+file_in+"' does not exist")
     131         print(" ")
     132         exit
     133      end if
     134   else   
     135      if (start_f .EQ. 0)then
     136         if (isfilepresent(file_in+".nc"))then
     137            files(0)=file_in+".nc"
     138            do i=1,end_f
     139               if (isfilepresent(file_in+"."+i+".nc"))then   
     140                  files(i)=file_in+"."+i+".nc"
     141               else
     142                  print(" ")
     143                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
     144                  print(" ")
     145                  exit 
     146               end if       
     147            end do         
     148         else
     149            print(" ")
     150            print("Input file: '"+file_in+".nc' does not exist")
     151            print(" ")
     152            exit
     153         end if
     154      else
     155         do i=start_f,end_f
     156            if (isfilepresent(file_in+"."+i+".nc"))then   
     157               files(i-start_f)=file_in+"."+i+".nc"
     158            else
     159               print(" ")
     160               print("Input file: '"+file_in+"."+i+".nc' does not exist")
     161               print(" ")
     162               exit 
     163            end if
     164         end do
     165      end if
     166   end if
     167
     168   f=addfiles(files,"r")
     169   f_att=addfile(files(0),"r")
     170   ListSetType(f,"cat")
     171   
     172   vNam = getfilevarnames(f_att)
    112173   print(" ")
    113174   print("Variables in input file:")
    … …  
    115176   print(" ")
    116177   dim = dimsizes(vNam)
    117    vDim = getfiledimsizes(f)
     178   vDim = getfiledimsizes(f_att)
    118179 
    119    t_all = f->time
     180   t_all = f[:]->time
    120181   nt    = dimsizes(t_all)
    121182   delta_t=t_all(nt-1)/nt
    122183   
    123    k_x=f->k_x
     184   k_x=f_att->k_x
    124185   dimx=dimsizes(k_x)
    125    k_y=f->k_y
     186   k_y=f_att->k_y
    126187   dimy=dimsizes(k_y)
    127188   
    … …  
    131192   do i=0,dim-1
    132193      if (vNam(i) .EQ. "zu_sp")then
    133          zu=f->zu_sp     
     194         zu=f_att->zu_sp         
    134195         if (height_level(0) .EQ. -1)then
    135196            dimz=dimsizes(zu)
    … …  
    150211      else
    151212         if (vNam(i) .EQ. "zw_sp")then
    152             zw=f->zw_sp
     213            zw=f_att->zw_sp
    153214            if (height_level(0) .EQ. -1)then             
    154215               dimz=dimsizes(zw)
    … …  
    182243         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")
    183244         print(" ")
    184          print("Please select another 'start_time_step'")
     245         print("Select another 'start_time_step'")
    185246         print(" ")
    186247         exit
    … …  
    204265      print("'start_time_step' = "+ start_time_step +"h is invalid")
    205266      print(" ")
    206       print("Please select another 'start_time_step'")
     267      print("Select another 'start_time_step'")
    207268      print(" ")
    208269      exit
    … …  
    220281         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")
    221282         print(" ")
    222          print("Please select another 'end_time_step'") 
     283         print("Select another 'end_time_step'") 
    223284         print(" ")
    224285         exit
    … …  
    228289         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
    229290         print(" ")
    230          print("Please select another 'start_time_step' or 'end_time_step'")
     291         print("Select another 'start_time_step' or 'end_time_step'")
    231292         print(" ")
    232293         exit
    … …  
    245306      print("'end_time_step' = "+ end_time_step +"h is invalid")
    246307      print(" ")
    247       print("Please select another 'end_time_step'")
     308      print("Select another 'end_time_step'")
    248309      print(" ")
    249310      exit
    … …  
    281342   res@lgLabelFont             = "helvetica"
    282343   res@tmLabelAutoStride       = True
    283    res@pmLegendDisplayMode     = "Always"
    284    res@pmLegendSide            = "Top"
    285    res@pmLegendParallelPosF    = 1.2
    286    res@pmLegendOrthogonalPosF  = -1.0
    287    res@pmLegendWidthF          = 0.12
    288    if (sort .EQ. "time")then
    289       res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
    290    else
    291       res@pmLegendHeightF         = 0.04*dimz
    292    end if
    293    res@lgLabelFontHeightF     = .025
    294    res@lgTitleFontHeightF     = .025
    295    res@txFontHeightF      = 0.025
    296    res@tiXAxisFontHeightF = 0.025
    297    res@tiYAxisFontHeightF = 0.025
    298    
    299    if (logx .EQ. 1) then
     344   
     345   res@lgLabelFontHeightF     = font_size_legend
     346   res@lgTitleFontHeightF     = font_size
     347   res@txFontHeightF      = font_size
     348   res@tiXAxisFontHeightF = font_size
     349   res@tiYAxisFontHeightF = font_size
     350   res@tmXBLabelFontHeightF = font_size
     351   res@tmYLLabelFontHeightF = font_size
     352   
     353   res@tmXBMinorPerMajor = 4
     354   res@tmYLMinorPerMajor = 4
     355   
     356   if (log_x .EQ. 1) then
    300357      res@trXLog = True
    301358   else
    302359      res@trXLog = False
    303360   end if
    304    if (logy .EQ. 1)then
     361   if (log_y .EQ. 1)then
    305362      res@trYLog = True
    306363   else 
    … …  
    326383      plot  = new(dim*dimt,graphic)
    327384      np=dimz
    328       res@lgTitleString = "Height [m]" 
     385     
    329386      do p=0,dimz-1
    330387         if (height_level(0) .EQ. -1)then
    … …  
    340397   if (black .eq. 0 ) then
    341398      if (np .EQ. 1)then
    342          res@xyLineColors = 237
     399         color = 237
    343400      else   
    344401         step=round(235/(np-1),3)
    345          res@xyLineColors = ispan(2,237,step)
     402         color = ispan(2,237,step)
    346403      end if
    347404   end if
    … …  
    352409   wks=gsn_open_wks(format_out,file_out)
    353410   gsn_define_colormap(wks,"rainbow+white")
    354 
    355    temp=new((/dimt,dimz,dimx/),float)
    356411
    357412   n=0
    … …  
    369424
    370425      if(check) then
    371      
    372          temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:)
    373            
    374          a=getvardims(temp)
     426
     427         temp = f[:]->$vNam(varn)$
     428         data = temp(start_time_step:end_time_step,0:dimz-1,:)
     429
     430         temp_att = f_att->$vNam(varn)$
     431         a=getvardims(temp_att)
    375432         b=dimsizes(a)
    376433
    377434         if (height_level(0) .NE. -1)then
    378435            do te=0,dimz-1
    379                temp(:,te,:) = f->$vNam(varn)$(start_time_step:end_time_step,height_level(te),:)       
     436            print(height_level(te))
     437               data(:,te,:) = temp(start_time_step:end_time_step,height_level(te),:)           
    380438            end do
    381439         end if
    382440
    383          temp=temp/(normy*normx)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
     441         data=data/(norm_y*norm_x)  ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)#
    384442           
    385443         do i=0,b-1           
    … …  
    402460            end if
    403461         end do 
    404 
     462         
     463         min_x=new(dimz,double)
     464         max_x=new(dimz,double)
     465         min_y=new(dimz,float)
     466         max_y=new(dimz,float)
     467         plot_h  = new(dimz,graphic)
     468         
    405469         if (isStrSubset(vNam(varn),"x"))then
    406             x_axis = f->k_x
    407             x_axis = x_axis/normx
    408             if (normx .NE. 1.)then
    409                res@tiXAxisString = "k_x / "+normx
     470            x_axis = new((/dimz,dimx/),double)
     471            do q=0,dimz-1
     472               x_axis(q,:) = f_att->k_x
     473               x_axis = x_axis/norm_x
     474            end do
     475            if (norm_x .NE. 1.)then
     476               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
    410477            else
    411                res@tiXAxisString = "k_x"
     478               if (norm_height .EQ. 1)then
     479                  res@tiXAxisString = "k~B~x~N~ x z [1]"
     480               else
     481                  res@tiXAxisString = "k~B~x~N~ [1/m]"
     482               end if
    412483            end if
     484            dim_r=dimx
    413485         else
    414             x_axis = f->k_y
    415             x_axis = x_axis/normx
    416             if (normx .NE. 1.)then
    417                res@tiXAxisString = "k_y / "+normx
     486            x_axis=new((/dimz,dimy/),double)
     487            do q=0,dimz-1
     488               x_axis(q,:) = f_att->k_y
     489               x_axis = x_axis/norm_x
     490            end do
     491            if (norm_x .NE. 1.)then
     492               res@tiXAxisString = "k~B~x~N~ ["+unit_x+"]"
    418493            else
    419                res@tiXAxisString = "k_y"
     494               if (norm_height .EQ. 1)then
     495                  res@tiXAxisString = "k~B~x~N~ x z [1]"
     496               else
     497                  res@tiXAxisString = "k~B~x~N~ [1/m]"
     498               end if
    420499            end if
     500            dim_r=dimy
    421501         end if
    422502       
    423503         if (sort .EQ. "time")
     504            res@xyLineColors = color
     505            res@pmLegendDisplayMode     = "Always"
     506            res@pmLegendSide            = "Top"
     507            res@pmLegendParallelPosF    = 1.2
     508            res@pmLegendOrthogonalPosF  = -1.0
     509            res@pmLegendWidthF          = 0.12
     510            res@pmLegendHeightF         = 0.04*(end_time_step-start_time_step+1)
    424511            do p=dimz-1,0,1
    425                if (logy .EQ. 1)then 
     512               if (log_y .EQ. 1)then 
    426513                  do q=0,dimt-1
    427                      do r=0,dimsizes(x_axis)-1
    428                         if (temp(q,p,r) .EQ. 0)then
     514                     do r=0,dim_r-1
     515                        if (data(q,p,r) .EQ. 0)then
    429516                           st=p+start_time_step
    430517                           print(" ")
    … …  
    436523                  end do
    437524               end if
    438                res@trXMinF = min(x_axis)
    439                res@trXMaxF = max(x_axis)
    440525               res@gsnLeftString      = vNam(varn)
    441                res@gsnRightString     = "Height = "+z(p)+"m"
    442                if (normy .NE. 1.)then
    443                   res@tiYAxisString      = vNam(varn)+" / "+normy
    444                else
    445                   res@tiYAxisString      = vNam(varn)
    446                end if
    447                res@xyExplicitLegendLabels  = legend_label
    448                plot(n)  = gsn_csm_xy(wks,x_axis,temp(:,p,:),res)
     526               res@gsnRightString     = "Height = "+z(p)+"m"               
     527               res@tiYAxisString      = "["+unit_y+"]"           
     528               res@xyExplicitLegendLabels  = legend_label             
     529               if (norm_height .EQ. 1)then
     530                  data(:,p,:)=data(:,p,:)*doubletofloat(z(p))
     531                  x_axis(p,:) = x_axis(p,:)*z(p)
     532               end if
     533               res@trXMinF = min(x_axis(p,:))
     534               res@trXMaxF = max(x_axis(p,:))
     535               plot(n)  = gsn_csm_xy(wks,x_axis(p,:),data(:,p,:),res)
    449536               n=n+1
    450537            end do
    … …  
    453540               do p=0,dimt-1           
    454541                  do q=0,dimz-1
    455                      do r=0,dimsizes(x_axis)-1
    456                         if (temp(p,q,r) .EQ. 0)then
     542                     do r=0,dim_r-1
     543                        if (data(p,q,r) .EQ. 0)then
    457544                           st=p+start_time_step
    458545                           print(" ")
    … …  
    462549                        end if
    463550                     end do
    464                   end do
    465                   res@trXMinF = min(x_axis)
    466                   res@trXMaxF = max(x_axis)
    467                   res@gsnLeftString      = vNam(varn)
    468                   res@gsnRightString     = "Time = "+legend_label(p)+"h"
    469                   if (normy .NE. 1.)then
    470                      res@tiYAxisString      = vNam(varn)+" / "+normy
    471                   else
    472                      res@tiYAxisString      = vNam(varn)
    473                   end if
    474                   res@xyExplicitLegendLabels  = legend_label_z
    475                   plot(n)  = gsn_csm_xy(wks,x_axis,temp(p,:,:),res)
    476                   n=n+1
     551                     if (norm_height .EQ. 1 .AND. p .EQ. 0)then
     552                        data(p,q,:) = data(p,q,:)*doubletofloat(legend_label_z(q))
     553                        x_axis(q,:) = x_axis(q,:)*doubletofloat(legend_label_z(q))
     554                     end if
     555                     max_y(q)=max(data(p,q,:))
     556                     min_y(q)=min(data(p,q,:))
     557                     min_x(q)=min(x_axis(q,:))
     558                     max_x(q)=max(x_axis(q,:))                                   
     559                  end do
     560                  do q=0,dimz-1
     561                     res@xyLineColor = color(q)
     562                     if (dash .EQ. 1)then
     563                        res@xyDashPattern = q
     564                     end if
     565                     if (q .EQ. 0)then
     566                        res@tiYAxisString      = "["+unit_y+"]"                 
     567                        res@gsnLeftString      = vNam(varn)
     568                        res@gsnRightString     = "Time = "+legend_label(p)+"h"
     569                        res@trXMinF = min(min_x)
     570                        res@trXMaxF = max(max_x)
     571                        res@trYMinF = min(min_y)
     572                        res@trYMaxF = max(max_y)
     573                       
     574                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
     575                       
     576                        lgres = True
     577                        if (dash .EQ. 0)then
     578                           lgres@lgMonoDashIndex = True
     579                        else
     580                           lgres@lgDashIndexes   = ispan(0,dimz-1,1)
     581                        end if
     582                        if (black .EQ. 1)then
     583                           lgres@lgMonoLineColors = True
     584                        else
     585                           lgres@lgLineColors = color
     586                        end if
     587                        lgres@lgTitleString = "Height [m]" 
     588                        lgres@lgLabelFont        = "helvetica"
     589                        lgres@lgLabelFontHeightF     = font_size_legend
     590                        lgres@lgTitleFontHeightF     = font_size     
     591                        lgres@vpWidthF           = 0.12           
     592                        lgres@vpHeightF          = font_size_legend/5*dimz           
     593 
     594                        lbid = gsn_create_legend(wks,dimz,legend_label_z,lgres)       
     595
     596                        amres = True
     597                        amres@amParallelPosF   = 0.75               
     598                        amres@amOrthogonalPosF = 0.15           
     599                        annoid1 = gsn_add_annotation(plot_h(q),lbid,amres)
     600                     else
     601                        plot_h(q)  = gsn_csm_xy(wks,x_axis(q,:),data(p,q,:),res)
     602                        overlay(plot_h(0),plot_h(q))
     603                     end if
     604                  end do             
     605                  plot(n)=plot_h(0)
     606                  n=n+1
    477607               end do
    478608            end if
    479609         end if
    480          delete(temp)
     610         delete(data)
     611         delete(temp)
    481612         delete(x_axis)
     613         delete(min_x)
     614         delete(max_x)
     615         delete(min_y)
     616         delete(max_y)
     617         delete(plot_h)
    482618      end if
    483619   end do
    … …  
    497633   resP                        = True
    498634   resP@txFont                 = "helvetica"
    499    resP@txString               = f@title
     635   resP@txString               = f_att@title
    500636   resP@txFuncCode             = "~"
    501    resP@txFontHeightF          = 0.014
     637   resP@txFontHeightF          = 0.015
    502638
    503639   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
    504640      gsn_panel(wks,plot,(/n,1/),resP)
    505641   else   
    506       do i = 0,n-1, no_lines*no_columns
    507          if( (i+no_lines*no_columns) .gt. (n-1)) then
    508             gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
     642      do i = 0,n-1, no_rows*no_columns
     643         if( (i+no_rows*no_columns) .gt. (n-1)) then
     644            gsn_panel(wks,plot(i:n-1),(/no_rows,no_columns/),resP)
    509645         else
    510            gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
     646           gsn_panel(wks,plot(i:i+no_rows*no_columns-1),(/no_rows,no_columns/),resP)
    511647         end if
    512648      end do
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |