undef("calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon") function calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon) begin ; printVarSummary(th_grid) ; printMinMax(th_grid,False) del_th = grad_latlon_cfd(th_regrid,dst_lat,dst_lon,False,False) dTdy = del_th[0] ; latitude gradient of var_regrid dTdx = del_th[1] ; longitude gradient dY_min = min(dTdy) dX_min = min(dTdx) print("dY_min = " + dY_min) print("dX_min = " + dX_min) dTdx@units = "K/m" dTdy@units = "K/m" ; barry del_u = grad_latlon_cfd(u_regrid,dst_lat,dst_lon,False,False) dudy = del_u[0] ; latitude gradient of var_regrid dudx = del_u[1] ; longitude gradient del_v = grad_latlon_cfd(v_regrid,dst_lat,dst_lon,False,False) dvdy = del_v[0] ; latitude gradient of var_regrid dvdx = del_v[1] ; longitude gradient ; Do the calculation with those gradient variables mag_del_theta = sqrt(dTdx^2 + dTdy^2) ; Magnitude of del theta mag_del_theta@_FillValue = 0 mag_min = min(mag_del_theta) print ("mag_min = " + mag_min) ;rh90 = where(rh_plane.gt.90,rh_plane,rh_plane@_FillValue) ; only want to shade RH values greater than 90% ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Other formulation of FGEN ; Copied from GEMPAK appendix B1, page 6 conversion = 10^5 * 10800 ; Scale to SPC Mesoanalysis units of K/100 km/3 hr divergence = dudx + dvdy shearing = dudy + dvdy stretching = dudx - dvdy dims_mag=dimsizes(mag_del_theta) print("dims mag = " + dims_mag) min_mag = min(mag_del_theta) print("min_mag = " + min_mag) stretching@_FillValue = 0 delta = 0.5 * atan(shearing/stretching) beta = asin((-cos(delta) * dTdx - sin(delta) * dTdy) / mag_del_theta) resultant_def = sqrt(shearing^2 + stretching^2) pett_fgen = 0.5 * conversion * mag_del_theta * (resultant_def * cos(2. * beta) - divergence) ; do i_smooth =1,10 ; pett_fgen=smth9(pett_fgen, 0.50, 0.25, False) ; Heavy local smoothing ; ; end do ; wrf_smooth_2d( pett_fgen, 3 ) pett_fgen@units = "K/100 km/3 hr" pett_fgen@short_name = "Petterssen Frontogenesis" ; delete (opts) return(pett_fgen) end undef("set_res1_list") function set_res1_list() begin res = True res@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res@mpDataSetName = "Earth..4" ; high resolution res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; National borders res@mpNationalLineColor = "Black" ; National borders color res@mpNationalLineThicknessF =16.0 res@mpGeophysicalLineThicknessF =12.0 res@mpUSStateLineThicknessF =12.0 ; res@gsnDraw = False ; don't draw ; res@gsnFrame = False ; don't advance frame res@cnFillOn = True ; turn on color res@cnLinesOn = False ; no contour lines ; res@gsnLeftString = "500 mb Wind Vectors" ; change left string res@gsnLeftString = "500 mb Heights" res@gsnRightString = "m" ; assign right string res@mpFillOn = False ; no map fill res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; res@cnLevels = ispan(5200,5800,50) ; set the contour levels cmap = read_colormap_file("NCV_rainbow2") ; cmap := cmap(50:, :) ; cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap ; set color map res@cnLevels = (/ -4.,-3.5,-3.,-2.5,-2.0,-1.5,-1.,\ -0.5,0.,0.5,1.0,1.5,2.0,2.5,\ 3.0,3.5,4.0/) res@cnFillColors = (/ 10,14,17,22,31,39,48,53,\ 61,67,89,99,108,115,123,\ 130,135,143,148,155,163,165/) ; set the colors to be used ; res1 label bar res@cnInfoLabelOn = False ; turn off contour info label res@lbAutoManage = False ; control label bar res@pmLabelBarDisplayMode = "Always" ; turns on label bar res@lbOrientation = "Vertical" res@pmLabelBarSide = "Right" res@lbLabelAutoStride = True res@pmLabelBarWidthF = 0.20 res@pmLabelBarHeightF = 0.7 res@lbLabelFontHeightF = .025 res@lbPerimOn = False res@lbLabelFont = "Helvetica" ; label font res@lbTitleOn = True ; turn on title ; res@lbTitleString = "Heights [~S~o~N~C]" ; title string res@lbTitleString = "Heights [m]" ; title string res@lbTitlePosition = "Top" ; title position res@lbTitleFontHeightF= .020 ; make title smaller res@lbTitleDirection = "Across" ; title direction res@gsnAddCyclic = False ; data already has cyclic point return(res) end begin plots = new(4,graphic) ;in_file = "heights.nc" ;in_file = "heights_b.nc" ;in_file = "heights_cu1.nc" ;in_file = "heights_b_nomp.nc" ;in_file = "heights_7.nc" in_file = "fgen_vars.nc" a = addfile(in_file,"r") xi = a->xi yi = a->yi xo = a->xo yo = a->yo xi_nam = a->xi_nam yi_nam = a->yi_nam z_gefs_ave = a->z_gefs_ave u_gefs_ave = a->u_gefs_ave v_gefs_ave = a->v_gefs_ave th_gefs_ave = a->th_gefs_ave z_wrf_ave = a->z_wrf_ave u_wrf_ave = a->u_wrf_ave v_wrf_ave = a->v_wrf_ave th_wrf_ave = a->th_wrf_ave z_wrf_oth = a->z_wrf_oth u_wrf_oth = a->u_wrf_oth v_wrf_oth = a->v_wrf_oth th_wrf_oth = a->th_wrf_oth z_nam = a->z_nam u_nam = a->u_nam v_nam = a->v_nam th_nam = a->th_nam ; fi = a->fi ; z_nam = a->z_nam ; fi_oth = a->fi_oth xo_min = xo xo_min = 180 xo = xo-xo_min xlat2d= a->xlat2d printVarSummary(xlat2d) printMinMax(z_gefs_ave,False) ; xoo is latitude xoo = xo(180:250) ; yoo is longitude yoo = yo(76:130:-1) ; print(xoo) outfile="fgen_maps_22_nocump" wtype = "png" ; wtype = "pdf" wtype@wkWidth = 7500 ; produces a better looking PNG wtype@wkHeight = 7500 wks = gsn_open_wks(wtype,outfile) gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") res1 = set_res1_list() res1@tiMainString = "GEFS Mean" ;; Map Projection (min and max Lon and Lat) lat_beg = 27 lat_end = 47 lon_beg = -82 + 0 lon_end = -58 + 00 res1@mpMinLatF = lat_beg res1@mpMaxLatF = lat_end res1@mpMinLonF = lon_beg res1@mpMaxLonF = lon_end res1@mpLimitMode="LatLon" ; dims2d = dimsizes(th_regrid) ; lat2d = conform_dims(dims2d,xoo,1) ; lon2d = conform_dims(dims2d,yoo,0) th_regrid =th_gefs_ave(76:130,180:250) u_regrid :=u_gefs_ave(76:130,180:250) v_regrid :=v_gefs_ave(76:130,180:250) dst_lat = yoo dst_lon = xoo fgen_gefs= calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon) printVarSummary(fgen_gefs) copy_VarCoords(th_regrid,fgen_gefs) printVarSummary(fgen_gefs) print(x) fgen_gefs@long_name = "height" fgen_gefs@units = "K/100 km/3 hr" printMinMax(fgen_gefs,False) plot_1 = gsn_csm_contour_map(wks,fgen_gefs,res1) ;print(x) plots(1)= plot_1 ;fi@long_name = "height" ;fi@units = "(m)" ; yi@lat_name = "observed latitudes" ; yi@units = "degrees_north" ; xi@long_name = "observed longitudes" ; xi@units = "degrees_east" printVarSummary(xoo) printVarSummary(yoo) ; printVarSummary(fi) ; z_grid = area_hi2lores_Wrap (xi,yi,z_wrf_ave, False, 1, xoo,yoo, False) ; z_wrf_oth = area_hi2lores_Wrap (xi,yi,z_wrf_oth, False, 1, xoo,yoo, False) th_grid = area_hi2lores_Wrap (xi,yi,th_wrf_ave, False, 1, xoo,yoo, False) u_grid = area_hi2lores_Wrap (xi,yi,u_wrf_ave, False, 1, xoo,yoo, False) v_grid = area_hi2lores_Wrap (xi,yi,v_wrf_ave, False, 1, xoo,yoo, False) th_regrid :=th_grid u_regrid :=u_grid v_regrid :=v_grid dst_lat = yoo dst_lon = xoo fgen_wrf= calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon) copy_VarCoords(th_regrid,fgen_wrf) fgen_wrf@long_name = "height" fgen_wrf@units = "K/100 km/3 hr" printMinMax(fgen_wrf,False) th_grid := area_hi2lores_Wrap (xi,yi,th_wrf_oth, False, 1, xoo,yoo, False) u_grid := area_hi2lores_Wrap (xi,yi,u_wrf_oth, False, 1, xoo,yoo, False) v_grid := area_hi2lores_Wrap (xi,yi,v_wrf_oth, False, 1, xoo,yoo, False) th_regrid :=th_grid u_regrid :=u_grid v_regrid :=v_grid dst_lat = yoo dst_lon = xoo fgen_wrf_oth= calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon) fgen_wrf_oth@long_name = "height" fgen_wrf_oth@units = "K/100 km/3 hr" printMinMax(fgen_wrf_oth,False) copy_VarCoords(th_regrid,fgen_wrf_oth) th_grid := area_hi2lores_Wrap (xi_nam,yi_nam,th_nam, False, 1, xoo,yoo, False) u_grid := area_hi2lores_Wrap (xi_nam,yi_nam,u_nam, False, 1, xoo,yoo, False) v_grid := area_hi2lores_Wrap (xi_nam,yi_nam,v_nam, False, 1, xoo,yoo, False) printMinMax(th_grid,False) th_regrid :=th_grid u_regrid :=u_grid v_regrid :=v_grid dst_lat = yoo dst_lon = xoo fgen_nam= calc_fgen(th_regrid,u_regrid,v_regrid,dst_lat,dst_lon) fgen_nam@long_name = "height" fgen_nam@units = "K/100 km/3 hr" printMinMax(fgen_nam,False) copy_VarCoords(th_regrid,fgen_nam) ; print(x) ;z_wrf_ave_o = area_hi2lores_Wrap (xi,yi,fi, False, 1, xoo,yoo, False) ;z_wrf_oth_o = area_hi2lores_Wrap (xi,yi,fi_oth, False, 1, xoo,yoo, False) ;z_nam_o = area_hi2lores_Wrap (xi_nam,yi_nam,z_nam, False, 1, xoo,yoo, False) ;printVarSummary(z_wrf_ave_o) ;printVarSummary(xlat2d) ;copy_VarCoords(xlat2d,fi) ; z_wrf_ave_o!0 = "yoo" ; z_wrf_ave_o!1 = "xoo" ; z_wrf_ave_o&fi@units = "degrees_north" ; z_wrf_ave_o&fi@units = "degrees_east" ;printVarSummary(z_wrf_ave_o) ; z_wrf_ave_o@long_name = "observations of some sort" ; z_wrf_ave_o@units = "not sure" ; z_wrf_ave_o!0 = "lat" ; name the dimension ; z_wrf_ave_o&lat = xi ; assign the coordinate variable ; z_wrf_ave_o!1 = "lon" ; name the dimension ; z_wrf_ave_o&lat = yi ; assign the coordinate variable ; lat_1@long_name = "observed latitudes" ; lat_1@units = "degrees_north" ; lon_1@long_name = "observed longitudes" ;printVarSummary(z_wrf_ave_o) ;printMinMax(z_wrf_ave_o,False) ;dims = dimsizes(z_wrf_ave_o) ;print(dims) ;do j= 0,dims(0)-1 ;do i= 0,dims(1)-1 ;check = ismissing(z_wrf_ave_o(j,i)) ;if (check.eq.False)then ;print("z, j,i = " + xoo(i) + " " + yoo(j)+" "+z_wrf_ave_o(j,i)) ;end if ;end do ;end do res2 = res1 res2@tiMainString = "WRF Mean (M-T)" res2@tiMainString = "WRF Mean" res2@gsnAddCyclic = False plot_2 = gsn_csm_contour_map(wks,fgen_wrf,res2) plots(2)= plot_2 res3 = res1 ;res3@tiMainString = "WRF Mean (K-F)" ;res3@tiMainString = "WRF No Micro Mean" res3@tiMainString = "WRF No Micro Mean" res3@gsnAddCyclic = False plot_3 = gsn_csm_contour_map(wks,fgen_wrf_oth,res3) plots(3)= plot_3 res4 = res1 res4@tiMainString = "Observed" res4@gsnAddCyclic = False plot_4 = gsn_csm_contour_map(wks,fgen_nam,res4) plots(0)= plot_4 resP = True ; panel only resources resP@gsnMaximize = True ; maximize plots gsn_panel(wks,plots,(/2,2/),resP) end