load "./shapefile_utils.ncl" ;dir="/Users/imoleayogbode/Documents/RRPrd/data/ECMWF/" dir="/glade/scratch/gbode/RNEW_SCRATCH/ECMWF/" YEAR = (/1983, 1984, 1987, 1988/);(/1973, 1976 ,1977, 1978,1979,1980,1984, 1988, 1996, 1997, 1998/) ;--1988, 1973, 1976 , 1984, 1977, 1978, 1979, 1980, 1996, 1997 and 1998 Dailydata = new((/dimsizes(YEAR)*3,366/),float) fields = new(dimsizes(YEAR)*3+1, string) fields(0) = "Days";tochar("Days") f1 = 1 f2 = 2 f3 = 3 dUl = 0 dUm = 1 sU = 2 ;--====================================== ; print(pkanoUm); wtype = "png" wtype@wkWidth = 1500 ; Set the pixel size of PNG image. wtype@wkHeight = 1800 wks = gsn_open_wks (wtype,"../figures/RROnset4rmERAUwindContourJan_Dec") ; send graphics to PNG file plot = new(dimsizes(YEAR), graphic) do inYr = 0,dimsizes(YEAR)-1 ;f= addfile(dir + "uwnd."+YEAR(inYr)+".nc","r") f= addfile(dir + "EraUwind"+YEAR(inYr)+".nc","r") print("===Processing Year"+" "+YEAR(inYr)+"...") ;print(f) ;exit if (YEAR(inYr) .ge. 1979) then uera:=short2flt(f->u) uera:=uera(:,::-1,::-1,:) else uera:=f->u uera:=uera(:,::-1,::-1,:) end if uera!1 = "level" uera!2 = "lat" uera!3 = "lon" printVarSummary(uera) print(uera&level) ;exit ;printVarSummary(uera) ;print(uera&level((/0,3,6/))) ;--Split time into year, month, day and hours time2 := uera&time ymd2 := ut_calendar(time2, -3) ; year-month-day year2 := ymd2/1000000 md2 := ymd2 - year2*1000000 mon2 := md2/10000 dy2 := (md2 - mon2*10000)/100 hr2 := md2%100 ;--define the first two months of the year lipYr=31+29 ;--Jan + Feb LeapYear nlipYr=31+28 ;--Non Leap Year ;--extract 12hourly data i12hr := ind(hr2 .eq. 12) ;print(ymd2(i12hr)) ;--Calculate daily and 12hourly values ueraDay := calculate_daily_values(uera, "avg", 0, False) uera12hr := uera(i12hr, :,:,:) ;printVarSummary(ueraDay) ;printVarSummary(uera12hr) time3 := uera12hr&time ymd3 := ut_calendar(time3, -3) ; year-month-day year3 := ymd3/1000000 md3 := ymd3 - year3*1000000 mon3 := md3/10000 dy3 := (md3 - mon3*10000)/100 hr3 := md3%100 ;print(year2(:10)) ;print(mon2(:10)) ;print(dy2(:10)) ;print(hr2(:10)) ;exit ;--Station subset uKano := ueraDay(:,:,{6},{5}) Uera := ueraDay(:,:,{0:20},{-20:20}) ;printVarSummary(uKano) ;--Select the first two months of the year (Jan-Feb) i:=ind(mon3.eq.1 .or. mon3.eq.2) ; dimz:=dimsizes(uKano(i,:)) dimz:= dimsizes(Uera(i,:,:,:)) ; gtvard = getvardims(uKano) ; print(gtvard) ; nD = 0 ; TNAME = gtvard(nD) ; typically TNAME="time" ; TNAME = "time" ; TIME = uKano&$TNAME$ ; extract for convenience ; print(TIME) print(dimz) ;exit if (dimz(0).eq.lipYr) then print("This is a leap year") uKanoMD := uKano(:,:) ;--Extract data from January-December ntime := uKano&time(:) uMD := Uera(:,:,:,:) ;--Extract data from January-December ntime := Uera&time(:) else if (dimz(0).eq.nlipYr) then print("This is a non-leap year") uKanoMD := uKano(:,:) ntime := uKano&time(:) uMD := Uera(:,:,:,:) ;--Extract data from January-December ntime := Uera&time(:) else print("Please check input data") end if end if ;-- print(uKanoMD&level) ;--NCL index 0, 3, and 6 are 1000, 700 and 400hPa, respectively Ul := uMD(:,1,:,:) - uMD(:,0,:,:) ;--U700hPa - Usfc(1000hPa) Um := uMD(:,2,:,:) - uMD(:,1,:,:) ;--U400hPa - U700hPa ;printVarSummary(Ul) ;exit kanoUl := uKanoMD(:,1) - uKanoMD(:,0) ;--U700hPa - Usfc(1000hPa) kanoUm := uKanoMD(:,2) - uKanoMD(:,1) ;--U400hPa - U700hPa ;---Create a table for Ul, Um and Ul+Um for each year ; Dailydata(dUl,:) = kanoUl ; Dailydata(dUm,:) = kanoUm ; Dailydata(sU,:) = (/kanoUl - kanoUm/) dUl = dUl + 3 dUm = dUm + 3 sU = sU + 3 ;--Create header of the fields fields(f1) = YEAR(inYr)+"Ul" fields(f2) = YEAR(inYr)+"Um" fields(f3) = YEAR(inYr)+"sUlUm" f1 = f1 + 3 f2 = f2 + 3 f3 = f3 + 3 kanoUl!0="time" kanoUl&time := ntime kanoUm!0="time" kanoUm&time := ntime Ul!0="time" Ul&time := ntime Um!0="time" Um&time := ntime ; printVarSummary(kanoUl) ; printVarSummary(kanoUm) ;exit Opt=True Opt@nval_crit = 4 ; require at least 4 values for a segment statistic to be calculated. ; default is one (1) Opt@segment_length = 5 ; Pentad averages are the default nDim = 0 pkanoUl := calculate_segment_values(kanoUl, "avg", nDim, Opt) pkanoUm := calculate_segment_values(kanoUm, "avg", nDim, Opt) pUl := calculate_segment_values(Ul, "avg", nDim, Opt) pUm := calculate_segment_values(Um, "avg", nDim, Opt) ;--Weekly averages Opt@segment_length = 7 ; weekly averages are the default wkanoUl := calculate_segment_values(kanoUl, "avg", nDim, Opt) wkanoUm := calculate_segment_values(kanoUm, "avg", nDim, Opt) aa:=dimsizes(pkanoUl) print(aa) ;---Alternative method to calculate the weekly (N=7) average of wind shear is by Looping ptest := new(aa,float) n=0 do ii=0,aa-1 ptest(ii) = dim_avg(kanoUl(n:(n+4))) n = n + 5 end do ;--0:6;7:7+6, ; printVarSummary(pkanoUl) ; printVarSummary(pkanoUm) ; printVarSummary(ptest) ;print(pkanoUl(38:42)) ;print(ptest(38:42)) ;---To plot multiple lines, you must put them into a mulidimensional array. ;--Daily data data := new((/2,dimsizes(kanoUl)/),float) xlimdy := dimsizes(kanoUl) xvaldy := ispan(1,xlimdy,1) data(0,:) = kanoUl data(1,:) = kanoUm ;--Pentad datap := new((/2,dimsizes(pkanoUl)/),float) xlimp := dimsizes(pkanoUl) xvalp := ispan(1,xlimp,1) datap(0,:) = pkanoUl datap(1,:) = pkanoUm ;--Weekly dataw := new((/2,dimsizes(wkanoUl)/),float) xlimw := dimsizes(wkanoUl) xvalw := ispan(1,xlimw,1) dataw(0,:) = wkanoUl dataw(1,:) = wkanoUm ;--Determine when condition -20 <= Ul <= -5 ;-- 0 <= Um <= 10 ;--Daily print("===Daily===") do it = 0,xlimdy - 1 if (kanoUl(it) .le. -5 .and. kanoUm(it) .ge. -5) then Uc = dimsizes(kanoUl(:it));num(pkanoUl(:it)) print("The wind shear condition was satisfied on day"+" "+Uc) ; typ= typeof(Uc) ; print(typ) break end if end do print("The wind shear condition was first satisfied on day"+" "+Uc) do it = 0,xlimdy - 1 if (kanoUl(it) .le. -5 .and. kanoUm(it) .ge. -5 .and. kanoUl(it+1) .le. -5 .and. kanoUm(it+1) .ge. -5 \ .and. kanoUl(it+2) .le. -5 .and. kanoUm(it+2) .ge. -5) then Uc3 = dimsizes(kanoUl(:it+2)) print("The wind shear condition was satisfied 3 consecutive times on day"+" "+Uc3) break end if end do print("The wind shear condition was satisfied 3 consecutive times on day"+" "+Uc3) print("==========") ;--Pentad print("===Pentad===") do it = 0,xlimp - 1 if (pkanoUl(it) .le. -2.5 .and. pkanoUm(it) .ge. -4) then Ucp = dimsizes(pkanoUl(:it));num(pkanoUl(:it)) print("The wind shear condition was satisfied on pentad"+" "+Ucp) ; typ= typeof(Ucp) ; print(typ) break end if end do print("The wind shear condition was first satisfied on pentad"+" "+Ucp) do it = 0,xlimp - 1 if (pkanoUl(it) .le. -2.5 .and. pkanoUm(it) .ge. -4 .and. pkanoUl(it+1) .le. -2.5 .and. pkanoUm(it+1) .ge. -4 \ .and. pkanoUl(it+2) .le. -2.5 .and. pkanoUm(it+2) .ge. -4 \ .and. pkanoUl(it+3) .le. -2.5 .and. pkanoUm(it+3) .ge. -4) then Ucp4 = dimsizes(pkanoUl(:it+3));num(pkanoUl(:it)) print("The wind shear condition was satisfied 3 consecutive times on pentad"+" "+Ucp4) break end if end do print("The wind shear condition was satisfied 3 consecutive times on pentad"+" "+Ucp4) print("==========") ;============Spatial=============================== print("===Pentad===") dimzg = dimsizes(pUl) print(dimzg) UcP = new((/dimzg(1),dimzg(2)/), float) copy_VarCoords(Uera(0,0,:,:), UcP) ;printVarSummary(UcP) ;exit do iy = 0,dimzg(1)-1 do ix = 0,dimzg(2)-1 do it = 0,dimzg(0) - 1 if (pUl(it,iy,ix) .le. -2.5 .and. pUm(it,iy,ix) .ge. -4) then UcP(iy,ix) = dimsizes(pUl(:it,iy,ix));num(pkanoUl(:it)) break end if end do end do end do print("The wind shear condition was first satisfied on pentad"+" "+UcP({6},{5})) UcP4 = new((/dimzg(1),dimzg(2)/), float) copy_VarCoords(Uera(0,0,:,:), UcP4) do iy = 0,dimzg(1)-1 do ix = 0,dimzg(2)-1 do it = 0,dimzg(0) - 4 if (pUl(it,iy,ix) .le. -2.5 .and. pUm(it,iy,ix) .ge. -4 .and. pUl(it+1,iy,ix) .le. -2.5 .and. pUm(it+1,iy,ix) .ge. -4 \ .and. pUl(it+2,iy,ix) .le. -2.5 .and. pUm(it+2,iy,ix) .ge. -4\ .and. pUl(it+3,iy,ix) .le. -2.5 .and. pUm(it+3,iy,ix) .ge. -4) then UcP4(iy,ix) = dimsizes(pUl(:it+3,iy,ix));num(pkanoUl(:it)) break end if end do end do end do print("The wind shear condition was satisfied 4 consecutive pentads on pentad"+" "+UcP4({6},{5})) print("==========") ;---Weekly print("===Weekly===") do it = 0,xlimw - 1 if (wkanoUl(it) .le. -5 .and. wkanoUm(it) .ge. -5) then Ucw = dimsizes(wkanoUl(:it));num(wkanoUl(:it)) ; typ= typeof(Ucw) ; print(typ) break end if end do print("The wind shear condition was first satisfied on week"+" "+Ucw) do it = 0,xlimw - 1 if (wkanoUl(it) .le. -5 .and. wkanoUm(it) .ge. -5 .and. wkanoUl(it+1) .le. -5 .and. wkanoUm(it+1) .ge. -5 \ .and. wkanoUl(it+2) .le. -5 .and. wkanoUm(it+2) .ge. -5) then Ucw3 = dimsizes(wkanoUl(:it+2));num(wkanoUl(:it)) print("The wind shear condition was satisfied 3 consecutive times on week"+" "+Ucw3) break end if end do print("The wind shear condition was satisfied 3 consecutive times on week"+" "+Ucw3) print("==========") ;---Set plotting parameters res = True ; plot mods desired res@cnFillOn = True res@gsnDraw = False res@gsnFrame = False res@gsnAddCyclic = False res@lbLabelBarOn = False res@cnLinesOn = True res@cnLineColor = "blue" res@mpMinLonF = -20 ;-- min longitude res@mpMaxLonF = 15 ;-- max longitude res@mpMinLatF = 4 ;-- min latitude res@mpMaxLatF = 16 ;-- max latitude ; res@mpOutlineBoundarySets = "National" res@pmTickMarkDisplayMode = "Always" res@mpGeophysicalLineThicknessF = 3 res@mpNationalLineThicknessF = 2 res@tiMainFont = 22 res@tiXAxisFont = 22 res@tiYAxisFont = 22 res@tmXBLabelFont = 22 res@tmYLLabelFont = 22 res@tmYRLabelFont = 22 ;res@cnFillPalette = "default" res@cnLineLabelPlacementMode = "constant" ; res@cnLevelSelectionMode = "ManualLevels" ; manually specify contour levels ; res@cnMinLevelValF =2 ;toint(min(t_mask)) ;-- Minimum contour level ; res@cnMaxLevelValF =28; toint(max(t_mask)) ;-- Maximum contour level ; res@cnLevelSpacingF = 2 ;-- contour level spacing res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/7, 8 ,9 , 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35/) ;res@cnLevels = (/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 21, 22, 23, 24, 25/) res@cnFillColors = (/"white","cyan", "green","yellow","darkorange","red","magenta","purple",\ "SlateBlue", "Khaki", "OliveDrab","BurlyWood", "LightSalmon", "Coral", \ "HotPink", "MediumTurquoise", "DarkSeaGreen", "Peru", "Tomato", \ "Orchid","PapayaWhip", "PeachPuff", "RoyalBlue","IndianRed","PaleGreen", "MintCream", "LemonChiffon", "AliceBlue", "LightGrey"/) ;res@cnFillColors = (/0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 21, 22, 23, 24, 25, 26/) ; res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; res@cnLevels = (/1.,2.,4.,6.,8.,10.,12.,14,16,18,20,22,24,26,28,30,32,34,36,38,40/) ; set levels ; res@cnFillOn = True ; turn on color fill ; res@cnLinesOn = False ; turn off the contour lines ; res@cnFillColors = (/1,3,5,9,11,13,23,30,36,41,45,47,59,63,68,74,81,91,93,96,99/) ; set the colors to be used ; res@cnLineLabelsOn = False ; turn the line labels off ;---------------------------------------------------------------------- ; Mask the data with the Nigeria outline. Set attribute ; "keep" to True, to indicate you want to *keep* values inside the ; shapefile outline. ;---------------------------------------------------------------------- ;shp_filename = "NGA_adm_shp/NGA_adm0.shp" shp_filename = "WestCentralAfricaShp/wca_adm0.shp" ; shp_filename = "NGA_adm_shp/NIR_outline_SHP-2/NIR_outline.shp" print("======================================================================") print("Masking data against " + shp_filename) opt2 = True opt2@keep = True ; Keep values inside this shape UcP_mask = shapefile_mask_data(UcP,shp_filename,opt2) UcP4_mask = shapefile_mask_data(UcP4,shp_filename,opt2) printMinMax(UcP_mask,0) printMinMax(UcP4_mask,0) print("# Valid values = " + num(.not.ismissing(UcP_mask))) print("# Msg values = " + num(ismissing(UcP_mask))) print("# Valid values = " + num(.not.ismissing(UcP4_mask))) print("# Msg values = " + num(ismissing(UcP4_mask))) ; plot = new(3, graphic) res@tiMainOffsetYF = -0.006 res@tiMainString = "Zonal wind shear (U~B~C~N~P) for"+" "+YEAR(inYr) ; plot(0) = gsn_csm_contour_map(wks,UcP_mask,res) res@tiMainString = "Zonal wind shear (U~B~C~N~P~B~4~N~) for"+" "+YEAR(inYr) ; plot(1) = gsn_csm_contour_map(wks,UcP4_mask,res) res@tiMainString = "Rainfall Onset (U~B~C~N~P~B~4~N~ plus 3 Pentads) for"+" "+YEAR(inYr) ROnD = UcP4_mask + 3 ;--Rainfall onset date is equal to UcP4 + 3Pentads copy_VarCoords(UcP4_mask, ROnD) plot(inYr) = gsn_csm_contour_map(wks,ROnD,res) ;---Add West African shapefile outlines to the plot. lnres = True lnres@gsLineColor = "gray25" lnres@gsLineThicknessF = 2.0 ; id_mask = gsn_add_shapefile_polylines(wks,plot(0),shp_filename,lnres) ; id_mask2 = gsn_add_shapefile_polylines(wks,plot(1),shp_filename,lnres) id_mask3 = gsn_add_shapefile_polylines(wks,plot(inYr),shp_filename,lnres) ;*************************************** ; Add stations as polymarkers * ;*************************************** res_mark = True res_mark@gsMarkerIndex = 1 ; polymarker style res_mark@gsMarkerSizeF = 0.018 ; polymarker size res_mark@gsMarkerColor = "black" ; change marker color slat=(/12, 6.7150, 6.58, 6.33,7.8,9.92,5.60,13.48/) slon=(/8.5, -1.5910, 3.33, 5.63,6.73,8.9,-0.17,2.17/) snam=(/"Kano", "Kumasi", "Ikeja", "Benin", "Lokoja","Jos","Accra","Niamey"/) ; dum = gsn_add_polymarker(wks,plot(0),slon,slat,res_mark) ; dum2 = gsn_add_polymarker(wks,plot(1),slon,slat,res_mark) dum3 = gsn_add_polymarker(wks,plot(inYr),slon,slat,res_mark) print("Now plotting text labels") ;*************************** ; Plot some text labels * ;*************************** txres = True ; Now draw the 4-letter ICAO station identifiers txres@txFontHeightF = 0.008 txres@txJust = "TopCenter" txres@txPerimOn = True txres@txPerimColor = "Black" txres@txPerimThicknessF = 0.3 txres@txPerimSpaceF = 0.4 txres@txBackgroundFillColor = "White" ; Calculate label offset as a function of map domain size ; dy = (maxlat - minlat)/125. ; print(dy) ; Put the station labels onto the map as text strings ; txt = gsn_add_text(wks,plot(0),snam,slon,slat,txres) ; this is really slow ; txt2 = gsn_add_text(wks,plot(1),snam,slon,slat,txres) txt3 = gsn_add_text(wks,plot(inYr),snam,slon,slat,txres) end do ;--Panel pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True gsn_panel(wks,plot(:2),(/3,1/),pres) exit ; res@gsnDraw = False ; res@gsnFrame = False res@gsnMaximize = True ; res@tiMainString = "Low and Mid-Level Zonal Wind shear over Kano (12N and 8.5E) for"+" "+YEAR(inYr) ; add title res@tiMainFontHeightF = 0.012 res@tiXAxisString = "Pentads (March-December)" res@tiYAxisString = "Vertical Zonal wind shear (m/s)" res@tiXAxisFontHeightF = 0.010 res@tiYAxisFontHeightF = 0.010 ; Similiar resources are xyLineThicknessF and xyLineColor, ; which will affect all lines in the array. ; res@xyLineThicknesses = (/ 1.0, 2.0/) ; make second line thicker res@xyLineColors = (/"blue","red"/) ; change line color res@gsnYRefLine = (/-5., 0.,10./) ; four X reference lines res@gsnYRefLineDashPatterns = (/2, 0, 2/) ; four X reference lines res@gsnXRefLine = (/Ucp, Ucp3/) ; four X reference lines res@gsnXRefLineColor = (/"Pink", "Green"/) ; four X reference lines res@vpHeightF = 0.3 res@vpWidthF = 0.6 res@tmYROn = False res@tmYRBorderOn = False res@tmXTOn = False res@tmXTBorderOn = False res@trYMinF = -20 res@trYMaxF = 20 ;---Make plot plotp = gsn_csm_xy (wks,xvalp,datap,res) ; create plot ;--Add arrow txres = True txres@txFontHeightF = 0.010 txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotp,"~F34~.",Ucp,pkanoUl(Ucp-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotp,"Uc",Ucp,pkanoUl(Ucp-1)-3,txres) ; . = draw arrow pointing txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotp,"~F34~.",Ucp3,pkanoUl(Ucp3-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotp,"UcP4",Ucp3,pkanoUl(Ucp3-1)-3,txres) ; . = draw arrow pointing ;---Make plots for weekly res@tiXAxisString = "Weekly (March-December)" res@gsnYRefLine = (/-5., 0.,10./) ; four X reference lines res@gsnXRefLine = (/Ucw, Ucw3/) ; four X reference lines ;---Make plot plotw = gsn_csm_xy (wks,xvalw,dataw,res) ; create plot ;--Add arrow txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotw,"~F34~.",Ucw,wkanoUl(Ucw-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotw,"Uc",Ucw,wkanoUl(Ucw-1)-3,txres) ; . = draw arrow pointing txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotw,"~F34~.",Ucw3,wkanoUl(Ucw3-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotw,"UcW3",Ucw3,wkanoUl(Ucw3-1)-3,txres) ; . = draw arrow pointing ;---Make plots for daily res@tiXAxisString = "Days (March-December)" res@gsnYRefLine = (/-5., 0.,10./) ; four X reference lines res@gsnXRefLine = (/Uc, Uc3/) ; four X reference lines ;---Make plot plotd = gsn_csm_xy (wks,xvaldy,data,res) ; create plot ;--Add arrow txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotd,"~F34~.",Uc,kanoUl(Uc-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotd,"Uc",Uc,kanoUl(Uc-1)-3,txres) ; . = draw arrow pointing txres@txAngleF = 90. ; change angle of arrow (from right) textA=gsn_add_text(wks,plotd,"~F34~.",Uc3,kanoUl(Uc3-1)-1,txres) ; . = draw arrow pointing txres@txAngleF = 0. ; change angle of arrow (from right) textUc=gsn_add_text(wks,plotd,"Uc3",Uc3,kanoUl(Uc3-1)-3,txres) ; . = draw arrow pointing ;;;;;;----- ;draw(plot) ;frame(wks) ;--Panel resources pres = True pres@gsnMaximize = True pres@gsnPanelMainFontHeightF = 0.015 pres@gsnPanelMainString = "Low and Mid-Level Zonal Wind shear over Kano (12N and 8.5E) for"+" "+YEAR(inYr) gsn_panel(wks, (/plotd,plotp,plotw/), (/3,1/), pres) end do ;--End loop ;---WRITE DATA INTO CSV ;print(Dailydata(:,32)) ;print(Dailydata(305,1)) ;print(Dailydata(305,2)) ;---Set up CSV file and header information for the file csv_filename = "ZonalWindShearTable2.csv" system("rm -f " + csv_filename) ; Remove file in case it exists. ;---Create a header line for CSV file fields2 = (/fields(0),fields(1),fields(2),fields(3),fields(4),fields(5),fields(6),fields(7),fields(8),fields(9),fields(10),fields(11),\ fields(12),fields(13),fields(14),fields(15),fields(16),fields(17),fields(18),fields(19),fields(20),fields(21),fields(22),\ fields(23),fields(24),fields(25),fields(26),fields(27),fields(28),fields(29),fields(30),fields(31),fields(32),fields(33)/) ; dq = str_get_dq() ; fields2 = dq + fields2 + dq ; Pre/append quotes to field names header = [/str_join(fields2,",")/] ; Header is field names separated print(fields2) ; by commas. ; ; Format to use for writing each variable to CSV file. ; If you don't want spaces in CSV file, use the following ; format string: ; format = "%s,%g,%g,%g" ; ;format = "%i,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2f,%6.2,%6.2f,%6.2f,%6.2f" format = "%i,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g" ; ; Loop through each time step and desired list of lat/lon values, ; and write a single line of data to CSV file. ; write_table(csv_filename, "w", header, "%s") ; Write header to CSV file. iday =ispan(1,306,1) ;do ij =0,305 alist = [/iday(:),Dailydata(0,:),Dailydata(1,:), Dailydata(2,:), Dailydata(3,:), Dailydata(4,:)\ ,Dailydata(5,:),Dailydata(6,:), Dailydata(7,:), Dailydata(8,:), Dailydata(9,:)\ ,Dailydata(10,:),Dailydata(11,:), Dailydata(12,:), Dailydata(13,:), Dailydata(14,:)\ ,Dailydata(15,:),Dailydata(16,:), Dailydata(17,:), Dailydata(18,:), Dailydata(19,:)\ ,Dailydata(20,:),Dailydata(21,:), Dailydata(22,:), Dailydata(23,:), Dailydata(24,:)\ ,Dailydata(25,:),Dailydata(26,:), Dailydata(27,:), Dailydata(28,:), Dailydata(29,:)\ ,Dailydata(30,:),Dailydata(31,:), Dailydata(32,:)/] ; Store data to be written in a list. write_table(csv_filename, "a", alist, format) ; Write list to CSV file. ;end do exit