anchor.diff:            end do
           end do
   
  -        do j=2,jmm1
  -          do i=2,imm1
  -            zflux(i,j,1)=0.d0
  -            if(itera.eq.1) zflux(i,j,1)=w(i,j,1)*f(i,j,1)*art(i,j)
  -            zflux(i,j,kb)=0.d0
  -          end do
  -        end do
  +        zflux(2:imm1,2:jmm1,1) = 0.
  +        if(itera.eq.1) then
  +          zflux(2:imm1,2:jmm1,1) = w(2:imm1,2:jmm1,1)*f(2:imm1,2:jmm1,1)
  +     $                            *art(2:imm1,2:jmm1)
  +        end if
  +        zflux(2:imm1,2:jmm1,kb) = 0.
   
           do k=2,kbm1
             do j=2,jmm1
  
positive.code:         end do
        end do
  
        do k=1,kb
          do j=1,jm
            do i=1,im
              a(i,j,k)=0.d0
              c(i,j,k)=0.d0
              ee(i,j,k)=0.d0
              gg(i,j,k)=0.d0
            end do
          end do
        end do
  
        do k=2,kbm1
          do j=1,jm
  
negative.code: 
  ! initialize fields that are not calculated on all boundaries
  ! but are later used there
        l0    = 0.
        boygr = 0.
        prod  = 0.
  
        do j=1,jmm1
          do i=1,imm1
  