NCL geostrophic: Difference between revisions
(Created page with '<pre>;This is a NCL file with custom functions for the GFI dynamic meteorology group load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "WRFUserARW_new.ncl" ;a…') |
(No difference)
|
Latest revision as of 15:20, 11 January 2013
;This is a NCL file with custom functions for the GFI dynamic meteorology group load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "WRFUserARW_new.ncl" ;a function for calculating the geostrophic and ageostrophic wind along and across a given cross section undef("wrf_geostrophic") function wrf_geostrophic(nc_file:file,it:numeric,lat_diff_f:numeric,lon_diff_f:numeric, \ WE:logical,NS:logical) begin a=nc_file P_total = wrf_user_getvar(a,"pres",it) P_base = wrf_user_getvar(a,"PB",it) G_pert = wrf_user_getvar(a,"PH",it) G_base = wrf_user_getvar(a,"PHB",it) q_vapor = wrf_user_getvar(a,"QVAPOR",it) q_cloud = wrf_user_getvar(a,"QCLOUD",it) q_rain = wrf_user_getvar(a,"QRAIN",it) mu = wrf_user_getvar(a,"MU",it) mub = wrf_user_getvar(a,"MUB",it) eta_full = wrf_user_getvar(a,"ZNW",it) eta_half = wrf_user_getvar(a,"ZNU",it) f = wrf_user_getvar(a,"F",it) ;P_total = P_pert+P_base G_total = G_pert+G_base mu_total = mu+mub delta_y = 3000. delta_x = 3000. dimG = dimsizes(G_pert) dimP = dimsizes(P_total) ;;;;;;Calculating all finite differences;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; delta_G_eta = G_total(1:dimG(0)-1,:,:)-G_total(0:dimG(0)-2,:,:) delta_G_x = G_total(:,:,1:dimG(2)-1)-G_total(:,:,0:dimG(2)-2) delta_G_y = G_total(:,1:dimG(1)-1,:)-G_total(:,0:dimG(1)-2,:) delta_P_eta = P_total(1:dimP(0)-1,:,:)-P_total(0:dimP(0)-2,:,:) delta_P_x = P_total(:,:,1:dimP(2)-1)-P_total(:,:,0:dimP(2)-2) delta_P_y = P_total(:,1:dimP(1)-1,:)-P_total(:,0:dimP(1)-2,:) delta_eta_full = eta_full(1:dimG(0)-1)-eta_full(0:dimG(0)-2) delta_eta_half = eta_half(1:dimP(0)-1)-eta_half(0:dimP(0)-2) delta_eta_full_conform_P = conform_dims((/dimP(0),dimP(1),dimP(2)/),delta_eta_full,0) delta_eta_half_conform_G = conform_dims((/dimG(0)-2,dimG(1),dimG(2)/),delta_eta_half,0) delta_G_eta = delta_G_eta/delta_eta_full_conform_P delta_P_eta = delta_P_eta/delta_eta_half_conform_G ;;;;;;;;Calculating alpha and alpha_D;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; mu_total_conform_P = conform_dims((/dimP(0),dimP(1),dimP(2)/),mu_total,(/1,2/)) alpha_D = -delta_G_eta/mu_total_conform_P alpha = alpha_D*(1+q_vapor+q_cloud+q_rain) ;;;;;;;;getting right dimensions and interpolating to vg;;;;;;;;;;;;; mu_total_interpol_x = (mu_total(:,0:dimP(2)-2)+mu_total(:,1:dimP(2)-1))/2 mu_total_interpol_y = (mu_total(0:dimP(1)-2,:)+mu_total(1:dimP(1)-1,:))/2 alpha_interpol_x = (alpha(:,:,0:dimP(2)-2)+alpha(:,:,1:dimP(2)-1))/2 alpha_interpol_y = (alpha(:,0:dimP(1)-2,:)+alpha(:,1:dimP(1)-1,:))/2 alpha_D_interpol_x = (alpha_D(:,:,0:dimP(2)-2)+alpha_D(:,:,1:dimP(2)-1))/2 alpha_D_interpol_y = (alpha_D(:,0:dimP(1)-2,:)+alpha_D(:,1:dimP(1)-1,:))/2 f_interpol_x = (f(:,0:dimP(2)-2)+f(:,1:dimP(2)-1))/2 f_interpol_y = (f(0:dimP(1)-2,:)+f(1:dimP(1)-1,:))/2 delta_P_eta_interpol_vg = (delta_P_eta(0:dimG(0)-4,:,0:dimG(2)-2)+delta_P_eta(1:dimG(0)-3,:,0:dimG(2)-2)+delta_P_eta(1:dimG(0)-3,:,1:dimG(2)-1)+delta_P_eta(0:dimG(0)-4,:,1:dimG(2)-1))/4 delta_P_eta_interpol_ug = (delta_P_eta(0:dimG(0)-4,0:dimG(1)-2,:)+delta_P_eta(1:dimG(0)-3,0:dimG(1)-2,:)+delta_P_eta(1:dimG(0)-3,1:dimG(1)-1,:)+delta_P_eta(0:dimG(0)-4,1:dimG(1)-1,:))/4 delta_G_x_interpol_vg = (delta_G_x(1:dimG(0)-3,:,0:dimG(2)-2)+delta_G_x(2:dimG(0)-2,:,0:dimG(2)-2))/2 delta_G_y_interpol_ug = (delta_G_y(1:dimG(0)-3,0:dimG(1)-2,:)+delta_G_y(2:dimG(0)-2,0:dimG(1)-2,:))/2 mu_total_interpol_x_conform_P = conform_dims((/dimP(0),dimP(1),dimP(2)-1/),mu_total_interpol_x,(/1,2/)) mu_total_interpol_y_conform_P = conform_dims((/dimP(0),dimP(1)-1,dimP(2)/),mu_total_interpol_y,(/1,2/)) f_interpol_x_conform_P = conform_dims((/dimP(0),dimP(1),dimP(2)-1/),f_interpol_x,(/1,2/)) f_interpol_y_conform_P = conform_dims((/dimP(0),dimP(1)-1,dimP(2)/),f_interpol_y,(/1,2/)) delta_eta_half_conform_Gx = conform_dims((/dimG(0)-2,dimG(1),dimG(2)-1/),delta_eta_half,0) delta_eta_half_conform_Gy = conform_dims((/dimG(0)-2,dimG(1)-1,dimG(2)/),delta_eta_half,0) ;;;;;;split it up in two terms and adding afterwards;;;;;;;;;;;;;;;;; vg_term1 = mu_total_interpol_x_conform_P(1:dimP(0)-2,:,:)*alpha_interpol_x(1:dimG(0)-3,:,:)/f_interpol_x_conform_P(1:dimP(0)-2,:,:)*delta_P_x(1:dimP(0)-2,:,:)/delta_x ug_term1 = -mu_total_interpol_y_conform_P(1:dimP(0)-2,:,:)*alpha_interpol_y(1:dimG(0)-3,:,:)/f_interpol_y_conform_P(1:dimP(0)-2,:,:)*delta_P_y(1:dimP(0)-2,:,:)/delta_y vg_term2 = delta_P_eta_interpol_vg*delta_G_x_interpol_vg/delta_x*alpha_interpol_x(1:dimG(0)-3,:,:)/(alpha_D_interpol_x(1:dimG(0)-3,:,:)*f_interpol_x_conform_P(1:dimP(0)-2,:,:)) ug_term2 = -delta_P_eta_interpol_ug*delta_G_y_interpol_ug/delta_y*alpha_interpol_y(1:dimG(0)-3,:,:)/(alpha_D_interpol_y(1:dimG(0)-3,:,:)*f_interpol_y_conform_P(1:dimP(0)-2,:,:)) vg=(vg_term1+vg_term2)/mu_total_interpol_x_conform_P(1:dimP(0)-2,:,:) ug=(ug_term1+ug_term2)/mu_total_interpol_y_conform_P(1:dimP(0)-2,:,:) ;interpolate the geostrophic winds to mass points vga=(vg(:,:,0:dimP(2)-3)+vg(:,:,1:dimP(2)-2))/2 uga=(ug(:,0:dimP(1)-3,:)+ug(:,1:dimP(1)-2,:))/2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Calculate the angle between th cross section and the regular grid if(WE .and. lat_diff_f.ne.0.) then gamma = atan(lat_diff_f/lon_diff_f) else if(WE .and. lat_diff_f.eq.0.) then gamma=0. end if end if if(NS .and. lon_diff_f.lt.0.)then lat_diff_f=-lat_diff_f gamma=3.14/2.0-atan(lat_diff_f/lon_diff_f) else if(NS .and. lon_diff_f.gt.0.)then gamma=atan(lat_diff_f/lon_diff_f)+3.14/2. else if(NS .and. lon_diff_f.eq.0)then gamma=0. end if end if end if ;;;;Make use of vector algebra and use a rotation matrix :-);;;;;;;;; print(gamma*180/3.14) wwg_along = cos(gamma)*uga(:,:,1:dimP(2)-2)+sin(gamma)*vga(:,1:dimP(1)-2,:) wwg_across = -sin(gamma)*uga(:,:,1:dimP(2)-2)+cos(gamma)*vga(:,1:dimP(1)-2,:) wwg_along@description = "Along geostrophic Wind" wwg_across@description = "Across geostrophic Wind" u = wrf_user_getvar(a,"U",it) ; 3D U at mass points v = wrf_user_getvar(a,"V",it) ; 3D V at mass points theta = wrf_user_getvar(a,"theta",it) ;;;;interpolate it to masspoints myself, since i don't know what wrf_user_getvar is really doing;;;;;;; ua=(u(1:dimP(0)-2,1:dimP(1)-2,0:dimP(2)-1)+u(1:dimP(0)-2,1:dimP(1)-2,1:dimP(2)))/2 va=(v(1:dimP(0)-2,0:dimP(1)-1,1:dimP(2)-2)+v(1:dimP(0)-2,1:dimP(1),1:dimP(2)-2))/2 ww_along = cos(gamma)*ua(:,:,1:dimP(2)-2)+sin(gamma)*va(:,1:dimP(1)-2,:) ww_across = -sin(gamma)*ua(:,:,1:dimP(2)-2)+cos(gamma)*va(:,1:dimP(1)-2,:) ;calculating the ageostrophic components along and across the cross section wwageo_along = ww_along-wwg_along wwageo_across = ww_across-wwg_across wwageo_along@description = "Along ageostrophic wind" wwageo_across@description= "Across ageostrophic wind" all_winds=new((/dimP(0)-2,dimP(1)-2,dimP(2)-2,6/),typeof(ww_along)) all_winds(:,:,:,0)=ww_along all_winds(:,:,:,0)=ww_along all_winds(:,:,:,1)=ww_across all_winds(:,:,:,2)=wwg_along all_winds(:,:,:,3)=wwg_across all_winds(:,:,:,4)=wwageo_along all_winds(:,:,:,5)=wwageo_across return(all_winds) end