source: trunk/src/interp_erai_t2m_1989_2009.pro @ 66

;+
;
; .. _interp_erai_t2m_1989_2009.pro:
;
; ==================================================================================
; interp_erai_t2m_1989_2009.pro - Interpolation of t2 from ERA-I grid to OAFLUX grid
; ==================================================================================
;
; Interpolation of t2 from ERA-I grid to OAFLUX grid
;
; :file:`${PROJECT_ID}/20c3m_erai_t2_TROP_1989_2009.nc` containing t2 from ERA-I have been produced
; by :ref:`compute_erai_daily_region_2d.sh`.
;
; :file:`${PROJECT_ID}/mask_oaflux_30N30S.nc` containing OAFLUX grid have been produced by :ref:`oaflux_mask_30N30S.pro`.
;
; Interpolated t2 is written in
; :file:`${PROJECT_OD}/erai_t2m_19890101_20091231_oafluxgrid.nc` if this file not already exists.
;
; This output file :file:`${PROJECT_OD}/erai_t2m_19890101_20091231_oafluxgrid.nc` must be processed after by :ref:`t2m_correction_ncdf.pro`.
;
;     .. graphviz::
;
;        digraph interp_erai_t2m_1989_2009 {
;           graph [
;           rankdir="TB",
;           ]
;           file_in [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/20c3m_erai_t2_TROP_1989_2009.nc"];
;           mask [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/mask_oaflux_30N30S.nc"];
;
;           ncfile [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_t2m_19890101_20091231_oafluxgrid.nc"];
;
;           interp_erai_t2m_1989_2009 [shape=box,
;           fontname=Courier,
;           color=blue,
;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/interp_erai_t2m_1989_2009.pro"",
;           label="${PROJECT}/src/interp_erai_t2m_1989_2009.pro"];
;
;           {file_in mask} -> {interp_erai_t2m_1989_2009} -> {ncfile}
;
;          }
;
;
; SEE ALSO
; ========
;
; :ref:`interpolate_data`
;
; :ref:`project_profile.sh`
;
; :ref:`compute_erai_daily_region_2d.sh`
;
; :ref:`oaflux_mask_30N30S.pro`
;
; :func:`report <saxo:report>`
; :func:`isadirectory <saxo:isadirectory>`
; :func:`isafile <saxo:isafile>`
; :func:`initncdf <saxo:initncdf>`
; :func:`read_ncdf <saxo:read_ncdf>`
; :func:`ncdf_lec <saxo:ncdf_lec>`
; :func:`domdef <saxo:domdef>`
; :func:`call_interp2d <saxo:call_interp2d>`
; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
; :func:`ncdf_getatt <saxo:ncdf_getatt>`
;
; :ref:`t2m_correction_ncdf.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> .compile file_interp
;  IDL> interp_erai_t2m_1989_2009
;
;
; TODO
; ====
;
; strange view (lat and lon shift with ncview) : check grid init
;
; remove useless netcdf_read
;
; check time values
;
; why use :func:`call_interp2d <saxo:call_interp2d>` : this is an hidden function of SAXO
; see :func:`file_interp <saxo:file_interp>`
;
; use as input a file produced by compute_erai_daily_region_2d.sh
;
; coding rules
;
; hard coded time in module name and in output filename
;
; why two "initncdf, fullfilename_msk"
;
; hard coded attributes for t2m (missing value, short name, axis) and time (origin) : use ncdf_getatt 
;
; CF conventions
;
; KNOWN ISSUES
; ============
;
; test of existence of fullfilename_msk not very efficient because
; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
;
; EVOLUTIONS
; ==========
;
; - fplod 20110103T153734Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * computed first and last dates
;   * start using ncdf_getatt instead of hard coded attributes
;
; - fplod 20101223T130406Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * replace /Volumes/PRAVEEN/TropFlux/input_uncor/ by ${TROPFLUX_OD}
;   * complete documentation
;   * get time dimension in data file
;   * replace read_ncdf by ncdf_lec
;
;     usage of timestep keyword in read_ncdf was not efficient because of  memory problem :
;     can not read the whole time range (ok with 7580, not ok with 7591 while need is 7670
;     ::
;
;       % Unable to allocate memory: to make array.
;       Cannot allocate memory
;
;
;     When timestep keyword is not used we can see these messages::
;
;        /usr/work/incas/fplod/tropflux_d/20c3m_erai_t2_TROP_1989_2009.nc as no time axis variable
;        given by Praveen
;
;        Value of Julian date is out of allowed range
;
;        % Array used to subscript array contains out of range subscript: DATAIN
;
;
; - fplod 20101222T163216Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * replace /Volumes/PRAVEEN/ERAI_global by ${TROPFLUX_ID}
;
; - fplod 20101217T140745Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coded directory for mask_oaflux_30N30S.nc
;
; - fplod 20101215T112140Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add graph in header
;
; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
; - pbk 2008
;
;   * creation
;
;-
pro interp_erai_t2m_1989_2009
;
@cm_4cal
@cm_4data
@cm_4mesh
@cm_4data
@cm_project
;
; check for input directory
;
; test if ${PROJECT_ID} defined
CASE project_id_env OF
    ''  :  BEGIN
     msg = 'eee : ${PROJECT_ID} is not defined'
     ras = report(msg)
     STOP
           END
 ELSE: BEGIN
     msg = 'iii : ${PROJECT_ID} is ' + project_id_env
     ras = report(msg)
       END
ENDCASE
;
iodirin = isadirectory(project_id_env)
;
; existence and protection of ${PROJECT_ID}
IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
   msg = 'eee : the directory' + iodirin  + ' is not accessible.'
   ras = report(msg)
   STOP
ENDIF
;
; build mask filename
filename_msk='mask_oaflux_30N30S.nc'
;
; check if this file exists
fullfilename_msk = isafile(iodirin + filename_msk, NEW=0, /MUST_EXIST)
IF fullfilename_msk[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfilename_msk + ' was not found.'
   ras = report(msg)
   STOP
ENDIF
;
; build t2 filename
filename_t2='20c3m_erai_t2_TROP_1989_2009.nc'
;
; check if this file exists
fullfilename_t2 = isafile(iodirin + filename_t2, NEW=0, /MUST_EXIST)
IF fullfilename_t2[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfilename_t2 + ' was not found.'
   ras = report(msg)
   STOP
ENDIF
;
; test if ${PROJECT_OD} defined
CASE project_od_env OF
  '' : BEGIN
         msg = 'eee : ${PROJECT_OD} is not defined'
         ras = report(msg)
       STOP
       END
  ELSE: BEGIN
          msg = 'iii : ${PROJECT_OD} is ' + project_od_env
          ras = report(msg)
        END
 ENDCASE
;
; check if output data will be possible
iodirout = isadirectory(project_od_env)
;
; existence and protection
IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
    msg = 'eee : the directory' + iodirout  + ' was not found.'
    ras = report(msg)
    STOP
ENDIF
;
; build output filename
filename_out = 'erai_t2m_19890101_20091231_oafluxgrid.nc'
fullfilename_out = iodirout + filename_out
; in order to avoid unexpected overwritten
IF (FILE_TEST(fullfilename_out) EQ 1) THEN BEGIN
   msg = 'eee : the file ' + fullfilename_out  + ' already exists.'
   ras = report(msg)
   STOP
ENDIF
;
; define grid parameters with t2 file
initncdf, fullfilename_t2
domdef
latin=reform(gphit(0,*)) 
lonin=reform(glamt(*,0))
print, 'lat grid from data',min(latin),max(latin),latin(1)-latin(0)
print, 'lon grid from data',min(lonin),max(lonin),lonin(1)-lonin(0)
;
; get time in t2 file
timein=ncdf_lec(fullfilename_t2,var='time')
jptin=n_elements(timein)
print, 'time steps from data ', jptin
print, 'The first 10 time values (variable timein) = ', timein[0:9]
;
; find first and last dates yyyymmdd
; they will be written in global attributes of output file
da=jul2date(julday(01, 01, 1957,timein[0]))
cda0=string(da,format='(i8.8)')
da=jul2date(julday(01, 01, 1957,timein[jptin-1]))
cda1=string(da,format='(i8.8)')
print, 'first date ', cda0
print, 'last date ' , cda1

; read t2 data
;++ pb memory t2min=read_ncdf("t2",0,jptin-1,/timestep,file=fullfilename_t2,/nostr)
; the following line is here just to prepare replacement of read_ncdf by ncdf_lec
t2min_read_ncdf=read_ncdf("t2",0,10,/timestep,file=fullfilename_t2,/nostr)
help, t2min_read_ncdf
;++print, 'The first 10 time values (variable time) = ', time[0:9]
;++print, 'time steps after read_ncdf (variable jpt) ', jpt
t2min=ncdf_lec(fullfilename_t2,var='t2')
;
mskin=glamt*0.+1.

; define grid parameters with mask file
initncdf, fullfilename_msk
domdef
latout=reform(gphit(0,*)) 
lonout=reform(glamt(*,0))
print, 'lat grid from mask ',min(latout),max(latout),latout(1)-latout(0)
print, 'lon grid from mask ',min(lonout),max(lonout),lonout(1)-lonout(0)
mskout=read_ncdf("msk", file=fullfilename_msk,/nostr)
;
help, t2min,lonin,latin,mskin,lonout,latout,mskout
;
; interpolation and mask
t2mout=fltarr(jpi,jpj,jptin)
for jt=0,jptin-1 do begin
  tab=reform(t2min(*,*,jt))
  t2mout(*,*,jt)=call_interp2d(tab,lonin,latin,mskin,lonout,latout,method='bilinear')
  t2mout(*,*,jt)=t2mout(*,*,jt)*mskout+(1.-mskout)*1.e20
endfor

initncdf, fullfilename_msk

lat=latout
lon=lonout
ncfile='!' + fullfilename_out
print, 'ncfile ', ncfile
ncdf_getatt, fullfilename_msk, 'longitude', units=units
ncdf_getatt, fullfilename_msk, 'longitude', long_name=long_name
lon_attr={units:units, long_name:long_name}
ncdf_getatt, fullfilename_msk, 'latitude', units=units
ncdf_getatt, fullfilename_msk, 'latitude', long_name=long_name
lat_attr={units:units, long_name:long_name}
ncdf_getatt, fullfilename_t2, 'time', units=units
ncdf_getatt, fullfilename_t2, 'time', long_name=long_name
time_attr={units:units, long_name:long_name, time_origin:' 1957-JAN-01 00:00:00'}
ncdf_getatt, fullfilename_t2, 't2', units=units
ncdf_getatt, fullfilename_t2, 't2', long_name=long_name
t2m_attr={units:units, long_name:long_name, missing_value:1e20, short_name:'t2m',axis:'TYX'}
globattr={source:'Data are from ECMWF ERA-Interim reanalysis', timerange:cda0+' - '+cda1}


help, t2mout
help, timein
;ncfields = 't2m[longitude,latitude,*time]=t2mout:t2m_attr; ' $
;                      + 'longitude[]=lon:lon_attr; ' $
;                      + 'latitude[]=lat:lat_attr; ' $
;                     + 'time[]=timein:time_attr ' $
;                      + ' @globattr'

; ok ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr'
; ok ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timein:time_attr'
; ok ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[*]=timein:time_attr'
; ok ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[*]=timein:time_attr@globattr'
; pas ok ncfields = 't2m[longitude,latitude,*time]=t2mout:t2m_attr;longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[*]=timein:time_attr@globattr'
; pas ok  ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[*]=timein:time_attr;t2m[longitude,latitude,*time]=t2mout:t2m_attr@globattr'

; the simpliest one
nlat=73
nlon=96
latitudetp=findgen(nlat)
longitudetp=findgen(nlon)
timetp=findgen(2)
pressure3d=findgen(nlon,nlat,2)
patts = {units:'pascal',missing_value:1e10}
globatts={source:'my program',version:2}
NCFIELDS = 'pressure3d[longitudetp,latitudetp,*timetp]=pressure3d:patts;longitudetp[];latitudetp[];timetp[]@globatts'
ncfile='my4.nc'
@ncdf_quickwrite

timepb=findgen(2)
t2m=findgen(350,60,2)
help, lat
help, lon
help, t2m 
help, time
help, timein
help, timepb
ncfile='!' + fullfilename_out
; pas ok ncfields = 'longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timein:time_attr;t2m[longitude,latitude,*time]'
; pas ok ncfields = 't2m[longitude,latitude,*time];longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timein:time_attr'
; ok avec timepb=findgen(2) ncfields = 't2m[longitude,latitude,*time];longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timepb:time_attr'
; ok
;timepb=findgen(2)
;timepb[0]=280512.
;timepb[1]=280536.
;ncfields = 't2m[longitude,latitude,*time];longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timepb:time_attr'

timepb=findgen(7670)
ncfields = 't2m[longitude,latitude,*time];longitude[]=lon:lon_attr;latitude[]=lat:lat_attr;time[]=timepb:time_attr'
@ncdf_quickwrite


end