;load "/usr/lib/ncarg/nclscripts/csm/gsn_code.ncl" ;load "/usr/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ;load "/usr/lib/ncarg/nclscripts/csm/contributed.ncl" 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" begin ;=================================================; ; open file and read in data ;=================================================; time_plot = (/ 11.0, 13.0, 15.0 /) * 3600.0 height_plot = (/ 0.1, 0.5, 0.75 /) ;115.0, 1150.0, 1800.0 /) ; 2070 height_plot = (/ 50.0, 200.0, 400.0 /) ; in meters upper_level = 1000.0 tinterval = 1800.0 j_start = 0 ; 800 running_av = 10 writeascii = True ;False ;True readascii = True if ( writeascii .eq. True ) then readascii = False end if ; load two files, the original one and one with the correct latitudes f3d = addfile( "/scratch/usr/nikmsueh/test_ches/OUTPUT/test_ches_av_3d.000.nc", "r") fts = addfile( "/scratch/usr/nikmsueh/test_ches/OUTPUT/test_ches_ts.000.nc", "r") time = f3d->time z = f3d->zu_3d y = f3d->y x = f3d->x dim_x = dimsizes(x) dim_y = dimsizes(y) z_ind_up = closest_val( upper_level, z ) wks = gsn_open_wks("eps","adjustment_tke" ) gsn_define_colormap(wks,"BlGrYeOrReVi200") plot = new(3,graphic) res = True res@gsnDraw = False ;False res@gsnFrame = False ;False res@vpWidthF = 0.7 res@vpHeightF = 0.2 res@tiMainOn = False res@tiYAxisString = "E (m~S~2~N~s~S~-2~N~)" ;/ w~B~*~N~~S~2~N~" ;~F33~s~F7~~B~u~B~H~N~~N~ , ~F33~s~F7~~B~w~N~ (m~S~2~N~s~S~-2~N~)" res@tiXAxisString = "~F7~d (m)" res@tiXAxisOn = False res@tiXAxisFontHeightF = 0.02 res@tiYAxisFontHeightF = 0.02 res@xyLineThicknessF = 2 res@xyDashPatterns = (/0, 4, 16 /) res@trXMinF = 0.0 res@trXMaxF = 50000.0 ; ;-Define top x-axis resources do t = 0, dimsizes(time_plot)-1 t_ind = closest_val( time_plot(t), time ) z_ind = new( dimsizes( height_plot), "integer" ) ; first, calculate the mean wind direction averaged over ABL depth time_ts = fts->time vars_ts = getfilevarnames(fts) tsts = closest_val( time_plot(t) - tinterval, time_ts ) tets = closest_val( time_plot(t), time_ts ) zi = fts->zi_theta ;zi_stg ; zi_theta ;$vars_ts(37)$ ; ;-Calculate mean boundary-layer depth and index of mean boundary-layer depth in profile data zi_av = 400.0 ;avg( zi(tsts:tets) ) zi_ind = closest_val( zi_av, z ) do k = 0, dimsizes(z_ind)-1 z_ind(k) = closest_val( height_plot(k), z ) end do print("plots at :" + z(z_ind)) ; Calculate variances ww_av = f3d->ww(t_ind,:z_ind_up,:,:) w_av = f3d->w(t_ind,:z_ind_up,:,:) uu_av = f3d->uu(t_ind,:z_ind_up,:,:) u_av = f3d->u(t_ind,:z_ind_up,:,:) vv_av = f3d->vv(t_ind,:z_ind_up,:,:) v_av = f3d->v(t_ind,:z_ind_up,:,:) ww = ww_av - w_av * w_av uu = uu_av - u_av * u_av vv = vv_av - v_av * v_av delete( u_av ) delete( v_av ) delete( w_av ) delete( uu_av ) delete( vv_av ) delete( ww_av ) tke = 0.5 * ( ww + uu + vv ) ; --------------------------------------- ; determine adjustment u_av = f3d->u(t_ind,:,:,:) v_av = f3d->v(t_ind,:,:,:) avg_u = avg(u_av(1:zi_ind,:,:)) avg_v = avg(v_av(1:zi_ind,:,:)) u_mean = sqrt( avg_u^2 + avg_v^2 ) print( "zi" + " " + zi_av ) ; ;- Calculate wind direction in degrees and absolute wind speed print( avg_u + " " + avg_v ) delete ( u_av ) delete ( v_av ) r2d = 45.0/atan(1.0) ; conversion factor dir = atan2(avg_u,avg_v) * r2d + 180.0 ; compute wind direction vel_abs = sqrt( avg_u*avg_u + avg_v*avg_v ) if ( avg_u .ge. 0.0 ) then x_bound = 0.0 x_bound_out = x(dim_x-1) else x_bound = x(dim_x-1) x_bound_out = 0.0 end if if ( avg_v .ge. 0.0 ) then y_bound = 0.0 y_bound_out = y(dim_y-1) else y_bound = y(dim_y-1) y_bound_out = 0.0 end if distances_precalculated = new( (/dim_y, dim_x/), "float" ) distances_precalculated = 0.0 u_vec = -avg_u v_vec = -avg_v ; compute distance to inflow boundary for each i,j grid point - ; therefore compute interception point between line equation (setup by wind direction) and the inflow boundary do i = 0, dim_x-1 do j = j_start, dim_y-1 line2_beta = y(j) + ( x_bound - x(i) ) / u_vec line1_alpha = x(i) - y(j) * u_vec / v_vec intersection_x2 = x_bound intersection_y2 = line2_beta intersection_x1 = line1_alpha intersection_y1 = 0.0 distances_precalculated(j,i) = min( (/ doubletofloat( sqrt( ( x(i) - intersection_x1 )^2 + ( y(j) - intersection_y1 )^2 ) ), \ doubletofloat( sqrt( ( x(i) - intersection_x2 )^2 + ( y(j) - intersection_y2 )^2 ) ) /) ) end do end do print( " distances calculated " ) ; 2D array into 1D array distances = ndtooned( distances_precalculated ) dim_d = dimsizes(distances) tke_1d = new((/z_ind_up+1,dim_d/),"float") tke_1d = 0.0 do k = 0, z_ind_up tke_1d(k,:) = ndtooned( tke(k,:,:) ) end do ; ;- Get the unique values, i.e. filter out equal distances. dim_d = count_unique_values( distances(:) ) d2 = get_unique_values( distances(:) ) max_d2 = max(d2) min_d2 = 0.0 num_bins = floattointeger(max_d2 / 100.0) bin = new(num_bins, "float") ; compute bin sizes (100 bins seem to give a smooth graph) do b = 0, num_bins-1 bin(b) = b * 100.0 end do ; print( d2 ) ; print( bin ) ; Array for tke values along distance tke_adj = new((/dimsizes(z_ind), num_bins/),float) tke_adj = 0.0 ; average over similar distances to the inflow boundary. Distances are binned. do b = 0, num_bins-2 inds = ind(distances .ge. bin(b) .and. distances .lt. bin(b+1)) if ( dimsizes(inds) .gt. 1 ) then tke_adj(:,b) = dim_avg( tke_1d(z_ind,inds) ) end if delete( inds ) end do print( "plot now" ) ; ;- Plot the variances if ( t .eq. 2 ) then res@tiXAxisOn = True res@tiXAxisString = "~F7~d (m)" end if res@trYMinF = 0.0 res@trYMaxF = max( runave(tke_adj(:,:),running_av,0) ) * 1.1 res@xyLineColor = 1 plot(t) = gsn_csm_xy( wks,bin(:),runave(tke_adj(:,:),running_av,0),res ) delete( tke_adj ) delete( bin ) delete( tke_1d ) delete( d2 ) delete( tke ) delete( distances_precalculated ) delete( x_bound ) delete( x_bound_out ) delete( y_bound ) delete( y_bound_out ) end do ; Add legends lab = (/ " 0.1 z~B~i~N~", " 0.5 z~B~i~N~", " 0.75 z~B~i~N~" /) lgres = True lgres@lgLineThicknessF = 2 lgres@lgDashIndexes = res@xyDashPatterns lgres@lgLineColors = (/ 1, 1, 1 /) lgres@lgLabelFontHeightF = .1 lgres@lgPerimOn = False lgres@vpWidthF = 0.16 lgres@vpHeightF = 0.1 legend = gsn_create_legend(wks,dimsizes(lab),lab,lgres) amres = True amres@amJust = "TopRight" ;"TopLeft" ;-- Define the point of the box ;-- for determining its location. amres@amParallelPosF = 0.4 ;-- Move legend to right amres@amOrthogonalPosF = -0.4 ;4 ;-- Move legend down. annoid = gsn_add_annotation(plot(0),legend, amres) lab2 = (/ " H~B~0~N~", " LE~B~0~N~"/) lgres2 = True lgres2@lgLineThicknessF = 2 lgres2@lgDashIndexes = (/ 0, 0 /) lgres2@lgLineColors = (/ 1, 150 /) lgres2@lgLabelFontHeightF = .1 lgres2@lgPerimOn = False lgres2@vpWidthF = 0.11 lgres2@vpHeightF = 0.07 legend2 = gsn_create_legend(wks,dimsizes(lab2),lab2,lgres2) amres2 = True amres2@amJust = "TopRight" ;"TopLeft" ;-- Define the point of the box ;-- for determining its location. amres2@amParallelPosF = 0.48 ;-- Move legend to right amres2@amOrthogonalPosF = -0.43 ;4 ;-- Move legend down. ; annoid2 = gsn_add_annotation(plot(1),legend2, amres) ; Setup the panel plot ;--------------------------------------------------------------------------------- resP = True ; modify the panel plot resP@gsnFrame = False resP@gsnPaperOrientation = "portrait" resP@gsnPaperWidth = 8.27 resP@gsnPaperHeight = 11.69 resP@gsnPaperMargin = 0.79 ;resP@cnLinesOn = False resP@gsnPanelFigureStrings= (/"a","b", "c"/) resP@amJust = "Bottom" resP@amOrthogonalPosF = -0.31 resP@amParallelPosF = -0.0 ;41 resP@gsnPanelFigureStringsPerimOn = False resP@gsnPanelLabelBar = False ; move lable bar up ; resP@pmLabelBarOrthogonalPosF = 0.25 ; resP@pmLabelBarParallelPosF = 0.0 resP@gsnMaximize = False ;True resP@gsnPanelBottom = 0.1 resP@gsnBoxMargin = 0.0 gsn_panel(wks,plot(0:2),(/3,1/),resP) end