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/contrib/cd_string.ncl" ;*************************************** begin ; define number of panel plots plot = new(2,graphic) n_t_vars=5 n_l_times = 0 ; ascii file a_asc = (/"gefs_time.file"/) ;---Read the data in as a small 1D array values = asciiread(a_asc,-1,"float") nvals = dimsizes(values) print("nvals = " + nvals) ; nx=dims(1) ; ny=dims(0) ;---Calculate ntime based on the number of values we read in. n_l_times = 39 print("n_times = " + n_l_times) ; ntxy = nblk+1 nvars = 5 nxy = 1 nblk = nxy * nvars ; create variables to hold the ascii file variables. ; ; temporary ny=1 nx=1 n_t_vars=5 time_array = new((/n_t_vars/),integer) year = new((/n_l_times,ny,nx/),integer) month = new((/n_l_times,ny,nx/),integer) day = new((/n_l_times,ny,nx/),integer) hour = new((/n_l_times,ny,nx/),integer) minute = new((/n_l_times,ny,nx/),integer) do it=0,n_l_times-1 nstart = it*(nblk+1)+1 nend = nstart+nblk-1 year(it,:,:) = toint(reshape(values(nstart:nend:nvars),(/ny,nx/))) month(it,:,:) = toint(reshape(values(nstart+1:nend:nvars),(/ny,nx/))) day(it,:,:) = toint(reshape(values(nstart+2:nend:nvars),(/ny,nx/))) hour(it,:,:) = toint(reshape(values(nstart+3:nend:nvars),(/ny,nx/))) minute(it,:,:) = toint(reshape(values(nstart+4:nend:nvars),(/ny,nx/))) end do print("year = " + year(:,0,0)) print("day = " + day(:,0,0)) nx = 720 ny = 361 ; slp = new((/ny,nx/),float) ; slp = new((/ny,nx/),float) ; u_ave = new((/ny,nx/),float) ; v_ave = new((/ny,nx/),float) DATADir = "/home/cust100021_vol1/barry/cust100021_vol2/GEFS" ; find directories with this date, etc cmd = "find " + DATADir + "/GEFS* -type d -print" dirWRF = systemfunc (cmd) print(dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "GEFS_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) print("n_dirWRF" + n_dirWRF) outfile="gefs_surface_map" wks = gsn_open_wks("pdf",outfile) ; open wk station ; http://www.ncl.ucar.edu/Document/Functions/Built-in/systemfunc.shtml ; to retrieve a list of files from another directory ;Date = "/home/cust100021_vol1/barry/cust100021_vol2/GEFS/GEFS_00" ;diri = DATADir+Date diri = dirWRF(0) + "/" all_files = systemfunc("ls " + diri + "ge*pgrb2a*grb") n_files = dimsizes(all_files) print("n_files = " + n_files) ;;; Now that we have identified the number of files and number of directories ptot = 0 ;do n=0,n_files -1,1 do n=0,1,1 print("n times steps = " + n) do i_dir = 0,n_dirWRF-1 ; do i_dir = 0,1,1 print("i_dir = " + i_dir) diri = dirWRF(i_dir) + "/" all_files = systemfunc("ls " + diri + "ge*pgrb2a*grb") print("all_files(n) = " + all_files(n)) filename = all_files(n) print ("filename = " + filename) a = addfile(filename,"r") fin = a ; print(fin) ;;; get variables: temperature, longitude, latitude t2m = fin->TMP_P1_L103_GLL0(:,:) dims2d = dimsizes(t2m) print("dims2d(0) = " + dims2d(0)) print("dims2d(1) = " + dims2d(1)) slp = fin->PRMSL_P1_L101_GLL0(:,:) p_wat = fin->PWAT_P1_L200_GLL0(:,:) if (i_dir.eq.0)then slp = 0 end if u = fin->UGRD_P1_L103_GLL0(:,:) v = fin->VGRD_P1_L103_GLL0(:,:) lon = fin->lon_0 lat = fin->lat_0 t2m = t2m-273.15 ; degrees if (i_dir.eq.0)then p_wat = t2m/n_dirWRF t2m_ave = t2m/n_dirWRF slp = slp/n_dirWRF u_ave = u/n_dirWRF v_ave = v/n_dirWRF end if if (i_dir.ne.0)then print("here") printMinMax (t2m, False) printMinMax (slp, False) slp = slp/n_dirWRF + slp p_wat = t2m/n_dirWRF + p_wat t2m_ave = t2m/n_dirWRF + t2m_ave ; to attach coordinate variables` t2m_ave!0="lat" t2m_ave!1="lon" u_ave = u/n_dirWRF + u_ave v_ave = v/n_dirWRF + v_ave end if print("t2m = " + t2m(100,100)) print("p_wat = " + p_wat(100,100)) ;;;; ;;; prepare ;;;; if (i_dir.eq.n_dirWRF-1)then ;;;; ;;; create plot (pdf) ;;;; gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") ; load color table ;;; Font Style res1 = True ; plots modification on res1@txFont = "Helvetica" res1@txFontQuality = "High" ;;; Title res1@tiMainString = "www.meteo-blog.net" ; add title ; res1@gsnLeftString = "T2m, SLP" res1@gsnLeftString = "" res1@gsnRightString = "" ;;; Map res1@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res1@mpDataSetName = "Earth..4" ; high resolution res1@mpOutlineBoundarySets = "National" ; National borders res1@mpNationalLineColor = "Black" ; National borders color ;;; T2m res1@cnFillOn = True ; turn on color fill res1@cnLinesOn = False ; turn of contour lines res1@gsnDraw = False ; do not draw the plot res1@gsnFrame = False ; do not advance the frame ;;; color fill tresholds and colors res1@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res1@cnLevels = (/ -26.,-24.,-22.,-20.,-18.,-16.,-14.,\ -12.,-10.,-8.,-6.,-4.,-2.,0.,\ 2.,4.,6.,8.,10.,12.,14.,16.,18.,\ 20.,22.,24.,26.,28.,30.,32.,34.,\ 36.,38.,40./) ; set levels res1@cnFillColors = (/ 10,12,14,15,17,19,22,25,31,35,39,43,48,53,\ 61,64,67,71,89,95,99,104,108,111,115,123,\ 130,135,143,148,150,155,159,163,165/) ; set the colors to be used ;;; label bars style ;;; if not set standard values are used. ;;; so you normally can delete these lines without any problems res1@cnInfoLabelOn = False ; turn off contour info label res1@lbAutoManage = False ; control label bar res1@pmLabelBarDisplayMode = "Always" ; turns on label bar res1@lbOrientation = "Vertical" res1@pmLabelBarSide = "Right" res1@lbLabelAutoStride = True res1@pmLabelBarWidthF = 0.1 res1@pmLabelBarHeightF = 0.5 res1@lbLabelFontHeightF = .010 res1@lbPerimOn = False res1@lbLabelFont = "Helvetica" ; label font res1@lbTitleOn = True ; turn on title res1@lbTitleString = "[~S~o~N~C]" ; title string res1@lbTitlePosition = "Top" ; title position res1@lbTitleFontHeightF= .010 ; make title smaller res1@lbTitleDirection = "Across" ; title direction res2 = True res2@gsnDraw = False ; do not draw the plot res2@gsnFrame = False ; do not advance the frame res2@cnConstFLabelPerimOn = False res2@cnLevelSelectionMode = "ManualLevels" res2@cnMinLevelValF = 960 res2@cnMaxLevelValF = 1040 res2@cnLevelSpacingF = 10 res2@cnLineLabelAngleF = 0.0 res2@cnLineLabelFontHeightF = 0.01 res2@cnLineLabelBackgroundColor = -1 res2@gsnContourLineThicknessesScale = 1 res2@tfDoNDCOverlay = True slp = slp*0.01 ; Plotting options for Wind Vectors res3 = True ; vector only resources res3@gsnDraw = False ; don't draw res3@gsnFrame = False ; don't advance frame res3@vcGlyphStyle = "CurlyVector" ; curly vectors res3@vcRefMagnitudeF = 20 ; define vector ref mag res3@vcMinDistanceF = 0.025 ; thin out windbarbs res3@vcRefLengthF = 0.045 ; define length of vec ref res3@gsnRightString = " " ; turn off right string res3@gsnLeftString = " " ; turn off left string res3@tiXAxisString = " " ; turn off axis label res3@vcRefAnnoOrthogonalPosF = -.535 ; move ref vector into plot vector_map = gsn_csm_vector (wks,u_ave,v_ave,res3) ;;; Map Projection (min and max Lon and Lat) res1@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom res1@mpProjection = "LambertConformal" ; res1@mpMinLatF = 20 ; res1@mpMaxLatF = 50 ; res1@mpMinLonF = 15 ; res1@mpMaxLonF = 50 res1@mpMinLatF = 10 res1@mpMaxLatF = 60 res1@mpMinLonF = 05 res1@mpMaxLonF = 60 res1@mpLimitMode="LatLon" printMinMax (slp, False) print("got to here 1") ; contour_t2 = gsn_csm_contour_map(wks,p_wat,res1) contour_t2 = gsn_csm_contour_map(wks,t2m_ave,res1) ; contour_t2 = gsn_csm_contour_map(wks,t2m,res1) contour_slp = gsn_csm_contour(wks,slp,res2) print("got to here 2") ;plot = gsn_csm_contour_map(wks,slp,res1) ; create plot print("finished drawing z") overlay(contour_t2,contour_slp) overlay(contour_t2,vector_map) print("finished overlay call") draw(contour_t2) ; plot = gsn_csm_contour_map(wks,slp,res1) ; create plot print("finished drawing overlay plot") frame(wks) end if end do ; end directory loop end do ; end file (or time step loop) end