;;	IDL script to perform remapping using SCRIP-computed
;;	interpolation weights.
;;
;;	For non-conservative remappings.
;;
;;	(template -- adjust as necessary)
;;

dirin = '/fis/cgd/cseg/csm/inputdata/ocn/pop/gx1v5/forcing/'
datain = dirin+'sfwf_mm_PHC2_salx_flxio_20061230.ieeer8

dataout = '/ptmp/yeager/sfwf_tx1v1.ieeer8'
ncout = '/ptmp/yeager/sfwf_tx1v1.nc'
mapfile = '/ccsm/ocn/yeager/scrip/mappings/map_gx1v5_to_tx1v1_conserv_da.nc'

nx_i = 320
ny_i = 384
nx_o = 360
ny_o = 240

datai = dblarr(nx_i,ny_i,12)	
datao = dblarr(nx_o,ny_o,12)

openr,1,datain,/SWAP_END
readu,1,datai
readu,1,datai
close,1

ncid = ncdf_open(mapfile)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Get Source grid information
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
varid = ncdf_varid(ncid,"src_grid_dims")
ncdf_varget,ncid,varid,src_grid_dims
nxs = src_grid_dims[0]
nys = src_grid_dims[1]
dimid = ncdf_dimid(ncid,'n_a')
ncdf_diminq,ncid,dimid,name,grdsz1
varid = ncdf_varid(ncid,'col')
ncdf_varget,ncid,varid,src_address

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Get Destination grid information
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
varid = ncdf_varid(ncid,"dst_grid_dims")
ncdf_varget,ncid,varid,dst_grid_dims
nxd = dst_grid_dims[0]
nyd = dst_grid_dims[1]
dimid = ncdf_dimid(ncid,'n_b')
ncdf_diminq,ncid,dimid,name,grdsz2
varid = ncdf_varid(ncid,'row')
ncdf_varget,ncid,varid,dst_address
varid = ncdf_varid(ncid,'xc_b')
ncdf_varget,ncid,varid,dst_x
varid = ncdf_varid(ncid,'yc_b')
ncdf_varget,ncid,varid,dst_y
dst_x = reform(dst_x,nxd,nyd)
dst_y = reform(dst_y,nxd,nyd)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Get Mapping information
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
dimid = ncdf_dimid(ncid,'n_s')
ncdf_diminq,ncid,dimid,name,nlnk
varid = ncdf_varid(ncid,'S')
ncdf_varget,ncid,varid,remap_matrix
ncdf_close,ncid

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Compute Destination values
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

for im=0,11 do begin

missval = -99.0
dstd = fltarr(nxd,nyd)
tmp1 = reform(datai[*,*,im],grdsz1)
;tmp2 = replicate(missval,grdsz2)
tmp2 = replicate(0.,grdsz2)
for in = 0L,nlnk-1 do begin
    if (src_address(in) NE 0 AND $
      tmp1(src_address(in)-1) NE missval) then begin

           if (tmp2(dst_address(in)-1) EQ missval) then begin
                tmp2(dst_address(in)-1) = $
                tmp1(src_address(in)-1)*remap_matrix(in)
           endif else begin
                tmp2(dst_address(in)-1) = tmp2(dst_address(in)-1) +  $
                tmp1(src_address(in)-1)*remap_matrix(in)
           endelse

    endif
endfor
dstd = reform(tmp2,nxd,nyd)
datao[*,*,im] = dstd
endfor

ncid = ncdf_create(ncout,/clobber)
dim1 = ncdf_dimdef(ncid,'lat',nyd)
dim2 = ncdf_dimdef(ncid,'lon',nxd)
dim3 = ncdf_dimdef(ncid,'month',12)
var1 = ncdf_vardef(ncid,'flxio',[dim2,dim1,dim3],/double)
ncdf_attput,ncid,var1,'missing_value',missval
ncdf_control,ncid,/ENDEF
ncdf_varput,ncid,var1,datao
ncdf_close,ncid

openw,1,dataout,/SWAP_END
writeu,1,datao
close,1

end