Changeset 190 for palm/trunk/SCRIPTS/NCL/spectra.ncl
 Timestamp:
 Aug 16, 2008 12:18:08 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SCRIPTS/NCL/spectra.ncl
r183 r190 2 2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 3 3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 4 4 5 ;*************************************************** 6 ; load ncl_preferences.ncl 7 ;*************************************************** 8 9 if (isfilepresent("~/ncl_preferences.ncl")) then 10 loadscript("~/ncl_preferences.ncl") 11 else 12 if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl")) then 13 loadscript( "~/palm/current_version/trunk/SCRIPTS/NCL/ncl_preferences.ncl") 14 else 15 print(" ") 16 print("'ncl_preferences.ncl' does not exist in $home or $home/palm/current_version/trunk/SCRIPTS/NCL/") 17 print(" ") 18 exit 19 end if 20 end if 21 5 22 begin 6 7 ; *************************************************** 8 ; read parameter_list 9 ; *************************************************** 10 11 if (isfilepresent("~/.ncl_preferences")) then 12 parameter = asciiread("~/.ncl_preferences",129,"string") 13 delete(parameter@_FillValue) 14 else 15 if (isfilepresent("~/palm/current_version/trunk/SCRIPTS/NCL/.ncl_preferences")) then 16 parameter = asciiread("~/palm/current_version/trunk/SCRIPTS/NCL/.ncl_preferences",129,"string") 17 delete(parameter@_FillValue) 18 else 19 print(" ") 20 print("'.ncl_preferences' does not exist in '~/palm/current_version/trunk/SCRIPTS/NCL/'") 21 print(" ") 22 exit 23 end if 24 end if 25 26 if ( .not. isvar("file_1") ) then 27 if (parameter(7) .EQ. "File in") then 28 print(" ") 29 print("Please provide 1st input file 'file_1=' either in prompt or parameter_list") 30 print(" ") 31 exit 32 else 33 file_in = parameter(7) 34 end if 23 24 if (cross_sections .NE. 0 .OR. profiles .NE. 0 .OR. timeseries .NE. 0 .OR. spectra .NE. 1)then 25 print(" ") 26 print("Please specify the used script in 'ncl_preferences.ncl' (Line 710)") 27 print(" ") 28 print("Set 'spectra' equal 1 and the other variables equal 0") 29 print(" ") 30 exit 31 end if 32 33 ;*************************************************** 34 ; set up default parameter values and strings 35 ;*************************************************** 36 37 if (file_1 .EQ. "File in") then 38 print(" ") 39 print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'") 40 print(" ") 41 exit 35 42 else 36 43 file_in = file_1 … … 38 45 if (.not. isfilepresent(file_in)) then 39 46 print(" ") 40 print(" Your1st input file: '"+file_in+"' does not exist")47 print("1st input file: '"+file_in+"' does not exist") 41 48 print(" ") 42 49 exit 43 50 end if 44 51 45 if ( .not. isvar("format_out") ) then 46 format_out = "x11" 47 if (parameter(9) .NE. "x11") then 48 format_out = parameter(9) 49 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 50 print(" ") 51 print("Your 'format_out = "+format_out+"' is invalid and set to'x11'") 52 print(" ") 53 end if 54 end if 55 else 56 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 57 print(" ") 58 print("Your 'format_out = "+format_out+"' is invalid and set to'x11'") 59 print(" ") 60 end if 61 end if 62 63 if ( .not. isvar("file_out") ) then 64 file_out = "test" 65 if (parameter(11) .NE. "test_ts") then 66 file_out = parameter(11) 67 end if 68 end if 69 if ( .not. isvar("no_columns") ) then 70 no_columns = 1 71 if (parameter(17) .NE. "1") then 72 no_columns = stringtointeger(parameter(17)) 73 end if 74 end if 75 if ( .not. isvar("no_lines") ) then 76 no_lines = 2 77 if (parameter(19) .NE. "2") then 78 no_lines = stringtointeger(parameter(19)) 79 end if 52 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 53 print(" ") 54 print("'format_out = "+format_out+"' is invalid and set to'x11'") 55 print(" ") 56 format_out="x11" 80 57 end if 81 82 if (.not. isvar("logy"))then 83 logy = 0 84 if (stringtointeger(parameter(77)) .NE. 0) then 85 logy = stringtointeger(parameter(77)) 86 if (stringtointeger(parameter(77)) .NE. 1) then 87 print(" ") 88 print("Your 'logy'= "+logy+" is invalid and set to 0") 89 print(" ") 90 logy = 0 91 end if 92 end if 93 else 94 if (logy .NE. 0 .AND. logy .NE. 1)then 95 print(" ") 96 print("Your 'logy'= "+logy+" is invalid and set to 0") 97 print(" ") 98 logy = 0 99 end if 100 end if 101 102 if ( .not. isvar("sort") ) then 103 sort = "layer" 104 if (parameter(51) .NE. "layer") then 105 sort = parameter(51) 106 if (sort .NE. "time") then 107 print(" ") 108 print("Your 'sort'= "+sort+" is invalid and set to 'layer'") 109 print(" ") 110 sort = "layer" 111 end if 112 end if 113 else 114 if (sort .NE. "time" .OR. sort .NE. "layer")then 115 print(" ") 116 print("Your 'sort'= "+sort+" is invalid and set to 'layer'") 117 print(" ") 118 sort = "layer" 119 end if 120 end if 121 122 if ( .not. isvar("black") ) then 58 59 if (logx .NE. 0 .AND. logx .NE. 1)then 60 print(" ") 61 print("'logx'= "+logx+" is invalid and set to 1") 62 print(" ") 63 logx = 1 64 end if 65 66 if (logy .NE. 0 .AND. logy .NE. 1)then 67 print(" ") 68 print("'logy'= "+logy+" is invalid and set to 1") 69 print(" ") 70 logy = 1 71 end if 72 73 if (normy .EQ. 0.) then 74 print(" ") 75 print("You cannot normalise the yaxis with 0, 'normy' is set to 1.0") 76 print(" ") 77 normy = 1.0 78 end if 79 if (normx .EQ. 0.) then 80 print(" ") 81 print("You cannot normalise the xaxis with 0, 'normx' is set to 1.0") 82 print(" ") 83 normx= 1.0 84 end if 85 86 if (sort .NE. "height" .AND. sort .NE. "time") then 87 print(" ") 88 print("'sort'= "+sort+" is invalid and set to 'height'") 89 print(" ") 90 sort = "height" 91 end if 92 93 if (black .NE. 0 .AND. black .NE. 1)then 94 print(" ") 95 print("'black'= "+black+" is invalid and set to 0") 96 print(" ") 123 97 black = 0 124 if (parameter(31) .NE. "0") then 125 black = stringtointeger(parameter(31)) 126 if (stringtointeger(parameter(31)) .NE. 1) then 127 print(" ") 128 print("Your 'black'= "+black+" is invalid and set to 0") 129 print(" ") 130 black = 0 131 end if 132 end if 133 else 134 if (black .NE. 0 .AND. black .NE. 1)then 135 print(" ") 136 print("Your 'black'= "+black+" is invalid and set to 0") 137 print(" ") 138 black = 0 139 end if 140 end if 141 if ( .not. isvar("dash") ) then 98 end if 99 100 if (dash .NE. 0 .AND. dash .NE. 1)then 101 print(" ") 102 print("'dash'= "+dash+" is invalid and set to 0") 103 print(" ") 142 104 dash = 0 143 if (parameter(29) .NE. "0") then 144 dash = stringtointeger(parameter(29)) 145 if (stringtointeger(parameter(29)) .NE. 1) then 146 print(" ") 147 print("Your 'dash'= "+dash+" is invalid and set to 0") 148 print(" ") 149 dash = 0 150 end if 151 end if 152 else 153 if (dash .NE. 0 .AND. dash .NE. 1)then 154 print(" ") 155 print("Your 'dash'= "+dash+" is invalid and set to 0") 156 print(" ") 157 dash = 0 158 end if 159 end if 160 161 if ( .not. isvar("norm") ) then 162 norm = 1.0 163 if (parameter(79) .NE. "1") then 164 norm = stringtofloat(parameter(79)) 165 if (stringtofloat(parameter(79)) .EQ. 0) then 166 print(" ") 167 print("You cannot normalise with 0, 'norm' is set to 1") 168 print(" ") 169 norm = 1.0 170 end if 171 end if 172 else 173 if (norm .EQ. 0) then 174 print(" ") 175 print("You cannot normalise with 0, 'norm' is set to 1") 176 print(" ") 177 norm = 1.0 178 end if 179 end if 105 end if 106 107 ;*************************************************** 108 ; open input file 109 ;*************************************************** 180 110 181 111 f=addfile(file_in,"r") … … 183 113 vNam = getfilevarnames(f) 184 114 print(" ") 185 print("Variable in input file: " + vNam) 115 print("Variables in input file:") 116 print(" "+ vNam) 186 117 print(" ") 187 118 dim = dimsizes(vNam) … … 192 123 delta_t=t_all(nt1)/nt 193 124 125 k_x=f>k_x 126 dimx=dimsizes(k_x) 127 k_y=f>k_y 128 dimy=dimsizes(k_y) 129 130 131 dim_level=dimsizes(height_level) 132 194 133 do i=0,dim1 195 134 if (vNam(i) .EQ. "zu_sp")then 196 zu=f>zu_sp 197 z=zu 198 dimz=dimsizes(zu) 135 zu=f>zu_sp 136 if (height_level(0) .EQ. 1)then 137 dimz=dimsizes(zu) 138 else 139 if (dim_level .GT. dimsizes(zu))then 140 print(" ") 141 print("'height_level' has more elements than available height levels in input file (= "+dimz+")") 142 print(" ") 143 exit 144 else 145 zuh=new(dim_level,double) 146 do le=0,dim_level1 147 zuh(le)=zu(height_level(le)) 148 end do 149 dimz=dim_level 150 end if 151 end if 199 152 else 200 153 if (vNam(i) .EQ. "zw_sp")then 201 zw=f>zw_sp 202 z=zw 203 dimz=dimsizes(zw) 154 zw=f>zw_sp 155 if (height_level(0) .EQ. 1)then 156 dimz=dimsizes(zw) 157 else 158 if (dim_level .GT. dimsizes(zw))then 159 print(" ") 160 print("'height_level' has more elements than available height levels in input file (= "+dimz+")") 161 print(" ") 162 exit 163 else 164 zwh=new(dim_level,double) 165 do le=0,dim_level1 166 zwh(le)=zw(height_level(le)) 167 end do 168 dimz=dim_level 169 end if 170 end if 204 171 end if 205 172 end if 206 173 end do 207 208 ; 174 175 ;**************************************************** 209 176 ; start of time step and different types of mistakes that could be done 210 ; 211 212 if ( .not. isvar("start_time_step")) then177 ;**************************************************** 178 179 if (start_time_step .EQ. 1.d) then 213 180 start_time_step=t_all(0)/3600 214 if (parameter(13) .NE. "t(0)") then 215 if (stringtodouble(parameter(13)) .GT. t_all(nt1)/3600)then 216 print(" ") 217 print("'start_time_step' = "+ parameter(13) +"h is greater than last time step = " + t_all(nt1)+"s = "+t_all(nt1)/3600+"h") 218 print(" ") 219 print("Please select another 'start_time_step'") 220 print(" ") 221 exit 222 end if 223 if (stringtofloat(parameter(13)) .LT. t_all(0)/3600)then 224 print(" ") 225 print("'start_time_step' = "+ parameter(13) +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h") 226 exit 227 end if 228 start_time_step=stringtodouble(parameter(13)) 229 end if 230 else 181 else 231 182 if (start_time_step .GT. t_all(nt1)/3600)then 232 183 print(" ") … … 243 194 end if 244 195 end if 245 start_time_step = start_time_step*3600246 196 247 197 do i=0,nt1 248 if (start_time_step .GE. t_all(i)delta_t/2 .AND. start_time_step .LT. t_all(i)+delta_t/2)then198 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 249 199 st=i 250 200 break … … 252 202 end do 253 203 254 ; **************************************************** 204 if (.not. isvar("st"))then 205 print(" ") 206 print("'start_time_step' = "+ start_time_step +"h is invalid") 207 print(" ") 208 print("Please select another 'start_time_step'") 209 print(" ") 210 exit 211 end if 212 213 ;**************************************************** 255 214 ; end of time step and different types of mistakes that could be done 256 ; 257 258 if ( .not. isvar("end_time_step")) then215 ;**************************************************** 216 217 if (end_time_step .EQ. 1.d) then 259 218 end_time_step = t_all(nt1)/3600 260 if (parameter(15) .NE. "t(end)") then261 if (stringtodouble(parameter(15)) .GT. t_all(nt1)/3600)then262 print(" ")263 print("'end_time_step' = "+ parameter(15) +"h is greater than last time step = " + t_all(nt1)+"s = "+t_all(nt1)/3600+"h")264 print(" ")265 print("Please select another 'end_time_step'")266 print(" ")267 exit268 end if269 if (stringtodouble(parameter(15)) .LT. start_time_step/3600)then270 print(" ")271 print("'end_time_step' = "+ parameter(15) +"h is lower than 'start_time_step' = "+start_time_step/3600+"h")272 print(" ")273 print("Please select another 'start_time_step' or 'end_time_step'")274 print(" ")275 exit276 end if277 end_time_step = stringtodouble(parameter(15))278 end if279 219 else 280 220 if (end_time_step .GT. t_all(nt1)/3600)then … … 295 235 end if 296 236 end if 297 end_time_step = end_time_step*3600298 237 299 238 do i=0,nt1 300 if (end_time_step .GE. t_all(i)delta_t/2 .AND. end_time_step .LT. t_all(i)+delta_t/2)then239 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 301 240 et=i 302 241 break 303 242 end if 304 243 end do 244 245 if (.not. isvar("et"))then 246 print(" ") 247 print("'end_time_step' = "+ end_time_step +"h is invalid") 248 print(" ") 249 print("Please select another 'end_time_step'") 250 print(" ") 251 exit 252 end if 305 253 306 254 delete(start_time_step) … … 316 264 dimt = end_time_stepstart_time_step+1 317 265 318 ; 266 ;*************************************************** 319 267 ; set up recourses 320 ; 268 ;*************************************************** 321 269 322 270 res = True … … 337 285 res@pmLegendDisplayMode = "Always" 338 286 res@pmLegendSide = "Top" 339 res@pmLegendParallelPosF = 1. 15287 res@pmLegendParallelPosF = 1.2 340 288 res@pmLegendOrthogonalPosF = 1.0 341 289 res@pmLegendWidthF = 0.12 … … 345 293 res@pmLegendHeightF = 0.04*dimz 346 294 end if 347 res@lgLabelFontHeightF = .02 348 res@lgTitleFontHeightF = .02 349 res@txFontHeightF = 0.02 350 res@tiXAxisFontHeightF = 0.02 351 res@tiYAxisFontHeightF = 0.02 352 353 if (log y.EQ. 1) then295 res@lgLabelFontHeightF = .025 296 res@lgTitleFontHeightF = .025 297 res@txFontHeightF = 0.025 298 res@tiXAxisFontHeightF = 0.025 299 res@tiYAxisFontHeightF = 0.025 300 301 if (logx .EQ. 1) then 354 302 res@trXLog = True 303 else 304 res@trXLog = False 305 end if 306 if (logy .EQ. 1)then 355 307 res@trYLog = True 356 else 357 res@trXLog = False 308 else 358 309 res@trYLog = False 359 310 end if … … 379 330 res@lgTitleString = "Height [m]" 380 331 do p=0,dimz1 381 legend_label_zu(p)=round(zu(p),3) 382 legend_label_zw(p)=round(zw(p),3) 332 if (height_level(0) .EQ. 1)then 333 legend_label_zu(p)=round(zu(p),3) 334 legend_label_zw(p)=round(zw(p),3) 335 else 336 legend_label_zu(p)=round(zuh(p),3) 337 legend_label_zw(p)=round(zwh(p),3) 338 end if 383 339 end do 384 340 end if 385 386 if ( black .eq. 0 ) then 387 res@xyLineColors = ispan(2,237,235/np) 341 342 step=round(235/(np1),3) 343 if (black .eq. 0 ) then 344 res@xyLineColors = ispan(2,237,step) 388 345 end if 389 346 if ( dash .eq. 0 ) then 390 res@xyMonoDashPattern 347 res@xyMonoDashPattern = True 391 348 end if 392 349 393 350 wks=gsn_open_wks(format_out,file_out) 394 351 gsn_define_colormap(wks,"rainbow+white") 352 353 temp=new((/dimt,dimz,dimx/),float) 395 354 396 355 n=0 … … 403 362 end if 404 363 405 if (.not. isvar("var")) then 406 if (parameter(21) .NE. "variables") then 407 var=parameter(21) 408 check = isStrSubset( var,","+vNam(varn)+"," ) 409 end if 410 else 364 if (var .NE. "all") then 411 365 check = isStrSubset( var,","+vNam(varn)+"," ) 412 366 end if 413 367 414 if(check) then 415 368 if(check) then 369 416 370 temp = f>$vNam(varn)$(start_time_step:end_time_step,0:dimz1,:) 417 418 temp=temp/norm ;SMOOTHING: #temp=smth9(temp/norm, 0.50, 0.25, False)# 419 371 420 372 a=getvardims(temp) 421 b=dimsizes(a) 373 b=dimsizes(a) 374 375 if (height_level(0) .NE. 1)then 376 do te=0,dimz1 377 temp(:,te,:) = f>$vNam(varn)$(start_time_step:end_time_step,height_level(te),:) 378 end do 379 end if 380 381 temp=temp/(normy*normx) ;SMOOTHING: #temp=smth9(temp/norm, 0.50, 0.25, False)# 382 422 383 do i=0,b1 423 384 if (isStrSubset( a(i),"zu_sp" ))then 424 385 legend_label_z=legend_label_zu 386 if (height_level(0) .NE. 1)then 387 z=zuh 388 else 389 z=zu 390 end if 425 391 else 426 392 if (isStrSubset( a(i),"zw_sp" ))then 427 legend_label_z=legend_label_zw 393 legend_label_z=legend_label_zw 394 if (height_level(0) .NE. 1)then 395 z=zwh 396 else 397 z=zw 398 end if 428 399 end if 429 400 end if 430 401 end do 402 431 403 if (isStrSubset(vNam(varn),"x"))then 432 404 x_axis = f>k_x 433 res@tiXAxisString = "k_x" 405 x_axis = x_axis/normx 406 if (normx .NE. 1.)then 407 res@tiXAxisString = "k_x / "+normx 408 else 409 res@tiXAxisString = "k_x" 410 end if 434 411 else 435 412 x_axis = f>k_y 436 res@tiXAxisString = "k_y" 413 x_axis = x_axis/normx 414 if (normx .NE. 1.)then 415 res@tiXAxisString = "k_y / "+normx 416 else 417 res@tiXAxisString = "k_y" 418 end if 437 419 end if 438 420 439 421 if (sort .EQ. "time") 440 do p=dimz1,0,1 441 do q=0,dimt1 442 do r=0,dimsizes(x_axis)1 443 if (temp(q,p,r) .EQ. 0)then 444 st=p+start_time_step 445 print(" ") 446 print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for yaxis cannot be used") 447 print(" ") 448 res@trYLog = False 449 end if 422 do p=dimz1,0,1 423 if (logy .EQ. 1)then 424 do q=0,dimt1 425 do r=0,dimsizes(x_axis)1 426 if (temp(q,p,r) .EQ. 0)then 427 st=p+start_time_step 428 print(" ") 429 print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for yaxis and height "+z(p)+" cannot be used") 430 print(" ") 431 res@trYLog = False 432 end if 433 end do 450 434 end do 451 end do 435 end if 452 436 res@trXMinF = min(x_axis) 453 437 res@trXMaxF = max(x_axis) 454 438 res@gsnLeftString = vNam(varn) 455 439 res@gsnRightString = "Height = "+z(p)+"m" 456 if (norm .NE. 1)then457 res@tiYAxisString = vNam(varn)+" / "+norm 440 if (normy .NE. 1.)then 441 res@tiYAxisString = vNam(varn)+" / "+normy 458 442 else 459 443 res@tiYAxisString = vNam(varn) … … 464 448 end do 465 449 else 466 if (sort .EQ. " layer")467 do p= dimt1,0,1450 if (sort .EQ. "height") 451 do p=0,dimt1 468 452 do q=0,dimz1 469 453 do r=0,dimsizes(x_axis)1 … … 471 455 st=p+start_time_step 472 456 print(" ") 473 print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for yaxis cannot be used")457 print("'"+vNam(varn)+"("+st+","+q+","+r+")' is equal 0; Logarithmic scale for yaxis and time "+legend_label(p)+" cannot be used") 474 458 print(" ") 475 459 res@trYLog = False … … 481 465 res@gsnLeftString = vNam(varn) 482 466 res@gsnRightString = "Time = "+legend_label(p)+"h" 483 if (norm .NE. 1)then484 res@tiYAxisString = vNam(varn)+" / "+norm 467 if (normy .NE. 1.)then 468 res@tiYAxisString = vNam(varn)+" / "+normy 485 469 else 486 470 res@tiYAxisString = vNam(varn) … … 490 474 n=n+1 491 475 end do 492 else493 print(" ")494 print("Please choose 'sort' either equal 'time' or 'height'")495 print(" ")496 exit497 476 end if 498 477 end if 499 478 delete(temp) 500 479 delete(x_axis) … … 504 483 if (n .EQ. 0) then 505 484 print(" ") 506 print("The variables 'var=Â°"+var+"Â°' do not exist on your input file") 485 print("The variables 'var="+var+"' do not exist on your input file;") 486 print("be sure to have one comma berfore and after each variable") 507 487 print(" ") 508 488 exit
Note: See TracChangeset
for help on using the changeset viewer.