text:
submodule (Mod_sld_BaseElmope) sldsup_BaseElmope

   contains

   !Gathers and interpolates
   module subroutine sigmaGather
      implicit none
      integer(ip) :: nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2

      call e%gather(tn,elsigma(:,:),sup%sigma(:,:,1))
   end subroutine

   module subroutine InterpolateGpSigma
      implicit none
      integer(ip) :: nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2

      call e%interpg(tn,elsigma(:,:),gpsigma(:))
   end subroutine

   module subroutine calculatePressGradientsAndStressDiv
       implicit none
       integer(ip) :: nd,tn,inode,idime,jdime
       real(rp)    :: istr(e%ndime,e%ndime,e%mnode)

       call sup_NH%Mesh%GetNdime(nd)
       tn=(nd*(nd+1))/2

       divdisp = 0.0_rp
       grdpre  = 0.0_rp
       divstr  = 0.0_rp
       grsig   = 0.0_rp

       !Pressure gradient
       call e%gradient(1,elpress(1,:),grdpre)

       !Stress divergence
       istr = 0.0_rp
       do inode=1,e%pnode
           call getStressTensor(tn,nd,elsigma(:,inode),istr(:,:,inode))
       end do
       do idime=1,e%ndime
           call e%divergence(istr(idime,:,:),divstr(idime))
       end do

       call e%divergence(eldisp(:,:,1),divdisp)

       call e%gradient(tn,elsigma,grsig)

   end subroutine calculatePressGradientsAndStressDiv

   module subroutine calculateGradients
       implicit none
       integer(ip) :: kfl_nonlinear

       call ProcHook%PhysicalProp
       call calculateGradientsAndDeter
       call calculatePressGradientsAndStressDiv

   end subroutine calculateGradients

   module subroutine AssembleStiffnessMatrixSUP
      implicit none
      integer(ip) :: pn,nd,tn,u1,uf,s1,sf,p1,bc

      call a%Mesh%GetNdime(nd)
      tn  =(nd*(nd+1))/2
      pn  = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      ! Momentum equation
      ! Assembly elmsv to elmat
      elmat(u1:uf,1:pn,s1:sf,1:pn) = elmat(u1:uf,1:pn,s1:sf,1:pn)  &
                                   + elmsv(1:nd,1:pn,1:tn,1:pn)
      ! Assembly elmpv to elmat
      elmat(u1:uf,1:pn,p1,1:pn)    = elmat(u1:uf,1:pn,p1,1:pn) &
                                   + elmpv(1:nd,1:pn,1,1:pn)

      !Constitutive equation
      ! Assembly elmut to elmat
      elmat(s1:sf,1:pn,u1:uf,1:pn) = elmat(s1:sf,1:pn,u1:uf,1:pn) &
                                   + elmut(1:tn,1:pn,1:nd,1:pn)
      ! Assembly elmst to elmat
      elmat(s1:sf,1:pn,s1:sf,1:pn) = elmat(s1:sf,1:pn,s1:sf,1:pn) &
                                   + elmst(1:tn,1:pn,1:tn,1:pn)

      ! Continuity Equation
      ! Assembly elmuq to elmat
      elmat(p1,1:pn,u1:uf,1:pn)    = elmat(p1,1:pn,u1:uf,1:pn) &
                                   + elmuq(1,1:pn,1:nd,1:pn)
      ! Assembly elmpq to elmat
      elmat(p1,1:pn,p1,1:pn)       = elmat(p1,1:pn,p1,1:pn) &
                                   + elmpq(1,1:pn,1,1:pn)

   end subroutine AssembleStiffnessMatrixSUP

   module subroutine AssembleRhsSUP
      implicit none
      integer(ip) :: pn,nd,u1,uf,s1,sf,tn,p1,bc

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      pn  = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      ! Assembly elrhu to elrhs
      elrhs(u1:uf,1:pn) = elrhs(u1:uf,1:pn) + elrhu(1:nd,1:pn)
      elrhs(s1:sf,1:pn) = elrhs(s1:sf,1:pn) + elrhd(1:tn,1:pn)
      elrhs(p1,1:pn)    = elrhs(p1,1:pn)    + elrhp(1,1:pn)

   end subroutine AssembleRhsSUP

   module subroutine AssembleRhsForcesSUP
      implicit none
      integer(ip) :: pn,nd,u1,uf,s1,sf,tn,p1,bc

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      pn  = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      ! Assembly elrhu to elrhs
      elrhs(u1:uf,1:pn) = elrhs(u1:uf,1:pn) + a%extForce(ielem)%a(u1:uf,:)
      elrhs(s1:sf,1:pn) = elrhs(s1:sf,1:pn) + a%extForce(ielem)%a(s1:sf,:)
      elrhs(p1,1:pn)    = elrhs(p1,1:pn)    + a%extForce(ielem)%a(p1,:)

   end subroutine AssembleRhsForcesSUP

   module subroutine getTauParameters_linear(tau_u,tau_s,tau_p)
      implicit none
      real(rp), intent(inout):: tau_u,tau_s,tau_p
      integer(ip)            :: pn,nd,tn,mn
      real(rp)               :: G2,cl,hl

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode
      pn  = e%pnode

      cl  = sup%kfl_PrbCharL
      hl  = e%hleng(nd)

      !Use secant shear modulus
      !It spoils convergence, so not used generally
      if(sup%kfl_useSecantModulus) then
          if(sup%devnorm .le. zeroc) then
              G2  = 2.0_rp*sup%mu
          else
              G2 = sup%signorm/sup%devnorm
          endif
      else
          G2  = 2.0_rp*sup%mu
      endif

      tau_u = sup%tau_u*((hl*hl)/(G2))
      tau_s = sup%tau_s*(hl/cl)
      tau_p = sup%tau_p*(hl/cl)

   end subroutine getTauParameters_linear

   module subroutine getTauParameters(tau_u,tau_s,tau_p)
      implicit none
      real(rp), intent(inout):: tau_u,tau_s,tau_p
      integer(ip)            :: pn,nd,tn,mn,idime
      real(rp)               :: hl,freq,tau_ut,tau_ux

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode
      pn  = e%pnode

      hl=e%hleng(1)
      do idime=2,nd
         hl=hl*e%hleng(idime)
      end do
      hl=hl**(1.0_rp/real(nd))

      freq = (sup%tau_u*sup%mu)/(hl*hl)

      if(freq>zeroc) then
         tau_ux= 1.0_rp/freq
      else
         tau_ux= 1.0e12_rp
      end if

      tau_ut  = tau_ux

      if(sup%kfl_trasg==1) then
        tau_ut = 1/(1/tau_ux + densi*a%dtinv2*tsch_deriv)
        tautau = (1.0_rp/tau_ux)*tau_ut
      endif

      tau_u  = tau_ut

      tau_s = sup%tau_s*sup%mu*2.0_rp
      tau_p = sup%tau_p*sup%mu*2.0_rp

   end subroutine getTauParameters

   module subroutine AssembleStiffnessMatrixSUP_SGS_dynsgs
      implicit none
      integer(ip)          :: pn,nd,tn,u1,uf,s1,sf,p1,bc

      call a%Mesh%GetNdime(nd)
      tn = (nd*(nd+1))/2
      pn = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      !-------------TU------------------
      ! Assembly elmvu to elmat
      elmat(u1:uf,1:pn,u1:uf,1:pn)    = elmat(u1:uf,1:pn,u1:uf,1:pn) &
                                   - elm_usgs_VU(1:nd,1:pn,1:nd,1:pn)
      ! Assembly elmvu to elmat
      elmat(u1:uf,1:pn,s1:sf,1:pn)    = elmat(u1:uf,1:pn,s1:sf,1:pn) &
                                   - elm_usgs_VS(1:nd,1:pn,1:tn,1:pn)

      ! Assembly elmqp to elmat
      elmat(u1:uf,1:pn,1,1:pn)       = elmat(u1:uf,1:pn,1,1:pn) &
                                   - elm_usgs_VP(1:nd,1:pn,1,1:pn)

   end subroutine AssembleStiffnessMatrixSUP_SGS_dynsgs

   module subroutine AssembleStiffnessMatrixSUP_SGS_dyn
      implicit none
      integer(ip)          :: pn,nd,tn,u1,uf,s1,sf,p1,bc
      real(rp)             :: tau_u,tau_s,tau_p

      call a%Mesh%GetNdime(nd)
      tn = (nd*(nd+1))/2
      pn = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      call ProcPointer%getTauParameters(tau_u,tau_s,tau_p)

      !-------------TU------------------
      ! Assembly elmsu to elmat
      elmat(s1:sf,1:pn,u1:uf,1:pn)    = elmat(s1:sf,1:pn,u1:uf,1:pn) &
                                   - tau_u*elm_usgs_EU(1:tn,1:pn,1:nd,1:pn)

      ! Assembly elmqp to elmat
      elmat(p1,1:pn,u1:uf,1:pn)       = elmat(p1,1:pn,u1:uf,1:pn) &
                                   - tau_u*elm_usgs_QU(1,1:pn,1:nd,1:pn)

   end subroutine AssembleStiffnessMatrixSUP_SGS_dyn

   module subroutine AssembleStiffnessMatrixSUP_SGS
      implicit none
      integer(ip)          :: pn,nd,tn,u1,uf,s1,sf,p1,bc
      real(rp)             :: tau_u,tau_s,tau_p

      call a%Mesh%GetNdime(nd)
      tn = (nd*(nd+1))/2
      pn = e%pnode

      call sup%GetMatrixOrganization(u1,uf,s1,sf,p1,bc)

      call ProcPointer%getTauParameters(tau_u,tau_s,tau_p)

      !-------------TU------------------
      ! Assembly elmes to elmat
      elmat(s1:sf,1:pn,s1:sf,1:pn) = elmat(s1:sf,1:pn,s1:sf,1:pn) &
                                   - tau_u*elm_usgs_ES(1:tn,1:pn,1:tn,1:pn)
      ! Assembly elmep to elmat
      elmat(s1:sf,1:pn,p1,1:pn)    = elmat(s1:sf,1:pn,p1,1:pn) &
                                   - tau_u*elm_usgs_EP(1:tn,1:pn,1,1:pn)
      ! Assembly elmqs to elmat
      elmat(p1,1:pn,s1:sf,1:pn)    = elmat(p1,1:pn,s1:sf,1:pn)  &
                                   - tau_u*elm_usgs_QS(1,1:pn,1:tn,1:pn)
      ! Assembly elmqp to elmat
      elmat(p1,1:pn,p1,1:pn)       = elmat(p1,1:pn,p1,1:pn) &
                                   - tau_u*elm_usgs_QP(1,1:pn,1,1:pn)

      !-------------TS------------------
      ! Momentum equation
      ! Assembly elmvu to elmat
      elmat(u1:uf,1:pn,u1:uf,1:pn) = elmat(u1:uf,1:pn,u1:uf,1:pn) &
                                   - tau_s*elm_ssgs_VU(1:nd,1:pn,1:nd,1:pn)
      !Assembly elmvs to elmat
      elmat(u1:uf,1:pn,s1:sf,1:pn) = elmat(u1:uf,1:pn,s1:sf,1:pn)  &
                                   - tau_s*elm_ssgs_VS(1:nd,1:pn,1:tn,1:pn)
      !Constitutive equation
      ! Assembly elmeu to elmat
      elmat(s1:sf,1:pn,u1:uf,1:pn) = elmat(s1:sf,1:pn,u1:uf,1:pn) &
                                   - tau_s*elm_ssgs_EU(1:tn,1:pn,1:nd,1:pn)
      ! Assembly elmes to elmat
      elmat(s1:sf,1:pn,s1:sf,1:pn) = elmat(s1:sf,1:pn,s1:sf,1:pn) &
                                   - tau_s*elm_ssgs_ES(1:tn,1:pn,1:tn,1:pn)

      !-------------TP------------------
      ! Momentum equation
      ! Assembly elmvu to elmat
      elmat(u1:uf,1:pn,u1:uf,1:pn) = elmat(u1:uf,1:pn,u1:uf,1:pn) &
                                   - tau_p*elm_psgs_VU(1:nd,1:pn,1:nd,1:pn)
      ! Assembly elmvp to elmat
      elmat(u1:uf,1:pn,p1,1:pn)    = elmat(u1:uf,1:pn,p1,1:pn) &
                                   - tau_p*elm_psgs_VP(1:nd,1:pn,1,1:pn)
      ! Continuity Equation
      ! Assembly elmqu to elmat
      elmat(p1,1:pn,u1:uf,1:pn)    = elmat(p1,1:pn,u1:uf,1:pn) &
                                   - tau_p*elm_psgs_QU(1,1:pn,1:nd,1:pn)
      ! Assembly elmqp to elmat
      elmat(p1,1:pn,p1,1:pn)       = elmat(p1,1:pn,p1,1:pn) &
                                   - tau_p*elm_psgs_QP(1,1:pn,1,1:pn)

      !---------------------------------END SGS----------------------------
   end subroutine AssembleStiffnessMatrixSUP_SGS

   module subroutine ResetSUPBaseArrays
      implicit none

      !-----------------Assembly-----------------
      elrhs       = 0.0_rp
      elmat       = 0.0_rp

      !-------------------RHS---------------
      elrhu       = 0.0_rp
      elrhd       = 0.0_rp
      elrhp       = 0.0_rp

      !-------------------u---------------
      elmst       = 0.0_rp
      elmsv       = 0.0_rp
      !-------------------s---------------
      elmuq       = 0.0_rp
      elmut       = 0.0_rp
      !-------------------p---------------
      elmpv       = 0.0_rp
      elmpq       = 0.0_rp

   end subroutine ResetSUPBaseArrays

   module subroutine ResetSUPBaseArrays_SGS_dynsgs
      implicit none
      !-------------------Stabilization---------------
      !-------------------Tau_u---------------
      elm_usgs_VU = 0.0_rp
      elm_usgs_VS = 0.0_rp
      elm_usgs_VP = 0.0_rp

   end subroutine ResetSUPBaseArrays_SGS_dynsgs

   module subroutine ResetSUPBaseArrays_SGS_dyn
      implicit none
      !-------------------Stabilization---------------
      !-------------------Tau_u---------------
      elm_usgs_EU = 0.0_rp
      elm_usgs_QU = 0.0_rp

   end subroutine ResetSUPBaseArrays_SGS_dyn

   module subroutine ResetSUPBaseArrays_SGS
      implicit none
      !-------------------Stabilization---------------
      !-------------------Tau_u---------------
      elm_usgs_ES = 0.0_rp
      elm_usgs_EP = 0.0_rp
      elm_usgs_QS = 0.0_rp
      elm_usgs_QP = 0.0_rp
      !-------------------Tau_s---------------
      elm_ssgs_VU = 0.0_rp
      elm_ssgs_VS = 0.0_rp
      elm_ssgs_EU = 0.0_rp
      elm_ssgs_ES = 0.0_rp
      !-------------------Tau_p---------------
      elm_psgs_VU = 0.0_rp
      elm_psgs_VP = 0.0_rp
      elm_psgs_QU = 0.0_rp
      elm_psgs_QP = 0.0_rp

   end subroutine ResetSUPBaseArrays_SGS

   module subroutine ResetSUPNonlinearDerTerms
      implicit none
      !-------------------Stabilization---------------
      !-------------------Tau_u---------------
      derJFki = 0.0_rp
      grdFmat = 0.0_rp
      grdFspt = 0.0_rp
      divFmat = 0.0_rp
      divFspt = 0.0_rp

   end subroutine ResetSUPNonlinearDerTerms

   module subroutine AllocateBaseElmopeMatricesSUP_SGS_dynsgs
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Displ stabilization
      call sup%Memor%alloc(nd,mn,nd,mn,elm_usgs_VU,'elm_usgs_VU','sldsup_elmope')
      !-----------Stress stabilization
      call sup%Memor%alloc(nd,mn,tn,mn,elm_usgs_VS,'elm_usgs_VS','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%alloc(nd,mn,1 ,mn,elm_usgs_VP,'elm_usgs_VP','sldsup_elmope')

   end subroutine

   module subroutine DeallocateBaseElmopeMatricesSUP_SGS_dynsgs
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Displ stabilization
      call sup%Memor%dealloc(nd,mn,nd,mn,elm_usgs_VU,'elm_usgs_VU','sldsup_elmope')
      !-----------Stress stabilization
      call sup%Memor%dealloc(nd,mn,tn,mn,elm_usgs_VS,'elm_usgs_VS','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%dealloc(nd,mn,1 ,mn,elm_usgs_VP,'elm_usgs_VP','sldsup_elmope')

   end subroutine

   module subroutine AllocateBaseElmopeMatricesSUP_SGS_dyn
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Stress stabilization
      call sup%Memor%alloc(tn,mn,nd,mn,elm_usgs_EU,'elm_usgs_EU','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%alloc(1 ,mn,nd,mn,elm_usgs_QU,'elm_usgs_QU','sldsup_elmope')

   end subroutine

   module subroutine DeallocateBaseElmopeMatricesSUP_SGS_dyn
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Stress stabilization
      call sup%Memor%dealloc(tn,mn,nd,mn,elm_usgs_EU,'elm_usgs_EU','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%dealloc(1 ,mn,nd,mn,elm_usgs_QU,'elm_usgs_QU','sldsup_elmope')

   end subroutine

   module subroutine AllocateBaseElmopeMatricesSUP_SGS
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Stress stabilization
      call sup%Memor%alloc(tn,mn,tn,mn,elm_usgs_ES,'elm_usgs_ES','sldsup_elmope')
      call sup%Memor%alloc(tn,mn,1 ,mn,elm_usgs_EP,'elm_usgs_EP','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%alloc(1 ,mn,tn,mn,elm_usgs_QS,'elm_usgs_QS','sldsup_elmope')
      call sup%Memor%alloc(1 ,mn,1 ,mn,elm_usgs_QP,'elm_usgs_QP','sldsup_elmope')

      !-------------------Tau_s matrix calculation----------
      !-----------Displacement stabilization
      call sup%Memor%alloc(nd,mn,nd,mn,elm_ssgs_VU,'elm_ssgs_VU','sldsup_elmope')
      call sup%Memor%alloc(nd,mn,tn,mn,elm_ssgs_VS,'elm_ssgs_VS','sldsup_elmope')
      !-----------Stress stabilization
      call sup%Memor%alloc(tn,mn,tn,mn,elm_ssgs_ES,'elm_ssgs_ES','sldsup_elmope')
      call sup%Memor%alloc(tn,mn,nd,mn,elm_ssgs_EU,'elm_ssgs_EU','sldsup_elmope')

      !-------------------Tau_p matrix calculation----------
      !-----------Displacement stabilization
      call sup%Memor%alloc(nd,mn,nd,mn,elm_psgs_VU,'elm_psgs_VU','sldsup_elmope')
      call sup%Memor%alloc(nd,mn,1 ,mn,elm_psgs_VP,'elm_psgs_VP','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%alloc(1 ,mn,nd,mn,elm_psgs_QU,'elm_psgs_QU','sldsup_elmope')
      call sup%Memor%alloc(1 ,mn,1 ,mn,elm_psgs_QP,'elm_psgs_QP','sldsup_elmope')

   end subroutine

   module subroutine DeallocateBaseElmopeMatricesSUP_SGS
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      !-------------------Tau_u matrix calculation----------
      !-----------Stress stabilization
      call sup%Memor%dealloc(tn,mn,tn,mn,elm_usgs_ES,'elm_usgs_ES','sldsup_elmope')
      call sup%Memor%dealloc(tn,mn,1 ,mn,elm_usgs_EP,'elm_usgs_EP','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%dealloc(1 ,mn,tn,mn,elm_usgs_QS,'elm_usgs_QS','sldsup_elmope')
      call sup%Memor%dealloc(1 ,mn,1 ,mn,elm_usgs_QP,'elm_usgs_QP','sldsup_elmope')

      !-------------------Tau_s matrix calculation----------
      !-----------Displacement stabilization
      call sup%Memor%dealloc(nd,mn,nd,mn,elm_ssgs_VU,'elm_ssgs_VU','sldsup_elmope')
      call sup%Memor%dealloc(nd,mn,tn,mn,elm_ssgs_VS,'elm_ssgs_VS','sldsup_elmope')
      !-----------Stress stabilization
      call sup%Memor%dealloc(tn,mn,tn,mn,elm_ssgs_ES,'elm_ssgs_ES','sldsup_elmope')
      call sup%Memor%dealloc(tn,mn,nd,mn,elm_ssgs_EU,'elm_ssgs_EU','sldsup_elmope')

      !-------------------Tau_p matrix calculation----------
      !-----------Displacement stabilization
      call sup%Memor%dealloc(nd,mn,nd,mn,elm_psgs_VU,'elm_psgs_VU','sldsup_elmope')
      call sup%Memor%dealloc(nd,mn,1 ,mn,elm_psgs_VP,'elm_psgs_VP','sldsup_elmope')
      !-----------Pressure stabilization
      call sup%Memor%dealloc(1 ,mn,nd,mn,elm_psgs_QU,'elm_psgs_QU','sldsup_elmope')
      call sup%Memor%dealloc(1 ,mn,1 ,mn,elm_psgs_QP,'elm_psgs_QP','sldsup_elmope')

   end subroutine

   module subroutine AllocateSUPAssemblyMatrices
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn  =(nd*(nd+1))/2
      mn  = e%mnode

      call sup%Memor%alloc(tn,mn,      elrhd,'elrhd','sldsup_Elmope')
      call sup%Memor%alloc(1,mn,       elrhp,'elrhp','sldsup_Elmope')
      !-----------Displacement testing
      call sup%Memor%alloc(nd,mn,tn,mn,elmsv,'elmsv','sldsup_elmope')
      call sup%Memor%alloc(nd,mn,1 ,mn,elmpv,'elmpv','sldsup_elmope')
      !-----------Stress testing
      call sup%Memor%alloc(tn,mn,tn,mn,elmst,'elmst','sldsup_elmope')
      call sup%Memor%alloc(tn,mn,nd,mn,elmut,'elmut','sldsup_elmope')
      !-----------Pressure testing
      call sup%Memor%alloc(1 ,mn,nd,mn,elmuq,'elmuq','sldsup_elmope')
      call sup%Memor%alloc(1 ,mn,1 ,mn,elmpq,'elmpq','sldsup_elmope')

   end subroutine

   module subroutine DeallocateSUPAssemblyMatrices
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn  =(nd*(nd+1))/2
      mn  = e%mnode

      call sup%Memor%dealloc(tn,mn,      elrhd,'elrhd','sldsup_Elmope')
      call sup%Memor%dealloc(1,mn,       elrhp,'elrhp','sldsup_Elmope')
      !-----------Displacement testing
      call sup%Memor%dealloc(nd,mn,tn,mn,elmsv,'elmsv','sldsup_elmope')
      call sup%Memor%dealloc(nd,mn,1 ,mn,elmpv,'elmpv','sldsup_elmope')
      !-----------Stress testing
      call sup%Memor%dealloc(tn,mn,tn,mn,elmst,'elmst','sldsup_elmope')
      call sup%Memor%dealloc(tn,mn,nd,mn,elmut,'elmut','sldsup_elmope')
      !-----------Pressure testing
      call sup%Memor%dealloc(1 ,mn,nd,mn,elmuq,'elmuq','sldsup_elmope')
      call sup%Memor%dealloc(1 ,mn,1 ,mn,elmpq,'elmpq','sldsup_elmope')

   end subroutine

   module subroutine AllocateBaseElmopeMatricesSUPLinear
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn  =(nd*(nd+1))/2
      mn  = e%mnode

      call sup%Memor%alloc(tn,tn,      P,      'P','sldsup_Elmope')
      call sup%Memor%alloc(tn,tn,      D,      'D','sldsup_Elmope')
      call sup%Memor%alloc(tn,tn,  C_dev,  'C_dev','sldsup_Elmope')
      call sup%Memor%alloc(tn,tn,  D_dev,  'D_dev','sldsup_Elmope')
   end subroutine

   module subroutine DeallocateBaseElmopeMatricesSUPLinear
      implicit none
      integer(ip) :: mn,nd,tn

      call a%Mesh%GetNdime(nd)
      tn=(nd*(nd+1))/2
      mn  = e%mnode

      call sup%Memor%dealloc(tn,tn,      P,      'P','sldsup_Elmope')
      call sup%Memor%dealloc(tn,tn,      D,      'D','sldsup_Elmope')
      call sup%Memor%dealloc(tn,tn,  C_dev,  'C_dev','sldsup_Elmope')
      call sup%Memor%dealloc(tn,tn,  D_dev,  'D_dev','sldsup_Elmope')
   end subroutine

end submodule sldsup_BaseElmope

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
C
C     Driver routine for REACTOR code
C
      PROGRAM DRIVER

      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      PARAMETER (LR=2000000, LI=500000, LC=5000)
      PARAMETER (LRDSL=500000, LIDSL=200000)
      PARAMETER (LIN=9,LPAR=10,LDAT=11,LINKCK=12,LOUT=6)
      DIMENSION RWORK(LR), IWORK(LI) !, TEMP(10), PRES(10)
      CHARACTER CWORK(LC)*16
	LOGICAL XML

      COMMON /ICONS/  NSPC,NRXN,NVAR,NT,NY,NP
      COMMON /IPOINT/ NIPAR,NICK,NIDSL
      COMMON /RPOINT/ NRPAR,NREAC,NWT,NU,NWDOT,NCVMS,NC,NQ,NDQ,NRCK,
     >                NRDSL

      DATA RWORK/LR*0.0D0/, IWORK/LI*0/, XML/.TRUE./

C     Open input and output files
	IF(XML) THEN
		OPEN (LIN,    STATUS='UNKNOWN',FILE='reactorInput.xml')
		OPEN (LDAT,   STATUS='UNKNOWN',FILE='reactorOutput.xml')
	ELSE
		OPEN (LIN,    STATUS='UNKNOWN',FILE='reactor.inp')
	    OPEN (LDAT,   STATUS='UNKNOWN',FILE='reactor.out')
	ENDIF

      OPEN (LINKCK, FORM='UNFORMATTED',STATUS='UNKNOWN',FILE='chem.bin')
      OPEN (LPAR,   STATUS='UNKNOWN',FILE='parameter.inp')
      OPEN (LOUT,   STATUS='UNKNOWN',FILE='diagnostics.out')

C     Get size of workspaces
      CALL CKLEN(LINKCK, LOUT, LENI, LENR, LENC)

C     Determine mechanism size
      IF (LENI.LT.LI .AND. LENR.LT.LR .AND. LENC.LT.LC) THEN
         CALL CKINIT(LI, LR, LC, LINKCK, LOUT, IWORK, RWORK, CWORK)
         CALL CKINDX(IWORK, RWORK, NELM, NSPC, NRXN, NFIT)
      END IF

C     Number of states and number of equations
      NVAR = NSPC + 2
      NEQ  = NVAR

C     Pointers to Q,DQ arrays
      NT    = 1               ! Pointer to temperature variable
      NY    = NT + 1          ! Pointer to mass fraction variable
      NP    = NY + NSPC       ! Pointer to pressure variable

C     Apportion integer workspace
      NNPAR = 1               ! Number of integer parameters (NRXN)
      NIPAR = NNPAR + 1       ! Integer parameters (point to position of A factors)
      NICK  = NIPAR + NRXN    ! Pointer to integer chemkin workspace
      NIDSL = NICK  + LENI    ! DSL48S integer work space
      LITOT = NIDSL + LIDSL-1 ! Length of IWORK

C     Apportion real workspace
      NRPAR = 1               ! Real parameters (A factors)
      NREAC = NRPAR + NRXN    ! Pointer to REAC array
      NWT   = NREAC + NSPC    ! Pointer to molecular weights
      NU    = NWT   + NSPC    ! Pointer to internal energies
      NWDOT = NU    + NSPC    ! Pointer to molar production rates
      NCVMS = NWDOT + NSPC    ! Pointer to Cv's
      NC    = NCVMS + NSPC    ! Pointer to concentrations
      NQ    = NC    + NSPC    ! Pointer to states and sensitivities
      NDQ   = NQ    + NEQ     ! Pointer to derivatives of states and sens
      NRCK  = NDQ   + NEQ     ! Pointer to real chemkin workspace
      NRDSL = NRCK  + LENR    ! DSL48S real work space
      LRTOT = NRDSL + LRDSL-1 ! Length of RWORK

C     Apportion character workspace
      NCCK  = 1              ! Pointer to character chemkin workspace
      NKSYM = NCCK  + LENC   ! Pointer to species names
      LCTOT = NKSYM + NSPC-1 ! Length of CWORK

C     Check if the alloted workspace is large enough
      IF (LITOT.GT.LI .OR. LRTOT.GT.LR .OR. LCTOT.GT.LC) THEN

         WRITE (LOUT, 7000) LI, LITOT, LR, LRTOT, LC, LCTOT
         WRITE (LOUT, '(A)') 'ERROR: NOT ENOUGH WORKSPACE PROVIDED'

 7000    FORMAT (/,'                WORKING SPACE REQUIREMENTS',
     1           /,'                 PROVIDED        REQUIRED ',
     2           /,' INTEGER  ' , 2I15,
     3           /,' REAL     ' , 2I15,
     4           /,' CHARACTER' , 2I15,/)

	   CALL ERRORXML(LDAT,'ERROR: NOT ENOUGH WORKSPACE PROVIDED')
         STOP
      END IF

C     Initialize CHEMKIN workspace
      DO I=1,LR
         RWORK(I) = 0.0D0
      END DO

      DO I=1,LI
         IWORK(I) = 0
      END DO

      REWIND(LINKCK)
      CALL CKINIT(LI, LR, LC, LINKCK, LOUT, IWORK(NICK), RWORK(NRCK),
     >           CWORK(NCCK))
      CLOSE(LINKCK)

c	begin test
c	do i=1,10
c		temp(i) = 300.0 + float(i-1)*20.0
c		pres(i) = 1.01325e6*2 + float(i-1)*1.01325e6
c	enddo
c	rwork(nc) = 1.0d0
c	do i=1,10
c		do j=1,10
c			call ckwyp(pres(i), temp(j), rwork(nc), iwork(nick),
c     >                   rwork(nrck), rwork(nwdot))
c	    enddo
c      enddo
c	stop
c	end test

C     Call reactor subroutine
      CALL REACTOR(NEQ,LIN,LPAR,LDAT,LOUT,RWORK,IWORK,CWORK,
     >             RWORK(NQ),RWORK(NDQ),CWORK(NKSYM),RWORK(NREAC),
     >             LRDSL,LIDSL,XML)

      CLOSE(LOUT)
	CLOSE(LDAT)

      STOP
      END

ext: f

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
C    %W% %G%
      subroutine wgvcal(to,vfdint)
C * * *
C * * * THIS SUBROUTINE CALCULATES GENERATOR FIELD VOLTAGE DEVIATION
C * * * AT EACH TIME STEP AND SAVES THE LOWEST 20 DEVIATIONS,
C * * * BUS NAMES WHERE THEY OCCUR, AND THE TIME WHEN THEY
C * * * OCCURRED.  IT IS CALLED BY NOUT2.
C * * *
      include 'tspinc/params.inc'
      include 'tspinc/comn56.inc'
      include 'tspinc/link56.inc'
      include 'tspinc/room.inc'
      include 'tspinc/wstequ.inc'
      include 'tspinc/wgv.inc'
      include 'tspinc/nwgntn.inc'
      dimension vfdint(MAXGEN)
      do 500 i = 1,isg
      if(vfldt(1,i) .eq. 0.0)go to 500
      vltdev = abs(vfldt(1,i) - vfdint(i))
      ibkv = nwgntn(i)
      do 200 j = 1,20
      if(nwgntc(1,i).eq.wgvnam(j).and.basekv(ibkv) .eq. wgvkv(j)
     1  .and. nwgntc(2,i) .eq. wgvid(j))then
          if(vltdev .gt. wgvval(j))then
            if(j .eq. 20)go to 250
            do 175 m = j,19
            wgvval(m) = wgvval(m+1)
            wgvnam(m) = wgvnam(m+1)
            wgvid(m) = wgvid(m+1)
            wgvkv(m) = wgvkv(m+1)
 175        wgvtim(m) = wgvtim(m+1)
            go to 250
         endif
         go to 500
      endif
 200  continue
 250  do 400 j= 1,20
      if(vltdev .gt. wgvval(j))then
         l = 21
         do 300 k = j+1,20
         l = l - 1
         wgvval(l) = wgvval(l-1)
         wgvtim(l) = wgvtim(l-1)
         wgvnam(l) = wgvnam(l-1)
         wgvid(l) = wgvid(l-1)
 300     wgvkv(l) = wgvkv(l-1)
         wgvval(j) = vltdev
         wgvtim(j) = to
         wgvnam(j) = nwgntc(1,i)
         wgvid(j) = nwgntc(2,i)
         wgvkv(j) = basekv(ibkv)
         go to 500
      endif
 400  continue
 500  continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


SUBROUTINE phi_linear_3D_standing(H,angle,Lx,Ly,xc,yc,z,d,omega,time,grav,phix,phiy,phiz,phi,phit)
!
USE Precision
USE Constants
IMPLICIT NONE
!
REAL(KIND=long), INTENT(IN)  :: H,Lx,Ly,xc,yc,z,d,omega,time,grav
REAL(KIND=long), INTENT(OUT) :: phix,phiy,phiz
REAL(KIND=long), INTENT(OUT), OPTIONAL :: phi,phit
! REAL(KIND=long), OPTIONAL :: y,phi,phiyy !Boundary fitted ?
! Local variables
REAL(KIND=long) :: kx, ky, k_omegat, sink_omegat, cosk_omegat, temp1, temp2
!
REAL(KIND=long) :: HH,k,w,g,c,angle,x,y,tmp
!
x = xc*COS(-angle)-yc*SIN(-angle)
y = xc*SIN(-angle)+yc*COS(-angle)

g = grav
kx    = two*pi/Lx; ky = two*pi/Ly
k = SQRT(kx**2+ky**2)
w = omega
c = SQRT(g/k*TANH(k*d))
HH = H

phix    =  kx*half*HH*c*COSH(k*(z+d))/SINH(k*d)*SIN(w*time)*SIN(kx*X)*COS(ky*Y)
phiy    =  ky*half*HH*c*COSH(k*(z+d))/SINH(k*d)*SIN(w*time)*COS(kx*X)*SIN(ky*Y)
phiz    = -half*k*HH*c*SINH(k*(z+d))/SINH(k*d)*SIN(w*time)*COS(kx*X)*COS(ky*Y)
!
IF(present(phit)) THEN
    phit   = -half*w*HH*c*COSH(k*(z+d))/SINH(k*d)*COS(w*time)*COS(kx*X)*COS(ky*Y)
ENDIF
IF(present(phi)) THEN
    phi    = -half*HH*c*COSH(k*(z+d))/SINH(k*d)*SIN(w*time)*COS(kx*X)*COS(ky*Y)
ENDIF
!
tmp = phix
phix = tmp*COS(angle)-phiy*SIN(angle)
phiy = tmp*SIN(angle)+phiy*COS(angle)
!
!xc = x*COS(angle)-y*SIN(angle)
!yc = x*SIN(angle)+y*COS(angle)

END SUBROUTINE phi_linear_3D_standing


text:
Describe Users/lreyes here.

20080317 18:22:24 nbsp Welcome to the Wiki. There is also a Wiki for wiki:Winters. You might want to list your business there as well. Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program main

integer f(30)
do i=1,30
f(i)=1
enddo

do i=1,3
call update()
print*,'global f',f
enddo
contains

subroutine update()

use mpi

integer i,j,istr,iend,avgtask,extra
integer f(30)
integer, allocatable :: lf(:),kf(:)
integer code,np,world,ierr
character*8 name
!================================
!do i=1,30
!   do j=1,30
!     f(i,j)=i+j
!  enddo
!enddo

!do i=1,30
!   do j=1,30
!        g(i,j)=i*j
!          enddo
!          enddo

!============================

!single
!do i=1,3
! do j=1,3
 !  k(i,j)=f(i,j)+g(i,j)
 ! enddo
!enddo

!print*,k
!=============================

!parallel
call mpi_init(ierr)
call mpi_comm_dup(mpi_comm_world,world,ierr)
call mpi_comm_rank(world,code,ierr)
call mpi_comm_size(world,np,ierr)

status=hostnm(name)
!write(*,*)'i can',code
!write(*,*)'core number:',code,'from','  ',name

!lf=f(1:3,:)
!lg=g(1:3,:)
!lk=lf+lg
!print*,'fuck'
avgtask=30/np
extra=mod(30,np)
if(code .le. extra)then
        istr=(avgtask+1)*code+1
        iend=(avgtask+1)*(code+1)
else
        istr=avgtask*code+1+extra
        iend=avgtask*(code+1)+extra
endif
allocate(lf(iend-istr+1))
if(code.eq.0)then
allocate(kf(30))
endif

do i=1,iend-istr+1
lf(i)=10*code
enddo
call mpi_gather(lf,iend-istr+1,mpi_integer,kf,iend-istr+1,mpi_integer,0,world,ierr)

if(code.eq.0)then
print*,'code:',code,'kf:',kf
do i=1,30
f(i)=kf(i)
enddo
endif

deallocate(lf)

!print*,'Glowbal f',f
call mpi_comm_free(world,ierr)
call mpi_finalize(ierr)
end subroutine update

end program

max_stars_repo_name: CaptainYin/ParLTRANS

max_issues_repo_name: CaptainYin/ParLTRANS

max_forks_repo_name: CaptainYin/ParLTRANS


text:
MODULE define_dimensions
  INTEGER            :: mp      ! # points in this simulation
  INTEGER, PARAMETER :: mf = 2  ! # leaves (sunlit, shaded)
  INTEGER, PARAMETER :: nrb = 3 ! # radiation bands
  INTEGER, PARAMETER :: ms = 6  ! # soil layers
  INTEGER, PARAMETER :: ncp = 3 ! # vegetation carbon stores
  INTEGER, PARAMETER :: ncs = 2 ! # soil carbon stores
END MODULE define_dimensions

MODULE define_types
  ! i_d is default kind for representing integer values.
!  INTEGER, PARAMETER :: i_d = KIND(9)
  ! r_1 is default kind for representing REAL values (typically 32 or 64 bits).
!  INTEGER, PARAMETER :: r_1  = KIND(1.0)
  ! r_2 is kind for representing REAL values with at least 10-digit precision
  ! (typically 64 bits).
!  INTEGER, PARAMETER :: r_2  = SELECTED_REAL_KIND(12, 50)

 TYPE veg_parameter_type
      integer, dimension(:), pointer   :: iveg
      real,    dimension(:,:), pointer :: froot
 END type veg_parameter_type

 TYPE soil_parameter_type
      integer, dimension(:), pointer  :: isoilm
      real,    dimension(:), pointer  :: sfc, swilt
      real,    dimension(:), pointer  :: clay, silt, sand, ssat
      real,    dimension(:), pointer  :: dzsoil
 END type soil_parameter_type

 TYPE (veg_parameter_type)  :: veg
 TYPE (soil_parameter_type) :: soil

Contains

   SUBROUTINE alloc_casavegsoil(mp,ms)
     IMPLICIT none

! TYPE (veg_parameter_type)  :: veg
! TYPE (soil_parameter_type) :: soil

     integer, INTENT(IN) :: mp,ms

     ALLOCATE(veg%iveg(mp))
     ALLOCATE(veg%froot(mp,ms))
     ALLOCATE(soil%isoilm(mp))
     ALLOCATE(soil%sfc(mp), soil%swilt(mp), soil%ssat(mp))
     ALLOCATE(soil%clay(mp), soil%silt(mp), soil%sand(mp))
     ALLOCATE(soil%dzsoil(ms))
   END subroutine alloc_casavegsoil

END MODULE define_types

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
Cast forth a twelve times table, suitable for chanting at school.
      INTEGER I,J	!Steppers.
      CHARACTER*52 ALINE	!Scratchpad.
      WRITE(6,1) (I,I = 1,12)	!Present the heading.
    1 FORMAT ("  ×|",12I4,/," --+",12("----"))	!Alas, can't do overprinting with underlines now.
      DO 3 I = 1,12		!Step down the lines.
        WRITE (ALINE,2) I,(I*J, J = 1,12)	!Prepare one line.
    2   FORMAT (I3,"|",12I4)		!Aligned with the heading.
        ALINE(5:1 + 4*I) = ""		!Scrub the unwanted part.
    3   WRITE (6,"(A)") ALINE		!Print the text.
      END	!"One one is one! One two is two! One three is three!...

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      SUBROUTINE MB03BC( K, AMAP, S, SINV, A, LDA1, LDA2, MACPAR, CV,
     $                   SV, DWORK )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute the product singular value decomposition of the K-1
C     triangular factors corresponding to a 2-by-2 product of K
C     factors in upper Hessenberg-triangular form.
C     For a general product of 2-by-2 triangular matrices
C
C                        S(2)        S(3)            S(K)
C            A = A(:,:,2)    A(:,:,3)    ... A(:,:,K),
C
C     Givens rotations are computed so that
C                                                          S(i)
C       [  CV(i-1) SV(i-1) ] [ A(1,1,i)(in)  A(1,2,i)(in) ]
C       [ -SV(i-1) CV(i-1) ] [     0         A(2,2,i)(in) ]
C                                      S(i)
C       [ A(1,1,i)(out) A(1,2,i)(out) ]    [  CV(i) SV(i) ]
C     = [     0         A(2,2,i)(out) ]    [ -SV(i) CV(i) ]
C
C     stays upper triangular and
C
C       [  CV(1) SV(1) ]       [ CV(K) -SV(K) ]
C       [ -SV(1) CV(1) ] * A * [ SV(K)  CV(K) ]
C
C     is diagonal.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     K       (input)  INTEGER
C             The number of factors.  K >= 1.
C
C     AMAP    (input) INTEGER array, dimension (K)
C             The map for accessing the factors, i.e., if AMAP(I) = J,
C             then the factor A_I is stored at the J-th position in A.
C
C     S       (input)  INTEGER array, dimension (K)
C             The signature array. Each entry of S must be 1 or -1.
C
C     SINV    (input) INTEGER
C             Signature multiplier. Entries of S are virtually
C             multiplied by SINV.
C
C     A       (input/output)  DOUBLE PRECISION array, dimension
C                             (LDA1,LDA2,K)
C             On entry, the leading 2-by-2-by-K part of this array must
C             contain a 2-by-2 product (implicitly represented by its K
C             factors) in upper Hessenberg-triangular form.
C             On exit, the leading 2-by-2-by-K part of this array
C             contains modified triangular factors such that their
C             product is diagonal.
C
C     LDA1    INTEGER
C             The first leading dimension of the array A.  LDA1 >= 2.
C
C     LDA2    INTEGER
C             The second leading dimension of the array A.  LDA2 >= 2.
C
C     MACPAR  (input)  DOUBLE PRECISION array, dimension (5)
C             Machine parameters:
C             MACPAR(1)  overflow threshold,         DLAMCH( 'O' );
C             MACPAR(2)  underflow threshold,        DLAMCH( 'U' );
C             MACPAR(3)  safe minimum,               DLAMCH( 'S' );
C             MACPAR(4)  relative machine precision, DLAMCH( 'E' );
C             MACPAR(5)  base of the machine,        DLAMCH( 'B' ).
C
C     CV      (output)  DOUBLE PRECISION array, dimension (K)
C             On exit, the first K elements of this array contain the
C             cosines of the Givens rotations.
C
C     SV      (output)  DOUBLE PRECISION array, dimension (K)
C             On exit, the first K elements of this array contain the
C             sines of the Givens rotations.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (3*(K-1))
C
C     METHOD
C
C     The product singular value decomposition of the K-1
C     triangular factors are computed as described in [1].
C
C     REFERENCES
C
C     [1] Bojanczyk, A. and Van Dooren, P.
C         On propagating orthogonal transformations in a product of 2x2
C         triangular matrices.
C         In Reichel, Ruttan and Varga: 'Numerical Linear Algebra',
C         pp. 1-9, 1993.
C
C     CONTRIBUTOR
C
C     D. Kressner, Technical Univ. Berlin, Germany, June 2001.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
C     July 2009, SLICOT Library version of the routine PLAPST.
C     V. Sima, Nov. 2010, Aug. 2011.
C
C     KEYWORDS
C
C     Eigenvalues, orthogonal transformation, singular values,
C     singular value decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE, HALF, TWO
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0,
     $                    TWO = 2.0D0 )
C     .. Scalar Arguments ..
      INTEGER           K, LDA1, LDA2, SINV
C     .. Array Arguments ..
      INTEGER           AMAP(*), S(*)
      DOUBLE PRECISION  A(LDA1,LDA2,*), CV(*), DWORK(*), MACPAR(*),
     $                  SV(*)
C     .. Local Scalars ..
      INTEGER           AI, I, PW, SCL
      DOUBLE PRECISION  A11, A12, A22, B11, B12, B22, BASE, CC, CL, CR,
     $                  EPS, MX, MX2, RMAX, RMIN, RMNS, RMXS, S11, S22,
     $                  SC, SFMN, SL, SR, SSMAX, SSMIN, T11, T12, T22,
     $                  TEMP, TWOS
C     .. External Subroutines ..
      EXTERNAL          DLARTG, DLASV2
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SQRT
C
C     .. Executable Statements ..
C
      RMAX = MACPAR(1)
      RMXS = SQRT( RMAX )
      RMIN = MACPAR(2)
      RMNS = SQRT( RMIN )
      SFMN = MACPAR(3)
      EPS  = MACPAR(4)
      BASE = MACPAR(5)
      TWOS = SQRT( TWO )
C
C     Compute the product of the 2-by-2 triangular matrices.
C
      PW  = 1
      T11 = ONE
      T12 = ZERO
      T22 = ONE
C
      DO 60  I = 2, K
         AI  = AMAP(I)
         A11 = A(1,1,AI)
         A12 = A(1,2,AI)
         A22 = A(2,2,AI)
         IF ( S(AI).NE.SINV ) THEN
            TEMP =  A11
            A11  =  A22
            A22  =  TEMP
            A12  = -A12
         END IF
C
C        A and T are scaled so that the elements of the resulting
C        product do not overflow.
C
         MX  = ABS( A11 ) / RMXS
         MX2 = ABS( T11 ) / RMXS
   10    CONTINUE
         IF ( MX*MX2.GE.ONE ) THEN
            IF (  MX.GE.ONE ) THEN
               MX  =  MX / BASE
               A11 = A11 / BASE
               A22 = A22 / BASE
               A12 = A12 / BASE
            END IF
            IF ( MX2.GE.ONE ) THEN
               MX2 = MX2 / BASE
               T11 = T11 / BASE
               T22 = T22 / BASE
               T12 = T12 / BASE
            END IF
            GOTO 10
         END IF
C
         MX  = ABS( A22 ) / RMXS
         MX2 = ABS( T22 ) / RMXS
   20    CONTINUE
         IF ( MX*MX2.GE.ONE ) THEN
            IF (  MX.GE.ONE ) THEN
               MX  = MX  / BASE
               A11 = A11 / BASE
               A22 = A22 / BASE
               A12 = A12 / BASE
            END IF
            IF ( MX2.GE.ONE ) THEN
               MX2 = MX2 / BASE
               T11 = T11 / BASE
               T22 = T22 / BASE
               T12 = T12 / BASE
            END IF
            GOTO 20
         END IF
C
         MX  = ABS( A12 ) / RMXS
         MX2 = ABS( T11 ) / RMXS
   30    CONTINUE
         IF ( MX*MX2.GE.HALF ) THEN
            IF (  MX.GE.HALF ) THEN
               MX  = MX  / BASE
               A11 = A11 / BASE
               A22 = A22 / BASE
               A12 = A12 / BASE
            END IF
            IF ( MX2.GE.HALF ) THEN
               MX2 = MX2 / BASE
               T11 = T11 / BASE
               T22 = T22 / BASE
               T12 = T12 / BASE
            END IF
            GOTO 30
         END IF
C
         MX  = ABS( A22 ) / RMXS
         MX2 = ABS( T12 ) / RMXS
   40    CONTINUE
         IF ( MX*MX2.GE.HALF ) THEN
            IF (  MX.GE.HALF ) THEN
               MX  = MX  / BASE
               A11 = A11 / BASE
               A22 = A22 / BASE
               A12 = A12 / BASE
            END IF
            IF ( MX2.GE.HALF ) THEN
               MX2 = MX2 / BASE
               T11 = T11 / BASE
               T22 = T22 / BASE
               T12 = T12 / BASE
            END IF
            GOTO 40
         END IF
C
C        Avoid underflow if possible.
C
         MX  = MAX( ABS( A11 ), ABS( A22 ), ABS( A12 ) )
         MX2 = MAX( ABS( T11 ), ABS( T22 ), ABS( T12 ) )
         IF ( MX.NE.ZERO .AND. MX2.NE.ZERO ) THEN
   50       CONTINUE
            IF ( ( MX.LE.( ONE/RMNS ) .AND. MX2.LE.RMNS ) .OR.
     $           ( MX.LE.RMNS .AND. MX2.LE.( ONE/RMNS ) ) )
     $      THEN
               IF ( MX.LE.MX2 ) THEN
                  MX  = MX  * BASE
                  A11 = A11 * BASE
                  A22 = A22 * BASE
                  A12 = A12 * BASE
               ELSE
                  MX2 = MX2 * BASE
                  T11 = T11 * BASE
                  T22 = T22 * BASE
                  T12 = T12 * BASE
               END IF
               GOTO 50
            END IF
         END IF
         T12 = T11 * A12 + T12 * A22
         T11 = T11 * A11
         T22 = T22 * A22
         IF ( I.LT.K ) THEN
            DWORK(PW)   = T11
            DWORK(PW+1) = T12
            DWORK(PW+2) = T22
            PW = PW + 3
         END IF
   60 CONTINUE
C
C     Compute the SVD of this product avoiding unnecessary
C     overflow/underflow in the singular values.
C
      TEMP = MAX( ABS( T11 / TWO ) + ABS( T12 / TWO ),
     $            ABS( T22 / TWO ) )
      IF ( TEMP.GT.( RMAX/( TWO * TWOS ) ) ) THEN
         TEMP = TEMP / BASE
         T11  =  T11 / BASE
         T12  =  T12 / BASE
         T22  =  T22 / BASE
      END IF
   70 CONTINUE
      IF ( TEMP.LT.( RMAX/( TWO * BASE * TWOS ) ) .AND.
     $     T11.NE.ZERO .AND. T22.NE.ZERO ) THEN
         SCL = 0
         IF ( ABS( T22 ).LE.TWOS * RMIN ) THEN
            SCL = 1
         ELSE IF ( EPS * ABS( T12 ).GT.ABS( T22 ) ) THEN
            IF ( SQRT( ABS( T11 ) ) * SQRT( ABS( T22 ) ).LE.
     $           ( SQRT( TWOS ) * RMNS ) * SQRT( ABS( T12 ) ) )
     $          SCL = 1
         ELSE
            IF ( ABS( T11 ).LE.TWOS * RMIN *
     $            ( ONE + ABS( T12 / T22 ) ) )
     $         SCL = 1
         END IF
         IF ( SCL.EQ.1 ) THEN
            TEMP = TEMP * BASE
            T11  =  T11 * BASE
            T12  =  T12 * BASE
            T22  =  T22 * BASE
            GOTO 70
         END IF
      END IF
C
      CALL DLASV2( T11, T12, T22, SSMIN, SSMAX, SR, CR, SL, CL )
C
C     Now, the last transformation is propagated to the front as
C     described in [1].
C
      S11 = T11
      S22 = T22
C
      CV(K) = CR
      SV(K) = SR
C
      DO 80  I = K, 2, -1
         AI = AMAP(I)
         IF ( S(AI).EQ.SINV ) THEN
            A11 = A(1,1,AI)
            A12 = A(1,2,AI)
            A22 = A(2,2,AI)
         ELSE
            A11 =  A(2,2,AI)
            A12 = -A(1,2,AI)
            A22 =  A(1,1,AI)
         END IF
         IF ( I.GT.2 ) THEN
            PW = PW - 3
            T11 = DWORK(PW)
            T12 = DWORK(PW+1)
            T22 = DWORK(PW+2)
            IF ( ABS( SR * CL * S22 ).LT.ABS( SL * CR * S11 ) ) THEN
               B11 =  T22
               B22 =  T11
               B12 = -T12
               CC  =  CL
               SC  =  SL
            ELSE
               B11 = A11
               B12 = A12
               B22 = A22
               CC  = CR
               SC  = SR
            END IF
            MX = MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ) )
            IF ( MX.GT.RMAX / TWO ) THEN
               B11 = B11 / TWO
               B22 = B22 / TWO
               B12 = B12 / TWO
            END IF
            CALL DLARTG( B11 * CC + B12 * SC, SC * B22, CC, SC, TEMP )
         ELSE
            CC = CL
            SC = SL
         END IF
         IF ( ABS( SC ).LT.SFMN * ABS( A22 ) ) THEN
            A(1,1,AI) = SC * SR * A22 + CC * ( CR * A11 + SR * A12 )
         ELSE
            A(1,1,AI) = ( A22 / SC ) * SR
         END IF
         IF ( ABS( SR ).LT.SFMN * ABS( A11 ) ) THEN
            A(2,2,AI) = SC * SR * A11 + CR * ( CC * A22 - SC * A12 )
         ELSE
            A(2,2,AI) = ( A11 / SR ) * SC
         END IF
         A(1,2,AI) = ( A12 * CR - A11 * SR ) * CC + A22 * CR * SC
         IF ( S(AI).NE.SINV ) THEN
            TEMP      =  A(1,1,AI)
            A(1,1,AI) =  A(2,2,AI)
            A(2,2,AI) =  TEMP
            A(1,2,AI) = -A(1,2,AI)
         END IF
         CR = CC
         SR = SC
         CV(I-1) = CR
         SV(I-1) = SR
         S11 = T11
         S22 = T22
   80 CONTINUE
C
      CV(1) = CL
      SV(1) = SL
C
      RETURN
C *** Last line of MB03BC ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


!***********************************/
!	Name:計算領域の設定を制御するためのプログラム
!	Alias:SetRegion
!	Description:
!	Type:Geom
!	Input:
!	Output:Geom
!	Note:1~3次元すべてに対応済み
!	Author:Akitaka Toyota
!	Date:2017.10.26
!	Update:2017.11.09
!	Other:
!***********************************/
subroutine SetRegion(Conf,Geom)
    use StructVar_Mod
    use LoopVar_Mod
    implicit none

    type(Configulation),intent(in) :: Conf
    type(Geometry),intent(inout) :: Geom

    if(Conf%UseReadRegion == 1) then
        call ReadRegionData(Geom)
    else if(Conf%UseReadRegion == 0) then
        call SetRegionData(Geom)
    end if

    allocate(Geom%Area(Geom%Dimension))
    allocate(Geom%Width(1,1,Geom%Dimension))
    allocate(Geom%Volume(1))

        iXmax = Geom%CellNumber(1)
        iYmax = max(1,Geom%CellNumber(2)) !セル総数と1の大きい方がiYmax
        iZmax = max(1,Geom%CellNumber(3))
    if(Conf%UseOverSet == 1) then
        !allocate(Geom%CellType(iXmax,iYmax,iZmax))
        allocate(Geom%Interpolated(iXmax,iYmax,iZmax))
    end if

!    allocate(Geom%Vector(Geom%Dimension,Geom%Dimension))

    call CalcWidth(Geom)
    call CalcVolume(Geom)
    call CalcEdgeArea(Geom)
    !call CalcVector(Geom)

return
end subroutine SetRegion