Subroutine fill
Use comm

Implicit None
Real(Kind=8), Parameter :: kte=3d2*Boltzmann/elec_charge
Real(Kind=8) :: kte2
Real(Kind=8), Parameter :: dettach=1d-3
Real(Kind=8) :: bp(4),cp(4),dp(4), temp0, temp1, temp2, temp3
Integer :: j,k,l,m

kte2=electron_temperature*Boltzmann/elec_charge
phi(0,:) = phi(1,:)
phi(ny+1,:) = phi(ny,:)
phi(:,0) = 0d0
phi(:,nz+1) = pc

rne(0,:) = rne(1,:)
rne(ny+1,:) = rne(ny,:)
rne(:,0) = rne(:,1)
rne(:,nz+1) = rne(:,nz)

rnp(0,:) = rnp(1,:)
rnp(ny+1,:) = rnp(ny,:)
rnp(:,0) = rnp(:,1)
rnp(:,nz+1) = rnp(:,nz)

rnn(0,:) = rnn(1,:)
rnn(ny+1,:) = rnn(ny,:)
rnn(:,0) = rnn(:,1)
rnn(:,nz+1) = rnn(:,nz)

Call amf

ap=0d0
af=0d0
Do j=1,ny
  Do k=1,nz
    Do l=1,nv
      m=(l-1)*nn
      bp=0d0
      cp=0d0
      Select Case(l)
      Case(1)
        bp(1) = (2d0*art(j)/(dz(k)+dz(k-1)))*(epsz(j,k-1))
        bp(2) = (2d0*arn(j-1,k)/(dy(j)+dy(j-1)))*(epsy(j-1,k))
        bp(3) = (2d0*arn(j,k)/(dy(j)+dy(j+1)))*(epsy(j,k))
        bp(4) = (2d0*art(j)/(dz(k)+dz(k+1)))*(epsz(j,k))
        If (k>1) Then
          ap(ijk(j,k),1) = ap(ijk(j,k),1) + bp(1)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(1)
        Else
          af(ijk(j,1)) = af(ijk(j,1)) - phi(j,0)*bp(1)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(1)
      End If
        If (j>1) Then
          ap(ijk(j,k),2) = ap(ijk(j,k),2) + bp(2)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(2)
        Else
          phi(0,k) = phi(1,k)
        End If
        If (j<ny) Then
          ap(ijk(j,k),4) = ap(ijk(j,k),4) + bp(3)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(3)
        Else
          phi(ny+1,k) = phi(ny,k)
        End If
        If (k<nz) Then
          ap(ijk(j,k),5) = ap(ijk(j,k),5) + bp(4)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(4)
        Else
          af(ijk(j,nz)) = af(ijk(j,nz)) - phi(j,nz+1)*bp(4)
          ap(ijk(j,k),3) = ap(ijk(j,k),3) - bp(4)
        End If
        temp0 = elec_charge*vol(j,k)/eps_0
        ap(ijk(j,k),8) = ap(ijk(j,k),8) - temp0
        ap(ijk(j,k),13) = ap(ijk(j,k),13) + temp0
        ap(ijk(j,k),18) = ap(ijk(j,k),18) - temp0
      Case(2)
        bp(1) = 2d0*art(j)*emobz(j,k-1)*rnez(j,k-1)/(dz(k)+dz(k-1))
        bp(2) = 2d0*arn(j-1,k)*emoby(j-1,k)*rney(j-1,k)/(dy(j)+dy(j-1))
        bp(3) = 2d0*arn(j,k)*emoby(j,k)*rney(j,k)/(dy(j)+dy(j+1))
        bp(4) = 2d0*art(j)*emobz(j,k)*rnez(j,k)/(dz(k)+dz(k+1))

        If (ian(j,k)==1.AND.high_charge_dens) Then
          cp=0d0
          dp=0d0
        Else
          Select Case (nezcase(j,k-1))
          Case(0)
            cp(1) = -art(j)*emobz(j,k-1)*ez(j,k-1)*rnez(j,k-1)**2
          Case(1)
            cp(1) = -art(j)*emobz(j,k-1)*ez(j,k-1)*rnez(j,k-1)
          Case(2,3)
            cp(1) = -art(j)*emobz(j,k-1)*ez(j,k-1)
          End Select
          Select Case (neycase(j-1,k))
          Case(0)
            cp(2) = -arn(j-1,k)*emoby(j-1,k)*ey(j-1,k)*rney(j-1,k)**2
          Case(1)
            cp(2) = -arn(j-1,k)*emoby(j-1,k)*ey(j-1,k)*rney(j-1,k)
          Case(2,3)
            cp(2) = -arn(j-1,k)*emoby(j-1,k)*ey(j-1,k)
          End Select
          Select Case (neycase(j,k))
          Case(0)
            cp(3) = -arn(j,k)*emoby(j,k)*ey(j,k)*rney(j,k)**2
          Case(1)
            cp(3) = -arn(j,k)*emoby(j,k)*ey(j,k)*rney(j,k)
          Case(2,3)
            cp(3) = -arn(j,k)*emoby(j,k)*ey(j,k)
          End Select
          Select Case (nezcase(j,k))
          Case(0)
            cp(4) = -art(j)*emobz(j,k)*ez(j,k)*rnez(j,k)**2
          Case(1)
            cp(4) = -art(j)*emobz(j,k)*ez(j,k)*rnez(j,k)
          Case(2,3)
            cp(4) = -art(j)*emobz(j,k)*ez(j,k)
          End Select
          dp(1) = (2d0*art(j)/(dz(k)+dz(k-1)))*kte2*Abs(emobz(j,k-1))
          dp(2) = (2d0*arn(j-1,k)/(dy(j)+dy(j-1)))*kte2*Abs(emoby(j-1,k))
          dp(3) = (2d0*arn(j,k)/(dy(j)+dy(j+1)))*kte2*Abs(emoby(j,k))
          dp(4) = (2d0*art(j)/(dz(k)+dz(k+1)))*kte2*Abs(emobz(j,k))
        End If
          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) + bp(1)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(1)
            If (ian(j,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(1)/=0d0) Then
                If (nezcase(j,k-1)==0) Then
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - h(k-1)*cp(1)/rne(j,k)**2
                  If (ian(j,k-1)==1.AND.high_charge_dens) Then
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*(h(k-1)/rne(j,k))
                  Else
                    ap(ijk(j,k)+m,6) = ap(ijk(j,k)+m,6) - (1d0-h(k-1))*cp(1)/rne(j,k-1)**2
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))/rne(j,k-1)+h(k-1)/rne(j,k))
                  End If
                Else If (nezcase(j,k-1)==1) Then
                  ap(ijk(j,k)+m,6) = ap(ijk(j,k)+m,6) - (1d0-h(k-1))*cp(1)/rne(j,k-1)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - h(k-1)*cp(1)/rne(j,k)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)
                Else If (nezcase(j,k-1)==2) Then
                  ap(ijk(j,k)+m,6) = ap(ijk(j,k)+m,6) - (1d0-h(k-1))*cp(1)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - h(k-1)*cp(1)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))*rne(j,k-1)+h(k-1)*rne(j,k))
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind z lower)", nezcase(j,k-1),j,k
                  Stop
                End If
              End If
              ap(ijk(j,k)+m,6) = ap(ijk(j,k)+m,6) + dp(1)
              ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - dp(1)
            End If
          Else
            af(ijk(j,1)+m) = af(ijk(j,1)+m) - phi(j,0)*bp(1)
            ap(ijk(j,1)+m,3) = ap(ijk(j,1)+m,3) - bp(1)
            rne(j,0)=rne(j,1)
            If (ian(j,0)==0.OR..NOT.high_charge_dens) Then
              If (cp(1)/=0d0) Then
                If (nezcase(j,0)==0) Then
                  ap(ijk(j,1)+m,8) = ap(ijk(j,1)+m,8) - cp(1)/rne(j,1)**2
                  af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)/rne(j,1)
                Else If (nezcase(j,0)==1) Then
                  ap(ijk(j,1)+m,8) = ap(ijk(j,1)+m,8) - cp(1)/rne(j,1)
                  af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)
                Else If (nezcase(j,0)==2) Then
                  ap(ijk(j,1)+m,8) = ap(ijk(j,1)+m,8) - cp(1)
                  af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)*rne(j,1)
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind z lower b)"
                  Stop
                End If
              End If
            End If
          End If
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) + bp(2)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(2)
            If (ian(j,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(2)/=0d0) Then
                If (neycase(j-1,k)==0) Then
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - g(j-1)*cp(2)/rne(j,k)**2
                  If (ian(j-1,k)==1.AND.high_charge_dens) Then
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*(g(j-1)/rne(j,k))
                  Else
                    ap(ijk(j,k)+m,7) = ap(ijk(j,k)+m,7) - (1d0-g(j-1))*cp(2)/rne(j-1,k)**2
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))/rne(j-1,k)+g(j-1)/rne(j,k))
                  End If
                Else If (neycase(j-1,k)==1) Then
                  ap(ijk(j,k)+m,7) = ap(ijk(j,k)+m,7) - (1d0-g(j-1))*cp(2)/rne(j-1,k)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - g(j-1)*cp(2)/rne(j,k)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)
                Else If (neycase(j-1,k)==2) Then
                  ap(ijk(j,k)+m,7) = ap(ijk(j,k)+m,7) - (1d0-g(j-1))*cp(2)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - g(j-1)*cp(2)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))*rne(j-1,k)+g(j-1)*rne(j,k))
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind y left)", neycase(j-1,k),j,k
                  Stop
                End If
              End If
              ap(ijk(j,k)+m,7) = ap(ijk(j,k)+m,7) + dp(2)
              ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - dp(2)
            End If
          Else
            rne(0,k)=rne(1,k)
            If (ian(0,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(2)/=0d0) Then
                If (neycase(0,k)==0) Then
                  ap(ijk(1,k)+m,8) = ap(ijk(1,k)+m,8) - cp(2)/rne(1,k)**2
                  af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)/rne(1,k)
                Else If (neycase(0,k)==1) Then
                  ap(ijk(1,k)+m,8) = ap(ijk(1,k)+m,8) - cp(2)/rne(1,k)
                  af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)
                Else If (neycase(0,k)==2) Then
                  ap(ijk(1,k)+m,8) = ap(ijk(1,k)+m,8) - cp(2)
                  af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)*rne(1,k)
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind y left b)"
                  Stop
                End If
              End If
            End If
          End If
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) + bp(3)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(3)
            If (ian(j,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(3)/=0d0) Then
                If (neycase(j,k)==0) Then
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-g(j))*cp(3)/rne(j,k)**2
                  If (ian(j+1,k)==1.AND.high_charge_dens) Then
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))/rne(j,k))
                  Else
                    ap(ijk(j,k)+m,9) = ap(ijk(j,k)+m,9) + g(j)*cp(3)/rne(j+1,k)**2
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))/rne(j,k)+g(j)/rne(j+1,k))
                  End If
                Else If (neycase(j,k)==1) Then
                  ap(ijk(j,k)+m,9) = ap(ijk(j,k)+m,9) + g(j)*cp(3)/rne(j+1,k)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-g(j))*cp(3)/rne(j,k)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)
                Else If (neycase(j,k)==2) Then
                  ap(ijk(j,k)+m,9) = ap(ijk(j,k)+m,9) + g(j)*cp(3)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-g(j))*cp(3)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))*rne(j,k)+g(j)*rne(j+1,k))
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind y right)",neycase(j,k),j,k
                  Stop
                End If
              End If
              ap(ijk(j,k)+m,9) = ap(ijk(j,k)+m,9) + dp(3)
              ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - dp(3)
            End If
          Else
            rne(ny+1,k)=rne(ny,k)
            If (ian(ny+1,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(3)/=0d0) Then
                If (neycase(ny,k)==0) Then
                  ap(ijk(ny,k)+m,8) = ap(ijk(ny,k)+m,8) + cp(3)/rne(ny,k)**2
                  af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)/rne(ny,k)
                Else If (neycase(ny,k)==1) Then
                  ap(ijk(ny,k)+m,8) = ap(ijk(ny,k)+m,8) + cp(3)/rne(ny,k)
                  af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)
                Else If (neycase(ny,k)==2) Then
                  ap(ijk(ny,k)+m,8) = ap(ijk(ny,k)+m,8) + cp(3)
                  af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)*rne(ny,k)
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind y right b)"
                  Stop
                End If
              End If
            End If
          End If

          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) + bp(4)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(4)
            If (ian(j,k)==0.OR..NOT.high_charge_dens) Then
              If (cp(4)/=0d0) Then
                If (nezcase(j,k)==0) Then
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-h(k))*cp(4)/rne(j,k)**2
                  If (ian(j,k+1)==1.AND.high_charge_dens) Then
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))/rne(j,k))
                  Else
                    ap(ijk(j,k)+m,10) = ap(ijk(j,k)+m,10) + h(k)*cp(4)/rne(j,k+1)**2
                    af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))/rne(j,k)+h(k)/rne(j,k+1))
                  End If
                Else If (nezcase(j,k)==1) Then
                  ap(ijk(j,k)+m,10) = ap(ijk(j,k)+m,10) + h(k)*cp(4)/rne(j,k+1)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-h(k))*cp(4)/rne(j,k)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)
                Else If (nezcase(j,k)==2) Then
                  ap(ijk(j,k)+m,10) = ap(ijk(j,k)+m,10) + h(k)*cp(4)
                  ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (1d0-h(k))*cp(4)
                  af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))*rne(j,k)+h(k)*rne(j,k+1))
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind z top)",nezcase(j,k),j,k
                  Stop
                End If
              End If
              ap(ijk(j,k)+m,10) = ap(ijk(j,k)+m,10) + dp(4)
              ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) - dp(4)
            End If
          Else
            af(ijk(j,nz)+m) = af(ijk(j,nz)+m) - phi(j,nz+1)*bp(4)
            ap(ijk(j,nz)+m,3) = ap(ijk(j,nz)+m,3) - bp(4)
            rne(j,nz+1)=rne(j,nz)
            If (ian(j,k)==0.OR..NOT.high_charge_dens) Then
            If (cp(4)/=0d0) Then
                If (nezcase(j,nz)==0) Then
                  ap(ijk(j,nz)+m,8) = ap(ijk(j,nz)+m,8) + cp(4)/rne(j,nz)**2
                  af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)/rne(j,nz)
                Else If (nezcase(j,nz)==1) Then
                  ap(ijk(j,nz)+m,8) = ap(ijk(j,nz)+m,8) + cp(4)/rne(j,nz)
                  af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)
                Else If (nezcase(j,nz)==2) Then
                  ap(ijk(j,nz)+m,8) = ap(ijk(j,nz)+m,8) + cp(4)
                  af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)*rne(j,nz)
                Else
                  Write(*,*) "fill: Error! This section is incomplete... (upwind z top b)"
                  Stop
                End If
              End If
            End If
          End If
        If (ian(j,k)==0) Then
          ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + ((alpha(j,k)-eta(j,k)) - gamma*rnp(j,k))*vol(j,k)
          ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - gamma*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - gamma*rne(j,k)*rnp(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - gamma*vol(j,k)*cp_ion2**2
          ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + dettach*vol(j,k)
          temp0 = -2d0*(dy(j)+dy(j+1))/((dy(j-1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp1 = 2d0*(dy(j+1)-dy(j-1))/((dy(j-1)+dy(j))*(dy(j+1)+dy(j)))
          temp2 = 2d0*(dy(j)+dy(j-1))/((dy(j+1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp3 = (detay(j,k)-dalphay(j,k))*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j-1,k)+temp1*phi(j,k)+temp2*phi(j+1,k))
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
          temp0 = -2d0*(dz(k)+dz(k+1))/((dz(k-1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp1 = 2d0*(dz(k+1)-dz(k-1))/((dz(k-1)+dz(k))*(dz(k+1)+dz(k)))
          temp2 = 2d0*(dz(k)+dz(k-1))/((dz(k+1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp3 = (detaz(j,k)-dalphaz(j,k))*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j,k-1)+temp1*phi(j,k)+temp2*phi(j,k+1))
          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
        End If
      Case(3)
        If (ian(j,k)==1) Then
          ap(ijk(j,k)+m,13) = -1d0
          If (high_charge_dens) Then
            af(ijk(j,k)+m) = 0d0
          Else
            af(ijk(j,k)+m) = -cp_ion
          End If
        Else
          bp(1) = 2d0*art(j)*pmobz(j,k-1)*rnpz(j,k-1)/(dz(k)+dz(k-1))
          bp(2) = 2d0*arn(j-1,k)*pmoby(j-1,k)*rnpy(j-1,k)/(dy(j)+dy(j-1))
          bp(3) = 2d0*arn(j,k)*pmoby(j,k)*rnpy(j,k)/(dy(j)+dy(j+1))
          bp(4) = 2d0*art(j)*pmobz(j,k)*rnpz(j,k)/(dz(k)+dz(k+1))

        Select Case (npzcase(j,k-1))
        Case(0)
          cp(1) = -art(j)*pmobz(j,k-1)*ez(j,k-1)*rnpz(j,k-1)**2
        Case(1)
          cp(1) = -art(j)*pmobz(j,k-1)*ez(j,k-1)*rnpz(j,k-1)
        Case(2,3)
          cp(1) = -art(j)*pmobz(j,k-1)*ez(j,k-1)
        End Select
        Select Case (npycase(j-1,k))
        Case(0)
          cp(2) = -arn(j-1,k)*pmoby(j-1,k)*ey(j-1,k)*rnpy(j-1,k)**2
        Case(1)
          cp(2) = -arn(j-1,k)*pmoby(j-1,k)*ey(j-1,k)*rnpy(j-1,k)
        Case(2,3)
          cp(2) = -arn(j-1,k)*pmoby(j-1,k)*ey(j-1,k)
        End Select
        Select Case (npycase(j,k))
        Case(0)
          cp(3) = -arn(j,k)*pmoby(j,k)*ey(j,k)*rnpy(j,k)**2
        Case(1)
          cp(3) = -arn(j,k)*pmoby(j,k)*ey(j,k)*rnpy(j,k)
        Case(2,3)
          cp(3) = -arn(j,k)*pmoby(j,k)*ey(j,k)
        End Select
        Select Case (npzcase(j,k))
        Case(0)
          cp(4) = -art(j)*pmobz(j,k)*ez(j,k)*rnpz(j,k)**2
        Case(1)
          cp(4) = -art(j)*pmobz(j,k)*ez(j,k)*rnpz(j,k)
        Case(2,3)
          cp(4) = -art(j)*pmobz(j,k)*ez(j,k)
        End Select
          dp(1) = (2d0*art(j)/(dz(k)+dz(k-1)))*kte*Abs(pmobz(j,k-1))
          dp(2) = (2d0*arn(j-1,k)/(dy(j)+dy(j-1)))*kte*Abs(pmoby(j-1,k))
          dp(3) = (2d0*arn(j,k)/(dy(j)+dy(j+1)))*kte*Abs(pmoby(j,k))
          dp(4) = (2d0*art(j)/(dz(k)+dz(k+1)))*kte*Abs(pmobz(j,k))

          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) + bp(1)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(1)
            If (cp(1)/=0d0) Then
              If (npzcase(j,k-1)==0) Then
                ap(ijk(j,k)+m,11) = ap(ijk(j,k)+m,11) - (1d0-h(k-1))*cp(1)/rnp(j,k-1)**2
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - h(k-1)*cp(1)/rnp(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))/rnp(j,k-1)+h(k-1)/rnp(j,k))
              Else If (npzcase(j,k-1)==1) Then
                ap(ijk(j,k)+m,11) = ap(ijk(j,k)+m,11) - (1d0-h(k-1))*cp(1)/rnp(j,k-1)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - h(k-1)*cp(1)/rnp(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)
              Else If (npzcase(j,k-2)==1) Then
                ap(ijk(j,k)+m,11) = ap(ijk(j,k)+m,11) - (1d0-h(k-1))*cp(1)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - h(k-1)*cp(1)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))*rnp(j,k-1)+h(k-1)*rnp(j,k))
              End If
            End If
            ap(ijk(j,k)+m,11) = ap(ijk(j,k)+m,11) + dp(1)
            ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - dp(1)
          Else
            af(ijk(j,1)+m) = af(ijk(j,1)+m) - phi(j,0)*bp(1)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(1)
            rnp(j,0)=rnp(j,1)
            If (cp(1)/=0d0) Then
              If (npzcase(j,0)==0) Then
                ap(ijk(j,1)+m,13) = ap(ijk(j,1)+m,13) - cp(1)/rnp(j,1)**2
                af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)/rnp(j,1)
              Else If (npzcase(j,0)==1) Then
                ap(ijk(j,1)+m,13) = ap(ijk(j,1)+m,13) - cp(1)/rnp(j,1)
                af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)
              Else If (npzcase(j,0)==2) Then
                ap(ijk(j,1)+m,13) = ap(ijk(j,1)+m,13) - cp(1)
                af(ijk(j,1)+m) = af(ijk(j,1)+m) - cp(1)*rnp(j,1)
              End If
            End If
          End If
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) + bp(2)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(2)
            If (cp(2)/=0d0) Then
              If (npycase(j-1,k)==0) Then
                ap(ijk(j,k)+m,12) = ap(ijk(j,k)+m,12) - (1d0-g(j-1))*cp(2)/rnp(j-1,k)**2
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - g(j-1)*cp(2)/rnp(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))/rnp(j-1,k)+g(j-1)/rnp(j,k))
              Else If (npycase(j-1,k)==1) Then
                ap(ijk(j,k)+m,12) = ap(ijk(j,k)+m,12) - (1d0-g(j-1))*cp(2)/rnp(j-1,k)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - g(j-1)*cp(2)/rnp(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)
              Else If (npycase(j-1,k)==2) Then
                ap(ijk(j,k)+m,12) = ap(ijk(j,k)+m,12) - (1d0-g(j-1))*cp(2)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - g(j-1)*cp(2)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))*rnp(j-1,k)+g(j-1)*rnp(j,k))
              End If
            End If
            ap(ijk(j,k)+m,12) = ap(ijk(j,k)+m,12) + dp(2)
            ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - dp(2)
          Else
            rnp(0,k)=rnp(1,k)
            If (cp(2)/=0d0) Then
              If (npycase(0,k)==0) Then
                ap(ijk(1,k)+m,13) = ap(ijk(1,k)+m,13) - cp(2)/rnp(1,k)**2
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)/rnp(1,k)
              Else If (npycase(0,k)==1) Then
                ap(ijk(1,k)+m,13) = ap(ijk(1,k)+m,13) - cp(2)/rnp(1,k)
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)
              Else If (npycase(0,k)==2) Then
                ap(ijk(1,k)+m,13) = ap(ijk(1,k)+m,13) - cp(2)
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)*rnp(1,k)
              End If
            End If
          End If
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) + bp(3)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(3)
            If (cp(3)/=0d0) Then
              If (npycase(j,k)==0) Then
                ap(ijk(j,k)+m,14) = ap(ijk(j,k)+m,14) + g(j)*cp(3)/rnp(j+1,k)**2
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-g(j))*cp(3)/rnp(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))/rnp(j,k)+g(j)/rnp(j+1,k))
              Else If (npycase(j,k)==1) Then
                ap(ijk(j,k)+m,14) = ap(ijk(j,k)+m,14) + g(j)*cp(3)/rnp(j+1,k)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-g(j))*cp(3)/rnp(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)
              Else If (npycase(j,k)==2) Then
                ap(ijk(j,k)+m,14) = ap(ijk(j,k)+m,14) + g(j)*cp(3)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-g(j))*cp(3)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))*rnp(j,k)+g(j)*rnp(j+1,k))
              End If
            End If
            ap(ijk(j,k)+m,14) = ap(ijk(j,k)+m,14) + dp(3)
            ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - dp(3)
          Else
            rnp(ny+1,k)=rnp(ny,k)
            If (cp(3)/=0d0) Then
              If (npycase(ny,k)==0) Then
                ap(ijk(ny,k)+m,13) = ap(ijk(ny,k)+m,13) + cp(3)/rnp(ny,k)**2
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)/rnp(ny,k)
              Else If (npycase(ny,k)==1) Then
                ap(ijk(ny,k)+m,13) = ap(ijk(ny,k)+m,13) + cp(3)/rnp(ny,k)
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)
              Else If (npycase(ny,k)==2) Then
                ap(ijk(ny,k)+m,13) = ap(ijk(ny,k)+m,13) + cp(3)
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)*rnp(ny,k)
              End If
            End If
          End If
          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) + bp(4)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(4)
            If (cp(4)/=0d0) Then
              If (npzcase(j,k)==0) Then
                ap(ijk(j,k)+m,15) = ap(ijk(j,k)+m,15) + h(k)*cp(4)/rnp(j,k+1)**2
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-h(k))*cp(4)/rnp(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))/rnp(j,k)+h(k)/rnp(j,k+1))
              Else If (npzcase(j,k)==1) Then
                ap(ijk(j,k)+m,15) = ap(ijk(j,k)+m,15) + h(k)*cp(4)/rnp(j,k+1)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-h(k))*cp(4)/rnp(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)
              Else If (npzcase(j,k)==2) Then
                ap(ijk(j,k)+m,15) = ap(ijk(j,k)+m,15) + h(k)*cp(4)
                ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) + (1d0-h(k))*cp(4)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))*rnp(j,k)+h(k)*rnp(j,k+1))
              End If
            End If
            ap(ijk(j,k)+m,15) = ap(ijk(j,k)+m,15) + dp(4)
            ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - dp(4)
          Else
            af(ijk(j,nz)+m) = af(ijk(j,nz)+m) - phi(j,nz+1)*bp(4)
            ap(ijk(j,nz)+m,3) = ap(ijk(j,nz)+m,3) - bp(4)
            rnp(j,nz+1)=rnp(j,nz)
            If (cp(4)/=0d0) Then
              If (npzcase(j,nz)==0) Then
                ap(ijk(j,nz)+m,13) = ap(ijk(j,nz)+m,13) + cp(4)/rnp(j,nz)**2
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)/rnp(j,nz)
              Else If (npzcase(j,nz)==1) Then
                ap(ijk(j,nz)+m,13) = ap(ijk(j,nz)+m,13) + cp(4)/rnp(j,nz)
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)
              Else If (npzcase(j,nz)==2) Then
                ap(ijk(j,nz)+m,13) = ap(ijk(j,nz)+m,13) + cp(4)
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)*rnp(j,nz)
              End If
            End If
          End If
          ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + (alpha(j,k) - gamma*rnp(j,k))*vol(j,k)
          ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - gamma*(rne(j,k)+rnn(j,k))*vol(j,k)
          ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - gamma*rnp(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - gamma*(rne(j,k)+rnn(j,k))*rnp(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - gamma*vol(j,k)*cp_ion2**2
          temp0 = -2d0*(dy(j)+dy(j+1))/((dy(j-1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp1 = 2d0*(dy(j+1)-dy(j-1))/((dy(j-1)+dy(j))*(dy(j+1)+dy(j)))
          temp2 = 2d0*(dy(j)+dy(j-1))/((dy(j+1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp3 = -dalphay(j,k)*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j-1,k)+temp1*phi(j,k)+temp2*phi(j+1,k))
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
          temp0 = -2d0*(dz(k)+dz(k+1))/((dz(k-1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp1 = 2d0*(dz(k+1)-dz(k-1))/((dz(k-1)+dz(k))*(dz(k+1)+dz(k)))
          temp2 = 2d0*(dz(k)+dz(k-1))/((dz(k+1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp3 = -dalphaz(j,k)*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j,k-1)+temp1*phi(j,k)+temp2*phi(j,k+1))
          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
        End If
      Case(4)
        If (ian(j,k)==1) Then
          ap(ijk(j,k)+m,18) = -1d30
          af(ijk(j,k)+m) = 0d0
          rnn(j,k)=0d0
        Else
          bp(1) = 2d0*art(j)*nmobz(j,k-1)*rnnz(j,k-1)/(dz(k)+dz(k-1))
          bp(2) = 2d0*arn(j-1,k)*nmoby(j-1,k)*rnny(j-1,k)/(dy(j)+dy(j-1))
          bp(3) = 2d0*arn(j,k)*nmoby(j,k)*rnny(j,k)/(dy(j)+dy(j+1))
          bp(4) = 2d0*art(j)*nmobz(j,k)*rnnz(j,k)/(dz(k)+dz(k+1))

        Select Case (nnzcase(j,k-1))
        Case(0)
          cp(1) = -art(j)*nmobz(j,k-1)*ez(j,k-1)*rnnz(j,k-1)**2
        Case(1)
          cp(1) = -art(j)*nmobz(j,k-1)*ez(j,k-1)*rnnz(j,k-1)
        Case(2,3)
          cp(1) = -art(j)*nmobz(j,k-1)*ez(j,k-1)
        End Select
        Select Case (nnycase(j-1,k))
        Case(0)
          cp(2) = -arn(j-1,k)*nmoby(j-1,k)*ey(j-1,k)*rnny(j-1,k)**2
        Case(1)
          cp(2) = -arn(j-1,k)*nmoby(j-1,k)*ey(j-1,k)*rnny(j-1,k)
        Case(2,3)
          cp(2) = -arn(j-1,k)*nmoby(j-1,k)*ey(j-1,k)
        End Select
        Select Case (nnycase(j,k))
        Case(0)
          cp(3) = -arn(j,k)*nmoby(j,k)*ey(j,k)*rnny(j,k)**2
        Case(1)
          cp(3) = -arn(j,k)*nmoby(j,k)*ey(j,k)*rnny(j,k)
        Case(2,3)
          cp(3) = -arn(j,k)*nmoby(j,k)*ey(j,k)
        End Select
        Select Case (nnzcase(j,k))
        Case(0)
          cp(4) = -art(j)*nmobz(j,k)*ez(j,k)*rnnz(j,k)**2
        Case(1)
          cp(4) = -art(j)*nmobz(j,k)*ez(j,k)*rnnz(j,k)
        Case(2,3)
          cp(4) = -art(j)*nmobz(j,k)*ez(j,k)
        End Select
          dp(1) = (2d0*art(j)/(dz(k)+dz(k-1)))*kte*Abs(nmobz(j,k-1))
          dp(2) = (2d0*arn(j-1,k)/(dy(j)+dy(j-1)))*kte*Abs(nmoby(j-1,k))
          dp(3) = (2d0*arn(j,k)/(dy(j)+dy(j+1)))*kte*Abs(nmoby(j,k))
          dp(4) = (2d0*art(j)/(dz(k)+dz(k+1)))*kte*Abs(nmobz(j,k))
          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) + bp(1)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(1)
            If (cp(1)/=0d0) Then
              If (nnzcase(j,k-1)==0) Then
                ap(ijk(j,k)+m,16) = ap(ijk(j,k)+m,16) - (1d0-h(k-1))*cp(1)/rnn(j,k-1)**2
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - h(k-1)*cp(1)/rnn(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))/rnn(j,k-1)+h(k-1)/rnn(j,k))
              Else If (nnzcase(j,k-1)==1) Then
                ap(ijk(j,k)+m,16) = ap(ijk(j,k)+m,16) - (1d0-h(k-1))*cp(1)/rnn(j,k-1)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - h(k-1)*cp(1)/rnn(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)
              Else If (nnzcase(j,k-1)==1) Then
                ap(ijk(j,k)+m,16) = ap(ijk(j,k)+m,16) - (1d0-h(k-1))*cp(1)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - h(k-1)*cp(1)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(1)*((1d0-h(k-1))*rnn(j,k-1)+h(k-1)*rnn(j,k))
              End If
            End If
            ap(ijk(j,k)+m,16) = ap(ijk(j,k)+m,16) + dp(1)
            ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - dp(1)
          Else
            af(ijk(j,1)+m) = af(ijk(j,1)+m) - phi(j,0)*bp(1)
            ap(ijk(j,1)+m,3) = ap(ijk(j,1)+m,3) - bp(1)
            rnn(j,0)=rnn(j,1)
            If (cp(1)/=0d0) Then
              If (nnzcase(j,0)==0) Then
                ap(ijk(j,1)+m,18) = ap(ijk(j,1)+m,18) - cp(1)/rnn(j,1)**2
                af(ijk(j,1)+m) = af(ijk(j,1)+m) + cp(1)/rnn(j,1)
              Else If (nnzcase(j,0)==1) Then
                ap(ijk(j,1)+m,18) = ap(ijk(j,1)+m,18) - cp(1)/rnn(j,1)
                af(ijk(j,1)+m) = af(ijk(j,1)+m) + cp(1)
              Else If (nnzcase(j,0)==2) Then
                ap(ijk(j,1)+m,18) = ap(ijk(j,1)+m,18) - cp(1)
                af(ijk(j,1)+m) = af(ijk(j,1)+m) + cp(1)*rnn(j,1)
              End If
            End If
          End If
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) + bp(2)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(2)
            If (cp(2)/=0d0) Then
              If (nnycase(j-1,k)==0) Then
                ap(ijk(j,k)+m,17) = ap(ijk(j,k)+m,17) - (1d0-g(j-1))*cp(2)/rnn(j-1,k)**2
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - g(j-1)*cp(2)/rnn(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))/rnn(j-1,k)+g(j-1)/rnn(j,k))
              Else If (nnycase(j-1,k)==1) Then
                ap(ijk(j,k)+m,17) = ap(ijk(j,k)+m,17) - (1d0-g(j-1))*cp(2)/rnn(j-1,k)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - g(j-1)*cp(2)/rnn(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)
              Else If (nnycase(j-1,k)==2) Then
                ap(ijk(j,k)+m,17) = ap(ijk(j,k)+m,17) - (1d0-g(j-1))*cp(2)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - g(j-1)*cp(2)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) - cp(2)*((1d0-g(j-1))*rnn(j-1,k)+g(j-1)*rnn(j,k))
              End If
            End If
            ap(ijk(j,k)+m,17) = ap(ijk(j,k)+m,17) + dp(2)
            ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - dp(2)
          Else
            rnn(0,k)=rnn(1,k)
            If (cp(2)/=0d0) Then
              If (nnycase(0,k)==0) Then
                ap(ijk(1,k)+m,18) = ap(ijk(1,k)+m,18) - cp(2)/rnn(1,k)**2
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)/rnn(1,k)
              Else If (nnycase(0,k)==1) Then
                ap(ijk(1,k)+m,18) = ap(ijk(1,k)+m,18) - cp(2)/rnn(1,k)
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)
              Else If (nnycase(0,k)==2) Then
                ap(ijk(1,k)+m,18) = ap(ijk(1,k)+m,18) - cp(2)
                af(ijk(1,k)+m) = af(ijk(1,k)+m) - cp(2)*rnn(1,k)
              End If
            End If
          End If
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) + bp(3)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(3)
            If (cp(3)/=0d0) Then
              If (nnycase(j,k)==0) Then
                ap(ijk(j,k)+m,19) = ap(ijk(j,k)+m,19) + g(j)*cp(3)/rnn(j+1,k)**2
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-g(j))*cp(3)/rnn(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))/rnn(j,k)+g(j)/rnn(j+1,k))
              Else If (nnycase(j,k)==1) Then
                ap(ijk(j,k)+m,19) = ap(ijk(j,k)+m,19) + g(j)*cp(3)/rnn(j+1,k)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-g(j))*cp(3)/rnn(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)
              Else If (nnycase(j,k)==2) Then
                ap(ijk(j,k)+m,19) = ap(ijk(j,k)+m,19) + g(j)*cp(3)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-g(j))*cp(3)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(3)*((1d0-g(j))*rnn(j,k)+g(j)*rnn(j+1,k))
              End If
            End If
            ap(ijk(j,k)+m,19) = ap(ijk(j,k)+m,19) + dp(3)
            ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - dp(3)
          Else
            rnn(ny+1,k)=rnn(ny,k)
            If (cp(3)/=0d0) Then
              If (nnycase(ny,k)==0) Then
                ap(ijk(ny,k)+m,18) = ap(ijk(ny,k)+m,18) + cp(3)/rnn(ny,k)**2
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)/rnn(ny,k)
              Else If (nnycase(ny,k)==1) Then
                ap(ijk(ny,k)+m,18) = ap(ijk(ny,k)+m,18) + cp(3)/rnn(ny,k)
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)
              Else If (nnycase(ny,k)==2) Then
                ap(ijk(ny,k)+m,18) = ap(ijk(ny,k)+m,18) + cp(3)
                af(ijk(ny,k)+m) = af(ijk(ny,k)+m) + cp(3)*rnn(ny,k)
              End If
            End If
          End If
          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) + bp(4)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(4)
            If (cp(4)/=0d0) Then
              If (nnzcase(j,k)==0) Then
                ap(ijk(j,k)+m,20) = ap(ijk(j,k)+m,20) + h(k)*cp(4)/rnn(j,k+1)**2
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-h(k))*cp(4)/rnn(j,k)**2
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))/rnn(j,k)+h(k)/rnn(j,k+1))
              Else If (nnzcase(j,k)==1) Then
                ap(ijk(j,k)+m,20) = ap(ijk(j,k)+m,20) + h(k)*cp(4)/rnn(j,k+1)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-h(k))*cp(4)/rnn(j,k)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)
              Else If (nnzcase(j,k)==2) Then
                ap(ijk(j,k)+m,20) = ap(ijk(j,k)+m,20) + h(k)*cp(4)
                ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) + (1d0-h(k))*cp(4)
                af(ijk(j,k)+m) = af(ijk(j,k)+m) + cp(4)*((1d0-h(k))*rnn(j,k)+h(k)*rnn(j,k+1))
              End If
            End If
            ap(ijk(j,k)+m,20) = ap(ijk(j,k)+m,20) + dp(4)
            ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - dp(4)
          Else
            af(ijk(j,nz)+m) = af(ijk(j,nz)+m) - phi(j,nz+1)*bp(4)
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - bp(4)
            rnn(j,nz+1)=rnn(j,nz)
            If (cp(4)/=0d0) Then
              If (nnzcase(j,nz)==0) Then
                ap(ijk(j,nz)+m,18) = ap(ijk(j,nz)+m,18) + cp(4)/rnn(j,nz)**2
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)/rnn(j,nz)
              Else If (nnzcase(j,nz)==1) Then
                ap(ijk(j,nz)+m,18) = ap(ijk(j,nz)+m,18) + cp(4)/rnn(j,nz)
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)
              Else If (nnzcase(j,nz)==2) Then
                ap(ijk(j,nz)+m,18) = ap(ijk(j,nz)+m,18) + cp(4)
                af(ijk(j,nz)+m) = af(ijk(j,nz)+m) + cp(4)*rnn(j,nz)
              End If
            End If
          End If
          ap(ijk(j,k)+m,8) = ap(ijk(j,k)+m,8) + eta(j,k)*vol(j,k)
          ap(ijk(j,k)+m,13) = ap(ijk(j,k)+m,13) - gamma*rnn(j,k)*vol(j,k)
          ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - gamma*rnp(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - gamma*rnn(j,k)*rnp(j,k)*vol(j,k)
          ap(ijk(j,k)+m,18) = ap(ijk(j,k)+m,18) - dettach*vol(j,k)
          temp0 = -2d0*(dy(j)+dy(j+1))/((dy(j-1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp1 = 2d0*(dy(j+1)-dy(j-1))/((dy(j-1)+dy(j))*(dy(j+1)+dy(j)))
          temp2 = 2d0*(dy(j)+dy(j-1))/((dy(j+1)+dy(j))*(dy(j-1)+2d0*dy(j)+dy(j+1)))
          temp3 = -detay(j,k)*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j-1,k)+temp1*phi(j,k)+temp2*phi(j+1,k))
          If (j>1) Then
            ap(ijk(j,k)+m,2) = ap(ijk(j,k)+m,2) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (j<ny) Then
            ap(ijk(j,k)+m,4) = ap(ijk(j,k)+m,4) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
          temp0 = -2d0*(dz(k)+dz(k+1))/((dz(k-1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp1 = 2d0*(dz(k+1)-dz(k-1))/((dz(k-1)+dz(k))*(dz(k+1)+dz(k)))
          temp2 = 2d0*(dz(k)+dz(k-1))/((dz(k+1)+dz(k))*(dz(k-1)+2d0*dz(k)+dz(k+1)))
          temp3 = -detaz(j,k)*rne(j,k)*vol(j,k)
          af(ijk(j,k)+m) = af(ijk(j,k)+m) - temp3*(temp0*phi(j,k-1)+temp1*phi(j,k)+temp2*phi(j,k+1))
          If (k>1) Then
            ap(ijk(j,k)+m,1) = ap(ijk(j,k)+m,1) - temp3*temp0
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp0
          End If
          ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp1
          If (k<nz) Then
            ap(ijk(j,k)+m,5) = ap(ijk(j,k)+m,5) - temp3*temp2
          Else
            ap(ijk(j,k)+m,3) = ap(ijk(j,k)+m,3) - temp3*temp2
          End If
        End If
      End Select
    End Do
  End Do
End Do

End Subroutine fill
