pro ssmem_get_input, vsl, cmap, cln, uv, vis, wgt, nmap, misc, $
 model = mmap, flux = flux, ferr = ferr, spwgt = spwgt, tmin = tmin

    mdldir = '../Model/'
    inp ={inuvfile:mdldir + 'a3.min', $
     tbegstr:'180000',tendstr:'220000', $
     bif:0, eif:39, $
     fif:[-1], $    ;tab added for SSMEM
;     fif:[      1,      3,      5,  6,  7,      9, $
;           10, 11,     13, 14, 15,     17, 18, 19, $
;               21, 22, 23,     25, 26, 27,     29, $
;           30, 31,     33, 34, 35,     37, 38, 39], $
     poln:'R', imsiz:256, $
     niter:2000, clgain:0.02, clstop:2., $
     method:'CLEAN + SSMEM', outfile:'tmp.mou'}
    spwgt = 1.
    cln_restore = mdldir + 'a3_4h2m_cr.mou'   ;set to restore CLEAN map
;    cln_save = mdldir + 'a1_2m_cr.mou'   ;set to save CLEAN map
    trufile = mdldir + 'a3_tru.mou'  ;set to retrieve true map flux
;    model = 'dg_1r.mou' ;set to retrieve model map
;    conv = 1b   ;set to convolve the model with beam.

    ssmem_vsel, inp, uv, vis, wgt, tp, vsl
    if vsl.chk eq 0 then return

    m = inp.imsiz
    m2 = m / 2
    f = vsl.f_all   ;frequency array
    nf = n_elements(f)
    nf1 = nf - 1l

    if keyword_set(cln_restore) then begin
        mou = readmou(cln_restore)
        f_mou = mou.f_ghz
        flag = bytarr(n_elements(f_mou))
        for k = 0l, nf1 do begin
            match = where(abs(f_mou - f[k]) lt 0.1, nmatch)
            if nmatch eq 1l then flag[match] = 1b
        endfor
        ksl = where(flag, nsl)
        sz = size(mou.tb_xy)
        if nsl eq nf and sz[1] eq m then begin
            cmap = mou.tb_xy[*, *, ksl]
            xyint = mou.xyint[ksl]
            cln = {xyint : xyint[nf1], bmin : mou.bmin * xyint[nf1], $
             bmaj : mou.bmaj * xyint[nf1], pa : mou.pa, $
             box : [-m2, -m2, m2 - 1, m2 - 1], f_ghz : f_mou[ksl]}
        endif else begin
            vsl.chk = 0
            print, 'Wrong CLEAN map.'
            return
        endelse
    endif else begin
        ssmem_clean, uv, vis, inp, vsl, cmap, cln;, flux = flux;, ferr = ferr

        xyint = cln.xyint * f[nf1] / f
        cbm = [cln.bmin / xyint[nf1], cln.bmaj / xyint[nf1], cln.pa]
        nul = fltarr(nf)

        if keyword_set(cln_save) then savemou, cln_save, f, cmap, xyint, cbm, nul, nul
    endelse

    TB2flux = 3.6E-5 * (xyint * f)^2

    if keyword_set(trufile) then begin
        mou = readmou(trufile)
        f_mou = mou.f_ghz
        flag = bytarr(n_elements(f_mou))
        for k = 0l, nf1 do begin
            match = where(abs(f_mou - f[k]) lt 0.1, nmatch)
            if nmatch eq 1l then flag[match] = 1b
        endfor
        ksl = where(flag, nsl)
        if nsl eq nf then begin
            f_mou = f_mou[ksl]
            xyint_mou = mou.xyint[ksl]
            tb2flux_mou = 3.6e-5 * (xyint_mou[0] * f_mou[0])^2
            sz = size(mou.tb_xy)
            m_mou = sz[1]
            m2_mou = m_mou / 2l
            x0 = round((-m2 * xyint[0] / xyint_mou[0]) + m2_mou) > 0l
            x1 = round(((m2 - 1l) * xyint[0] / xyint_mou[0]) + m2_mou) < (m_mou - 1l)
            flux = total(total(mou.tb_xy[x0 : x1, x0 : x1, ksl], 1), 1) * tb2flux_mou
        endif else begin
            vsl.chk = 0
            print, 'Wrong true map.'
            return
        endelse
    endif else begin
;        sz = size(tp)
;        flux = total(total(tp[*, [0, 1], *], 1), 1) / (sz[1] * 2l)

        flux = total(total((cmap + 0.01) > 0., 1), 1) * TB2flux
    endelse

    if keyword_set(model) then begin
        mou = readmou(model)
        f_mou = mou.f_ghz
        flag = bytarr(n_elements(f_mou))
        for k = 0l, nf1 do begin
            match = where(abs(f_mou - f[k]) lt 0.1, nmatch)
            if nmatch eq 1l then flag[match] = 1b
        endfor
        ksl = where(flag, nsl)
        if nsl eq nf then begin
            f_mou = f_mou[ksl]
            xyint_mou = mou.xyint[ksl]
            tb2flux_mou = 3.6e-5 * (xyint_mou[0] * f_mou[0])^2
            sz = size(mou.tb_xy)
            m_mou = sz[1]
            m2_mou = m_mou / 2l
            x = (findgen(m) - m2) * (xyint[0] / xyint_mou[0]) + m2_mou
            mmap = fltarr(m, m, nf)
            mcb2 = long(mou.bmaj * 4.) < (m2_mou - 1l)
            mcb = 2l * mcb2 + 1l
            xcb = findgen(mcb) - mcb2
            one = replicate(1., mcb)
            ycb = one # xcb
            xcb = xcb # one
            ker = cvd_ellipse(xcb, ycb, [1., 0., 0., mou.bmin, mou.bmaj, -mou.pa, 0.]) $
             / (!pi * mou.bmin * mou.bmaj)
            for k = 0l, nf1 do begin
                if conv then mmapk = convol(mou.tb_xy[*, *, ksl[k]], ker, /edge_wrap) $
                 else mmapk = mou.tb_xy[*, *, ksl[k]]
                mmap[*, *, k] = interpolate(mmapk, x, x, /grid)
            endfor
        endif else begin
            vsl.chk = 0
            print, 'Wrong model map.'
            return
        endelse
    endif else begin
        mmap = fltarr(m, m, nf)
        for k = 0l, nf1 do mmap[*, *, k] = flux[k] / tb2flux[k] / float(m)^2
    endelse
end