;+
; NAME:
;      NEI2005_MASK
;
;
; PURPOSE: To build the mask used for EPA 2005 emissions, with options
;           to keep the Mexican emission and/or the Canadian regions
;
; KEYWORDS:
; 
;       /SAVEIT: set to save the mask in a bpch file
;       
;       MASKFILE: name of the file where the mask is saved (default
;                  should be fine)
;       
;       /CANADA : to zero the mask where we have canadian CAC data
;
;       /MEXICO : to zero the mask where we have mexican BRAVO emissions
;       
;       /FT : fine tune when doing /CANADA (to add a couple of boxes
;              to the mnask)
;
;       
;
; CATEGORY:
;       Mask, emissions
;
; CALLING SEQUENCE:
;
;  !! NOTE : pass /USE_SAVED after the first call, when calling
;  !!        the routine several times in a row, so you do not 
;  !!        have to recompute the regridding weights each time.
;
;
;  ; To see what's up:
;       NEI2005_MASK, limit=[18., -140., 60., -55], gylab=0, gxlab=0, /isotro, rectangle=0, min_valid=0.5
;       NEI2005_MASK, limit=[18., -140., 60., -55], gylab=0, gxlab=0, /isotro, rectangle=0, min_valid=0.5,/use_saved,/canada
;       NEI2005_MASK, limit=[18., -140., 60., -55], gylab=0, gxlab=0, /isotro, rectangle=0, min_valid=0.5,/use_saved,/canada,/ft
;
;  ; to save:
;      NEI2005_MASK, /saveit,/use_saved
;      NEI2005_MASK, /saveit,/use_saved,/ft,/mex
;      NEI2005_MASK, /saveit,/use_saved,/ft,/can
;      NEI2005_MASK, /saveit,/use_saved,/ft,/mex,/can
;
;  ; to check
;     test_nie05mask
;
;
; MODIFICATION HISTORY:
;
;       Mon Oct 26 11:56:41 2009, Philippe Le Sager
;       <phs@sol.as.harvard.edu>
;
;		initial version
;
;-


pro test_nie05mask

   f = ['usa.mex.mask.nei2005.geos.1x1', 'usa.can.mask.nei2005.geos.1x1', $
        'usa.can.mex.mask.nei2005.geos.1x1', 'usa.mask.nei2005.geos.1x1']


   multipanel, rows=2, cols=2
   
   for k=0, 3 do begin

      ctm_plot, filename=f[k], 'LANDMAP',  Tracer=2, tau0=0., $
                limit=[18., -140., 60., -55], dlon=20, dlat=20., /sample, $
                title=StrUpCase(neim_str(f[k], '.mask.nei2005.geos.1x1')),  $
                /nocbar, /countries,  cbunit=''

   endfor

   multipanel,  /off

   return
end

function neim_str, str, regex

   ; little function to take off a substring

   p = stregex(str, regex, len=l)
   
   return, p eq -1 ? str : strmid(str,0,p) + strmid(str,p+l,strlen(str)-1)

end




pro NEI2005_MASK, _extra=e, saveit=saveit, maskfile=maskfile, canada=canada, $
                  mexico=mexico, ft=ft


   ; ========== FILES

   cacF = '/as/data/geos/GEOS_1x1/CAC_200801/CanadaMask.geos.1x1'

   bravoF = '/as/data/geos/GEOS_1x1/BRAVO_200607/BRAVO.MexicoMask.generic.1x1'

   nei2005F = '/home/phs/data/NEI2005/data/version2/NEI2005.GEOS5.1x1.AVG.bpch'


   ; ========== DATA

   ctm_get_data, cac,     'landmap', file=cacF,     tracer=2
   ctm_get_data, nei2005,            file=nei2005F, tracer=1

   ctm_get_data, bravo,   'landmap', file=bravoF,   tracer=2


   ; ========== GRIDS

   ; geos 1x1
   getModelandGridinfo, cac, model, grid, lon=lon, lat=lat

   ; generic 1x1
   getModelandGridinfo, bravo, model2, grid2

   
   ; ========== REGRID generic to GC

   bravo2   = CTM_RegridH( *bravo.data, Grid2, Grid, _extra=e)


   
   ; ========== DEFINE Mask

   ; where we have data
   newmask = (*nei2005[0].data gt 0)[*, *, 0]
   
   ; to zero where we have CAC
   if keyword_set(canada) then begin

      keeper = where(*cac.data ne 0.)
   
      k2d = array_indices(*cac.data, keeper)

      ; just add 1 to latitude index to keep the border
      newmask[ k2d[0, *],  k2d[1, *]+1 ] = 0

      ; fine tuning: add a couple of boxes
      if keyword_set(ft) then newmask[97, 133:134] = 1
      
   endif

   ; to zero where we have BRAVO
   if keyword_set(mexico) then begin

      keeper = where(bravo2 ne 0.)
   
      k2d = array_indices(bravo2, keeper)
      
      ; just substract 1 to latitude index to keep the border with NEI 05
      newmask[ k2d[0, *],  k2d[1, *]-2 ] = 0

      ; fine tuning: add a couple of boxes
      if keyword_set(ft) then newmask[[80, 81],[118, 117]] = 1
      
   endif



   ; ======= SAVE MASK on the nested grid
   if keyword_set(saveit) then begin


      Success = CTM_Make_DataInfo( float(newmask),   $
                                   ThisDataInfo,     ThisFileInfo,   $
                                   ModelInfo=model,  GridInfo=Grid,  $
                                   DiagN='LANDMAP',  Tracer=2,       $
                                   Tau0=0d0,         Tau1=0d0,       $
                                   Unit='unitless',                  $
                                   /No_Global )

      ; output filename if you save the new mask
      midname = 'usa.can.mex.'
      if keyword_set(canada) then midname = neim_str(midname, 'can.')
      if keyword_set(mexico) then midname = neim_str(midname, 'mex.')

      if n_elements(maskfile) eq 0 then $
         maskfile = midname+'mask.nei2005.geos.1x1'


      CTM_WriteBpch, ThisDataInfo, ThisFileInfo, FileName=maskfile

      ; clean No_Global pointers (delete the copy of data made in
      ; ctm_make_datainfo, and still alive because pointed to)
      ptr_free, ThisDataInfo.data
      ptr_free, Thisfileinfo.gridinfo

   endif




   ; ==============  PLOT

   myct, 120 ; load 3-colors accent ct
   
   tvmap, newmask, lon, lat, SAMPLE=1, _extra=e, $
          /contin, /usa, title='EPA mask & data',  $
             
          ; to use TVMAP features (call to MAP_GRID and MAP_LABELS):
          /grid, $
          lats=grid.yedge, lons=grid.xedge


   return
end