;Written by Xiaoming Hu, CAPS, OU, Sept. 28, 2011 ; contact yuanfangcan@gmail.com load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" ;*********************************************************** begin ; filename = "hf005-04-soil-temp_XMHu.csv" filename = "hf005-04-soil-temp_XMHu_2.csv" ; expends to 2016 ;---------------------------------------------------------------------- ; This method uses "asciicread" to read the data ;---------------------------------------------------------------------- ;---Read the values in as 1D array of strings to get rows and columns. values_1d = asciiread(filename,-1,"string") ncols = dimsizes(str_split(values_1d(1),",")) nrows = dimsizes(values_1d)-1 print("This file has " + nrows + " rows and " + ncols + " columns.") ;---Reread values as integers and reshape at the same time. values_2d = onedtond(asciiread(filename,-1,"double"),(/nrows,ncols/)) opt = True opt@title = values_1d(0) opt@fout = filename+".txt" ; write_matrix(values_2d,"4f6.0,3f15.9",opt) res1= True res1@gsnFrame = False ; don't draw yet res1@gsnDraw = False ; don't advance frame ; res1@mpShapeMode = "FreeAspect" res1@vpHeightF =0.2 res1@vpWidthF = 0.9 res1@vpXF = 0.2 ; res1@vpYF =-0.9 res1@gsnPaperOrientation = "portrait" res1@gsnMaximize = True res1@gsnPaperMargin = 0.01 ; res1@tiYAxisString = "Soil temperature, ~S~o~N~C" res1@tiYAxisString = "Soil temperature, K" res1@tiXAxisString = "Day of Month";;"DOY"; "PM~B~2.5~N~, ~F33~m~F~g m~S~-3" ; label bottom axis with units ; res1@xyDashPatterns = (/16,4,16,2,15,3/) ; 0 is default (solid) res1@xyDashPatterns = (/0,4,16,2,15,3/) ; 0 is default (solid) ; res1@xyMarkLineModes = (/"Markers","MarkLines","MarkLines","MarkLines","MarkLines","MarkLines" /) res1@xyMarkLineMode = "Lines" res1@xyMarkers = (/16,16,16,16,16,16/) ; choose type of marker ; res1@xyMarkerColor = "Turquoise3" res1@xyLineColors = (/"red","Green","blue","Purple4","black","brown"/) res1@xyMarkerSizeF = 0.01 res1@xyLineThicknessF = 3 res1@tiXAxisFuncCode = "~" res1@pmTickMarkDisplayMode = "Always" ; turn on tickmarks ; res1@tmXTOn = False ; turn off top labels ; res1@tmYROn = False ; turn off right labels res1@tmXTLabelsOn = False res1@tmYRLabelsOn = False res1@lgPerimOn = False res1@lgLabelFontHeightF = 0.018 res1@tiXAxisFontHeightF= 0.018 res1@tiYAxisFontHeightF= 0.018 res1@tmXBLabelFontHeightF = 0.018 res1@tmYLLabelFontHeightF = 0.018 res1@tmYRMajorOutwardLengthF = 0 res1@tmYLMajorOutwardLengthF = 0 res1@tmXBMajorOutwardLengthF = 0 res1@tmXBMinorOutwardLengthF = 0 res1@pmLegendSide = "Top" ; Change location of res1@pmLegendParallelPosF = .82 ; move units right res1@pmLegendOrthogonalPosF = -.3 res1@pmLegendWidthF = 0.08 ; Change width and res1@pmLegendHeightF = 0.18 ; height of legend. ; res1@tiYAxisOffsetXF = 0.055 res1@tmYLLabelDeltaF = -0.8 res1@tmXBMode = "Explicit" res1@tmXBLabels = ispan(1,33,1) res1@tmXBValues = ispan(1,33,1) res1@trXMaxF = 32.5 res1@trXMinF =0.5 res1@tmYLFormat = "f" res1@tmYLLabelStride =1 res1@tmXBLabelStride = 2 res1@tmYLMinorOn = False ; res1@tmYLValues = ispan(330,370,5) ;(/6,12,18,24,30,36,42,42+6,42+12,18,24,30,36,42/) ; location of explicit labels ; res1@pmLegendDisplayMode = "Always" ; create a legend res1@tmYMajorGrid = True ; implement x grid res1@tmYMajorGridThicknessF = 1.0 ; 2.0 is default res1@tmYMajorGridLineDashPattern = 2 ; select short dash lines Year = values_2d(:,0) Month = values_2d(:,1) Day = values_2d(:,2) ; do iyear = 2009, 2010 do iyear = 2011, 2016 do imonth = 1, 12 ; do iyear = 1991, 1991 ; do imonth = 7, 12 if (imonth.lt.10) then CFSR_data_name = "soilt1.gdas."+iyear+"0"+imonth+".grb2_Tsoil0mhourly.txt" CFSR_data_name_outObs = "soilt1.gdas."+iyear+"0"+imonth+".grb2" else CFSR_data_name = "soilt1.gdas."+iyear+""+imonth+".grb2_Tsoil0mhourly.txt" CFSR_data_name_outObs = "soilt1.gdas."+iyear+""+imonth+".grb2" end if print("going to open /nsftor/xhu/public_html/WRF-UCM/CFSR/"+CFSR_data_name) CFSR_data = readAsciiTable("/nsftor/xhu/public_html/WRF-UCM/CFSR/"+CFSR_data_name, 2 , "double", 1) dims= dimsizes(CFSR_data) CFSR_data_2d = onedtond(CFSR_data(:,1), (/dimsizes(CFSR_data(:,1))/24, 24/)) CFSR_data_2d_dailyave = dim_avg(CFSR_data_2d) delete(CFSR_data_2d) ; CFSR_data_2d_dailyave_expand = new ((/dimsizes(CFSR_data_2d_dailyave),24/),double) ; do icolumn = 0, 23 ; CFSR_data_2d_dailyave_expand(:,icolumn)= CFSR_data_2d_dailyave ; end do ; CFSR_data_2d_only_dirunalVariation = CFSR_data(:,1) - ndtooned(CFSR_data_2d_dailyave_expand) CFSR_data_2d_only_dirunalVariation = CFSR_data(:,1) - dim_avg(CFSR_data(:,1)) ifile=(iyear-1992)*12+imonth-1 ; figurename = "TS_soil-temp_addDiurnal_"+ifile figurename = CFSR_data_name+"_basedONobs" res1@gsnRightString = iyear+"/"+imonth ;files(ifiles) wks = gsn_open_wks("eps" ,figurename) ; eps,pdf,x11,ncgm,eeps setvalues wks ; "wkColorMap" : "BlAqGrYeOrRevi200" "wkForegroundColor" : "Black";"Blue" end setvalues index= ind(Year.eq.iyear.and.Month.eq.imonth) if( mod(dims(0), dimsizes(index)).gt. 0) then print("something wrong!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ") end if do i3method = 4,6 ; values_expanded = new ((/dimsizes(index),24/),double) ; do icolumn = 0, 23 ; values_expanded(:,icolumn)= values_2d(index,i3method) ; end do ; values_expanded_1d = ndtooned(values_expanded) values_expanded_1d = linint1 (ispan(1,dimsizes(index),1),values_2d(index,i3method), True , fspan(1,dimsizes(index)+1,dimsizes(index)*24), 0) scaling = dim_stddev_Wrap(values_2d(index,i3method))/ dim_stddev_Wrap(CFSR_data_2d_dailyave) print("scaling factor="+scaling) ; obs_plus_CFSRV = values_expanded_1d + CFSR_data_2d_only_dirunalVariation*scaling + 273.15 obs_plus_CFSRV = (1-scaling)*(values_expanded_1d-dim_avg(values_expanded_1d)) + CFSR_data_2d_only_dirunalVariation*scaling + 273.15+ dim_avg(values_expanded_1d) ; obs_plus_CFSRV = values_expanded_1d+ 273.15 ; obs_plus_CFSRV = CFSR_data_2d_only_dirunalVariation if (i3method.eq.4) then obs_plus_CFSRV_3method = new((/3,dimsizes(obs_plus_CFSRV)/),double) obs_plus_CFSRV_3method_output = new((/dimsizes(obs_plus_CFSRV),4/),double) obs_plus_CFSRV_3method_output(:,0) = (/CFSR_data(:,0)/) end if obs_plus_CFSRV_3method(i3method-4,:) = (/obs_plus_CFSRV/) obs_plus_CFSRV_3method_output(:,i3method-4+1) = (/obs_plus_CFSRV/) end do ; i3method delete(CFSR_data_2d_dailyave) res1@tmXBMode = "Explicit" ; plot = gsn_csm_xy (wks,Day(index),transpose(values_2d(index,4:6))+273.15,res1) ; print(fspan(1,dimsizes(index)+1,dimsizes(index)*24)) plot = gsn_csm_xy (wks,fspan(1,dimsizes(index)+1,dimsizes(index)*24), obs_plus_CFSRV_3method,res1) opt = True ; opt@title = values_1d(0) opt@title = "date warming control controlD" opt@fout = "/nsftor/xhu/public_html/WRF-UCM/CFSR/"+CFSR_data_name_outObs+"_basedONobs.txt" write_matrix(obs_plus_CFSRV_3method_output,"f15.0,3f20.9",opt) delete(index) draw(plot) frame(wks) ; system("convert -trim -density 300 -resize 1000x400 "+figurename+".eps /nsftor/xhu/public_html/WRF-UCM/CFSR/"+figurename+".png") system("convert -trim -density 300 -resize 1000x400 "+figurename+".eps /nsftor/xhu/public_html/WRF-UCM/CFSR/"+figurename+".png") system("rm "+figurename+".eps") print("finish "+figurename+".eps") delete(CFSR_data) delete(obs_plus_CFSRV) delete(obs_plus_CFSRV_3method_output) delete(obs_plus_CFSRV_3method) delete(values_expanded_1d) delete(CFSR_data_2d_only_dirunalVariation) end do ; imonth end do ; iyear end