load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ;*************************************************** ; load ncl_preferences.ncl ;*************************************************** function which_script() local script begin script="spectra" return(script) end if (isfilepresent("~/ncl_preferences.ncl")) then loadscript("~/ncl_preferences.ncl") else if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl")) then loadscript( "~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl") else print(" ") print("'ncl_preferences.ncl' does not exist in $home or $home/palm/current_version/trunk/SCRIPTS/NCL/") print(" ") exit end if end if begin ;*************************************************** ; set up default parameter values and strings ;*************************************************** if (file_1 .EQ. "File in") then print(" ") print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'") print(" ") exit else file_in = file_1 end if if (.not. isfilepresent(file_in)) then print(" ") print("1st input file: '"+file_in+"' does not exist") print(" ") exit end if 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 print(" ") print("'format_out = "+format_out+"' is invalid and set to'x11'") print(" ") format_out="x11" end if if (logx .NE. 0 .AND. logx .NE. 1)then print(" ") print("'logx'= "+logx+" is invalid and set to 1") print(" ") logx = 1 end if if (logy .NE. 0 .AND. logy .NE. 1)then print(" ") print("'logy'= "+logy+" is invalid and set to 1") print(" ") logy = 1 end if if (normy .EQ. 0.) then print(" ") print("You cannot normalise the y-axis with 0, 'normy' is set to 1.0") print(" ") normy = 1.0 end if if (normx .EQ. 0.) then print(" ") print("You cannot normalise the x-axis with 0, 'normx' is set to 1.0") print(" ") normx= 1.0 end if if (sort .NE. "height" .AND. sort .NE. "time") then print(" ") print("'sort'= "+sort+" is invalid and set to 'height'") print(" ") sort = "height" end if if (black .NE. 0 .AND. black .NE. 1)then print(" ") print("'black'= "+black+" is invalid and set to 0") print(" ") black = 0 end if if (dash .NE. 0 .AND. dash .NE. 1)then print(" ") print("'dash'= "+dash+" is invalid and set to 0") print(" ") dash = 0 end if ;*************************************************** ; open input file ;*************************************************** f=addfile(file_in,"r") vNam = getfilevarnames(f) print(" ") print("Variables in input file:") print("- "+ vNam) print(" ") dim = dimsizes(vNam) vDim = getfiledimsizes(f) t_all = f->time nt = dimsizes(t_all) delta_t=t_all(nt-1)/nt k_x=f->k_x dimx=dimsizes(k_x) k_y=f->k_y dimy=dimsizes(k_y) dim_level=dimsizes(height_level) do i=0,dim-1 if (vNam(i) .EQ. "zu_sp")then zu=f->zu_sp if (height_level(0) .EQ. -1)then dimz=dimsizes(zu) else if (dim_level .GT. dimsizes(zu))then print(" ") print("'height_level' has more elements than available height levels in input file (= "+dimz+")") print(" ") exit else zuh=new(dim_level,double) do le=0,dim_level-1 zuh(le)=zu(height_level(le)) end do dimz=dim_level end if end if else if (vNam(i) .EQ. "zw_sp")then zw=f->zw_sp if (height_level(0) .EQ. -1)then dimz=dimsizes(zw) else if (dim_level .GT. dimsizes(zw))then print(" ") print("'height_level' has more elements than available height levels in input file (= "+dimz+")") print(" ") exit else zwh=new(dim_level,double) do le=0,dim_level-1 zwh(le)=zw(height_level(le)) end do dimz=dim_level end if end if end if end if end do ;**************************************************** ; start of time step and different types of mistakes that could be done ;**************************************************** if (start_time_step .EQ. -1.d) then start_time_step=t_all(0)/3600 else if (start_time_step .GT. t_all(nt-1)/3600)then print(" ") 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") print(" ") print("Please select another 'start_time_step'") print(" ") exit end if if (start_time_step .LT. t_all(0)/3600)then print(" ") print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h") exit end if end if do i=0,nt-1 if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. start_time_step .LT. (t_all(i)+delta_t/2)/3600)then st=i break end if end do if (.not. isvar("st"))then print(" ") print("'start_time_step' = "+ start_time_step +"h is invalid") print(" ") print("Please select another 'start_time_step'") print(" ") exit end if ;**************************************************** ; end of time step and different types of mistakes that could be done ;**************************************************** if (end_time_step .EQ. -1.d) then end_time_step = t_all(nt-1)/3600 else if (end_time_step .GT. t_all(nt-1)/3600)then print(" ") 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") print(" ") print("Please select another 'end_time_step'") print(" ") exit end if if (end_time_step .LT. start_time_step/3600)then print(" ") print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h") print(" ") print("Please select another 'start_time_step' or 'end_time_step'") print(" ") exit end if end if do i=0,nt-1 if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. end_time_step .LT. (t_all(i)+delta_t/2)/3600)then et=i break end if end do if (.not. isvar("et"))then print(" ") print("'end_time_step' = "+ end_time_step +"h is invalid") print(" ") print("Please select another 'end_time_step'") print(" ") exit end if delete(start_time_step) start_time_step=round(st,3) delete(end_time_step) end_time_step=round(et,3) print(" ") print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+t_all(start_time_step)+" s => index = "+start_time_step) print(" till "+t_all(end_time_step)/3600+" h = "+t_all(end_time_step)+" s => index = "+end_time_step) print(" ") dimt = end_time_step-start_time_step+1 ;*************************************************** ; set up recourses ;*************************************************** res = True res@gsnDraw = False res@gsnFrame = False res@gsnPaperOrientation = "portrait" res@gsnPaperWidth = 8.27 res@gsnPaperHeight = 11.69 res@gsnPaperMargin = 0.79 res@txFont = "helvetica" res@tiMainFont = "helvetica" res@tiXAxisFont = "helvetica" res@tiYAxisFont = "helvetica" res@tmXBLabelFont = "helvetica" res@tmYLLabelFont = "helvetica" res@lgLabelFont = "helvetica" res@tmLabelAutoStride = True res@pmLegendDisplayMode = "Always" res@pmLegendSide = "Top" res@pmLegendParallelPosF = 1.2 res@pmLegendOrthogonalPosF = -1.0 res@pmLegendWidthF = 0.12 if (sort .EQ. "time")then res@pmLegendHeightF = 0.04*(end_time_step-start_time_step+1) else res@pmLegendHeightF = 0.04*dimz end if res@lgLabelFontHeightF = .025 res@lgTitleFontHeightF = .025 res@txFontHeightF = 0.025 res@tiXAxisFontHeightF = 0.025 res@tiYAxisFontHeightF = 0.025 if (logx .EQ. 1) then res@trXLog = True else res@trXLog = False end if if (logy .EQ. 1)then res@trYLog = True else res@trYLog = False end if legend_label=new(dimt,double) legend_label_zu=new(dimz,double) legend_label_zw=new(dimz,double) legend_label_z=new(dimz,double) do p=start_time_step,end_time_step if (t_all(p)/3600 .LT. 1) then legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,2,True) else legend_label(p-start_time_step)=decimalPlaces(t_all(p)/3600,0,True) end if end do if (sort .EQ. "time") plot = new(dim*dimz,graphic) np=dimt res@lgTitleString = "Time [h]" else plot = new(dim*dimt,graphic) np=dimz res@lgTitleString = "Height [m]" do p=0,dimz-1 if (height_level(0) .EQ. -1)then legend_label_zu(p)=round(zu(p),3) legend_label_zw(p)=round(zw(p),3) else legend_label_zu(p)=round(zuh(p),3) legend_label_zw(p)=round(zwh(p),3) end if end do end if if (black .eq. 0 ) then if (np .EQ. 1)then res@xyLineColors = 237 else step=round(235/(np-1),3) res@xyLineColors = ispan(2,237,step) end if end if if ( dash .eq. 0 ) then res@xyMonoDashPattern = True end if wks=gsn_open_wks(format_out,file_out) gsn_define_colormap(wks,"rainbow+white") temp=new((/dimt,dimz,dimx/),float) n=0 do varn =dim-1,0,1 check = True if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "zu_sp") .OR. isStrSubset( vNam(varn), "zw_sp") .OR. isStrSubset( vNam(varn), "k_x") .OR. isStrSubset( vNam(varn), "k_y")) then check = False end if if (var .NE. "all") then check = isStrSubset( var,","+vNam(varn)+"," ) end if if(check) then temp = f->$vNam(varn)$(start_time_step:end_time_step,0:dimz-1,:) a=getvardims(temp) b=dimsizes(a) if (height_level(0) .NE. -1)then do te=0,dimz-1 temp(:,te,:) = f->$vNam(varn)$(start_time_step:end_time_step,height_level(te),:) end do end if temp=temp/(normy*normx) ;SMOOTHING: #temp=smth9(temp/norm, 0.50, -0.25, False)# do i=0,b-1 if (isStrSubset( a(i),"zu_sp" ))then legend_label_z=legend_label_zu if (height_level(0) .NE. -1)then z=zuh else z=zu end if else if (isStrSubset( a(i),"zw_sp" ))then legend_label_z=legend_label_zw if (height_level(0) .NE. -1)then z=zwh else z=zw end if end if end if end do if (isStrSubset(vNam(varn),"x"))then x_axis = f->k_x x_axis = x_axis/normx if (normx .NE. 1.)then res@tiXAxisString = "k_x / "+normx else res@tiXAxisString = "k_x" end if else x_axis = f->k_y x_axis = x_axis/normx if (normx .NE. 1.)then res@tiXAxisString = "k_y / "+normx else res@tiXAxisString = "k_y" end if end if if (sort .EQ. "time") do p=dimz-1,0,1 if (logy .EQ. 1)then do q=0,dimt-1 do r=0,dimsizes(x_axis)-1 if (temp(q,p,r) .EQ. 0)then st=p+start_time_step print(" ") print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and height "+z(p)+" cannot be used") print(" ") res@trYLog = False end if end do end do end if res@trXMinF = min(x_axis) res@trXMaxF = max(x_axis) res@gsnLeftString = vNam(varn) res@gsnRightString = "Height = "+z(p)+"m" if (normy .NE. 1.)then res@tiYAxisString = vNam(varn)+" / "+normy else res@tiYAxisString = vNam(varn) end if res@xyExplicitLegendLabels = legend_label plot(n) = gsn_csm_xy(wks,x_axis,temp(:,p,:),res) n=n+1 end do else if (sort .EQ. "height") do p=0,dimt-1 do q=0,dimz-1 do r=0,dimsizes(x_axis)-1 if (temp(p,q,r) .EQ. 0)then st=p+start_time_step print(" ") print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for y-axis and time "+legend_label(p)+" cannot be used") print(" ") res@trYLog = False end if end do end do res@trXMinF = min(x_axis) res@trXMaxF = max(x_axis) res@gsnLeftString = vNam(varn) res@gsnRightString = "Time = "+legend_label(p)+"h" if (normy .NE. 1.)then res@tiYAxisString = vNam(varn)+" / "+normy else res@tiYAxisString = vNam(varn) end if res@xyExplicitLegendLabels = legend_label_z plot(n) = gsn_csm_xy(wks,x_axis,temp(p,:,:),res) n=n+1 end do end if end if delete(temp) delete(x_axis) end if end do if (n .EQ. 0) then print(" ") print("The variables 'var="+var+"' do not exist on your input file;") print("be sure to have one comma berfore and after each variable") print(" ") exit end if ; *************************************************** ; merge plots onto one page ; *************************************************** resP = True resP@txFont = "helvetica" resP@txString = f@title resP@txFuncCode = "~" resP@txFontHeightF = 0.014 if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then gsn_panel(wks,plot,(/n,1/),resP) else do i = 0,n-1, no_lines*no_columns if( (i+no_lines*no_columns) .gt. (n-1)) then gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP) else gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP) end if end do end if print(" ") print("Output to: " + file_out +"."+ format_out) print(" ") end