; BEFORE YOU RUN THIS SCRIPT, you need to set things up. ; Search for "<<< SET THINGS HERE" to see where to do this. ; (12/30/2004) ; ERA40 eight single level fields (ds117.0). ; The fields are: CI,SSTK,STL1,SD,SDOR,STL2,LSM,STL3,STL4. ; Does conversion from GRIB to NetCDF. This code derived from ds117.2.0.ncl. ; ; NOTES ; ----- ; 1. LSM is constant in time. We keep the time dimension in the output. ; 2. Fields CI and SSTK have missing values over land and are not ; interpolated correctly by NCL. Hence, we use a more round-about, ; but correct, method of getting the fields. This involves wgrib and ; calls to EMOSLIB. (Kudos to Dave Stepaniak in DSS for coming up ; with this method and the Fortran source code for red2reg -- see ; /home/tempest/mai/bin/src/prog_redn80_regn80.f.) ; 3. 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. ; 4. This version runs on tempest only! load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "~mai/ncl/functions/getEndDate.ncl" ; --------------- Subroutines --------------- undef ("date2flnm") procedure date2flnm(i:integer, msn2d:string, ln2d:string, indx2: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.0. ; See the following URL for the name tables: ; ; http://dss.ucar.edu/datasets/ds117.0/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 surface dataset (ds117.0) msn2d = "U0"+(12*(iyr-1833)+imo-7) ; works for Sep 1957 -> Aug 2002 ; ..Local file name for surface dataset ln2d = "e4oper.an.sfc."+cyr+cmo+"01" ; ..Index for 0Z of this date for surface dataset indx2 = 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,"level")) then delete(x@level) end if if (isatt(x,"initial_time0")) then delete(x@initial_time0) end if if (isatt(x,"gwt")) then delete(x@gwt) 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 ; ..Absolute paths of executables RED2REG = "/home/"+HOST+"/mai/bin/red2reg" WGRIB = "/usr/local/dss/bin/wgrib" ; <<< 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") msn2d = new( 1, "string") indx2 = new(nfils, "integer") do nf = 0, nfils-1 date2flnm(idate(nf), msn2d, fils(nf), indx2(nf)) system("msrcp -n mss:/DSS/"+msn2d+" "+diri+fils(nf)) end do print(fils) print(indx2) ; ..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 files grib_in = addfile(diri+fils(nf)+".grb","r") wallClockElapseTime(wcStrt, "Read/open file: "+fils(nf), 0) ; ..Read and create time ntSt2 = indx2(nf)+ntStrt ntLa2 = indx2(nf)+ntLast+4*(ndays(nf)-1) stime = grib_in->initial_time0(ntSt2:ntLa2:Stride) time = grib_in->initial_time0_hours(ntSt2:ntLa2:Stride) ntim = dimsizes(time) DATE = doubletointeger(grib_in->initial_time0_encoded(ntSt2:ntLa2:Stride)) date = DATE/100 print(date) if (date(ntStrt) .eq. idate(nf)) then print("Date in surface file matches specified date") else print("STOP -- date in surface 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" ; ..Spectral truncation parameters ntrm = trunc ntrn = trunc ntrk = trunc ntrm@long_name = "spectral truncation parameter M" ntrn@long_name = "spectral truncation parameter N" ntrk@long_name = "spectral truncation parameter K" if (nf.eq.0) then print (stime) print (time) print (datesec) print (ntrm) print (ntrn) print (ntrk) end if ; ..Read standard scalar variables: GRIB lats run N->S filp = fils(nf) ; CIQ = grib_in->CI_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) system("cd "+diri+"; "+WGRIB+" -s "+filp+" | egrep ':CI:' | "\ +WGRIB+" -i "+filp+" -grib -o "+filp+".CI.redn80.grb > /dev/null ; "\ +RED2REG+" "+filp+".CI") grib_ic = addfile(diri+filp+".CI.regn80.grb","r") CIQ = grib_ic->CI_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) CIQ = CIQ(:,::-1,:) CIQ@_FillValue = 1.e20 printVarSummary(CIQ) system("cd "+diri+"; rm -f "+filp+".CI.redn80.grb "+filp+".CI.regn80.grb") ; SSTKQ = grib_in->SSTK_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) system("cd "+diri+"; "+WGRIB+" -s "+filp+" | egrep ':SSTK:' | "\ +WGRIB+" -i "+filp+" -grib -o "+filp+".SSTK.redn80.grb > /dev/null ; "\ +RED2REG+" "+filp+".SSTK") grib_is = addfile(diri+filp+".SSTK.regn80.grb","r") SSTKQ = grib_is->SSTK_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) SSTKQ = SSTKQ(:,::-1,:) SSTKQ@_FillValue = 1.e20 printVarSummary(SSTKQ) system("cd "+diri+"; rm -f "+filp+".SSTK.redn80.grb "+filp+".SSTK.regn80.grb") STL1Q = grib_in->STL1_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL1Q = STL1Q(:,::-1,:) STL1Q@_FillValue = 1.e20 printVarSummary(STL1Q) SDQ = grib_in->SD_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) SDQ = SDQ(:,::-1,:) SDQ@_FillValue = 1.e20 printVarSummary(SDQ) SDORQ = grib_in->SDOR_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) SDORQ = SDORQ(:,::-1,:) SDORQ@_FillValue = 1.e20 printVarSummary(SDORQ) STL2Q = grib_in->STL2_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL2Q = STL2Q(:,::-1,:) STL2Q@_FillValue = 1.e20 printVarSummary(STL2Q) LSMQ = grib_in->LSM_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) LSMQ = LSMQ(:,::-1,:) LSMQ@_FillValue = 1.e20 printVarSummary(LSMQ) STL3Q = grib_in->STL3_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL3Q = STL3Q(:,::-1,:) STL3Q@_FillValue = 1.e20 printVarSummary(STL3Q) STL4Q = grib_in->STL4_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL4Q = STL4Q(:,::-1,:) STL4Q@_FillValue = 1.e20 printVarSummary(STL4Q) delete(filp) ; ..Regrid standard scalar variables to target grid wcQ = systemfunc("date") CIt = linint2_Wrap(lonin,latin,CIQ ,True,lon,lat,0) ; have missing values CI = linint1_Wrap(lon,CIt,True,lon,0) SSTt = linint2_Wrap(lonin,latin,SSTKQ,True,lon,lat,0) ; SSTK = linint1_Wrap(lon,SSTt,True,lon,0) STL1 = g2gsh_Wrap(STL1Q, (/nlat,mlon/), trunc) SD = g2gsh_Wrap(SDQ, (/nlat,mlon/), trunc) SDOR = g2gsh_Wrap(SDORQ, (/nlat,mlon/), trunc) STL2 = g2gsh_Wrap(STL2Q, (/nlat,mlon/), trunc) LSM = linint2_Wrap(lonin,latin,LSMQ,True,lon,lat,0) ; not a good field to interp spectrally STL3 = g2gsh_Wrap(STL3Q, (/nlat,mlon/), trunc) STL4 = g2gsh_Wrap(STL4Q, (/nlat,mlon/), trunc) delete(CIQ) delete(CIt) delete(SSTKQ) delete(SSTt) delete(STL1Q) delete(SDQ) delete(SDORQ) delete(STL2Q) delete(LSMQ) delete(STL3Q) delete(STL4Q) wallClockElapseTime(wcQ, "Scalars: g2gsh_Wrap", 0) ; ..Fix up dimensions and metadata CI@_FillValue = 1.e20 CI!0 = "time" CI!1 = "lat" CI!2 = "lon" grib_meta(CI) printVarSummary(CI) SSTK!0 = "time" SSTK!1 = "lat" SSTK!2 = "lon" grib_meta(SSTK) printVarSummary(SSTK) STL1!0 = "time" STL1!1 = "lat" STL1!2 = "lon" grib_meta(STL1) printVarSummary(STL1) SD!0 = "time" SD!1 = "lat" SD!2 = "lon" grib_meta(SD) printVarSummary(SD) SDOR!0 = "time" SDOR!1 = "lat" SDOR!2 = "lon" SDOR@units = "m" grib_meta(SDOR) printVarSummary(SDOR) STL2!0 = "time" STL2!1 = "lat" STL2!2 = "lon" grib_meta(STL2) printVarSummary(STL2) LSM!0 = "time" LSM!1 = "lat" LSM!2 = "lon" grib_meta(LSM) printVarSummary(LSM) STL3!0 = "time" STL3!1 = "lat" STL3!2 = "lon" grib_meta(STL3) printVarSummary(STL3) STL4!0 = "time" STL4!1 = "lat" STL4!2 = "lon" grib_meta(STL4) printVarSummary(STL4) ; ..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: "+nl+\ " Some variables [eg: specific humidity] are directly available on a full "+nl+\ " grid while other variables [eg: sea-ice cover and SST] are on an N80 "+nl+\ " reduced grid. The expansion of this data leads to an nlat=160 / mlon=320"+nl+\ " T159 Gaussian grid. "+nl+\ " (1) Read file via NCL "+nl+\ " (2) Read/unpack Q directly (60,160,320) "+nl+\ " Use: g2gsh(Q , (/64,128/), 42) to get target grid "+nl+\ " (3) Read/unpack VAR directly (160,320), where "+nl+\ " VAR = STL1,SD,SDOR,STL2,LSM,STL3,STL4 "+nl+\ " Use: g2gsh(VAR , (/64,128/), 42) to get target grid "+nl+\ " (linint2 for LSM -- land fraction) "+nl+\ " (4) Read/unpack VAR (160,320) using wgrib and emoslib, where "+nl+\ " VAR = CI,SSTK (missing values not handled correctly by NCL) "+nl+\ " Use: linint2 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", "lat", "lon" /) dimSizes = (/ ntim , nlat , mlon /) dimUnlim = (/ True , False, False /) filedimdef(netcdf_out, dimNames, dimSizes, dimUnlim) filevardef (netcdf_out, "time", typeof(time), "time") filevarattdef(netcdf_out, "time", time) 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, "ntrm" , typeof(trunc) , "ncl_scalar") filevarattdef(netcdf_out, "ntrm" , ntrm) filevardef (netcdf_out, "ntrn" , typeof(trunc) , "ncl_scalar") filevarattdef(netcdf_out, "ntrn" , ntrn) filevardef (netcdf_out, "ntrk" , typeof(trunc) , "ncl_scalar") filevarattdef(netcdf_out, "ntrk" , ntrk) filevardef (netcdf_out, "CI" , typeof(CI) , getvardims(CI)) filevarattdef(netcdf_out, "CI" , CI) filevardef (netcdf_out, "SSTK" , typeof(SSTK) , getvardims(SSTK)) filevarattdef(netcdf_out, "SSTK" , SSTK) filevardef (netcdf_out, "STL1" , typeof(STL1) , getvardims(STL1)) filevarattdef(netcdf_out, "STL1" , STL1) filevardef (netcdf_out, "SD" , typeof(SD) , getvardims(SD)) filevarattdef(netcdf_out, "SD" , SD) filevardef (netcdf_out, "SDOR" , typeof(SDOR) , getvardims(SDOR)) filevarattdef(netcdf_out, "SDOR" , SDOR) filevardef (netcdf_out, "STL2" , typeof(STL2) , getvardims(STL2)) filevarattdef(netcdf_out, "STL2" , STL2) filevardef (netcdf_out, "LSM" , typeof(LSM) , getvardims(LSM)) filevarattdef(netcdf_out, "LSM" , LSM) filevardef (netcdf_out, "STL3" , typeof(STL3) , getvardims(STL3)) filevarattdef(netcdf_out, "STL3" , STL3) filevardef (netcdf_out, "STL4" , typeof(STL4) , getvardims(STL4)) filevarattdef(netcdf_out, "STL4" , STL4) wallClockElapseTime(wcDef, "Definition phase", 0) wcWrite = systemfunc("date") ; ..Write data values to predefined locations netcdf_out->time = (/ time /) netcdf_out->lat = (/ lat /) netcdf_out->lon = (/ lon /) netcdf_out->gw = (/ gwt /) netcdf_out->ntrm = (/ ntrm /) netcdf_out->ntrn = (/ ntrn /) netcdf_out->ntrk = (/ ntrk /) netcdf_out->date = (/ date /) netcdf_out->datesec = (/ datesec /) netcdf_out->CI = (/ CI /) netcdf_out->SSTK = (/ SSTK /) netcdf_out->STL1 = (/ STL1 /) netcdf_out->SD = (/ SD /) netcdf_out->SDOR = (/ SDOR /) netcdf_out->STL2 = (/ STL2 /) netcdf_out->LSM = (/ LSM /) netcdf_out->STL3 = (/ STL3 /) netcdf_out->STL4 = (/ STL4 /) wallClockElapseTime(wcWrite, "Write phase", 0) delete(stime) ; may change in next file delete(time) delete(DATE) delete(date) delete(datesec) delete(hh) delete(CI) delete(SSTK) delete(STL1) delete(SD) delete(SDOR) delete(STL2) delete(LSM) delete(STL3) delete(STL4) end do ; loop on files print("Normal termination of program ds117.0.ncl") end