View source for NCL geostrophic

<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 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
</pre>

[[Category:NCL]]
[[Category:WRF]]