; 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 on model levels (ds117.2) + eight single level fields (ds117.0). ; The added single level fields are: CI,SSTK,STL1,SD,SDOR,STL2,LSM,STL3,STL4. ; Does conversion from GRIB to NetCDF. This code derived from ds117.2.ncl. ; ; NOTES ; ----- ; 1. Files from ds117.2 each contain about 10 days while the files from ; ds117.0 each contain exactly one month. To be more specific, the ; files in the ds117.2 dataset are arranged so that each month is in ; three files. The first two are always days 1-10 and 11-20 while the ; last file contains the remaining days of the month. ; 2. LSM is constant in time. We keep the time dimension in the output. ; 3. 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.) ; 4. 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. ; 5. 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, msn2d:string,\ ln3d:string, ln2d:string, indx3:integer, 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.2 and ; ds117.0. See the following two URLs for the name tables: ; ; http://dss.ucar.edu/datasets/ds117.2/MSS-file-list.html ; 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) ; ..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 ; ..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 model level dataset ln3d = "e4oper.an.ml."+cyr+cmo+cdy ; ..Local file name for surface dataset ln2d = "e4oper.an.sfc."+cyr+cmo+"01" ; ..Index for 0Z of this date for model level dataset indx3 = 4*(idy-kdy-1) ; ..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") filv = new(nfils, "string") msn3d = new( 1, "string") msn2d = new( 1, "string") indx3 = new(nfils, "integer") indx2 = new(nfils, "integer") do nf = 0, nfils-1 date2flnm(idate(nf), msn3d, msn2d, fils(nf), filv(nf), indx3(nf), indx2(nf)) system("msrcp -n mss:/DSS/"+msn3d+" "+diri+fils(nf)) system("msrcp -n mss:/DSS/"+msn2d+" "+diri+filv(nf)) end do print(fils) print(filv) print(indx3) 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", "") ; ..Vertical coordinates and attributes (layer midpoints) klev = 60 hyam = new (klev, "float") delete(hyam@_FillValue) hyam!0 = "lev" hyam@long_name = "hybrid A coefficient at layer midpoints" hybm = new (klev, "float") delete(hybm@_FillValue) hybm!0 = "lev" hybm@long_name = "hybrid B coefficient at layer midpoints" lev = new (klev, "float") ; want type float delete(lev@_FillValue) ; coord should not have _FillValue attribute lev!0 = "lev" lev@long_name = "hybrid level at midpoints (1000*(A+B))" lev@units = "level" lev@positive = "down" lev@standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" lev@formula_terms = "a: hyam b: hybm p0: P0 ps: PS" ; ..Vertical coordinates and attributes (layer interfaces) hyai = new (klev+1, "float") delete(hyai@_FillValue) hyai!0 = "ilev" hyai@long_name = "hybrid A coefficient at layer interfaces" hybi = new (klev+1, "float") delete(hybi@_FillValue) hybi!0 = "ilev" hybi@long_name = "hybrid B coefficient at layer interfaces" ilev = new (klev+1, "float") delete(ilev@_FillValue) ilev!0 = "ilev" ilev@long_name = "hybrid level at interfaces (1000*(A+B))" ilev@units = "level" ilev@positive = "down" ilev@standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ilev@formula_terms = "a: hyai b: hybi p0: P0 ps: PS" ; ..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") grib_iv = addfile(diri+filv(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) ntSt2 = indx2(nf)+ntStrt ntLa2 = indx2(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" DATEv = doubletointeger(grib_iv->initial_time0_encoded(ntSt2:ntLa2:Stride)) datev = DATEv/100 print(datev) if (datev(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 hh = (DATE-date*100) datesec = hh*3600 datesec!0 = "time" datesec@long_name = "current seconds of current date" datesec@units = "s" ; ..Read and fill in the CAM standard vertical coordinate variables hyam = (/ grib_in->lv_HYBL4_a /) hybm = (/ grib_in->lv_HYBL4_b /) hyai = (/ grib_in->lv_HYBL_i5_a /) hybi = (/ grib_in->lv_HYBL_i5_b /) P0 = grib_in->P0 lev = 1000.*(hyam+hybm) ilev = 1000.*(hyai+hybi) ; ..Make sure P0 has correct attributes P0@long_name = "reference pressure" P0@units = "Pa" ; ..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 (lev) print (ilev) print (hyam) print (hybm) print (hyai) print (hybi) print (P0) print (ntrm) print (ntrn) print (ntrk) end if ; ..Read standard scalar variables: GRIB lats run N->S QQ = grib_in->Q_GDS4_HYBL(ntSt3:ntLa3:Stride,:,:,:) QQ = QQ(:,:,::-1,:) ; make lat coords S->N QQ@_FillValue = 1.e20 printVarSummary(QQ) filp = filv(nf) ; CIQ = grib_iv->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_iv->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_iv->STL1_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL1Q = STL1Q(:,::-1,:) STL1Q@_FillValue = 1.e20 printVarSummary(STL1Q) SDQ = grib_iv->SD_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) SDQ = SDQ(:,::-1,:) SDQ@_FillValue = 1.e20 printVarSummary(SDQ) SDORQ = grib_iv->SDOR_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) SDORQ = SDORQ(:,::-1,:) SDORQ@_FillValue = 1.e20 printVarSummary(SDORQ) STL2Q = grib_iv->STL2_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL2Q = STL2Q(:,::-1,:) STL2Q@_FillValue = 1.e20 printVarSummary(STL2Q) LSMQ = grib_iv->LSM_GDS4_SFC(ntSt2:ntLa2:Stride,:,:) LSMQ = LSMQ(:,::-1,:) LSMQ@_FillValue = 1.e20 printVarSummary(LSMQ) STL3Q = grib_iv->STL3_GDS4_DBLY(ntSt2:ntLa2:Stride,:,:) STL3Q = STL3Q(:,::-1,:) STL3Q@_FillValue = 1.e20 printVarSummary(STL3Q) STL4Q = grib_iv->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") Q = g2gsh_Wrap(QQ, (/nlat,mlon/), trunc) 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(QQ) 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 Q!0 = "time" Q!1 = "lev" Q!2 = "lat" Q!3 = "lon" grib_meta(Q) ; remove/add GRIB meta-data printVarSummary(Q) 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) ; ..Read complex scalar variables; convert to target grid wcT = systemfunc("date") PS = cc2g_ntStrt_ntLast (grib_in,"LNSP_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlat,mlon/),trunc) PS = exp(PS) PHIS = cc2g_ntStrt_ntLast (grib_in,"Z_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlat,mlon/),trunc) T = cc2g_ntStrt_ntLast (grib_in,"T_GDS50_HYBL",ntSt3,ntLa3,Stride,(/nlat,mlon/),trunc) wallClockElapseTime(wcT, "T: scalar: cc2g_ntStrt_ntLast", 0) PS@_FillValue = 1.e20 PHIS@_FillValue = 1.e20 T@_FillValue = 1.e20 printVarSummary(T) grib_meta(PS) grib_meta(PHIS) 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_HYBL","VO_GDS50_HYBL"/)) wallClockElapseTime(wcUV, "UV", 0) U = UV(0,:,:,:,:) ; separate the two components V = UV(1,:,:,:,:) ; ..Meta info for each variable PS@long_name = "Surface pressure" PS@units = "Pa" 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 = "e4oper.an."+idate(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. Still "+nl+\ " other variables [eg: sea-ice cover and sea surface temperature] are on "+nl+\ " an N80 reduced grid. The expansion of this data leads to an nlat=160 / "+nl+\ " mlon=320 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+\ " (5) 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+\ " (6) 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", "ilev", "lat", "lon" /) dimSizes = (/ ntim , klev , klev+1, nlat , mlon /) dimUnlim = (/ True , False, False , False, False /) filedimdef(netcdf_out, dimNames, dimSizes, dimUnlim) filevardef (netcdf_out, "time", typeof(time), "time") filevarattdef(netcdf_out, "time", time) filevardef (netcdf_out, "hyam", typeof(hyam), "lev") filevarattdef(netcdf_out, "hyam", hyam) filevardef (netcdf_out, "hybm", typeof(hybm), "lev") filevarattdef(netcdf_out, "hybm", hybm) filevardef (netcdf_out, "lev" , typeof(lev) , "lev") filevarattdef(netcdf_out, "lev" , lev) filevardef (netcdf_out, "hyai", typeof(hyai), "ilev") filevarattdef(netcdf_out, "hyai", hyai) filevardef (netcdf_out, "hybi", typeof(hybi), "ilev") filevarattdef(netcdf_out, "hybi", hybi) filevardef (netcdf_out, "ilev", typeof(ilev), "ilev") filevarattdef(netcdf_out, "ilev", ilev) 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, "P0" , typeof(P0) , "ncl_scalar") filevarattdef(netcdf_out, "P0" , P0) filevardef (netcdf_out, "PS" , typeof(PS) , getvardims(PS)) filevarattdef(netcdf_out, "PS" , PS) filevardef (netcdf_out, "PHIS" , typeof(PHIS) , getvardims(PHIS)) filevarattdef(netcdf_out, "PHIS" , PHIS) 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) 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->hyam = (/ hyam /) netcdf_out->hybm = (/ hybm /) netcdf_out->hyai = (/ hyai /) netcdf_out->hybi = (/ hybi /) netcdf_out->lev = (/ lev /) netcdf_out->ilev = (/ ilev /) 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->P0 = (/ P0 /) netcdf_out->PS = (/ PS /) netcdf_out->PHIS = (/ PHIS /) netcdf_out->Q = (/ Q /) netcdf_out->T = (/ T /) netcdf_out->U = (/ U /) netcdf_out->V = (/ V /) 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(DATEv) delete(datev) delete(datesec) delete(hh) delete(P0) delete(PS) delete(PHIS) delete(Q) delete(T) delete(UV) delete(U) delete(V) 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.2.0.ncl") end