! This file was automatically generated by SWIG (http://www.swig.org).
! Version 4.0.0
!
! Do not make changes to this file unless you know what you are doing--modify
! the SWIG interface file instead.

! ---------------------------------------------------------------
! Programmer(s): Auto-generated by swig.
! ---------------------------------------------------------------
! SUNDIALS Copyright Start
! Copyright (c) 2002-2022, Lawrence Livermore National Security
! and Southern Methodist University.
! All rights reserved.
!
! See the top-level LICENSE and NOTICE files for details.
!
! SPDX-License-Identifier: BSD-3-Clause
! SUNDIALS Copyright End
! ---------------------------------------------------------------

module fnvector_pthreads_mod
 use, intrinsic :: ISO_C_BINDING
 use fsundials_nvector_mod
 use fsundials_context_mod
 use fsundials_types_mod
 implicit none
 private

 ! DECLARATION CONSTRUCTS
 public :: FN_VNew_Pthreads
 public :: FN_VNewEmpty_Pthreads
 public :: FN_VMake_Pthreads
 public :: FN_VGetLength_Pthreads
 public :: FN_VPrint_Pthreads
 public :: FN_VPrintFile_Pthreads
 public :: FN_VGetVectorID_Pthreads
 public :: FN_VCloneEmpty_Pthreads
 public :: FN_VClone_Pthreads
 public :: FN_VDestroy_Pthreads
 public :: FN_VSpace_Pthreads
 public :: FN_VGetArrayPointer_Pthreads
 public :: FN_VSetArrayPointer_Pthreads
 public :: FN_VLinearSum_Pthreads
 public :: FN_VConst_Pthreads
 public :: FN_VProd_Pthreads
 public :: FN_VDiv_Pthreads
 public :: FN_VScale_Pthreads
 public :: FN_VAbs_Pthreads
 public :: FN_VInv_Pthreads
 public :: FN_VAddConst_Pthreads
 public :: FN_VDotProd_Pthreads
 public :: FN_VMaxNorm_Pthreads
 public :: FN_VWrmsNorm_Pthreads
 public :: FN_VWrmsNormMask_Pthreads
 public :: FN_VMin_Pthreads
 public :: FN_VWL2Norm_Pthreads
 public :: FN_VL1Norm_Pthreads
 public :: FN_VCompare_Pthreads
 public :: FN_VInvTest_Pthreads
 public :: FN_VConstrMask_Pthreads
 public :: FN_VMinQuotient_Pthreads
 public :: FN_VLinearCombination_Pthreads
 public :: FN_VScaleAddMulti_Pthreads
 public :: FN_VDotProdMulti_Pthreads
 public :: FN_VLinearSumVectorArray_Pthreads
 public :: FN_VScaleVectorArray_Pthreads
 public :: FN_VConstVectorArray_Pthreads
 public :: FN_VWrmsNormVectorArray_Pthreads
 public :: FN_VWrmsNormMaskVectorArray_Pthreads
 public :: FN_VWSqrSumLocal_Pthreads
 public :: FN_VWSqrSumMaskLocal_Pthreads
 public :: FN_VBufSize_Pthreads
 public :: FN_VBufPack_Pthreads
 public :: FN_VBufUnpack_Pthreads
 public :: FN_VEnableFusedOps_Pthreads
 public :: FN_VEnableLinearCombination_Pthreads
 public :: FN_VEnableScaleAddMulti_Pthreads
 public :: FN_VEnableDotProdMulti_Pthreads
 public :: FN_VEnableLinearSumVectorArray_Pthreads
 public :: FN_VEnableScaleVectorArray_Pthreads
 public :: FN_VEnableConstVectorArray_Pthreads
 public :: FN_VEnableWrmsNormVectorArray_Pthreads
 public :: FN_VEnableWrmsNormMaskVectorArray_Pthreads
 public :: FN_VCloneVectorArray_Pthreads
 public :: FN_VCloneVectorArrayEmpty_Pthreads
 public :: FN_VDestroyVectorArray_Pthreads

! WRAPPER DECLARATIONS
interface
function swigc_FN_VNew_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VNew_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), intent(in) :: farg1
integer(C_INT), intent(in) :: farg2
type(C_PTR), value :: farg3
type(C_PTR) :: fresult
end function

function swigc_FN_VNewEmpty_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VNewEmpty_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), intent(in) :: farg1
integer(C_INT), intent(in) :: farg2
type(C_PTR), value :: farg3
type(C_PTR) :: fresult
end function

function swigc_FN_VMake_Pthreads(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FN_VMake_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T), intent(in) :: farg1
integer(C_INT), intent(in) :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
type(C_PTR) :: fresult
end function

function swigc_FN_VGetLength_Pthreads(farg1) &
bind(C, name="_wrap_FN_VGetLength_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT64_T) :: fresult
end function

subroutine swigc_FN_VPrint_Pthreads(farg1) &
bind(C, name="_wrap_FN_VPrint_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
end subroutine

subroutine swigc_FN_VPrintFile_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VPrintFile_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
end subroutine

function swigc_FN_VGetVectorID_Pthreads(farg1) &
bind(C, name="_wrap_FN_VGetVectorID_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT) :: fresult
end function

function swigc_FN_VCloneEmpty_Pthreads(farg1) &
bind(C, name="_wrap_FN_VCloneEmpty_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

function swigc_FN_VClone_Pthreads(farg1) &
bind(C, name="_wrap_FN_VClone_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

subroutine swigc_FN_VDestroy_Pthreads(farg1) &
bind(C, name="_wrap_FN_VDestroy_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
end subroutine

subroutine swigc_FN_VSpace_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VSpace_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
end subroutine

function swigc_FN_VGetArrayPointer_Pthreads(farg1) &
bind(C, name="_wrap_FN_VGetArrayPointer_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR) :: fresult
end function

subroutine swigc_FN_VSetArrayPointer_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VSetArrayPointer_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
end subroutine

subroutine swigc_FN_VLinearSum_Pthreads(farg1, farg2, farg3, farg4, farg5) &
bind(C, name="_wrap_FN_VLinearSum_Pthreads")
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE), intent(in) :: farg3
type(C_PTR), value :: farg4
type(C_PTR), value :: farg5
end subroutine

subroutine swigc_FN_VConst_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VConst_Pthreads")
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
end subroutine

subroutine swigc_FN_VProd_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VProd_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
end subroutine

subroutine swigc_FN_VDiv_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VDiv_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
end subroutine

subroutine swigc_FN_VScale_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VScale_Pthreads")
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
end subroutine

subroutine swigc_FN_VAbs_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VAbs_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
end subroutine

subroutine swigc_FN_VInv_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VInv_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
end subroutine

subroutine swigc_FN_VAddConst_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VAddConst_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE), intent(in) :: farg2
type(C_PTR), value :: farg3
end subroutine

function swigc_FN_VDotProd_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VDotProd_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VMaxNorm_Pthreads(farg1) &
bind(C, name="_wrap_FN_VMaxNorm_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VWrmsNorm_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VWrmsNorm_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VWrmsNormMask_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VWrmsNormMask_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VMin_Pthreads(farg1) &
bind(C, name="_wrap_FN_VMin_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VWL2Norm_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VWL2Norm_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VL1Norm_Pthreads(farg1) &
bind(C, name="_wrap_FN_VL1Norm_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
real(C_DOUBLE) :: fresult
end function

subroutine swigc_FN_VCompare_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VCompare_Pthreads")
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
end subroutine

function swigc_FN_VInvTest_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VInvTest_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VConstrMask_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VConstrMask_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
integer(C_INT) :: fresult
end function

function swigc_FN_VMinQuotient_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VMinQuotient_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VLinearCombination_Pthreads(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FN_VLinearCombination_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
integer(C_INT) :: fresult
end function

function swigc_FN_VScaleAddMulti_Pthreads(farg1, farg2, farg3, farg4, farg5) &
bind(C, name="_wrap_FN_VScaleAddMulti_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
type(C_PTR), value :: farg5
integer(C_INT) :: fresult
end function

function swigc_FN_VDotProdMulti_Pthreads(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FN_VDotProdMulti_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
integer(C_INT) :: fresult
end function

function swigc_FN_VLinearSumVectorArray_Pthreads(farg1, farg2, farg3, farg4, farg5, farg6) &
bind(C, name="_wrap_FN_VLinearSumVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
real(C_DOUBLE), intent(in) :: farg2
type(C_PTR), value :: farg3
real(C_DOUBLE), intent(in) :: farg4
type(C_PTR), value :: farg5
type(C_PTR), value :: farg6
integer(C_INT) :: fresult
end function

function swigc_FN_VScaleVectorArray_Pthreads(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FN_VScaleVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
integer(C_INT) :: fresult
end function

function swigc_FN_VConstVectorArray_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VConstVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
real(C_DOUBLE), intent(in) :: farg2
type(C_PTR), value :: farg3
integer(C_INT) :: fresult
end function

function swigc_FN_VWrmsNormVectorArray_Pthreads(farg1, farg2, farg3, farg4) &
bind(C, name="_wrap_FN_VWrmsNormVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
integer(C_INT) :: fresult
end function

function swigc_FN_VWrmsNormMaskVectorArray_Pthreads(farg1, farg2, farg3, farg4, farg5) &
bind(C, name="_wrap_FN_VWrmsNormMaskVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
type(C_PTR), value :: farg4
type(C_PTR), value :: farg5
integer(C_INT) :: fresult
end function

function swigc_FN_VWSqrSumLocal_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VWSqrSumLocal_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VWSqrSumMaskLocal_Pthreads(farg1, farg2, farg3) &
bind(C, name="_wrap_FN_VWSqrSumMaskLocal_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
type(C_PTR), value :: farg3
real(C_DOUBLE) :: fresult
end function

function swigc_FN_VBufSize_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VBufSize_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VBufPack_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VBufPack_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VBufUnpack_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VBufUnpack_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
type(C_PTR), value :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableFusedOps_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableFusedOps_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableLinearCombination_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableLinearCombination_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableScaleAddMulti_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableScaleAddMulti_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableDotProdMulti_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableDotProdMulti_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableLinearSumVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableLinearSumVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableScaleVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableScaleVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableConstVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableConstVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableWrmsNormVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableWrmsNormVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VEnableWrmsNormMaskVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VEnableWrmsNormMaskVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
integer(C_INT) :: fresult
end function

function swigc_FN_VCloneVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VCloneVectorArray_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR) :: fresult
end function

function swigc_FN_VCloneVectorArrayEmpty_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VCloneVectorArrayEmpty_Pthreads") &
result(fresult)
use, intrinsic :: ISO_C_BINDING
integer(C_INT), intent(in) :: farg1
type(C_PTR), value :: farg2
type(C_PTR) :: fresult
end function

subroutine swigc_FN_VDestroyVectorArray_Pthreads(farg1, farg2) &
bind(C, name="_wrap_FN_VDestroyVectorArray_Pthreads")
use, intrinsic :: ISO_C_BINDING
type(C_PTR), value :: farg1
integer(C_INT), intent(in) :: farg2
end subroutine

end interface

contains
 ! MODULE SUBPROGRAMS
function FN_VNew_Pthreads(vec_length, n_threads, sunctx) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), pointer :: swig_result
integer(C_INT64_T), intent(in) :: vec_length
integer(C_INT), intent(in) :: n_threads
type(C_PTR) :: sunctx
type(C_PTR) :: fresult
integer(C_INT64_T) :: farg1
integer(C_INT) :: farg2
type(C_PTR) :: farg3

farg1 = vec_length
farg2 = n_threads
farg3 = sunctx
fresult = swigc_FN_VNew_Pthreads(farg1, farg2, farg3)
call c_f_pointer(fresult, swig_result)
end function

function FN_VNewEmpty_Pthreads(vec_length, n_threads, sunctx) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), pointer :: swig_result
integer(C_INT64_T), intent(in) :: vec_length
integer(C_INT), intent(in) :: n_threads
type(C_PTR) :: sunctx
type(C_PTR) :: fresult
integer(C_INT64_T) :: farg1
integer(C_INT) :: farg2
type(C_PTR) :: farg3

farg1 = vec_length
farg2 = n_threads
farg3 = sunctx
fresult = swigc_FN_VNewEmpty_Pthreads(farg1, farg2, farg3)
call c_f_pointer(fresult, swig_result)
end function

function FN_VMake_Pthreads(vec_length, n_threads, v_data, sunctx) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), pointer :: swig_result
integer(C_INT64_T), intent(in) :: vec_length
integer(C_INT), intent(in) :: n_threads
real(C_DOUBLE), dimension(*), target, intent(inout) :: v_data
type(C_PTR) :: sunctx
type(C_PTR) :: fresult
integer(C_INT64_T) :: farg1
integer(C_INT) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4

farg1 = vec_length
farg2 = n_threads
farg3 = c_loc(v_data(1))
farg4 = sunctx
fresult = swigc_FN_VMake_Pthreads(farg1, farg2, farg3, farg4)
call c_f_pointer(fresult, swig_result)
end function

function FN_VGetLength_Pthreads(v) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT64_T) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT64_T) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(v)
fresult = swigc_FN_VGetLength_Pthreads(farg1)
swig_result = fresult
end function

subroutine FN_VPrint_Pthreads(v)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: v
type(C_PTR) :: farg1

farg1 = c_loc(v)
call swigc_FN_VPrint_Pthreads(farg1)
end subroutine

subroutine FN_VPrintFile_Pthreads(v, outfile)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: v
type(C_PTR) :: outfile
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(v)
farg2 = outfile
call swigc_FN_VPrintFile_Pthreads(farg1, farg2)
end subroutine

function FN_VGetVectorID_Pthreads(v) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(N_Vector_ID) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(v)
fresult = swigc_FN_VGetVectorID_Pthreads(farg1)
swig_result = fresult
end function

function FN_VCloneEmpty_Pthreads(w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), pointer :: swig_result
type(N_Vector), target, intent(inout) :: w
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(w)
fresult = swigc_FN_VCloneEmpty_Pthreads(farg1)
call c_f_pointer(fresult, swig_result)
end function

function FN_VClone_Pthreads(w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), pointer :: swig_result
type(N_Vector), target, intent(inout) :: w
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(w)
fresult = swigc_FN_VClone_Pthreads(farg1)
call c_f_pointer(fresult, swig_result)
end function

subroutine FN_VDestroy_Pthreads(v)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: v
type(C_PTR) :: farg1

farg1 = c_loc(v)
call swigc_FN_VDestroy_Pthreads(farg1)
end subroutine

subroutine FN_VSpace_Pthreads(v, lrw, liw)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: v
integer(C_INT64_T), dimension(*), target, intent(inout) :: lrw
integer(C_INT64_T), dimension(*), target, intent(inout) :: liw
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(v)
farg2 = c_loc(lrw(1))
farg3 = c_loc(liw(1))
call swigc_FN_VSpace_Pthreads(farg1, farg2, farg3)
end subroutine

function FN_VGetArrayPointer_Pthreads(v) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), dimension(:), pointer :: swig_result
type(N_Vector), target, intent(inout) :: v
type(C_PTR) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(v)
fresult = swigc_FN_VGetArrayPointer_Pthreads(farg1)
call c_f_pointer(fresult, swig_result, [1])
end function

subroutine FN_VSetArrayPointer_Pthreads(v_data, v)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), dimension(*), target, intent(inout) :: v_data
type(N_Vector), target, intent(inout) :: v
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(v_data(1))
farg2 = c_loc(v)
call swigc_FN_VSetArrayPointer_Pthreads(farg1, farg2)
end subroutine

subroutine FN_VLinearSum_Pthreads(a, x, b, y, z)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: a
type(N_Vector), target, intent(inout) :: x
real(C_DOUBLE), intent(in) :: b
type(N_Vector), target, intent(inout) :: y
type(N_Vector), target, intent(inout) :: z
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2
real(C_DOUBLE) :: farg3
type(C_PTR) :: farg4
type(C_PTR) :: farg5

farg1 = a
farg2 = c_loc(x)
farg3 = b
farg4 = c_loc(y)
farg5 = c_loc(z)
call swigc_FN_VLinearSum_Pthreads(farg1, farg2, farg3, farg4, farg5)
end subroutine

subroutine FN_VConst_Pthreads(c, z)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: c
type(N_Vector), target, intent(inout) :: z
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2

farg1 = c
farg2 = c_loc(z)
call swigc_FN_VConst_Pthreads(farg1, farg2)
end subroutine

subroutine FN_VProd_Pthreads(x, y, z)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: y
type(N_Vector), target, intent(inout) :: z
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(x)
farg2 = c_loc(y)
farg3 = c_loc(z)
call swigc_FN_VProd_Pthreads(farg1, farg2, farg3)
end subroutine

subroutine FN_VDiv_Pthreads(x, y, z)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: y
type(N_Vector), target, intent(inout) :: z
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(x)
farg2 = c_loc(y)
farg3 = c_loc(z)
call swigc_FN_VDiv_Pthreads(farg1, farg2, farg3)
end subroutine

subroutine FN_VScale_Pthreads(c, x, z)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: c
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: z
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c
farg2 = c_loc(x)
farg3 = c_loc(z)
call swigc_FN_VScale_Pthreads(farg1, farg2, farg3)
end subroutine

subroutine FN_VAbs_Pthreads(x, z)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: z
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(z)
call swigc_FN_VAbs_Pthreads(farg1, farg2)
end subroutine

subroutine FN_VInv_Pthreads(x, z)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: z
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(z)
call swigc_FN_VInv_Pthreads(farg1, farg2)
end subroutine

subroutine FN_VAddConst_Pthreads(x, b, z)
use, intrinsic :: ISO_C_BINDING
type(N_Vector), target, intent(inout) :: x
real(C_DOUBLE), intent(in) :: b
type(N_Vector), target, intent(inout) :: z
type(C_PTR) :: farg1
real(C_DOUBLE) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(x)
farg2 = b
farg3 = c_loc(z)
call swigc_FN_VAddConst_Pthreads(farg1, farg2, farg3)
end subroutine

function FN_VDotProd_Pthreads(x, y) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: y
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(y)
fresult = swigc_FN_VDotProd_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VMaxNorm_Pthreads(x) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(x)
fresult = swigc_FN_VMaxNorm_Pthreads(farg1)
swig_result = fresult
end function

function FN_VWrmsNorm_Pthreads(x, w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: w
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(w)
fresult = swigc_FN_VWrmsNorm_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VWrmsNormMask_Pthreads(x, w, id) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: w
type(N_Vector), target, intent(inout) :: id
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(x)
farg2 = c_loc(w)
farg3 = c_loc(id)
fresult = swigc_FN_VWrmsNormMask_Pthreads(farg1, farg2, farg3)
swig_result = fresult
end function

function FN_VMin_Pthreads(x) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(x)
fresult = swigc_FN_VMin_Pthreads(farg1)
swig_result = fresult
end function

function FN_VWL2Norm_Pthreads(x, w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: w
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(w)
fresult = swigc_FN_VWL2Norm_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VL1Norm_Pthreads(x) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1

farg1 = c_loc(x)
fresult = swigc_FN_VL1Norm_Pthreads(farg1)
swig_result = fresult
end function

subroutine FN_VCompare_Pthreads(c, x, z)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE), intent(in) :: c
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: z
real(C_DOUBLE) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c
farg2 = c_loc(x)
farg3 = c_loc(z)
call swigc_FN_VCompare_Pthreads(farg1, farg2, farg3)
end subroutine

function FN_VInvTest_Pthreads(x, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: z
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(z)
fresult = swigc_FN_VInvTest_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VConstrMask_Pthreads(c, x, m) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: c
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: m
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(c)
farg2 = c_loc(x)
farg3 = c_loc(m)
fresult = swigc_FN_VConstrMask_Pthreads(farg1, farg2, farg3)
swig_result = fresult
end function

function FN_VMinQuotient_Pthreads(num, denom) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: num
type(N_Vector), target, intent(inout) :: denom
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(num)
farg2 = c_loc(denom)
fresult = swigc_FN_VMinQuotient_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VLinearCombination_Pthreads(nvec, c, x, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
real(C_DOUBLE), dimension(*), target, intent(inout) :: c
type(C_PTR) :: x
type(N_Vector), target, intent(inout) :: z
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4

farg1 = nvec
farg2 = c_loc(c(1))
farg3 = x
farg4 = c_loc(z)
fresult = swigc_FN_VLinearCombination_Pthreads(farg1, farg2, farg3, farg4)
swig_result = fresult
end function

function FN_VScaleAddMulti_Pthreads(nvec, a, x, y, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
real(C_DOUBLE), dimension(*), target, intent(inout) :: a
type(N_Vector), target, intent(inout) :: x
type(C_PTR) :: y
type(C_PTR) :: z
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4
type(C_PTR) :: farg5

farg1 = nvec
farg2 = c_loc(a(1))
farg3 = c_loc(x)
farg4 = y
farg5 = z
fresult = swigc_FN_VScaleAddMulti_Pthreads(farg1, farg2, farg3, farg4, farg5)
swig_result = fresult
end function

function FN_VDotProdMulti_Pthreads(nvec, x, y, dotprods) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
type(N_Vector), target, intent(inout) :: x
type(C_PTR) :: y
real(C_DOUBLE), dimension(*), target, intent(inout) :: dotprods
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4

farg1 = nvec
farg2 = c_loc(x)
farg3 = y
farg4 = c_loc(dotprods(1))
fresult = swigc_FN_VDotProdMulti_Pthreads(farg1, farg2, farg3, farg4)
swig_result = fresult
end function

function FN_VLinearSumVectorArray_Pthreads(nvec, a, x, b, y, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
real(C_DOUBLE), intent(in) :: a
type(C_PTR) :: x
real(C_DOUBLE), intent(in) :: b
type(C_PTR) :: y
type(C_PTR) :: z
integer(C_INT) :: fresult
integer(C_INT) :: farg1
real(C_DOUBLE) :: farg2
type(C_PTR) :: farg3
real(C_DOUBLE) :: farg4
type(C_PTR) :: farg5
type(C_PTR) :: farg6

farg1 = nvec
farg2 = a
farg3 = x
farg4 = b
farg5 = y
farg6 = z
fresult = swigc_FN_VLinearSumVectorArray_Pthreads(farg1, farg2, farg3, farg4, farg5, farg6)
swig_result = fresult
end function

function FN_VScaleVectorArray_Pthreads(nvec, c, x, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
real(C_DOUBLE), dimension(*), target, intent(inout) :: c
type(C_PTR) :: x
type(C_PTR) :: z
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4

farg1 = nvec
farg2 = c_loc(c(1))
farg3 = x
farg4 = z
fresult = swigc_FN_VScaleVectorArray_Pthreads(farg1, farg2, farg3, farg4)
swig_result = fresult
end function

function FN_VConstVectorArray_Pthreads(nvec, c, z) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
real(C_DOUBLE), intent(in) :: c
type(C_PTR) :: z
integer(C_INT) :: fresult
integer(C_INT) :: farg1
real(C_DOUBLE) :: farg2
type(C_PTR) :: farg3

farg1 = nvec
farg2 = c
farg3 = z
fresult = swigc_FN_VConstVectorArray_Pthreads(farg1, farg2, farg3)
swig_result = fresult
end function

function FN_VWrmsNormVectorArray_Pthreads(nvec, x, w, nrm) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
type(C_PTR) :: x
type(C_PTR) :: w
real(C_DOUBLE), dimension(*), target, intent(inout) :: nrm
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4

farg1 = nvec
farg2 = x
farg3 = w
farg4 = c_loc(nrm(1))
fresult = swigc_FN_VWrmsNormVectorArray_Pthreads(farg1, farg2, farg3, farg4)
swig_result = fresult
end function

function FN_VWrmsNormMaskVectorArray_Pthreads(nvec, x, w, id, nrm) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
integer(C_INT), intent(in) :: nvec
type(C_PTR) :: x
type(C_PTR) :: w
type(N_Vector), target, intent(inout) :: id
real(C_DOUBLE), dimension(*), target, intent(inout) :: nrm
integer(C_INT) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3
type(C_PTR) :: farg4
type(C_PTR) :: farg5

farg1 = nvec
farg2 = x
farg3 = w
farg4 = c_loc(id)
farg5 = c_loc(nrm(1))
fresult = swigc_FN_VWrmsNormMaskVectorArray_Pthreads(farg1, farg2, farg3, farg4, farg5)
swig_result = fresult
end function

function FN_VWSqrSumLocal_Pthreads(x, w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: w
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(w)
fresult = swigc_FN_VWSqrSumLocal_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VWSqrSumMaskLocal_Pthreads(x, w, id) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
real(C_DOUBLE) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(N_Vector), target, intent(inout) :: w
type(N_Vector), target, intent(inout) :: id
real(C_DOUBLE) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2
type(C_PTR) :: farg3

farg1 = c_loc(x)
farg2 = c_loc(w)
farg3 = c_loc(id)
fresult = swigc_FN_VWSqrSumMaskLocal_Pthreads(farg1, farg2, farg3)
swig_result = fresult
end function

function FN_VBufSize_Pthreads(x, size) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: x
integer(C_INT64_T), dimension(*), target, intent(inout) :: size
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = c_loc(size(1))
fresult = swigc_FN_VBufSize_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VBufPack_Pthreads(x, buf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(C_PTR) :: buf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = buf
fresult = swigc_FN_VBufPack_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VBufUnpack_Pthreads(x, buf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: x
type(C_PTR) :: buf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
type(C_PTR) :: farg2

farg1 = c_loc(x)
farg2 = buf
fresult = swigc_FN_VBufUnpack_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableFusedOps_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableFusedOps_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableLinearCombination_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableLinearCombination_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableScaleAddMulti_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableScaleAddMulti_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableDotProdMulti_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableDotProdMulti_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableLinearSumVectorArray_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableLinearSumVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableScaleVectorArray_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableScaleVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableConstVectorArray_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableConstVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableWrmsNormVectorArray_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableWrmsNormVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VEnableWrmsNormMaskVectorArray_Pthreads(v, tf) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
integer(C_INT) :: swig_result
type(N_Vector), target, intent(inout) :: v
integer(C_INT), intent(in) :: tf
integer(C_INT) :: fresult
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = c_loc(v)
farg2 = tf
fresult = swigc_FN_VEnableWrmsNormMaskVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VCloneVectorArray_Pthreads(count, w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(C_PTR) :: swig_result
integer(C_INT), intent(in) :: count
type(N_Vector), target, intent(inout) :: w
type(C_PTR) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2

farg1 = count
farg2 = c_loc(w)
fresult = swigc_FN_VCloneVectorArray_Pthreads(farg1, farg2)
swig_result = fresult
end function

function FN_VCloneVectorArrayEmpty_Pthreads(count, w) &
result(swig_result)
use, intrinsic :: ISO_C_BINDING
type(C_PTR) :: swig_result
integer(C_INT), intent(in) :: count
type(N_Vector), target, intent(inout) :: w
type(C_PTR) :: fresult
integer(C_INT) :: farg1
type(C_PTR) :: farg2

farg1 = count
farg2 = c_loc(w)
fresult = swigc_FN_VCloneVectorArrayEmpty_Pthreads(farg1, farg2)
swig_result = fresult
end function

subroutine FN_VDestroyVectorArray_Pthreads(vs, count)
use, intrinsic :: ISO_C_BINDING
type(C_PTR) :: vs
integer(C_INT), intent(in) :: count
type(C_PTR) :: farg1
integer(C_INT) :: farg2

farg1 = vs
farg2 = count
call swigc_FN_VDestroyVectorArray_Pthreads(farg1, farg2)
end subroutine

end module


text:
cls

.( Demo 'create', 'does>', 'here', and , ) cr cr

\  0        1         2         3         4         5         6         7         8
\  12345678901234567890123456789012345678901234567890123456789012345678901234567890
.( The 'create' - 'does>' pair is the hardest part of forth language. Sir yap      )
.( taught me this way " forth 'create' - 'does>' is simply to define an OO class   )
.( with any number of properties but with only one method."                        ) cr

.( Consider forth was appeared around the same time as C language .. or earlier.   )
.( So early that Object Oriented concept was too aggressive to most people.        )
.( However, we can understand it now because we already have some knowledge of     )
.( Object Oriented programming.                                                    ) cr

.( The 'constant' command word is actually created by 'create' - 'does>' pair, its )
.( definition is,                                                                  ) cr
.(      : constant create ,  does>  @  ; 										   ) cr

.( The new word CONSTANT itself is a OO class. How to use it?                      )

.( 123 constant one-two-three )
123 constant one-two-three             \ now one-two-three is a new word created by the 'create' command in 'constant'
                                       \ it's value is given by the ',' command in 'constant' that compiles the TOS which was 123 into
                                       \ one-two-three's PFA. So this is the class constant's property for one-two-three.
                                       \ When one-two-three executed, actually executes the '@' command defined by does>.
                                       \ So does> defines the method of the 'constant' class. In this case the method is
                                       \ to fetch '@' the PFA of one-two-three. That's 123.

cr one-two-three . cr

bye

ext: f

max_stars_repo_name: hcchengithub/weforth

max_issues_repo_name: hcchengithub/weforth

max_forks_repo_name: hcchengithub/weforth


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE FLUXES (XS,YS,ICON,LISTEL,LISTND,UN,USAVE)
C
C     ******************************************************************
C
C     PROGRAM TO EVALUATE FLUX QUANTITIES (DEVIATORIC STRESSES,
C     VORTICITY, HEAT FLUX, ETC.) ON ELEMENT BOUNDARIES OR AT THE
C     ELEMENT INTEGRATION POINTS
C
C     ******************************************************************
C
      CHARACTER*10 CDATA,CODE
C
      COMMON /INDATR/ RDATA(150)
      COMMON /INDATI/ IDATA(150)
      COMMON /INDATC/ CDATA(150)
      COMMON /TAPES/  NIN,NOUT,NTP0,NTP1,NTP2,NTP3,NTP4,NTP5,NTP6,NTP7,
     1                NTP8,NTP9,NTP10,NTP11,NTP12,NTP13
      COMMON /SZDAT/  NUMEL,NUMNOD,NUMVAR,NUMDOF,NODSOL
      COMMON /MATDAT/ NMAT,NPROP,PROP(15,10),NXMAT,NXPROP,XPROP(15,10)
      COMMON /PRBDAT/ IAXSYM,ITMDEP,IFORCE,IFREE,IVAR1,IVAR2,IPFUNC,
     1                IPNLTY,PNLTY
      COMMON /ELMDAT/ NNELM(6),NNCOR(6),NNMID(6),NNCTR(6),NSIDET(5),
     1                NSIDEQ(7),NNSIDE(6,4,3)
      COMMON /TRI6/   F6Q(6,7),F6L(3,7),DF6QDS(6,7),DF6QDT(6,7)
      COMMON /QUAD8/  F8Q(8,9),F8L(4,9),DF8QDS(8,9),DF8QDT(8,9)
      COMMON /QUAD9/  F9Q(9,9),F9L(4,9),DF9QDS(9,9),DF9QDT(9,9)
      COMMON /FLXDAT/ NFLUX
C
      DIMENSION ICON(NUMEL,*), LISTEL(*), LISTND(*), XS(*), YS(*)
      DIMENSION UN(NUMNOD,*), USAVE(NUMNOD,*)
      DIMENSION X(9), Y(9), U(9), V(9), P(9), T(9), V1(9), V2(9)
      DIMENSION SQD(8), TQD(8), STRI(6), TTRI(6)
      DIMENSION SPTQA(8), TPTQA(8), SPTQB(4), TPTQB(4)
      DIMENSION SPTTA(6), TPTTA(6), SPTTB(3), TPTTB(3)
      DIMENSION TXX(8), TYY(8), TXY(8), TZZ(8), VORT(8)
      DIMENSION TN(8), TS(8), TSN(8), FN(4), FS(4)
      DIMENSION QXX(8), QYY(8), QN(8), QA(4)
      DIMENSION F1XX(8), F1YY(8), F2XX(8), F2YY(8)
      DIMENSION F1N(8), F2N(8), F1A(4), F2A(4)
      DIMENSION CODE(6), NEL(50)
C
      DATA (CODE(I),I=1,6)/'TRI6/3','TRI6/6','QUAD8/4','QUAD8/8',
     1     'QUAD9/4','QUAD9/9'/
      DATA (SPTQA(I),I=1,8)/-.5,.5,1.,1.,.5,-.5,-1.,-1./
      DATA (TPTQA(I),I=1,8)/-1.,-1.,-.5,.5,1.,1.,.5,-.5/
      DATA (SPTQB(I),I=1,4)/-.57735026,.57735026,.57735026,-.57735026/
      DATA (TPTQB(I),I=1,4)/-.57735026,-.57735026,.57735026,.57735026/
      DATA (SPTTA(I),I=1,6)/.75,.25,0.,0.,.25,.75/
      DATA (TPTTA(I),I=1,6)/.25,.75,.75,.25,0.,0./
      DATA (SPTTB(I),I=1,3)/.73333333,.13333333,.13333333/
      DATA (TPTTB(I),I=1,3)/.13333333,.73333333,.13333333/
C
      PARAMETER (PI=3.1415926535)
C
C     ******************************************************************
C
C     STATEMENT FUNCTIONS FOR DETERMINANTS AND ELEMENT QUANTITIES
C
      F11Q(S,R)=.25*(D1+S*D2+R*D3+S*R*D4+S**2*D5+S**2*R*D6)
      F12Q(S,R)=-.25*(B1+R*B2+S*B3+S*R*B4+R**2*B5+S*R**2*B6)
      F21Q(S,R)=-.25*(C1+S*C2+R*C3+S*R*C4+S**2*C5+S**2*R*C6)
      F22Q(S,R)=.25*(A1+R*A2+S*A3+S*R*A4+R**2*A5+S*R**2*A6)
      DETJQ(S,R)=F11Q(S,R)*F22Q(S,R)-F12Q(S,R)*F21Q(S,R)
C
      F11T(S,R)=S*A1+R*A2+A3
      F12T(S,R)=-(S*C1+R*C2+C3)
      F21T(S,R)=-(S*B1+R*B2+B3)
      F22T(S,R)=S*D1+R*D2+D3
      DETJT(S,R)=F11T(S,R)*F22T(S,R)-F12T(S,R)*F21T(S,R)
C
      U6Q(L)=F6Q(1,L)*U(1)+F6Q(2,L)*U(2)+F6Q(3,L)*U(3)+F6Q(4,L)*U(4)+
     1F6Q(5,L)*U(5)+F6Q(6,L)*U(6)
      U8Q(L)=F8Q(1,L)*U(1)+F8Q(2,L)*U(2)+F8Q(3,L)*U(3)+F8Q(4,L)*U(4)+
     1F8Q(5,L)*U(5)+F8Q(6,L)*U(6)+F8Q(7,L)*U(7)+F8Q(8,L)*U(8)
      U9Q(L)=F9Q(1,L)*U(1)+F9Q(2,L)*U(2)+F9Q(3,L)*U(3)+F9Q(4,L)*U(4)+
     1F9Q(5,L)*U(5)+F9Q(6,L)*U(6)+F9Q(7,L)*U(7)+F9Q(8,L)*U(8)+
     2F9Q(9,L)*U(9)
C
      T6Q(L)=F6Q(1,L)*T(1)+F6Q(2,L)*T(2)+F6Q(3,L)*T(3)+F6Q(4,L)*T(4)+
     1F6Q(5,L)*T(5)+F6Q(6,L)*T(6)
      T8Q(L)=F8Q(1,L)*T(1)+F8Q(2,L)*T(2)+F8Q(3,L)*T(3)+F8Q(4,L)*T(4)+
     1F8Q(5,L)*T(5)+F8Q(6,L)*T(6)+F8Q(7,L)*T(7)+F8Q(8,L)*T(8)
      T9Q(L)=F9Q(1,L)*T(1)+F9Q(2,L)*T(2)+F9Q(3,L)*T(3)+F9Q(4,L)*T(4)+
     1F9Q(5,L)*T(5)+F9Q(6,L)*T(6)+F9Q(7,L)*T(7)+F9Q(8,L)*T(8)+
     2F9Q(9,L)*T(9)
C
      V16Q(L)=F6Q(1,L)*V1(1)+F6Q(2,L)*V1(2)+F6Q(3,L)*V1(3)+
     1F6Q(4,L)*V1(4)+F6Q(5,L)*V1(5)+F6Q(6,L)*V1(6)
      V18Q(L)=F8Q(1,L)*V1(1)+F8Q(2,L)*V1(2)+F8Q(3,L)*V1(3)+
     1F8Q(4,L)*V1(4)+F8Q(5,L)*V1(5)+F8Q(6,L)*V1(6)+F8Q(7,L)*V1(7)+
     2F8Q(8,L)*V1(8)
      V19Q(L)=F9Q(1,L)*V1(1)+F9Q(2,L)*V1(2)+F9Q(3,L)*V1(3)+
     1F9Q(4,L)*V1(4)+F9Q(5,L)*V1(5)+F9Q(6,L)*V1(6)+F9Q(7,L)*V1(7)+
     2F9Q(8,L)*V1(8)+F9Q(9,L)*V1(9)
C
      V26Q(L)=F6Q(1,L)*V2(1)+F6Q(2,L)*V2(2)+F6Q(3,L)*V2(3)+
     1F6Q(4,L)*V2(4)+F6Q(5,L)*V2(5)+F6Q(6,L)*V2(6)
      V28Q(L)=F8Q(1,L)*V2(1)+F8Q(2,L)*V2(2)+F8Q(3,L)*V2(3)+
     1F8Q(4,L)*V2(4)+F8Q(5,L)*V2(5)+F8Q(6,L)*V2(6)+F8Q(7,L)*V2(7)+
     2F8Q(8,L)*V2(8)
      V29Q(L)=F9Q(1,L)*V2(1)+F9Q(2,L)*V2(2)+F9Q(3,L)*V2(3)+
     1F9Q(4,L)*V2(4)+F9Q(5,L)*V2(5)+F9Q(6,L)*V2(6)+F9Q(7,L)*V2(7)+
     2F9Q(8,L)*V2(8)+F9Q(9,L)*V2(9)
C
      RAD6L(L)=F6L(1,L)*X(1)+F6L(2,L)*X(2)+F6L(3,L)*X(3)
      RAD6Q(L)=F6Q(1,L)*X(1)+F6Q(2,L)*X(2)+F6Q(3,L)*X(3)+F6Q(4,L)*X(4)+
     1F6Q(5,L)*X(5)+F6Q(6,L)*X(6)
C
      RAD8L(L)=F8L(1,L)*X(1)+F8L(2,L)*X(2)+F8L(3,L)*X(3)+F8L(4,L)*X(4)
      RAD8Q(L)=F8Q(1,L)*X(1)+F8Q(2,L)*X(2)+F8Q(3,L)*X(3)+F8Q(4,L)*X(4)+
     1F8Q(5,L)*X(5)+F8Q(6,L)*X(6)+F8Q(7,L)*X(7)+F8Q(8,L)*X(8)
C
      RAD9L(L)=F9L(1,L)*X(1)+F9L(2,L)*X(2)+F9L(3,L)*X(3)+F9L(4,L)*X(4)
      RAD9Q(L)=F9Q(1,L)*X(1)+F9Q(2,L)*X(2)+F9Q(3,L)*X(3)+F9Q(4,L)*X(4)+
     1F9Q(5,L)*X(5)+F9Q(6,L)*X(6)+F9Q(7,L)*X(7)+F9Q(8,L)*X(8)+
     2F9Q(9,L)*X(9)
C
      P6L(L)=F6L(1,L)*P(1)+F6L(2,L)*P(2)+F6L(3,L)*P(3)
      P8L(L)=F8L(1,L)*P(1)+F8L(2,L)*P(2)+F8L(3,L)*P(3)+F8L(4,L)*P(4)
      P9L(L)=F9L(1,L)*P(1)+F9L(2,L)*P(2)+F9L(3,L)*P(3)+F9L(4,L)*P(4)
C
      DX(S)=.5*(2.*S-1.)*X1-2.*S*X2+.5*(2.*S+1.)*X3
      DY(S)=.5*(2.*S-1.)*Y1-2.*S*Y2+.5*(2.*S+1.)*Y3
      DS(S)=SQRT(DX(S)**2+DY(S)**2)
C
C     SET CONTROL PARAMETERS
C
      IWRTP=0
      ILOCM=0
      ILOCE=0
      IF (CDATA(3).EQ.'SAVE      ') THEN
      IWRTP=1
      ILOCE=1
      ILOCM=1
      WRITE (NOUT, 1430)
       ELSE
      IF (CDATA(2).EQ.'GAUSS     ') THEN
      ILOCE=1
      CDATA(4)='FULL'
      END IF
      IF (CDATA(4).EQ.'FULL'.OR.IDATA(4).EQ.0) THEN
      ILOCM=1
      WRITE (NOUT, 1440)
       ELSE
      DO 10 I=1,50
      IF (IDATA(I+3).EQ.0) GO TO 20
      NEL(I)=IDATA(I+3)
   10 CONTINUE
      I=51
   20 CONTINUE
      NUM=I-1
      WRITE (NOUT, 1160) (NEL(I),I=1,NUM)
      END IF
      END IF
      IF (ILOCE.EQ.0) WRITE (NOUT, 1140)
      IF (ILOCE.EQ.1) WRITE (NOUT, 1150)
C
C     EVALUATE SHAPE FUNCTIONS AND  DERIVATIVES AT THE FLUX
C     COMPUTATION POINTS
C
      IF (ILOCE.EQ.0) THEN
      NSTRSQ=8
      NSTRST=6
      DO 30 I=1,NSTRSQ
      SQD(I)=SPTQA(I)
      TQD(I)=TPTQA(I)
   30 CONTINUE
      DO 40 I=1,NSTRST
      STRI(I)=SPTTA(I)
      TTRI(I)=TPTTA(I)
   40 CONTINUE
       ELSE
      NSTRSQ=4
      NSTRST=3
      DO 50 I=1,NSTRSQ
      SQD(I)=SPTQB(I)
      TQD(I)=TPTQB(I)
   50 CONTINUE
      DO 60 I=1,NSTRST
      STRI(I)=SPTTB(I)
      TTRI(I)=TPTTB(I)
   60 CONTINUE
      END IF
      CALL FLUX6 (NSTRST,STRI,TTRI)
      CALL FLUX8 (NSTRSQ,SQD,TQD)
      CALL FLUX9 (NSTRSQ,SQD,TQD)
C
C     SET PARAMETERS FOR FLUX TYPE, TIME INCREMENTS, ETC.
C
      REWIND (NTP8)
   70 CONTINUE
      KOUNT=0
      NSTEPS=0
      NINC=0
      ITIME=0
      ISTRES=0
      IHEAT=0
      IFLUX1=0
      IFLUX2=0
      REWIND (NTP9)
      IF (IWRTP.EQ.1) THEN
      CDATA (1)='FLUXLOOP'
      CDATA (2)='ALLTIMES'
      GO TO 90
      END IF
   80 CONTINUE
      CALL RDFFLD (NIN,RDATA,IDATA,CDATA)
   90 CONTINUE
      IF (CDATA(1).EQ.'END       ') GO TO 1130
      IF (CDATA(1).EQ.'FLUXEND   ') GO TO 70
      IF (CDATA(1).EQ.'FLUXLOOP  ') GO TO 100
      IF (CDATA(1).EQ.'ALLFLUX   ') GO TO 110
      IF (CDATA(1).EQ.'STRESS    ') GO TO 120
      IF (CDATA(1).EQ.'HEATFLUX  ') GO TO 130
      IF (CDATA(1).EQ.'AUXFLUX1  ') GO TO 140
      IF (CDATA(1).EQ.'AUXFLUX2  ') GO TO 150
      CALL ERROR ('FLUXES','UNRECOGNIZED FLUX COMMAND' ,' ',0,' ',0,
     1'WORD',CDATA(1),1)
  100 CONTINUE
      IF (CDATA(2).EQ.'ALLTIMES  ') THEN
      NSTEPS=5000
      NINC=1
      ITIME=1
       ELSE
      NSTEPS=IDATA(2)
      NINC=IDATA(3)
      END IF
      GO TO 80
  110 CONTINUE
      ISTRES=1
      IF (IFORCE+IFREE.GT.0) IHEAT=1
      IF (IVAR1.EQ.1) IFLUX1=1
      IF (IVAR2.EQ.1) IFLUX2=1
      KSTEP=IDATA(2)
      IF (ITIME.EQ.1) KSTEP=1
      GO TO 160
  120 CONTINUE
      ISTRES=1
      KSTEP=IDATA(2)
      GO TO 160
  130 CONTINUE
      IF (IFORCE+IFREE.GT.0) IHEAT=1
      KSTEP=IDATA(2)
      GO TO 160
  140 CONTINUE
      IF (IVAR1.EQ.1) IFLUX1=1
      KSTEP=IDATA(2)
      GO TO 160
  150 CONTINUE
      IF (IVAR2.EQ.1) IFLUX2=1
      KSTEP=IDATA(2)
  160 CONTINUE
      IF (IWRTP.EQ.1) KSTEP=1
      IF (ITMDEP.EQ.0) THEN
      NSTEPS=1
      KSTEP=1
      NINC=0
      END IF
      WRITE (NOUT, 1460)
      IF (ISTRES.EQ.1) WRITE (NOUT, 1480)
      IF (IHEAT.EQ.1)  WRITE (NOUT, 1490)
      IF (IFLUX1.EQ.1) WRITE (NOUT, 1500)
      IF (IFLUX2.EQ.1) WRITE (NOUT, 1510)
      IF (IWRTP.EQ.1) GO TO 170
      IF (ITIME.EQ.1) THEN
      WRITE (NOUT, 1520)
       ELSE
      WRITE (NOUT, 1530) NSTEPS,KSTEP,NINC
      END IF
  170 CONTINUE
C
C     READ SOLUTION FIELD
C
      DO 180 K=1,KSTEP
      READ (NTP9, END=190) TIME,NUMNOD,NUMVAR,((UN(I,J),I=1,NUMNOD),
     1 J=1,NUMVAR)
  180 CONTINUE
      GO TO 200
  190 CONTINUE
      WRITE (NOUT, 1450)
      GO TO 80
  200 CONTINUE
      KOUNT=KOUNT+1
      IF (ILOCE.EQ.0) WRITE (NOUT, 1470) TIME
      DO 210 I=1,NUMNOD
      LISTND(I)=0
  210 CONTINUE
      DO 220 I=1,NUMNOD
      DO 220 J=1,NFLUX
      USAVE(I,J)=0.
  220 CONTINUE
C
C     LOOP ON ELEMENTS AND COMPUTE FLUXES
C
      DO 980 N=1,NUMEL
      IF (ILOCM.EQ.1) GO TO 240
      DO 230 I=1,NUM
      IF (N.EQ.NEL(I)) GO TO 240
  230 CONTINUE
      GO TO 980
  240 CONTINUE
      MAT=ABS(LISTEL(N))/100
      ISOLID=0
      IF (PROP(15,MAT).EQ.3.0) ISOLID=1
      NN=ABS(LISTEL(N))
      KIND=MOD(NN,100)
      NN=NNELM(KIND)
      DO 250 I=1,NN
      INODE=ICON(N,I)
      X(I)=XS(INODE)
      Y(I)=YS(INODE)
      U(I)=UN(INODE,1)
      V(I)=UN(INODE,2)
      P(I)=UN(INODE,3)
  250 CONTINUE
      IF (IFORCE+IFREE.GT.0) THEN
      DO 260 I=1,NN
      INODE=ICON(N,I)
      T(I)=UN(INODE,4)
  260 CONTINUE
      END IF
      IF (IVAR1.EQ.1) THEN
      DO 270 I=1,NN
      INODE=ICON(N,I)
      V1(I)=UN(INODE,5)
  270 CONTINUE
      END IF
      IF (IVAR2.EQ.1) THEN
      DO 280 I=1,NN
      INODE=ICON(N,I)
      V2(I)=UN(INODE,6)
  280 CONTINUE
      END IF
C
      GO TO (650, 660, 290, 300, 290, 310), KIND
C
C     SUBPARAMETRIC AND ISOPARAMETRIC QUADRILATERALS
C
  290 CONTINUE
      A1=-X(1)+X(2)+X(3)-X(4)
      A2=X(1)-X(2)+X(3)-X(4)
      A3=0.
      A4=0.
      A5=0.
      A6=0.
      B1=-Y(1)+Y(2)+Y(3)-Y(4)
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=0.
      B4=0.
      B5=0.
      B6=0.
      C1=-X(1)-X(2)+X(3)+X(4)
      C2=X(1)-X(2)+X(3)-X(4)
      C3=0.
      C4=0.
      C5=0.
      C6=0.
      D1=-Y(1)-Y(2)+Y(3)+Y(4)
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=0.
      D4=0.
      D5=0.
      D6=0.
      GO TO 320
  300 CONTINUE
      A1=2.*(X(6)-X(8))
      A2=X(1)-X(2)+X(3)-X(4)
      A3=2.*(X(1)+X(2)+X(3)+X(4)-2.*X(5)-2.*X(7))
      A4=2.*(-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7))
      A5=-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8)
      A6=0.
      B1=2.*(Y(6)-Y(8))
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*Y(5)-2.*Y(7))
      B4=2.*(-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7))
      B5=-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8)
      B6=0.
      C1=2.*(X(7)-X(5))
      C2=X(1)-X(2)+X(3)-X(4)
      C3=2.*(X(1)+X(2)+X(3)+X(4)-2.*X(6)-2.*X(8))
      C4=2.*(-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8))
      C5=-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7)
      C6=0.
      D1=2.*(Y(7)-Y(5))
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=2.*(Y(1)+Y(2)+Y(3)+Y(4)-2.*Y(6)-2.*Y(8))
      D4=2.*(-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8))
      D5=-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7)
      D6=0.
      GO TO 320
  310 CONTINUE
      A1=2.*(X(6)-X(8))
      A2=X(1)-X(2)+X(3)-X(4)
      A3=4.*(X(6)+X(8)-2.*X(9))
      A4=2.*(-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7))
      A5=-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8)
      A6=2.*(X(1)+X(2)+X(3)+X(4)-2.*(X(5)+X(6)+X(7)+X(8))+4.*X(9))
      B1=2.*(Y(6)-Y(8))
      B2=Y(1)-Y(2)+Y(3)-Y(4)
      B3=4.*(Y(6)+Y(8)-2.*Y(9))
      B4=2.*(-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7))
      B5=-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8)
      B6=2*(Y(1)+Y(2)+Y(3)+Y(4)-2.*(Y(5)+Y(6)+Y(7)+Y(8))+4.*Y(9))
      C1=2.*(X(7)-X(5))
      C2=X(1)-X(2)+X(3)-X(4)
      C3=4.*(X(5)+X(7)-2.*X(9))
      C4=2.*(-X(1)+X(2)+X(3)-X(4)-2.*X(6)+2.*X(8))
      C5=-X(1)-X(2)+X(3)+X(4)+2.*X(5)-2.*X(7)
      C6=2.*(X(1)+X(2)+X(3)+X(4)-2.*(X(5)+X(6)+X(7)+X(8))+4.*X(9))
      D1=2.*(Y(7)-Y(5))
      D2=Y(1)-Y(2)+Y(3)-Y(4)
      D3=4.*(Y(5)+Y(7)-2.*Y(9))
      D4=2.*(-Y(1)+Y(2)+Y(3)-Y(4)-2.*Y(6)+2.*Y(8))
      D5=-Y(1)-Y(2)+Y(3)+Y(4)+2.*Y(5)-2.*Y(7)
      D6=2*(Y(1)+Y(2)+Y(3)+Y(4)-2.*(Y(5)+Y(6)+Y(7)+Y(8))+4.*Y(9))
  320 CONTINUE
C
C     STRESS AND VORTICITY EVALUATION
C
      IF (ISTRES.EQ.0) GO TO 430
      DO 380 I=1,NSTRSQ
      QNT1=0.
      QNT2=0.
      QNT3=0.
      QNT4=0.
      A=F11Q(SQD(I),TQD(I))
      B=F21Q(SQD(I),TQD(I))
      C=F12Q(SQD(I),TQD(I))
      D=F22Q(SQD(I),TQD(I))
      E=1./DETJQ(SQD(I),TQD(I))
      IF (KIND.GE.5) GO TO 340
C
C     EIGHT-NODE ELEMENTS
C
      RG=RAD8L(I)
      IF (KIND.EQ.4) RG=RAD8Q(I)
      UG=U8Q(I)
      DO 330 J=1,NN
      AA=DF8QDS(J,I)*A+DF8QDT(J,I)*C
      BB=DF8QDS(J,I)*B+DF8QDT(J,I)*D
      QNT1=QNT1+AA*U(J)*E
      QNT2=QNT2+BB*V(J)*E
      QNT3=QNT3+AA*V(J)*E+BB*U(J)*E
      QNT4=QNT4+(BB*U(J)-AA*V(J))*E
  330 CONTINUE
      PRESS=P8L(I)
      PARAM1=PROP(2,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(3,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 370
      TLOC=T8Q(I)
      SLOC=2.*(QNT1**2+QNT2**2)+(QNT3**2)+IAXSYM*2.*(UG/RG)**2
      V1LOC=V18Q(I)
      V2LOC=V28Q(I)
      GO TO 360
C
C     NINE-NODE ELEMENTS
C
  340 CONTINUE
      RG=RAD9L(I)
      IF (KIND.EQ.6) RG=RAD9Q(I)
      UG=U9Q(I)
      DO 350 J=1,NN
      AA=DF9QDS(J,I)*A+DF9QDT(J,I)*C
      BB=DF9QDS(J,I)*B+DF9QDT(J,I)*D
      QNT1=QNT1+AA*U(J)*E
      QNT2=QNT2+BB*V(J)*E
      QNT3=QNT3+AA*V(J)*E+BB*U(J)*E
      QNT4=QNT4+(BB*U(J)-AA*V(J))*E
  350 CONTINUE
      PRESS=P9L(I)
      PARAM1=PROP(2,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(3,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 370
      TLOC=T9Q(I)
      SLOC=2.*(QNT1**2+QNT2**2)+(QNT3**2)+IAXSYM*2.*(UG/RG)**2
      V1LOC=V19Q(I)
      V2LOC=V29Q(I)
  360 CONTINUE
C
      CALL USRVIS (PARAM1,TLOC,SLOC,V1LOC,V2LOC,1,MAT,0)
  370 CONTINUE
      TXX(I)=2.*PARAM1*QNT1
      TYY(I)=2.*PARAM1*QNT2
      TXY(I)=PARAM1*QNT3
      VORT(I)=QNT4
      IF (IAXSYM.EQ.1) TZZ(I)=2.*PARAM1*(UG/RG)
  380 CONTINUE
      NSTRS=NSTRSQ
C
      IF (ILOCE.EQ.1) GO TO 420
      NSIDE=NSTRSQ/2
      DO 410 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X1=X(1)
      IF (I.EQ.NSIDE) Y1=Y(1)
      GO TO (390, 390, 390, 400, 390, 400), KIND
  390 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(II)=TXX(II)*(COS(ANG)**2)+TYY(II)*(SIN(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TS(II)=TXX(II)*(SIN(ANG)**2)+TYY(II)*(COS(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TSN(II)=(-TXX(II)+TYY(II))*SIN(ANG)*COS(ANG)+
     1TXY(II)*(COS(ANG)**2-SIN(ANG)**2)
      TN(JJ)=TXX(JJ)*(COS(ANG)**2)+TYY(JJ)*(SIN(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TS(JJ)=TXX(JJ)*(SIN(ANG)**2)+TYY(JJ)*(COS(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TSN(JJ)=(-TXX(JJ)+TYY(JJ))*SIN(ANG)*COS(ANG)+
     1TXY(JJ)*(COS(ANG)**2-SIN(ANG)**2)
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      FN(I)=(TN(II)+TN(JJ))*DSP*FACTOR
      FS(I)=(TSN(II)+TSN(JJ))*DSP*FACTOR
      GO TO 410
  400 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(II)=TXX(II)*(COS(ANG)**2)+TYY(II)*(SIN(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TS(II)=TXX(II)*(SIN(ANG)**2)+TYY(II)*(COS(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TSN(II)=(-TXX(II)+TYY(II))*SIN(ANG)*COS(ANG)+
     1TXY(II)*(COS(ANG)**2-SIN(ANG)**2)
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(JJ)=TXX(JJ)*(COS(ANG)**2)+TYY(JJ)*(SIN(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TS(JJ)=TXX(JJ)*(SIN(ANG)**2)+TYY(JJ)*(COS(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TSN(JJ)=(-TXX(JJ)+TYY(JJ))*SIN(ANG)*COS(ANG)+
     1TXY(JJ)*(COS(ANG)**2-SIN(ANG)**2)
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      FN(I)=TN(II)*XSIDE+TN(JJ)*YSIDE
      FS(I)=TSN(II)*XSIDE+TSN(JJ)*YSIDE
  410 CONTINUE
      GO TO 430
  420 CONTINUE
      CALL EXTRAP (KIND,TXX)
      CALL EXTRAP (KIND,TYY)
      CALL EXTRAP (KIND,TXY)
      IF (IAXSYM.EQ.1) CALL EXTRAP (KIND,TZZ)
      CALL EXTRAP (KIND,VORT)
C
C     HEAT FLUX EVALUATION
C
  430 CONTINUE
      IF (IHEAT.EQ.0) GO TO 540
      DO 490 I=1,NSTRSQ
      QNT1=0.
      QNT2=0.
      A=F11Q(SQD(I),TQD(I))
      B=F21Q(SQD(I),TQD(I))
      C=F12Q(SQD(I),TQD(I))
      D=F22Q(SQD(I),TQD(I))
      E=1./DETJQ(SQD(I),TQD(I))
      IF (KIND.GE.5) GO TO 450
C
C     EIGHT-NODE ELEMENTS
C
      DO 440 J=1,NN
      AA=DF8QDS(J,I)*A+DF8QDT(J,I)*C
      BB=DF8QDS(J,I)*B+DF8QDT(J,I)*D
      QNT1=QNT1+AA*T(J)*E
      QNT2=QNT2+BB*T(J)*E
  440 CONTINUE
      PARAM1=PROP(4,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(7,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 480
      TLOC=T8Q(I)
      V1LOC=V18Q(I)
      V2LOC=V28Q(I)
      GO TO 470
C
C     NINE-NODE ELEMENTS
C
  450 CONTINUE
      DO 460 J=1,NN
      AA=DF9QDS(J,I)*A+DF9QDT(J,I)*C
      BB=DF9QDS(J,I)*B+DF9QDT(J,I)*D
      QNT1=QNT1+AA*T(J)*E
      QNT2=QNT2+BB*T(J)*E
  460 CONTINUE
      PARAM1=PROP(4,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(7,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 480
      TLOC=T9Q(I)
      SLOC=1.0
      V1LOC=V19Q(I)
      V2LOC=V29Q(I)
  470 CONTINUE
C
      CALL USRCON (PARAM1,TLOC,SLOC,V1LOC,V2LOC,1,MAT)
  480 CONTINUE
      QXX(I)=-PARAM1*QNT1
      QYY(I)=-PARAM1*QNT2
  490 CONTINUE
      NSTRS=NSTRSQ
C
      IF (ILOCE.EQ.1) GO TO 530
      NSIDE=NSTRSQ/2
      DO 520 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X3=X(1)
      IF (I.EQ.NSIDE) Y3=Y(1)
      GO TO (500, 500, 500, 510, 500, 510), KIND
  500 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      QN(II)=(QXX(II)*YSIDE-QYY(II)*XSIDE)/DSP
      QN(JJ)=(QXX(JJ)*YSIDE-QYY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      QA(I)=(QN(II)+QN(JJ))*DSP*FACTOR
      GO TO 520
  510 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      QN(II)=(QXX(II)*YSIDE-QYY(II)*XSIDE)/DSP
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      QN(JJ)=(QXX(JJ)*YSIDE-QYY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      QA(I)=QN(II)*XSIDE+QN(JJ)*YSIDE
  520 CONTINUE
      GO TO 540
  530 CONTINUE
      CALL EXTRAP (KIND,QXX)
      CALL EXTRAP (KIND,QYY)
C
C     VARIABLE 1 AND VARIABLE 2 FLUX EVALUATION
C
  540 CONTINUE
      IF (IFLUX1.EQ.0.AND.IFLUX2.EQ.0) GO TO 910
      DO 600 I=1,NSTRSQ
      QNT1=0.
      QNT2=0.
      QNT3=0.
      QNT4=0.
      A=F11Q(SQD(I),TQD(I))
      B=F21Q(SQD(I),TQD(I))
      C=F12Q(SQD(I),TQD(I))
      D=F22Q(SQD(I),TQD(I))
      E=1./DETJQ(SQD(I),TQD(I))
      IF (KIND.GE.5) GO TO 560
C
C     EIGHT-NODE ELEMENTS
C
      DO 550 J=1,NN
      AA=DF8QDS(J,I)*A+DF8QDT(J,I)*C
      BB=DF8QDS(J,I)*B+DF8QDT(J,I)*D
      QNT1=QNT1+AA*V1(J)*E
      QNT2=QNT2+BB*V1(J)*E
      QNT3=QNT3+AA*V2(J)*E
      QNT4=QNT4+BB*V2(J)*E
  550 CONTINUE
      PARAM1=XPROP(2,MAT)
      PARAM2=XPROP(8,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 590
      TLOC=T8Q(I)
      V1LOC=V18Q(I)
      V2LOC=V28Q(I)
      GO TO 580
C
C     NINE-NODE ELEMENTS
C
  560 CONTINUE
      DO 570 J=1,NN
      AA=DF9QDS(J,I)*A+DF9QDT(J,I)*C
      BB=DF9QDS(J,I)*B+DF9QDT(J,I)*D
      QNT1=QNT1+AA*V1(J)*E
      QNT2=QNT2+BB*V1(J)*E
      QNT3=QNT3+AA*V2(J)*E
      QNT4=QNT4+BB*V2(J)*E
  570 CONTINUE
      PARAM1=XPROP(2,MAT)
      PARAM2=XPROP(8,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 590
      TLOC=T9Q(I)
      V1LOC=V19Q(I)
      V2LOC=V29Q(I)
  580 CONTINUE
C
      IF (IVAR1.EQ.1) CALL USRDIF (PARAM1,TLOC,V1LOC,V2LOC,1,MAT,1)
      IF (IVAR2.EQ.1) CALL USRDIF (PARAM2,TLOC,V1LOC,V2LOC,1,MAT,2)
  590 CONTINUE
      F1XX(I)=-PARAM1*QNT1
      F1YY(I)=-PARAM1*QNT2
      F2XX(I)=-PARAM2*QNT3
      F2YY(I)=-PARAM2*QNT4
  600 CONTINUE
      NSTRS=NSTRSQ
C
      IF (ILOCE.EQ.1) GO TO 640
      NSIDE=NSTRSQ/2
      DO 630 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X3=X(1)
      IF (I.EQ.NSIDE) Y3=Y(1)
      GO TO (610, 620, 610, 620, 610, 620), KIND
  610 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      F1N(II)=(F1XX(II)*YSIDE-F1YY(II)*XSIDE)/DSP
      F1N(JJ)=(F1XX(JJ)*YSIDE-F1YY(JJ)*XSIDE)/DSP
      F2N(II)=(F2XX(II)*YSIDE-F2YY(II)*XSIDE)/DSP
      F2N(JJ)=(F2XX(JJ)*YSIDE-F2YY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      F1A(I)=(F1N(II)+F1N(JJ))*DSP*FACTOR
      F2A(I)=(F2N(II)+F2N(JJ))*DSP*FACTOR
      GO TO 630
  620 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      F1N(II)=(F1XX(II)*YSIDE-F1YY(II)*XSIDE)/DSP
      F2N(II)=(F2XX(II)*YSIDE-F2YY(II)*XSIDE)/DSP
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      F1N(JJ)=(F1XX(JJ)*YSIDE-F1YY(JJ)*XSIDE)/DSP
      F2N(JJ)=(F2XX(JJ)*YSIDE-F2YY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      F1A(I)=F1N(II)*XSIDE+F1N(JJ)*YSIDE
      F2A(I)=F2N(II)*XSIDE+F2N(JJ)*YSIDE
  630 CONTINUE
      GO TO 910
  640 CONTINUE
      CALL EXTRAP (KIND,F1XX)
      CALL EXTRAP (KIND,F1YY)
      IF (IFLUX2.EQ.1) THEN
      CALL EXTRAP (KIND,F2XX)
      CALL EXTRAP (KIND,F2YY)
      END IF
      GO TO 910
C
C     SUBPARAMETRIC AND ISOPARAMETRIC TRIANGLE
C
  650 CONTINUE
      A1=0.
      A2=0.
      A3=Y(2)-Y(3)
      B1=0.
      B2=0.
      B3=-X(2)+X(3)
      C1=0.
      C2=0.
      C3=-Y(1)+Y(3)
      D1=0.
      D2=0.
      D3=X(1)-X(3)
      GO TO 670
  660 CONTINUE
      A1=4.*(Y(3)+Y(4)-Y(5)-Y(6))
      A2=4.*(Y(2)+Y(3)-2.*Y(5))
      A3=4.*Y(5)-3.*Y(3)-Y(2)
      B1=4.*(X(3)+X(4)-X(5)-X(6))
      B2=4.*(X(2)+X(3)-2.*X(5))
      B3=4.*X(5)-3.*X(3)-X(2)
      C1=4.*(Y(1)+Y(3)-2.*Y(6))
      C2=4.*(Y(3)+Y(4)-Y(5)-Y(6))
      C3=4.*Y(6)-3.*Y(3)-Y(1)
      D1=4.*(X(1)+X(3)-2.*X(6))
      D2=4.*(X(3)+X(4)-X(5)-X(6))
      D3=4.*X(6)-3.*X(3)-X(1)
  670 CONTINUE
      IF (ISTRES.EQ.0) GO TO 750
      DO 700 I=1,NSTRST
      QNT1=0.
      QNT2=0.
      QNT3=0.
      QNT4=0.
      A=F11T(STRI(I),TTRI(I))
      B=F21T(STRI(I),TTRI(I))
      C=F12T(STRI(I),TTRI(I))
      D=F22T(STRI(I),TTRI(I))
      E=1./DETJT(STRI(I),TTRI(I))
      RG=RAD6L(I)
      IF (KIND.EQ.2) RG=RAD6Q(I)
      UG=U6Q(I)
      DO 680 J=1,NN
      AA=DF6QDS(J,I)*A+DF6QDT(J,I)*C
      BB=DF6QDS(J,I)*B+DF6QDT(J,I)*D
      QNT1=QNT1+AA*U(J)*E
      QNT2=QNT2+BB*V(J)*E
      QNT3=QNT3+AA*V(J)*E+BB*U(J)*E
      QNT4=QNT4+(BB*U(J)-AA*V(J))*E
  680 CONTINUE
      PRESS=P6L(I)
      PARAM1=PROP(2,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(3,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 690
      TLOC=T6Q(I)
      SLOC=2.*(QNT1**2+QNT2**2)+(QNT3**2)+IAXSYM*2.*(UG/RG)**2
      V1LOC=V16Q(I)
      V2LOC=V26Q(I)
C
      CALL USRVIS (PARAM1,TLOC,SLOC,V1LOC,V2LOC,1,MAT,0)
  690 CONTINUE
      TXX(I)=2.*PARAM1*QNT1
      TYY(I)=2.*PARAM1*QNT2
      TXY(I)=PARAM1*QNT3
      VORT(I)=QNT4
      IF (IAXSYM.EQ.1) TZZ(I)=2.*PARAM1*(UG/RG)
  700 CONTINUE
      NSTRS=NSTRST
C
      IF (ILOCE.EQ.1) GO TO 740
      NSIDE=NSTRST/2
      DO 730 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X1=X(1)
      IF (I.EQ.NSIDE) Y1=Y(1)
      GO TO (710, 720, 720, 720, 720, 720), KIND
  710 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(II)=TXX(II)*(COS(ANG)**2)+TYY(II)*(SIN(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TS(II)=TXX(II)*(SIN(ANG)**2)+TYY(II)*(COS(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TSN(II)=(-TXX(II)+TYY(II))*SIN(ANG)*COS(ANG)+
     1TXY(II)*(COS(ANG)**2-SIN(ANG)**2)
      TN(JJ)=TXX(JJ)*(COS(ANG)**2)+TYY(JJ)*(SIN(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TS(JJ)=TXX(JJ)*(SIN(ANG)**2)+TYY(JJ)*(COS(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TSN(JJ)=(-TXX(JJ)+TYY(JJ))*SIN(ANG)*COS(ANG)+
     1TXY(JJ)*(COS(ANG)**2-SIN(ANG)**2)
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      FN(I)=(TN(II)+TN(JJ))*DSP*FACTOR
      FS(I)=(TSN(II)+TSN(JJ))*DSP*FACTOR
      GO TO 730
  720 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(II)=TXX(II)*(COS(ANG)**2)+TYY(II)*(SIN(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TS(II)=TXX(II)*(SIN(ANG)**2)+TYY(II)*(COS(ANG)**2)+
     12.*TXY(II)*SIN(ANG)*COS(ANG)
      TSN(II)=(-TXX(II)+TYY(II))*SIN(ANG)*COS(ANG)+
     1TXY(II)*(COS(ANG)**2-SIN(ANG)**2)
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      A=XSIDE/DSP
      B=YSIDE/DSP
      IF (B.EQ.0.) B=1.0E-15
      ANG=ATAN2(A,B)
      TN(JJ)=TXX(JJ)*(COS(ANG)**2)+TYY(JJ)*(SIN(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TS(JJ)=TXX(JJ)*(SIN(ANG)**2)+TYY(JJ)*(COS(ANG)**2)+
     12.*TXY(JJ)*SIN(ANG)*COS(ANG)
      TSN(JJ)=(-TXX(JJ)+TYY(JJ))*SIN(ANG)*COS(ANG)+
     1TXY(JJ)*(COS(ANG)**2-SIN(ANG)**2)
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      FN(I)=TN(II)*XSIDE+TN(JJ)*YSIDE
      FS(I)=TSN(II)*XSIDE+TSN(JJ)*YSIDE
  730 CONTINUE
      GO TO 750
  740 CONTINUE
      CALL EXTRAP (KIND,TXX)
      CALL EXTRAP (KIND,TYY)
      CALL EXTRAP (KIND,TXY)
      IF (IAXSYM.EQ.1) CALL EXTRAP (KIND,TZZ)
      CALL EXTRAP (KIND,VORT)
C
C     HEAT FLUX EVALUATION
C
  750 CONTINUE
      IF (IHEAT.EQ.0) GO TO 830
      DO 780 I=1,NSTRST
      QNT1=0.
      QNT2=0.
      A=F11T(STRI(I),TTRI(I))
      B=F21T(STRI(I),TTRI(I))
      C=F12T(STRI(I),TTRI(I))
      D=F22T(STRI(I),TTRI(I))
      E=1./DETJT(STRI(I),TTRI(I))
      DO 760 J=1,NN
      AA=DF6QDS(J,I)*A+DF6QDT(J,I)*C
      BB=DF6QDS(J,I)*B+DF6QDT(J,I)*D
      QNT1=QNT1+AA*T(J)*E
      QNT2=QNT2+BB*T(J)*E
  760 CONTINUE
      PARAM1=PROP(4,MAT)
      IF (PROP(15,MAT).EQ.4.0) PARAM1=PROP(7,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 770
      TLOC=T6Q(I)
      SLOC=1.0
      V1LOC=V16Q(I)
      V2LOC=V26Q(I)
C
      CALL USRCON (PARAM1,TLOC,SLOC,V1LOC,V2LOC,1,MAT)
  770 CONTINUE
      QXX(I)=-PARAM1*QNT1
      QYY(I)=-PARAM1*QNT2
  780 CONTINUE
      NSTRS=NSTRST
C
      IF (ILOCE.EQ.1) GO TO 820
      NSIDE=NSTRST/2
      DO 810 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X3=X(1)
      IF (I.EQ.NSIDE) Y3=Y(1)
      GO TO (790, 800, 800, 800, 800, 800), KIND
  790 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      QN(II)=(QXX(II)*YSIDE-QYY(II)*XSIDE)/DSP
      QN(JJ)=(QXX(JJ)*YSIDE-QYY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      QA(I)=(QN(II)+QN(JJ))*DSP*FACTOR
      GO TO 810
  800 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      QN(II)=(QXX(II)*YSIDE-QYY(II)*XSIDE)/DSP
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      QN(JJ)=(QXX(JJ)*YSIDE-QYY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      QA(I)=QN(II)*XSIDE+QN(JJ)*YSIDE
  810 CONTINUE
      GO TO 830
  820 CONTINUE
      CALL EXTRAP (KIND,QXX)
      CALL EXTRAP (KIND,QYY)
C
C     VARIABLE 1 AND VARIABLE 2 FLUX EVALUATION
C
  830 CONTINUE
      IF (IFLUX1.EQ.0.AND.IFLUX2.EQ.0) GO TO 910
      DO 860 I=1,NSTRST
      QNT1=0.
      QNT2=0.
      QNT3=0.
      QNT4=0.
      A=F11T(STRI(I),TTRI(I))
      B=F21T(STRI(I),TTRI(I))
      C=F12T(STRI(I),TTRI(I))
      D=F22T(STRI(I),TTRI(I))
      E=1./DETJT(STRI(I),TTRI(I))
      DO 840 J=1,NN
      AA=DF6QDS(J,I)*A+DF6QDT(J,I)*C
      BB=DF6QDS(J,I)*B+DF6QDT(J,I)*D
      QNT1=QNT1+AA*V1(J)*E
      QNT2=QNT2+BB*V1(J)*E
      QNT3=QNT3+AA*V2(J)*E
      QNT4=QNT4+BB*V2(J)*E
  840 CONTINUE
      PARAM1=XPROP(2,MAT)
      PARAM2=XPROP(8,MAT)
      IF (PROP(8,MAT).EQ.1.0) GO TO 850
      TLOC=T6Q(I)
      V1LOC=V16Q(I)
      V2LOC=V26Q(I)
C
      IF (IVAR1.EQ.1) CALL USRDIF (PARAM1,TLOC,V1LOC,V2LOC,1,MAT,1)
      IF (IVAR2.EQ.1) CALL USRDIF (PARAM2,TLOC,V1LOC,V2LOC,1,MAT,2)
  850 CONTINUE
      F1XX(I)=-PARAM1*QNT1
      F1YY(I)=-PARAM1*QNT2
      F2XX(I)=-PARAM2*QNT3
      F2YY(I)=-PARAM2*QNT4
  860 CONTINUE
      NSTRS=NSTRST
C
      IF (ILOCE.EQ.1) GO TO 900
      NSIDE=NSTRST/2
      DO 890 I=1,NSIDE
      X1=X(I)
      X2=X(I+NSIDE)
      X3=X(I+1)
      Y1=Y(I)
      Y2=Y(I+NSIDE)
      Y3=Y(I+1)
      IF (I.EQ.NSIDE) X3=X(1)
      IF (I.EQ.NSIDE) Y3=Y(1)
      GO TO (870, 880, 880, 880, 880, 880), KIND
  870 CONTINUE
      XSIDE=(X3-X1)/2.
      YSIDE=(Y3-Y1)/2.
      DSP=SQRT(XSIDE**2+YSIDE**2)
      II=2*I-1
      JJ=2*I
      F1N(II)=(F1XX(II)*YSIDE-F1YY(II)*XSIDE)/DSP
      F1N(JJ)=(F1XX(JJ)*YSIDE-F1YY(JJ)*XSIDE)/DSP
      F2N(II)=(F2XX(II)*YSIDE-F2YY(II)*XSIDE)/DSP
      F2N(JJ)=(F2XX(JJ)*YSIDE-F2YY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(X3+X1)
      F1A(I)=(F1N(II)+F1N(JJ))*DSP*FACTOR
      F2A(I)=(F2N(II)+F2N(JJ))*DSP*FACTOR
      GO TO 890
  880 CONTINUE
      II=2*I-1
      JJ=2*I
      XSIDE=DX(-.5)
      YSIDE=DY(-.5)
      DSP=DS(-.5)
      F1N(II)=(F1XX(II)*YSIDE-F1YY(II)*XSIDE)/DSP
      F2N(II)=(F2XX(II)*YSIDE-F2YY(II)*XSIDE)/DSP
      XSIDE=DX(.5)
      YSIDE=DY(.5)
      DSP=DS(.5)
      F1N(JJ)=(F1XX(JJ)*YSIDE-F1YY(JJ)*XSIDE)/DSP
      F2N(JJ)=(F2XX(JJ)*YSIDE-F2YY(JJ)*XSIDE)/DSP
      FACTOR=1.0
      IF (IAXSYM.EQ.1) FACTOR=PI*(.75*X1+1.5*X2-.25*X3)
      XSIDE=(X2-X1)**2+(Y2-Y1)**2
      XSIDE=SQRT(XSIDE)*FACTOR
      IF (IAXSYM.EQ.1) FACTOR=PI*(-.25*X1+1.5*X2+.75*X3)
      YSIDE=(X3-X2)**2+(Y3-Y2)**2
      YSIDE=SQRT(YSIDE)*FACTOR
      F1A(I)=F1N(II)*XSIDE+F1N(JJ)*YSIDE
      F2A(I)=F2N(II)*XSIDE+F2N(JJ)*YSIDE
  890 CONTINUE
      GO TO 910
  900 CONTINUE
      CALL EXTRAP (KIND,F1XX)
      CALL EXTRAP (KIND,F1YY)
      IF (IFLUX2.EQ.1) THEN
      CALL EXTRAP (KIND,F2XX)
      CALL EXTRAP (KIND,F2YY)
      END IF
C
C     OUTPUT BOUNDARY POINT FLUXES
C
  910 CONTINUE
      IF (ILOCE.EQ.1) GO TO 950
      WRITE (NOUT, 1170) N,CODE(KIND)
      IF (ISTRES.EQ.0) GO TO 920
      IF (ISOLID.EQ.1) GO TO 920
      WRITE (NOUT, 1180)
      WRITE (NOUT, 1190) (TXX(I),I=1,NSTRS)
      WRITE (NOUT, 1200) (TYY(I),I=1,NSTRS)
      WRITE (NOUT, 1210) (TXY(I),I=1,NSTRS)
      IF (IAXSYM.EQ.1) WRITE (NOUT, 1220) (TZZ(I),I=1,NSTRS)
      WRITE (NOUT, 1230) (VORT(I),I=1,NSTRS)
      WRITE (NOUT, 1240) (TN(I),I=1,NSTRS)
      WRITE (NOUT, 1250) (TS(I),I=1,NSTRS)
      WRITE (NOUT, 1260) (TSN(I),I=1,NSTRS)
      WRITE (NOUT, 1270) (FN(I),I=1,NSIDE)
      WRITE (NOUT, 1280) (FS(I),I=1,NSIDE)
  920 CONTINUE
      IF (IHEAT.EQ.0) GO TO 930
      WRITE (NOUT, 1290)
      WRITE (NOUT, 1300) (QXX(I),I=1,NSTRS)
      WRITE (NOUT, 1310) (QYY(I),I=1,NSTRS)
      WRITE (NOUT, 1320) (QN(I),I=1,NSTRS)
      WRITE (NOUT, 1330) (QA(I),I=1,NSIDE)
  930 CONTINUE
      IF (IFLUX1+IFLUX2.GT.0) WRITE (NOUT, 1340)
      IF (IFLUX1.EQ.0) GO TO 940
      WRITE (NOUT, 1350) (F1XX(I),I=1,NSTRS)
      WRITE (NOUT, 1360) (F1YY(I),I=1,NSTRS)
      WRITE (NOUT, 1370) (F1N(I),I=1,NSTRS)
      WRITE (NOUT, 1380) (F1A(I),I=1,NSIDE)
  940 CONTINUE
      IF (IFLUX2.EQ.0) GO TO 980
      WRITE (NOUT, 1390) (F2XX(I),I=1,NSTRS)
      WRITE (NOUT, 1400) (F2YY(I),I=1,NSTRS)
      WRITE (NOUT, 1410) (F2N(I),I=1,NSTRS)
      WRITE (NOUT, 1420) (F2A(I),I=1,NSIDE)
      GO TO 980
C
C     ACCUMULATE EXTRAPOLATED INTEGRATION POINT FLUXES FOR
C     LATER AVERAGING AT NODES
C
  950 CONTINUE
      DO 970 I=1,NSTRS
      INODE=ICON(N,I)
      IF (ISOLID.EQ.1) GO TO 960
      IF (ISTRES.EQ.1) THEN
      USAVE(INODE,1)=USAVE(INODE,1)+TXX(I)
      USAVE(INODE,2)=USAVE(INODE,2)+TYY(I)
      USAVE(INODE,3)=USAVE(INODE,3)+TXY(I)
      USAVE(INODE,4)=USAVE(INODE,4)+TZZ(I)
      USAVE(INODE,5)=USAVE(INODE,5)+VORT(I)
      END IF
  960 CONTINUE
      IF (IHEAT.EQ.1) THEN
      USAVE(INODE,6)=USAVE(INODE,6)+QXX(I)
      USAVE(INODE,7)=USAVE(INODE,7)+QYY(I)
      END IF
      IF (IFLUX1.EQ.1) THEN
      USAVE(INODE,8)=USAVE(INODE,8)+F1XX(I)
      USAVE(INODE,9)=USAVE(INODE,9)+F1YY(I)
      END IF
      IF (IFLUX2.EQ.1) THEN
      USAVE(INODE,10)=USAVE(INODE,10)+F2XX(I)
      USAVE(INODE,11)=USAVE(INODE,11)+F2YY(I)
      END IF
      LISTND(INODE)=LISTND(INODE)+1
  970 CONTINUE
  980 CONTINUE
C
C     END LOOP ON ELEMENTS
C
C     AVERAGE NODAL POINT FLUXES
C
      IF (ILOCE.EQ.0) GO TO 1120
      DO 1000 I=1,NUMNOD
      NN=LISTND(I)
      IF (NN.EQ.0) GO TO 1000
      DO 990 J=1,NFLUX
      USAVE(I,J)=USAVE(I,J)/NN
  990 CONTINUE
 1000 CONTINUE
      IF (IWRTP.EQ.1) GO TO 1090
C
C     OUTPUT FLUXES
C
      WRITE (NOUT, 1470) TIME
      DO 1080 N=1,NUMEL
      IF (ILOCM.EQ.1) GO TO 1020
      DO 1010 I=1,NUM
      IF (N.EQ.NEL(I)) GO TO 1020
 1010 CONTINUE
      GO TO 1080
 1020 CONTINUE
      NN=ABS(LISTEL(N))
      KIND=MOD(NN,100)
      MAT=NN/100
      ISOLID=0
      IF (PROP(15,MAT).EQ.3.0) ISOLID=1
      NSTRS=NNCOR(KIND)
      WRITE (NOUT, 1170) N,CODE(KIND)
      IF (ISOLID.EQ.1) GO TO 1040
      IF (ISTRES.EQ.1) THEN
      DO 1030 I=1,NSTRS
      INODE=ICON(N,I)
      TXX(I)=USAVE(INODE,1)
      TYY(I)=USAVE(INODE,2)
      TXY(I)=USAVE(INODE,3)
      TZZ(I)=USAVE(INODE,4)
      VORT(I)=USAVE(INODE,5)
 1030 CONTINUE
      WRITE (NOUT, 1180)
      WRITE (NOUT, 1190) (TXX(I),I=1,NSTRS)
      WRITE (NOUT, 1200) (TYY(I),I=1,NSTRS)
      WRITE (NOUT, 1210) (TXY(I),I=1,NSTRS)
      IF (IAXSYM.EQ.1) WRITE (NOUT, 1220) (TZZ(I),I=1,NSTRS)
      WRITE (NOUT, 1230) (VORT(I),I=1,NSTRS)
      END IF
C
 1040 CONTINUE
      IF (IHEAT.EQ.1) THEN
      DO 1050 I=1,NSTRS
      INODE=ICON(N,I)
      QXX(I)=USAVE(INODE,6)
      QYY(I)=USAVE(INODE,7)
 1050 CONTINUE
      WRITE (NOUT, 1290)
      WRITE (NOUT, 1300) (QXX(I),I=1,NSTRS)
      WRITE (NOUT, 1310) (QYY(I),I=1,NSTRS)
      END IF
C
      IF (IFLUX1+IFLUX2.GT.0) WRITE (NOUT, 1340)
      IF (IFLUX1.EQ.1) THEN
      DO 1060 I=1,NSTRS
      INODE=ICON(N,I)
      F1XX(I)=USAVE(INODE,8)
      F1YY(I)=USAVE(INODE,9)
 1060 CONTINUE
      WRITE (NOUT, 1350) (F1XX(I),I=1,NSTRS)
      WRITE (NOUT, 1360) (F1YY(I),I=1,NSTRS)
      END IF
C
      IF (IFLUX2.EQ.1) THEN
      DO 1070 I=1,NSTRS
      INODE=ICON(N,I)
      F2XX(I)=USAVE(INODE,10)
      F2YY(I)=USAVE(INODE,11)
 1070 CONTINUE
      WRITE (NOUT, 1390) (F2XX(I),I=1,NSTRS)
      WRITE (NOUT, 1400) (F2YY(I),I=1,NSTRS)
      END IF
 1080 CONTINUE
      GO TO 1120
C
C     GENERATE FLUX DATA AT MIDSIDE NODES AND WRITE TO TAPE
C
 1090 CONTINUE
      IF (ISTRES.EQ.1) THEN
      DO 1100 I=1,5
      CALL MIDSID (USAVE(1,I),ICON,LISTEL)
 1100 CONTINUE
      NVR=5
      WRITE (NTP8) TIME,NUMNOD,NVR,((USAVE(I,J),I=1,NUMNOD),J=1,NVR)
      END IF
      IF (IHEAT+IFLUX1+IFLUX2.EQ.0) GO TO 1120
      DO 1110 I=6,11
      CALL MIDSID (USAVE(1,I),ICON,LISTEL)
 1110 CONTINUE
      NVR=6
      WRITE (NTP7) TIME,NUMNOD,NVR,((USAVE(I,J),I=1,NUMNOD),J=6,11)
 1120 CONTINUE
C
C     CHECK FOR COMPLETION OF LOOPING
C
      IF (NSTEPS.EQ.0) GO TO 70
      IF (KOUNT.EQ.NSTEPS) GO TO 80
      KSTEP=NINC
      GO TO 170
 1130 CONTINUE
      RETURN
C
 1140 FORMAT (//,10X,'FLUX VALUES ARE CALCULATED ON THE ELEMENT BOUNDARI
     1ES' ,/,10X,'MIDWAY BETWEEN NODAL POINTS (1/4 AND 3/4 POINTS)',
     2/,10X,'REPORTED FLUXES CONSIST OF THE FOLLOWING :',/,10X,
     3'A) BASIC X,Y (OR R,Z) FLUX COMPONENTS',/,10X,'B) FLUX COMPONENT N
     4ORMAL (AND TANGENTIAL, IF APPROPRIATE) TO ELEMENT BOUNDARY',/,10X,
     5'C) INTEGRATED (NORMAL) FLUX OVER ELEMENT BOUNDARY')
 1150 FORMAT (//,10X,'FLUX VALUES ARE CALCULATED AT THE ELEMENT INTEGRAT
     1ION POINTS' ,/,10X,'(2X2 GAUSS POINTS AND 3 HAMMER POINTS)' ,
     2/,10X,'AND (BILINEARLY) EXTRAPOLATED TO THE NODES' ,/,10X,'REPORTE
     3D NODAL FLUX VALUES ARE AVERAGED FOR ADJOINING ELEMENTS',/,10X,
     4'REPORTED FLUXES CONSIST OF THE FOLLOWING :',/,10X,'A) BASIC X,Y (
     5OR R,Z) FLUX COMPONENTS')
 1160 FORMAT (///,10X,'FLUX CALCULATIONS WERE REQUESTED FOR THE FOLLOWIN
     1G ELEMENTS--' ,/,10X,20I5,/,10X,20I5,/,10X,10I5)
 1170 FORMAT (//,10X,'ELEMENT NO.' ,I4,' - ',A8)
 1180 FORMAT (/,5X,'STRESSES AND FORCES :')
 1190 FORMAT (5X,'TXX',2X,8E15.7)
 1200 FORMAT (5X,'TYY',2X,8E15.7)
 1210 FORMAT (5X,'TXY',2X,8E15.7)
 1220 FORMAT (5X,'TZZ',2X,8E15.7)
 1230 FORMAT (4X,'VORT',2X,8E15.7)
 1240 FORMAT (5X,'TN ',2X,8E15.7)
 1250 FORMAT (5X,'TS ',2X,8E15.7)
 1260 FORMAT (5X,'TNS',2X,8E15.7)
 1270 FORMAT (5X,'FN ',2X,4(E15.7,15X))
 1280 FORMAT (5X,'FS ',2X,4(E15.7,15X))
 1290 FORMAT (/,5X,'HEAT FLUXES :')
 1300 FORMAT (5X,'QX ',2X,8E15.7)
 1310 FORMAT (5X,'QY ',2X,8E15.7)
 1320 FORMAT (5X,'QN ',2X,8E15.7)
 1330 FORMAT (5X,'QA ',2X,4(E15.7,15X))
 1340 FORMAT (/,5X,'AUXILIARY VARIABLE FLUXES :')
 1350 FORMAT (5X,'F1X',2X,8E15.7)
 1360 FORMAT (5X,'F1Y',2X,8E15.7)
 1370 FORMAT (5X,'F1N',2X,8E15.7)
 1380 FORMAT (5X,'F1A',2X,4(E15.7,15X))
 1390 FORMAT (5X,'F2X',2X,8E15.7)
 1400 FORMAT (5X,'F2Y',2X,8E15.7)
 1410 FORMAT (5X,'F2N',2X,8E15.7)
 1420 FORMAT (5X,'F2A',2X,4(E15.7,15X))
 1430 FORMAT (///,10X,'FLUX CALCULATIONS WERE REQUESTED FOR ALL ELEMENTS
     1, FOR ALL TIME STEPS' ,/,10X,'FLUX VALUES ARE STORED ON TAPE FOR F
     2OR PLOTTING',/,10X,'NO TABULATED VALUES ARE PRINTED FOR THIS OPTIO
     3N')
 1440 FORMAT (///,10X,'FLUX CALCULATIONS WERE REQUESTED FOR ALL ELEMENTS
     1 IN THE MESH')
 1450 FORMAT (//,10X,'END OF FILE ENCOUNTERED ON SOLUTION FILE, PROCESSI
     1NG CONTINUES')
 1460 FORMAT (//,10X,'CALCULATIONS WERE REQUESTED FOR THE FOLLOWING QUAN
     1TITIES-' )
 1470 FORMAT (///,10X,'CURRENT TIME =' ,E15.7)
 1480 FORMAT (15X,'DEVIATORIC STRESSES',/,15X,'VORTICITY')
 1490 FORMAT (15X,'HEAT FLUX')
 1500 FORMAT (15X,'FLUX-VARIABLE 1')
 1510 FORMAT (15X,'FLUX-VARIABLE 2')
 1520 FORMAT (//,10X,'FLUX VALUES WERE REQUESTED FOR ALL TIME PLANES' )
 1530 FORMAT (//,10X,'FLUX VALUES WERE REQUESTED FOR ' ,I3,' TIME PLANES
     1' ,/,10X,'BEGINNING WITH TIME PLANE ' ,I3,' AND A TIME PLANE INCRE
     2MENT OF ' ,I3)
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


text:
subroutine CilmPlus(cilm, gridin, lmax, nmax, mass, d, rho, gridtype, w, zero, plx, n, dref)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This routine will compute the potential coefficients associated
!	with the input gridded relief using the method of Wieczorek and
!	Phillips (1998). The input grid must contain the degree-0 term
!	and the computed coefficients will be referenced to the corresonding
!	degree-0 radius. Note that the array plx is optional, and should not
!	be precomputed when memory is an issue (i.e., lmax>360).
!
!	Calling Parameters
!		IN
!			gridin:		Input grid to be transformed to spherical
!					harmonics.
!			lmax:		Maximum spherical harmonic degree to compute. For gridtype 3 and 4,
!					this must be less than or equal to N/2 - 1.
!			nmax:		Order of expansion.
!			mass:		Mass of planet.
!			rho:		density of the relief.
!			gridtype	1 = Gauss-Legendre quadrature grid corresponding to LMAX.
!					2 = N by N Driscoll and Healy grid corresponding to LMAX.
!					3 = N by 2N Driscoll and Healy grid corresponding to LMAX.
!					(4 = 2D Cartesian using MakeGrid2D is not implemented).
!		OUT
!			cilm:		Output spherical harmonic coefficients with dimensions (2, lmax+1, lmax+1).
!			d:	 	The radius that the coefficients are referenced
!					to. This parameter corresponds to the degree zero term of the data.
!		OPTIONAL
!			w:		Gauss-Legendre points used in the integrations of dimension lmax+1.
!			zero:		Array of dimension lmax+1 that contains the latitudinal
!					gridpoints used in the Gauss-Legendre quadrature integration
!					scheme. Only needed when PLX is not included.
!					(Determined from a call to PreCompute).
!			plx:		Input array of Associated Legendre Polnomials computed
!					at the Gauss points (determined from a call to
!					PreCompute). If this is not included, then the optional
!					array zero MUST be inlcuded.
!			N:		Number of latitude points in the Driscoll and Healy grid. Required for Gridtype 2 and 3
!			dref:		The reference radius used to be used when calculating the spherical
!					harmonic coefficients. If not specified, this will be set equal to the
!					mean radius of GRIDIN.
!
!	All units assumed to be SI.
!
!	Dependencies:		NGLQSH, SHExpandGLQ, SHExpandDH
!
!	Written by Mark Wieczorek 2003
!		September 3, 2005. Modifed so that the array plx is now optional.
!		April 2009. Added the option to use N by N and N by 2N Driscoll and Healy (1994) grids.
!			GRIDTYPE must now be specified, and W is now defined as an optional parameter.
!			Added the optional parameter DREF, which sets the reference radius of the expansion.
!
!	Copyright (c) 2005-2009, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	use SHTOOLS, only: NGLQSH, SHExpandGLQ, SHExpandDH

	implicit none

	real*8, intent(in) :: 	gridin(:,:), mass, rho
	real*8, intent(in), optional :: w(:), zero(:), plx(:,:), dref
	real*8, intent(out) ::	cilm(:,:,:), d
	integer, intent(in) ::	lmax, nmax, gridtype
	integer, intent(in), optional :: n
	real*8 ::		prod, pi
	real*8, allocatable ::	cilmn(:, :, :), grid(:, :)
	integer ::		j, l, k, nlat, nlong, astat(2), lmax_dh

	pi = acos(-1.0d0)

	if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax+1 .or. size(cilm(1,1,:)) < lmax+1) then
		print*, "Error --- CilmPlus"
		print*, "CILM must be dimensioned as (2, LMAX+1, LMAX+1) where LMAX is ", lmax
		print*, "Input dimension is ", size(cilm(:,1,1)), size(cilm(1,:,1)), size(cilm(1,1,:))
		stop
	endif

	if (gridtype == 4) then
		print*, "Error --- CilmPlus"
		print*, "GRIDTYPE 4 (Cartesian obtained from MakeGrid2D) is not allowed."
		stop

	elseif (gridtype == 1) then
		if (present(n)) then
			print*, "Error --- CilmPlus"
			print*, "N can not be present when using GLQ grids."
			stop
		endif

		if (present(w) .and. present(zero) .and. present(plx)) then
			print*, "Error --- CilmPlus"
			print*, "For GLQ grids, either W and ZERO or W and PLX must be present, but not all three."
			stop
		elseif (present(w) .and. present(zero)) then
			if (size(w) < lmax + 1) then
				print*, "Error --- CilmPlus"
				print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
				print*, "Input dimension is ", size(w)
				stop
			endif

			if (size(zero) < lmax + 1) then
				print*, "Error --- CilmPlus"
				print*, "ZERO must be dimensioned as (LMAX+1) where LMAX is ", lmax
				print*, "Input dimension is ", size(zero)
				stop
			endif

		elseif (present(plx) .and. present(w)) then
			if (size(w) < lmax + 1) then
				print*, "Error --- CilmPlus"
				print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
				print*, "Input dimension is ", size(w)
				stop
			endif

			if (size(plx(:,1)) < lmax+1 .or. size(plx(1,:)) < (lmax+1)*(lmax+2)/2) then
				print*, "Error --- CilmPlus"
				print*, "PLX must be dimensioned as (LMAX+1, (LMAX+1)*(LMAX+2)/2) where LMAX is ", lmax
				print*, "Input dimension is ", size(plx(:,1)), size(plx(1,:))
				stop
			endif

		else
			print*, "Error --- CilmPlus"
			print*, "For GLQ grids, either W and ZERO or W and PLX must be present"
			stop
		endif

	elseif (gridtype == 2) then
		if (present(w) .or. present(zero) .or. present(plx)) then
			print*, "Error --- CilmPlus"
			print*, "W, ZERO and PLX can not be present for Driscoll-Healy grids."
			stop
		elseif (.not.present(N)) then
			print*, "Error --- CilmPlus"
			print*, "N must be present when GRIDTYPE is 2 or 3."
			stop
		endif

	elseif (gridtype == 3) then
		if (present(w) .or. present(zero) .or. present(plx)) then
			print*, "Error --- CilmPlus"
			print*, "W, ZERO and PLX can not be present for Driscoll-Healy grids."
			stop
		elseif (.not.present(N)) then
			print*, "Error --- CilmPlus"
			print*, "N must be present when GRIDTYPE is 2 or 3."
			stop
		endif

	else
		print*, "Error --- CilmPlus"
		print*, "GRIDTYPE must be 2 (GLQ), 3 (NxN) or 4 (Nx2N)"
		print*, "Input value is ", gridtype
		stop
	endif

	if (gridtype == 1) then
		nlat = NGLQSH(lmax) !  lmax+1
		nlong = 2*lmax+1
	elseif (gridtype == 2) then
		nlat = N
		nlong = N
		lmax_dh = N/2 - 1
		if (lmax > lmax_dh) then
			print*, "Error --- CilmPlus"
			print*, "For Driscoll-Healy grids, LMAX must be less than or equal to N/2 -1, where N is ", N
			print*, "Input value of LMAX is ", lmax
			stop
		endif
	elseif (gridtype == 3) then
		nlat = N
		nlong = 2*N
		lmax_dh = N/2 - 1
		if (lmax > lmax_dh) then
			print*, "Error --- CilmPlus"
			print*, "For Driscoll-Healy grids, LMAX must be less than or equal to N/2 -1, where N is ", N
			print*, "Input value of LMAX is ", lmax
			stop
		endif
	endif

	if (size(gridin(1,:)) < nlong .or. size(gridin(:,1)) < nlat) then
		print*, "Error --- CilmPlus"
		if (gridtype == 2) then
			print*, "GRIDIN must be dimensioned as (LMAX+1, 2*LMAX+1) where LMAX is ", lmax
			print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
			stop
		elseif (gridtype == 3) then
			print*, "GRIDIN must be dimensioned as (N, N) where N is ", n
			print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
			stop
		elseif (gridtype == 4) then
			print*, "GRIDIN must be dimensioned as (N, 2N) where N is ", n
			print*, "Input dimension is ", size(gridin(1,:)), size(gridin(:,1))
			stop
		endif
	endif

	allocate(cilmn(2, lmax+1, lmax+1), stat=astat(1))
	allocate(grid(nlat, nlong), stat=astat(2))
	if (astat(1) /= 0 .or. astat(2) /= 0) then
		print*, "Error --- CilmPlus"
		print*, "Problem allocating arrays CILMN and GRID", astat(1), astat(2)
		stop
	endif

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!
	!	Do the expansion.
	!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	grid(1:nlat, 1:nlong) = gridin(1:nlat,1:nlong)
	cilm = 0.0d0
	cilmn = 0.0d0

	! Do k = 1 terms first
	select case(gridtype)
		case(1)
			if (present(plx)) then
				call SHExpandGLQ(cilmn, lmax, grid, w, plx = plx, norm = 1, csphase = 1)
			else
				call SHExpandGLQ(cilmn, lmax, grid, w, zero=zero, norm = 1, csphase = 1)
			endif
		case(2); call SHExpandDH(grid, n, cilmn, lmax_dh, norm = 1, sampling = 1, csphase = 1, lmax_calc = lmax)
		case(3); call SHExpandDH(grid, n, cilmn, lmax_dh, norm = 1, sampling = 2, csphase = 1, lmax_calc = lmax)
	end select

	if (present(dref)) then
		d = dref
	else
		d = cilmn(1,1,1)	! use mean radius of relief for reference sphere
	endif

	grid = grid - d
	cilmn(1,1,1) = cilmn(1,1,1) - d
	cilmn = cilmn/d

	do l=0, lmax
		cilm(:,l+1,1:l+1) = 4.0d0*pi*rho*(d**3)*cilmn(:,l+1,1:l+1) /mass /dble(2*l+1)
	enddo

	do k=2, nmax
		select case(gridtype)
			case(1)
				if (present(plx)) then
					call SHExpandGLQ(cilmn, lmax, (grid/d)**k, w, plx=plx, norm = 1, csphase = 1)
				else
					call SHExpandGLQ(cilmn, lmax, (grid/d)**k, w, zero=zero, norm = 1, csphase = 1)
				endif
			case(2); call SHExpandDH((grid/d)**k, n, cilmn, lmax_dh, norm = 1, sampling = 1, csphase = 1, lmax_calc = lmax)
			case(3); call SHExpandDH((grid/d)**k, n, cilmn, lmax_dh, norm = 1, sampling = 2, csphase = 1, lmax_calc = lmax)
		end select

		do l = 0, lmax
			prod = 4.0d0*pi*rho*(d**3)/mass
			do j=1, k
				prod = prod * dble(l+4-j)
			enddo
			prod = prod/( dble(2*l+1) * dble(l+3) * dble(fact(k)) )

			cilm(:,l+1,1:l+1) = cilm(:,l+1,1:l+1) + cilmn(:,l+1,1:l+1)*prod
		enddo
	enddo

	deallocate(cilmn)
	deallocate(grid)

	contains

		function fact(i)
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!
		!	This function computes the factorial of an integer
		!
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			implicit none
			integer ::	i, j, fact

			if (i == 0) then
				fact = 1
			elseif (i .lt. 0) then
				print*, "Argument to FACT must be positive"
				stop
			else
				fact = 1
				do j = 1, i
					fact = fact * j
				enddo
			endif

		end function fact

end subroutine CilmPlus

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
      SUBROUTINE SYMSLV(A,C,N,NA,NC)
C
C   PURPOSE:
C      Solve the real matrix equation A'X + XA = C, where C=C' and A
C      is n x n and in upper real Schur form.
C
C   REFERENCES:
C      Bartels, R.H.; and Stewart, G.W.: Algorithm 432 - Solution of
C        the Matrix Equation AX + XB = C.  Commun. ACM, vol. 15, no. 9,
C        Sept. 1972, pp. 820-826.
C
C   Subroutines employed by SYMSLV: SYSSLV
C   Subroutines employing SYMSLV: ATXPXA
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8
     1A(NA,1),C(NC,1),T,P
      INTEGER
     1N,NA,NC,K,KK,DK,KM1,L,LL,DL,LDL,I,IA,J,NSYS
      COMMON/SLVBLK/T(5,5),P(5),NSYS
      L = 1
   10   DL = 1
        IF(L .EQ. N) GO TO 20
        IF(A(L+1,L) .NE. 0.) DL = 2
   20   LL = L+DL-1
        K = L
   30     KM1 = K-1
          DK = 1
          IF(K .EQ. N) GO TO 35
          IF(A(K+1,K) .NE. 0.) DK = 2
   35     KK = K+DK-1
          IF(K .EQ. L) GO TO 45
          DO 40 I=K,KK
            DO 40 J=L,LL
              DO 40 IA=L,KM1
                C(I,J) = C(I,J) - A(IA,I)*C(IA,J)
   40     CONTINUE
   45     IF(DL .EQ. 2) GO TO 60
          IF(DK .EQ. 2 ) GO TO 50
          T(1,1) = A(K,K) + A(L,L)
          IF(T(1,1) .EQ. 0.) STOP
          C(K,L) = C(K,L)/T(1,1)
          GO TO 90
  50      T(1,1) = A(K,K) + A(L,L)
          T(1,2) = A(KK,K)
          T(2,1) = A(K,KK)
          T(2,2) = A(KK,KK) + A(L,L)
          P(1) = C(K,L)
          P(2) = C(KK,L)
          NSYS = 2
          CALL SYSSLV
          C(K,L) = P(1)
        C(KK,L) = P(2)
          GO TO 90
  60      IF(DK .EQ. 2) GO TO 70
          T(1,1) = A(K,K) + A(L,L)
          T(1,2) = A(LL,L)
          T(2,1) = A(L,LL)
          T(2,2) = A(K,K) + A(LL,LL)
          P(1) = C(K,L)
          P(2) = C(K,LL)
          NSYS = 2
          CALL SYSSLV
          C(K,L) = P(1)
          C(K,LL) = P(2)
          GO TO 90
  70      IF(K .NE. L) GO TO 80
          T(1,1) = A(L,L)
          T(1,2) = A(LL,L)
          T(1,3) = 0.
          T(2,1) = A(L,LL)
          T(2,2) = A(L,L) + A(LL,LL)
          T(2,3) = T(1,2)
          T(3,1) = 0.
          T(3,2) = T(2,1)
          T(3,3) = A(LL,LL)
          P(1) = C(L,L)/2.
          P(2) = C(LL,L)
          P(3) = C(LL,LL)/2.
          NSYS = 3
          CALL SYSSLV
          C(L,L) = P(1)
          C(LL,L) = P(2)
          C(L,LL) = P(2)
          C(LL,LL) = P(3)
          GO TO 90
  80      T(1,1) = A(K,K) + A(L,L)
          T(1,2) = A(KK,K)
          T(1,3) = A(LL,L)
          T(1,4) = 0.
          T(2,1) = A(K,KK)
          T(2,2) = A(KK,KK) + A(L,L)
          T(2,3) = 0.
          T(2,4) = T(1,3)
          T(3,1) = A(L,LL)
          T(3,2) = 0.
          T(3,3) = A(K,K) + A(LL,LL)
          T(3,4) = T(1,2)
          T(4,1) = 0.
          T(4,2) = T(3,1)
          T(4,3) = T(2,1)
          T(4,4) = A(KK,KK) + A(LL,LL)
          P(1) = C(K,L)
          P(2) = C(KK,L)
          P(3) = C(K,LL)
          P(4) = C(KK,LL)
          NSYS = 4
          CALL SYSSLV
          C(K,L) = P(1)
          C(KK,L) = P(2)
          C(K,LL) = P(3)
          C(KK,LL) = P(4)
  90    K = K + DK
        IF(K .LE. N) GO TO 30
        LDL = L + DL
        IF(LDL .GT. N) RETURN
        DO 120 J=LDL,N
          DO 100 I=L,LL
            C(I,J) = C(J,I)
  100     CONTINUE
          DO 120 I=J,N
            DO 110 K=L,LL
              C(I,J) = C(I,J) - C(I,K)*A(K,J) - A(K,I)*C(K,J)
  110     CONTINUE
          C(J,I) = C(I,J)
  120 CONTINUE
      L = LDL
      GO TO 10
      END

ext: for