text:
      subroutine fovrg (il, ihard, rmt, xmt, jri, e, nr, dx, ri, v, dny,
     1                  pu, qu, p, q, ps, qs, vm)
      implicit double precision (a-h, o-z)

c     Input:
c        il      ang mom number + 1
c        ihard   number of times convergence test fails
c        rmt     muffin tin radius
c        xmt     x such that rmt = exp ((x-1)*dx - 8.8)
c        jri     first interstitial grid point (imt + 1)
c        e       current complex energy
c        nr      number of points in r grid
c        dx      dx in Loucks' grid (usually .05)
c        ri(nr)  Loucks' position grid, r = exp ((i-1)*dx - 8.8)
c        v(nr)   total complex potential including energy dep xc
c                v is in the form  pot*r**2
c
c     Work space:
c        complex*16 p(nr), q(nr), ps(nr), qs(nr), vm(nr)
c        Must be dimensioned in calling program.  Coded like this
c        to make using different r-grids with different nrmax easy.
c
c     Output:
c        ihard   incremented each time convergence test fails
c        dny     r*g'/g, see loucks (4-85), q/p = cf/g (eq 4-86)
c        pu, qu  upper and lower components at muffin tin
c        q and q arrays  upper and lower components (see comments)

      complex*16 v(nr), e
      dimension ri(nr)
      complex*16 dny, pu, qu
      complex*16 p(nr), q(nr), ps(nr), qs(nr), vm(nr)

      include 'const.h'
      parameter (c = clight)
      parameter (csq = c**2)

      double precision lp1, ldcsq
      complex*16 c1,c2,c3,pc,qc,dp1,dq1,dp2,dq2,dp3,dq3,dp4,dq4
      complex*16 vh,vmh,vmnp1,psn,qsn,psnm1,qsnm1,psnm2,qsnm2
      complex*16 psnm3,qsnm3,psnm4,qsnm4,pp,qp,psnp1,qsnp1,prel,qrel
      complex*16 psu,vu,dummy
      complex*16 vn,vmn

c     test=1.e+04 value in loucks
      test=1.e+05
      nrk=6

      expdxh=exp(dx/2.0)
      dxd4=dx/4.0
      dxd8=dx/8.0
      a1=dx*3.30
      a2=-dx*4.20
      a3=dx*7.80
      a4=dx*14.0/45.0
      a5=dx*64.0/45.0
      a6=dx*24.0/45.0
      call diff (v,dx,jri,vm)
      twoz=-dble (v(1))/ri(1)
      l=il-1
      lp1=l+1.0
      ldcsq=l/csq
      ie=1
      r=ri(1)
      vn=v(1)
      vmn=vm(1)
cv    p(1)=1.0
      p(1)=1.e-20
      q(1)=-e/(2.0*l+3.0)*r*p(1)
      beta=lp1
      if (twoz.eq.0.0) go to 10
      beta=sqrt(lp1*l+1.0-(twoz/c)**2)
      sb0=(beta-lp1)*csq/twoz
      sa1=(3.0*beta-(twoz/c)**2)/(2.0*beta+1.0)
      sb1=csq/twoz*((beta-l)*sa1-1.0)-sb0
      sa2=((beta+3.0*lp1)*sa1-3.0*l+twoz/csq*(beta+lp1+3.0)*sb1)/
     1 (beta+1.0)/4.0
      sb2=(csq/twoz*(2.0*l*(beta+2.0-lp1)-l-(twoz/c)**2)*sa1-3.0*l
     1 *csq/twoz*(beta+2.0-lp1)+(beta+3.0-2.0*lp1-(twoz/c)**2)*sb1)/
     2 (beta+1.0)/4.0
      delta=r*csq/twoz
      q(1)=(sb0+delta*(sb1+delta*sb2))/(1.0+delta*(sa1+delta*sa2))*p(1)
   10 continue
c     runge kutta method  (see loucks)
      c1=vn/r**2-e
      c2=1.0-c1/csq
      c3=(vmn-2.0*vn)/c2/c2*ldcsq
      ps(1)=r*c2*q(1)+lp1*p(1)
      qs(1)=-lp1*q(1)+(r*c1-c3/r**3)*p(1)
      n=1
   20 continue
      pc=p(n)
      qc=q(n)
      dp1=dx*(r*c2*qc+lp1*pc)
      dq1=dx*(-lp1*qc+(r*c1-c3/r**3)*pc)
      pc=pc+0.50*dp1
      qc=qc+0.50*dq1
      r=r*expdxh
      vnp1=v(n+1)
      vmnp1=vm(n+1)
      vh=(vn+vnp1)*.50+(vmn-vmnp1)*dxd8
      vmh=(1.50*(vnp1-vn)-(vmn+vmnp1)*dxd4)/dx
      c1=vh/r/r-e
      c2=1.0-c1/csq
      c3=(vmh-2.0*vh)/c2/c2*ldcsq
      dp2=dx*(r*c2*qc+lp1*pc)
      dq2=dx*(-lp1*qc+(r*c1-c3/r**3)*pc)
      pc=pc+0.50*(dp2-dp1)
      qc=qc+0.50*(dq2-dq1)
      dp3=dx*(r*c2*qc+lp1*pc)
      dq3=dx*(-lp1*qc+(r*c1-c3/r**3)*pc)
      pc=pc+dp3-0.50*dp2
      qc=qc+dq3-0.50*dq2
      n=n+1
      r=ri(n)
      c1=vnp1/r/r-e
      c2=1.0-c1/csq
      c3=(vmnp1-2.0*vnp1)/c2/c2*ldcsq
      dp4=dx*(r*c2*qc+lp1*pc)
      dq4=dx*(-lp1*qc+(r*c1-c3/r**3)*pc)
      p(n)=p(n-1)+(dp1+2.0*(dp2+dp3)+dp4)/6.0
      q(n)=q(n-1)+(dq1+2.0*(dq2+dq3)+dq4)/6.0
      ps(n)=r*c2*q(n)+lp1*p(n)
      qs(n)=-lp1*q(n)+(r*c1-c3/r**3)*p(n)
      vn=vnp1
      vmn=vmnp1
      if (n-nrk) 20,30,30
   30 if (n.ge.jri) go to 120
      psn=ps(nrk)
      qsn=qs(nrk)
      psnm1=ps(nrk-1)
      qsnm1=qs(nrk-1)
      psnm2=ps(nrk-2)
      qsnm2=qs(nrk-2)
      psnm3=ps(nrk-3)
      qsnm3=qs(nrk-3)
      psnm4=ps(nrk-4)
      qsnm4=qs(nrk-4)
c     milne method
   40 r=ri(n+1)
      c1=v(n+1)/r/r-e
      c2=1.0-c1/csq
      c3=(vm(n+1)-2.0*v(n+1))/c2/c2*ldcsq
      pp=p(n-5)+a1*(psn+psnm4)+a2*(psnm1+psnm3)+a3*psnm2
      qp=q(n-5)+a1*(qsn+qsnm4)+a2*(qsnm1+qsnm3)+a3*qsnm2
      nit=0
   50 psnp1=r*c2*qp+lp1*pp
      qsnp1=-lp1*qp+(r*c1-c3/r**3)*pp
      pc=p(n-3)+a4*(psnp1+psnm3)+a5*(psn+psnm2)+a6*psnm1
      qc=q(n-3)+a4*(qsnp1+qsnm3)+a5*(qsn+qsnm2)+a6*qsnm1
      if (abs(test*(pc-pp))-abs(pc)) 60,60,70
   60 if (abs(test*(qc-qp))-abs(qc)) 110,110,70
   70 if (nit-40) 100,80,100
c  70 if (nit-5) 100,80,100 value in loucks
   80 prel=(pc-pp)/pc
      qrel=(qc-qp)/qc
c     count times hard test fails
      ihard = ihard + 1
c     print90, il,ie,n,prel,qrel
   90 format (' hard test in fovrg il=',i2,' ie=',i1,' n=',i3,' prel='
     1 ,e16.8,' qrel=',e16.8,' **********')
      go to 110
  100 nit=nit+1
      pp=pc
      qp=qc
      go to 50
  110 n=n+1
      p(n)=pc
      q(n)=qc
      ps(n)=psnp1
      qs(n)=qsnp1
      psnm4=psnm3
      psnm3=psnm2
      psnm2=psnm1
      psnm1=psn
      psn=psnp1
      qsnm4=qsnm3
      qsnm3=qsnm2
      qsnm2=qsnm1
      qsnm1=qsn
      qsn=qsnp1
c     introduce scale factor to prevent overflow on vax jjr
      if(abs(pc).lt.1.e+20) go to 119
      scale=1.e-20
      do 112 mm=1,6
      nm=n-mm+1
      p(nm)=scale*p(nm)
      q(nm)=scale*q(nm)
      ps(nm)=scale*ps(nm)
      qs(nm)=scale*qs(nm)
  112 continue
      psnm4=scale*psnm4
      psnm3=scale*psnm3
      psnm2=scale*psnm2
      psnm1=scale*psnm1
      psn=scale*psn
      qsnm4=scale*qsnm4
      qsnm3=scale*qsnm3
      qsnm2=scale*qsnm2
      qsnm1=scale*qsnm1
      qsn=scale*qsn
  119 if (n-jri) 40,120,120
  120 jm=jri-1
      x=dx*(xmt-jm)
      call intpol (zero,dx,p(jm),p(jri),ps(jm),ps(jri),x,pu,psu)
      call intpol (zero,dx,q(jm),q(jri),qs(jm),qs(jri),x,qu,dummy)
      call intpol (zero,dx,v(jm),v(jri),vm(jm),vm(jri),x,vu,dummy)
      dny=rmt*(1.0-(vu/rmt**2-e)/csq)*qu/pu+l
c dny is r*g'/g, see loucks (4-85), q/p = cf/g (eq 4-86)
c (watch for factors of rmt)
      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
      subroutine ZZmassivebox(j1,j2,j3,j4,j5,j6,za,zb,mt,box,drat)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'ggZZintegrals.f'
      include 'ZZdlabels.f'
      integer j1,j2,j3,j4,j5,j6,h1,h2,h3,h4,j
      double precision mt
      double complex Xpp(2,2),Xmp(2,2),Xpm(2,2),Xmm(2,2),
     & box(2,2,2,2,-2:0),d(2,2,2,2,3),
     & Xrat(2,2,2,2),drat(2,2,2,2,3)

      call ZZbox2LL(j1,j2,j3,j4,j5,j6,za,zb,mt,Xpp,Xmp,Xpm,Xmm,Xrat)
      d(2,2,:,:,d2_1_34)=Xpp(:,:)
      d(1,2,:,:,d2_1_34)=Xmp(:,:)
      d(2,1,:,:,d2_1_34)=Xpm(:,:)
      d(1,1,:,:,d2_1_34)=Xmm(:,:)
      drat(:,:,:,:,d2_1_34)=Xrat(:,:,:,:)
c      write(6,*) 'Xmp(1,1)',Xmp(1,1)
c      write(6,*) 'Xmp(1,2)',Xmp(1,2)
c      write(6,*) 'Xmp(2,1)',Xmp(2,1)
c      write(6,*) 'Xmp(2,2)',Xmp(2,2)
c      write(6,*) 'Xrat(1,1,1,1)',Xrat(1,1,1,1)

      call ZZbox2LL(j1,j2,j5,j6,j3,j4,za,zb,mt,Xpp,Xmp,Xpm,Xmm,Xrat)
      d(2,2,:,:,d1_2_34)=Xpp(:,:)
      d(1,2,:,:,d1_2_34)=Xmp(:,:)
      d(2,1,:,:,d1_2_34)=Xpm(:,:)
      d(1,1,:,:,d1_2_34)=Xmm(:,:)
      drat(:,:,:,:,d1_2_34)=Xrat(:,:,:,:)
c      write(6,*) 'Xmp(1,1)',Xmp(1,1)
c      write(6,*) 'Xmp(1,2)',Xmp(1,2)
c      write(6,*) 'Xmp(2,1)',Xmp(2,1)
c      write(6,*) 'Xmp(2,2)',Xmp(2,2)
c      write(6,*) 'Xrat(1,1,1,1)',Xrat(1,1,1,1)

      call ZZbox1LL(j1,j2,j3,j4,j5,j6,za,zb,mt,Xpp,Xmp,Xpm,Xmm,Xrat)
      d(2,2,:,:,d1_34_2)=Xpp(:,:)
      d(1,2,:,:,d1_34_2)=Xmp(:,:)
      d(2,1,:,:,d1_34_2)=Xpm(:,:)
      d(1,1,:,:,d1_34_2)=Xmm(:,:)
      drat(:,:,:,:,d1_34_2)=Xrat(:,:,:,:)
c      write(6,*) 'Xpp(1,1)',Xpp(1,1)
c      write(6,*) 'Xpp(1,2)',Xpp(1,2)
c      write(6,*) 'Xpp(2,1)',Xpp(2,1)
c      write(6,*) 'Xpp(2,2)',Xpp(2,2)
c      write(6,*) 'Xrat(1,1,1,1)',Xrat(1,1,1,1)

c--- These integrals are now initialized in ZZintegraleval
c      do e=-2,0
c      Dint(1,e)
c     & =qlI4(s34,0d0,0d0,s56,s134,s12,mtsq,mtsq,mtsq,mtsq,musq,e)
c      enddo
c      do e=-2,0
c      Dint(2,e)
c     & =qlI4(s56,0d0,0d0,s34,s156,s12,mtsq,mtsq,mtsq,mtsq,musq,e)
c      enddo
c      do e=-2,0
c      Dint(3,e)
c     & =qlI4(s56,0d0,s34,0d0,s156,s134,mtsq,mtsq,mtsq,mtsq,musq,e)
c      enddo

c      write(6,*) 'check',D0(d2_1_34)/Dint(1,0)
c      write(6,*) 'check',D0(d1_2_34)/Dint(2,0)
c      write(6,*) 'check',D0(d1_34_2)/Dint(3,0)
c      pause

c--- DEBUG
c      write(6,*) 'd(1,1,1,1,d1_2_34)/s12/s234,D0(d1_2_34)*s12*s234',
c     & d(1,1,1,1,d1_2_34)/s12/s156,D0(d1_2_34)*s12*s156
c      write(6,*) 'd(1,1,1,1,d2_1_34)/s12/s134,D0(d2_1_34)*s12*s134',
c     & d(1,1,1,1,d2_1_34)/s12/s134,D0(d2_1_34)*s12*s134
c--- DEBUG

      box(:,:,:,:,:)=czip

c--- only compute finite part here
      do h1=1,2
      do h2=1,2
      do h3=1,2
      do h4=1,2
      do j=1,3
      box(h1,h2,h3,h4,0)=box(h1,h2,h3,h4,0)+d(h1,h2,h3,h4,j)*D0(j)
      enddo
      enddo
      enddo
      enddo
      enddo

c      write(6,*) 'box(1,1,1,1,0)',box(1,1,1,1,0)

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C    %W% %G%
      subroutine savfil
c
c     CREATES A NEW SAVED DATA FILE (FOR015) FROM DATA ENTERED
c     ON THE INPUT FILE (FOR005).  CALLED BY SDATA.
c
      include 'tspinc/params.inc'
      include 'tspinc/blkcom1.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/reread.inc'
      include 'tspinc/cntrl2.inc'
      include 'tspinc/in1n.inc'
      include 'tspinc/in1c.inc'
      include 'tspinc/lnet.inc'
      include 'tspinc/param.inc'
      include 'tspinc/contrl.inc'
      include 'tspinc/comn34.inc'
      include 'tspinc/pointr.inc'
      include 'tspinc/titles.inc'

      character*10 machin(8,100),lsrcp(8,100),
     &             lrelay(8,100),lrrod(8,100),
     &             lrepdt(8,100),
     &             lshed(8,100),dccard(8,150)
      character*8 plotd
      logical debug                                                     !dem

      data plotd /'VER01'/
c
c     Begin     Begin     Begin     Begin     Begin     Begin

      call mpost('SAVFIL')
      debug = .true.                                                    !dem
      blnk10 = '          '

      if (debug)                                                        !dem
     &  call dbgeko ('SAVFIL - at start - fetching file 1 rec nums')    !dem

      call whrec1 ('MAC', ilo, ihi, isz)                                !dem
      if (debug)                                                        !dem
     &  call dbgwri ('  for gens, count rec # = ',ihi)                  !dem
      read (l1,rec=ihi) mac                                              !dem
c     read (l1,rec=101) MAC
      call whrec1 ('LDA', ilo, ihi, isz)                                !dem
      if (debug)                                                        !dem
     &  call dbgwri ('  for area load rep, count rec # = ',ihi)         !dem
      read (l1,rec=ihi) ldar                                             !dem
c     read (l1,rec=139) LDAR
      call whrec1 ('LDZ', ilo, ihi, isz)                                !dem
      if (debug)                                                        !dem
     &  call dbgwri ('  for zone load rep, count rec # = ',ihi)         !dem
      read (l1,rec=ihi) ldz                                              !dem
c     read (l1,rec=141) LDZ
      call whrec1 ('LDB', ilo, ihi, isz)                                !dem
      if (debug)                                                        !dem
     &  call dbgwri ('  for bus  load rep, count rec # = ',ihi)         !dem
      read (l1,rec=ihi) lrep                                             !dem
c     read (l1,rec=145) LREP
      call whrec1 ('LSH', ilo, ihi, isz)                                !dem
      if (debug)                                                        !dem
     &  call dbgwri ('  for load shedding, count rec # = ',ihi)         !dem
      read (l1,rec=ihi) ls                                               !dem
c     read (l1,rec=170) LS
      call whrec1 ('RG ', ilo, ihi, isz)                                !dem
      read (l1,rec=ihi) lsc                                              !dem
c     read (l1,rec=177) LSC
      call whrec1 ('RD ', ilo, ihi, isz)                                !dem
      read (l1,rec=ihi) lrd                                              !dem
c     read (l1,rec=180) LRD
      call whrec1 ('RR ', ilo, ihi, isz)                                !dem
      read (l1,rec=ihi) lrr                                              !dem
c     read (l1,rec=184) LRR
      call whrec1 ('DC ', ilo, ihi, isz)                                !dem
      read (l1,rec=ihi) jdctot,jdc                                       !dem
c     read (l1,rec=187) JDCTOT,JDC
      call opents(4)

c     ! Write header
      write (l15) swdata,savout,swdate,swdate

c     ! Write number of base kV's
      write (l15) ibxyz,(basekv(i),i=1,ibxyz)

c     ! Write LOCAL RELAY header and number of records
      write (l15) relay,blnk10,datx,lrd
      if (lrd .gt. 0) then
         call whrec1 ('RD ', ilo, ihi, isz)                                !dem
c        MSLRD=181
         mslrd = ihi + 1                                                   !dem
         if (debug)                                                        !dem
     &     call dbgwri ('  for local relays, data rec # = ',mslrd)         !dem
         numrec = lrd
         do while (numrec .gt. 0)
            num = min (100, numrec)
            read (l1,rec=mslrd) ((lrelay(i,j),i=1,8),j=1,num)
            write (l15) ((lrelay(i,j),i=1,8),j=1,num)
            mslrd = mslrd + 1
            numrec = numrec - num
         enddo
      endif

c     ! Write REMOTE RELAY header and number of records
      write (l15) remote,relay,datx,lrr
      if (lrr .gt. 0) then
         call whrec1 ('RR ', ilo, ihi, isz)                                !dem
         mslrr = ihi + 1                                                   !dem
c        MSLRR=185
         if (debug)                                                        !dem
     &     call dbgwri ('  for remote relays, data rec # = ',mslrr)        !dem
         read (l1,rec=mslrr) lrrod
         write (l15) ((lrrod(i,j),i=1,8),j=1,lrr)
      endif

c     ! Write SERIES CAPACITOR header and number of records
      write (l15) series,capac,datx,lsc
      if (lsc .gt. 0) then
         call whrec1 ('RG ', ilo, ihi, isz)                                !dem
         mslsc = ihi + 1                                                   !dem
c        MSLSC=178
         if (debug)                                                        !dem
     &     call dbgwri ('  for series cap relays, data rec # = ',mslsc)         !dem
         read (l1,rec=mslsc) lsrcp
         write (l15) ((lsrcp(i,j),i=1,8),j=1,lsc)
      endif

c     ! Write LOAD REPRESENTATION BY AREA header and number of records
      write (l15) load,repr,datx,ldar
      if (ldar .gt. 0) then
         call whrec1 ('LDA', ilo, ihi, isz)                                !dem
         mslrep = ihi + 1                                                  !dem
c        MSLREP=140
         if (debug)                                                        !dem
     &     call dbgwri ('  for area load rep, data rec # = ',mslrep)       !dem
         read (l1,rec=mslrep) (( lrepdt(i,j),i=1,8),j=1,ldar)
         write (l15) ((lrepdt(i,j),i=1,8),j=1,ldar)
      endif

c     ! Write LOAD REPRESENTATION BY ZONE header and number of records
      write (l15) load,repr,datx,ldz
      if (ldz .gt. 0) then
         call whrec1 ('LDZ', ilo, ihi, isz)                                !dem
         mslrep = ihi + 1                                                  !dem
c        MSLREP=142
         if (debug)                                                        !dem
     &      call dbgwri ('  for zone load rep, data rec # = ',mslrep)       !dem
         numrec = ldz
         do while (numrec .gt. 0)
            num = min (100, numrec)
            read (l1,rec=mslrep) ((lrepdt(i,j),i=1,8),j=1,num)
            write (l15) ((lrepdt(i,j),i=1,8),j=1,num)
            mslrep = mslrep + 1
            numrec = numrec - num
         enddo
      endif

c     ! Write LOAD REPRESENTATION BY BUS header and number of records
      write (l15) load,repr,datx,lrep
      if (lrep .gt. 0) then
         call whrec1 ('LDB', ilo, ihi, isz)                                !dem
         mslrep = ihi + 1                                                  !dem
c        MSLREP=146
         if (debug)                                                        !dem
     &     call dbgwri ('  for bus load rep, data rec # = ',mslrep)        !dem
         numrec = lrep
         do while (numrec .gt. 0)
            num = min (100, numrec)
            read (l1,rec=mslrep) ((lrepdt(i,j),i=1,8),j=1,num)
            write (l15) ((lrepdt(i,j),i=1,8),j=1,num)
            mslrep = mslrep + 1
            numrec = numrec - num
         enddo
      endif

c     ! Write LOAD SHEDDING header and number of records
      write (l15) load,shed,datx,ls
      if (ls .gt. 0) then
         call whrec1 ('LSH', ilo, ihi, isz)                                !dem
         msls = ihi + 1                                                    !dem
c        MSLS=171
         if (debug)                                                        !dem
     &     call dbgwri ('  for load shedding, data rec # = ',msls)         !dem
         numrec = ls
         do while (numrec .gt. 0)
            num = min (100, numrec)
            read (l1,rec=msls) ((lshed(i,j),i=1,8),j=1,num)
            write (l15) ((lshed(i,j),i=1,8),j=1,num)
            msls = msls + 1
            numrec = numrec - num
         enddo
      endif

c     ! Write LOAD NETTING header and number of records
      write (l15) load,net,datx,ln,ibxyz
      if (ln .gt. 0) then
         call whrec1 ('LN ', ilo, ihi, isz)                                !dem
         msln = ihi + 1                                                    !dem
c        MSLN=168
         if (debug)                                                        !dem
     &     call dbgwri ('  for load netting, data rec # = ',msln)          !dem
         write (l15) (lnetc(i),i=1,ln),(lnetn(i),i=1,ln),
     &               (basekv(ii),ii=1,ibxyz)
1060     write (l15) machn,blnk10,datx,mac
      endif

c     ! Write MACHINE (PLANT DATA) header and number of records
      if (mac .gt. 0) then
         call whrec1 ('MAC', ilo, ihi, isz)                                !dem
         msmac = ihi + 1                                                   !dem
c        MSMAC=102
         if (debug)                                                        !dem
     &      call dbgwri ('  for machines, data rec # = ',msmac)             !dem
         numrec = mac
         do while (numrec .gt. 0)
            num = min (100, numrec)
            read (l1,rec=msmac) ((machin(i,j),i=1,8),j=1,num)
            write (l15) ((machin(i,j),i=1,8),j=1,num)
            msmac = msmac + 1
            numrec = numrec - num
         enddo
      endif

c     ! Write DIRECT CURRENT header and number of records
      write (l15) direct,currnt,datx,jdc
      if (jdc .gt. 0) then
         call whrec1 ('DC ', ilo, ihi, isz)                                !dem
         msjdc = ihi + 1                                                   !dem
c        MSJDC=188
         if (debug)                                                        !dem
     &      call dbgwri ('  for dc convertors, data rec # = ',msjdc)        !dem
        read (l1,rec=msjdc) (( dccard(i,j),i=1,8),j=1,jdc)
        write (l15) ((dccard(i,j),i=1,8),j=1,jdc)
      endif

      write (l15) swdata,xend,blnk10,blnk10
      call closts (l15,0)

      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


text:
!  ------------------------------------------------------------------------------------------------------
!
!    Program HAMS for the diffraction and radiation of waves
!    by 3D structures.
!
!  License:
!
!    This routine is part of HAMS.
!
!    HAMS is a free software framework: you can redistribute it and/or modify it
!    under the terms of the Apache License, Version 2.0 (the "License"); you may
!    not use these subroutines except in compliance with the License. The software
!    is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND.
!
!    You should have received a copy of the Apache License, Version 2.0, along with
!    HAMS. If not, see <http://www.apache.org/licenses/LICENSE-2.0>.
!
!  Code Original Author:
!
!    Yingyi Liu
!
!  ------------------------------------------------------------------------------------------------------
MODULE PrintOutput

   USE Const_mod

   IMPLICIT NONE

      ! ..... Public Subroutines ...................................................................................................

   PUBLIC :: PrintHeading
   PUBLIC :: PrintBody_RealVal

   PUBLIC :: PrintBody_CmplxVal
   PUBLIC :: PrintBody_Exfc
   PUBLIC :: PrintEnd

CONTAINS
!   --------------------------------------------------------
!       Print the headings in HydroStar output files
!   --------------------------------------------------------

      SUBROUTINE PrintHeading(NFILE,NBETA,REFL,RAOType,MD1,MD2,H,XW,XC,WVHD)
      USE Const_mod
      IMPLICIT NONE

      INTEGER  NFILE,NBETA,MD1,MD2,II
      REAL*8   REFL,H,XW(2),XC(3),WVHD(NBETA),NWVHD(NBETA)
      CHARACTER(*)::RAOType

      WRITE(NFILE, '(A11)')  '# Project :'
      WRITE(NFILE, '(A11)')  '# User    :'
      IF (RAOType.EQ.'AddedMass') THEN
	   WRITE(NFILE, '(A8,A10,A1,I1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,MD2,'.rao'
      ELSEIF (adjustl(trim(RAOType)).EQ.'WaveDamping') THEN
	   WRITE(NFILE, '(A8,A12,A1,I1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,MD2,'.rao'
      ELSEIF (RAOType.EQ.'Cm') THEN
	   WRITE(NFILE, '(A8,A3,A1,I1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,MD2,'.dat'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Ca') THEN
	   WRITE(NFILE, '(A8,A3,A1,I1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,MD2,'.dat'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Excitation') THEN
	   WRITE(NFILE, '(A8,A11,A1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,'.rao'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Motion') THEN
	   WRITE(NFILE, '(A8,A7,A1,I1,A4)')  '# File :',adjustl(trim(RAOType)),'_',MD1,'.rao'
      ENDIF
	  WRITE(NFILE, '(A1)')  '#'
	  WRITE(NFILE, '(A34)')  '# Constants used in computations :'
	  WRITE(NFILE, '(A28,F11.4)')  '#     Reference length     :', REFL
	  WRITE(NFILE, '(A28,F11.4)')  '#     Water density (rho)  :', RHO
	  WRITE(NFILE, '(A28,F11.4)')  '#     Gravity acceleration :', G
      IF (H.LT.0.D0) THEN
	   WRITE(NFILE, '(A35)')  '#     Waterdepth           :  Inf. '
      ELSE
	   WRITE(NFILE, '(A28,F11.4)')  '#     Waterdepth           :',H
      ENDIF
	  WRITE(NFILE, '(A28,A2,2F9.4,A1)')  '#     Ref.pt incident wave :',' (',XW(1),XW(2),')'
	  WRITE(NFILE, '(A43)')  '#            Forward speed :   0.0000  m/s '
	  WRITE(NFILE, '(A1)')  '#'
	  WRITE(NFILE, '(A28,A2,3F9.4,A1)')  '# Reference point of body 1:' ,' (',XC(1),XC(2),XC(3),')'
	  WRITE(NFILE, '(A25)')  '# MEANVALUE :   0.0000E+00'
	  WRITE(NFILE, '(A13)')  '#   AMP/PHASE'
	  WRITE(NFILE, '(A72)')  '#------------------------------------------------------------------------'
      IF (RAOType.EQ.'AddedMass') THEN
	   WRITE(NFILE, '(A15,A9)')  '#RAOTYPE    :  ',RAOType
      ELSEIF (adjustl(trim(RAOType)).EQ.'WaveDamping') THEN
	   WRITE(NFILE, '(A15,A11)')  '#RAOTYPE    :  ',RAOType
      ELSEIF (adjustl(trim(RAOType)).EQ.'Cm') THEN
	   WRITE(NFILE, '(A15,A9)')  '#RAOTYPE    :  ',RAOType
      ELSEIF (adjustl(trim(RAOType)).EQ.'Ca') THEN
	   WRITE(NFILE, '(A15,A11)')  '#RAOTYPE    :  ',RAOType
      ELSEIF (adjustl(trim(RAOType)).EQ.'Excitation') THEN
	   WRITE(NFILE, '(A15,A10)')  '#RAOTYPE    :  ',RAOType
      ELSEIF (adjustl(trim(RAOType)).EQ.'Motion') THEN
	   WRITE(NFILE, '(A15,A6)')  '#RAOTYPE    :  ',RAOType
      ENDIF
	  WRITE(NFILE, '(A15,I1,I1)')  '#COMPONENT  :  ',MD1,MD2
      IF (adjustl(trim(RAOType)).EQ.'AddedMass') THEN
	   WRITE(NFILE, '(A15,A2)')  '#UNIT       :  ','kg'
      ELSEIF (adjustl(trim(RAOType)).EQ.'WaveDamping') THEN
	   WRITE(NFILE, '(A15,A4)')  '#UNIT       :  ','kg/s'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Cm') THEN
	   WRITE(NFILE, '(A15,A2)')  '#UNIT       :  ','kg'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Ca') THEN
	   WRITE(NFILE, '(A15,A4)')  '#UNIT       :  ','kg/s'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Excitation') THEN
	   WRITE(NFILE, '(A15,A9)')  '#UNIT       :  ','N/m, Nm/m'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Motion'.AND.MD1.LE.3) THEN
	   WRITE(NFILE, '(A15,A3)')  '#UNIT       :  ','m/m'
      ELSEIF (adjustl(trim(RAOType)).EQ.'Motion'.AND.MD1.LE.6) THEN
	   WRITE(NFILE, '(A15,A5)')  '#UNIT       :  ','deg/m'
      ENDIF

      WRITE(NFILE, '(A10,I6)')  '#NBHEADING',NBETA
      IF (RAOType.EQ.'Motion') THEN
       DO II=1,NBETA
        IF (WVHD(II).LT.0.D0) THEN
         NWVHD(II)=WVHD(II)+360.D0
        ELSE
         NWVHD(II)=WVHD(II)
        ENDIF
       ENDDO
	   WRITE(NFILE, '(A8,<NBETA>(7X,F7.2))')  '#HEADING  ',(NWVHD(II),II=1,NBETA)
      ELSE
	   WRITE(NFILE, '(A8,<NBETA>(7X,F7.2))')  '#HEADING  ',(WVHD(II),II=1,NBETA)
      ENDIF
	  WRITE(NFILE, '(A63)')  '#---w(r/s)-----------------------------------------------------'

      RETURN
      END SUBROUTINE PrintHeading

!   --------------------------------------------------------
!       Print the body in HydroStar output files
!   --------------------------------------------------------

      SUBROUTINE PrintBody_RealVal(NFILE,W1,NBETA,RAOType,VAB)
      USE Const_mod
      IMPLICIT NONE

      INTEGER  NFILE,NBETA,II
      REAL*8   W1,VAB,REL(NBETA),IMG(NBETA)
      CHARACTER(*)::RAOType
      CHARACTER(LEN=100)::FMT

      DO II=1,NBETA
       REL(II)=VAB
       IMG(II)=0.D0
      ENDDO

      WRITE(FMT,*) '(F8.4,',NBETA,'(ES14.6),',NBETA,'(F12.4))'
      WRITE(NFILE,FMT) W1,(REL(II),II=1,NBETA),(IMG(II),II=1,NBETA)

      RETURN
      END SUBROUTINE PrintBody_RealVal

!   --------------------------------------------------------
!       Print the body in HydroStar output files
!   --------------------------------------------------------

      SUBROUTINE PrintBody_CmplxVal(NFILE,W1,NBETA,RAOType,CVAB)
      USE Const_mod
      IMPLICIT NONE

      INTEGER  NFILE,NBETA,II
      REAL*8   W1,NREL(NBETA),NIMG(NBETA),REL(NBETA),IMG(NBETA),MDL(NBETA),PHS(NBETA)
      COMPLEX*16 CVAB(NBETA)
      CHARACTER(*)::RAOType
      CHARACTER(LEN=100)::FMT

      DO II=1,NBETA
        REL(II)=DREAL(CVAB(II))
        IMG(II)=DIMAG(CVAB(II))
        MDL(II)=CDABS(CVAB(II))
        NREL(II)=-IMG(II)
        NIMG(II)=-REL(II)
        PHS(II)=ATAN2D(NIMG(II),NREL(II))
        IF (PHS(II).LT.0.D0) PHS(II)=PHS(II)+360.D0
      ENDDO

      WRITE(FMT,*) '(F8.4,',NBETA,'(ES14.6),',NBETA,'(F12.4))'
      WRITE(NFILE,FMT) W1,(MDL(II),II=1,NBETA),(PHS(II),II=1,NBETA)
!	  WRITE(NFILE,200) W1,(NREL(II),II=1,NBETA),(NIMG(II),II=1,NBETA)

100	  FORMAT(F8.4,<NBETA>(ES14.6),<NBETA>(F12.4))
200	  FORMAT(F8.4,<NBETA>(ES14.6),<NBETA>(ES14.6))
      RETURN
      END SUBROUTINE PrintBody_CmplxVal

!   --------------------------------------------------------
!       Print the body in HydroStar output files
!   --------------------------------------------------------

      SUBROUTINE PrintBody_Exfc(NFILE,W1,NBETA,RAOType,CVAB)
      USE Const_mod
      IMPLICIT NONE

      INTEGER  NFILE,NBETA,II
      REAL*8   W1,NREL(NBETA),NIMG(NBETA),REL(NBETA),IMG(NBETA),MDL(NBETA),PHS(NBETA)
      COMPLEX*16 CVAB(NBETA)
      CHARACTER(*)::RAOType
      CHARACTER(LEN=100)::FMT

      DO II=1,NBETA
        REL(II)=DREAL(CVAB(II))
        IMG(II)=DIMAG(CVAB(II))
        MDL(II)=CDABS(CVAB(II))
        PHS(II)=ATAN2D(REL(II),IMG(II))
      ENDDO

      WRITE(FMT,*) '(F8.4,',NBETA,'(ES14.6),',NBETA,'(F12.4))'
      WRITE(NFILE,FMT) W1,(MDL(II),II=1,NBETA),(PHS(II),II=1,NBETA)

100	  FORMAT(F8.4,<NBETA>(ES14.6),<NBETA>(F12.4))
200	  FORMAT(F8.4,<NBETA>(ES14.6),<NBETA>(ES14.6))
      RETURN
      END SUBROUTINE PrintBody_Exfc

!   --------------------------------------------------------
!       Print the end in HydroStar output files
!   --------------------------------------------------------

      SUBROUTINE PrintEnd(NFILE)
      IMPLICIT NONE

      INTEGER NFILE

      WRITE(NFILE, '(A61)')  '#------------------------------------------------------------'
      WRITE(NFILE, '(A8)')  '#ENDFILE'

      RETURN
      END SUBROUTINE PrintEnd
!-------------------------------------------------------------------------------
END MODULE PrintOutput
!*******************************************************************************

max_stars_repo_name: Greentp/HAMS

max_issues_repo_name: Greentp/HAMS

max_forks_repo_name: Greentp/HAMS


text:
      SUBROUTINE zfver (CFNAME, CVER, IVER)
C
C
C     Determine a DSS file's version number
C     CVER is returned with the 4 character version (6-FD)
C     IVER is returned with:
C        -3:  Not a DSS file
C        -2:  Unable to access file (but it exists)
C        -1:  File does not exist
C         4:  Version 4 file
C         5:  Version 5 file
C         6:  Version 6 file
C
      INTEGER*4 IREC(5)
      CHARACTER CFNAME*(*), CVER*(*), CNAME*64
      LOGICAL LEXIST, LOPEN
C
C
      CVER = ' '
C
      CALL zfname6 (CFNAME, CNAME, ILEN, LEXIST)
      IF (.NOT.LEXIST) GO TO 910
C
C     Try to find a unit to attach to the DSS file
      JUNIT = 80
      DO 20 I=1,10
      INQUIRE (UNIT=JUNIT, OPENED=LOPEN)
      IF (.NOT.LOPEN) GO TO 40
      JUNIT = JUNIT + 1
 20   CONTINUE
C
      JUNIT = 70
      INQUIRE (UNIT=JUNIT, OPENED=LOPEN)
      IF (.NOT.LOPEN) GO TO 40
      GO TO 920
C
 40   CONTINUE
C
      OPEN (UNIT=JUNIT, FILE=CNAME(1:ILEN), ACCESS='DIRECT', RECL=20,
     * IOSTAT=ISTAT)
      IF (ISTAT.NE.0) GO TO 920
*     READ (UNIT=JUNIT, REC=1, END=930, ERR=930) IREC
      READ (UNIT=JUNIT, REC=1, ERR=930) IREC
      CLOSE (UNIT=JUNIT)
C
      CNAME = ' '
      CALL HOLCHR (IREC(5), 1, 4, CNAME, 1)
C
      IF (CNAME(2:2).NE.'-') THEN
          IF (CNAME(3:3).EQ.'-') THEN
             CNAME(1:1) = CNAME(4:4)
          ELSE
             GO TO 930
          ENDIF
      ENDIF
C
      CVER = CNAME
C
      IF (CNAME(1:1).EQ.'4') THEN
      IVER = 4
      ELSE IF (CNAME(1:1).EQ.'5') THEN
      IVER = 5
      ELSE IF (CNAME(1:1).EQ.'6') THEN
      IVER = 6
      ELSE IF (CNAME(1:1).EQ.'7') THEN
      IVER = 7
      ELSE
      GO TO 930
      ENDIF
C
C
 800  CONTINUE
      RETURN
C
 910  CONTINUE
      IVER = -1
      GO TO 800
C
 920  CONTINUE
      IVER = -2
      GO TO 800
C
 930  CONTINUE
      IVER = -3
      GO TO 800
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib