; 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 2a.) ; ERA40 on model levels (ds117.2) conversion from GRIB to NetCDF with ; vertical interpolation to pressure surfaces; output 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.2. See the ; following URL for the name tables: ; ; http://dss.ucar.edu/datasets/ds117.2/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) ; ..Is this date in the first, second or third file of the model level data? jdy = (idy-1)/10 + 1 if (jdy .gt. 3) then jdy = 3 end if ; ..Offset in days for index in the model level data kdy = 10*(jdy-1) cdy = sprinti("%0.2i",(kdy+1)) ; ..MSS file name for model level dataset (ds117.2) msn3d = "U1"+(36*(iyr-1918)+3*(imo-6)+jdy) ; works for Sep 1957 -> Aug 2002 ; ..Local file name for model level dataset ln3d = "e4oper.an.ml."+cyr+cmo+cdy ; ..Index for 0Z of this date for model level dataset indx3 = 4*(idy-kdy-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_HYBL4")) then delete(x@lv_HYBL4) end if if (isatt(x,"lv_HYBL_i5")) then delete(x@lv_HYBL_i5) 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 >>> ; <<< SET THINGS HERE when changing output levels >>> ; ..Specify pressure levels levp = (/ 1., 2., 3., 5., 7., 10., 20., 30., 50., 70., 100., 150., \ 200., 250., 300., 400., 500., 600., 700., 775., 850., 925., 1000. /) ; <<< END SET THINGS HERE when changing output levels >>> ; ..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", "") ; ..Other inits for vertical interpolation to pressure surfaces intyp = 1 ; 1=linear, 2=log, 3=log-log kxtrp = True ; True=extrapolate p0 = 1. ; do not use P0 on input volume ; ..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" nlev = dimsizes(grib_in->lv_HYBL4) nlevp = dimsizes(levp) lev = new(nlevp,"float") ; want type float delete(lev@_FillValue) ; coord should not have _FillValue attribute lev = levp lev!0 = "lev" lev@positive = "down" if (nf.eq.0) then print (stime) print (time) print (datesec) print (lev) end if ; ..Read standard scalar variables: GRIB lats run N->S Qml = grib_in->Q_GDS4_HYBL(ntSt3:ntLa3:Stride,:,:,:) hyam = grib_in->lv_HYBL4_a hybm = grib_in->lv_HYBL4_b psfc = cc2g_ntStrt_ntLast (grib_in,"LNSP_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlati,mloni/),0) psfc = exp(psfc) Tml = cc2g_ntStrt_ntLast (grib_in,"T_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlati,mloni/),0) tbot = Tml(:,nlev-1,:,:) phis = cc2g_ntStrt_ntLast (grib_in,"Z_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlati,mloni/),0) varflg = 0 ; not T or Z QQ = vinth2p_ecmwf(Qml,hyam,hybm,lev,psfc,intyp,p0,1,kxtrp,varflg,tbot,phis) QQ = QQ(:,:,::-1,:) ; make lat coords S->N QQ@_FillValue = 1.e20 printVarSummary(QQ) delete(Qml) ; ..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") varflg = 1 ; 1 for T, -1 for Z TT = vinth2p_ecmwf(Tml,hyam,hybm,lev,psfc,intyp,p0,1,kxtrp,varflg,tbot,phis) TT = TT(:,:,::-1,:) TT@_FillValue = 1.e20 printVarSummary(TT) delete(Tml) ; ..Regrid standard scalar variable to target grid wcT = systemfunc("date") T = g2gsh_Wrap(TT, (/nlat,mlon/), trunc) delete(TT) wallClockElapseTime(wcT, "T: scalar: g2gsh_Wrap", 0) ; ..Fix up dimensions and metadata T!0 = "time" T!1 = "lev" T!2 = "lat" T!3 = "lon" grib_meta(T) printVarSummary(T) ; ..Read scalar variables (div and rel vort) and convert to vectors wcUV = systemfunc("date") UVml = ccvd2uvg_ntStrt_ntLast(grib_in,ntSt3,ntLa3,Stride,(/nlati,mloni/),0,(/"D_GDS50_HYBL","VO_GDS50_HYBL"/)) wallClockElapseTime(wcUV, "UV: vector: ccvd2uvg_ntStrt_ntLast", 0) Uml = UVml(0,:,:,:,:) ; separate variables Vml = UVml(1,:,:,:,:) varflg = 0 ; not T or Z UU = vinth2p_ecmwf(Uml,hyam,hybm,lev,psfc,intyp,p0,1,kxtrp,varflg,tbot,phis) VV = vinth2p_ecmwf(Vml,hyam,hybm,lev,psfc,intyp,p0,1,kxtrp,varflg,tbot,phis) UU = UU(:,:,::-1,:) VV = VV(:,:,::-1,:) ; ..Meta info for each variable UU@long_name = "zonal wind component" UU@units = "m/s" UU@_FillValue = 1.e20 UU@missing_value = UU@_FillValue VV@long_name = "meridional wind component" VV@units = "m/s" VV@_FillValue = 1.e20 VV@missing_value = VV@_FillValue printVarSummary(UU) printVarSummary(VV) delete(UVml) delete(Uml) delete(Vml) ; ..Regrid vector variable to target grid wcUV = systemfunc("date") U = new((/ntim,nlevp,nlat,mlon/),typeof(UU)) V = new((/ntim,nlevp,nlat,mlon/),typeof(VV)) g2gshv_Wrap(UU,VV,U,V,trunc) delete(UU) delete(VV) wallClockElapseTime(wcUV, "UV: vector: g2gshv_Wrap", 0) ; ..Fix up dimensions and metadata U!0 = "time" U!1 = "lev" U!2 = "lat" U!3 = "lon" V!0 = "time" V!1 = "lev" V!2 = "lat" V!3 = "lon" grib_meta(U) grib_meta(V) printVarSummary(U) printVarSummary(V) ; ..Open output NetCDF file; create parameters and output file flnmo = fils(nf)+"_vintp.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 (60,160,320) "+nl+\ " Use: vinth2p_ecmwf(Qml,...) to get pressure levels "+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: vinth2p_ecmwf(Tml,...) to get pressure levels "+nl+\ " Use: g2gsh(Tml,(/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: vinth2p_ecmwf(Uml,...) to get U onto pressure levels "+nl+\ " (d) Use: vinth2p_ecmwf(Vml,...) to get V onto pressure levels "+nl+\ " (e) 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 , nlevp, 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(lev) delete(stime) ; may change in next file delete(time) delete(DATE) delete(date) delete(datesec) delete(hh) delete(Q) delete(T) delete(U) delete(V) delete(hyam) delete(hybm) delete(psfc) delete(tbot) delete(phis) end do ; loop on files print("Normal termination of program ds117.2_vintp.ncl") end