module lomis2 {
  use mask1 + helicon + polydiv + cgstep_mod + solver_mod + patch + mkwallwt
contains
# fill in missing data by minimizing power out of a given filter.
# by helix magic works in any number of dimensions
  subroutine lomis1( niter, xx,yy,aa, mask, doprec, npatch, nwall, nwind)
 {
    logical,                      intent( in)  :: doprec
    integer,                      intent( in)  :: niter
    type( filter), dimension (:), intent( in)  :: aa
    real,          dimension( :), intent( in)  :: mask, xx
    real,          dimension( :), intent( out) :: yy	# fitting variables
    real,          dimension( :), allocatable  :: dd, wall
    integer,       dimension( :), pointer      :: npatch, nwall, nwind
    real,          dimension( :), pointer      :: windata, winmask
    logical,       dimension( :), pointer      :: kk
    integer                                    :: nw, stat1, stat2, stat3, ip
    nw = product( nwind); allocate (windata (nw), winmask (nw), kk (nw))
    if( doprec)                          #  preconditioned
	call mask1_init (kk)
    else {
       allocate( dd( nw)); dd = 0.
    }
    call patch_init( npatch, nwall, nwind); yy = 0.
    do ip = 1, product( npatch) {
	stat1 = patch_lop (.false.,.false., xx,   windata)
        stat2 = patch_lop (.false.,.false., mask, winmask)
        kk = (winmask /= 0.)	
	if (count (kk) < nw .and. count (kk) > 0) {
	   if (doprec) {
        	call polydiv_init( nw, aa (ip))
            	call solver_prec( mask1_lop, cgstep, niter= niter, x= windata, 
                              dat= windata, 
                              prec= polydiv_lop, nprec= nw, eps= 0.)
            } else {                               #  regularized
                call helicon_init( aa (ip))
                call solver( helicon_lop, cgstep, niter= niter, x= windata, 
                             dat= dd, known = kk, x0= windata)
            }
            call cgstep_close( )
        }
        stat3 = patch_lop (.true.,.true., yy, windata)
        call patch_close ()
#	call snap ("winout.H",nwind(1), nwind(2), windata)
#	call snap ("winmask.H",nwind(1), nwind(2), winmask)
    }
    allocate (wall (size (yy))); windata = 1.
    call wallwtn(npatch, nwall, nwind, windata, wall)
    yy = yy*wall

    if (doprec)
       call polydiv_close()
    else
       deallocate (dd)
    deallocate (windata, winmask, kk, wall)  
  }
}