;+
;
; ======================
; oaflux_mask_30n30s.pro
; ======================
;
; .. function:: oaflux_mask_30n30s()
;
; DESCRIPTION
; ===========
;
; Produce a NetCDF file with OAFLUX mask over global tropical oceans (30N-30S)
;
; .. note::
;
;    Longitude range [30.5,379.5] was choosen in order not to make any cut in
;    the Atlantic Ocean.
;    Otherwise if we choose longitude range to [0.5,359.5], Atlantic Ocean
;    will be split to two and will appear on right and left sides of the map.
;
; .. graphviz::
;
;    digraph oaflux_mask_30n30s {
;
;       file_oaflux [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/lh_oaflux_2004.nc"];
;       mask [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/mask_oaflux_30N30S.nc"];
;
;       oaflux_mask_30n30s [shape=box,
;       fontname=Courier,
;       color=blue,
;       URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/oaflux_mask_30n30s.pro",
;       label="${PROJECT}/src/oaflux_mask_30n30s.pro"];
;
;       {file_oaflux} -> {oaflux_mask_30n30s} -> {mask}
;
;      }
;
; :returns: 0 if ok, -1 if error
; :rtype: integer
;
; SEE ALSO
; ========
;
; :ref:`guide Gather data <gather_data>`
;
; :ref:`guide data OAFLUX <data_in_oaflux>`
;
; Environement :
; :ref:`project_profile.sh`,
; :ref:`project_init.pro`,
; :ref:`cm_project.pro`
;
; Previous step : :ref:`get_oaflux.sh`
;
; Use :
; :func:`report <saxo:report>`,
; :func:`initncdf <saxo:initncdf>`,
; :func:`ncdf_lec <saxo:ncdf_lec>`,
; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`,
; :func:`isadirectory <saxo:isadirectory>`,
; :func:`isafile <saxo:isafile>`
;
; Following program use :func:`oaflux_mask_30n30s` output:
; :ref:`TropFlux_19890101_20091231.pro`,
; :func:`interp_erai_dewt`,
; :func:`interp_erai_lwr`,
; :ref:`interp_erai_sst_1989_2009.pro`,
; :func:`interp_erai_t2m`,
; :ref:`interp_erai_ws_1989_2009.pro`,
; :ref:`interp_olr_30n30s_1989_2009.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> result = oaflux_mask_30n30s()
;
; TODO
; ====
;
; SAXO new (keep compatibility false)
;
; coding rules
;
; NetCDF CF (may be area_type)
;
; global attributes corrections :
;
;  - written now : OAFLUX mask over the Indian Ocean: 40E-120E, 30S-30N
;  - while data latitude=[-29.5,+29.5] longitude=[30.5,379.5]
;
; 30n30s in name vs future global output ... parameters lat and long max
;
; KNOWN ISSUES
; ============
;
; test of existence of fullfilename_in not very efficient because
; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20120319
;
;   * taking project_overwite into account
;
; - fplod 20120306
;
;   * try to add compile_opt
;
; - fplod 20120305
;
;   * pro -> function
;
; - fplod 20110411T140133Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * usage of tropflux_init and cm_project
;   * replace 30N30S by 30n30s
;
; - fplod 20101217T081915Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coded directory - usage of ${PROJECT_ID} and ${PROJECT_OD}
;   * add IO test
;
; - fplod 20101216T141137Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
; - jv 2010
;
;   * creation
;
;-
;
FUNCTION oaflux_mask_30n30s
;
;++ compile_opt idl2, strictarrsubs, logical_predicate
;
@cm_4cal
@cm_4data
@cm_4mesh
@cm_4data
@cm_project

; Return to caller if errors
on_error, 2
;
result = -1
;
usage = 'result=oaflux_mask_30n30s()'
nparam = N_PARAMS()
IF (nparam NE 0) THEN BEGIN
   ras = report(['Incorrect number of arguments.' $
         + '!C' $
         + 'Usage : ' + usage])
   return, result
ENDIF
;
; check for input file
; build input filename
yyyy=2004
cy=string(yyyy,format='(i4.4)')
filename_in='lh_oaflux_'+cy+'.nc'
;
; check if this file exists
fullfilename_in = isafile(project_id_env + filename_in, NEW=0, /MUST_EXIST)
IF fullfilename_in[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfilename_in + ' was not found.'
   ras = report(msg)
   return, result
ENDIF
;
fullfilename_out=project_od_env+'mask_oaflux_30N30S.nc'
; in order to avoid unexpected overwritten
PRINT, 'project_overwrite : ',project_overwrite
IF ((FILE_TEST(fullfilename_out) EQ 1) AND (project_overwrite EQ 0)) THEN BEGIN
   msg = 'eee : the file ' + fullfilename_out  + ' already exists.'
   ras = report(msg)
   return, result
ENDIF
;
; Read oaflux grid and mask
initncdf, fullfilename_in
lh=ncdf_lec(fullfilename_in,var='lhtfl',count=[jpi,jpj,1])
lh=float(lh)
valmask=1.e20
ind=where(lh ge 32000,compl=nind)
lh[ind]=valmask
lh[nind]=-0.1*lh[nind]
mask_out=(lh ne valmask)
domdef, 0.,360.,-30.,30.

lon=reform(glamt[*,0])
ind1=where(lon le 20.)
ind2=where(lon ge 30.)

lat_out=reform(gphit[0,firstyt:lastyt])

lon_out=[reform(glamt[ind2,0]),reform(glamt[ind1,0])+360]

; ++ this line is incompatible with 
; compile_opt idl2, strictarrsubs, logical_predicate
mask_out=mask_out([ind2,ind1],firstyt:lastyt)*1.

ncfile='!'+fullfilename_out
lon_attr={units:'degrees_east',long_name:'Longitude'}
lat_attr={units:'degrees_north',long_name:'Latitude'}
msk_attr={long_name:'land sea mask',short_name:'msk',axis:'YX'}
globattr={description:'OAFLUX mask over the Indian Ocean: 40E-120E, 30S-30N'}

ncfields = 'msk[longitude,latitude]=mask_out:msk_attr; ' $
                      + 'longitude[]=lon_out:lon_attr; ' $
                      + 'latitude[]=lat_out:lat_attr; ' $
                      + ' @ globattr'

@ncdf_quickwrite
;
result = 0
return, result
;
END