纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b

!*************************************************************
!*************************************************************
! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b
!*************************************************************
! SINOMACH
!
! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之b
!
! 周勇 20221125
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Subroutine hcell_dof_set(region_density_x, region_density_y, maxnd,&
node_num, node_dof_index)

!*********************************************************************72
!
!c HCELL_DOF_SET assigns degrees of freedom to each node.
!
Implicit None

Integer maxnd

Integer col
Integer dof
Logical found
Integer i
Integer j
Integer node
Integer node_dof_index(maxnd, 3)
Integer node_num
Integer region_density_x(3)
Integer region_density_y(3)
Integer row

node = 0
dof = 0
j = 0
!
! Working in column 1, and NOT handling nodes in extreme right.
!
Do col = 1, 2*region_density_x(1)
!
! Working in row 1.
!
j = j + 1
i = 0

Do row = 1, 2*region_density_y(1) + 1

  i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1

Do row = 1, 2*region_density_y(3) + 1

  i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

End Do
!
! Working in column 2, including nodes on extreme right and left.
!
Do col = 1, 2*region_density_x(2) + 1
!
! Working in row 1.
!
j = j + 1
i = 0

Do row = 1, 2*region_density_y(1) + 1

  i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

!
! Working in row 2.
!
Do row = 2, 2*region_density_y(2)

i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1

i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

End Do
!
! Working in column 3, and NOT handling nodes in extreme left.
!
Do col = 2, 2*region_density_x(3) + 1
!
! Working in row 1.
!
j = j + 1
i = 0

Do row = 1, 2*region_density_y(1) + 1

  i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1

Do row = 1, 2*region_density_y(3) + 1

  i = i + 1
  node = node + 1

  dof = dof + 1
  node_dof_index(node, 1) = dof

  dof = dof + 1
  node_dof_index(node, 2) = dof

  If (mod(j,2)==1 .And. mod(i,2)==1) Then
    dof = dof + 1
    node_dof_index(node, 3) = dof
  Else
    node_dof_index(node, 3) = 0
  End If

End Do

End Do
!
! Replace one continuity equation by a pressure = 0 equation.
!
found = .False.

Do node = 1, node_num

If (0<node_dof_index(node,3)) Then
  node_dof_index(node, 3) = -1
  found = .True.
  Exit
End If

End Do

If (.Not. found) Then
Write (, '(a)') ' '
Write (, '(a)') 'HCELL_DOF_SET - Fatal errorc'
Write (, '(a)') ' Could not find a degree of freedom'
Write (, '(a)') ' associated with the continuity equation.'
Stop
End If

Return
End Subroutine hcell_dof_set
Subroutine hcell_element_count(region_density_x, region_density_y, element_num)

!*********************************************************************72
!
!c HCELL_ELEMENT_COUNT determines the number of elements in the region.
!
Implicit None

Integer element_num
Integer region_density_x(3)
Integer region_density_y(3)

element_num = 2*region_density_x(1)region_density_y(1) + &
2region_density_x(1)region_density_y(3) + &
2region_density_x(2)region_density_y(1) +&
2region_density_x(2)region_density_y(2) + &
2region_density_x(2)region_density_y(3) +&
2region_density_x(3)region_density_y(1) + &
2region_density_x(3)*region_density_y(3)

Write (, '(a)') ' '
Write (, '(a)') 'HCELL_ELEMENT_COUNT:'
Write (*, '(a,i6)') ' Number of elements = ', element_num

Return
End Subroutine hcell_element_count
Subroutine hcell_element_node(region_density_x, region_density_y,&
maxel, element_num, nnodes, element_node)

!*********************************************************************72
!
!c HCELL_ELEMENT_NODE determines the nodes that make up each element.
!
Implicit None

Integer maxel
Integer nnodes

Integer col
Integer col2
Integer element
Integer element_node(maxel, nnodes)
Integer element_num
Integer inc1
Integer inc2
Integer node_sw
Integer region_density_x(3)
Integer region_density_y(3)
Integer row
Integer row2

element = 0

Do col = 1, 3
Do col2 = 1, region_density_x(col)

Do row = 1, 3

    If (row/=2 .Or. col==2) Then

      If (col==1) Then

        If (col2<region_density_x(1)) Then

          If (row==1) Then

            If (col2==1) Then
              node_sw = 1
            Else
              node_sw = node_sw + inc1 + 1
            End If

          Else

            node_sw = node_sw + 1

          End If

          inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

        Else If (row==1) Then

          node_sw = node_sw + inc1 + 1

          inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

        Else If (row==3) Then

          node_sw = node_sw + 1

          inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
          + (2*region_density_y(3)+1)

        End If

      Else If (col==2) Then

        If (row==1) Then

          node_sw = node_sw + inc1 + 1

        End If

        inc1 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
        + (2*region_density_y(3)+1)

        inc2 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
        + (2*region_density_y(3)+1)

      Else If (col==3) Then

        If (1<col2) Then

          If (row==1) Then
            node_sw = node_sw + inc1 + 1
          Else
            node_sw = node_sw + 1
          End If

          inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

        Else If (row==1) Then

          node_sw = node_sw + inc1 + 1

          inc1 = (2*region_density_y(1)+1) +&
          (2*region_density_y(2)-1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

        Else If (row==3) Then

          node_sw = node_sw + (2*region_density_y(2)-1) + 1

          inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

          inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)

        End If

      End If

      Do row2 = 1, region_density_y(row)

        element = element + 1
        element_node(element, 1) = node_sw
        element_node(element, 2) = node_sw + inc1 + inc2
        element_node(element, 3) = node_sw + 2
        element_node(element, 4) = node_sw + inc1
        element_node(element, 5) = node_sw + inc1 + 1
        element_node(element, 6) = node_sw + 1

        element = element + 1
        element_node(element, 1) = node_sw + inc1 + inc2 + 2
        element_node(element, 2) = node_sw + 2
        element_node(element, 3) = node_sw + inc1 + inc2
        element_node(element, 4) = node_sw + inc1 + 2
        element_node(element, 5) = node_sw + inc1 + 1
        element_node(element, 6) = node_sw + inc1 + inc2 + 1

        node_sw = node_sw + 2

      End Do

    End If

  End Do
End Do

End Do

Return
End Subroutine hcell_element_node
Subroutine hcell_node_count(region_density_x, region_density_y, node_num)

!*********************************************************************72
!
!c HCELL_NODE_NUM determines the number of nodes in the region.
!
Implicit None

Integer node_num
Integer region_density_x(3)
Integer region_density_y(3)
!
! Count the nodes.
!
node_num = (2region_density_x(1)+1)(2*region_density_y(1)+1)&

• (2region_density_x(1)+1)(2*region_density_y(3)+1) &
• (2region_density_x(2)-1)(2*region_density_y(1)+1)&
• (2region_density_x(2)+1)(2*region_density_y(2)-1)&
• (2region_density_x(2)-1)(2*region_density_y(3)+1)&
• (2region_density_x(3)+1)(2*region_density_y(1)+1) &
• (2region_density_x(3)+1)(2*region_density_y(3)+1)

Write (, '(a)') ' '
Write (, '(a)') 'HCELL_NODE_COUNT:'
Write (*, '(a,i6)') ' Number of nodes = ', node_num

Return
End Subroutine hcell_node_count
Subroutine hcell_node_xy(region_density_x, region_density_y,&
node_num, region_x, region_y, node_x, node_y)

!*********************************************************************72
!
!c HCELL_NODE_XY assigns coordinates to each node.
!
Implicit None

Integer node_num

Integer col
Integer i
Integer j
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Integer region_density_x(3)
Integer region_density_y(3)
Double Precision region_x(4)
Double Precision region_y(4)
Integer row

node = 0
j = 0
!
! Working in column 1, except for the extreme right.
!
Do col = 1, 2*region_density_x(1)

j = j + 1
i = 0

!
! Working in row 1.
!
Do row = 1, 2*region_density_y(1) + 1

i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(1)+1-col)*region_x(1)+&
  dble(-1+col)*region_x(2))/dble(2*region_density_x(1))

  node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
  dble(-1+row)*region_y(2))/dble(2*region_density_y(1))

End Do

!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1

Do row = 1, 2*region_density_y(3) + 1

  i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(1)+1-col)*region_x(1)+&
  dble(-1+col)*region_x(2))/dble(2*region_density_x(1))

  node_y(node) = (dble(2*region_density_y(3)+1-row)*region_y(3)+&
  dble(-1+row)*region_y(4))/dble(2*region_density_y(3))

End Do

End Do
!
! Working in column 2, including extreme left and right.
!
Do col = 1, 2*region_density_x(2) + 1
!
! Working in row 1.
!
j = j + 1
i = 0

Do row = 1, 2*region_density_y(1) + 1

  i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
  dble(-1+col)*region_x(3))/dble(2*region_density_x(2))

  node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
  dble(-1+row)*region_y(2))/dble(2*region_density_y(1))

End Do

!
! Working in row 2.
!
Do row = 2, 2*region_density_y(2)

i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
  dble(-1+col)*region_x(3))/dble(2*region_density_x(2))

  node_y(node) = (dble(2*region_density_y(2)+1-row)*region_y(2)+&
  dble(-1+row)*region_y(3))/dble(2*region_density_y(2))

End Do

!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1

i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
  dble(-1+col)*region_x(3))/dble(2*region_density_x(2))

  node_y(node) = (dble(2*region_density_y(3)+1-row)*region_y(3)+&
  dble(-1+row)*region_y(4))/dble(2*region_density_y(3))

End Do

End Do
!
! Working in column 3, except for extreme left.
!
Do col = 2, 2*region_density_x(3) + 1
!
! Working in row 1.
!
j = j + 1
i = 0

Do row = 1, 2*region_density_y(1) + 1

  i = i + 1
  node = node + 1

  node_x(node) = (dble(2*region_density_x(3)+1-col)*region_x(3)+&
  dble(-1+col)*region_x(4))/dble(2*region_density_x(3))

  node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
  dble(-1+row)*region_y(2))/dble(2*region_density_y(1))

End Do

i = i + 2*region_density_y(2) - 1

!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1

node_x(node) = (dble(2*region_density_x(3)+1-col)&
  *region_x(3)+dble(-1+col)*region_x(4))/dble(2*region_density_x(3))

  node_y(node) = (dble(2*region_density_y(3)+1-row)&
  *region_y(3)+dble(-1+row)*region_y(4))/dble(2*region_density_y(3))
End Do

End Do

Return
End Subroutine hcell_node_xy
Integer Function idamax(n, dx, incx)

!*********************************************************************72
!
!c IDAMAX finds the vector element of largest magnitude.
!
Double Precision dx(1), dmax
Integer i, incx, ix, n

idamax = 0
If (n<1) Return
idamax = 1
If (n==1) Return
If (incx==1) Goto 20
!
! CODE FOR INCREMENT NOT EQUAL TO 1
!
ix = 1
dmax = dabs(dx(1))
ix = ix + incx
Do i = 2, n
If (dabs(dx(ix))<=dmax) Goto 5
idamax = i
dmax = dabs(dx(ix))
5 ix = ix + incx
End Do
Return
!
! CODE FOR INCREMENT EQUAL TO 1
!
20 dmax = dabs(dx(1))
Do i = 2, n
If (dabs(dx(i))<=dmax) Goto 30
idamax = i
dmax = dabs(dx(i))
30 End Do
Return
End Function idamax
Subroutine node_eps(node_file_name, np, node_mask, xc, yc, title)

!*********************************************************************72
!
!c NODE_EPS creates an EPS file containing an image of the nodes.
!
Implicit None

Integer np

Double Precision ave_x
Double Precision ave_y
Integer circle_size
Double Precision dif
Integer eps_unit
Integer eps_x
Integer eps_y
Character (Len=*) node_file_name
Integer i
Integer ios
Integer ip1
Integer j
Integer k
Integer node
Logical node_mask(np)
Integer node_unmasked_num
Double Precision node_x_max
Double Precision node_x_min
Double Precision node_y_max
Double Precision node_y_min
Double Precision r8_huge
Double Precision scale
Character (Len=40) string
Character (Len=*) title
Double Precision xc(np)
Double Precision yc(np)
!
! Determine the range and number of the unmasked nodes.
!
node_x_min = r8_huge()
node_x_max = -r8_huge()
node_y_min = r8_huge()
node_y_max = -r8_huge()

node_unmasked_num = 0

Do node = 1, np
If (node_mask(node)) Then
node_unmasked_num = node_unmasked_num + 1
node_x_min = min(node_x_min, xc(node))
node_x_max = max(node_x_max, xc(node))
node_y_min = min(node_y_min, yc(node))
node_y_max = max(node_y_max, yc(node))
End If
End Do

If (node_unmasked_num<=200) Then
circle_size = 3
Else If (node_unmasked_num<=800) Then
circle_size = 2
Else
circle_size = 1
End If

If (node_y_max-node_y_min<node_x_max-node_x_min) Then
scale = node_x_max - node_x_min
dif = (node_x_max-node_x_min) - (node_y_max-node_y_min)
node_y_max = node_y_max + 0.5dif
node_y_min = node_y_min - 0.5dif
Else
scale = node_y_max - node_y_min
dif = (node_y_max-node_y_min) - (node_x_max-node_x_min)
node_x_max = node_x_max + 0.5dif
node_x_min = node_x_min - 0.5dif
End If

eps_unit = 1

Open (Unit=eps_unit, File=node_file_name, Status='unknown')

Write (eps_unit, '(a)') '%cPS-Adobe-3.0 EPSF-3.0'
Write (eps_unit, '(a)') '%%Creator: node_eps(hcell.f)'
Write (eps_unit, '(a,a)') '%% : ', node_file_name
Write (eps_unit, '(a)') '%%Pages: 1'
Write (eps_unit, '(a)') '%%BoundingBox: 36 36 576 756'
Write (eps_unit, '(a)') '%%Document-Fonts: Times-Roman'
Write (eps_unit, '(a)') '%%LanguageLevel: 1'
Write (eps_unit, '(a)') '%%EndComments'
Write (eps_unit, '(a)') '%%BeginProlog'
Write (eps_unit, '(a)') '/inch {72 mul} def'
Write (eps_unit, '(a)') '%%EndProlog'
Write (eps_unit, '(a)') '%%Page: 1 1'
Write (eps_unit, '(a)') 'save'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Set RGB line color.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.9000 0.9000 0.9000 setrgbcolor'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Draw a gray border around the page.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') 'newpath'
Write (eps_unit, '(a)') ' 36 126 moveto'
Write (eps_unit, '(a)') ' 576 126 lineto'
Write (eps_unit, '(a)') ' 576 666 lineto'
Write (eps_unit, '(a)') ' 36 666 lineto'
Write (eps_unit, '(a)') ' 36 126 lineto'
Write (eps_unit, '(a)') 'stroke'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Set RGB line color.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'

Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Label the plot:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'
Write (eps_unit, '(a)') '/Times-Roman findfont 0.50 inch scalefont setfont'
Write (eps_unit, '(a)') ' 36 666 moveto'
Write (eps_unit, '(a,a,a)') '(',title, ') show'

Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Define a clipping polygon'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 36 126 moveto'
Write (eps_unit, '(a)') ' 576 126 lineto'
Write (eps_unit, '(a)') ' 576 666 lineto'
Write (eps_unit, '(a)') ' 36 666 lineto'
Write (eps_unit, '(a)') ' 36 126 lineto'
Write (eps_unit, '(a)') 'clip newpath'

Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Draw filled dots at each node:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 1.0000 setrgbcolor'

Do node = 1, np

If (node_mask(node)) Then

  eps_x = int((node_x_max-xc(node))*61.0+(+xc(node)-node_x_min)*551.0)/scale

  eps_y = int((node_y_max-yc(node))*151.0+(yc(node)-node_y_min)*641.0)/scale

  Write (eps_unit, '(a,i4,2x,i4,2x,i4,a)') 'newpath  ',&
  eps_x, eps_y, circle_size, ' 0 360 arc closepath fill'

End If

End Do
!
! Label each node, but only if there aren't too many of themc
!
If (node_unmasked_num<=200) Then

Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '%  Label the nodes:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'
Write (eps_unit, '(a)') '/Times-Roman findfont 0.20 inch scalefont setfont'

Do node = 1, np

  If (node_mask(node)) Then

    eps_x = int((node_x_max-xc(node))*61.0+(+xc(node)-node_x_min)*551.0)/scale

    eps_y = int((node_y_max-yc(node))*151.0+(yc(node)-node_y_min)*641.0)/scale

    Write (eps_unit, '(i4,2x,i4,a,i4,a)') eps_x, eps_y + 5,&
    ' moveto (', node, ') show'

  End If

End Do

End If

Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') 'restore showpage'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% End of page'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '%%Trailer'
Write (eps_unit, '(a)') '%%EOF'

Close (Unit=eps_unit)

Write (, '(a)') ' '
Write (, '(a)') 'NODE_EPS:'
Write (, '(a)') ' An encapsulated Post file was created'
Write (, '(a)') ' containing an image of the nodes.'
Write (, '(a)') ' The file name is'
Write (, '(4x,a)') node_file_name

Return
End Subroutine node_eps
Subroutine nstoke(xc, yc, area, xm, ym, a, f, g, uold, reynld,&
tolns, xlngth, ylngth, node, indx, ipivot, mrow1, nlband, &
nuband, nband, nrow1, ncol1, nelemn, np, nnodes, nuk, nquad,&
neqn, nsteps, nsim, maxnd, maxel, rdel, alpha)

!*********************************************************************72
!
!c NSTOKE solves the Navier-Stokes equations using Taylor-Hood elements.
!
Implicit Double Precision (A-H, O-Z)

Double Precision a(mrow1, )
Double Precision area()
Double Precision f()
Double Precision g()
Integer indx(maxnd, 3)
Integer ipivot()
Integer it
Integer iter
Integer node(maxel, )
Double Precision reynld
Double Precision un(2)
Double Precision unx(2)
Double Precision uny(2)
Double Precision uold()
Double Precision uold_qp
Double Precision visc
Double Precision vold_qp
Double Precision x
Double Precision xc()
Double Precision xm(maxel, )
Double Precision y
Double Precision yc()
Double Precision ym(maxel, *)
!
! zero arrays
! g array contains the previous iterate
! f array contains the right hand side initially and then the current
! iterate is overwritten on f
!
visc = 1.D0/reynld
!
! matrix assembly triangle by triangle
! nsim is the number of simple iterations performed
!
csim = 0.D0
Do iter = 1, nsteps

niter = iter
If (iter>nsim) csim = 1.D0

Do i = 1, nrow1
  Do j = 1, ncol1
    a(i, j) = 0.D0
  End Do
End Do

Do it = 1, nelemn

  arr = area(it)/3.D0

  Do iquad = 1, nquad

    y = ym(it, iquad)
    x = xm(it, iquad)
    Call trans(it, x, y, det, xix, xiy, etax, etay, xc, yc, node, maxel)
    ar = arr*det

    Do kk = 1, 2
      un(kk) = 0.D0
      uny(kk) = 0.D0
      unx(kk) = 0.D0
    End Do
    uold_qp = 0.0D0
    vold_qp = 0.0D0

    Do iq = 1, nnodes

      Call refqbf(x, y, iq, bb, tbx, tby)
      bx = tbx*xix + tby*etax
      by = tbx*xiy + tby*etay
      ip = node(it, iq)

      Do iuk = 1, 2

        iun = indx(ip, iuk)

        If (0<iun) Then
          un(iuk) = un(iuk) + bb*g(iun)
          unx(iuk) = unx(iuk) + bx*g(iun)
          uny(iuk) = uny(iuk) + by*g(iun)
          If (iuk==1) uold_qp = uold_qp + bb*uold(iun)
          If (iuk==2) vold_qp = vold_qp + bb*uold(iun)
        Else If (iun<0) Then
          ubc = alpha*ubdry(iuk, ip, xc, yc)
          un(iuk) = un(iuk) + bb*ubc
          unx(iuk) = unx(iuk) + bx*ubc
          uny(iuk) = uny(iuk) + by*ubc
          If (iuk==1) uold_qp = uold_qp + bb*ubc
          If (iuk==2) vold_qp = vold_qp + bb*ubc
        End If

      End Do

    End Do

    Do iq = 1, nnodes
      ip = node(it, iq)
      Call refqbf(x, y, iq, bb, tbx, tby)
      bx = tbx*xix + tby*etax
      by = tbx*xiy + tby*etay
      bbl = refbsp(x, y, iq)
      Do iuk = 1, nuk
        i = indx(ip, iuk)
        If (i<=0) Goto 210
        If (iuk==1) f(i) = f(i) + csim*((un(1)*unx(1)+un(2)*uny(1))&
        *bb)*ar + rdel*uold_qp*bb*ar
        If (iuk==2) f(i) = f(i) + csim*((un(1)*unx(2)+un(2)*uny(2))&
        *bb)*ar + rdel*vold_qp*bb*ar
        Do iqq = 1, nnodes
          ipp = node(it, iqq)
          Call refqbf(x, y, iqq, bbb, tbbx, tbby)
          bbx = tbbx*xix + tbby*etax
          bby = tbbx*xiy + tbby*etay
          bbbl = refbsp(x, y, iqq)
          Do iukk = 1, nuk
            j = indx(ipp, iukk)
            If (j==0) Goto 190
            aij = 0.D0
            If (i==neqn) Goto 190
            If (iuk==1) Then
              If (iukk==1) aij = visc*(by*bby+bx*bbx)&
              + (bbb*unx(1)*bb)*csim + bb*bbx*un(1) +&
              bb*bby*un(2) + rdel*(bb*bbb)
              If (iukk==2) aij = csim*(bb*bbb*uny(1))
              If (iukk==3) aij = -bx*bbbl
            Else If (iuk==2) Then
              If (iukk==2) aij = (visc*(by*bby+bx*bbx)+&
              (bb*bbb*uny(2))*csim+bb*bby*un(2)+bb*bbx*un(1)) + rdel*(bb*bbb)
              If (iukk==1) aij = csim*(bb*bbb*unx(2))
              If (iukk==3) aij = -by*bbbl
            Else
              If (iukk==1) aij = bbx*bbl
              If (iukk==2) aij = bby*bbl
            End If
            If (j<0) Goto 185
            iuse = i - j + nband
            a(iuse, j) = a(iuse, j) + aij*ar
            Goto 190

!
! add terms to rhs for inhomogeneous boundary condition
!
185 Continue
f(i) = f(i) - aralphaubdry(iukk, ipp, xc, yc)*aij
190 End Do
End Do
210 End Do
End Do
End Do
End Do
f(neqn) = 0.D0
Do j = neqn - nlband, neqn - 1
i = neqn - j + nband
a(i, j) = 0.D0
End Do
a(nband, neqn) = 1.D0
!
! solve system
!
job = 0
Call dgbfa(a, mrow1, neqn, nlband, nuband, ipivot, info)
Call dgbsl(a, mrow1, neqn, nlband, nuband, ipivot, f, job)
!
! check for convergence
!
diff = 0.D0
Do i = 1, neqn
diff = diff + (g(i)-f(i))**2
End Do
diff = sqrt(diff)

Write (*, *) '  Iteration ', iter, '  Difference is ', diff

If (diff<=tolns) Then
  Return
End If

Do i = 1, neqn
  g(i) = f(i)
  f(i) = 0.D0
End Do

End Do

Return
End Subroutine nstoke