added first version of spectral method, changed corrected functions to get correct path to files

code/mpie_spectral.f90

@@ -92,22 +92,24 @@ function getSolverJobName()
  posExt = scan(outName,'.',back=.true.)
  posSep = scan(outName,pathSep,back=.true.)
  
- if (posExt <= posSep) posExt = len_trim(outName)       ! no extension present
+ if (posExt <= posSep) posExt = len_trim(outName)+1       ! no extension present
-
+ getSolverJobName = outName(1:posExt-1)                   ! path to mesh file (excl. extension)
- getSolverJobName = outName(posSep+1:posExt-1)          ! path to mesh file (excl. extension)
  
  if (scan(getSolverJobName,pathSep) /= 1) then            ! relative path given as command line argument
    call getcwd(cwd)
-   getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),&
+   getSolverJobName = rectifyPath(trim(cwd)//'/'//getSolverJobName)
-                                       rectifyPath(trim(cwd)//'/'//getSolverJobName))
+ else
+   getSolverJobName = rectifyPath(getSolverJobName)
  endif
  
+ getSolverJobName = makeRelativePath(getSolverWorkingDirectoryName(),&
+                                    getSolverJobName)
  return
 end function
 
 
 !********************************************************************
-! realive path of loadcase from command line arguments
+! relative path of loadcase from command line arguments
 !
 !********************************************************************
 function getLoadcaseName()
@@ -119,17 +121,22 @@ function getLoadcaseName()
  character(len=1024) getLoadcaseName, outName, cwd
  character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \
  integer(pInt) posExt,posSep
+ posExt=0                                                 !not sure if its needed
 
  call getarg(2,getLoadcaseName)
+ posExt = scan(getLoadcaseName,'.',back=.true.)
+ posSep = scan(getLoadcaseName,pathSep,back=.true.)
  
+ if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load')   ! no extension present
  if (scan(getLoadcaseName,pathSep) /= 1) then          ! relative path given as command line argument
    call getcwd(cwd)
    getLoadcaseName = rectifyPath(trim(cwd)//'/'//getLoadcaseName)
+ else
+   getLoadcaseName = rectifyPath(getLoadcaseName)
  endif
  
  getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
                                     getLoadcaseName)
-  
  return
 end function
 
@@ -167,6 +174,7 @@ function rectifyPath(path)
     rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j)
     i = j+index(rectifyPath(j+1:l),'../')
  end do
+ if(len_trim(rectifyPath)==0) rectifyPath='/'
  return
  
  end function rectifyPath
@@ -188,7 +196,6 @@ function makeRelativePath(a,b)
 
  posLastCommonSlash = 0
  remainingSlashes = 0
-
  do i = 1,min(1024,len_trim(a),len_trim(b))
    if (a(i:i) /= b(i:i)) exit
    if (a(i:i) == '/') posLastCommonSlash = i
@@ -197,7 +204,6 @@ function makeRelativePath(a,b)
    if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1
  enddo
  makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
- 
  return
 end function makeRelativePath
 
@@ -233,7 +239,8 @@ program mpie_spectral
  use mpie_interface
  use prec, only: pInt, pReal
  use IO
- use math, only: math_I3,math_transpose3x3
+! use math, only: math_I3,math_transpose3x3,math_Mandel66to3333
+  use math
  use CPFEM, only: CPFEM_general
 
  implicit none
@@ -262,6 +269,43 @@ program mpie_spectral
  integer(pInt) unit, N_l, N_s, N_t, N_n, N, i, j, k, l            ! numbers of identifiers, loop variables
  integer(pInt) a, b, c, e, homog
  real(pReal) x, y, z
+!-------------------------
+!begin RL
+!-------------------------
+
+ ! begin change by m.diehl
+ !integer(pInt), parameter :: npts1 = 32                            
+ !integer(pInt), parameter :: npts2 = 32
+ !integer(pInt), parameter :: npts3 = 32
+ 
+ integer(pInt) npts1, npts2, npts3 ! will be allocated later
+ ! end change by m.diehl
+
+ real(pReal), dimension(:), allocatable :: datafft
+ real(pReal), dimension(:,:,:,:,:), allocatable :: workfft,workfftim,sg,disgrad,defgradold
+
+ integer(pInt), dimension (3,3) :: iudot,iscau
+  
+ real(pReal), dimension(3,3) :: disgradmacro,disgradmacroactual
+ real(pReal), dimension(3,3) :: ddisgradmacro,ddisgradmacroacum,ddisgrad,ddisgradim
+ real(pReal), dimension(3,3) :: defgrad0,defgrad
+ real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33
+
+ real(pReal), dimension(3) :: delt,nn,xk,xk2
+ real(pReal), dimension(2) :: nn2
+ real(pReal), dimension(6) :: aux6
+ real(pReal), dimension(3,3,3,3) :: c0,s0,g1
+ real(pReal), dimension(6,6) :: c066,s066
+
+ integer(pInt) itmax,nsteps,jload,ielem,ii,jj,k1,kxx,kyy,kzz,kx,ky,kz,idum,iter,imicro,m1,n1,p,q
+
+ real(pReal) prodnn,wgt,error,tdot,erre,errs,evm,svm,det,xknorm
+
+ logical errmatinv
+
+!-------------------------
+!end RL
+!-------------------------
  
  if (IargC() < 2) call IO_error(102)
 
@@ -273,6 +317,9 @@ program mpie_spectral
  N_t = 0_pInt
  N_n = 0_pInt
  
+ print*,'Loadcase: ',trim(path)
+ print*,'Workingdir: ',trim(getSolverWorkingDirectoryName())
+
  if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg=path)
 
  rewind(unit)
@@ -366,7 +413,8 @@ program mpie_spectral
  gotDimension =       .false.
  gotHomogenization =  .false.
  path = getSolverJobName()
- if (.not. IO_open_file(unit,path)) call IO_error(101,ext_msg=path)
+ print*,'JobName: ',trim(path)
+ if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg=path)
 
  rewind(unit)
  do 
@@ -405,9 +453,491 @@ program mpie_spectral
    end select
    if (gotDimension .and. gotHomogenization .and. gotResolution) exit
  enddo
+ 100 close(unit)
+ 
+   print '(a,/,i3,i3,i3)','resolution a b c',a,b,c
+   print '(a,/,f6.2,f6.2,f6.2)','dimension x y z',x, y, z
+   print *,'homogenization',homog
+   print *, ''
 
  temperature = 300.0_pReal
- call CPFEM_general(2,math_i3,math_i3,temperature,0.0_pReal,1_pInt,1_pInt,stress,dsde)
+ !-------------------------
+!RL
+!-------------------------
+!begin change m.diehl
+ npts1=a
+ npts2=b
+ npts3=c
+!end change m.diehl
+
+ allocate (datafft(2*npts1*npts2*npts3))
+ allocate (workfft(3,3,npts1,npts2,npts3))
+ allocate (workfftim(3,3,npts1,npts2,npts3))
+ allocate (sg(3,3,npts1,npts2,npts3))
+ allocate (disgrad(3,3,npts1,npts2,npts3))
+ allocate (defgradold(3,3,npts1,npts2,npts3))
+
+ error=0.00001
+ itmax=100
+
+ delt(1)=1.
+ delt(2)=1.
+ delt(3)=1.
+
+ nn(1)=npts1
+ nn(2)=npts2
+ nn(3)=npts3
+
+ nn2(1)=npts1
+ nn2(2)=npts2
+
+ prodnn=nn(1)*nn(2)*nn(3)
+ wgt=1./prodnn
+
+! C_0 and S_0 CALCULATION
+
+!!
+!! PHILIP: FE_exec_elem?
+!!
+
+ ielem=0
+ c0=0.
+ c066=0.
+
+ do k=1,npts3
+ do j=1,npts2
+ do i=1,npts1
+
+  ielem=ielem+1
+
+  call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
+
+  c066=c066+dsde   
+  c0=c0+math_Mandel66to3333(dsde)
+
+ enddo
+ enddo
+ enddo
+
+  c066=c066*wgt   
+  c0=c0*wgt
+  
+ call math_invert(6,c066,s066,idum,errmatinv)
+
+ if(errmatinv) then
+  write(*,*) 'ERROR IN C0 INVERSION'
+  stop
+ endif
+
+ s0=math_Mandel66to3333(s066)
+
+! INITIALIZATION BEFORE STARTING WITH LOADINGS
+
+  disgradmacro=0.
+
+  do k=1,npts3
+  do j=1,npts2
+  do i=1,npts1
+
+   defgradold(:,:,i,j,k)=math_i3(:,:)
+
+  enddo
+  enddo
+  enddo
+
+ do jload = 1,N
+
+  udot(:,:)=bc_velocityGrad(:,:,jload)
+  scauchy(:,:)=bc_stress(:,:,jload)
+
+  iudot=0
+  iscau=0
+  do i=1,3
+  do j=1,3
+   if(bc_mask(i,j,1,jload)) iudot(i,j)=1
+   if(bc_mask(i,j,2,jload)) iscau(i,j)=1
+ enddo
+ enddo
+ 
+  nsteps=bc_steps(jload)
+  tdot=bc_timeIncrement(jload)/bc_steps(jload)
+
+! evm=dvm*tdot ?
+
+  do imicro=1,nsteps
+
+   write(*,*) '***************************************************'
+   write(*,*) 'STEP = ',imicro
+!  write(21,*) 'STEP = ',imicro
+
+! INITIALIZATION BEFORE NEW TIME STEP
+
+   disgradmacro=disgradmacro+udot*tdot
+   ddisgradmacro=0.
+
+   ielem=0  
+
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    ielem=ielem+1
+
+    disgrad(:,:,i,j,k)=disgradmacro(:,:)
+
+    defgrad0(:,:)=defgradold(:,:,i,j,k)
+    defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
+
+    call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
+
+!    sg(:,:,i,j,k)=math_Mandel6to33(stress) ?
+
+   enddo
+   enddo
+   enddo
+
+   ielem=0  
+
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    ielem=ielem+1
+
+    disgrad(:,:,i,j,k)=disgradmacro(:,:)
+
+    defgrad0(:,:)=defgradold(:,:,i,j,k)
+    defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
+
+    call CPFEM_general(1,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
+
+    sg(:,:,i,j,k)=math_Mandel6to33(stress)
+
+   enddo
+   enddo
+   enddo
+  
+   ddisgradmacroacum=0.
+
+   iter=0
+   erre=2.*error
+   errs=2.*error
+
+   do while(iter.lt.itmax.and.(errs.gt.error.or.erre.gt.error))
+
+    iter=iter+1
+
+    write(*,*)'ITER = ',iter
+    write(*,*) 'DIRECT FFT OF STRESS FIELD'
+
+    do ii=1,3
+    do jj=1,3
+
+     k1=0
+
+     do k=1,npts3
+
+!      write(*,'(1H+,a,i2,2(a,i4))') 
+!     #   'STRESS - COMPONENT',ii,jj,'  -   Z = ',k,'   OUT OF ',npts
+
+      do j=1,npts2
+      do i=1,npts1
+
+      k1=k1+1
+      datafft(k1)=sg(ii,jj,i,j,k)
+      k1=k1+1
+      datafft(k1)=0.
+
+     enddo
+     enddo
+     enddo
+
+     if(npts3.gt.1) then
+      call fourn(datafft,nn,3,1)
+     else
+      call fourn(datafft,nn2,2,1)
+     endif
+
+     k1=0
+
+     do k=1,npts3
+     do j=1,npts2
+     do i=1,npts1
+
+      k1=k1+1
+      workfft(ii,jj,i,j,k)=datafft(k1)
+
+      k1=k1+1
+      workfftim(ii,jj,i,j,k)=datafft(k1)
+
+     enddo
+     enddo
+     enddo
+
+    enddo
+    enddo
+
+    write(*,*) 'CALCULATING G^pq,ij : SG^ij ...'
+
+    do kzz=1,npts3
+    do kyy=1,npts2
+    do kxx=1,npts1
+
+     if(kxx.le.npts1/2) kx=kxx-1
+     if(kxx.gt.npts1/2) kx=kxx-npts1-1
+
+     if(kyy.le.npts2/2) ky=kyy-1
+     if(kyy.gt.npts2/2) ky=kyy-npts2-1
+
+     if(kzz.le.npts3/2) kz=kzz-1
+     if(kzz.gt.npts3/2) kz=kzz-npts3-1
+
+     xk(1)=kx/(delt(1)*nn(1))
+     xk(2)=ky/(delt(2)*nn(2))
+
+     if(npts3.gt.1) then
+      xk(3)=kz/(delt(3)*nn(3))
+     else
+      xk(3)=0.
+     endif
+
+     xknorm=sqrt(xk(1)**2+xk(2)**2+xk(3)**2)
+
+     if (xknorm.ne.0.) then
+      do i=1,3
+       xk2(i)=xk(i)/(xknorm*xknorm*2.*pi)
+       xk(i)=xk(i)/xknorm
+      enddo
+     endif
+
+     do i=1,3
+     do k=1,3
+      aux33(i,k)=0.
+     do j=1,3
+     do l=1,3
+      aux33(i,k)=aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l)
+     enddo
+     enddo
+     enddo
+     enddo
+
+     aux33=math_inv3x3(aux33)
+
+!     call minv(aux33,3,det,minv1,minv2)
+!     if(det.eq.0) then
+!      write(*,*) kx,ky,kz,'  --> SINGULAR SYSTEM'
+!      stop
+!      pause
+!     endif
+
+     do p=1,3
+     do q=1,3
+     do i=1,3
+     do j=1,3
+      g1(p,q,i,j)=-aux33(p,i)*xk(q)*xk(j)
+     enddo
+     enddo
+     enddo
+     enddo
+
+     do i=1,3
+     do j=1,3
+
+      ddisgrad(i,j)=0.
+      ddisgradim(i,j)=0.
+
+!      if(kx.eq.0.and.ky.eq.0.and.kz.eq.0.) goto 4
+      if(.not.(kx.eq.0.and.ky.eq.0.and.kz.eq.0.)) then
+
+       do k=1,3
+       do l=1,3
+
+        ddisgrad(i,j)=ddisgrad(i,j)+g1(i,j,k,l)*workfft(k,l,kxx,kyy,kzz)
+        ddisgradim(i,j)=ddisgradim(i,j)+g1(i,j,k,l)*workfftim(k,l,kxx,kyy,kzz)
+
+       enddo
+       enddo
+
+      endif
+
+     enddo
+     enddo
+
+      workfft(:,:,kxx,kyy,kzz)=ddisgrad(:,:)
+      workfftim(:,:,kxx,kyy,kzz)=ddisgradim(:,:)
+
+    enddo
+    enddo
+    enddo
+
+    write(*,*) 'INVERSE FFT TO GET DISPLACEMENT GRADIENT FIELD'
+
+    do m1=1,3
+    do n1=1,3
+
+     k1=0
+
+     do k=1,npts3
+     do j=1,npts2
+     do i=1,npts1
+
+      k1=k1+1
+      datafft(k1)=workfft(m1,n1,i,j,k)
+
+      k1=k1+1
+      datafft(k1)=workfftim(m1,n1,i,j,k)
+
+     enddo
+     enddo
+     enddo
+
+     if(npts3.gt.1) then
+      call fourn(datafft,nn,3,-1)
+     else
+      call fourn(datafft,nn2,2,-1)
+     endif
+
+      datafft=datafft*wgt
+
+     k1=0
+     do k=1,npts3
+     do j=1,npts2
+     do i=1,npts1
+
+      k1=k1+1
+      disgrad(m1,n1,i,j,k)=disgrad(m1,n1,i,j,k)+ddisgradmacro(m1,n1)+datafft(k1)
+      k1=k1+1
+
+     enddo
+     enddo
+     enddo
+
+    enddo
+    enddo
+
+    write(*,*) 'UPDATE STRESS FIELD'
+
+!!!    call evpal
+
+   ielem=0  
+
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    ielem=ielem+1
+
+    defgrad0(:,:)=defgradold(:,:,i,j,k)
+    defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
+
+    call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
+
+    sg(:,:,i,j,k)=math_Mandel6to33(stress)
+
+   enddo
+   enddo
+   enddo
+
+   ielem=0  
+
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    ielem=ielem+1
+
+    defgrad0(:,:)=defgradold(:,:,i,j,k)
+    defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k)
+
+    call CPFEM_general(2,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
+
+    sg(:,:,i,j,k)=math_Mandel6to33(stress)
+
+   enddo
+   enddo
+   enddo
+
+!!!    call get_smacro
+
+!  AVERAGE STRESS
+
+   scauav=0.
+   
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    scauav(:,:)=scauav(:,:)+sg(:,:,i,j,k)*wgt
+
+   enddo
+   enddo
+   enddo
+
+!  MIXED BC
+
+   do i=1,3
+   do j=1,3
+
+    ddisgradmacro(i,j)=0.
+   
+    if(iudot(i,j).eq.0) then
+
+     do k=1,3
+     do l=1,3
+ 
+      ddisgradmacro(i,j)=ddisgradmacro(i,j)+s0(i,j,k,l)*iscau(k,l)*(scauchy(k,l)-scauav(k,l))
+
+     enddo
+     enddo
+
+    endif
+
+    enddo
+    enddo
+
+    ddisgradmacroacum=ddisgradmacroacum+ddisgradmacro
+ 
+!      write(*,*) 'DDEFGRADMACRO(1,1),(2,2)=',ddefgradmacro(1,1),ddefgradmacro(2,2)
+
+!!  svm= ?
+
+!    erre=erre/evm
+!    errs=errs/svm
+
+    write(*,*) 'STRESS FIELD ERROR =',errs
+    write(*,*) 'STRAIN FIELD ERROR =',erre
+
+!      write(21,101) iter,erre,errs,svm
+!101   format(i3,4(1x,e10.4),10(1x,F7.4))
+
+   enddo    ! WHILE  ENDDO
+
+   disgradmacroactual=disgradmacro+ddisgradmacroacum
+
+!!      write(*,*) 'defgradmacro(1,1),defgradmacro(2,2),defgradmacro(3,3)'
+!!      write(*,*) defgradmacroactual(1,1),defgradmacroactual(2,2),defgradmacroactual(3,3)
+!!      write(*,*) 'defgradmacro(1,1)/defgradmacro(3,3)'
+!!      write(*,*) defgradmacroactual(1,1)/defgradmacroactual(3,3)
+!!      write(*,*) 'scauav(1,1),scauav(2,2),scauav(3,3)'
+!!      write(*,*) scauav(1,1),scauav(2,2),scauav(3,3)
+
+   do k=1,npts3
+   do j=1,npts2
+   do i=1,npts1
+
+    defgradold(:,:,i,j,k)=math_i3(:,:)+disgrad(:,:,i,j,k)
+
+   enddo
+   enddo
+   enddo
+    
+  enddo  ! IMICRO ENDDO
+
+ enddo  ! JLOAD ENDDO
+
+!-------------------------
+!RL
+!-------------------------
 
 end program mpie_spectral
 
@@ -424,3 +954,110 @@ subroutine quit(id)
 
  stop
 end subroutine
+
+   SUBROUTINE fourn(data,nn,ndim,isign)
+!
+!     FROM NUMERICAL RECIPES IN F77 (FIXED FORMAT), 
+!     CONVERTED INTO FREE FORMAT (RL @ MPIE, JUNE 2010)
+!
+      INTEGER isign,ndim,nn(ndim)
+      REAL data(*)
+!      INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,
+!     *k2,n,nprev,nrem,ntot
+ INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,k2,n,nprev,nrem,ntot
+      REAL tempi,tempr
+      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
+      ntot=1
+!      do 11 idim=1,ndim
+ do idim=1,ndim  ! 11
+!
+        ntot=ntot*nn(idim)
+!11    continue
+ enddo ! 11
+!
+      nprev=1
+!      do 18 idim=1,ndim
+ do idim=1,ndim  ! 18
+        n=nn(idim)
+        nrem=ntot/(n*nprev)
+        ip1=2*nprev
+        ip2=ip1*n
+        ip3=ip2*nrem
+        i2rev=1
+!        do 14 i2=1,ip2,ip1
+   do i2=1,ip2,ip1  !  14
+          if(i2.lt.i2rev)then
+
+!            do 13 i1=i2,i2+ip1-2,2
+!              do 12 i3=i1,ip3,ip2
+     do i1=i2,i2+ip1-2,2  ! 13
+       do i3=i1,ip3,ip2  ! 12
+                i3rev=i2rev+i3-i2
+                tempr=data(i3)
+                tempi=data(i3+1)
+                data(i3)=data(i3rev)
+                data(i3+1)=data(i3rev+1)
+                data(i3rev)=tempr
+                data(i3rev+1)=tempi
+!12            continue
+!13          continue
+     enddo  !  12
+       enddo  !  13
+          endif
+          ibit=ip2/2
+
+!1         if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
+     do while ((ibit.ge.ip1).and.(i2rev.gt.ibit))  ! if 1
+            i2rev=i2rev-ibit
+            ibit=ibit/2
+!          goto 1
+!          endif
+     enddo  ! do while (if 1)
+
+          i2rev=i2rev+ibit
+!14      continue
+   enddo  ! 14
+        ifp1=ip1
+
+!2       if(ifp1.lt.ip2)then
+    do while (ifp1.lt.ip2)  !  if 2
+          ifp2=2*ifp1
+          theta=isign*6.28318530717959d0/(ifp2/ip1)
+          wpr=-2.d0*sin(0.5d0*theta)**2
+          wpi=sin(theta)
+          wr=1.d0
+          wi=0.d0
+!          do 17 i3=1,ifp1,ip1
+!            do 16 i1=i3,i3+ip1-2,2
+!              do 15 i2=i1,ip3,ifp2
+     do i3=1,ifp1,ip1  !  17
+       do i1=i3,i3+ip1-2,2  !  16
+         do i2=i1,ip3,ifp2  !  15
+                k1=i2
+                k2=k1+ifp1
+                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
+                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
+                data(k2)=data(k1)-tempr
+
+                data(k2+1)=data(k1+1)-tempi
+                data(k1)=data(k1)+tempr
+                data(k1+1)=data(k1+1)+tempi
+!15            continue
+!16          continue
+         enddo  !  15
+       enddo  !  16
+            wtemp=wr
+            wr=wr*wpr-wi*wpi+wr
+            wi=wi*wpr+wtemp*wpi+wi
+!17        continue
+     enddo  !  17
+          ifp1=ifp2
+!        goto 2
+!        endif
+    enddo  !  do while (if 2)
+
+        nprev=n*nprev
+!18    continue
+ enddo  ! 18
+      return
+      END subroutine