!    do k=1,nwrk
!        do i=2,nz
!            Z=0;X=0;T=0
!            alfa=velw(i,k)*dt
!!
!! 8th order extrapolator in space: Fornberg, B - Mathematics of Computation, Vol. 51, No. 184. (Oct., 1988)
!!
!            if (k>4.and.k<(nwrk-3)) then
!                X=kx2*velw(i,k)**2*(-0.00178571428571*(u(i,k+4,it)+u(i,k-4,it))+ & 
!                                0.02539682539683*(u(i,k+3,it)+u(i,k-3,it))- &
!                                0.20000000000000*(u(i,k+2,it)+u(i,k-2,it))+ &
!                                1.60000000000000*(u(i,k+1,it)+u(i,k-1,it))- &
!                                2.84722222222222*(u(i,k,it)))
!            end if

!            if (i==2) then
!                Z=kz2*velw(i,k)**2*(u(i+1,k,it)+u(i-1,k,it)-2*u(i,k,it))
!            else if (i==3.or.i==4) then
!                Z=kz2/12*velw(i,k)**2*(-u(i+2,k,it)-u(i-2,k,it)+16*u(i+1,k,it)+16*u(i-1,k,it)-30*u(i,k,it))
!            else if (i>4.and.i<(nz-3)) then
!                Z=kz2*velw(i,k)**2*(-0.00178571428571*(u(i+4,k,it)+u(i-4,k,it))+ & 
!                                0.02539682539683*(u(i+3,k,it)+u(i-3,k,it))- &
!                                0.20000000000000*(u(i+2,k,it)+u(i-2,k,it))+ &
!                                1.60000000000000*(u(i+1,k,it)+u(i-1,k,it))- &
!                                2.84722222222222*(u(i,k,it)))
!            end if

!            T=2*u(i,k,it)-u(i,k,it-1)

!            u(i,k,it+1)=(Z+X+T)!*taper(i,k)

!!
!! Boundary conditions: Reynolds, A. - GEOPHYSICS. VOL. 43, NO. 6 (OCTOBER 1978)
!!
!            if (k<=4) then
!               u(i,k,it+1)=u(i,k,it)+u(i,k+1,it)-u(i,k+1,it-1)+ &
!                   alfa/dwrk*(u(i,k+1,it)-u(i,k,it)-(u(i,k+2,it-1)-u(i,k+1,it-1)))
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            else if (k>=nwrk-3) then
!               u(i,k,it+1)=u(i,k,it)+u(i,k-1,it)-u(i,k-1,it-1)- &
!                   alfa/dwrk*(u(i,k,it)-u(i,k-1,it)-(u(i,k-1,it-1)-u(i,k-2,it-1)))
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            end if
            
!            if (i>=nz-3) then
!               u(i,k,it+1)=u(i,k,it)+u(i-1,k,it)-u(i-1,k,it-1)- &
!                   alfa/dz*(u(i,k,it)-u(i-1,k,it)-(u(i-1,k,it-1)-u(i-2,k,it-1)))   
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            end if

!            if (it==nt) then
!                T=u(i,k,it)
!            else
!                T=u(i,k,it+1)-2*u(i,k,it)
!            end if
            
!            if (k>4 .and. k<(nwrk-3)) then
!                X=dwrk**(-2)*(-0.00178571428571*(u(i,k+4,it)+u(i,k-4,it))+ & 
!                            0.02539682539683*(u(i,k+3,it)+u(i,k-3,it))- &
!                            0.20000000000000*(u(i,k+2,it)+u(i,k-2,it))+ &
!                            1.60000000000000*(u(i,k+1,it)+u(i,k-1,it))- &
!                            2.84722222222222*(u(i,k,it)))
!            end if
                        
!            if (i==1) then
!                Z=dz**(-2)*(u(i+1,k,it)-2*u(i,k,it))
!            else if (i==2) then
!                Z=dz**(-2)*(u(i+1,k,it)+u(i-1,k,it)-2*u(i,k,it))
!            else if (i==3 .or. i==4) then
!                Z=dz**(-2)/12*(-u(i+2,k,it)-u(i-2,k,it)+16*u(i+1,k,it)+16*u(i-1,k,it)-30*u(i,k,it))
!            else if (i>4 .and. i<(nz-3)) then
!                Z=dz**(-2)*(-0.00178571428571*(u(i+4,k,it)+u(i-4,k,it))+ & 
!                            0.02539682539683*(u(i+3,k,it)+u(i-3,k,it))- &
!                            0.20000000000000*(u(i+2,k,it)+u(i-2,k,it))+ &
!                            1.60000000000000*(u(i+1,k,it)+u(i-1,k,it))- &
!                            2.84722222222222*(u(i,k,it)))
!            end if
            
!            u(i,k,it-1)=-T+(Z+X)*(velw(i,k)*dt)**(2)

!            if (k<=4) then
!               u(i,k,it-1)=-(-u(i,k,it)+u(i,k+1,it)-u(i,k+1,it-1)+ &
!                   alfa/dwrk*(u(i,k+1,it)-u(i,k,it)-(u(i,k+2,it-1)-u(i,k+1,it-1))))
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            else if (k>=nwrk-3) then
!               u(i,k,it+1)=u(i,k,it)+u(i,k-1,it)-u(i,k-1,it-1)- &
!                   alfa/dwrk*(u(i,k,it)-u(i,k-1,it)-(u(i,k-1,it-1)-u(i,k-2,it-1)))
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            end if
!            
!            if (i>=nz-3) then
!               u(i,k,it+1)=u(i,k,it)+u(i-1,k,it)-u(i-1,k,it-1)- &
!                   alfa/dz*(u(i,k,it)-u(i-1,k,it)-(u(i-1,k,it-1)-u(i-2,k,it-1)))
!               u(i,k,it+1)=taper(i,k)*u(i,k,it+1)
!            end if
