[spam] suggestions for reducing the run time for the script (ascii writing)i

Dear NCL users,

I have written a NCL script in version 5.1.1. It uses Cygwin in windows.

The script reads IPCC data (for the whole world, at daily time step) in netcdf format and writes a region (Eastern North America bounded by lat=25 to 50 deg. N; lon=275 to 300 deg E) to ascii file.

The ascii file has matrix comprising of coumns the time-step, lat, lon & value of the variable. The sample output file is

time lat lon var. val(in this case precipitation)
"0000058804.5 26.5 275.6 0.06
0000058805.5 26.5 275.6 0.38"

The script runs without error., however it takes a very long time (~10 hrs) for a 10 yr daily data (for 1322*9*9 lines). I think the 4 do loops is reducing the speed. The 4 loops: one loop for the 145 files i need to read, the rest 3 are for extracting the value lat, lon & time.

The script is pasted below. Can this script be improved for efficiency to run faster? I am new to ncl . Any suggestion in this direction would be very useful.

;------program to subset a netcdf file for ENA----
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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin

 ;=========================================;
 ; step 1: get the list of netcdf files from a directory
 ;=========================================;

 dir = "E:/IPCC_data_copy/20c3m_atm_da/pr"

 dblquote = integertochar(34)

 fili = systemfunc("find " + dir + "/* -type f -name " + dblquote + "*.nc" + dblquote + " -print")

 ;t(fili)
   nfili = dimsizes(fili)
 ;print(nfili)

 ;=========================================;
 ; step 2: loop to read one netcdf file at a time from a list of files
 ;=========================================;

 do nf=0,nfili-1
      in= addfile (fili(nf), "r")
      orig = in->pr
      time = in->time
      lat = in->lat
      lon = in->lon
      print(fili(nf))

 ;=========================================;
 ; step 3: create output file name(fnameo);
 ;=========================================;

  fout=new(3,"string") ; defing fout ass an array of strings
  fout(0)="E:/IPCC_data_extract_ENA/20c3m_atm_da" ; adding a folder name
  fex=stringtochar(fili(nf));
  fout(1)=chartostring(fex(30:dimsizes(fex)-5)) ; extracting a portion of the file name from input file
  fout(2)="_ENA_out.txt"
  fNameO=str_concat(fout)
  print(fNameO)

 ;=========================================;
 ; step 5: for a specified region get the list of lat (lat_v),
 ; lon(lon_v); its dimension values(lon_d),(lat_d)
 ;=========================================;
 lats = 25 ; lower limit of lat
 latn = 50 ; upper limit of lon
 lonw = 275 ; left limit of lon
 lone = 300 ; right limit of lon
 lat_v = lat(ind(lat.gt.lats .and. lat.lt.latn))
 lat_d = ind(lat.gt.lats .and. lat.lt.latn)
 nlat= dimsizes(lat_d)

 lon_v = lon(ind(lon.gt.lonw .and. lon.lt.lone))
 lon_d = ind(lon.gt.lonw .and. lon.lt.lone)
 nlon=dimsizes(lon_d)
 ntim=dimsizes(time)

 ;=========================================;
 ; step 6: subset the value of variable for a specified region get
 ; the list of lat (lat_v),
 ; lon(lon_v); its dimension values(lon_d),(lat_d)
 ;=========================================;

 orig_ss = orig(:,{lat_v(0):lat_v(nlat-1)},{lon_v(0):lon_v(nlon-1)})

 ;============================================================================;
 ; step 7: write the lat, lon, time, var values for each output file (fName0)
 ; 3 loops to read every value of lat, lon, time & print to ascii file
 ;============================================================================;

    ; printVarSummary(in) ; prints the header information
 data=new(1,"string")
     do nl=0,nlat-1
      do ml=0,nlon-1
       do nt=0,ntim-1
        npt =0
        data(npt) = sprintf("%0.5i", (npt+1) )
        data(npt) = data(npt) + sprintf("%7.1f ",time(nt))
        data(npt) = data(npt) + sprintf("%7.1f ",lat_v(nl))
        data(npt) = data(npt) + sprintf("%4.1f ",lon_v(ml))
        data(npt) = data(npt) + sprintf("%6.2f ", orig_ss(nt,nl,ml)*86400)
        data_app= data(npt)
        ;filename = "C:/anandhi/Part2_GCM_evaluation_methods/analysis_GCM_evaluation/foo_ra1.txt"
        system("echo " + data_app+ ">>" + fNameO)

       end do
      end do
     end do
     ;asciiwrite (fNameO , data)

     end do

end
;-----------------------------
thank you all in advance for the support.

regards
anandhi
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sat Sep 05 2009 - 10:01:10 MDT