; BEFORE YOU RUN THIS SCRIPT, you need to set things up. ; Search for "<<< SET THINGS HERE" to see where to do this. ; (11/19/2004 1.) ; ERA40 on pressure levels (ds117.1) conversion from GRIB to NetCDF; ; variables are Q,T,U,V ; ; NOTES ; ----- ; 1. The input is (1) an 8-digit integer in the form of yyyymmdd (the ; starting date), (2) the number of days and (3) the hour code -- ; a single-digit integer in the range [0-4]. The meaning of the ; hour code is: ; 0->0Z ; 1->6Z ; 2->12Z ; 3->18Z ; 4->all four times ; So the user can specify all times or one time. No other options are ; provided for specifying the hour of the day. ; 2. This version runs on tempest only! load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "~mai/svn/datasets/ECMWF/cc2g_ntStrt_ntLast_Stride.ncl" ; spec to grid load "~mai/svn/datasets/ECMWF/ccvd2uvg_ntStrt_ntLast_Stride.ncl" ; vort/div to u/v load "~mai/ncl/functions/getEndDate.ncl" ; --------------- Subroutines --------------- undef ("date2flnm") procedure date2flnm(i:integer, msn3d:string, ln3d:string, indx3:integer) ; ; Convert 8-digit integer-encoded date (yyyymmdd) to corresponding MSS ; and local file names. Also return the indices (zero-based!) for 0Z of ; the date in these files. This routine only works for ds117.1. See the ; following URL for the name tables: ; ; http://dss.ucar.edu/datasets/ds117.1/MSS-file-list.html ; begin ; ..Extract year, month, day from 8-digit integer-encoded date iyr = i/10000 imo = i/100 % 100 idy = i % 100 cyr = sprinti("%0.4i",iyr) cmo = sprinti("%0.2i",imo) ; ..MSS file name for model level dataset (ds117.1) msn3d = "U0"+(12*(iyr-1788)+(imo-7)) ; works for Sep 1957 -> Aug 2002 ; ..Local file name for model level dataset ln3d = "e4oper.an.pl."+cyr+cmo+"01" ; ..Index for 0Z of this date for model level dataset indx3 = 4*(idy-1) end undef ("grib_meta") procedure grib_meta(x) ; Bookkeeping: remove/add some GRIB meta-data begin if (isatt(x,"center")) then delete(x@center) end if if (isatt(x,"forecast_time")) then delete(x@forecast_time) end if if (isatt(x,"parameter_number")) then delete(x@parameter_number) end if if (isatt(x,"grid_number")) then delete(x@grid_number) end if if (isatt(x,"level_indicator")) then delete(x@level_indicator) end if if (isatt(x,"initial_time0")) then delete(x@initial_time0) end if if (isatt(x,"gwt")) then delete(x@gwt) end if if (isatt(x,"lv_ISBL1")) then delete(x@lv_ISBL1) end if x@missing_value = x@_FillValue end ; --------------- Main block --------------- begin ; <<< SET THINGS HERE when porting to another machine >>> ; ..Environment variables HOST = getenv("HOST") USER = getenv("USER") PTMP = getenv("PTMP") if (ismissing(PTMP)) then PTMP = "/ptmp/" end if ; <<< END SET THINGS when porting to another machine >>> ; ..Read starting date, number of days, hour and working directory from ; environment variable: setenv DS117 "${idate};${ndays};${ihour};${wdir}" cla = getenv("DS117") clc = stringtochar(cla) lst = dimsizes(clc)-2 if (clc(lst) .eq. "/") then lst = lst-1 end if ; ..Find the three semicolons ib = 0 i = 0 do while (clc(i) .ne. ";") i = i+1 end do ie = i-1 i = i+1 do while (clc(i) .ne. ";") i = i+1 end do jb = ie+2 je = i-1 i = i+1 do while (clc(i) .ne. ";") i = i+1 end do kb = je+2 ke = i-1 lb = ke+2 le = lst ; ..Create the three needed integer arrays nfils = 1 do i = ib,ie if (clc(i) .eq. ",") then nfils = nfils+1 end if end do idate = new(nfils, "integer") ndays = new(nfils, "integer") ihour = new(nfils, "integer") ; ..Recover the list of dates and place in the integer array idate inb = ib i = ib knt = 0 do while (clc(i) .ne. ";") i = i+1 if (clc(i) .eq. ",") then ine = i-1 knt = knt+1 idate(knt-1) = stringtointeger(chartostring(clc(inb:ine))) inb = i+1 end if end do ine = i-1 knt = knt+1 idate(knt-1) = stringtointeger(chartostring(clc(inb:ine))) ; ..Recover the list of numbers of days and place in the integer array ndays inb = jb i = jb knt = 0 do while (clc(i) .ne. ";") i = i+1 if (clc(i) .eq. ",") then ine = i-1 knt = knt+1 ndays(knt-1) = stringtointeger(chartostring(clc(inb:ine))) inb = i+1 end if end do ine = i-1 knt = knt+1 ndays(knt-1) = stringtointeger(chartostring(clc(inb:ine))) ; ..Recover the list of hours and place in the integer array ihour inb = kb i = kb knt = 0 do while (clc(i) .ne. ";") i = i+1 if (clc(i) .eq. ",") then ine = i-1 knt = knt+1 ihour(knt-1) = stringtointeger(chartostring(clc(inb:ine))) inb = i+1 end if end do ine = i-1 knt = knt+1 ihour(knt-1) = stringtointeger(chartostring(clc(inb:ine))) ; ..Check that specified dates are actually present in the available data stp = False do ifil = 0, nfils-1 istrt = idate(ifil) iend = getEndDate(istrt, ndays(ifil)) if (istrt .lt. 19570901 .or. istrt .gt. 20020801) then print("START DATE: "+istrt+" Must have 19570901 <= start date <= 20020801") stp = True end if if (iend .lt. 19570901 .or. iend .gt. 20020801) then print("END DATE : "+iend +" Must have 19570901 <= end date <= 20020801") stp = True end if end do if (stp) then exit end if ; ..Set up working directory tdir = chartostring(clc(lb:le)) system ("mkdir -p "+tdir) wdir = tdir+"/" ; ..Load correct input files to disk (for a given time spec, all needed ; time slices must be contained within a single MSS file) diri = wdir if (isfilepresent(diri)) then print("Using existing working directory: "+diri) else print("Making new working directory: "+diri) system("mkdir -p "+diri) end if diro = diri fils = new(nfils, "string") msn3d = new( 1, "string") indx3 = new(nfils, "integer") do nf = 0, nfils-1 date2flnm(idate(nf), msn3d, fils(nf), indx3(nf)) system("msrcp -n mss:/DSS/"+msn3d+" "+diri+fils(nf)) end do print(fils) print(indx3) ; ..T159 input grid nlati = 160 mloni = 320 latin = latGau (nlati, "lat", "latitude", "degrees_north") lonin = lonGlobeF (mloni, "lon", "longitude", "degrees_east") ; ..Target Gaussian grid nlat = 64 mlon = 128 trunc = 42 lat = latGau (nlat, "lat", "latitude", "degrees_north") lon = lonGlobeF (mlon, "lon", "longitude", "degrees_east") gwt = latGauWgt (nlat, "lat", "gauss weights", "") ; ..Loop over all files do nf=0,nfils-1 ; ..Set index and stride of starting and ending times (0-3): ; 0 --> 0Z ; 1 --> 6Z ; 2 --> 12Z ; 3 --> 18Z if (ihour(nf) .eq. 4) then ; include all four hours of the day ntStrt = 0 ntLast = 3 Stride = 1 else if (ihour(nf) .ge. 0 .and. ihour(nf) .le. 3) then ntStrt = ihour(nf) ntLast = ihour(nf) Stride = 4 else print("STOP -- ihour out of range (0-4): "+ihour(nf)) exit() end if end if print ("==========> nf="+nf+" <==========") wcStrt = systemfunc("date") ; wall clock start ; ..Open grib file grib_in = addfile(diri+fils(nf)+".grb","r") wallClockElapseTime(wcStrt, "Read/open file: "+fils(nf), 0) ; ..Read and create time ntSt3 = indx3(nf)+ntStrt ntLa3 = indx3(nf)+ntLast+4*(ndays(nf)-1) stime = grib_in->initial_time0(ntSt3:ntLa3:Stride) time = grib_in->initial_time0_hours(ntSt3:ntLa3:Stride) ntim = dimsizes(time) DATE = doubletointeger(grib_in->initial_time0_encoded(ntSt3:ntLa3:Stride)) date = DATE/100 print(date) if (date(ntStrt) .eq. idate(nf)) then print("Date in model level file matches specified date") else print("STOP -- date in model level file does not match specified date") exit() end if date@long_name = "current date as 8 digit integer (YYYYMMDD)" date!0= "time" hh = (DATE-date*100) datesec = hh*3600 datesec!0 = "time" datesec@long_name = "current seconds of current date" datesec@units = "s" levi = grib_in->lv_ISBL1 ; this is integer on file klev = dimsizes(levi) lev = new ( klev, "float") ; want type float delete(lev@_FillValue) ; coord should not have _FillValue attribute lev = levi lev!0 = "lev" lev@positive = "down" if (nf.eq.0) then ; 1st file only [debug] print (stime) print (time) print (datesec) print (lev) end if ; ..Read standard scalar variables: GRIB lats run N->S QQ = grib_in->Q_GDS4_ISBL(ntSt3:ntLa3:Stride,:,:,:) QQ = QQ(:,:,::-1,:) ; make lat coords S->N QQ@_FillValue = 1.e20 printVarSummary(QQ) ; ..Regrid standard scalar variable to target grid wcQ = systemfunc("date") Q = g2gsh_Wrap(QQ, (/nlat,mlon/), trunc) delete(QQ) wallClockElapseTime(wcQ, "Q: scalar: g2gsh_Wrap", 0) ; ..Fix up dimensions and metadata Q!0 = "time" Q!1 = "lev" Q!2 = "lat" Q!3 = "lon" grib_meta(Q) ; remove/add GRIB meta-data printVarSummary(Q) ; ..Read complex scalar variables; convert to target grid wcT = systemfunc("date") T = cc2g_ntStrt_ntLast (grib_in,"T_GDS50_ISBL",ntSt3,ntLa3,Stride,(/nlat,mlon/),trunc) wallClockElapseTime(wcT, "T: scalar: cc2g_ntStrt_ntLast", 0) T@_FillValue = 1.e20 printVarSummary(T) grib_meta(T) printVarSummary(T) ; ..Read scalar variables (div and rel vort) and convert to vectors wcUV = systemfunc("date") UV = ccvd2uvg_ntStrt_ntLast(grib_in,ntSt3,ntLa3,Stride,(/nlat,mlon/),trunc,(/"D_GDS50_ISBL","VO_GDS50_ISBL"/)) wallClockElapseTime(wcUV, "UV", 0) U = UV(0,:,:,:,:) ; separate the two components V = UV(1,:,:,:,:) delete(UV) ; ..Meta info for each variable U@long_name = "zonal wind component" U@units = "m/s" U@_FillValue = 1.e20 U@missing_value = U@_FillValue V@long_name = "meridional wind component" V@units = "m/s" V@_FillValue = 1.e20 V@missing_value = V@_FillValue printVarSummary(U) printVarSummary(V) ; ..Open output NetCDF file; create parameters and output file flnmo = fils(nf)+".nc" system ("rm -f "+diro+flnmo) netcdf_out = addfile(diro+flnmo,"c") ; ..Assign file attributes wcDef = systemfunc("date") nl = inttochar(10) ; carriage return fAtt = True fAtt@title = "ERA40 GRIB-to-NetCDF" fAtt@story = nl+\ " ERA40 GRIB file: Some variables are represented by complex spherical "+nl+\ " harmonic coefs [eg: divergence, relative vorticity and temperature]. "+nl+\ " Other variables [eg: specific humidity] are directly available. This "+nl+\ " data is on an nlat=160/mlon=320 Gaussian grid. "+nl+\ " (1) Read file via NCL "+nl+\ " (2) Read/unpack Q directly (23,160,320) "+nl+\ " Use: g2gsh(Q , (/64,128/), 42) to get target grid "+nl+\ " (3) Read/unpack T complex coeffficients "+nl+\ " Perform spherical harmonic synthesis => T(160,320) "+nl+\ " Use: g2gsh(T, (/64,128/), 42) to get target grid "+nl+\ " (4) Read/unpack div/rvort complex coef => DV/RV(160,320) "+nl+\ " (a) Compute divergent (ud,vd) and rotational (ur,vr) wind "+nl+\ " components. "+nl+\ " (b) u = ur+ud ; v = vr+vd "+nl+\ " (c) Use g2gshv(u,v, U,V, 42) to get target grid "+nl fAtt@Conventions = "None" fAtt@source_file = fils(nf) fAtt@source_center = "European Center for Medium-Range Weather Forecasts - Reading" fAtt@creation_date = systemfunc ("date") fileattdef( netcdf_out, fAtt ) ; ..Predefine coordinate information dimNames = (/ "time", "lev", "lat", "lon" /) dimSizes = (/ ntim , klev , nlat , mlon /) dimUnlim = (/ True , False, False, False /) filedimdef(netcdf_out, dimNames, dimSizes, dimUnlim) filevardef (netcdf_out, "time", typeof(time), "time") filevarattdef(netcdf_out, "time", time) filevardef (netcdf_out, "lev" , typeof(lev) , "lev") filevarattdef(netcdf_out, "lev" , lev) filevardef (netcdf_out, "lat" , typeof(lat) , "lat") filevarattdef(netcdf_out, "lat" , lat) filevardef (netcdf_out, "lon" , typeof(lon) , "lon") filevarattdef(netcdf_out, "lon" , lon) ; ..Predefine variable sizes filevardef (netcdf_out, "date" , typeof(date) , getvardims(date)) filevarattdef(netcdf_out, "date" , date) filevardef (netcdf_out, "datesec", typeof(datesec), getvardims(datesec)) filevarattdef(netcdf_out, "datesec", datesec) filevardef (netcdf_out, "gw" , typeof(gwt) , getvardims(gwt)) filevarattdef(netcdf_out, "gw" , gwt) filevardef (netcdf_out, "Q" , typeof(Q) , getvardims(Q)) filevarattdef(netcdf_out, "Q" , Q) filevardef (netcdf_out, "T" , typeof(T) , getvardims(T)) filevarattdef(netcdf_out, "T" , T) filevardef (netcdf_out, "U" , typeof(U) , getvardims(U)) filevarattdef(netcdf_out, "U" , U) filevardef (netcdf_out, "V" , typeof(V) , getvardims(V)) filevarattdef(netcdf_out, "V" , V) wallClockElapseTime(wcDef, "Definition phase took ", 0) wcWrite = systemfunc("date") ; ..Write data values to predefined locations wcStrt = systemfunc("date") netcdf_out->time = (/ time /) netcdf_out->lev = (/ lev /) netcdf_out->lat = (/ lat /) netcdf_out->lon = (/ lon /) netcdf_out->gw = (/ gwt /) netcdf_out->date = (/ date /) netcdf_out->datesec = (/ datesec /) netcdf_out->Q = (/ Q /) netcdf_out->T = (/ T /) netcdf_out->U = (/ U /) netcdf_out->V = (/ V /) wallClockElapseTime(wcWrite, "Write phase took ", 0) delete(levi) delete(lev) delete(stime) delete(time) delete(DATE) delete(date) delete(datesec) delete(hh) delete(Q) delete(T) delete(U) delete(V) end do ; loop on files print("Normal termination of program ds117.1.ncl") end