text:
Program stitch
! This code stitches the Low- and High-frequency synthetic
! ground motions together in wavelet domain.
!
! Here assumed that the inputted ground motions are velocities.
!
! Written by Pengcheng Liu
! Copyright (c) 2007 by Pengcheng Liu
! modified by J. Schmedes to read in 3comp files
!
 implicit NONE
 character(len=256):: fileinp,fileRec_list,fileVM_list,fileVM, &
                     st_name,path1D,path3D
 integer:: num_sm,id_motion,nst,ks,ism,getlen
 real:: freq_join,freq_max,depth_hyp
 real:: st_x,st_y,st_z,stan_h1,stan_h2,xx00,yy00,afa
 real:: dis_x,dis_y,dist,time0
!
 fileinp='stitchBBP.inp'
 open(unit=8,file=fileinp,status='old')

! write(*,*) 'Enter the name of file listing the stations'
 read(8,'(1a)') fileRec_list

! write(*,*) 'Enter the name of file listng 1D VM for each stations'
 read(8,'(1a)') fileVM_list

! write(*,*) 'Enter the name of folder(and Path) with 1D synthetics'
 read(8,'(1a)') path1D

! write(*,*) 'Enter the name of folder(and Path) with nonlinear synthetics'

 read(8,'(1a)') path3D

! write(*,*) 'Enter the joint frequency and fmax '
 read(8,*) freq_join, freq_max

! write(*,*) 'Enter the depth of hypocenter (km)'
 read(8,*) depth_hyp

! write(*,*) 'Enter the number of source model'
 read(8,*) num_sm

! write(*,*) 'Outptut Displacement (=1), Vel. (2), or Acc (3)'
 read(8,*) id_motion

 close(8)
!============================================================

 open(10,file=fileVM_list,status='old', position='rewind')
 open(9,file=fileRec_list,status='old', position='rewind')
 read(9,*) nst,xx00,yy00,afa
 do ks=1,nst
   read(10,'(1a)') fileVM
   call cmodel(fileVM,depth_hyp)
!
   read(9,'(1a)') st_name
   read(9,*) st_x,st_y,st_z,stan_h1,stan_h2
   dis_x=(st_x-xx00)/1000.0
   dis_y=(st_y-yy00)/1000.0
   dist=sqrt(dis_x**2+dis_y**2)
   call trav(dist,time0)
   !write(*,*) ks,time0
   !write(*,*) st_name
   do ism=1,num_sm
     call ReadVel(ism,stan_h1,stan_h2,st_name,path1D,path3D)
     call align(freq_join,time0)
     call stitch_wt(freq_join)
     call time_history(id_motion,freq_max)
     call output3(ism,stan_h1,stan_h2,st_name)
   enddo  ! end loop of source model (ism)
 enddo  ! end loop of station (ks)
 close(9)
 close(10)
 stop
end Program stitch
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
MODULE ccmodel
 implicit NONE
 save
 integer:: ns,lmax,nlayer,nb
 real:: freq_ref
 real, allocatable, dimension(:):: afa,beta,rho,thick,qp,qs,dep
 real(kind=8), allocatable, dimension(:):: th,v
END MODULE ccmodel
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine cmodel(fileVM,depth)
 use ccmodel
 implicit none
 character(len=*), intent(IN):: fileVM
 real, intent(IN):: depth
 integer:: jo,j1,j
!
 open(unit=19,file=fileVM,status='old',position='rewind')
 read(19,*) jo, freq_ref
 if(allocated(beta)) deallocate(afa,beta,rho,thick,qp,qs,dep,th,v)
 j1=jo+1
 allocate(afa(j1),beta(j1),rho(j1),thick(j1),qp(j1),qs(j1))
 allocate(th(j1),v(j1),dep(j1+1))
 do j=1,jo
   read (19,*) afa(j),beta(j),rho(j),thick(j),qp(j),qs(j)
 enddo
 close(19)
!
 dep(1)=0.0
 do j=1,jo
   dep(j+1)=dep(j)+thick(j)
 enddo
 dep(jo+1)=amax1(dep(jo+1),depth+10.0)
 nb=2
 do j=1,jo
   if(depth>dep(j)) nb=j+1
 enddo
 if(abs(depth-dep(nb)) < 1.e-5) then
   nlayer=jo
 else if(abs(depth-dep(nb-1)) < 1.e-5) then
   nlayer=jo;   nb=nb-1
 else
   do j=jo,nb-1,-1
     j1=j+1
     afa(j1)  =afa(j)
     beta(j1) =beta(j)
     rho(j1)  =rho(j)
     thick(j1)=thick(j)
     qp(j1)   =qp(j)
     qs(j1)   =qs(j)
   enddo
   thick(nb-1)=depth-dep(nb-1)
   thick(nb)  =dep(nb)-depth
   nlayer=jo+1
 endif
end subroutine cmodel
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
MODULE synths_comm
 implicit NONE
 save
 integer:: nt1,nt3,npt,npt2
 real:: dt1,dt3,dt
 real, allocatable, dimension(:,:):: vel1d,vel3d
END MODULE synths_comm
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ReadVel(ism,stan_h1,stan_h2,st_name,path1D,path3D)
 use synths_comm
 implicit NONE
 character(len=*), intent(INOUT):: st_name,path1D,path3D
 integer,intent(IN):: ism
 real,intent(IN):: stan_h1,stan_h2
 character(len=256):: fh1,fh2,fup,chsm,smidx,prefix,suffix,junk
 character(len=3):: ch1,ch2,ch3
 integer:: nc_st,nc_p1,nc_p3,nc_sm,getlen,ntt,i,k
 real, allocatable, dimension(:,:):: syn3d
 integer:: comment, status, npt1
!
 call num2ch(float(ism)+0.01,smidx)
 call num2ch(stan_h1,ch1)
 call num2ch(stan_h2,ch2)
 ch3='ver'
 nc_p1=getlen(path1D)
 nc_p3=getlen(path3D)
 nc_st=getlen(st_name)
 nc_sm=getlen(smidx)
! Input 1D synthetics (linear)
!BBP format
   fh1=path1D(1:nc_p1)//st_name(1:nc_st)//'.3comp'
   open(unit=13,file=fh1,status='old',position='rewind',IOSTAT=status)
   comment=0;
   read(13,'(1a)') junk
   do while ((junk(1:1).eq.'%').or.(junk(1:1).eq.'#'))
	comment=comment+1;
	 read(13,'(1a)') junk
   end do
   npt1=1
   status=0
   do while(status.eq.0)
	 read(13,'(1a)',IOSTAT=status) junk
	 npt1=npt1+1
    end do
    npt1=npt1-1;
    nt3=npt1
    if(allocated(vel3d)) deallocate(vel3d)
    allocate(vel3d(nt3,4))
   close(13)
   open(unit=13,file=fh1,status='old',position='rewind',IOSTAT=status)
   do k=1,comment
	 read(13,'(1a)') junk
   end do
   read(13,*) (vel3d(i,4),vel3d(i,1),vel3d(i,2),vel3d(i,3), i=1,npt1)
   close(13)
   dt3=vel3d(2,4)-vel3d(1,4)
   write(*,*) dt3, nt3, comment
!
! Input nonlinear synthetics
 chsm='.nl1D.'//smidx(1:nc_sm)
 k=getlen(chsm)
 fh1=''; fh2=''; fup=''
 fh1=path3D(1:nc_p3)//st_name(1:nc_st)//'.'//ch1//chsm(1:k)
 fh2=path3D(1:nc_p3)//st_name(1:nc_st)//'.'//ch2//chsm(1:k)
 fup=path3D(1:nc_p3)//st_name(1:nc_st)//'.'//ch3//chsm(1:k)
 open(21,file=fh1,status='old', position='rewind')
 open(22,file=fh2,status='old', position='rewind')
 open(23,file=fup,status='old', position='rewind')
 read(21,*) nt1,dt1
 read(22,*) ntt,dt1
 nt1=min0(nt1,ntt)
 read(23,*) ntt,dt1
 nt1=min0(nt1,ntt)
 if(allocated(vel1d)) deallocate(vel1d)
 allocate(vel1d(nt1,3))
 read(21,*) (vel1d(i,1), i=1,nt1)
 read(22,*) (vel1d(i,2), i=1,nt1)
 read(23,*) (vel1d(i,3), i=1,nt1)
 close(21);  close(22); close(23)
!
! Interpolating to the same time step
 npt=ifix((nt3-1)*dt3/dt1)
 ntt=max0(npt,nt3)
 allocate(syn3d(ntt,3))
 do k=1,3
   do i=1,ntt
     syn3d(i,k) = 0.0
   enddo
   syn3d(1:nt3,k)=vel3d(1:nt3,k)
   if(abs(dt1-dt3).gt.0.01*dt1) then
     if(dt1.gt.dt3) call xapiir(syn3d(1,k),nt3,'BU',0.0,0.0,4, &
                             'LP', 0.001, 0.48/dt1, dt3, 2, ntt)
     call interpol(0.0,dt3,nt3,syn3d(1,k),dt1,npt, ntt)
   endif
 enddo
 npt=min0(npt, nt1)
 deallocate(vel3d);   allocate(vel3d(npt,3))
 vel3d(1:npt, :)=syn3d(1:npt, :)
 nt3=npt;  dt3=dt1
 npt2=2
 do while (npt2<npt*2)
   npt2=npt2*2
 enddo
 deallocate(syn3d)
end subroutine ReadVel
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine align(freq_join,ts_start)
 use synths_comm
 implicit NONE
 real, intent(IN):: freq_join,ts_start
 integer:: i,j,k,ntm1,ntm2,nwd1,nwd2,ish
 real:: vmx1,vmx2,crr,cross,sum1,sum2,sum3,shift_half
 real, dimension(npt):: vh1,vh2

 vh1=vel1d(1:npt,1);  vh2=vel1d(1:npt,2)
 shift_half=amin1(3.0,0.9*ts_start,amax1(0.1*ts_start,0.5))
 ntm1=-int(shift_half/dt1+0.5)
 nwd1=max0(int(amax1(0.0,(ts_start-1.0/freq_join))/dt1+1.5),1-ntm1)
 nwd2=min0(int((ts_start+2.0/freq_join)/dt1+1.5),npt)
 ntm2=min0(int(shift_half/dt1+0.5), npt-nwd2)
 call xapiir(vh1, npt, 'BU', 0.0, 0.0, 4, &
             'BP',freq_join/2.0,freq_join,dt1,2,npt)
 call xapiir(vh2, npt, 'BU', 0.0, 0.0, 4, &
             'BP',freq_join/2.0,freq_join,dt1,2,npt)
 vmx1=0.0;  vmx2=0.0
 do i=nwd1,nwd2
   vmx1=amax1(vmx1,abs(vh1(i)))
   vmx2=amax1(vmx2,abs(vh2(i)))
 enddo
 k=1;  if(vmx2 > vmx1) k=2
 vh1=vel1d(1:npt,k);  vh2=vel3d(1:npt,k)
 call xapiir(vh1, npt, 'BU', 0.0, 0.0, 4, &
             'BP',freq_join/2.0,freq_join,dt1,2,npt)
 call xapiir(vh2, npt, 'BU', 0.0, 0.0, 4, &
             'BP',freq_join/2.0,freq_join,dt1,2,npt)
 cross=0.0;  ish=0
 do j=ntm1,ntm2
   sum1=0.0;  sum2=0.0;  sum3=0.0
   do i=nwd1,nwd2
     sum1=sum1+vh1(i)*vh2(i+j)
     sum2=sum2+vh1(i)**2
     sum3=sum3+vh2(i+j)**2
   enddo
   crr=sum1/(sqrt(sum2)*sqrt(sum3))
   if(cross < crr) then
     ish=j
     cross = crr
   endif
 enddo

 do k=1,3
   if(ish < 0) then
     do i=npt,1-ish,-1
       vel3d(i,k)=vel3d(i+ish,k)
     enddo
     do i=1,-ish
       vel3d(i,k)=0.0
     enddo
   else
     do i=1,npt-ish
       vel3d(i,k)=vel3d(i+ish,k)
     enddo
     do i=npt-ish+1,npt
       vel3d(i,k)=0.0
     enddo
   endif
 enddo
end subroutine align
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine stitch_wt(freq_join)
 use synths_comm
 implicit NONE
 real, intent(IN):: freq_join
 integer:: k,kj,nb,ne,i
 real:: factor
 real, dimension(npt2):: syn1d,syn3d
!
 factor=real(npt2)*freq_join
 kj=1;  dt=2.0/factor
 do while(dt <= dt1)
   kj=kj+1
   dt=2.0**kj/factor
 enddo
 kj=kj-1
 dt=2**kj/factor
 ne=2**(kj+1);  nb=ne/2+1
 do k=1,3
   syn1d=0.0;   syn3d=0.0
   do i=1,npt2
      syn1d(i) = 0.0
      syn3d(i) = 0.0
   enddo
   syn1d(1:npt)=vel1d(1:npt,k)
   syn3d(1:npt)=vel3d(1:npt,k)
   call interpol(0.0, dt1,npt2,syn1d,dt,npt2,npt2)
   call interpol(0.0, dt3,npt2,syn3d,dt,npt2,npt2)
   call fwd_wt(npt2,syn1d)
   call fwd_wt(npt2,syn3d)
   syn1d(1:ne)=syn3d(1:ne)
   call inv_wt(npt2,syn1d)
   call xapiir(syn1d,npt2,'BU',0.0,0.0,4, &
               'LP',0.001,0.48/dt1,dt,2,npt2)
   call interpol(0.0,dt,npt2,syn1d,dt1,npt,npt2)
   vel1d(1:npt, k)=syn1d(1:npt)
 enddo
end subroutine stitch_wt
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine time_history(id_motion)
 use synths_comm
 implicit NONE
 integer, intent(IN):: id_motion
 integer:: k
 do k=1,3
   if(id_motion == 1) then
     call vel2dis(nt1,dt1,vel1d(1,k))
     call vel2dis(nt3,dt3,vel3d(1,k))
   else if(id_motion == 3) then
     call vel2acc(nt1,dt1,vel1d(1,k))
     call vel2acc(nt3,dt3,vel3d(1,k))
   endif
 enddo
end subroutine time_history
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine output3(ism,stan_h1,stan_h2,st_name)
 use synths_comm
 implicit NONE
 character(len=*), intent(IN):: st_name
 integer, intent(IN):: ism
 real, intent(IN):: stan_h1,stan_h2
 character(len=256):: fout,chsm
 character(len=3):: ch1,ch2,ch3
 integer:: i,getlen,nc_st

 call num2ch(float(ism)+0.01,ch3)
 chsm='.gmBB.'//ch3
 ch3='ver'
 call num2ch(stan_h1,ch1)
 call num2ch(stan_h2,ch2)
 nc_st=getlen(st_name)
 fout=''; fout=st_name(1:nc_st)//'.'//ch1//chsm(1:9)
 open(19, file=fout,status='replace',position='rewind')
 write(19,*) npt,dt1
 write(19,925) (vel1d(i,1), i=1,npt)
 close(19)
 fout=''; fout=st_name(1:nc_st)//'.'//ch2//chsm(1:9)
 open(19, file=fout,status='replace',position='rewind')
 write(19,*) npt,dt1
 write(19,925) (vel1d(i,2), i=1,npt)
 close(19)
 fout=''; fout=st_name(1:nc_st)//'.'//ch3//chsm(1:9)
 open(19, file=fout,status='replace',position='rewind')
 write(19,*) npt,dt1
 write(19,925) (vel1d(i,3), i=1,npt)
 close(19)
 925   format(5e14.5)
end subroutine output3
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine vel2acc(nt,dt,acc)
 implicit NONE
 integer, intent(IN) :: nt
 real, intent(IN) :: dt
 real, intent(INOUT), dimension(nt) :: acc
 integer :: i
 do i=1,nt-1
   acc(i)=(acc(i+1)-acc(i))/dt
 enddo
 acc(nt)=0.0
 do i=nt,2,-1
   acc(i)=0.5*(acc(i)+acc(i-1))
 enddo
 acc(1)=0.0
end subroutine vel2acc
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine vel2dis(nt,dt,acc)
 implicit NONE
 integer, intent(IN) :: nt
 real, intent(IN) :: dt
 real, intent(INOUT), dimension(nt) :: acc
 real :: v1,v2,dt05
 integer :: i
 dt05=0.5*dt;  v1=acc(1);  acc(1)=0.0
 do i=2,nt
   v2=v1;   v1=acc(i)
   acc(i)=acc(i-1)+(v2+v1)*dt05
 enddo
 return
end subroutine vel2dis
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
function getlen(string)
 implicit NONE
 character(len=*), intent(INOUT):: string
 character(len=256):: chatmp*256
 integer:: nl,getlen

 chatmp=adjustl(string)
 nl=len_trim(chatmp)
 string=''
 string=chatmp(1:nl)
 getlen=nl
end function getlen
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine num2ch(r,ch)
 implicit NONE
 real, intent(IN):: r
 character(len=*), intent(OUT):: ch
 integer:: i,n,n2,nl,nc,nzero

 ch='';  n=int(r+0.01)
 nc=3; if(n > 99) nc=3
 nl=10**(nc-1)
 nzero=ichar('0')
 do i=1,nc
   n2=n/nl
   ch(i:i)=char(nzero+n2)
   n=n-n2*nl
   nl=nl/10
 enddo
end subroutine num2ch
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
subroutine trav(dis,tmin)
! calculate travel time for horizontal layered model
! it outputs both times for first arrival and direct arrival
!
 use ccmodel
 implicit none
 real, intent(IN):: dis
 real, intent(OUT):: tmin
 real(kind=8):: t, t0, td, x
 complex(kind=8):: pd, p0, p
 integer:: i

 ns=nb-1
 do i=1, nlayer
   th(i)=dble(thick(i))
   v(i)=1.0d0/dble(beta(i))/dble(beta(i))
 enddo

 x=dble(dis)
 lmax = ns
 call find2(x,td,pd)
 t0 = td
 p0 = pd
 do lmax=ns+1, nlayer-1
   call find2(x,t,p)
   if (t < t0) then
     t0=t
     p0=p
   endif
 enddo
 tmin=sngl(t0)
end subroutine trav
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine find2(x,t0,p0)
! solve d(tau)/dp = 0 for p0, tau=p*x+eta*z
! input:  x distance range
! output: p0, t0
!
 use ccmodel
 implicit none
 real(kind=8), intent(IN):: x
 real(kind=8), intent(OUT):: t0
 complex(kind=8), intent(OUT):: p0
 complex(kind=8):: taup,dtdp,p1,p2
 real(kind=8):: aminm,dtdp0,pm
 integer:: i

 aminm = 1.d-6
 p1 = dcmplx(0.d0,1.0d-20)    ! make p slightly above the real axis
 pm = 9999999.d0
 do i=1,lmax
   pm = dmin1(pm,v(lmax))
 enddo
 p2 = dsqrt(pm) + p1
 p0 = 0.5d0*(p1+p2)          ! the actual p0 lies betw. p1 and p2
 do while( dble(p2-p1) > aminm )
   dtdp0 = dtdp(x,p0)
   if( dabs(dtdp0) < aminm ) exit
   if( dtdp0 > 0.0d0 ) then
     p1 = p0
   else
     p2 = p0
   endif
   p0 = 0.5d0*(p1+p2)
 enddo
 pm = dmin1(pm,v(lmax+1))
 pm = dsqrt(pm)
 if (lmax>ns .and. pm<dble(p0)) p0 = pm
 t0 = taup(p0,x)
end subroutine find2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function taup(p,x)
! define function tau(p) = p x + eta h
 use ccmodel
 implicit none
 complex(kind=8), intent(IN):: p
 real(kind=8), intent(IN):: x
 complex(kind=8):: pp,taup
 integer:: i

 taup = p*x
 pp = p*p
 do i = 1, ns
   taup=taup+cdsqrt(v(i)-pp)*th(i)
 enddo
 do i = ns+1,lmax
   taup=taup+cdsqrt(v(i)-pp)*th(i)*2.d0
 enddo
end function taup
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function dtdp(x,p)
! define d(tau)/dp
 use ccmodel
 implicit none
 complex(kind=8), intent(IN):: p
 real(kind=8), intent(IN):: x
 complex(kind=8):: pp,dtdp
 integer:: j
 pp = p*p
 dtdp = cmplx(0.0d0,0.0d0)
 do j = 1, ns
   dtdp=dtdp-th(j)/cdsqrt(v(j)-pp)
 enddo
 do j = ns+1, lmax
   dtdp=dtdp-2*th(j)/cdsqrt(v(j)-pp)
 enddo
 dtdp = x + p*dtdp
end function dtdp
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
module wavelet_base
 implicit NONE
 SAVE
 integer :: num_pt=0, iwt_save=0
 real, allocatable, dimension(:) :: rwt
end module wavelet_base
!
subroutine inv_wt(nlen,u)
! nlen: the size of u. nlen must be a power of 2 (nlen=2**N)
! u: inputting values are wavelet coefficients from fwd_wt
!    outputs are time histories.
! The inputted wavelet coefficients are not multiplied sqrt(dt)
! level(j=0,1, N-1),  input-wavelet-coefficients(u),
! #-of-wavelet-coefficients
!     0                     u(2)                               1
!     1                     u(3), u(4)                         2
!     j                     u(2**j+1 : 2**(j+1))               2**j
!
! USES: real_fft.f90
!
 use wavelet_base
 implicit NONE
 integer, intent(IN) :: nlen
 real, intent(INOUT), dimension(nlen):: u
 real, dimension(2*nlen) :: wcf,fu
 integer :: i1,i2,k1,k2,kmax,kstep,i
 real :: cff
 real :: cp1,cp2, pi2
 complex :: wave, cc
 if(num_pt /= nlen .or. iwt_save/=1) then
   if(allocated(rwt)) deallocate(rwt)
   allocate(rwt(4*nlen))
   do i=1,4*nlen
      rwt(i) = 0.0
   enddo
   do i=1,2*nlen
      wcf(i) = 0.0
      fu(i) = 0.0
   enddo
   num_pt = nlen
   iwt_save = 1
   kmax=1
   i2=0
   do while( kmax < nlen)
     do k1=1,2*kmax
       i1=i2+2*k1-1
       cc=wave(nlen,kmax,k1,1)
       rwt(i1)  =real(cc)
       rwt(i1+1)=aimag(cc)
     enddo
     i2=i2+4*kmax
     kmax=kmax*2
   enddo
 endif
 wcf=0.0
 k1=nlen+1
 k2=k1+1
 i1=nlen/2+1
 i2=i1+1
! j=0
 wcf(1)=u(1)
! j=1, kmax=1
 wcf(3)= u(2)*rwt(3)
 wcf(4)= u(2)*rwt(4)
! j=2, kmax=2
 wcf(3)= wcf(3)+rwt(7)*(u(3)-u(4))
 wcf(4)= wcf(4)+rwt(8)*(u(3)-u(4))
 wcf(5)=        rwt(9)*(u(3)+u(4))
 wcf(6)=        rwt(10)*(u(3)+u(4))
! j=3
 kmax=4
 do
   do i1=1,kmax
     fu(i1)= u(i1+kmax)
   enddo
   call realft(fu,kmax,-1)
   fu(kmax+1)=fu(2)
   fu(kmax+2)=0.0
   fu(2)=0.0
   k2=2*kmax+2
   do i1=3,kmax,2
     k1=k2-i1
     fu(k1)  = fu(i1)
     fu(k1+1)=-fu(i1+1)
   enddo
   kmax=kmax*2
   kstep=2*kmax-4
   do i2=4,kmax,2
     i1=i2-1
     k1=i1+kstep
     wcf(i1)=wcf(i1)+fu(i1)*rwt(k1)  -fu(i2)*rwt(k1+1)
     wcf(i2)=wcf(i2)+fu(i1)*rwt(k1+1)+fu(i2)*rwt(k1)
   enddo
!
   k2=kmax/2
   kstep=3*kmax-4
   do i2=2,k2,2
     i1=i2-1
     k1=i1+kstep
     wcf(i1+kmax)=fu(i1)*rwt(k1)  -fu(i2)*rwt(k1+1)
     wcf(i2+kmax)=fu(i1)*rwt(k1+1)+fu(i2)*rwt(k1)
   enddo
   if(kmax == nlen) EXIT
 enddo
 call four1(wcf,nlen,1)
 cff=2./float(nlen)
 do i1=1,nlen
   u(i1)=cff*wcf(2*i1-1)
 enddo
 return
end
!
subroutine fwd_wt(nlen,u)
! nlen: the size of u. nlen must be a power of 2
! u: inputting values are time histories.
!    outputs are wavelet coefficients at different level
!    (frequencies band).
! The computed wavelet coefficients are not multiplied sqrt(dt)
! level(j=0,1, N-1),  coefficients(u),    #-of-wavelet-coefficients
!     0                u(2)                               1
!     1                u(3), u(4)                         2
!     j                u(2**j+1 : 2**(j+1))               2**j
!
 use wavelet_base
 implicit NONE
 integer, intent(IN) :: nlen
 real, intent(INOUT), dimension(nlen):: u
 real, dimension(2*nlen) :: fu
 real, dimension(nlen) :: wcf
 integer :: j,i1,i2,k1,k2,kmax,kstep
 real :: cff
 complex:: wave, cc
 if(num_pt /= nlen .or. iwt_save/=2) then
   if(allocated(rwt)) deallocate(rwt)
   allocate(rwt(4*nlen))
   num_pt = nlen
   do j=1,nlen*4
      rwt(j) = 0.0
   enddo
   do j=1,2*nlen
      fu(j) = 0.0
   enddo
   do j=1,nlen
      wcf(j) = 0.0
   enddo
   iwt_save = 2
   kmax=1
   i2=0
   do while( kmax < nlen)
     do k1=1,2*kmax
       i1=i2+2*k1-1
       cc=wave(nlen,kmax,k1,2)
       rwt(i1)  =real(cc)
       rwt(i1+1)=aimag(cc)
     enddo
     i2=i2+4*kmax
     kmax=kmax*2
   enddo
 endif
!
 cff=1./real(nlen)
 do i1=1,nlen
   i2=2*i1
   fu(i2-1) = cff*u(i1)
   fu(i2)   = 0.0
 enddo
 call four1(fu,nlen,-1)
!
! Zero frequency
 u(1)=fu(1)
!j=0, kmax=1
 u(2)=2.*(fu(3)*rwt(3)-fu(4)*rwt(4))
!j=1, kmax=2
 u(3)=2.*(fu(5)*rwt(9)-fu(6)*rwt(10)+fu(3)*rwt(7)-fu(4)*rwt(8))
 u(4)=2.*(fu(5)*rwt(9)-fu(6)*rwt(10)-fu(3)*rwt(7)+fu(4)*rwt(8))
! j=2
 kmax=4
 do while(kmax < nlen)
   kmax=2*kmax
   kstep=2*kmax-4
   do i1=1,kmax,2
     i2=i1+kmax
     k1=kstep+i1
     k2=kstep+i2
     wcf(i1)  =fu(i1)*rwt(k1)-fu(i1+1)*rwt(k1+1)+ &
               fu(i2)*rwt(k2)-fu(i2+1)*rwt(k2+1)
     wcf(i1+1)=fu(i1)*rwt(k1+1)+fu(i1+1)*rwt(k1)+ &
               fu(i2)*rwt(k2+1)+fu(i2+1)*rwt(k2)
   enddo
   k2=kmax/2
   call four1(wcf,k2,1)
   do k1=1,kmax,2
     k2=k2+1
     u(k2)=2.*wcf(k1)
   enddo
 enddo
 return
end
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function wave(nlen,kmax,iw,ic)
 implicit NONE
 integer, parameter :: double=kind(0.0d0)
 real(kind=double), parameter :: pi=3.1415926535897d0, &
                   p23=2.0d0*pi/3.0d0,p43=4.0d0*pi/3.0d0
 integer, intent(IN) :: nlen,kmax,iw,ic
 real(kind=double), dimension(2):: g,p
 real(kind=double):: w,hw,ww,wd,wf,wg,wp,fw
 real:: r1,r2
 integer:: i,j,k
 complex:: wave
 wd=0.0d0
 w=2.0d0*pi*real(iw-1,double)/real(kmax,double)
 hw=w*0.5d0
 ww=dabs(w)
 wave=cmplx(0.0,0.0)
 if(ww >(2.0d0*p43).or.ww < p23) return
 do i=1,2
   wp=hw*i
   fw=1.0d0
   do j=1,2
     if(j == 1)then
       wf=-wp
     else
       wf=wp
     end if
     do k=1,2
       if(k == 1) wd=p43-wf
       if(k == 2) wd=wf-p23
       if(wd <= 0.0d0)then
         g(k)=0.0d0
       else
        g(k)=exp(-1.0/wd/wd)
       end if
     enddo
     wg=g(1)/(g(1)+g(2))
     fw=fw*wg
   enddo
   p(i)=fw
 enddo
 r1=real(hw)
 r2=real(p(1)-p(2))*float(nlen)/float(kmax)
 wave = cmplx(cos(r1),-sin(r1))*csqrt(cmplx(r2, 0.0))
 if(ic == 1) then
!   w= 2.d0*pi*real(iw-1,double)/real(nlen*dt,double)
!   wave=wave*cmplx(-sngl(w*w),0.0)
 else
   wave = conjg(wave)
 endif
 return
end function wave

max_stars_repo_name: kevinmilner/bbp

max_issues_repo_name: kevinmilner/bbp

max_forks_repo_name: kevinmilner/bbp


text:
*
*     -------------------------------------------------------------
*      T W O P A R T I C L E 3 1
*     -------------------------------------------------------------
*                                                                  *
*     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
*     OF TWO PARTICLE OPERATOR IN CASE :       N'1 = N1 - 1        *
*                                              N'2 = N2 + 1        *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University,  Nashville             October 1996   *
*
      SUBROUTINE TWOPARTICLE31(IA,IB,IIA,IIB,IIC,IID)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/MEDEFN/IHSH,NJ(16),LJ(16),NOSH1(16),NOSH2(16),J1QN1(31,3),
     :     J1QN2(31,3),IJFUL(16)
      COMMON/TRK/BD1(3),BD2(3),BK1(3),BK2(3),
     *ID1(7),ID2(7),IK1(7),IK2(7)
      COMMON/PERMAT/IRS(2,2),IRL(20,20),RS(2,2),RL(20,20)
      COMMON/KAMPAS/ IW1(2,20),IW2(2,20),IWAA(2,2,20,20),
     :RW1(2,20),RW2(2,20),RWAA(2,2,20,20)
      COMMON /CASEOP/ IOCASE
      EXTERNAL TWO31
      IF(IHSH.LE.1)RETURN
      CALL HIBFF(IA,IB,IA,IA,2)
      IF(ID1(4).LT.2)RETURN
      IF(ITTK(1,IK2(6),ID2(6)).EQ.0) RETURN
      IF(ITTK(2*IK2(3),IK2(5),ID2(5)).EQ.0) RETURN
      IF(IABS(IK1(6)-ID1(6)).GT.3) RETURN
      IF(IABS(IK1(5)-ID1(5)).GT.6*IK1(3)) RETURN
      IF(IOCASE.EQ.1) THEN
* Spin-Spin operator
        IP=ITREXG2(LJ(IA),LJ(IA),LJ(IA),LJ(IB),IKK)+1
      ELSEIF(IOCASE.EQ.2) THEN
* Spin-other-orbit operator
        IP=ITREXG(LJ(IA),LJ(IA),LJ(IA),LJ(IB),IKK)+1
      ENDIF
      IF(IKK.LE.0)RETURN
      IG=IP+IKK-1
      DO 2 I1=1,2
        IRS(I1,1)=0
        DO 3 I2=1,20
          IRL(I2,1)=0
          DO 4 I3=1,2
            DO 5 I4=1,20
	      IWAA(I1,I3,I2,I4)=0
    5       CONTINUE
    4     CONTINUE
    3   CONTINUE
    2 CONTINUE
      DO 1 I=IP,IG,2
        KL=I-1
C   COULOMB
C        CALL TWO31(KL,0,KL,0,0,IA,IB,IIA,IIB,IIC,IID,C1,C2)
C   COULOMB
*
        IF(IOCASE.EQ.1) THEN
* Spin-spin operator
          CALL SSC(IG,KL,IA,IB,IB,IB,IB,IA,IA,IA,KL,TWO31)
        ELSEIF(IOCASE.EQ.2) THEN
* Spin-other-orbit operator
          CALL SOOC(IG,KL,IA,IB,IB,IB,IB,IA,IA,IA,KL,TWO31)
        ENDIF
*
    1 CONTINUE
      RETURN
      END

ext: f


text:
! Copyright (c) 1999, 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.
!
! Target of a pointer assignment is a derived type array section with member.
! reported by a user in pgf90 3.0

module m
  type et
   integer:: v
   complex::xxyy
   double precision :: d(2)
   integer:: v2(2)
  end type
end module

  use m
  integer,pointer,dimension(:) :: p,q,r,s,t,u,w

  integer, dimension(37) :: results, expect
  data expect/100,200,300,400,500,600,700,800,900,1000, &
             &100,200,300,400,500,600,700,800,900,1000, &
             &200,400,600,800,1000,100,300,500,700,900, &
             &200,400,600,800,1000,400,800/

  type(et),dimension(10),target :: e

  do i = 1,10
   e(i)%v = 100*i
   e(i)%xxyy = (-1.0,-10.0)*i
   e(i)%d(1) = 99
   e(i)%d(2) = 98
   e(i)%v2(1) = 1000*i
   e(i)%v2(2) = 999*i
  enddo

 5 format(a,i2,a,i2,a)
10 format( 10i6 )

  q => e%v
  !print 5,'q(',lbound(q,1),':',ubound(q,1),') => e%v         = (100,200,...,1000)'
  !print 10,q
  results(1:10) = q

  r => q
  !print 5,'r(',lbound(r,1),':',ubound(r,1),') => q           = (100,200,...,1000)'
  !print 10,r
  results(11:20) = r

  s => r(2:10:2)
  !print 5,'s(',lbound(s,1),':',ubound(s,1),') => r(2:10:2)   = (200,400,...,1000)'
  !print 10,s
  results(21:25) = s

  t => s(1:5:2)
  !print 5,'t(',lbound(t,1),':',ubound(t,1),') => s(1:5:2)    = (200,600,1000)'
  !print 10,t
  results(21:25) = s

  u => e(1:10:2)%v
  !print 5,'u(',lbound(u,1),':',ubound(u,1),') => e(1:10:2)%v = (100,300,...,900)'
  !print 10,u
  results(26:30) = u

  w => e(2:10:2)%v
  !print 5,'w(',lbound(u,1),':',ubound(u,1),') => e(2:10:2)%v = (200,400,...,1000)'
  !print 10,w
  results(31:35) = w

  p => w(2:5:2)
  !print 5,'p(',lbound(u,1),':',ubound(u,1),') => w(2:5:2)    = (400,800)'
  !print 10,p
  results(36:37) = p
  call check(results,expect,size(expect))
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


! { dg-do run }
!
! PR fortran/52022
!

module check
  integer, save :: icheck = 0
end module check

module t
implicit none
      contains
subroutine  sol(cost)
   use check
   interface
        function cost(p) result(y)
                double precision,dimension(:) :: p
                double precision,dimension(:),allocatable :: y
        end function cost
   end interface

   if (any (cost([1d0,2d0]) /= [2.d0, 4.d0])) STOP 1
   icheck = icheck + 1
end subroutine

end module t

module tt
   procedure(cost1),pointer :: pcost
contains
  subroutine init()
        pcost=>cost1
  end subroutine

  function cost1(x) result(y)
        double precision,dimension(:) :: x
        double precision,dimension(:),allocatable :: y
        allocate(y(2))
        y=2d0*x
  end function cost1

  function cost(x) result(y)
        double precision,dimension(:) :: x
        double precision,dimension(:),allocatable :: y
        allocate(y(2))
        y=pcost(x)
  end function cost
end module

program test
        use tt
        use t
        use check
        implicit none

        call init()
        if (any (cost([3.d0,7.d0]) /= [6.d0, 14.d0])) STOP 2
        if (icheck /= 0) STOP 3
        call sol(cost)
        if (icheck /= 1) STOP 4
end program test


text:
CHARACTER xyz(3)
INTEGER(4) idata(35)
DATA idata/35*-1/,xyz/'x', 'y', 'z'/
OPEN(3,FILE='USEQ',FORM='UNFORMATTED')
WRITE(3) idata
WRITE(3) xyz
CLOSE(3)
END

max_stars_repo_name: quchunguang/test


text:
*> \brief \b CLAHRD
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CLAHRD + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
*
*       .. Scalar Arguments ..
*       INTEGER            K, LDA, LDT, LDY, N, NB
*       ..
*       .. Array Arguments ..
*       COMPLEX            A( LDA, * ), T( LDT, NB ), TAU( NB ),
*      $                   Y( LDY, NB )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CLAHRD reduces the first NB columns of a complex general n-by-(n-k+1)
*> matrix A so that elements below the k-th subdiagonal are zero. The
*> reduction is performed by a unitary similarity transformation
*> Q**H * A * Q. The routine returns the matrices V and T which determine
*> Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T.
*>
*> This is an OBSOLETE auxiliary routine.
*> This routine will be 'deprecated' in a  future release.
*> Please use the new routine CLAHR2 instead.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*>          The offset for the reduction. Elements below the k-th
*>          subdiagonal in the first NB columns are reduced to zero.
*> \endverbatim
*>
*> \param[in] NB
*> \verbatim
*>          NB is INTEGER
*>          The number of columns to be reduced.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX array, dimension (LDA,N-K+1)
*>          On entry, the n-by-(n-k+1) general matrix A.
*>          On exit, the elements on and above the k-th subdiagonal in
*>          the first NB columns are overwritten with the corresponding
*>          elements of the reduced matrix; the elements below the k-th
*>          subdiagonal, with the array TAU, represent the matrix Q as a
*>          product of elementary reflectors. The other columns of A are
*>          unchanged. See Further Details.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*>          TAU is COMPLEX array, dimension (NB)
*>          The scalar factors of the elementary reflectors. See Further
*>          Details.
*> \endverbatim
*>
*> \param[out] T
*> \verbatim
*>          T is COMPLEX array, dimension (LDT,NB)
*>          The upper triangular matrix T.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T.  LDT >= NB.
*> \endverbatim
*>
*> \param[out] Y
*> \verbatim
*>          Y is COMPLEX array, dimension (LDY,NB)
*>          The n-by-nb matrix Y.
*> \endverbatim
*>
*> \param[in] LDY
*> \verbatim
*>          LDY is INTEGER
*>          The leading dimension of the array Y. LDY >= max(1,N).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup complexOTHERauxiliary
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  The matrix Q is represented as a product of nb elementary reflectors
*>
*>     Q = H(1) H(2) . . . H(nb).
*>
*>  Each H(i) has the form
*>
*>     H(i) = I - tau * v * v**H
*>
*>  where tau is a complex scalar, and v is a complex vector with
*>  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
*>  A(i+k+1:n,i), and tau in TAU(i).
*>
*>  The elements of the vectors v together form the (n-k+1)-by-nb matrix
*>  V which is needed, with T and Y, to apply the transformation to the
*>  unreduced part of the matrix, using an update of the form:
*>  A := (I - V*T*V**H) * (A - Y*V**H).
*>
*>  The contents of A on exit are illustrated by the following example
*>  with n = 7, k = 3 and nb = 2:
*>
*>     ( a   h   a   a   a )
*>     ( a   h   a   a   a )
*>     ( a   h   a   a   a )
*>     ( h   h   a   a   a )
*>     ( v1  h   a   a   a )
*>     ( v1  v2  a   a   a )
*>     ( v1  v2  a   a   a )
*>
*>  where a denotes an element of the original matrix A, h denotes a
*>  modified element of the upper Hessenberg matrix H, and vi denotes an
*>  element of the vector defining H(i).
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE CLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
*
*  -- LAPACK auxiliary routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            K, LDA, LDT, LDY, N, NB
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), T( LDT, NB ), TAU( NB ),
     $                   Y( LDY, NB )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ),
     $                   ONE = ( 1.0E+0, 0.0E+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      COMPLEX            EI
*     ..
*     .. External Subroutines ..
      EXTERNAL           CAXPY, CCOPY, CGEMV, CLACGV, CLARFG, CSCAL,
     $                   CTRMV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
*
      DO 10 I = 1, NB
         IF( I.GT.1 ) THEN
*
*           Update A(1:n,i)
*
*           Compute i-th column of A - Y * V**H
*
            CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
            CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
     $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
            CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
*
*           Apply I - V * T**H * V**H to this column (call it b) from the
*           left, using the last column of T as workspace
*
*           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
*                    ( V2 )             ( b2 )
*
*           where V1 is unit lower triangular
*
*           w := V1**H * b1
*
            CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
            CALL CTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
     $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
*
*           w := w + V2**H *b2
*
            CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
     $                  A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
     $                  T( 1, NB ), 1 )
*
*           w := T**H *w
*
            CALL CTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
     $                  T, LDT, T( 1, NB ), 1 )
*
*           b2 := b2 - V2*w
*
            CALL CGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
     $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
*
*           b1 := b1 - V1*w
*
            CALL CTRMV( 'Lower', 'No transpose', 'Unit', I-1,
     $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
            CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
*
            A( K+I-1, I-1 ) = EI
         END IF
*
*        Generate the elementary reflector H(i) to annihilate
*        A(k+i+1:n,i)
*
         EI = A( K+I, I )
         CALL CLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
     $                TAU( I ) )
         A( K+I, I ) = ONE
*
*        Compute  Y(1:n,i)
*
         CALL CGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
     $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
         CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
     $               A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
     $               1 )
         CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
     $               ONE, Y( 1, I ), 1 )
         CALL CSCAL( N, TAU( I ), Y( 1, I ), 1 )
*
*        Compute T(1:i,i)
*
         CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
         CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
     $               T( 1, I ), 1 )
         T( I, I ) = TAU( I )
*
   10 CONTINUE
      A( K+NB, NB ) = EI
*
      RETURN
*
*     End of CLAHRD
*
      END

max_stars_repo_name: LuaDist/lapack

max_issues_repo_name: LuaDist/lapack

max_forks_repo_name: LuaDist/lapack


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!   C                                                                      C
!   C  References:                                                         C
!   C          [1] Duff, I., Marrone, M., Radicati, G., and Vittoli, C.    C
!   C              Level 3 basic linear algebra subprograms for sparse     C
!   C              matrices: a user level interface                        C
!   C              ACM Trans. Math. Softw., 23(3), 379-401, 1997.          C
!   C                                                                      C
!   C                                                                      C
!   C         [2]  S. Filippone, M. Colajanni                              C
!   C              PSBLAS: A library for parallel linear algebra           C
!   C              computation on sparse matrices                          C
!   C              ACM Trans. on Math. Softw., 26(4), 527-550, Dec. 2000.  C
!   C                                                                      C
!   C         [3] M. Arioli, I. Duff, M. Ruiz                              C
!   C             Stopping criteria for iterative solvers                  C
!   C             SIAM J. Matrix Anal. Appl., Vol. 13, pp. 138-144, 1992   C
!   C                                                                      C
!   C                                                                      C
!   C         [4] R. Barrett et al                                         C
!   C             Templates for the solution of linear systems             C
!   C             SIAM, 1993
!   C                                                                      C
!   C                                                                      C
!   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! File:  psb_scgstab.f90
!
! Subroutine: psb_scgstab
!    This subroutine implements the BiCG Stabilized method.
!
!
! Arguments:
!
!    a      -  type(psb_sspmat_type)      Input: sparse matrix containing A.
!    prec   -  class(psb_sprec_type)       Input: preconditioner
!    b      -  real,dimension(:)       Input: vector containing the
!                                         right hand side B
!    x      -  real,dimension(:)       Input/Output: vector containing the
!                                         initial guess and final solution X.
!    eps    -  real                       Input: Stopping tolerance; the iteration is
!                                         stopped when the error estimate |err| <= eps
!    desc_a -  type(psb_desc_type).       Input: The communication descriptor.
!    info   -  integer.                   Output: Return code
!
!    itmax  -  integer(optional)          Input: maximum number of iterations to be
!                                         performed.
!    iter   -  integer(optional)          Output: how many iterations have been
!                                         performed.
!                                         performed.
!    err    -  real   (optional)          Output: error estimate on exit. If the
!                                         denominator of the estimate is exactly
!                                         0, it is changed into 1.
!    itrace -  integer(optional)          Input: print an informational message
!                                         with the error estimate every itrace
!                                         iterations
!    istop  -  integer(optional)          Input: stopping criterion, or how
!                                         to estimate the error.
!                                         1: err =  |r|/(|a||x|+|b|);  here the iteration is
!                                            stopped when  |r| <= eps * (|a||x|+|b|)
!                                         2: err =  |r|/|b|; here the iteration is
!                                            stopped when  |r| <= eps * |b|
!                                         where r is the (preconditioned, recursive
!                                         estimate of) residual.
!
Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
  use psb_base_mod
  use psb_prec_mod
  use psb_s_krylov_conv_mod
  use psb_krylov_mod
  implicit none
  type(psb_sspmat_type), intent(in)  :: a
  class(psb_sprec_type), Intent(inout)  :: prec
  Type(psb_desc_type), Intent(in)    :: desc_a
  type(psb_s_vect_type), Intent(inout)   :: b
  type(psb_s_vect_type), Intent(inout)   :: x
  Real(psb_spk_), Intent(in)         :: eps
  integer(psb_ipk_), intent(out)               :: info
  integer(psb_ipk_), Optional, Intent(in)      :: itmax, itrace, istop
  integer(psb_ipk_), Optional, Intent(out)     :: iter
  Real(psb_spk_), Optional, Intent(out) :: err
! =   Local data
  real(psb_spk_), allocatable, target   :: aux(:),wwrk(:,:)
  type(psb_s_vect_type) :: q, r, p, v, s, t, z, f

  integer(psb_ipk_) :: itmax_, naux, it,itrace_,&
       & n_row, n_col
  integer(psb_lpk_) :: mglob
  integer(psb_ipk_) :: debug_level, debug_unit
  Logical, Parameter :: exchange=.True., noexchange=.False., debug1 = .False.
  integer(psb_ipk_), Parameter :: irmax = 8
  integer(psb_ipk_) :: itx, err_act, i
  integer(psb_ipk_) :: istop_
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: np, me
  real(psb_dpk_)     :: derr
  real(psb_spk_)     :: alpha, beta, rho, rho_old, sigma, omega, tau
  type(psb_itconv_type) :: stopdat

  character(len=20)           :: name
  character(len=*), parameter :: methdname='BiCGStab'

  info = psb_success_
  name = 'psb_scgstab'
  call psb_erractionsave(err_act)
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()
  ctxt = desc_a%get_context()
  call psb_info(ctxt, me, np)
  if (debug_level >= psb_debug_ext_)&
       & write(debug_unit,*) me,' ',trim(name),': from psb_info',np
  if (.not.allocated(x%v)) then
    info = psb_err_invalid_vect_state_
    call psb_errpush(info,name)
    goto 9999
  endif
  if (.not.allocated(b%v)) then
    info = psb_err_invalid_vect_state_
    call psb_errpush(info,name)
    goto 9999
  endif

  mglob = desc_a%get_global_rows()
  n_row = desc_a%get_local_rows()
  n_col = desc_a%get_local_cols()

  If (Present(istop)) Then
    istop_ = istop
  Else
    istop_ = 2
  Endif
  !
  !  ISTOP_ = 1:  Normwise backward error, infinity norm
  !  ISTOP_ = 2:  ||r||/||b||   norm 2
  !
  ! =  if (.not.same_type_as(x,b)) then
  ! =    write(0,*) 'Warning: different dynamic types for X and B '
  ! =  end if

  call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
  if (psb_errstatus_fatal()) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='psb_chkvect on X')
    goto 9999
  end if
  call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
  if (psb_errstatus_fatal()) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='psb_chkvect on B')
    goto 9999
  end if
  naux=6*n_col
  if (info == psb_success_) allocate(aux(naux),stat=info)
  if (info /= psb_success_) then
    info=psb_err_alloc_dealloc_
    call psb_errpush(info,name)
    goto 9999
  End If

  call psb_geasb(q,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(r,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(p,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(v,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(s,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(t,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(z,desc_a,info,mold=x%v,scratch=.true.)
  call psb_geasb(f,desc_a,info,mold=x%v,scratch=.true.)

  if (psb_errstatus_fatal()) then
    info=psb_err_from_subroutine_non_
    call psb_errpush(info,name)
    goto 9999
  End If

  If (Present(itmax)) Then
    itmax_ = itmax
  Else
    itmax_ = 1000
  Endif

  If (Present(itrace)) Then
    itrace_ = itrace
  Else
    itrace_ = 0
  End If

  itx   = 0
  call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
  if (psb_errstatus_fatal()) Then
    call psb_errpush(psb_err_from_subroutine_non_,name)
    goto 9999
  End If

  restart: Do

    if (itx >= itmax_) exit restart

    it = 0
    call psb_geaxpby(sone,b,szero,r,desc_a,info)

    call psb_spmm(-sone,a,x,sone,r,desc_a,info,work=aux)
    call psb_geaxpby(sone,r,szero,q,desc_a,info)

    ! Perhaps we already satisfy the convergence criterion...
    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart

    if (psb_errstatus_fatal()) then
      info=psb_err_from_subroutine_
      call psb_errpush(info,name,a_err='Init residual chk')
      goto 9999
    end if

    rho = szero

    iteration:  Do
      it   = it + 1
      itx = itx + 1

      if (debug_level >= psb_debug_ext_)&
           & write(debug_unit,*) me,' ',trim(name),&
           & ' Iteration: ',itx

      rho_old = rho
      rho     = psb_gedot(q,r,desc_a,info)

      if (debug_level >= psb_debug_ext_) then
        call q%sync()
        call r%sync()
      end if
      if (debug_level >= psb_debug_ext_)&
           & write(debug_unit,*) me,' ',trim(name),&
           & ' Rho: ',rho

      if (rho == szero) then
        if (debug_level >= psb_debug_ext_) &
             & write(debug_unit,*) me,' ',trim(name),&
             & ' Iteration breakdown R',rho
        exit iteration
      endif

      if (it == 1) then
        call psb_geaxpby(sone,r,szero,p,desc_a,info)
      else
        beta = (rho/rho_old)*(alpha/omega)
        call psb_geaxpby(-omega,v,sone,p,desc_a,info)
        call psb_geaxpby(sone,r,beta,p,desc_a,info)
      End If

      call prec%apply(p,f,desc_a,info,work=aux)

      call psb_spmm(sone,a,f,szero,v,desc_a,info,&
           & work=aux)

      sigma = psb_gedot(q,v,desc_a,info)

      if (debug_level >= psb_debug_ext_) then
        call q%sync()
        call v%sync()
      end if
      if (debug_level >= psb_debug_ext_)&
           & write(debug_unit,*) me,' ',trim(name),&
           & ' Sigma: ',sigma

      if (sigma == szero) then
        if (debug_level >= psb_debug_ext_) &
             & write(debug_unit,*) me,' ',trim(name),&
             & ' Iteration breakdown S1', sigma
        exit iteration
      endif

      alpha = rho/sigma
      call psb_geaxpby(sone,r,szero,s,desc_a,info)
      call psb_geaxpby(-alpha,v,sone,s,desc_a,info)
      if (debug_level >= psb_debug_ext_)&
           & write(debug_unit,*) me,' ',trim(name),&
           & ' alpha: ',alpha

      if (psb_errstatus_fatal()) then
        call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_geaxpby')
        goto 9999
      end if

      call prec%apply(s,z,desc_a,info,work=aux)
      Call psb_spmm(sone,a,z,szero,t,desc_a,info,work=aux)

      if(psb_errstatus_fatal()) then
        call psb_errpush(psb_err_from_subroutine_,name,a_err='precaply/spmm')
        goto 9999
      end if

      sigma = psb_gedot(t,t,desc_a,info)
      if (sigma == szero) then
        if (debug_level >= psb_debug_ext_) &
             & write(debug_unit,*) me,' ',trim(name),&
             & ' Iteration breakdown S2', sigma
        exit iteration
      endif

      tau   = psb_gedot(t,s,desc_a,info)
      omega = tau/sigma

      if (debug_level >= psb_debug_ext_) then
        call t%sync()
        call s%sync()
      end if
      if (debug_level >= psb_debug_ext_)&
           & write(debug_unit,*) me,' ',trim(name),&
           & ' sigma, tau, omega: ',sigma, tau, omega

      if (omega == szero) then
        if (debug_level >= psb_debug_ext_) &
             & write(debug_unit,*) me,' ',trim(name),&
             & ' Iteration breakdown O',omega
        exit iteration
      endif

      call psb_geaxpby(alpha,f,sone,x,desc_a,info)
      call psb_geaxpby(omega,z,sone,x,desc_a,info)
      call psb_geaxpby(sone,s,szero,r,desc_a,info)
      call psb_geaxpby(-omega,t,sone,r,desc_a,info)
      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart

      if (psb_errstatus_fatal()) Then
        call psb_errpush(psb_err_from_subroutine_,name,a_err='X/R update ')
        goto 9999
      End If

    end do iteration
  end do restart

  call psb_end_conv(methdname,itx,desc_a,stopdat,info,derr,iter)
  if (present(err)) err = derr

  deallocate(aux,stat=info)

  call x%sync()
  call psb_gefree(q,desc_a,info)
  call psb_gefree(r,desc_a,info)
  call psb_gefree(p,desc_a,info)
  call psb_gefree(v,desc_a,info)
  call psb_gefree(s,desc_a,info)
  call psb_gefree(t,desc_a,info)
  call psb_gefree(z,desc_a,info)
  call psb_gefree(f,desc_a,info)

  if(psb_errstatus_fatal()) then
    call psb_errpush(info,name)
    goto 9999
  end if

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)
  return

End Subroutine psb_scgstab_vect


text:
      SUBROUTINE DL7TSQ(N, A, L)
C
C  ***  SET A TO LOWER TRIANGLE OF (L**T) * L  ***
C
C  ***  L = N X N LOWER TRIANG. MATRIX STORED ROWWISE.  ***
C  ***  A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L.  ***
C
      INTEGER N
      DOUBLE PRECISION A(1), L(1)
C     DIMENSION A(N*(N+1)/2), L(N*(N+1)/2)
C
      INTEGER I, II, IIM1, I1, J, K, M
      DOUBLE PRECISION LII, LJ
C
      II = 0
      DO 50 I = 1, N
         I1 = II + 1
         II = II + I
         M = 1
         IF (I .EQ. 1) GO TO 30
         IIM1 = II - 1
         DO 20 J = I1, IIM1
              LJ = L(J)
              DO 10 K = I1, J
                   A(M) = A(M) + LJ*L(K)
                   M = M + 1
 10                CONTINUE
 20           CONTINUE
 30      LII = L(II)
         DO 40 J = I1, II
 40           A(J) = LII * L(J)
 50      CONTINUE
C
 999  RETURN
C  ***  LAST CARD OF DL7TSQ FOLLOWS  ***
      END

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
      subroutine caprsh2(irho)
c--
c  Evaluate CAPR, CAPS, and H2.  Here, m is the poloidal mode number
c  and n is the toroidal mode number/periods.
c--LAB--changed calculation of H2 to use proper 1/b**2 Jacobian for
c  flux surface average
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      use parambs
      use vmec0
      implicit none
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      integer irho
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      real(rprec), parameter :: zero = 0, one = 1
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      integer :: nh, m, i, j
      real(rprec) :: h2top_th, h2top_phi, den, qn, e2
      real(rprec) ::  b2, h2top
      real(rprec) :: h2top_tmp, den_tmp, sin_eps
C-----------------------------------------------------------------------

c  the following code was added to evaluate H2 LAB
      h2top = zero
      den   = zero
      do i=1, nthetah
         do j=1, nzetah
            h2top_th = zero
            h2top_phi = zero
            den_tmp = zero
            do nh = 0, nbuse            !  evaluate i,j terms in num and denom of H2
               qn = qsafety(irho)*periods*nh*zetasign
               do m = -mbuse, mbuse
                  sin_eps = amnfit(irho,m,nh)*
     1            (sinmi(m,i)*cosnj(nh,j)+cosmi(m,i)*sinnj(nh,j))
                  h2top_th = h2top_th - m*sin_eps
                  h2top_phi = h2top_phi - qn*sin_eps
                  den_tmp = den_tmp - (m + qn)*sin_eps
               enddo
            enddo
            b2 = bfield(i,j)**2
            h2top = h2top + (h2top_th**2 - h2top_phi**2)/b2
            den  = den + den_tmp**2/b2
         enddo
      enddo

      ! Modified by SAL so code won't stop
      !if (den .eq. zero) stop 'den = 0 in caprsh2'
      IF (den .eq. zero) THEN
         IF (lscreen) THEN
            WRITE(6,'(A)')      'WARNING(BOOTSJ): den = 0 in caprsh2'
            WRITE(6,'(A,I3.3)') '                 irho: ',irho
         END IF
         h2(irho) = 0.0
      ELSE
         h2(irho) = h2top/den
      END IF
      capr(irho) = (one - h2(irho))/(2*qsafety(irho))
      caps(irho) = (one + h2(irho))/2

      end subroutine caprsh2

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
!
! This work (Common Community Physics Package), identified by NOAA, NCAR,
! CU/CIRES, is free of known copyright restrictions and is placed in the
! public domain.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
! IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
! CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!

!>
!! @brief A NO-OP physics modules.
!!
!
module check_noop

    use, intrinsic :: iso_c_binding,                                   &
                      only: c_f_pointer, c_ptr
    use            :: ccpp_types,                                      &
                      only: ccpp_t
    use            :: ccpp_fields,                                     &
                      only: ccpp_field_get
    implicit none

    private
    public :: noop_init_cap, noop_run_cap, noop_finalize_cap

    contains

    subroutine noop_init_cap(ptr) bind(c)
        implicit none
        type(c_ptr), intent(inout) :: ptr

        type(ccpp_t), pointer      :: cdata

        call c_f_pointer(ptr, cdata)

        print *, 'In noop_init_cap'
        print *, cdata%suite%groups(1)%subcycles(1)%schemes(1)%name

    end subroutine noop_init_cap

    subroutine noop_run_cap(ptr) bind(c)
        implicit none
        type(c_ptr), intent(inout) :: ptr

        type(ccpp_t), pointer      :: cdata

        call c_f_pointer(ptr, cdata)

        print *, 'In noop_run_cap'
        print *, cdata%suite%groups(1)%subcycles(1)%schemes(1)%name

    end subroutine noop_run_cap

    subroutine noop_finalize_cap(ptr) bind(c)
        implicit none
        type(c_ptr), intent(inout) :: ptr

        type(ccpp_t), pointer      :: cdata

        call c_f_pointer(ptr, cdata)

        print *, 'In noop_finalize_cap'
        print *, cdata%suite%groups(1)%subcycles(1)%schemes(1)%name

    end subroutine noop_finalize_cap

end module check_noop

max_stars_repo_name: panll/ccpp-framework

max_issues_repo_name: panll/ccpp-framework

max_forks_repo_name: panll/ccpp-framework


text:
C Copyright 1981-2012 ECMWF.
C
C This software is licensed under the terms of the Apache Licence
C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation
C nor does it submit to any jurisdiction.
C
C NIO/NCL
C Modified for the NIO library by David I. Brown
C Allocation is eliminated. All allocation is performed by the calling
C function and allocated space passed in via arguments.
C Removed dependence on grprs.h -- errors printed to standard out (for now)
C Lines removed are left in the code but commented out.
C Added lines are bracketed with
C NIO/NCL
C NIO/NCL END
C
C#ifndef USE_NO_POINTERS
C      SUBROUTINE QU2REG3(PFIELD,KPOINT,KLAT,KLON,KCODE,PMSVAL,KRET,
C     X                   OMISNG,OPERIO,OVEGGY)
      SUBROUTINE QU2REG3(PFIELD,KPOINT,KLAT,KLON,KCODE,PMSVAL,KRET,
     X                   OMISNG,OPERIO,OVEGGY,JPMAX,ZTEMP,ZLINE,ZWORK)
C
C---->
C**** QU2REG3 - Convert quasi-regular grid data to regular.
C
C     Purpose.
C     --------
C
C     Convert quasi-regular grid data to regular,
C     using either a linear or cubic interpolation.
C
C
C**   Interface.
C     ----------
C
C     CALL QU2REG3(PFIELD,KPOINT,KLAT,KLON,KCODE,PMSVAL,OMISNG,OPERIO,
C    X            OVEGGY)
C
C
C     Input Parameters.
C     -----------------
C
C     PFIELD     - Array containing quasi-regular grid
C                  data.
C
C     KPOINT     - Array containing list of the number of
C                  points on each latitude (or longitude) of
C                  the quasi-regular grid.
C
C     KLAT       - Number of latitude lines
C
C     KLON       - Number of longitude lines
C
C     KCODE      - Interpolation required.
C                  1 , linear - data quasi-regular on
C                               latitude lines.
C                  3 , cubic -  data quasi-regular on
C                               latitude lines.
C                  11, linear - data quasi-regular on
C                               longitude lines.
C                  13, cubic -  data quasi-regular on
C                               longitude lines.
C
C     PMSVAL     - Value used for missing data indicator.
C
C     OMISNG     - True if missing values are present in field.
C
C     OPERIO     - True if input field is periodic.
C
C     OVEGGY     - True if 'nearest neighbour' processing must be used
C                  for interpolation
C
C
C     Output Parameters.
C     ------------------
C
C     KRET       - return code
C                  0 = OK
C                  non-zero indicates fatal error
C
C
C     Output Parameters.
C     ------------------
C
C     PFIELD     - Array containing regular grid data.
C
C
C     Method.
C     -------
C
C     Data is interpolated and expanded into a temporary array,
C     which is then copied back into the user's array.
C     Returns an error code if an invalid interpolation is requested
C     or field size exceeds array dimensions.
C
C
C     Externals.
C     ----------
C
C     ROWINA3
C
C
C     Reference.
C     ----------
C
C     WMO Manual on Codes for GRIB code specifications of
C     quasi-regular grids.
C
C
C     Comments.
C     ---------
C
C     This routine is an adaptation of QU2REG to allow missing data
C     values, and hence bit mapped fields.
C
C
C     Author.
C     -------
C
C     J.D.Chambers     ECMWF      22.07.94
C
C
C     Modifications.
C     --------------
C
C     J.D.Chambers     ECMWF      13.09.94
C     Add return code KRET and remove calls to ABORT.
C
C     J.D.Chambers     ECMWF        Feb 1997
C     Allow for 64-bit pointers
C
C     J. Clochard, Meteo France, for ECMWF - January 1998.
C     Addition of OMISNG and OPERIO arguments.
C     Fix message for longitude number out of bounds, and routine
C     name in title and formats.
C
C
C----<
C     -----------------------------------------------------------------
C*    Section 0. Definition of variables. Data statements.
C     -----------------------------------------------------------------
C
      IMPLICIT NONE
C
C#include "grprs.h"
C
C     Parameters
C
      INTEGER JPBYTES, JPMAX
C#ifdef REAL_8
C      PARAMETER (JPBYTES = 8)
C#else
      PARAMETER (JPBYTES = 4)
C#endif
C      PARAMETER (JPMAX=3000)
C     Maximum number of latitudes (or longitudes), for which arrays
C     are dimensioned.
C
C     Subroutine arguments
C
      REAL PFIELD, PMSVAL
      DIMENSION PFIELD(*)
      INTEGER KPOINT, KLAT, KLON, KCODE, KRET
      DIMENSION KPOINT(*)
      REAL ZLINE, ZWORK, ZTEMP
      DIMENSION ZLINE(JPMAX*2)
      DIMENSION ZWORK(0:JPMAX*2+2,3)
      DIMENSION ZTEMP(JPMAX*JPMAX*2)
C
      LOGICAL OMISNG, OPERIO, OVEGGY
C
C     Local variables.
C
      INTEGER ICODE, ILII, ILIO, IQUANO, IREGNO
      INTEGER J210, J220, J225, J230, J240
C
C      REAL ZLINE, ZWORK
C
C#ifndef _CRAYFTN
C#ifdef POINTER_64
C      INTEGER*8 IZTEMP
C#endif
C#endif
C      REAL ZTEMP
C      DIMENSION ZTEMP(1)
C      POINTER ( IZTEMP, ZTEMP )
C      INTEGER ISIZE
C      SAVE ISIZE, IZTEMP
C
C     Externals
C
C#ifdef POINTER_64
C      INTEGER*8 JMALLOC
C#else
C      INTEGER JMALLOC
C#endif
C      EXTERNAL JMALLOC
C
C      DATA ISIZE/0/
C
C     ------------------------------------------------------------------
C*    Section 1. Set initial values.
C     ------------------------------------------------------------------
C
  100 CONTINUE
C
      KRET = 0
C
C     Check input parameters.
C
      IF (KCODE.NE.1.AND.KCODE.NE.3.AND.KCODE.NE.11.AND.KCODE.NE.13)
     X  THEN
C        WRITE(GRPRSM,9001) KCODE
        WRITE(*,9001) KCODE
        KRET = 1
        GOTO 900
      ELSEIF (KLAT.GT.JPMAX) THEN
C        WRITE(GRPRSM,9002) KLAT , JPMAX
        WRITE(*,9002) KLAT , JPMAX
        KRET = 2
        GOTO 900
      ELSEIF (KLON.GT.JPMAX*2) THEN
C        WRITE(GRPRSM,9003) KLON , JPMAX*2
        WRITE(*,9003) KLON , JPMAX*2
        KRET = 3
        GOTO 900
      ENDIF
C
C     Set array indices to 0.
C
      ILII  = 0
      ILIO  = 0
C
C     Establish values of loop parameters.
C
      IF (KCODE.GT.10) THEN
C
C       Quasi-regular along longitude lines.
C
        IQUANO = KLON
        IREGNO = KLAT
        ICODE  = KCODE - 10
      ELSE
C
C       Quasi-regular along latitude lines.
C
        IQUANO = KLAT
        IREGNO = KLON
        ICODE  = KCODE
      ENDIF
C
C     Allocate memory (first time only)
C      IF( ISIZE.EQ.0 ) THEN
C        ISIZE = (JPMAX*JPMAX*2)*JPBYTES
C        IZTEMP = JMALLOC(ISIZE)
C#ifdef hpR64
C        IZTEMP = IZTEMP/(1024*1024*1024*4)
C#endif
C        IF( IZTEMP.EQ.0 ) THEN
C          WRITE(GRPRSM,*) 'QU2REG3: Memory allocation failed.'
C          WRITE(GRPRSM,*) 'QU2REG3: Number of bytes required = ', ISIZE
C          KRET = 5
C          GOTO 900
C        ENDIF
C      ENDIF
C
C     ------------------------------------------------------------------
C*    Section 2. Interpolate field from quasi to regular grid.
C     ------------------------------------------------------------------
C
  200 CONTINUE
C
      DO 230 J230=1,IQUANO
C
        IF (IREGNO.NE.KPOINT(J230)) THEN
C
C         Line contains less values than required,so
C         extract quasi-regular grid values for a line
C
          DO 210 J210=1,KPOINT(J230)
            ILII        = ILII+1
            ZLINE(J210) = PFIELD(ILII)
  210     CONTINUE
C
C         and interpolate this line.
C
          CALL ROWINA3( ZLINE, IREGNO, KPOINT(J230), ZWORK, ICODE,
     X                  PMSVAL, KRET, OMISNG, OPERIO , OVEGGY)
          IF ( KRET .NE. 0 ) GOTO 900
C
C         Add regular grid values for this line to the temporary array.
C
          DO 220 J220=1,IREGNO
            ILIO        = ILIO+1
            ZTEMP(ILIO) = ZLINE(J220)
  220     CONTINUE
C
        ELSE
C
C         Line contains the required number of values, so add
C         this line to the temporary array.
C
          DO 225 J225=1,IREGNO
            ILIO        = ILIO+1
            ILII        = ILII+1
            ZTEMP(ILIO) = PFIELD(ILII)
  225     CONTINUE
C
        ENDIF
C
  230 CONTINUE
C
C     Copy temporary array to user array.
C
      DO 240 J240=1,KLON*KLAT
        PFIELD(J240) = ZTEMP(J240)
  240 CONTINUE
C
C     ------------------------------------------------------------------
C*    Section 9. Return to calling routine. Format statements.
C     ------------------------------------------------------------------
C
  900 CONTINUE
C
      RETURN
C
 9001 FORMAT (1H ,'QU2REG3 : Invalid interpolation type code = ',I3)
C
 9002 FORMAT (1H ,'QU2REG3 : Number of latitudes is ',I4,', maximum ',
     C                      'allowed is ',I4,'.')
C
 9003 FORMAT (1H ,'QU2REG3 : Number of longitudes is ',I4,', maximum ',
     C                      'allowed is ',I4,'.')
C
      END
C#endif

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


!
!  Copyright 2016 ARTED developers
!
!  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.
!
subroutine PSE_hpsi(ik)
  use global_variables
  use PSE_variables
  implicit none
  integer :: ik

  select case(type_spatial_difference)
  case('FD')
    call PSE_hpsi_grid(ik)
  case('FT')
    call PSE_hpsi_DFT(ik)
  case default
    err_message = 'invalid parameter in type_spatial_difference'
    call err_finalize
  end select

  return
end subroutine PSE_hpsi


text:
*DECK ZRATI
      SUBROUTINE ZRATI (ZR, ZI, FNU, N, CYR, CYI, TOL)
C***BEGIN PROLOGUE  ZRATI
C***SUBSIDIARY
C***PURPOSE  Subsidiary to ZBESH, ZBESI and ZBESK
C***LIBRARY   SLATEC
C***TYPE      ALL (CRATI-A, ZRATI-A)
C***AUTHOR  Amos, D. E., (SNL)
C***DESCRIPTION
C
C     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
C     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
C     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
C     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
C     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
C     BY D. J. SOOKNE.
C
C***SEE ALSO  ZBESH, ZBESI, ZBESK
C***ROUTINES CALLED  ZABS, ZDIV
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   910415  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  ZRATI
      DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR,
     * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU,
     * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI,
     * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS
      INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
      DIMENSION CYR(N), CYI(N)
      EXTERNAL ZABS
      DATA CZEROR,CZEROI,CONER,CONEI,RT2/
     1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 /
C***FIRST EXECUTABLE STATEMENT  ZRATI
      AZ = ZABS(ZR,ZI)
      INU = FNU
      IDNU = INU + N - 1
      MAGZ = AZ
      AMAGZ = MAGZ+1
      FDNU = IDNU
      FNUP = MAX(AMAGZ,FDNU)
      ID = IDNU - MAGZ - 1
      ITIME = 1
      K = 1
      PTR = 1.0D0/AZ
      RZR = PTR*(ZR+ZR)*PTR
      RZI = -PTR*(ZI+ZI)*PTR
      T1R = RZR*FNUP
      T1I = RZI*FNUP
      P2R = -T1R
      P2I = -T1I
      P1R = CONER
      P1I = CONEI
      T1R = T1R + RZR
      T1I = T1I + RZI
      IF (ID.GT.0) ID = 0
      AP2 = ZABS(P2R,P2I)
      AP1 = ZABS(P1R,P1I)
C-----------------------------------------------------------------------
C     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
C     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
C     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
C     PREMATURELY.
C-----------------------------------------------------------------------
      ARG = (AP2+AP2)/(AP1*TOL)
      TEST1 = SQRT(ARG)
      TEST = TEST1
      RAP1 = 1.0D0/AP1
      P1R = P1R*RAP1
      P1I = P1I*RAP1
      P2R = P2R*RAP1
      P2I = P2I*RAP1
      AP2 = AP2*RAP1
   10 CONTINUE
      K = K + 1
      AP1 = AP2
      PTR = P2R
      PTI = P2I
      P2R = P1R - (T1R*PTR-T1I*PTI)
      P2I = P1I - (T1R*PTI+T1I*PTR)
      P1R = PTR
      P1I = PTI
      T1R = T1R + RZR
      T1I = T1I + RZI
      AP2 = ZABS(P2R,P2I)
      IF (AP1.LE.TEST) GO TO 10
      IF (ITIME.EQ.2) GO TO 20
      AK = ZABS(T1R,T1I)*0.5D0
      FLAM = AK + SQRT(AK*AK-1.0D0)
      RHO = MIN(AP2/AP1,FLAM)
      TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0D0))
      ITIME = 2
      GO TO 10
   20 CONTINUE
      KK = K + 1 - ID
      AK = KK
      T1R = AK
      T1I = CZEROI
      DFNU = FNU + (N-1)
      P1R = 1.0D0/AP2
      P1I = CZEROI
      P2R = CZEROR
      P2I = CZEROI
      DO 30 I=1,KK
        PTR = P1R
        PTI = P1I
        RAP1 = DFNU + T1R
        TTR = RZR*RAP1
        TTI = RZI*RAP1
        P1R = (PTR*TTR-PTI*TTI) + P2R
        P1I = (PTR*TTI+PTI*TTR) + P2I
        P2R = PTR
        P2I = PTI
        T1R = T1R - CONER
   30 CONTINUE
      IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40
      P1R = TOL
      P1I = TOL
   40 CONTINUE
      CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N))
      IF (N.EQ.1) RETURN
      K = N - 1
      AK = K
      T1R = AK
      T1I = CZEROI
      CDFNUR = FNU*RZR
      CDFNUI = FNU*RZI
      DO 60 I=2,N
        PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1)
        PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1)
        AK = ZABS(PTR,PTI)
        IF (AK.NE.CZEROR) GO TO 50
        PTR = TOL
        PTI = TOL
        AK = TOL*RT2
   50   CONTINUE
        RAK = CONER/AK
        CYR(K) = RAK*PTR*RAK
        CYI(K) = -RAK*PTI*RAK
        T1R = T1R - CONER
        K = K - 1
   60 CONTINUE
      RETURN
      END


text:
! ======================================================================
! NAME
!
!     io [ io ] - Input-Output Module
!
! SYNOPSIS
!
!     USE :: IO
!
! DESCRIPTION
!
!     IO is a custom Fortran module written to perform input-output
!     operations for gpse.
!
! OPTIONS
!
! SEE ALSO
!
! BUGS
!
! HISTORY
!
! AUTHOR
!
!     Marty Kandes, Ph.D.
!     Computational & Data Science Research Specialist
!     High-Performance Computing User Services Group
!     San Diego Supercomputer Center
!     University of California, San Diego
!
! COPYRIGHT
!
!     Copyright (c) 2014, 2015, 2016, 2017, 2018, 2019, 2020 Martin Charles Kandes
!
! LAST UPDATED
!
!     Monday, May 13th, 2019
!
! ----------------------------------------------------------------------

      MODULE IO

! --- MODULE DECLARATIONS ----------------------------------------------

      USE, INTRINSIC :: ISO_FORTRAN_ENV

! --- MODULE DEFINITIONS -----------------------------------------------
!
!     ISO_FORTRAN_ENV is the intrinsic Fortran module that provides
!        information about the run-time environment.
!
! ----------------------------------------------------------------------

      IMPLICIT NONE
      PRIVATE

! --- SUBROUTINE DECLARATIONS ------------------------------------------

      PUBLIC :: io_read_bin_psi
      PUBLIC :: io_read_bin_vex
      PUBLIC :: io_write_bin_psi
      PUBLIC :: io_write_bin_vex
      PUBLIC :: io_write_splot_psi2
      PUBLIC :: io_write_vtk_header
      PUBLIC :: io_write_vtk_xcoordinates
      PUBLIC :: io_write_vtk_ycoordinates
      PUBLIC :: io_write_vtk_zcoordinates
      PUBLIC :: io_write_vtk_repsi
      PUBLIC :: io_write_vtk_impsi
      PUBLIC :: io_read_vtk
      PUBLIC :: io_write_vtk

! --- SUBROUTINE DEFINITIONS -------------------------------------------
!
!     io_read_bin_psi is a PUBLIC SUBROUTINE that performs unformatted
!        stream I/O to read a wave function from a binary file.
!
!     io_read_bin_vex is a PUBLIC SUBROUTINE that performs unformatted
!        stream I/O to read an external potential from an binary file.
!
!     io_write_bin_psi is a PUBLIC SUBROUTINE that performs unformatted
!         stream I/O to write a wave function to a binary file.
!
!     io_write_bin_vex is a PUBLIC SUBROUTINE that performs unformatted
!        stream I/O to write an external potential to a binary file.
!
!     io_write_vtk_header is a PUBLIC SUBROUTINE that performs formatted
!        stream I/O to write the standard legacy header of a vtk file
!        for an output wave function.
!
!     io_write_vtk_xcoordinates is a PUBLIC SUBROUTINE that performs
!        formatted stream I/O to write the x-axis coordinates of the
!        grid points to a vtk file.
!
!     io_write_vtk_ycoordinates is a PUBLIC SUBROUTINE that performs
!        formatted stream I/O to write the y-axis coordinates of the
!        grid points to a vtk file.
!
!     io_write_vtk_zcoordinates is a PUBLIC SUBROUTINE that performs
!        formatted stream I/O to write the z-axis coordinates of the
!        grid points to a vtk file.
!
!     io_write_vtk_repsi is a PUBLIC SUBROUTINE that performs formatted
!        stream I/O to write the real component of an output wave
!        function to a vtk file.
!
!     io_write_vtk_impsi is a PUBLIC SUBROUTINE that performs formatted
!        stream I/O to write the imaginary component of an output wave
!        function to a vtk file.
!
!     io_read_vtk is a PUBLIC SUBROUTINE that performs formatted
!        record-based I/O to read in a wave function from a vtk file,
!        which is primarily used to implement simple parallel input.
!
!     io_write_vtk is a PUBLIC SUBROUTINE that performs formatted
!        record-based I/O to write a wave function to a vtk file, which
!        is primarily used to implement simple parallel output.
!
! ----------------------------------------------------------------------

      CONTAINS

! ----------------------------------------------------------------------

      SUBROUTINE io_read_bin_psi(fileUnit, filePosition, nXa, nXb, &
         & nXbc, nYa, nYb, nYbc, nZa, nZb, nZbc, Psi3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM('psi-'//fileUnitChar//'.bin'), &
         & ACCESS='STREAM', ACTION='READ', FORM='UNFORMATTED', &
         & STATUS='UNKNOWN')

         DO l = nZa, nZb
            DO k = nYa, nYb
               DO j = nXa, nXb
                  READ(UNIT=fileUnit, POS=filePosition) Psi3(j,k,l)
                  INQUIRE(UNIT=fileUnit, POS=filePosition)
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_read_bin_vex(fileUnit, filePosition, nXa, nXb, &
         & nXbc, nYa, nYb, nYbc, nZa, nZb, nZbc, Vex3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(INOUT) :: Vex3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM('psi-'//fileUnitChar//'.bin'), &
         & ACCESS='STREAM', ACTION='READ', FORM='UNFORMATTED', &
         & STATUS='UNKNOWN')

         DO l = nZa, nZb
            DO k = nYa, nYb
               DO j = nXa, nXb
                  READ(UNIT=fileUnit, POS=filePosition) Vex3(j,k,l)
                  INQUIRE(UNIT=fileUnit, POS=filePosition)
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_bin_psi(fileUnit, filePosition, nXa, nXb, &
         & nXbc, nYa, nYb, nYbc, nZa, nZb, nZbc, Psi3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM('psi-'//fileUnitChar//'.bin'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='UNFORMATTED', &
         & STATUS='UNKNOWN')

         DO l = nZa, nZb
            DO k = nYa, nYb
               DO j = nXa, nXb
                  WRITE(UNIT=fileUnit, POS=filePosition) Psi3(j,k,l)
                  INQUIRE(UNIT=fileUnit, POS=filePosition)
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_bin_vex(fileUnit, filePosition, nXa, nXb, &
         & nXbc, nYa, nYb, nYbc, nZa, nZb, nZbc, Vex3)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc
      INTEGER, INTENT(IN) :: nZa, nZb, nZbc

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Vex3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM('vex-'//fileUnitChar//'.bin'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='UNFORMATTED', &
         & STATUS='UNKNOWN')

         DO l = nZa, nZb
            DO k = nYa, nYb
               DO j = nXa, nXb
                  WRITE(UNIT=fileUnit, POS=filePosition) Vex3(j,k,l)
                  INQUIRE(UNIT=fileUnit, POS=filePosition)
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_splot_psi2(fileName, fileUnit, nXa, nXb, &
         & nXbc, nYa, nYb, nYbc, X, Y, Psi2)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER, INTENT(IN) :: nXa, nXb, nXbc
      INTEGER, INTENT(IN) :: nYa, nYb, nYbc

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: X
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Y

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc), INTENT (IN) :: Psi2

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.splot'), &
         & ACCESS='SEQUENTIAL', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         DO k = nYa, nYb
            DO j = nXa, nXb
               WRITE(UNIT=fileUnit, FMT=*) Y(k), X(j), Psi2(j,k)
            END DO
            WRITE(UNIT=fileUnit , FMT=*) ' '
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_header(fileName, fileUnit, filePosition, &
         & nX, nY, nZ, dNx, dNy, dNz)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nX, nY, nZ
      INTEGER, INTENT(IN) :: dNx, dNy, dNz

      CHARACTER(LEN=4) :: fileUnitChar

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         WRITE(UNIT=fileUnit, FMT='(A26)') '# vtk DataFile Version 3.0' ! VTK file identifier and version number
         WRITE(UNIT=fileUnit, FMT='(A26)') 'STANDARD LEGACY VTK FORMAT' ! VTK file header; 256 characters maximum
         WRITE(UNIT=fileUnit, FMT='(A5)')  'ASCII'                      ! VTK file format: ASCII or BINARY
         WRITE(UNIT=fileUnit, FMT='(A24)') 'DATASET RECTILINEAR_GRID'   ! VTK dataset structure
         WRITE(UNIT=fileUnit, FMT='(A10,1X,I4.1,1X,I4.1,1X,I4.1)') &
            & 'DIMENSIONS', nX / dNx + MODULO(nX,dNx), &
            & nY / dNy + MODULO(nY,dNy), nZ / dNz + MODULO(nZ,dNz)
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE
! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_xcoordinates(fileName, fileUnit, &
         & filePosition, nX, nXa, nXb, nXbc, dNx, X)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nX, nXa, nXb, nXbc, dNx

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: X

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         WRITE(UNIT=fileUnit, POS=filePosition, &
            & FMT='(A13,1X,I4.1,1X,A6)' ) 'X_COORDINATES', &
            & nX / dNx + MODULO(nX,dNx) , 'double'
         DO j = nXa, nXb, dNx
            WRITE(UNIT=fileUnit, FMT=*) X(j)
         END DO
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_ycoordinates(fileName, fileUnit, &
         & filePosition, nY, nYa, nYb, nYbc, dNy, Y)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: nY, nYa, nYb, nYbc, dNy

      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Y

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: k

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         WRITE(UNIT=fileUnit, POS=filePosition, &
            & FMT='(A13,1X,I4.1,1X,A6)') 'Y_COORDINATES', &
            & nY / dNy + MODULO(nY,dNy), 'double'
         DO k = nYa, nYb, dNy
            WRITE(UNIT=fileUnit, FMT=*) Y(k)
         END DO
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_zcoordinates(fileName, fileUnit, &
         & filePosition, mpiSource, nZ, nZa, nZb, nZbc, dNz, Z)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: mpiSource
      INTEGER, INTENT(IN) :: nZ, nZa, nZb, nZbc, dNz

      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Z

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         IF (mpiSource == 0) THEN ! i.e., if mpiSource == MPI_MASTER

            WRITE(UNIT=fileUnit, POS=filePosition, &
               & FMT='(A13,1X,I4.1,1X,A6)') 'Z_COORDINATES', &
               & nZ / dNz + MODULO(nZ,dNz), 'double'
            DO l = nZa, nZb, dNz
               WRITE(UNIT=fileUnit, FMT=*) Z(l)
            END DO

         ELSE

            WRITE(UNIT=fileUnit, POS=filePosition, FMT='(F23.15)') Z(nZa)
            DO l = nZa + dNz, nZb, dNz
               WRITE(UNIT=fileUnit, FMT=*) Z(l)
            END DO

         END IF
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_repsi(fileName, fileUnit, filePosition, &
         & mpiSource, nX, nXa, nXb, nXbc, dNx, nY, nYa, nYb, nYbc, dNy,&
         & nZ, nZa, nZb, nZbc, dNz, Psi3)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: mpiSource
      INTEGER, INTENT(IN) :: nX, nXa, nXb, nXbc, dNx
      INTEGER, INTENT(IN) :: nY, nYa, nYb, nYbc, dNy
      INTEGER, INTENT(IN) :: nZ, nZa, nZb, nZbc, dNz

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         IF (mpiSource == 0) THEN

            WRITE(UNIT=fileUnit, POS=filePosition, &
               & FMT='(A10,1X,I19.1)') 'POINT_DATA', &
               & (nX / dNx + MODULO(nX, dNx)) &
               & * (nY / dNy + MODULO(nY, dNy)) &
               & * (nZ / dNz + MODULO(nZ, dNz))
            WRITE(UNIT=fileUnit, FMT='(A24)') 'SCALARS RePsi double 1'
            WRITE(UNIT=fileUnit, FMT='(A20)') 'LOOKUP_TABLE default'
            DO l = nZa, nZb, dNz
               DO k = nYa, nYb, dNy
                  DO j = nXa, nXb, dNx
                     WRITE(UNIT=fileUnit, FMT=*) REAL(Psi3(j,k,l))
                  END DO
               END DO
            END DO

         ELSE

            DO l = nZa, nZb, dNz
               DO k = nYa, nYb, dNy
                  DO j = nXa, nXb, dNx
                     WRITE(UNIT=fileUnit, POS=filePosition, FMT=*) &
                        & REAL(Psi3(j,k,l))
                     INQUIRE(UNIT=fileUnit, POS=filePosition)
                  END DO
               END DO
            END DO

         END IF
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk_impsi(fileName, fileUnit, filePosition, &
         & mpiSource, nX, nXa, nXb, nXbc, dNx, nY, nYa, nYb, nYbc, dNy,&
         & nZ, nZa, nZb, nZbc, dNz, Psi3)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER(KIND=8), INTENT(INOUT) :: filePosition ! hard coded KIND = 8 to avoid integer overflow when writing large files
      INTEGER, INTENT(IN) :: mpiSource
      INTEGER, INTENT(IN) :: nX, nXa, nXb, nXbc, dNx
      INTEGER, INTENT(IN) :: nY, nYa, nYb, nYbc, dNy
      INTEGER, INTENT(IN) :: nZ, nZa, nZb, nZbc, dNz

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'.vtk'), &
         & ACCESS='STREAM', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='UNKNOWN')

         IF (mpiSource == 0) THEN

            WRITE(UNIT=fileUnit, POS=filePosition, FMT='(A24)') &
               & 'SCALARS ImPsi double 1'
            WRITE(UNIT=fileUnit, FMT='(A20)') 'LOOKUP_TABLE default'
            DO l = nZa, nZb, dNz
               DO k = nYa, nYb, dNy
                  DO j = nXa, nXb, dNx
                     WRITE(UNIT=fileUnit, FMT=*) AIMAG(Psi3(j,k,l))
                  END DO
               END DO
            END DO

         ELSE

            DO l = nZa, nZb, dNz
               DO k = nYa, nYb, dNy
                  DO j = nXa, nXb, dNx
                     WRITE(UNIT=fileUnit, POS=filePosition, FMT=*) &
                        & AIMAG(Psi3(j,k,l))
                     INQUIRE(UNIT=fileUnit, POS=filePosition)
                  END DO
               END DO
            END DO

         END IF
         INQUIRE(UNIT=fileUnit, POS=filePosition)

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_read_vtk(fileName, fileUnit, mpiRank, nX, nXa, &
         & nXb, nXbc, dNx, nY, nYa, nYb, nYbc, dNy, nZ, nZa, nZb, nZbc,&
         & dNz, X, Y, Z, Psi3)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER, INTENT(IN) :: mpiRank
      INTEGER, INTENT(IN) :: nX, nXa, nXb, nXbc, dNx
      INTEGER, INTENT(IN) :: nY, nYa, nYb, nYbc, dNy
      INTEGER, INTENT(IN) :: nZ, nZa, nZb, nZbc, dNz

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(OUT) :: X
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(OUT) :: Y
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(OUT) :: Z

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(OUT) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar
      CHARACTER(LEN=4) :: mpiRankChar

      CHARACTER(LEN=26) :: vtkFileIdentifierVersionNumber
      CHARACTER(LEN=26) :: vtkFileHeader
      CHARACTER(LEN=5)  :: vtkFileFormat
      CHARACTER(LEN=24) :: vtkDatasetStructure
      CHARACTER(LEN=10) :: vtkDimensions
      CHARACTER(LEN=13) :: vtkXCoordinates
      CHARACTER(LEN=6)  :: vtkXCoordinatesType
      CHARACTER(LEN=13) :: vtkYCoordinates
      CHARACTER(LEN=6)  :: vtkYCoordinatesType
      CHARACTER(LEN=13) :: vtkZCoordinates
      CHARACTER(LEN=6)  :: vtkZCoordinatesType
      CHARACTER(LEN=10) :: vtkPointData
      CHARACTER(LEN=24) :: vtkScalarsRePsi
      CHARACTER(LEN=24) :: vtkScalarsImPsi
      CHARACTER(LEN=20) :: vtkLookUpTable

      INTEGER :: j, k, l

      INTEGER :: nXvtk, nYvtk, nZvtk, nXnYnZvtk

      REAL, ALLOCATABLE, DIMENSION(:, :, :) :: RePsi3
      REAL, ALLOCATABLE, DIMENSION(:, :, :) :: ImPsi3

      ALLOCATE(RePsi3(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc))
      ALLOCATE(ImPsi3(nXa - nXbc : nXb + nXbc, &
                    & nYa - nYbc : nYb + nYbc, &
                    & nZa - nZbc : nZb + nZbc))

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit
      WRITE(UNIT=mpiRankChar, FMT='(I4.4)') mpiRank

      OPEN(UNIT=fileUnit, FILE=TRIM(fileName//fileUnitChar//'-'//&
         &mpiRankChar//'.vtk'), ACCESS='SEQUENTIAL', ACTION='READ', &
         & FORM='FORMATTED', STATUS='UNKNOWN')

         READ(UNIT=fileUnit, FMT='(A26)') vtkFileIdentifierVersionNumber
         READ(UNIT=fileUnit, FMT='(A26)') vtkFileHeader
         READ(UNIT=fileUnit, FMT='(A5)') vtkFileFormat
         READ(UNIT=fileUnit, FMT='(A24)') vtkDatasetStructure
         READ(UNIT=fileUnit, FMT='(A10,1X,I4.1,1X,I4.1,1X,I4.1)') &
            & vtkDimensions, nXvtk, nYvtk, nZvtk

         READ(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
            & vtkXCoordinates, nXvtk, vtkXCoordinatesType
         DO j = nXa, nXb, dNx
            READ(UNIT=fileUnit, FMT=*) X(j)
         END DO

         READ(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
            & vtkYCoordinates, nYvtk, vtkYCoordinatesType
         DO k = nYa, nYb, dNy
            READ(UNIT=fileUnit, FMT=*) Y(k)
         END DO

         READ(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
            & vtkZCoordinates, nZvtk, vtkZCoordinatesType
         DO l = nZa, nZb, dNz
            READ(UNIT=fileUnit, FMT=*) Z(l)
         END DO

         Psi3 = CMPLX(0.0, 0.0)

         READ(UNIT=fileUnit, FMT='(A10,1X,I19.1)') vtkPointData, &
            & nXnYnZvtk
         READ(UNIT=fileUnit, FMT='(A24)') vtkScalarsRePsi
         READ(UNIT=fileUnit, FMT='(A20)') vtkLookUpTable
         DO l = nZa, nZb, dNz
            DO k = nYa, nYb, dNy
               DO j = nXa, nXb, dNx
                  READ(UNIT=fileUnit, FMT=*) RePsi3(j,k,l)
               END DO
            END DO
         END DO

         READ(UNIT=fileUnit, FMT='(A24)') vtkScalarsImPsi
         READ(UNIT=fileUnit, FMT='(A20)') vtkLookUpTable
         DO l = nZa, nZb, dNz
            DO k = nYa, nYb, dNy
               DO j = nXa, nXb, dNx
                  READ(UNIT=fileUnit, FMT=*) ImPsi3(j,k,l)
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      DO l = nZa, nZb, dNz
         DO k = nYa, nYb, dNy
            DO j = nXa, nXb, dNx
               Psi3(j,k,l) = CMPLX(RePsi3(j,k,l), ImPsi3(j,k,l))
            END DO
         END DO
      END DO

      DEALLOCATE(ImPsi3)
      DEALLOCATE(RePsi3)

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      SUBROUTINE io_write_vtk(fileName, fileUnit, mpiRank, nX, nXa, &
         & nXb, nXbc, dNx, nY, nYa, nYb, nYbc, dNy, nZ, nZa, nZb, nZbc,&
         & dNz, X, Y, Z, Psi3)
      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN) :: fileName

      INTEGER, INTENT(IN) :: fileUnit
      INTEGER, INTENT(IN) :: mpiRank
      INTEGER, INTENT(IN) :: nX, nXa, nXb, nXbc, dNx
      INTEGER, INTENT(IN) :: nY, nYa, nYb, nYbc, dNy
      INTEGER, INTENT(IN) :: nZ, nZa, nZb, nZbc, dNz

      REAL, DIMENSION(nXa - nXbc : nXb + nXbc), INTENT(IN) :: X
      REAL, DIMENSION(nYa - nYbc : nYb + nYbc), INTENT(IN) :: Y
      REAL, DIMENSION(nZa - nZbc : nZb + nZbc), INTENT(IN) :: Z

      COMPLEX, DIMENSION(nXa - nXbc : nXb + nXbc, &
                       & nYa - nYbc : nYb + nYbc, &
                       & nZa - nZbc : nZb + nZbc), INTENT(IN) :: Psi3

      CHARACTER(LEN=4) :: fileUnitChar
      CHARACTER(LEN=4) :: mpiRankChar

      INTEGER :: j, k, l

      WRITE(UNIT=fileUnitChar, FMT='(I4.4)') fileUnit
      WRITE(UNIT=mpiRankChar , FMT='(I4.4)') mpiRank

      OPEN(UNIT=fileUnit, &
         & FILE=TRIM(fileName//fileUnitChar//'-'//mpiRankChar//'.vtk'),&
         & ACCESS='SEQUENTIAL', ACTION='WRITE', FORM='FORMATTED', &
         & STATUS='NEW')

         WRITE(UNIT=fileUnit, FMT='(A26)') '# vtk DataFile Version 3.0' ! VTK file identifier and version number
         WRITE(UNIT=fileUnit, FMT='(A26)') 'STANDARD LEGACY VTK FORMAT' ! VTK file header; 256 characters maximum
         WRITE(UNIT=fileUnit, FMT='(A5)')  'ASCII'                      ! VTK file format: ASCII or BINARY
         WRITE(UNIT=fileUnit, FMT='(A24)') 'DATASET RECTILINEAR_GRID'   ! VTK dataset structure
         WRITE(UNIT=fileUnit, FMT='(A10,1X,I4.1,1X,I4.1,1X,I4.1)') &
            & 'DIMENSIONS', nX / dNx + MODULO(nX, dNx), &
            & nY / dNy + MODULO(nY, dNy), &
            & (nZb - nZa + 1) / dNz + MODULO(nZ , dNz)

         WRITE(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
            & 'X_COORDINATES', nX / dNx + MODULO(nX, dNx), 'double'
         DO j = nXa, nXb, dNx
            WRITE(UNIT=fileUnit, FMT=*) X(j)
         END DO

         WRITE(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
           & 'Y_COORDINATES', nY / dNy + MODULO(nY, dNy), 'double'
         DO k = nYa, nYb, dNy
            WRITE(UNIT=fileUnit, FMT=*) Y(k)
         END DO

         WRITE(UNIT=fileUnit, FMT='(A13,1X,I4.1,1X,A6)') &
           & 'Z_COORDINATES', &
           & (nZb - nZa + 1) / dNz + MODULO(nZb - nZb + 1, dNz), &
           & 'double'
         DO l = nZa, nZb, dNz
            WRITE(UNIT=fileUnit, FMT=*) Z(l)
         END DO

         WRITE(UNIT=fileUnit, FMT='(A10,1X,I19.1)') 'POINT_DATA', &
            & (nX / dNx + MODULO(nX, dNx)) &
            & * (nY / dNy + MODULO(nY, dNy)) &
            & * ((nZb - nZa + 1) / dNz + MODULO(nZb - nZa + 1, dNz))
         WRITE(UNIT=fileUnit, FMT='(A24)') 'SCALARS RePsi double 1'
         WRITE(UNIT=fileUnit, FMT='(A20)') 'LOOKUP_TABLE default'
         DO l = nZa, nZb, dNz
            DO k = nYa, nYb, dNy
               DO j = nXa, nXb, dNx
                  WRITE(UNIT=fileUnit, FMT=*) REAL(Psi3(j,k,l))
               END DO
            END DO
         END DO

         WRITE(UNIT=fileUnit, FMT='(A24)') 'SCALARS ImPsi double 1'
         WRITE(UNIT=fileUnit, FMT='(A20)') 'LOOKUP_TABLE default'
         DO l = nZa, nZb, dNz
            DO k = nYa, nYb, dNy
               DO j = nXa, nXb, dNx
                  WRITE(UNIT=fileUnit, FMT=*) AIMAG(Psi3(j,k,l))
               END DO
            END DO
         END DO

      CLOSE(UNIT=fileUnit, STATUS='KEEP')

      RETURN
      END SUBROUTINE

! ----------------------------------------------------------------------

      END MODULE

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

max_stars_repo_name: mkandes/gpse

max_issues_repo_name: mkandes/gpse

max_forks_repo_name: mkandes/gpse