PRO oaflux_mask_30N30S

@common

;
;; Read oaflux grid and mask
;
dir='/Volumes/vialardj/OAFLUX/'
dir='/Users/jv/data/OAFLUX/'

yy=2004 & cy=string(yy,format='(i4.4)')
fi=dir+'lh_oaflux_'+cy+'.nc'
initncdf, fi
lh=ncdf_lec(fi,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='!'+dir+'mask_oaflux_30N30S.nc'
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