text:
C *********************************************************************
C Programme de calcul de matrice de scattering
C Couplage Grill Plasma avec differents cas poloidaux
C *********************************************************************

      Common/err/epsabs
      Common/nlim/nlimit
      Common/dim/b,z
      Common/idem/bcte
      Common/onde/k0
      Common/plasma/X0,D0,pertes,S_dielectric,P_dielectric
      Common/integration/max_nz

      Integer *4 nbgp,i,j,q,intg,ni,nlimit
      Integer *4 T_grill,D_guide_max
      Integer *4 Pg,Nmh,Nme,Nmhgn,Nmegn,Nmheg,Nmhegn,prhe,prhegn
      Integer *4 Pgg(5),Nmhg(5),Nmeg(5),prheg(5),Nmhg1(5)
      Integer *4 prh(60),pre(60),prhg(5,60),preg(5,60)

      Real *8 pi,k0,f,cl,me,Eps0,qe,epsrel,epsabs
      Real *8 a,ac,b(60),bc,b0,ei,e(60),z(60),pc
      Real *8 X0,D0,pertes,max_nz,n0,nc,dn0,B0_version7
      Real *8 k0g(5),D0g(5),X0g(5)
      Real *8 ag(5),bg(5,60),eg(5,60),zg(5,60)
      Real *8 SMAmpli(2*10/2*60*5,2*10/2*60*5)
      Real *8 SMPhas(2*10/2*60*5,2*10/2*60*5)
      Real *8 SMAmplig(2*10*60*5,2*10*60*5)
      Real *8 ne,nD,nT,nH,nHe,wce,wcD,wcT,wcH,wcHe,wpe,wpD,wpT,wpH,wpHe
      Real *8 mD,mT,mH,mHe
      Real *8 S_dielectric,P_dielectric

      Complex *16 S(2*10*60,2*10*60),Spr(2*10/2*60,2*10/2*60)
      Complex *16 Sg(2*10*60*5,2*10*60*5),Sprg(10*60*5,10*60*5)
      Complex *16 RZheg(2*10*60,2*10*60),RZhegpr(2*10/2*60,2*10/2*60)
      Complex *16 RZhegg(2*10*60*5,2*10*60*5)
      Complex *16 RZhegprg(10*60*5,10*60*5)

      Character *1 repspc,repgril,reppr,sortie,bcte,rep

C Declaration des constantes

      pi=4.*ATAN(1.)
      cl=3.E8
      me=9.1091E-31
      Eps0=8.854E-12
      qe=1.6021E-19
      mD=3.345E-27
      mT=5.018E-27
      mH=1.673E-27
      mHe=5.018E-27

      pc=1.*1.E-3

      Open(10,File='par_grill.dat',Status='unknown',Form='formatted')

      Read(10,*) Nmh
      Read(10,*) Nme
      Read(10,*) f
      Read(10,*) n0
      Read(10,*) dn0
      Read(10,*) B0_version7
      Read(10,*) Pg
      Read(10,*) a
      Read(10,*) (b(i),i=1,Pg)
      Read(10,*) (z(i),i=1,Pg)
      Read(10,*) T_grill
      Read(10,*) D_guide_max
      Read(10,*) epsrel
      Read(10,*) pertes
      Read(10,*) max_nz

      Rewind 10
      Close(10)

      k0=2.*pi*f/cl
      nc=(2.*pi*f)**2*me*Eps0/qe**2

      X0=n0/nc
      D0=k0*n0/dn0

      ne=n0
      nD=0.9*ne
      nT=0.
      nH=0.1*ne
      nHe=0.

      wce=qe*B0_version7/me
      wcD=qe*B0_version7/mD
      wcT=qe*B0_version7/mT
      wcH=qe*B0_version7/mH
      wcHe=2*qe*B0_version7/mHe

      wpe=sqrt(ne*qe**2/(Eps0*me))
      wpD=sqrt(nD*qe**2/(Eps0*mD))
      wpT=sqrt(nT*qe**2/(Eps0*mT))
      wpH=sqrt(nH*qe**2/(Eps0*mH))
      wpHe=sqrt(nHe*(2*qe)**2/(Eps0*mHe))

      S_dielectric=1-( wpe**2/((2*pi*f)**2-wce**2) +
     &      wpD**2/((2*pi*f)**2-wcD**2) +
     &      wpT**2/((2*pi*f)**2-wcT**2) +
     &      wpH**2/((2*pi*f)**2-wcH**2) +
     &      wpHe**2/((2*pi*f)**2-wcHe**2) )

      P_dielectric=1-( (wpe/(2*pi*f))**2 +
     &      (wpD/(2*pi*f))**2 +
     &      (wpT/(2*pi*f))**2 +
     &      (wpH/(2*pi*f))**2 +
     &      (wpHe/(2*pi*f))**2 )

      S_dielectric=1.

      nlimit=100000
      ni = 32
      epsabs=0.

      Call SMGrill(intg,ni,epsrel,T_grill,D_guide_max,Nmh,Nme,Pg,a,
     &prh,pre,S,RZheg,Spr,RZhegpr)

C Sauvegarde des resultats pour matlab

      Open(10,File='S_plasma.dat',Status='unknown',Form='unformatted')

      Write(10) ((S(i,j),i=1,(Nme+Nmh)*Pg),j=1,(Nme+Nmh)*Pg)
      Write(10) ((RZheg(i,j),i=1,(Nme+Nmh)*Pg),j=1,(Nme+Nmh)*Pg)

      Rewind 10
      Close(10)

      Stop

      End

C Fin du programme principal

C ***
C Surotine de CONTrole de la SYMetrie du module d une mat. de scattering
C ***
      Subroutine ContSymg(Dim,Sg,Statsymgg)

      Integer *4 i,j,Dim
      Integer *4 Statsymgg(2*10*60*5,2*10*60*5)
      Real *8 A21r,A12r,Diffr
      Real *8 A21i,A12i,Diffi
      Complex *16 Sg(2*10*60*5,2*10*60*5)
      Character *1 rep

      Do 100 j=1,Dim
         Do 110 i=1,Dim
            If (i.NE.j) Then
               A12r=DBLE(Sg(i,j))
               A21r=DBLE(Sg(j,i))
               A12i=DIMAG(Sg(i,j))
               A21i=DIMAG(Sg(j,i))
               Diffr=DABS(A12r-A21r)
               Diffi=DABS(A12i-A21i)
               If ((Diffr.LT.0.001).AND.(Diffi.LT.0.001)) Then
                  Statsymgg(i,j)=0
               Else
C               If ((Diffr.GT.0.001).OR.(Diffi.GT.0.001)) Then
                  Statsymgg(i,j)=-1
C                  Stop 'Probleme de Symetrie de la matice '
               EndIf
            EndIf
 110     Continue
 100  Continue

   1  Return
      End

C ***
C Suroutine de CONTrole de l UNItarite d une matrice de scattering
C ***
      Subroutine ContUnig(Dim,Sprg,Statuni)

      Integer *4 i,j,Dim
      Integer *4 Statuni(10*60*5)
      Real *8 Diff,Somm,Som(10*60*5)
      Complex *16 Sprg(10*60*5,10*60*5)
      Character *1 rep

      Do 100 j=1,Dim
         Somm=0.
         Do 110 i=1,Dim
            Somm=Somm+(CDABS(Sprg(i,j)))**2
 110     Continue
         Som(j)=Somm
C         Diff=DABS(Somm-1)
C         If (Diff.LT.0.01) Then
         If (Somm.LT.1.) Then
            Statuni(i)=0
         Else
             Statuni(i)=-1
C             Write (*,*) 'Probleme d Unitarite '
         EndIf
 100  Continue

C      Write (*,1) (Som(j),j=1,Dim)
C   1  Format(F8.4)
      If (Dim.EQ.1) Then
         Write (*,1) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.2) Then
         Write (*,2) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.3) Then
         Write (*,3) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.4) Then
         Write (*,4) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.5) Then
         Write (*,5) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.6) Then
         Write (*,6) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.7) Then
         Write (*,7) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.8) Then
         Write (*,8) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.9) Then
         Write (*,9) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.10) Then
         Write (*,10) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.11) Then
         Write (*,11) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.12) Then
         Write (*,12) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.13) Then
         Write (*,13) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.14) Then
         Write (*,14) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.15) Then
         Write (*,15) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.16) Then
         Write (*,16) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.17) Then
         Write (*,17) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.18) Then
         Write (*,18) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.19) Then
         Write (*,19) (Som(j),j=1,Dim)
      ElseIf (Dim.EQ.20) Then
         Write (*,20) (Som(j),j=1,Dim)
      EndIf

   1  Format(1(F8.4))
   2  Format(2(F8.4))
   3  Format(3(F8.4))
   4  Format(4(F8.4))
   5  Format(5(F8.4))
   6  Format(6(F8.4))
   7  Format(7(F8.4))
   8  Format(8(F8.4))
   9  Format(9(F8.4))
  10  Format(10(F8.4))
  11  Format(11(F8.4))
  12  Format(12(F8.4))
  13  Format(13(F8.4))
  14  Format(14(F8.4))
  15  Format(15(F8.4))
  16  Format(16(F8.4))
  17  Format(17(F8.4))
  18  Format(18(F8.4))
  19  Format(19(F8.4))
  20  Format(20(F8.4))

      Return
      End

C ***
C Subroutine de CALcul de la MATrice d AMPLItude d une matrice complexe
C ***
      Subroutine CalMatAmplig(prhe,Sprg,Matamp)

      Integer *4 i,j,prhe
      Real *8 Matamp(2*10/2*60*5,2*10/2*60*5)
      Complex *16 Sprg(2*10/2*60*5,2*10/2*60*5)

      Do 100 j=1,prhe
         Do 110 i=1,prhe
            Matamp(i,j)=CDABS(Sprg(i,j))
 110     Continue
 100  Continue

      Return
      End
C ***
      Subroutine CalMatAmpligg(prhe,Sg,Matampg)

      Integer *4 i,j,prhe
      Real *8 Matampg(2*10*60*5,2*10*60*5)
      Complex *16 Sg(2*10*60*5,2*10*60*5)

      Do 100 j=1,prhe
         Do 110 i=1,prhe
            Matampg(i,j)=CDABS(Sg(i,j))
 110     Continue
 100  Continue

      Return
      End

C ***
C Subroutine de CALcul de la MATrice de PHASe d une matrice complexe

      Subroutine CalMatPhasg(prhe,Sprg,Matpha)

      Integer *4 i,j,prhe
      Real *8 pReel,pIma,phi,pi
      Real *8 Matpha(2*10/2*60*5,2*10/2*60*5)
      Complex *16 Sprg(2*10/2*60*5,2*10/2*60*5)

      pi=4.*ATAN(1.)

      Do 100 j=1,prhe
         Do 110 i=1,prhe
            pReel=DBLE(Sprg(i,j))
            pIma=DIMAG(Sprg(i,j))
            If (pReel.EQ.0) Then
               If (pIma.GT.0) Then
                  Matpha(i,j)=90.
               ElseIf (pIma.LT.0) Then
                  Matpha(i,j)=-90.
               Else
                  Matpha(i,j)=0.
               EndIf
            Else
               phi=pIma/pReel
               If (pReel.GT.0) Then
                  Matpha(i,j)=(180./pi)*DATAN(phi)
               Else
                  Matpha(i,j)=180.+(180./pi)*DATAN(phi)
                  If (Matpha(i,j).GT.180.) Then
                     Matpha(i,j)=Matpha(i,j)-360.
                  EndIf
               EndIf
            EndIf
 110     Continue
 100  Continue

      Return
      End

C ***
C Subroutine pour IMPrimer une MATrice Reelle d Amplitude
C ***
      Subroutine ImpMatRAg(lprhe,cprhe,mat)

      Integer *4 i,j,cprhe,lprhe
      Real *8 mat(2*10/2*60*5,2*10/2*60*5)

      If (cprhe.EQ.1) Then
         Write (*,1) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.2) Then
         Write (*,2) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.3) Then
         Write (*,3) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.4) Then
         Write (*,4) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.5) Then
         Write (*,5) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.6) Then
         Write (*,6) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.7) Then
         Write (*,7) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.8) Then
         Write (*,8) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.9) Then
         Write (*,9) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.10) Then
         Write (*,10) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.11) Then
         Write (*,11) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.12) Then
         Write (*,12) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.13) Then
         Write (*,13) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.14) Then
         Write (*,14) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.15) Then
         Write (*,15) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.16) Then
         Write (*,16) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.17) Then
         Write (*,17) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.18) Then
         Write (*,18) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.19) Then
         Write (*,19) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.20) Then
         Write (*,20) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.21) Then
         Write (*,21) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.22) Then
         Write (*,22) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.23) Then
         Write (*,23) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.24) Then
         Write (*,24) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.25) Then
         Write (*,25) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.26) Then
         Write (*,26) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.27) Then
         Write (*,27) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.28) Then
         Write (*,28) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.29) Then
         Write (*,29) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.30) Then
         Write (*,30) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      EndIf

   1  Format(1(F8.4))
   2  Format(2(F8.4))
   3  Format(3(F8.4))
   4  Format(4(F8.4))
   5  Format(5(F8.4))
   6  Format(6(F8.4))
   7  Format(7(F8.4))
   8  Format(8(F8.4))
   9  Format(9(F8.4))
  10  Format(10(F8.4))
  11  Format(11(F8.4))
  12  Format(12(F8.4))
  13  Format(13(F8.4))
  14  Format(14(F8.4))
  15  Format(15(F8.4))
  16  Format(16(F8.4))
  17  Format(17(F8.4))
  18  Format(18(F8.4))
  19  Format(19(F8.4))
  20  Format(20(F8.4))
  21  Format(21(F8.4))
  22  Format(22(F8.4))
  23  Format(23(F8.4))
  24  Format(24(F8.4))
  25  Format(25(F8.4))
  26  Format(26(F8.4))
  27  Format(27(F8.4))
  28  Format(28(F8.4))
  29  Format(29(F8.4))
  30  Format(30(F8.4))

      Return
      End
C ***
      Subroutine ImpMatRAgg(lprhe,cprhe,matg)

      Integer *4 i,j,cprhe,lprhe
      Real *8 matg(2*10*60*5,2*10*60*5)

      If (cprhe.EQ.1) Then
         Write (*,1) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.2) Then
         Write (*,2) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.3) Then
         Write (*,3) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.4) Then
         Write (*,4) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.5) Then
         Write (*,5) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.6) Then
         Write (*,6) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.7) Then
         Write (*,7) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.8) Then
         Write (*,8) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.9) Then
         Write (*,9) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.10) Then
         Write (*,10) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.11) Then
         Write (*,11) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.12) Then
         Write (*,12) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.13) Then
         Write (*,13) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.14) Then
         Write (*,14) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.15) Then
         Write (*,15) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.16) Then
         Write (*,16) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.17) Then
         Write (*,17) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.18) Then
         Write (*,18) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.19) Then
         Write (*,19) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.20) Then
         Write (*,20) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.21) Then
         Write (*,21) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.22) Then
         Write (*,22) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.23) Then
         Write (*,23) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.24) Then
         Write (*,24) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.25) Then
         Write (*,25) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.26) Then
         Write (*,26) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.27) Then
         Write (*,27) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.28) Then
         Write (*,28) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.29) Then
         Write (*,29) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.30) Then
         Write (*,30) ((matg(i,j),j=1,cprhe),i=1,lprhe)
      EndIf

   1  Format(1(F8.4))
   2  Format(2(F8.4))
   3  Format(3(F8.4))
   4  Format(4(F8.4))
   5  Format(5(F8.4))
   6  Format(6(F8.4))
   7  Format(7(F8.4))
   8  Format(8(F8.4))
   9  Format(9(F8.4))
  10  Format(10(F8.4))
  11  Format(11(F8.4))
  12  Format(12(F8.4))
  13  Format(13(F8.4))
  14  Format(14(F8.4))
  15  Format(15(F8.4))
  16  Format(16(F8.4))
  17  Format(17(F8.4))
  18  Format(18(F8.4))
  19  Format(19(F8.4))
  20  Format(20(F8.4))
  21  Format(21(F8.4))
  22  Format(22(F8.4))
  23  Format(23(F8.4))
  24  Format(24(F8.4))
  25  Format(25(F8.4))
  26  Format(26(F8.4))
  27  Format(27(F8.4))
  28  Format(28(F8.4))
  29  Format(29(F8.4))
  30  Format(30(F8.4))

      Return
      End

C Subroutine pour IMPrimer une MATrice Reelle de Phase
C ***
      Subroutine ImpMatRPg(lprhe,cprhe,mat)

      Integer *4 i,j,cprhe,lprhe
      Real *8 mat(2*10/2*60*5,2*10/2*60*5)

      If (cprhe.EQ.1) Then
         Write (*,1) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.2) Then
         Write (*,2) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.3) Then
         Write (*,3) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.4) Then
         Write (*,4) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.5) Then
         Write (*,5) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.6) Then
         Write (*,6) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.7) Then
         Write (*,7) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.8) Then
         Write (*,8) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.9) Then
         Write (*,9) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.10) Then
         Write (*,10) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.11) Then
         Write (*,11) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.12) Then
         Write (*,12) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.13) Then
         Write (*,13) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.14) Then
         Write (*,14) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.15) Then
         Write (*,15) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.16) Then
         Write (*,16) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.17) Then
         Write (*,17) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.18) Then
         Write (*,18) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.19) Then
         Write (*,19) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.20) Then
         Write (*,20) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.21) Then
         Write (*,21) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.22) Then
         Write (*,22) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.23) Then
         Write (*,23) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.24) Then
         Write (*,24) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.25) Then
         Write (*,25) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.26) Then
         Write (*,26) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.27) Then
         Write (*,27) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.28) Then
         Write (*,28) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.29) Then
         Write (*,29) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      ElseIf (cprhe.EQ.30) Then
         Write (*,30) ((mat(i,j),j=1,cprhe),i=1,lprhe)
      EndIf

   1  Format(1(F8.2))
   2  Format(2(F8.2))
   3  Format(3(F8.2))
   4  Format(4(F8.2))
   5  Format(5(F8.2))
   6  Format(6(F8.2))
   7  Format(7(F8.2))
   8  Format(8(F8.2))
   9  Format(9(F8.2))
  10  Format(10(F8.2))
  11  Format(11(F8.2))
  12  Format(12(F8.2))
  13  Format(13(F8.2))
  14  Format(14(F8.2))
  15  Format(15(F8.2))
  16  Format(16(F8.2))
  17  Format(17(F8.2))
  18  Format(18(F8.2))
  19  Format(19(F8.2))
  20  Format(20(F8.2))
  21  Format(21(F8.2))
  22  Format(22(F8.2))
  23  Format(23(F8.2))
  24  Format(24(F8.2))
  25  Format(25(F8.2))
  26  Format(26(F8.2))
  27  Format(27(F8.2))
  28  Format(28(F8.2))
  29  Format(29(F8.2))
  30  Format(30(F8.2))

      Return
      End

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


text:
C    @(#)bldxti.f	20.11 10/13/99
      subroutine bldxti (z_threshold)
      real z_threshold

C     Build "TXTIE" array of phase shifter bus ties and
C     optionally, of low impedance branches

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/beta.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/brtype.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/komps.inc'
      include 'ipfinc/optim1.inc'
      include 'ipfinc/phase.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/qksrt.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/snput.inc'
      include 'ipfinc/tran.inc'
      include 'ipfinc/lfiles.inc'

      common /is_batch / is_batch

      common /ran002/ brtypec(9)
      character brtypec*2

      complex * 16 y(2,2)
      double precision ratio, angle
      character id*1
      integer p, q, pold, error, ktxtie, intrst(MAXBUS)
      logical found, finished
      external kmptxt,swptxt

      ntxtie = 0
      if (option(9) .le. 0) go to 900

      do i = 1, ntot
         intrst(i) = 0
      enddo
c
c     Pass 1 - process phase shifters
c
      do 110 jt = 1,jphno
        k1 = jphid(1,jt)
        k2 = jphid(2,jt)
        kt = inp2opt(k1)
        mt = inp2opt(k2)
        ltc = 0
C
C       Find LTC (automatic transformer) and NBR (branch) index.
C
        do ksw = 0, 1
          p = kbsdta(16,k1)
          do while (p .gt. 0)
            if (ky(p) .eq. k2) then
               found = .true.
               do while (p .gt. 0)  ! Found branch k1-k2. Restrict loop
C                                   ! to this branch only.
                 ltyp = brtype(p)
C
                 if (ltyp .eq. BRTYP_R) then
C
C                  LTC phase shifters will be processed in a
C                  special LP mode. Store index in TXTIE(7,*)
C
                   do i = 1, ntota
                     if (ltran(10,i) .eq. 3) then
                       if (ltran(1,i) .eq. kt .and.
     &                     ltran(9,i) .eq. mt) then
                         ltc = i
                         go to 100
                       else if (ltran(1,i) .eq. mt .and.
     &                          ltran(9,i) .eq. kt) then
                         ltc = -i
                         go to 100
                       endif
                     endif
                   enddo
                   ltc = 0
  100              continue

                 else if (ltyp .eq. BRTYP_TP) then
                   q = brnch_ptr(p)
                   nbr = iabs(q)
                   if (brnch(5,nbr) .gt. brnch(6,nbr) .or.
     &                 brnch(5,nbr) .lt. 0.0 .or.
     &                 brnch(6,nbr) .lt. 0.0) then
                     write (errbuf(1), 10000) bus(k1), base(k1),
     &                 bus(k2),  base(k2), id, brtypec(ltyp),
     &                 brnch(5,nbr), brnch(6,nbr)
10000                format ('0 Illegal phase shifter impedance ',
     &                 a8, f7.1, 1x, a8, f7.1, ' id ', a, ' type ',
     &                 a2, ' R, X = ', 2f8.5)
                     if (is_batch .eq. 0) then
                       call prterx ('W',1)
                     else
                       call prterx ('E',1)
                     endif
                   endif
                   ratio = 1.0d0
                   if (abs(brnch(9,nbr)) .lt. option(9) .and.
     &                 ltc .eq. 0) go to 110
                   if (q .gt. 0) then
                      angle = 0.0174532925d0 * brnch(9,nbr)
                   else
                      angle = -0.0174532925d0 * brnch(9,nbr)
                   endif
c
c                  Consolidate parallels
c
                   ktxtie = 0
                   i = 1
                   do while (i .le. ntxtie .and. ktxtie .eq. 0)
                     if (txtie(1,i) .eq. kt .and.
     &                   txtie(2,i) .eq. mt) then
                       ktxtie = i
                     else
                       i = i + 1
                     endif
                   enddo
                   if (ktxtie .eq. 0) then
                     ntxtie = ntxtie + 1
                     if (ntxtie .gt. MAXBUSTIE) then
                       write (errbuf(1), 10010) MAXBUSTIE
10010                  format(' More than ', i4, ' bus ties.  ',
     &                        'Ideal bus tie scheme aborted.')
                       call prterx ('W',1)
                       ntxtie = 0
                       go to 900
                     endif
                     ktxtie = ntxtie
                     txtie(1,ntxtie) = kt
                     txtie(2,ntxtie) = mt
                     txtie(3,ntxtie) = ratio
                     if (ltc .eq. 0) then
                       txtie(4,ntxtie) = angle
                     else
                       x = isign (1,ltc)
                       i = iabs (ltc)
                       txtie(4,ntxtie) = x * tap(i)
                     endif
                     intrst(k1) = ntxtie
                     intrst(k2) = ntxtie
C
C                    The following is a temporary flag to indicate
C                    an LTC phase shifter.  It will be sorted last in
C                    TXTIE array.
C
                     txtie(7,ntxtie) = ltc
                     txtie(8,ntxtie) = 0

                     txtie(5,ntxtie) = 0.0
                     txtie(6,ntxtie) = 0.0
                     txtie(9,ntxtie) = 0.0
                     txtie(10,ntxtie) = 0.0
                   endif
C
C                  Obtain 2-port Y-matrix
C
                   call pieqiv (p,y,error)
                   txtie(5,ktxtie) = txtie(5,ktxtie) + dreal(y(1,2))
                   txtie(6,ktxtie) = txtie(6,ktxtie) + dimag(y(1,2))
                   txtie(9,ktxtie) = txtie(9,ktxtie) + dreal(y(1,1))
                   txtie(10,ktxtie) = txtie(10,ktxtie) + dimag(y(1,1))

                 else if ( ltyp .eq. BRTYP_PEQ ) then
c
c                  ignore phase shifter if in section
c
                   p = 0        ! force end of loop

                 else
C
C                  Illegal branch type - program error!!
C
                   id = brid(p)
                   write (errbuf(1), 10020) bus(k1), base(k1),
     &                bus(k2),  base(k2), id, brtypec(ltyp)
10020              format ('0 Illegal branch type ', a8, f7.1, 1x,
     &                a8,f7.1, ' id ', a, ' type ', a2,
     &                ' in parallel with a phase shifter ')
                   if (is_batch .eq. 0) then
                     call prterx ('W',1)
                   else
                     call prterx ('E',1)
                   endif
                 endif

                 if (p.gt.0) p = brnch_nxt(p)
                 if (p .gt. 0) then
                   if (ky(p) .ne. k2) p = 0  ! End of branch k1-k2
C                                            ! Force end of loop
                 endif
               enddo
            else
               p = brnch_nxt(p)
            endif
          enddo
          if (.not.found) then
C
C           Could not locate branch - program error!!
C
            write (errbuf(1), 10030) bus(k1), base(k1), bus(k2),
     &         base(k2)
10030       format ('0 Could not locate phase shifter ', a8, f7.1,
     &         1x, a8, f7.1)
            if (is_batch .eq. 0) then
               call prterx ('W',1)
            else
               call prterx ('E',1)
            endif
          endif
c
c         Process transpose k2-k1 using same logic.
c
          ktemp = k1
          k1 = k2
          k2 = ktemp
          kt = inp2opt(k1)
          mt = inp2opt(k2)
          ltc = -ltc
        enddo
  110 continue
c
c     Pass 2 - process low impedance branches
c
      if (z_threshold .le. 0.0) go to 900
      y_threshold = 1.0 / z_threshold
      do 140 ib = 1, ntot_alf
        nb = alf2inp(ib)
        p = kbsdta(16,nb)
        do while (p .gt. 0)
c
c         Process new branch, Skip type R, LD, LM, and TP
c         records (TP's have been processed in pass 1 above)
c
          k1 = kx(p)
          k2 = ky(p)
          kt = inp2opt(k1)
          mt = inp2opt(k2)
          if (brtype(p) .eq. 4 .or. brtype(p) .eq. 2 .or.
     &        brtype(p) .eq. 6 .or. brtype(p) .eq. 7) then
            p = brnch_nxt(p)
            go to 130
c
c         For this version only, supplement TP records with T record
c         bus ties.
c
          else if (brtype(p) .eq. 1) then
            do while (p .gt. 0 .and. ky(p) .eq. k2)
              p = brnch_nxt(p)
            enddo
            go to 130
          else if (brtype(p) .ne. 5) then
            p = brnch_nxt(p)
            go to 130
          else
            if (intrst(k1) .ne. 0 .and. intrst(k2) .ne. 0) then
              found = .false.
              jt = 1
              do while (jt .le. ntxtie .and. .not. found)
                if (txtie(1,jt) .eq. kt .and. txtie(2,jt) .eq. mt) then
                  found = .true.
                else
                  jt = jt + 1
                endif
              enddo
              if (found) then
                do while (p .gt. 0 .and. ky(p) .eq. k2)
                  p = brnch_nxt(p)
                enddo
                go to 130
              endif
            endif
          endif
c
c         Use Y-matrix to screen admittance
c
          do l = km(kt), km(kt)+kmlen(kt)-1
            if (ikmu(l) .eq. mt) then
              if (abs (bkmu(l)) .gt. y_threshold) then
                go to 120
              else
                p = brnch_nxt(p)
                go to 130
              endif
            endif
          enddo
          p = brnch_nxt(p)
          go to 130

  120     continue
          ltyp = brtype(p)
          if (ltyp .eq. 1 .or. ltyp .eq. 3 .or. ltyp .eq. 8) then
            ratio = 1.0d0
            angle = 1.0d0
          else if (ltyp .eq. 5) then
            q = brnch_ptr(p)
            nbr = iabs (q)
            if (q .gt. 0) then
              ratio = (brnch(9,nbr) / base(k1))
     &              / (brnch(10,nbr) / base(k2))
              angle = 0.0
            else
              ratio = (brnch(10,nbr) / base(k1))
     &              / (brnch(9,nbr) / base(k2))
              angle = 0.0
            endif
          else
            write (errbuf(1), 10040) bus(k1), base(k1),
     &        bus(k2),  base(k2), brid(p), brtypec(ltyp)
10040       format ('0 Program error 01 (bldxti) ', a8, f7.1, 1x,
     &        a8,f7.1, ' id ', a, ' type ', a2)
            if (is_batch .eq. 0) then
              call prterx ('W',1)
            else
              call prterx ('E',1)
            endif
          endif
C
C         Obtain 2-port Y-matrix
C
          call pieqiv (p, y, error)
          if (dimag(y(1,2)) .lt. y_threshold) then
            p = brnch_nxt(p)
            go to 130
          endif
          ntxtie = ntxtie + 1
          if (ntxtie .gt. MAXBUSTIE) then
            write (errbuf(1), 10010) MAXBUSTIE
            call prterx ('W',1)
            ntxtie = 0
            go to 900
          endif
          txtie(1,ntxtie) = kt
          txtie(2,ntxtie) = mt
          txtie(3,ntxtie) = ratio
          txtie(4,ntxtie) = angle
          intrst(k1) = ntxtie
          intrst(k2) = ntxtie
          txtie(5,ntxtie) = dreal(y(1,2))
          txtie(6,ntxtie) = dimag(y(1,2))
          txtie(7,ntxtie) = 0.0
          txtie(8,ntxtie) = 0
          txtie(9,ntxtie) = dreal(y(1,1))
          txtie(10,ntxtie) = dimag(y(1,1))
c
c         Loop through any parallel branches
c
          finished = .false.
          do while (.not. finished)
            if (ltyp .eq. 1) then
              pold = p
              p = brnch_nxt(p)
              do while (p .gt. 0 .and. ky(p) .eq. k2 .and.
     &                  brid(p) .eq. brid(pold))
                p = brnch_nxt(p)
              enddo
            else
              p = brnch_nxt(p)
            endif
            if (ky(p) .eq. k2) then
C
C             Obtain 2-port Y-matrix
C
              call pieqiv (p, y, error)
              txtie(5,ntxtie) = txtie(5,ntxtie) + dreal(y(1,2))
              txtie(6,ntxtie) = txtie(6,ntxtie) + dimag(y(1,2))
              txtie(9,ntxtie) = txtie(9,ntxtie) + dreal(y(1,1))
              txtie(10,ntxtie) = txtie(10,ntxtie) + dimag(y(1,1))

            else
              finished = .true.
            endif
          enddo
  130     continue
       enddo
  140  continue

       if (ntxtie .gt. 0) then
         call qiksrt(1, ntxtie, kmptxt, swptxt)
         if ( idswa .ne. 0 ) then
           write (dbug, 10050)
10050      format( '0 DUMP OF txtie(i,j) from BLDXTI' )
            do i = 1,ntxtie
              write (dbug, 10060) i, ifix(sngl(txtie(1,i))),
     &              ifix(sngl(txtie(2,i))), (txtie(j,i), j=3,6),
     &              ifix(sngl(txtie(7,i))), ifix(sngl(txtie(8,i))),
     &              (txtie(j,i), j = 9,10)
10060         format( i3, 2i5, 2f7.3, 1x, 2e10.3, 2i5, 2e10.3)
            enddo
          endif
        endif

  900   continue

        return
        end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE INITER_DEM
C***********************************************************************
C*               INITER_DEM Module of AERMAP Terrain Preprocessor
C*
C*       PURPOSE: Read Terrain Data from Raw DEM Data Files and Write Out
C*                the Data to Direct Access Binary Terrain Data Files.
C*                Also Create Record Index Files
C*
C*
C*       PROGRAMMER: Jayant Hardikar, Roger Brode
C*
C*       DATE:    September 29, 1995
C*
C*       MODIFIED: February 9, 2009
C*
C*                Modified output format for UTM Zone to I3 to
C*                accommodate Southern Hemisphere applications,
C*                since Southern zones are negative.  Also
C*                adjusted output format for XBASE and YBASE
C*                to accommodate double precision values.
C*                Adjusted XBASE (longitude) for LAT files to
C*                use negative for west longitude.
C*                Modified to account for no-DOMAIN option.
C*                Also modified handling of missing elevation
C*                data and elevation of local datum (LOCEL).
C*                Added information on direct access files to
C*                the DOMDETAIL.OUT debug file.
C*                Roger W. Brode, U.S. EPA, OAQPS, AQMG
C*
C*       MODIFIED: December 7, 2006
C*
C*                Corrected several problems related to NAD conversion
C*                process, procedure for optimizing critical hill height
C*                calculations for neighboring DEM files, and other issues.
C*                See header comments in AERMAP.FOR source file and
C*                AERMAP MCB#1 for more details.
C*                Roger W. Brode, U.S. EPA, OAQPS, AQMG
C*
C*       INPUTS:  Raw Terrain Data (DEM) Files
C*
C*       OUTPUTS: Direct Access Binary Terrain Files
C*
C*       CALLED FROM:   MAIN
C***********************************************************************

C*    Variable Declarations
      USE MAIN1

      IMPLICIT NONE

      SAVE

      CHARACTER*12 MODNAM

      LOGICAL NEWPRO
      DOUBLE PRECISION XNODEI,YNODEI
      DOUBLE PRECISION NADFLG
      DOUBLE PRECISION ARGE, ARGN, XARG, YARG
      DOUBLE PRECISION XBUTM, YBUTM, XNUTM, YNUTM
      DOUBLE PRECISION XBLON, YBLAT
      INTEGER IPROF, J, K, L, M, IERR, INOD

      INTEGER, ALLOCATABLE :: NUMNOD(:)

      INTEGER NSTART, N, JRECRD, RC, ZMINSHFT, ZMAXSHFT
      INTEGER IZONE, ISPHERE

      logical   :: L_FirstProf  ! flg to indicate first profile within domain

      DOUBLE PRECISION :: ELV

      MODNAM = 'INITER_DEM'

      L_FirstProf = .true.

      XBUTM = 0.0D0
      YBUTM = 0.0D0
      XNUTM = 0.0D0
      YNUTM = 0.0D0
      XBLON = 0.0D0
      YBLAT = 0.0D0

C*    Loop Over Terrain Files
      DEMLOOP: DO IDEM = 1, NUMDEM

C*       Assign J=IDEM for READ Statements
         J = IDEM

C*       Initialize Some Variables
         NUMPRF(IDEM) = 0

C*       Open the DEM File * Read Total Number of Elevation
C*       Profiles etc.
         REWIND (IDMUNT(IDEM))

         IF (FT(J) .NE. 1) THEN
C*          Open as standard formatted file
            OPEN(UNIT = IDMUNT(IDEM), FILE= DEMFIL(IDEM), ERR = 95)
            READ (IDMUNT(IDEM),51, ERR=97)
     &       MAPN(J), FREEF(J), FILR1(J), PROCODE(J), FILR2(J),
     &       SECTNL(J), MCTR(J), DEMLVL(J), ELEVPAT(J),
     &       IPLAN(J), IZO(J), (MPROJ(J,L), L = 1,15), CUNIT(J),
     &       ELUNIT(J), SIDZ(J), ((DMCNR(J,L,M),L=1,2),M=1,4),
     &       ELEVMN(J), ELEVMX(J), CNTRC(J),ACCUC(J),
     &       DXM(J), DYM(J),DCI(J), NROW(J), NPROF(J), LPRIM(J),
     &       LPINT(J), SPRIM(J), SPINT(J),
     &       DDATE(J), DINSP(J), INSPF(J), DVALD(J), SUSF(J),
     &       VDAT(J), NADD(J), EDITN(J), PVOID(J)
         ELSE
C*          For FT(J) = 1 (no record delimiters) open as direct access file
            OPEN(UNIT = IDMUNT(IDEM), FILE= DEMFIL(IDEM),
     &         RECL = 1024,
     &         ACCESS = 'DIRECT',
     &         FORM ='FORMATTED',
     &         ERR = 95)
            READ (IDMUNT(IDEM), FMT=51, REC = 1, ERR=97)
     &       MAPN(J), FREEF(J), FILR1(J), PROCODE(J), FILR2(J),
     &       SECTNL(J), MCTR(J), DEMLVL(J), ELEVPAT(J),
     &       IPLAN(J), IZO(J), (MPROJ(J,L), L = 1,15), CUNIT(J),
     &       ELUNIT(J), SIDZ(J), ((DMCNR(J,L,M),L=1,2),M=1,4),
     &       ELEVMN(J), ELEVMX(J), CNTRC(J), ACCUC(J),
     &       DXM(J),DYM(J),DCI(J), NROW(J), NPROF(J), LPRIM(J),
     &       LPINT(J), SPRIM(J), SPINT(J),
     &       DDATE(J), DINSP(J), INSPF(J), DVALD(J), SUSF(J),
     &       VDAT(J), NADD(J), EDITN(J), PVOID(J)
         END IF

  51     FORMAT(2A40, A55, 2A1, A3, A4, 2I6,
     &          2I6, 15D24.15, I6,
     &          2I6, 4(2D24.15), 2D24.15, D24.15,
     &          I6, 2E12.6, E12.6, 2I6, 2(I5,I1),
     &          2I4, A1, I1, I2,
     &          2I2, 2I4)

C*       Assign DEM zone from Record A to IZOND(NDEM) array
         IF (IZO(IDEM) .NE. 0) THEN
            IZOND(IDEM) = IZO(IDEM)
         END IF

C*       Open Temporary Direct Access Terrain Files and Index Files
         OPEN (IDRUNT(IDEM),FILE=DIRFIL(IDEM),
     &         ACCESS='DIRECT',RECL=LREC_DIR)
         OPEN (IDXUNT(IDEM),FILE=IDXFIL(IDEM))

C*       For Each Profile in the DEM File
         RC = 1

C*       Check for blank (0) NADD from DEM Record Type A
C*       This was already done in CHKADJ based on DEMCHK read,
C*       but needs to be checked again here for new read.
         IF (NADD(IDEM) .EQ. 0 .AND. IPLAN(IDEM) .EQ. 0) THEN
C*          Assume WGS72 for LAT files
            NADD(IDEM) = 2
         ELSE IF (NADD(IDEM) .EQ. 0 .AND. IPLAN(IDEM) .EQ. 1) THEN
C*          Assum NAD27 for UTM files
            NADD(IDEM) = 1
         END IF

         IF (.NOT. GOTDOMFLG) THEN
C*          No user-specified domain; skip domain shift calcs
            XDMNSHFT = 0.0D0
            YDMNSHFT = 0.0D0
            XDMXSHFT = 0.0D0
            YDMXSHFT = 0.0D0
            GO TO 777
         END IF
C*
C*       Determine whether to apply NAD shift to Domain
         IF (NADA .EQ. 0 .OR. NADA .EQ. NADD(IDEM)) THEN

            NADFLG = 0.0D0                ! No NAD shift needed

         ELSE IF( (NADA.EQ.1    .OR.  NADA.GE.5) .AND.
     &            (NADD(IDEM).GE.2 .AND. NADD(IDEM).LE.4) )THEN

            NADFLG = 1.0D0                ! Include NAD shift

         ELSE IF( (NADA.GE.2    .AND. NADA.LE.4) .AND.
     &            (NADD(IDEM).EQ.1 .OR.  NADD(IDEM).GE.5) )THEN

            NADFLG = 1.0D0                ! Include NAD shift

         ELSE

            NADFLG = 0.0D0                ! No NAD shift needed

         END IF

C*       Apply NAD shift if necessary, based on type of DEM data
         IF (TYPDEM .EQ. 'LAT') THEN
C*          Apply datum shift in arc-seconds for LAT files
            XDMNSHFT = DOMLL(2,1) + (XDMNDIFS*NADFLG)
            YDMNSHFT = DOMLL(1,1) + (YDMNDIFS*NADFLG)
            XDMXSHFT = DOMLL(2,3) + (XDMXDIFS*NADFLG)
            YDMXSHFT = DOMLL(1,3) + (YDMXDIFS*NADFLG)
         ELSE IF (TYPDEM .EQ. 'UTM') THEN
C*          Apply total shift in meters for UTM files
            XDMNSHFT = XDMIN + (XDMNDIFM*NADFLG)
            YDMNSHFT = YDMIN + (YDMNDIFM*NADFLG)
            XDMXSHFT = XDMAX + (XDMXDIFM*NADFLG)
            YDMXSHFT = YDMAX + (YDMXDIFM*NADFLG)
            IF (NADFLG .EQ. 0.0D0) THEN
C*             Use original domain zones if no NAD shift
               ZMINSHFT = ZONMIN
               ZMAXSHFT = ZONMAX
            ELSE
C*             Use UTM zones based on shifted domain
               ZMINSHFT = ZONMIN_SHFT
               ZMAXSHFT = ZONMAX_SHFT
            END IF
         END IF

C*       Loop through profiles to extract data within domain
C*       for this DEM file:

  777    CONTINUE

C*       Allocate array for number of nodes, NUMNOD(NPROF(IDEM))
C*       First deallocate if needed
         IF (ALLOCATED(NUMNOD)) DEALLOCATE(NUMNOD)
         ALLOCATE (NUMNOD(NPROF(IDEM)))

         DO IPROF = 1, NPROF(IDEM)

C*          Initialize Variables
            NEWPRO = .TRUE.
            XNODEI = 0.0D0
            YNODEI = 0.0D0
            NUMNOD(IPROF) = 0

C*          Read Elevation From Each Profile Record.
            CALL GETPRO(IDEM,IPROF,IERR,RC)

C*          Check for error condition from GETPRO
            IF (IERR .NE. 0) GO TO 999

C*          Check to see if this profile lies between the min
C*          and maximum horizontal coords (LONG or UTMX)
            IF (.NOT. GOTDOMFLG) THEN
C*             No user-specified domain; use all profiles in data files
               GO TO 888
            END IF

            IF (TYPDEM .EQ. 'LAT' .AND. IPLAN(IDEM) .EQ. 0) THEN
C*             DEM/DOMAIN Type matches data file type
               IF (XBASE(IPROF) .LT. XDMNSHFT) THEN
                  CYCLE                        ! Cycle to next profile
               ELSE IF (XBASE(IPROF) .GT. XDMXSHFT) THEN
                  EXIT                         ! Exit profile loop
               END IF
            ELSE IF (TYPDEM .EQ. 'UTM' .AND. IPLAN(IDEM) .EQ. 1) THEN
C*             DEM/DOMAIN Type matches data file type
               IF (ZMINSHFT .EQ. ZMAXSHFT) THEN
                  IF (XBASE(IPROF) .LT. XDMNSHFT) THEN
                     CYCLE                      ! Cycle to next profile
                  ELSE IF (XBASE(IPROF) .GT. XDMXSHFT) THEN
                     EXIT                       ! Exit profile loop
                  END IF
               ELSE IF (IZO(IDEM) .EQ. ZMINSHFT) THEN
                  IF (XBASE(IPROF) .LT. XDMNSHFT) THEN
                     CYCLE                      ! Cycle to next profile
                  END IF
               ELSE IF (IZO(IDEM) .EQ. ZMAXSHFT) THEN
                  IF (XBASE(IPROF) .GT. XDMXSHFT) THEN
                     EXIT                       ! Exit profile loop
                  END IF
               END IF
            ELSE IF (TYPDEM .EQ. 'UTM' .AND. IPLAN(IDEM) .EQ. 0) THEN
C*             DEM/DOMAIN Type does NOT match data file type, convert xbase
C*             Calculate the Coords in UTM relative to NADA ellipsoid

               SELECT CASE (NADD(IDEM))
                 CASE (0)
                   ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                 CASE (2:4)
                   ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                 CASE (1,5:6)
                   ISPHERE = 0   ! UTM/LL conversions based on Clarke 1866 ellipsoid
                 CASE DEFAULT
                   ISPHERE = 4   ! DEFAULT CASE shouldn't occur
                 END SELECT

               ARGE  = XBASE(IPROF)
               ARGN  = YBASE(IPROF)
               IZONE = 0
               XARG  = 0.0D0
               YARG  = 0.0D0
               CALL UTMGEO (555,0,IZONE,XARG,YARG,ARGE,ARGN,ISPHERE)
               XBUTM = XARG
               YBUTM = YARG
               IF (ZMINSHFT .EQ. ZMAXSHFT) THEN
                  IF (XBUTM .LT. XDMNSHFT) THEN
                     CYCLE                     ! Cycle to next profile
                  ELSE IF (XBUTM .GT. XDMXSHFT) THEN
                     EXIT                      ! Exit profile loop
                  END IF
               ELSE IF (IZONE .EQ. ZMINSHFT) THEN
                  IF (XBUTM .LT. XDMNSHFT) THEN
                     CYCLE                     ! Cycle to next profile
                  END IF
               ELSE IF (IZONE .EQ. ZMAXSHFT) THEN
                  IF (XBUTM .GT. XDMXSHFT) THEN
                     EXIT                      ! Exit profile loop
                  END IF
               END IF
            ELSE IF (TYPDEM .EQ. 'LAT' .AND. IPLAN(IDEM) .EQ. 1) THEN
C*             DEM/DOMAIN Type does NOT match data file type, convert xbase
C*             Calculate the Coords in Lat/Lon relative to NADA ellipsoid

               SELECT CASE (NADD(IDEM))
                 CASE (0)
                   ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                 CASE (2:4)
                   ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                 CASE (1,5:6)
                   ISPHERE = 0   ! UTM/LL conversions based on Clarke 1866 ellipsoid
                 CASE DEFAULT
                   ISPHERE = 4   ! DEFAULT CASE shouldn't occur
                 END SELECT

               XARG  = XBASE(IPROF)
               YARG  = YBASE(IPROF)
               IZDUM = 0
               ARGE  = 0.0D0
               ARGN  = 0.0D0
               CALL UTMGEO (333,IZO(IDEM),IZDUM,XARG,YARG,
     &                                          ARGE,ARGN,ISPHERE)
               XBLON = ARGE
               YBLAT = ARGN
               IF (XBLON .LT. XDMNSHFT) THEN
                  CYCLE                        ! Cycle to next profile
               ELSE IF (XBLON .GT. XDMXSHFT) THEN
                  EXIT                         ! Exit profile loop
               END IF
            END IF

  888       CONTINUE

C*          For Each Node in this Profile
            DO INOD = 1, NODES(IPROF)

C*             Determine the Coordinates of this Point (L/L)   ---  CALL GETNOD
               CALL GETNOD(IPROF,INOD)

               IF (.NOT. GOTDOMFLG) THEN
C*                No user-specified domain; use all nodes in each profile

C*                Increment the Number of Nodes in this Profile
C*                Inside the Domain
                  NUMNOD(IPROF) = NUMNOD(IPROF)+1

                  IF (NEWPRO) THEN

C*                   NUMPRF is the number of profiles written to
C*                   the direct access file.
                     NUMPRF(IDEM) = NUMPRF(IDEM) + 1

                     NSTART = INOD
                     NEWPRO = .FALSE.

                     XNODEI = XNODE
                     YNODEI = YNODE

                  END IF

               ELSE

C*                Check to See if This Node Lies Between the Min
C*                and Maximum vertical Coords (LAT or UTMY)

                  IF (TYPDEM .EQ. 'UTM' .AND. IPLAN(IDEM) .EQ. 0) THEN
C*                   DEM/DOMAIN Type does NOT match data file type, convert xbase
C*                   Calculate the Coords in UTM relative to NADA ellipsoid

                     SELECT CASE (NADD(IDEM))
                       CASE (0)
                         ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                       CASE (2:4)
                         ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                       CASE (1,5:6)
                         ISPHERE = 0   ! UTM/LL conversions based on Clarke 1866 ellipsoid
                       CASE DEFAULT
                         ISPHERE = 4   ! DEFAULT CASE shouldn't occur
                       END SELECT

                     ARGE  = XBASE(IPROF)
                     ARGN  = YNODE
                     IZONE = 0
                     XARG  = 0.0D0
                     YARG  = 0.0D0
                     CALL UTMGEO (555,0,IZONE,XARG,YARG,ARGE,ARGN,
     &                                                       ISPHERE)
                     XNUTM = XARG
                     YNUTM = YARG

                  ELSE IF (TYPDEM.EQ.'LAT' .AND. IPLAN(IDEM).EQ.1) THEN
C*                   DEM/DOMAIN Type does NOT match data file type, convert xbase
C*                   Calculate the Coords in Lat/Lon relative to NADA ellipsoid

                     SELECT CASE (NADD(IDEM))
                       CASE (0)
                         ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                       CASE (2:4)
                         ISPHERE = 4   ! UTM/LAT-LONG conversions based on GRS80 ellipsoid
                       CASE (1,5:6)
                         ISPHERE = 0   ! UTM/LL conversions based on Clarke 1866 ellipsoid
                       CASE DEFAULT
                         ISPHERE = 4   ! DEFAULT CASE shouldn't occur
                       END SELECT

                     XARG  = XBASE(IPROF)
                     YARG  = YNODE
                     IZDUM = 0
                     ARGE  = 0.0D0
                     ARGN  = 0.0D0
                     CALL UTMGEO (333,IZO(IDEM),IZDUM,XARG,YARG,
     &                                                ARGE,ARGN,ISPHERE)
                     XBLON = ARGE
                     YBLAT = ARGN
                  END IF

                  IF ((TYPDEM.EQ.'UTM'.AND.IPLAN(IDEM).EQ.0 .AND.
     &                 YNUTM .GE. YDMNSHFT .AND.
     &                 YNUTM .LE. YDMXSHFT) .OR.
     &
     &                (TYPDEM.EQ.'LAT'.AND.IPLAN(IDEM).EQ.1 .AND.
     &                 YBLAT .GE. YDMNSHFT .AND.
     &                 YBLAT .LE. YDMXSHFT) .OR.
     &
     &                (YNODE .GE. YDMNSHFT .AND.
     &                 YNODE .LE. YDMXSHFT) ) THEN

C*                   This node is within the domain

C*                   Increment the Number of Nodes in this Profile
C*                   Inside the Domain
                     NUMNOD(IPROF) = NUMNOD(IPROF)+1

                     IF (NEWPRO) THEN

C*                      NUMPRF is the number of profiles written to
C*                      the direct access file.
                        NUMPRF(IDEM) = NUMPRF(IDEM) + 1

                        NSTART = INOD
                        NEWPRO = .FALSE.

                        XNODEI = XNODE
                        YNODEI = YNODE

                     END IF

                  ELSE IF ((TYPDEM.EQ.'UTM'.AND.IPLAN(IDEM).EQ.0 .AND.
     &                      YNUTM .GT. YDMXSHFT) .OR.
     &
     &                     (TYPDEM.EQ.'LAT'.AND.IPLAN(IDEM).EQ.1 .AND.
     &                      YBLAT .GT. YDMXSHFT) .OR.
     &
     &                     (YNODE .GT. YDMXSHFT) ) THEN

C*                   This node is beyond the domain

                     EXIT                    ! Exit loop on nodes

                  END IF

               END IF

C           End Loop on Nodes
            END DO

C*          If there were any nodes in this profile that were inside
C*          the domain, determine record numbers and write out data
C*          to direct access and index files
            IF (NUMNOD(IPROF) .GT. 0) THEN
C*             Write the baseline coords and number of nodes
C*             in this profile to the master index file
               WRITE (IDXUNT(IDEM),IDXFRM) XNODEI, YNODEI,
     &                NUMNOD(IPROF), IZO(IDEM), MAXEL(IPROF)

C*             For each node in the profile inside the domain
               DO N = 1, NUMNOD(IPROF)

C*                Calculate record number for the node
                  JRECRD = JRECRD + 1

C*                Extract elevation for this node from array
                  ELV = DBLE(ZDEM(NSTART-1+N))

C ---             Convert elevation units and apply other adjustments, except for missing
                  IF (ELV .GT. -9000.0D0) THEN
C*                   Apply vertical resolution factor from header (DCI) to elevations
C                    from file, adjust for elevation of local datum (LOCEL), and
C                    convert from feet to meters, if necessary
                     IF (ELUNIT(IDEM) .EQ. 1) THEN
C                       Convert from feet to meters
                        ELV = ((ELV* DCI(IDEM))+ LOCEL(IPROF))* 0.3048D0
                     ELSE
C                       Elevations in meters
                        ELV = (ELV* DCI(IDEM)) + LOCEL(IPROF)
                     END IF
                  END IF

C*                Write the elevations to a direct access file.
                  CALL WRITEZ (IDEM,JRECRD,ELV)

C              End Loop on Nodes
               END DO

            END IF

C*       End loop on profiles
         END DO

         IF (GOTDOMFLG) THEN
C ---       Write information to DOMDETAIL.OUT debug file
C           regarding direct access files
            IF (IDEM .EQ. 1) THEN
               WRITE(DOMK,50) NUMDEM
  50           FORMAT('From INITER_DEM:',
     &              //'Information on Direct Access Files Within ',
     &                                    'User-specified Domain:',
     &              //'No. of DEM Files = ',I6/)
            END IF
            WRITE(DOMK,60) IDEM, NUMPRF(IDEM)
  60        FORMAT('Direct Access File for DEM File No. ',I6,
     &             ' Contains ',I7,' Profiles.')
C ---       If direct access is not empty, NUMPRF(IDEM) > 0,
C           loop through profiles to print out first and last
            IF (NUMPRF(IDEM) .GT. 0) THEN
               L_FirstProf = .TRUE.
               DO K = 1, NPROF(IDEM)
                  IF (L_FirstProf .AND. NUMNOD(K) .GT. 0) THEN
C ---                This is the first non-empty profile,
C                    or the first profile within the direct access file
                     WRITE(DOMK,70) XBASE(K),YBASE(K),NUMNOD(K)
  70                 FORMAT(/' XBASE, YBASE and NODES for First',
     &                                        ' Profile: ',2F16.4,I7)
                     L_FirstProf = .FALSE.
                  ELSE IF ((.NOT.L_FirstProf .AND. NUMNOD(K).EQ.0)) THEN
C ---                The last non-empty profile has just been passed,
C                    print out info for previous profile
                     WRITE(DOMK,80) XBASE(K-1),YBASE(K-1),NUMNOD(K-1)
  80                 FORMAT(' XBASE, YBASE and NODES for Last ',
     &                                        ' Profile: ',2F16.4,I7/)
                     EXIT
                  ELSE IF (K .EQ. NPROF(IDEM)) THEN
C ---                This is the last profile in the DEM file, this
C                    must also be the last profile in the direct access file
                     WRITE(DOMK,80) XBASE(K),YBASE(K),NUMNOD(K)
                     EXIT
                  END IF
               END DO
            ELSE
               WRITE(DOMK,*)
            END IF
         END IF

         CLOSE (IDXUNT(IDEM))
         CLOSE (IDRUNT(IDEM))
         CLOSE (IDMUNT(IDEM))
         JRECRD = 0

         CYCLE DEMLOOP

 95      CONTINUE
C        Error opening a DEM file

         RUNERR = .TRUE.
         WRITE(DUMMY,'("DEM#",I4)') MIN(IDEM,9999)
         CALL ERRHDL(PATH,MODNAM,'E','500',DUMMY)

         CLOSE (IDXUNT(IDEM))
         CLOSE (IDRUNT(IDEM))
         CLOSE (IDMUNT(IDEM))
         JRECRD = 0

         CYCLE DEMLOOP

 97      CONTINUE
C        Error reading a DEM file

         RUNERR = .TRUE.
         WRITE(DUMMY,'("DEM#",I4)') MIN(IDEM,9999)
         CALL ERRHDL(PATH,MODNAM,'E','510',DUMMY)

         CLOSE (IDXUNT(IDEM))
         CLOSE (IDRUNT(IDEM))
         CLOSE (IDMUNT(IDEM))
         JRECRD = 0

         CYCLE DEMLOOP

C     End loop on dem files
      END DO DEMLOOP

C     Close DOMDETAIL.OUT debug file
      CLOSE(DOMK)

      WRITE(*,*) 'Exiting INITER_DEM'

999   RETURN
      end subroutine

      SUBROUTINE GETPRO (IDM,IPROF,IERR, RC)
C***********************************************************************
C*       PURPOSE:  THIS SUBROUTINE WILL READ X,Y NODE OF THE BASELINE AND
C*              THE Z VALUES FROM A DEM PROFILE
C*
C*       PROGRAMMER: Jayant Hardikar, Roger Brode
C*
C*       DATE:    September 29, 1995
C*
C*       MODIFIED: June 30, 2008
C*
C*                Added error checking for maximum number of nodes,
C*                and modified read format to accept up to 5,000 nodes.
C*                Roger W. Brode, U.S. EPA, OAQPS, AQMG
C*
C*       INPUTS:  Raw Terrain Data (DEM) Files
C*
C*       OUTPUTS: Direct Access Binary Terrain Files
C*
C*       CALLED FROM:   MAIN
C***********************************************************************

C*    Variable Declarations
      USE MAIN1

      IMPLICIT NONE

      SAVE

      CHARACTER*12 MODNAM

      INTEGER IERR, IDM, IPROF, ROWNUM, COLNUM, COLZ
      INTEGER L, J1, J2, KK, RR, RC

      MODNAM = 'GETPRO'
      IERR   = 0

      NODES(IPROF) = 0

      IF (FT(IDM) .NE. 1) THEN
         READ (IDMUNT(IDM), 208, ERR = 99, END = 99)
     &         ROWNUM, COLNUM, NODES(IPROF)
  208    FORMAT (3I6)

C*       Deallocate and reallocate ZDEM array for this profile
         IF (ALLOCATED(ZDEM)) DEALLOCATE(ZDEM)
         ALLOCATE(ZDEM(NODES(IPROF)))
         BACKSPACE (IDMUNT(IDM))

         READ (IDMUNT(IDM),209,ERR=99,END=99)
     &         ROWNUM, COLNUM,
     &         NODES(IPROF), COLZ, XBASE(IPROF), YBASE(IPROF),
     &         LOCEL(IPROF), MINEL(IPROF), MAXEL(IPROF),
     &         (ZDEM(L), L=1,NODES(IPROF))
C*        Read format supports up to 100,000 nodes per profile
  209     FORMAT (2I6,2I6,2D24.15,D24.15,2D24.15,
     &            146I6,587(:/170I6),:/64I6,:)

      ELSE

         RC = RC + 1
         READ (IDMUNT(IDM),FMT=209, REC=RC, ERR = 99)
     &         ROWNUM, COLNUM, NODES(IPROF)

C*       Deallocate and reallocate ZDEM array for this profile
         IF (ALLOCATED(ZDEM)) DEALLOCATE(ZDEM)
         ALLOCATE(ZDEM(NODES(IPROF)))

         J1 = 1
         RR = 0
         IF (NODES(IPROF) .GT. 146) THEN
            J2 = 146
            RR = INT((NODES(IPROF) - 146)/170 + 1)
            IF (MODULO((NODES(IPROF) - 146),170) .EQ. 0) RR = RR - 1
         ELSE
            J2 = NODES(IPROF)
         END IF

         READ (IDMUNT(IDM),FMT=209, REC=RC, ERR = 99)
     &         ROWNUM, COLNUM,
     &         NODES(IPROF), COLZ, XBASE(IPROF), YBASE(IPROF),
     &         LOCEL(IPROF), MINEL(IPROF), MAXEL(IPROF),
     &         (ZDEM(L), L=J1,J2)

         IF (RR .NE. 0) THEN
            DO KK = 1 , RR
              RC = RC + 1
              J1 = J2 + 1
              J2 = J1 + 170 - 1
              IF (J2 .GT. NODES(IPROF)) J2 = NODES(IPROF)
              READ (IDMUNT(IDM), FMT=211, REC=RC, ERR = 99)
     &             (ZDEM(L), L=J1,J2)
  211          FORMAT(170I6)
            END DO
         END IF
      END IF

      GO TO 999

99    RUNERR = .TRUE.
      IERR   = -1
      WRITE(DUMMY,'("DEM#",I4)') MIN(IDM,9999)
      CALL ERRHDL(PATH,MODNAM,'E','510',DUMMY)

999   RETURN

      end subroutine

      SUBROUTINE GETNOD(IPR,INO)
C***********************************************************************
C*               GETNOD Module of AERMAP Terrain Preprocessor
C*
C*       PURPOSE: Determine the Coordinates of a DEM Point (L/L)
C*
C*       PROGRAMMER: Jayant Hardikar, Roger Brode
C*
C*       DATE:    September 29, 1995
C*
C*       MODIFIED: June 30, 2008
C*
C*                Corrected to use DYM, node spacing in 'Y' direction,
C*                to assign YNODE values in profile.
C*                Roger W. Brode, U.S. EPA, OAQPS, AQMG
C*
C*       INPUTS:  Profile Number, Node Number, Coordinates of Baseline
C*
C*       OUTPUTS: Coordinates of a DEM Point
C*
C*       CALLED FROM:   INITER
C***********************************************************************

C*    Variable Declarations
      USE MAIN1

      IMPLICIT NONE

      CHARACTER*12 MODNAM

      INTEGER INO, IPR

      MODNAM = 'GETNOD'

C*    For the baseline in a profile, the coords of a node are the same
c*    as the coords of the baseline
      IF (INO .EQ. 1) THEN
         XNODE = XBASE(IPR)
         YNODE = YBASE(IPR)

C*    Otherwise, calculate the coords based on the profile and node number
      ELSE
         XNODE = XBASE(IPR)
         YNODE = YNODE + DYM(IDEM)
      END IF

      RETURN
      end subroutine

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
module tensor_converter_interface
    use, intrinsic :: iso_fortran_env, only : int64

    use :: tensor_package_api, only : &
            tensor, &
            allocate_tensor_object

    implicit none
    private

    public :: allocate_and_reshape_tensor
    public :: reshape_tensor

contains
    subroutine allocate_and_reshape_tensor(dst, src, dims)
        class(tensor), allocatable, intent(inout) :: dst
        class(tensor), intent(in) :: src
        integer(int64), dimension(:), intent(in) :: dims

        call allocate_tensor_object(dst, dims)
        call reshape_tensor(dst, src, dims)
    end subroutine allocate_and_reshape_tensor

    subroutine reshape_tensor(dst, src, dims)
        class(tensor), intent(inout) :: dst
        class(tensor), intent(in) :: src
        integer(int64), dimension(:), intent(in) :: dims

        if ( dst%rank /= size(dims) ) &
                error stop "tensor_converter_interface::reshape_tensor:Error in rank."
        if ( src%number_of_elements /= product(dims) ) &
                error stop "tensor_converter_interface::reshape_tensor:Error in size"

        call dst%setup(src%storage, src%datatype, dims)
    end subroutine reshape_tensor
end module tensor_converter_interface

max_stars_repo_name: cheshyre/ntcl-tensor

max_issues_repo_name: cheshyre/ntcl-tensor

max_forks_repo_name: cheshyre/ntcl-tensor


program test_get_grid_x

  use bmif_2_0, only: BMI_FAILURE
  use bmiheatf
  use fixtures, only: config_file, status

  implicit none

  integer, parameter :: grid_id = 1
  integer, parameter :: nx = 1
  double precision, parameter, dimension(nx) :: expected_x = (/ 0.0 /)

  type (bmi_heat) :: m
  double precision, dimension(nx) :: grid_x
  integer :: i

  status = m%initialize(config_file)
  status = m%get_grid_x(grid_id, grid_x)
  status = m%finalize()

  print *, grid_x
  print *, expected_x

  do i = 1, nx
     if (grid_x(i).ne.expected_x(i)) then
        stop BMI_FAILURE
     end if
  end do
end program test_get_grid_x


text:
! P03-complicated-test-1.f90
! The ROSE front end will fail an assertion under the following
! detailed conditions. All of these must hold to provoke the bug:
!
! (1) A derived type t1 must be declared before a derived type t2.
! (2) Type t1 must have a component f of type 'pointer to array of t2'.
! (3) Types t1 and t2 must each have a component c with the same name,
!     where c != f.
! (4) Type t2 must have an additional component d.
! (5) Following the type declarations must be two statements, each
!     referencing a component of t2 by indexing t1's component f, where:
!     (a) the first data reference accesses t2's additional component d
!     (b) the second data reference accesses t2's same-named component c.
!
! NOTE: this is the same as P04, but with the data references in 5(a)
! and 5(b) interchanged.

program p

   type :: t1
      integer :: c
      type (t2), pointer, dimension(:)  :: f
   end type

   type :: t2
      integer :: c
      integer :: d
   end type

   type (t1) :: v

   v % f(1) % d = 0
   v % f(1) % c  = 0    ! assertion failure: expression stack not empty

end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program cellular_automata
    use :: cellular_automata_mod
    implicit none
    integer :: nr_cells, rule_nr, max_steps
    type(automaton_t) :: automaton
    integer :: step

    call get_parameters(rule_nr, nr_cells, max_steps)
    automaton = init_automaton(nr_cells, rule_nr)
    call print_rules(automaton)
    call print_automaton(automaton)
    do step = 1, max_steps
        call step_automaton(automaton)
        call print_automaton(automaton)
    end do

contains

    subroutine get_parameters(rule_nr, nr_cells, max_steps)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        integer, intent(out) :: rule_nr, nr_cells, max_steps
        integer :: istat
        character(len=1024) :: buffer, msg

        rule_nr = 90
        nr_cells = 50
        max_steps = 100

        if (command_argument_count() > 0) then
            call get_command_argument(1, buffer)
            read (buffer, fmt=*, iostat=istat, iomsg=msg) rule_nr
            if (istat /= 0) then
                print '(2A)', 'error: ', msg
                stop 1
            end if
        end if
        if (command_argument_count() > 1) then
            call get_command_argument(2, buffer)
            read (buffer, fmt=*, iostat=istat, iomsg=msg) nr_cells
            if (istat /= 0) then
                print '(2A)', 'error: ', msg
                stop 1
            end if
        end if
        if (command_argument_count() > 2) then
            call get_command_argument(3, buffer)
            read (buffer, fmt=*, iostat=istat, iomsg=msg) max_steps
            if (istat /= 0) then
                print '(2A)', 'error: ', msg
                stop 1
            end if
        end if
    end subroutine get_parameters

end program cellular_automata

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
Subroutine ReadLTMatrix(iU,iSkip0,N,LDA,A,nleft,ncol,lencol,iskip)

Implicit Real(8) (A-H,O-Z)

Real(8) A(LDA,N)

character(255) str

! Subroutine ReadLTMatrix reads low triangular matrix A from file iU
!   Reading statrts from the current position of file iU
!   iSkip0  - number of lines to be skipped before reading of first row (the row with element A(1,1))
!   NLEFT - number of symbols to be skipped from the beginning of each line.
!   NCOL  - number of columns.
!   LENCOL - nuber of symbols in each colums
!   ISKIP  - number of lines to be skipped after next portion of columns
!   A      - matrix to be read in. It will be symmetrized after reading. All values are read by format G<lencol>.0

Do i=1,iSkip0
    Read(iU,*)
Enddo

Nport=N/ncol
Nlast=N-Nport*ncol
If (Nlast>0) Nport=Nport+1

jend=0
Do ip=1,Nport
    jbeg=jend+1
    jend=jend+ncol
    If (jend>N) jend=N
    Do i=jbeg,N
        Read(iU,'(<nleft>x,<ncol>G<lencol>.0)')(A(i,j),j=jbeg,jend)
    Enddo
    Do k=1,iskip
        If (.not.EOF(iU)) Read(iU,*)
    Enddo
Enddo

Do i=1,N
    Do j=1,i-1
        A(j,i)=A(i,j)
    Enddo
Enddo

End

max_stars_repo_name: skignatov/ReadExcitedStates

max_issues_repo_name: skignatov/ReadExcitedStates

max_forks_repo_name: skignatov/ReadExcitedStates


text:
      SUBROUTINE LOADSU
C
C     LOADSU SETS UP LOAD INFOTMATION FOR PROLAT FROM NSLT.
C     Z(IST)IS THE STARTING POINT FOR OPEN CORE,Z(MCORE) IS THE LAST
C     AVAILABLE WORD, NTOT IS THE NUMBER OF WORDS PUT INTO OPEN CORE
C     BY THIS ROUTINE. LOAD IS THE LOAD ID.
C
      LOGICAL         REMFL
      INTEGER         SUBCAS,BUF2,SCR1,FILE,HEST,BGPDT
      DIMENSION       NWORDS(19),MCB(7),IZ(1),L(2),ZL(2),NAM(2)
      CHARACTER       UFM*23,UWM*25
      COMMON /XMSSG / UFM,UWM
      COMMON /ZZZZZZ/ Z(1)
      COMMON /SYSTEM/ SYSBUF,IOUT
      COMMON /BIOT  / NG1,NG2,IST,SUBCAS,X1,Y1,Z1,X2,Y2,Z2,BUF2,REMFL,
     1                MCORE,LOAD,NSLT,SCR1,HEST,NTOT
      EQUIVALENCE     (Z(1),IZ(1)),(L(1),ZL(1))
      DATA    NAM   / 4HLOAD,4HSU  /
      DATA    NWORDS/ 6,6,4,4,6,6,2,5,5,6,6,7,12,10,10,19,38,7,5/
C
      BGPDT  = 103
      MCB(1) = BGPDT
      CALL RDTRL (MCB)
      NROWSP = MCB(2)
      MCB(1) = HEST
      CALL RDTRL (MCB)
      NEL    = MCB(2)
      NSIMP  = 0
      FILE   = NSLT
      CALL OPEN (*1001,NSLT,Z(BUF2),0)
      CALL READ (*1002,*10,NSLT,Z(IST+1),MCORE,0,IWORDS)
      GO TO 1008
   10 NLOADS = IWORDS-2
C
C     CHECK LOAD SELECTION AGAINST SIMPLE LOAD ID-S
C
      IF (NLOADS .EQ. 0) GO TO 35
      DO 20 I = 1,NLOADS
      IF (IZ(IST+2+I) .EQ. LOAD) GO TO 80
   20 CONTINUE
C
C     NOT A SIMPLE LOAD-MUST BEA LOAD COMBINATION. SKIP NLOADS RECORDS
C     AND SEARCH FOR PROPER LOAD ID
C
      DO 30 I = 1,NLOADS
      CALL FWDREC (*1002,NSLT)
   30 CONTINUE
C
C     READ 2 WORDS AT A TIME -1,-1 SIGNIFIES END OF LOAD CARD
C
   35 ILOAD = IST + IWORDS
   40 CALL READ (*1002,*500,NSLT,L,2,0,IFLAG)
      IF (L(1) .EQ. LOAD) GO TO 60
C
C     NO MATCH-SKIP TO -1-S
C
   50 CALL FREAD (NSLT,L,2,0)
      IF (L(1).EQ.-1 .AND. L(2).EQ.-1) GO TO 40
      GO TO 50
C
C     MATCH
C
   60 ALLS  = ZL(2)
   70 CALL FREAD (NSLT,L,2,0)
      IF (L(1).EQ.-1 .AND. L(2).EQ.-1) GO TO 90
      NSIMP = NSIMP + 1
      IF (ILOAD+2*NSIMP .GT. MCORE) GO TO 1008
      ISUB  = 2*NSIMP - 1
      Z(ILOAD+ISUB) = ZL(1)
      IZ(ILOAD+ISUB+1) = L(2)
      GO TO 70
C
C     WE HAVE NSIMP SIMPLE LOADS. FOR ONE LOAD,SET PROPER PARAMETERS
C
   80 NSIMP = 1
      ALLS  = 1.
      ILOAD = IST + IWORDS
      Z(ILOAD+1) = 1.
      IZ(ILOAD+2) = LOAD
C
C     FOR EACH SIMPLE LOAD, FIND PROPER LOAD ID AND THEN POSITION TO
C     PROPER LOAD RECORD IN NSLT
C
   90 NTOT  = 0
      ISIMP = ILOAD + 2*NSIMP
      DO 270 NS = 1,NSIMP
C
      ISUB   = ILOAD + 2*NS - 1
      FACTOR = Z(ISUB)
      ID     = IZ(ISUB+1)
      NCARDS = 0
      CALL REWIND (NSLT)
      I = 1
      IF (NLOADS .EQ. 0) GO TO 110
      DO 100 I = 1,NLOADS
      IF (ID .EQ. IZ(IST+2+I)) GO TO 110
  100 CONTINUE
      GO TO 499
C
  110 DO 120 J = 1,I
      CALL FWDREC (*1002,NSLT)
  120 CONTINUE
C
  125 CALL READ  (*1002,*260,NSLT,NOBLD,1,0,IFLAG)
      CALL FREAD (NSLT,IDO,1,0)
      IF (ISIMP+2 .GT. MCORE) GO TO 1008
      IZ(ISIMP+1) = NOBLD
      IZ(ISIMP+2) = IDO
      ISIMP = ISIMP + 2
      NTOT  = NTOT + 2
C
C     SKIP NOBLD=-20. IF NOBLD=24(REMFLUX), STORE ONLY NOBLD AND IDO,
C     BUT SKIP REMFLUX INFO ON NSLT
C
      IF (NOBLD .EQ. -20) GO TO 250
      IF (NOBLD .LE.  19) GO TO 245
      KTYPE = NOBLD - 19
      GO TO (126,127,128,129,130), KTYPE
  126 MWORDS = 3*NROWSP
      GO TO 140
  127 MWORDS = 12
      GO TO 140
  128 MWORDS = 48
      GO TO 140
  129 MWORDS = 9
      GO TO 140
  130 MWORDS = 3*NEL
      MWORDS = -MWORDS
      GO TO 141
C
  140 IF(ISIMP+MWORDS*IDO .GT. MCORE) GO TO 1008
      NTOT = NTOT + MWORDS*IDO
  141 DO 240 J = 1,IDO
C
C     NCARDS TELLS HOW MANY SIMPLE LOAD CARDS HAVE THE PRESENT FACTOR
C     APPLIED TO IT
C
      NCARDS = NCARDS + 1
      CALL FREAD (NSLT,Z(ISIMP+1),MWORDS,0)
      IF (NOBLD .NE. 24) ISIMP = ISIMP + MWORDS
  240 CONTINUE
C
C     DONE WITH CARDS OF PRESENT TYPE-GET ANOTHER TYPE
C
      GO TO 125
C
C     TYPE=-20    SKIP IT
C
  250 CALL FREAD (NSLT,Z,-(3*NROWSP),0)
      GO TO 125
C
C     NOT A MAGNETICS TYPE OF LOAD. - SKIP IT
C
  245 WRITE  (IOUT,246) UWM,LOAD
  246 FORMAT (A25,', IN FUNCTIONAL MODULE PROLATE, LOAD SET',I8, /5X,
     1       'CONTAINS A NONMAGNETIC LOAD TYPE. IT WILL BE IGNORED.')
      DO 247 I = 1,IDO
      CALL FREAD (NSLT,Z,-NWORDS(NOBLD),0)
  247 CONTINUE
C
C     EOR ON NSLT-DONE WITH THIS SIMPLE LOAD-GET ANOTHER SIMPLE LOAD
C
C     SUBSTITUTE IN OPEN CORE NCARDS FOR THE SIMPLE LOAD ID. WE NO
C     LONGER NEED THE ID, BUT WE MUST SAVE NCARDS
C
  260 CONTINUE
      IZ(ISUB+1) = NCARDS
C
  270 CONTINUE
C
C     DONE
C
C     STORE ALL THIS INFO BACK AT Z(IST) AS FOLLOWS
C
C     ALLS,NSIMP,(LOAD FACTOR,NCARDS) FOR EACH SIMPLE LOAD ID,
C     ALL LOAD INFO FOR EACH SIMPLE LOAD STARTING WITH NOBLD AND IDO
C
      Z(IST+1)  = ALLS
      IZ(IST+2) = NSIMP
      NS2   = 2*NSIMP
      DO 280 I = 1,NS2
  280 Z(IST+2+I) = Z(ILOAD+I)
      ISUB1 = IST + NS2 + 2
      ISUB2 = ILOAD + 2*NSIMP
      DO 290 I = 1,NTOT
  290 Z(ISUB1+I) = Z(ISUB2+I)
      NTOT = NTOT + 2*NSIMP + 2
      CALL CLOSE (NSLT,1)
      RETURN
C
  499 LOAD = ID
  500 WRITE  (IOUT,501) UFM,LOAD
  501 FORMAT (A23,', CANNOT FIND LOAD',I8,' ON NSLT IN BIOTSV')
      CALL MESAGE (-61,0,0)
C
 1001 N =-1
      GO TO 1010
 1002 N =-2
      GO TO 1010
 1008 N =-8
      FILE = 0
 1010 CALL MESAGE (N,FILE,NAM)
      RETURN
      END


text:
C NCLFORTSTART
      SUBROUTINE DPOPREMAP(DST_ARRAY,MAP_WTS,DST_ADD,SRC_ADD,SRC_ARRAY,
     +                     NDST,NLINK,NW,NSRC,XMSG)

c written in f77 for the GNU f77 compiler for portability reasons

      IMPLICIT NONE
      INTEGER NLINK,NW,NDST,NSRC
      DOUBLE PRECISION MAP_WTS(NW,NLINK)
      DOUBLE PRECISION DST_ARRAY(NDST)
      DOUBLE PRECISION SRC_ARRAY(NSRC)
      INTEGER DST_ADD(NLINK)
      INTEGER SRC_ADD(NLINK)
      DOUBLE PRECISION XMSG
C NCLEND
      INTEGER N

      DO N = 1,NDST
c initialize
          DST_ARRAY(N) = XMSG
      END DO

      DO N = 1,NLINK
          IF (SRC_ARRAY(SRC_ADD(N)).NE.XMSG) THEN
              IF (DST_ARRAY(DST_ADD(N)).EQ.XMSG) THEN
                  DST_ARRAY(DST_ADD(N)) = SRC_ARRAY(SRC_ADD(N))*
     +                                    MAP_WTS(1,N)
              ELSE
                  DST_ARRAY(DST_ADD(N)) = DST_ARRAY(DST_ADD(N)) +
     +                                    SRC_ARRAY(SRC_ADD(N))*
     +                                    MAP_WTS(1,N)
              END IF
          END IF
      END DO

      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


   !################################################################################!
 !/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\!
!######################################################################################!
!##!                                                                                !##!
!##!                                                                                !##!
MODULE Coefficient_Vector_Module                                                    !##!
!##!                                                                                !##!
!##!________________________________________________________________________________!##!
!##!                                                                                !##!
!##!    Contains the functions and subroutines that concern the coefficient vector, !##!
!##!        most importantly, the main routine responsible for solving the linear   !##!
!##!        system of equations.                                                    !##!
!##!                                                                                !##!
!##!================================================================================!##!
!##!                                                                                !##!
!##!    Contains:                                                                   !##!
!##!                                                                                !##!
!##!    +101+   Allocate_Coefficient_Vector                                         !##!
!##!    +102+   Deallocate_Coefficeint_Vector                                       !##!
!##!                                                                                !##!
!##!    +201+   Calculate_Coefficient_Vector                                        !##!
!##!                                                                                !##!
!##!                                                                                !##!
!######################################################################################!
 !\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/!
   !################################################################################!

!*D*================================!
!                                   !
!           Dependencies            !
!                                   !
!===================================!

USE constants, &
            ONLY :  idp

USE Global_Variables_And_Parameters, &
            ONLY :  NUM_R_ELEMENTS, NUM_R_NODES,                                    &
                    DEGREE, L_LIMIT, MATRIX_FORMAT, LINEAR_SOLVER,                  &
                    Source_Vector,                                                  &
                    STF_MAT, STF_ELEM_VAL, STF_COL_PTR, STF_ROW_IND, STF_NNZ,       &
                    Coefficient_Vector,                                             &
                    Matrix_Cholesky_Factorized_Flag

USE Boundary_Conditions_Module, &
            ONLY :  DIRICHLET_BC, NEUMANN_BC, DIRICHLET_BC_CCS, NEUMANN_BC_CCS,     &
                    DIRICHLET_BC_CHOL

USE Linear_Solvers_And_Preconditioners, &
            ONLY :  PRECOND_CONJ_GRAD_FULL, PRECOND_CONJ_GRAD_CCS

USE Cholesky_Module,    &
            ONLY :  Cholesky_Factorization, CCS_Forward_Substitution, CCS_Back_Substitution

IMPLICIT NONE

CONTAINS
!+101+############################################################################!
!                                                                                !
!               ALLOCATE_COEFFS - Allocate the coefficent vectors                !
!                                                                                !
!################################################################################!
SUBROUTINE Allocate_Coefficient_Vector()

ALLOCATE( Coefficient_Vector(0:NUM_R_NODES-1, -L_LIMIT:L_LIMIT, 0:L_LIMIT))

END SUBROUTINE Allocate_Coefficient_Vector

!+102+###########################################################################!
!                                                                                !
!             DEALLOCATE_COEFFS - Allocate the coefficent vectors                !
!                                                                                !
!################################################################################!
SUBROUTINE Deallocate_Coefficient_Vector()

DEALLOCATE( Coefficient_Vector )

END SUBROUTINE Deallocate_Coefficient_Vector

!+201+##########################################################################!
!                                                                               !
!      Calculate_Coefficient_Vector - Solve linear system STF_MAT*Coeffs = Src  !
!                                                                               !
!###############################################################################!
SUBROUTINE Calculate_Coefficient_Vector()

INTEGER                                                                     ::  INFO, LDAB
INTEGER, DIMENSION(NUM_R_NODES)                                             ::  IPIV

REAL(KIND = idp)                                                            ::  SCALE_FACTOR
COMPLEX(KIND = idp), DIMENSION(0:NUM_R_NODES-1)                             ::  WORK_VEC, WORK_VECB

REAL(KIND = idp), ALLOCATABLE, DIMENSION(:,:)                               ::  WORK_MAT

REAL(KIND = idp), ALLOCATABLE, DIMENSION(:)                                 ::  WORK_ELEM_VAL

INTEGER                                                                     ::  NNZ
INTEGER                                                                     ::  l, m, k
INTEGER                                                                     ::  Guess_Flag

REAL(KIND = idp)       :: timea, timeb

IF (MATRIX_FORMAT =='FULL') THEN

    ALLOCATE (WORK_MAT(0:NUM_R_NODES-1, 0:NUM_R_NODES-1))

ELSE IF (MATRIX_FORMAT == 'CCS' ) THEN

    NNZ = NUM_R_ELEMENTS*(DEGREE + 1)*(DEGREE + 1) - NUM_R_ELEMENTS + 1
    ALLOCATE(WORK_ELEM_VAL(0:NNZ-1))

END IF

WORK_VECB = 0

IF (( LINEAR_SOLVER == "CHOL" ) .AND. (Matrix_Cholesky_Factorized_Flag .EQV. .FALSE.) ) THEN

    !
    !   This only needs to be done everytime the radial mesh is defined/redefined.
    !   This performs Cholesky factorization on the stiffness matrix and overwrites
    !   the stiffness matrix variables STF_ELEM_VAL, STF_ROW_IND, and STF_COL_PTR to
    !   represent the factorization matrix, L.  This matrix can then be reused to
    !   solve the linear system using forward and backward substitution.
    !

    CALL Cholesky_Factorization()
    Matrix_Cholesky_Factorized_Flag = .TRUE.

END IF

DO l = 0,L_LIMIT

    DO m = -l,l

        IF (LINEAR_SOLVER =='FULL') THEN
            !####################################!
            !
            !           Full Matrix Solver       !
            !
            !####################################!

        WORK_MAT = STF_MAT(:,:,l)
        WORK_VEC = Source_Vector(:,m,l)

        CALL DIRICHLET_BC(WORK_MAT, WORK_VEC, l, m)

        CALL NEUMANN_BC(l, WORK_VEC)

        CALL DGESV(NUM_R_NODES, 1, WORK_MAT, NUM_R_NODES, IPIV, WORK_VEC, NUM_R_NODES, INFO)
        IF (INFO > 0) THEN
            print*,"DGESV has failed with INFO = ",INFO
        END IF

        !CALL PRECOND_CONJ_GRAD_FULL(WORK_MAT, WORK_VEC)

        ELSE IF (LINEAR_SOLVER == 'CCS') THEN
            !#######################################################################!
            !                                                                       !
            !           CCS Preconditioned Conjugate Gradient Matrix Solver         !
            !                                                                       !
            !#######################################################################!

            WORK_ELEM_VAL = STF_ELEM_VAL(:,l)
            WORK_VEC = Source_Vector(:,m,l)

            CALL DIRICHLET_BC_CCS(  NUM_R_NODES,   NNZ, l, m,                               &
                                    WORK_ELEM_VAL, STF_COL_PTR, STF_ROW_IND, WORK_VEC)

            CALL NEUMANN_BC_CCS(    NUM_R_NODES, NNZ, l, m,                                 &
                                    WORK_ELEM_VAL, STF_COL_PTR, STF_ROW_IND, WORK_VEC)

            GUESS_FLAG = 0  !!! 0 Means no guess, 1 means guess

            IF ( 1 == 0) THEN

                WORK_VECB = WORK_VEC

                CALL PRECOND_CONJ_GRAD_CCS(NUM_R_NODES, NNZ, WORK_ELEM_VAL, STF_COL_PTR, STF_ROW_IND, &
                                            WORK_VECB, GUESS_FLAG, WORK_VECB,0)

                GUESS_FLAG = 1

            END IF

            CALL PRECOND_CONJ_GRAD_CCS(NUM_R_NODES, NNZ, WORK_ELEM_VAL, STF_COL_PTR, STF_ROW_IND,       &
                                        WORK_VEC,GUESS_FLAG, WORK_VECB,0)

        ELSE IF (LINEAR_SOLVER == "CHOL") THEN
            !#######################################################################!
            !                                                                       !
            !               CCS Cholesky Factorization Matrix Solver                !
            !                                                                       !
            !#######################################################################!

            WORK_VEC = Source_Vector(:,m,l)

            CALL DIRICHLET_BC_CHOL(  NUM_R_NODES, STF_NNZ, l, m,                                &
            STF_COL_PTR, STF_ROW_IND, WORK_VEC)

            CALL NEUMANN_BC_CCS(    NUM_R_NODES, STF_NNZ, l, m,                                 &
            WORK_ELEM_VAL, STF_COL_PTR, STF_ROW_IND, WORK_VEC)

            CALL CCS_Forward_Substitution(NUM_R_NODES, STF_NNZ, STF_ELEM_VAL(:,l), STF_COL_PTR, STF_ROW_IND, WORK_VEC )

            CALL CCS_Back_Substitution(NUM_R_NODES, STF_NNZ, STF_ELEM_VAL(:,l), STF_COL_PTR, STF_ROW_IND, WORK_VEC )

        END IF

        Do k = 0,NUM_R_NODES - 1
            Coefficient_Vector(k,m,l) = WORK_VEC(k)
        END DO

    END DO
END DO

IF (MATRIX_FORMAT =='FULL') THEN

    DEALLOCATE (WORK_MAT)

ELSE IF (MATRIX_FORMAT == 'CCS' ) THEN

    DEALLOCATE(WORK_ELEM_VAL)

END IF

END SUBROUTINE Calculate_Coefficient_Vector

END MODULE Coefficient_Vector_Module


text:
!
! Copyright (c) 2015, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!

program testnewline
parameter(N=20)
integer :: results(N)
integer :: expected(N)

character*5 c1, c2, c3
logical l1(3), l2(3)
integer*8 abc(20)
integer iosint(3)
integer myintread(20)
integer myintwrite(20)
       data results / -1, -1, 0, 0, -1, 0, 0, 0, &
                      -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0/
       data expected / -1, -1, 0, 0, -1, 0, 0, 0, &
                      -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0/

! new_line
open(28, file='t.txt', status='replace', access='stream', form='formatted')
write(28,'(a)') 'hello'//new_line('x')//'world'
close(28)
!

! end of file
open(28, file='t.txt', form='formatted')
read(28,100,iostat=ios) c1
l1(1) = is_iostat_end(ios)
l2(1) = is_iostat_eor(ios)
read(28,100,iostat=ios) c2
l1(2) = is_iostat_end(ios)
l2(2) = is_iostat_eor(ios)
read(28,100,iostat=ios) c3
l1(3) = is_iostat_end(ios)
l2(3) = is_iostat_eor(ios)
100 format(a5)
close(28)

results(1)=(c1 == 'hello')
results(2)=(c2 == 'world')
results(3)=l1(1)
results(4)=l1(2)
results(5)=l1(3)
results(6)=l2(1)
results(7)=l2(2)
results(8)=l2(3)

!end of record
open(28, file='t2.txt',status='replace',form='formatted')
write(28,'(i13)') myintwrite
close(28)

open(28, file='t2.txt', form='formatted')
read(28,200,advance='no',iostat=ios) myintread
200 format(i20)
close(28)

results(9)=ios
results(10)=is_iostat_end(ios)
results(11)=is_iostat_eor(ios)

open(28, file='t2.txt', form='formatted')
read(28,300,advance='no',iostat=ios) myintread
300 format(i13)
close(28)

results(12)=ios
results(13)=is_iostat_end(ios)
results(14)=is_iostat_eor(ios)

!testing array
iosint=0
iosint=is_iostat_end(iosint)
results(15:17)=iosint

iosint=is_iostat_eor(iosint)
results(18:20)=iosint

call check(results,expected,N)

end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


!======================================================================================
! PROGRAMA QUE COLOCA UN NUMERO N DE PARTICULAS Y LAS MUEVE CONCIDERANDO CONDICIONES
! PERIODICAS. SE HACE SEGUIMIENTO DEL MOVIMIENTO DE PARTICULAS (2 PARTICULAS).
! SE REALIZO PARA 3 DIMENSIONES (X,Y,Z). ARREGLO UNIFORME.
!
! AUTOR: MARTIN ALEJANDRO PAREDES SOSA
!======================================================================================
Program Main
  Use cte
  Implicit None
  Integer :: i, j, IStep, N3                                !CONTADORES
  Real :: VLRC, VI, V, VOld, VNew, DV, VN                   !ENERGIAS
  Real :: OldX, OldY,OldZ, NewX, NewY, NewZ                 !VALORES TEMP DE POSC
  Real :: RanX, RanY,RanZ, Dummy                            !VALORES ALEATORIOS
  Logical :: Ctrl                                           !CONTROL LOGICO
  Integer :: h, p, g                                        !VALOR ALEATORIO TRAZADORAS
  Real :: Dens2                                             !DENSIDAD  FINAL
  Real, Dimension(2) :: Temp                                !TEMPORAL PARA PARTICULAS ALEATORIAS

  !PEDIR DENSIDAD Y NUMERO DE PARTICULAS
  Write(*,*) "NUMERO DE PARTICULAS"
  Read(*,*) N
  Write(*,*) "CONCENTRACION REDUCIDA"
  Read(*,*) Dens
  Write(*,*) "NUMERO DE CICLOS"
  Read(*,*) NStep

  !Allocate( X(N), Y(N), Z(N)  )

  !GENERAR LA CONFIGURACION INICIAL
  Call ConfigIni(N3)
  Write(*,*) "CONFIGURACION INICIAL LISTA"
  RCut = BoxL / 2.0
  dRMax = 1.0

  !ABRIENDO ARCHIVOS PARA GUARDAR INFO DEL SISTEMA
  Open(2, File="ConFin.dat")

  !SELECCION DE PARTICULAS ALEATORIAMENTE
  Call Random_Number(Temp)
  Temp = Temp * N3
  h = Anint( Temp(1) )
  p = Anint( Temp(2) )
  !g = Anint( Temp(3) )
  Write(*,*) "SE SELECIONARON LAS PARTICULAS: ", h, p!, g
  Open(11, File= "Part1.dat")
  Open(12, File= "Part2.dat")
  !Open(13, File= "Part3.dat")

  Write(*,*) "INICIO MOVIMIENTO PARTICULAS"
  !MOVIMIENTO DE PARTICULAS ALEATORIA
  Configuracion: Do iStep = 1, NStep

     Particula: Do i = 1, N3

        OldX = X(i)
        OldY = Y(i)
        OldZ = Z(i)

        !GENERAR VALORES ALEATORIOS PARA MOV TENTATIVOS
        Call Random_Number(RanX)
        Call Random_Number(RanY)
        Call Random_Number(RanZ)

        !MOVIMIENTO DE PARTICULA
        NewX = OldX + (2.0*RanX - 1.0)*dRMax
        NewY = OldY + (2.0*RanY - 1.0)*dRMax
        NewZ = OldZ + (2.0*RanZ - 1.0)*dRMax

        !CONDICIONES PERIODICAS (MANTENER MISMA N EN TODA CONFIGURACION)
        NewX = NewX - BoxL*Anint(NewX/BoxL)
        NewY = Newy - BoxL*Anint(NewY/BoxL)
        Newz = NewZ - Boxl*Anint(NewZ/BoxL)

        X(i) = NewX
        Y(i) = NewY
        Z(i) = NewZ
     End Do Particula

     !SEGUIMIENTO DE PARTICULAS (TRAZADORAS)
     Write(11,*) X(h), Y(h), Z(h)
     Write(12,*) X(p), Y(p), Z(p)
     !Write(13,*) X(g), Y(g), Z(g)

  End Do Configuracion
  Write(*,*) "CONFIGURACIONES LISTAS"
  Close(11)
  Close(12)
  !Close(13)

  !CHECA LA DENSIDAD DE LA CELDA ORIGINAL (CONCENTRACION FINAL)
  !Call Conc(Dens2)

  !GUARDAR CONFIG FINAL
  ConfigFin: Do i=1, N3
     Write(2,*) X(i), Y(i), Z(i)
  End Do ConfigFin

  Deallocate( X, Y, Z )

End Program Main