;+
;
; .. _oaflux_mask_30N30S.pro:
;
; ======================
; oaflux_mask_30N30S.pro
; ======================
;
; Produce a NetCDF file with  OAFLUX mask over global tropical oceans 
;
;     .. graphviz::
;
;        digraph oaflux_mask_30N30S {
;           graph [
;           rankdir="LR",
;           ]
;           file_oaflux [shape=ellipse,fontname=Courier,label="${TROPFLUX_ID}/lh_oaflux_2004.nc"];
;           mask [shape=ellipse,fontname=Courier,label="${TROPFLUX_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="${TROPFLUX}/src/oaflux_mask_30N30S.pro"];
;
;           {file_oaflux} -> {oaflux_mask_30N30S} -> {mask}
;
;          }
;
; SEE ALSO
; ========
;
; :ref:`guide data OAFLUX <data_in_oaflux>`
;
; :ref:`tropflux_profile.sh`
;
; :func:`report`
; :func:`initncdf`
; :func:`ncdf_lec`
; :ref:`ncdf_quickwrite`
; :func:`isadirectory`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> oaflux_mask_30N30S
;
; TODO
; ====
;
; SAXO new (keep compatibility false)
;
; coding rules
;
; pro -> function
;
; 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]
;
; KNOWN ISSUES
; ============
;
; test of existence of fullfilename_in not very efficient because
; MUST_EXIST keyword of :func:`isafile` not yet implemented
;
; EVOLUTIONS
; ==========
;
; - fplod 20101217T081915Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coded directory - usage of ${TROPFLUX_ID} and ${TROPFLUX_OD}
;   * add IO test
;
; - fplod 20101216T141137Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
; - jv 2010
;
;   * creation
;
;-
;
PRO oaflux_mask_30N30S
;
@common
;
; check for input file
;
; test if ${TROPFLUX_ID} defined
tropflux_id_env=GETENV('TROPFLUX_ID')
CASE tropflux_id_env OF
    ''  :  BEGIN
     msg = 'eee : ${TROPFLUX_ID} is not defined'
     ras = report(msg)
     STOP
           END
 ELSE: BEGIN
     msg = 'iii : ${TROPFLUX_ID} is ' + tropflux_id_env
     ras = report(msg)
       END
ENDCASE
;
iodirin = isadirectory(tropflux_id_env)
;
; existence and protection of ${TROPFLUX_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 input filename
yyyy=2004
cy=string(yyyy,format='(i4.4)')
filename_in='lh_oaflux_'+cy+'.nc'
;
; check if this file exists
fullfilename_in = isafile(iodirin + 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)
   STOP
ENDIF
;
; test if ${TROPFLUX_OD} defined
tropflux_od_env=GETENV('TROPFLUX_OD')
CASE tropflux_od_env OF
  '' : BEGIN
   msg = 'eee : ${TROPFLUX_OD} is not defined'
   ras = report(msg)
   STOP
       END
ELSE: BEGIN
   msg = 'iii : ${TROPFLUX_OD} is ' + tropflux_od_env
   ras = report(msg)
        END
ENDCASE
;
; check if output data will be possible
iodirout = isadirectory(tropflux_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
;
fullfilename_out=iodirout+'mask_oaflux_30N30S.nc'
; 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
;
; 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]

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

END