pfarm/cpu/src_mpi_omp_atom/DEFS

+DEFS_intel_skl
\ No newline at end of file

pfarm/cpu/src_mpi_omp_atom/DEFS_cray

+LIBS = 
+STATIC_LIBS = -L../../lib -lfxdr
+
+# Default Archer2 environment (Cray) 
+F90 = ftn
+
+# Skylake options
+#F90FLAGS = -mtune=skylake -mkl=parallel -O3
+
+# AMD options
+F90FLAGS = -O3
+
+F90FLAGS1 = $(F90FLAGS)
+LDFLAGS=$(F90FLAGS)

pfarm/cpu/src_mpi_omp_atom/DEFS_intel_skl

+LIBS = 
+STATIC_LIBS = -L../../lib -lfxdr -ltirpc
+
+F90 = mpif90
+
+# KNL options
+#F90FLAGS = -mtune=knl -mkl=parallel -O3
+
+# Skylake options
+F90FLAGS = -mtune=skylake -mkl=parallel -Ofast
+
+F90FLAGS1 = $(F90FLAGS)
+LDFLAGS=$(F90FLAGS)

pfarm/cpu/src_mpi_omp_atom/Makefile

+PROG = ../bin/exdig_mpi_omp_atom
+
+include DEFS
+
+XDR = xdr_files
+
+SRCS =	amp.f90 def_sectors_bench.f90 pdg_ctl.f90 dist_mat.f90 \
+	energy_grid.f90 rmx1.f90 potl.f90 scaling.f90 \
+	potl_coefs.f90 precisn.f90 rd_Hfile.f90 rmx1_in.f90 \
+	sec_ham.f90 slp.f90 error_out.f90 io_units.f90 mpi_params.f90
+SRC1 = $(XDR).f90
+
+
+OBJS = $(SRCS:.f90=.o)
+OBJ1 = $(SRC1:.f90=.o)
+
+
+all: $(PROG)
+
+$(PROG): $(OBJ1) $(OBJS) $(ULIBS)
+	$(F90) $(LDFLAGS) -o $@ $(OBJS) $(OBJ1) $(STATIC_LIBS) $(LIBS)
+
+$(OBJ1): $(SRC1)
+	$(F90) $(F90FLAGS) -c  $<
+clean:
+	rm -f $(PROG) $(OBJS) $(OBJ1) *.mod *~
+very_clean:
+	rm -f $(PROG) $(OBJS) $(OBJ1) *.mod *~
+
+.SUFFIXES: $(SUFFIXES) .f90
+
+.f90.o:
+	$(F90) $(F90FLAGS) -c $<
+
+amp.o: amp.f90 def_sectors_bench.o dist_mat.o energy_grid.o \
+	io_units.o precisn.o slp.o $(XDR).o error_out.o mpi_params.o
+def_sectors_bench.o: def_sectors_bench.f90 energy_grid.o io_units.o precisn.o \
+	rd_Hfile.o rmx1_in.o error_out.o mpi_params.o
+pdg_ctl.o: pdg_ctl.f90 amp.o def_sectors_bench.o dist_mat.o io_units.o potl.o \
+	potl_coefs.o precisn.o rd_Hfile.o rmx1_in.o slp.o mpi_params.o
+dist_mat.o: dist_mat.f90 io_units.o precisn.o sec_ham.o error_out.o
+energy_grid.o: energy_grid.f90 io_units.o precisn.o rd_Hfile.o \
+	rmx1_in.o error_out.o mpi_params.o
+error_prt.o: error_prt.f90 io_units.o
+rmx1.o: rmx1.f90 pdg_ctl.o energy_grid.o io_units.o precisn.o rd_Hfile.o \
+	rmx1_in.o slp.o error_out.o
+potl.o: potl.f90 io_units.o potl_coefs.o precisn.o rd_Hfile.o rmx1_in.o
+potl_coefs.o: potl_coefs.f90 io_units.o precisn.o scaling.o slp.o $(XDR).o
+rd_Hfile.o: rd_Hfile.f90 io_units.o precisn.o rmx1_in.o scaling.o $(XDR).o error_out.o mpi_params.o
+rmx1_in.o: rmx1_in.f90 io_units.o precisn.o error_out.o mpi_params.o
+scaling.o: scaling.f90 io_units.o precisn.o
+sec_ham.o: sec_ham.f90 precisn.o io_units.o potl.o error_out.o
+$(XDR).o: $(XDR).f90 io_units.o precisn.o
+slp.o: slp.f90
+io_units.o: io_units.f90
+precisn.o: precisn.f90
+error_out.o: error_out.f90 io_units.o

pfarm/cpu/src_mpi_omp_atom/a_descriptor.f90

+module array_descriptor
+! define parameters used to descibe BLACS distributed arrays:
+! Time-stamp: "03/07/04 14:31:08 cjn"
+  implicit none
+
+  public
+
+  integer, parameter :: dlen_ = 9
+  integer, parameter :: block_cyclic_2d = 1
+  integer, parameter :: ctxt_ = 2
+  integer, parameter :: m_ = 3
+  integer, parameter :: n_ = 4
+  integer, parameter :: mb_ = 5
+  integer, parameter :: nb_ = 6
+  integer, parameter :: rsrc_ = 7
+  integer, parameter :: csrc_ = 8
+  integer, parameter :: lld_ = 9
+end module array_descriptor

pfarm/cpu/src_mpi_omp_atom/amp.f90

+module rw_amplitudes
+! calculated sector reduced-width amplitudes
+! Time-stamp: "2003-10-06 14:23:21 cjn"
+  use io_units, only: fo, fa
+  use precisn, only: wp
+  use def_sectors, only: nh, nsect
+  use mpi_params, only: taskid, numtasks, io_processor
+  implicit none
+
+  integer, save     :: ixdro   ! XDR unit number for rwa file
+
+  private
+  public start_sa_file, close_sa_file, ampltd
+
+contains
+
+
+  subroutine start_sa_file (sqno, st, sect, imesh)
+! open XDR file to hold sector reduced width amplitudes
+    use xdr_files, only: open_xdr, xdr_io
+    use energy_grid, only: fine_mesh
+    use slp, only: qno, val_qno
+    use rmx1_in, only: xdr, split_prop
+    type(qno), intent(in)   :: sqno        ! scattering q#s
+    integer, intent(in)     :: st         ! second target spin
+    integer, intent(in)     :: sect ! sector id
+    integer, intent(in)     :: imesh ! mesh id    
+    character(len=13)       :: hname
+    character(len=2)        :: sp, l
+    character(len=1)        :: p
+    character(len=4)        :: stem
+    character(len=4)        :: isect
+    integer                 :: ios, qn(3), status
+    integer                 :: nspn ! scattering spin.  
+
+! File naming convention:
+!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
+!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
+!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
+!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
+! For no-exchange LS files target spin multiplicity is used
+! For J files L is replaced by 2*J and total spin multiplicity by 0 
+!                                  or target spin multiplicity by 2*K
+! For no-exchange files L is replaced by 2*K 
+!                   and total spin multiplicity by 0
+
+    ixdro = fa
+    stem = 'AMPC'
+    if (fine_mesh) stem = 'AMPF'
+    qn = val_qno(sqno)
+    if (split_prop) then   ! channel splitting in use
+       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
+    else
+       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
+       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
+    end if
+    write (l,'(i2.2)') abs(qn(1))
+    write (p, '(i1.1)') qn(3)
+    write (isect, '(i4.4)') sect
+    hname = stem // sp // l // p // isect
+
+    if (xdr) then   ! open XDR output file for Hamiltonian
+       ixdro = open_xdr (file=TRIM(hname), action='write')
+       call xdr_io (ixdro, nh)
+       call xdr_io (ixdro, nsect)
+    else
+       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
+            status='replace', form='unformatted', action='write', &
+            iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
+               TRIM(hname) // ' iostat = ', ios
+          stop
+       end if
+       write (ixdro) nh, nsect
+    end if
+    if (xdr) then
+       write (fo,'(/,a,i0,a,i0,a,i0,a)') ' Mesh ', imesh, ', Sector ', sect, ': Opening &
+            &XDR Representation Amplitude file: ' // TRIM(hname) 
+    else
+       write (fo,'(/,a,i0,a,i0,a,i0,a)') ' Mesh ', imesh, ', Sector ', sect, ': Opening &
+            &Native Representation Amplitude file: ' // TRIM(hname)
+    end if
+    call flush(fo)
+  end subroutine start_sa_file
+
+  subroutine start_sa_file_id (sqno, st)
+! open XDR file to hold sector reduced width amplitudes
+    use xdr_files, only: open_xdr, xdr_io
+    use energy_grid, only: fine_mesh
+    use slp, only: qno, val_qno
+    use rmx1_in, only: xdr, split_prop
+    !use mpi
+    include 'mpif.h'
+
+    type(qno), intent(in)   :: sqno        ! scattering q#s
+    integer, intent(in)     :: st         ! second target spin
+    character(len=13)        :: hname
+    character(len=2)        :: sp, l
+    character(len=1)        :: p
+    character(len=4)        :: stem, s_mpi
+    integer                 :: ios, qn(3), status
+    integer                 :: nspn ! scattering spin.  
+    integer ( kind = 4 ) ierr
+! File naming convention:
+!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
+!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
+!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
+!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
+! For no-exchange LS files target spin multiplicity is used
+! For J files L is replaced by 2*J and total spin multiplicity by 0 
+!                                  or target spin multiplicity by 2*K
+! For no-exchange files L is replaced by 2*K 
+!                   and total spin multiplicity by 0
+
+    ixdro = fa
+    stem = 'AMPC'
+    if (fine_mesh) stem = 'AMPF'
+    qn = val_qno(sqno)
+    if (split_prop) then   ! channel splitting in use
+       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
+    else
+       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
+       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
+    end if
+    write (l,'(i2.2)') abs(qn(1))
+    write (p, '(i1.1)') qn(3)
+    write (s_mpi, '(i4.4)') taskid
+    hname = stem // sp // l // p // s_mpi
+
+
+    if (xdr) then   ! open XDR output file for Hamiltonian
+       ixdro = open_xdr (file=TRIM(hname), action='write')
+       call xdr_io (ixdro, nh)
+       call xdr_io (ixdro, nsect)
+    else
+       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
+            status='replace', form='unformatted', action='write', &
+            iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
+               TRIM(hname) // ' iostat = ', ios
+          stop
+       end if
+       write (ixdro) nh, nsect
+    end if
+    if (xdr) then
+       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
+            &TRIM(hname)
+    else
+       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
+            &TRIM(hname)
+    end if
+  end subroutine start_sa_file_id
+
+
+  subroutine close_sa_file ()
+   use xdr_files, only: close_xdr
+   call flush(ixdro)
+   call close_xdr(ixdro)
+   
+  end subroutine close_sa_file
+
+
+  subroutine ampltd (ra1, ra2, nc, nl, nh)
+! ampltd : computes surface amplitudes and writes to file
+    use xdr_files, only: xdr_io
+    use sec_ham, only: vbl, vbr
+    use dist_mat, only: evals
+    use rmx1_in, only: bug8, xdr
+!    use blacs, only: p0, q0, ctxt, io_processor, p_error, dlen_, &
+!         ctxt_, mb_
+    use error_out, only: error_check
+    use dist_mat, only: z
+    integer, intent(in)      :: nc         ! # channels
+    integer, intent(in)      :: nl         ! # Legendre functions
+    integer, intent(in)      :: nh         ! sec Hamiltonian dimension
+    real(wp), intent(in)     :: ra1, ra2   ! end pt radii of subrange
+    real(wp), allocatable    :: ampal(:), ampar(:) ! local amps
+    real(wp), allocatable    :: evec(:)
+    integer                  :: k, j, ik, i, jj, status
+    integer                  :: mm, nn, lwork, mb, nb, rsrc, csrc
+    integer                  :: ldd, isize, istart, iend
+    integer                  :: jsize, jstart, jend
+    real(wp)                 :: alpha, beta, sml, smr
+!    integer                  :: descx(dlen_)
+
+!    if (io_processor) then
+       allocate (ampal(nc), ampar(nc), stat=status)
+       if (status /= 0) then
+          write (fo,'(a,i6)') 'ampltd: allocation error = ', status
+          STOP
+       end if
+    
+       if (bug8 > 0) then
+          write (fo,'(a,e20.12)') 'Sector eigenvalue sum = ', SUM(evals)
+       else if (bug8 > 1) then
+          write (fo,'(a)') 'R-Matrix Eigenvalues (Ryd)'
+          write (fo,'(5f12.5)') evals
+       end if
+
+! write eigenvalues to disk
+       if (xdr) then
+          call xdr_io (ixdro, evals, nh)  ! write R-matrix eigenvalues
+       else
+          write (ixdro) evals
+       end if
+!    end if
+
+    allocate (evec(nh), stat=status)
+    call error_check (status, 'ampltd: evec allocation')
+    lwork = 96
+!    lwork = descz(mb_)   ! eigenvectgor block size
+!    ctxt = descz(ctxt_)  ! context
+
+! form descriptor for evec
+    mm = MAX(1, MIN(nh, lwork))
+    nn = MAX(1, INT(lwork / mm))
+    mb = mm
+    nb = nn
+!    rsrc = p0
+!    csrc = q0
+!    ldd = MAX(1, mm)
+!    call descset (descx, mm, nn, mb, nb, rsrc, csrc, ctxt, ldd)
+
+    alpha = 1.0_wp
+    beta = 0.0_wp
+    sml = 0.0_wp
+    smr = 0.0_wp
+    evec = z
+    jst: do jstart = 1, nh, nn
+       jend = MIN(nh, jstart+nn-1)
+       jsize = jend - jstart + 1
+ !      row_blocks: do istart = 1, nh, mm
+ !         iend = MIN(nh, istart+mm-1)
+ !         isize = iend - istart + 1
+ !         call pdgeadd ('notrans', isize, jsize, alpha, z, istart, &
+ !              jstart, descz, beta, evec(istart:iend), 1, 1, descx)
+ !      end do row_blocks
+ !      call BLACS_BARRIER (ctxt, 'A')
+! now have the full eigenvector on the i/o processor:
+
+!       if (io_processor) then
+          ampal = 0.0_wp
+          ampar = 0.0_wp
+          ik = 0
+          do i = 1, nc
+             do j = 1, nl
+                ik = ik + 1
+                ampal(i) = ampal(i) + evec(ik) * vbl(j)
+                ampar(i) = ampar(i) + evec(ik) * vbr(j)
+             end do
+          end do
+          
+          if (bug8 > 0) then
+             sml = sml + SUM(ampal)
+             smr = smr + SUM(ampar)
+          else if (bug8 > 2) then
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra1,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampal
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra2,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampar
+          end if
+
+! write sector reduced width amplitudes to disk
+          if (xdr) then
+             call xdr_io (ixdro, ampal, nc) ! write left amplitudes to disk
+             call xdr_io (ixdro, ampar, nc) ! write right amplitudes to disk
+          else
+             write (ixdro) ampal(1:nc)
+             write (ixdro) ampar(1:nc)
+          end if
+ !      end if
+    end do jst
+!    if (io_processor) then
+       if (bug8 > 0) write (fo,'(a,2e20.12)') 'RWAmplitude sums = ', &
+            sml, smr
+       deallocate (ampal, ampar, evec, stat=status)
+!    else
+!       deallocate (evec, stat=status)
+!    end if
+    call error_check (status, 'ampltd: deallocation')
+  end subroutine ampltd
+end module rw_amplitudes

pfarm/cpu/src_mpi_omp_atom/blacs.f90

+module blacs
+! store blacs context and processor data
+! Time-stamp: "2003-09-23 13:41:59 cjn"
+  implicit none
+
+! define parameters used to descibe BLACS distributed arrays:
+  integer, parameter :: dlen_ = 9
+  integer, parameter :: block_cyclic_2d = 1
+  integer, parameter :: ctxt_ = 2
+  integer, parameter :: m_ = 3
+  integer, parameter :: n_ = 4
+  integer, parameter :: mb_ = 5
+  integer, parameter :: nb_ = 6
+  integer, parameter :: rsrc_ = 7
+  integer, parameter :: csrc_ = 8
+  integer, parameter :: lld_ = 9
+
+! initial linear context:
+  integer, save      :: ictxt  ! initial blacs context handle
+  integer, save      :: iam    ! processor id in initial context
+  integer, save      :: nprocs ! total # processors
+
+! context for main calculation: p * q <= nprocs
+  integer, save      :: p, q   ! # processors in grid
+  integer, save      :: ctxt   ! blacs context handle
+  integer, save      :: mycol, myrow ! processor coordinates
+  logical, save      :: io_processor   ! i/o processor flag
+  integer, save      :: p0, q0 ! i/o processor coordinates
+  integer, save      :: nblock ! blacs blocking factor
+  integer, save      :: mynumrows   ! Hamiltonian local # rows
+  integer, save      :: mynumcols   ! Hamiltonian local # cols
+
+  integer, save      :: cur_ctxt=-5 ! current contxt for p_error
+
+  private
+  public p_error, set_cur_ctxt
+  public dlen_, block_cyclic_2d, ctxt_, m_, n_, mb_, nb_
+  public rsrc_, csrc_, lld_
+  public ictxt, nprocs, iam
+  public ctxt, io_processor, p, q, mycol, myrow, p0, q0
+  public nblock, mynumrows, mynumcols
+
+contains
+
+  subroutine set_cur_ctxt (n)
+! Define the Blacs context handle to be used by p_error
+    integer, intent(in)    :: n
+    select case (n)
+    case (1)
+       cur_ctxt = ictxt
+    case (2)
+       cur_ctxt = ctxt
+    case default
+       cur_ctxt = ictxt
+    end select
+  end subroutine set_cur_ctxt
+
+  subroutine p_error (err, msg)
+    use io_units, only: fo
+    integer, intent(in)          :: err ! flag
+    character(len=*), intent(in) :: msg   ! error message
+    integer                      :: no
+
+    no = 1
+    if (cur_ctxt /= -5) then
+       call igsum2d (cur_ctxt, 'all', ' ', 1, 1, err, 1, -1, -1)
+       if (err /= 0) then
+          if (iam == 0) write (fo, fmt='(a)') msg
+          call BLACS_ABORT (cur_ctxt, no)
+       end if
+    else
+       if (err /= 0) then
+          write (fo, fmt='(a)') msg
+          call BLACS_GET (-1, 0, cur_ctxt) ! find default context
+          call BLACS_ABORT (cur_ctxt, no)
+       end if
+    end if
+  end subroutine p_error
+
+end module blacs

pfarm/cpu/src_mpi_omp_atom/def_sectors.f90

+module def_sectors
+! Time-stamp: "2003-10-14 12:51:24 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use error_out, only: error_check
+  use mpi_params, only: io_processor
+  implicit none
+
+  private
+  public get_nsect_nl, hack_get_nsect_nl, reset_sec, sectors
+  public nsect, nh, nl, nx, asect, w_sect
+
+  integer, save               :: nsect   ! # sectors
+  integer, save               :: nh      ! Hamiltonian dimension (all channels)
+  integer, save               :: nl      ! Legendre order
+  integer, save               :: nx      ! # Gauss points
+  real(wp), allocatable, save :: asect(:)
+  real(wp)                    :: w_sect
+
+contains
+
+  subroutine reset_sec
+! deallocate memory
+    integer            :: status
+    deallocate (asect, stat=status)
+    call error_check(status, " Error in reset_sec ")
+  end subroutine reset_sec
+
+  subroutine sectors (nc)
+ ! calculate sector end; rlast lies on the end of last sector
+    use energy_grid, only: fine_mesh, ebig
+    use rmx1_in, only: nhmax, ncmax, rafin
+    use hfile_data, only: rmatr
+    integer, intent(in) :: nc    ! # channels
+    real(wp)            :: rfirst, rlast ! propagation range
+    integer             :: ibuf(4), status, i
+    real(wp)            :: dum(1)
+
+
+       rfirst = rmatr
+       rlast = rafin
+       if (rafin < rmatr) then
+          write (fo,'(a)') 'sectors: propagation range error'
+          write (fo,'(a,2f12.6)') 'rmatr, rafin = ', rmatr, rafin
+          STOP
+       end if
+       call hack_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)
+       
+       w_sect = (rlast - rfirst) / REAL(nsect,wp)   ! actual sector width
+       allocate (asect(nsect+1), stat=status)
+       asect(1) = rfirst
+       do i = 1, nsect
+          asect(i+1) = asect(1) + i * w_sect
+       end do
+       nh = nl * nc
+
+! output sector radius
+       if(io_processor) then
+         write (fo,'(/a,i4)') 'Number of basis functions *** HACKED for Fine & Coarse Regions ***, nl = ', nl
+         write (fo,'(a,i4)') 'Number of abscissae, nx       = ', nx
+         write (fo,'(/a,i8/)') 'Hamiltonian dimension = ', nh
+         write (fo,'(a,f12.6)') 'Sector width   = ', w_sect
+         do i = 1, nsect+1
+           write (fo,'(i4,5x,a,f12.6)') i, 'asect = ', asect(i)
+         end do
+       end if
+   end subroutine sectors
+
+  subroutine get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+       anlmax = nhmax / ncmax
+       ansect = (rlast - rfirst) * kmax / anlmax
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+    else ! coarse mesh, therefore set nl to 10
+       nl = 10
+       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+    end if
+    nx = nl + 5
+  end subroutine get_nsect_nl
+
+  subroutine hack_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! **** AGS Hacked to provide larger Hamiltonian matrices ****
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+    CHARACTER(len=255) :: c_nl
+    integer :: i_nl, status_val
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+!       anlmax = nhmax / ncmax
+!       ansect = (rlast - rfirst) * kmax / anlmax
+!       nsect = MAX(INT(ansect+1.0_wp), 2)
+!       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+!       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 4)
+       nsect = 5
+       i_nl=-1
+       CALL get_environment_variable("DSYEVD_NL_1", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if(io_processor) WRITE (*,*) "Environmental variable DSYEVD_NL_1 : ",i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 12
+       end if
+    else ! coarse mesh, therefore set nl to 10
+       i_nl=-1
+       CALL get_environment_variable("DSYEVD_NL_2", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if(io_processor) WRITE (*,*) "Environmental variable DSYEVD_NL_2 : " ,i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 6
+       end if
+!       nl = 10
+!       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+!       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nsect = 20
+    end if
+    nx = nl + 5
+  end subroutine hack_get_nsect_nl
+
+end module def_sectors

pfarm/cpu/src_mpi_omp_atom/def_sectors_bench.f90

+module def_sectors
+! Time-stamp: "2003-10-14 12:51:24 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use error_out, only: error_check
+  use mpi_params, only: io_processor
+  implicit none
+
+  private
+  public reset_sec, sectors
+  public nsect, nh, nl, nx, asect, w_sect
+
+  integer, save               :: nsect   ! # sectors
+  integer, save               :: nh      ! Hamiltonian dimension (all channels)
+  integer, save               :: nl      ! Legendre order
+  integer, save               :: nx      ! # Gauss points
+  real(wp), allocatable, save :: asect(:)
+  real(wp)                    :: w_sect
+
+contains
+
+  subroutine reset_sec
+! deallocate memory
+    integer            :: status
+    deallocate (asect, stat=status)
+    call error_check(status, " Error in reset_sec ")
+  end subroutine reset_sec
+
+  subroutine sectors (nc)
+ ! calculate sector end; rlast lies on the end of last sector
+    use energy_grid, only: fine_mesh, ebig
+    use rmx1_in, only: nhmax, ncmax, rafin
+    use hfile_data, only: rmatr
+    integer, intent(in) :: nc    ! # channels
+    real(wp)            :: rfirst, rlast ! propagation range
+    integer             :: ibuf(4), status, i
+    real(wp)            :: dum(1)
+
+
+       rfirst = rmatr
+       rlast = rafin
+       if (rafin < rmatr) then
+          write (fo,'(a)') 'sectors: propagation range error'
+          write (fo,'(a,2f12.6)') 'rmatr, rafin = ', rmatr, rafin
+          STOP
+       end if
+
+       ! Benchmark version of code allowing easier problem size variation by user
+       call bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)
+
+       w_sect = (rlast - rfirst) / REAL(nsect,wp)   ! actual sector width
+       allocate (asect(nsect+1), stat=status)
+       asect(1) = rfirst
+       do i = 1, nsect
+          asect(i+1) = asect(1) + i * w_sect
+       end do
+       nh = nl * nc
+
+! output sector radius
+       if(io_processor) then
+         write (fo,'(/a,i4)') 'Number of basis functions (Benchmark Case) = ', nl
+         write (fo,'(a,i4)') 'Number of abscissae, nx = ', nx
+         write (fo,'(a,i8)') 'Hamiltonian dimension = ', nh
+         write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
+         write (fo,'(a,f12.6)') 'Sector width   = ', w_sect
+         do i = 1, nsect+1
+           write (fo,'(i4,5x,a,f12.6)') i, 'asect = ', asect(i)
+         end do
+         write (fo,'(/)')
+         call flush(fo)
+       end if
+   end subroutine sectors
+
+  subroutine get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+       anlmax = nhmax / ncmax
+       ansect = (rlast - rfirst) * kmax / anlmax
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+    else ! coarse mesh, therefore set nl to 10
+       nl = 10
+       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+    end if
+    nx = nl + 5
+  end subroutine get_nsect_nl
+
+  subroutine bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+    ! **** Benchmark version of routine to allow greater user control of problem size ****
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+    CHARACTER(len=255) :: c_nl, c_nsect
+    integer :: i_nl, i_nsect, status_val
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+
+       nsect = -1
+       i_nl = -1
+       CALL get_environment_variable("RMX_NL_FINE", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 12
+       end if
+
+       CALL get_environment_variable("RMX_NSECT_FINE", c_nsect, STATUS=status_val)
+       read( c_nsect, '(i10)' ) i_nsect
+       if (i_nsect > 0 ) then
+          nsect = i_nsect
+       else
+          nsect = 5
+       end if
+
+       if(io_processor) then
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Fine Region = ", nsect
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Fine Region = ", nl
+       end if
+
+    else ! coarse mesh
+       i_nsect=-1
+       i_nl=-1
+       CALL get_environment_variable("RMX_NL_COARSE", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 6
+       end if
+
+       CALL get_environment_variable("RMX_NSECT_COARSE", c_nsect, STATUS=status_val)
+       read( c_nsect, '(i10)' ) i_nsect
+       if (i_nsect > 0 ) then
+          nsect = i_nsect
+       else
+          nsect = 20
+       end if
+
+       if(io_processor) then
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Coarse Region = ", nsect
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Coarse Region = ", nl
+       end if
+
+    end if
+
+    ! number of Gauss points
+    nx = nl + 5
+
+  end subroutine bench_get_nsect_nl
+
+end module def_sectors

pfarm/cpu/src_mpi_omp_atom/dist_mat.f90

+module dist_mat
+! Form diestributed sector Hamiltonian matrix and diagonalize
+! Time-stamp: "2003-09-23 17:08:23 cjn"
+  use io_units, only: fo
+  use precisn, only: wp, ip_long
+!  use blacs, only: ctxt, myrow, mycol, p, q, nblock, mynumrows,&
+!       mynumcols, p_error, dlen_, nb_, csrc_
+  use error_out, only: error_check
+  implicit none
+
+  integer, save               :: ny           ! a dimension
+  real(wp), allocatable, save :: a(:)         ! distributed matrix
+  real(wp), allocatable, save :: z(:)         ! eigenvectors
+  real(wp), allocatable, save :: evals(:)     ! eigenvalues
+!  integer, save               :: desca(dlen_) ! Hamiltonian descriptor
+!  integer, save               :: descz(dlen_) ! Eigenvector descriptor
+
+  private
+  public A_fill, sh_diag, kill_az
+  public evals, z
+
+contains
+
+  subroutine A_fill (n)
+! fill lower triangle of symmetric matrix, A
+    use sec_ham, only: h_el
+    integer, intent(in)   :: n       ! matrix dimension
+    integer               :: j, jj, jloc, jglob, i, ii, iloc, iglob
+    integer               :: status, info, numroc
+    integer :: ind
+
+    ny = n
+
+! Compute local dimensions of matrix A
+!    mynumrows = NUMROC(n, nblock, myrow, 0, p)
+!    mynumcols = NUMROC(n, nblock, mycol, 0, q)
+
+! Allocate local part of A
+    allocate (a(ny*ny), z(ny*ny), stat=status)
+    call error_check (status, ' Allocation Error, A_fill')
+
+    a = 0.0_wp        ! initialize A
+
+! Form array descriptor od destributed matrix A
+!    call descinit (desca, n, n, nblock, nblock, 0, 0, ctxt, mynumrows,&
+!         info)
+!    descz = desca
+
+!    do j = 1, mynumcols, nblock                   ! local column blocks
+!       do jj = 1, MIN(nblock, mynumcols-j+1)      ! all columns
+!          jloc = j - 1 + jj                       ! local column index
+!          jglob = (j-1) * q + mycol * nblock + jj ! global column index
+
+!          do i = 1, mynumrows, nblock
+!             do ii = 1, MIN(nblock, mynumrows-i+1)
+!                iloc = i - 1 + ii                       ! local index
+!                iglob = (i-1) * p + myrow * nblock + ii ! global index
+!                if (iglob > jglob) cycle
+                
+!                a(iloc+(jloc-1)*mynumrows) = h_el(iglob, jglob)
+
+!             end do
+!          end do
+!       end do
+!    end do
+
+! sequential version - full matrix
+  ind = 0
+  do j = 1, ny
+    do i = 1, j
+      ind = ind + 1
+      a(ind) = h_el(i,j)
+    end do
+  end do
+
+  end subroutine A_fill
+
+  subroutine sh_diag(imesh, sect)
+! diagonalize the matrix a
+    integer, intent(in)          :: imesh, sect
+    real(wp), allocatable        :: work(:)
+    integer, allocatable         :: iwork(:)
+    integer                      :: lwork, liwork
+    real(wp)                     :: alpha, beta
+    character(len=1)             :: jobz, uplo
+    logical                      :: ismycol
+    real(wp)                     :: dnorm2
+    integer                      :: status, ii, nstrt, trilwmin
+    integer                      :: np, nq, lda, iproc, info, numroc
+    integer                      :: indxg2p, lwork1, ip, iq, imyrow, imycol
+    integer :: rows, cols
+    real(wp)                     :: t0, t1
+    integer(ip_long)            :: c0, c1, cr
+    real(wp), external :: dnrm2
+
+! initialize eigenvector distributed array, z.
+     alpha = 0.0_wp
+     beta = 1.0_wp
+!     call pdlaset ('all', ny, ny, alpha, beta, z, 1, 1, descz)
+!     np = mynumrows
+!     nq = mynumcols5
+
+! perform diagonalization:
+     jobz = 'V'
+     uplo = 'U'
+     liwork = 3 + 5*ny
+     lwork = 1 + 6*ny + 2*ny*ny
+
+     allocate (work(lwork), iwork(liwork), evals(ny), stat=status)
+     call error_check (status, 'sh_diag: allocation error')
+
+     write (fo,'(/,a,i0,a,i0,a)') ' Mesh ',imesh,', Sector ',sect,': Starting DSYEVD '
+     call flush(fo)
+     
+     call cpu_time (t0)
+     call system_clock (count=c0)
+
+     call dsyevd (jobz, uplo, ny, a, ny, evals, work, lwork, iwork, liwork, info)
+
+     write (fo,'(/,a,i0,a,i0,a)') ' Mesh ',imesh,', Sector ',sect,': Ended DSYEVD '
+     call cpu_time (t1)
+     call system_clock (count=c1, count_rate=cr)
+
+     write (fo,'(/,a,i0,a,i0,a,f14.3,a)') ' Mesh ',imesh,', Sector ',sect,&
+          ': DSYEVD CPU time = ', t1 - t0, ' secs'
+     write (fo,'(a,i0,a,i0,a,f14.3,a,/)') ' Mesh ',imesh,', Sector ',sect,&
+          ': DSYEVD Elapsed time = ', REAL(c1-c0,wp) / REAL(cr,wp), ' secs'
+
+     call error_check (info, 'sh_diag: pdsyevd error')
+
+     deallocate (work, iwork, stat=status)
+     call error_check (status, 'sh_diag: deallocation error')
+
+     write (fo,'(/,a,i0,a,i0,a,5f14.4)') ' Mesh ',imesh,', Sector ',sect,': first five eigenvalues = ', evals(1:5)
+     write (fo,'(a,i0,a,i0,a,5f14.4,/)') ' Mesh ',imesh,', Sector ',sect,': final five eigenvalues = ', evals(ny-4:ny)
+     call flush(fo)
+     
+     z = a
+     evecs: do ii = 1, ny ! normalize eigenvectors
+! indxg2p computes process coord which posseses entry of a
+! distributed matrix specified by a global index INDXGLOB.
+!        iproc = indxg2p(ii, descz(nb_), mycol, descz(csrc_), q)
+!        ismycol = (iproc == mycol)
+!        if (ismycol) then
+!           call pdnrm2 (ny, dnorm2, z, 1, ii, descz, 1)
+            
+           dnorm2 = dnrm2(ny,z,1)    
+           alpha = 1.0_wp / REAL(dnorm2,wp)
+           call dscal(ny, alpha, z, 1)
+!        end if
+     end do evecs
+   end subroutine sh_diag
+
+   subroutine kill_az
+     integer       :: status
+     deallocate (a, z, evals, stat=status)
+     call error_check (status, 'kill_az: deallocation error')
+   end subroutine kill_az
+
+end module dist_mat

pfarm/cpu/src_mpi_omp_atom/energy_grid.f90

+module energy_grid
+! Define maximum scattering energy for energy grid region
+! Time-stamp: "2003-09-19 13:19:02 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use mpi_params, only: io_processor
+  implicit none
+
+  logical, save            :: fine_mesh ! fine/coarse energy grid
+  real(wp), save           :: ebig     ! largest scattering energy
+
+  private
+  public energy_max
+  public fine_mesh, ebig
+
+contains
+
+  subroutine energy_max (imesh, nmesh)
+! Define ebig, the maximum scattering energy for case
+    use rmx1_in, only:  emax, targ_high_sce, uniform_e_grid, belthr
+    use hfile_data, only: ntarg, etarg, tchl, nc
+    integer, intent(in)   :: imesh     ! current mesh
+    integer, intent(in)   :: nmesh     ! number of energy regions
+    real(wp)              :: e1
+
+    if (imesh == 1 .and. nmesh > 1) then ! resonance region (fine mesh)
+       fine_mesh = .true.
+       e1 = MINVAL(etarg(1:ntarg)) ! ground state energy threshold
+
+       if (targ_high_sce > 0) then ! sequence of final target
+          ebig = etarg(targ_high_sce) - e1 - belthr
+
+       else if (uniform_e_grid) then ! end on highest threshold
+          ebig = MAXVAL(etarg(1:ntarg)) - e1 - belthr
+          
+       else ! end on highest theshold for this symmetry
+          ebig = MAXVAL(etarg(tchl(1:nc))) - e1 - belthr
+       end if
+
+    else ! coarse energy mesh
+       fine_mesh = .false.
+       ebig = emax
+    end if
+    if(io_processor) write (fo,'(/a,es12.4,a/)') 'Max. energy, ebig = ', ebig, &
+         ' Scaled Rydbergs'
+  end subroutine energy_max
+
+end module energy_grid

pfarm/cpu/src_mpi_omp_atom/error_out.f90

+module error_out
+  implicit none
+
+  private
+  public error_check
+
+contains
+
+  subroutine error_check (err, msg)
+     use io_units, only: fo
+     integer, intent(in)          :: err ! flag
+     character(len=*), intent(in) :: msg   ! error message
+     integer                      :: no
+
+     if (err /= 0) then
+       write (fo, fmt='(a)') msg
+     end if
+   end subroutine error_check
+end module error_out
+
+
+

pfarm/cpu/src_mpi_omp_atom/io_units.f90

+module io_units
+! Time-stamp: "2003-06-12 07:37:43 cjn"
+  implicit none
+  private
+  public fi, fo, dsk, dskp1, dskp2
+  public set_dsk, set_dskp1, set_dskp2
+! ang values:
+  public jbufiv, jbuffv, jbufie, jbuffe
+  public jbufiv1, jbuffv1, jbufiv2, jbuffv2
+! ham values:
+  public jbuff1, jbuff2, jbuff4, jbuff5, itape1, itape3, jbufd
+  public jbufim, jbuffm, jbufev
+
+  public jibb, jfbb, jibc, jfbc, fl, ihunit, itm
+  public jbufivd, jbuffvd, jbufiv2d, jbuffv2d
+
+  public nfnml, fa
+
+  integer, parameter      :: fi = 5    ! standard input
+  integer, parameter      :: fo = 6    ! standard output
+  integer, parameter      :: nfnml = 7 ! namelist file
+  integer, parameter      :: fa = 9    ! sector data file
+
+
+  integer, save           :: dsk = 0   ! scratch disk unit
+  integer, save           :: dskp1 = 0   ! scratch disk unit
+  integer, save           :: dskp2 = 0   ! scratch disk unit
+
+  integer, save           :: jbufiv = 20
+  integer, save           :: jbuffv = 21
+  integer, save           :: jbuff1 = 24
+  integer, save           :: jbuff2 = 25
+  integer, save           :: itape1 = 26
+  integer, save           :: jbuff4 = 27
+  integer, save           :: jbufiv1 = 24    ! -> jbuff1
+  integer, save           :: jbuffv1 = 25    ! -> jbuff2
+  integer, save           :: jbufiv2 = 26    ! -> itape1
+  integer, save           :: jbuffv2 = 27    ! -> jbuff4
+  integer, save           :: jbuff5 = 28
+  integer, save           :: jbufd = 29
+  integer, save           :: itape3 = 10
+
+  integer, save           :: jbufie = 22
+  integer, save           :: jbuffe = 23
+
+  integer, save           :: jbufim = 45
+  integer, save           :: jbuffm = 46
+
+  integer, save           :: jbufev = 47
+
+! dipole units used in stmmat:
+  integer, save           :: jbufivd = 41
+  integer, save           :: jbuffvd = 42
+  integer, save           :: jbufiv2d = 43
+  integer, save           :: jbuffv2d = 44
+
+  integer, save           :: jibb = 30    ! abbi
+  integer, save           :: jfbb = 31    ! abbr
+  integer, save           :: jibc = 32    ! abci
+  integer, save           :: jfbc = 33    ! abcr
+
+  integer, save           :: fl = 36      ! tgt dataset unit
+  integer, save           :: ihunit = 37  ! unit for Hamiltonian output
+  integer, save           :: itm = 38     ! unit for tarmom file
+  
+contains
+  subroutine set_dsk (dsk_unit)
+    integer, intent(in)   :: dsk_unit
+    dsk = dsk_unit
+    return
+  end subroutine set_dsk
+
+  subroutine set_dskp1 (dsk_unit)
+    integer, intent(in)   :: dsk_unit
+    dskp1 = dsk_unit
+    return
+  end subroutine set_dskp1
+
+  subroutine set_dskp2 (dsk_unit)
+    integer, intent(in)   :: dsk_unit
+    dskp2 = dsk_unit
+    return
+  end subroutine set_dskp2
+
+end module io_units

pfarm/cpu/src_mpi_omp_atom/mpi_params.f90

+module mpi_params
+ implicit none
+ private
+ public taskid, numtasks
+ public io_processor
+
+ integer, save :: taskid, numtasks
+ logical, save :: io_processor=.false.
+
+end module mpi_params
+  

pfarm/cpu/src_mpi_omp_atom/pdg_ctl.f90

+module pdg_ctl
+! control diagonalization of sector Hamiltonians and RWA calculation
+! Time-stamp: "2003-10-14 12:49:30 cjn"
+! modified for jpi_coupling april to dec 2006 vmb
+  use precisn, only: wp
+  use io_units, only: fo
+  implicit none
+
+  private
+  public diag
+
+contains
+
+  subroutine diag (nc, ncp, spins, sqno, imesh)
+   use hfile_data, only: lchl, kschl, tchl
+   use rw_amplitudes, only: start_sa_file, close_sa_file, ampltd
+   use potl_cofs, only: cfs, cfsj, del_cf
+   use rmx1_in, only: inc_lrp_prop
+   use slp, only: qno
+   use def_sectors, only: nsect, nl, nx, w_sect, asect, reset_sec, &
+        sectors
+   use sector_potl, only: potl, potl_pams, reset_potl
+   use dist_mat, only: a_fill, sh_diag, kill_az
+!   use blacs, only: ctxt, p, q, myrow, mycol, io_processor, p0, &
+!        q0, p_error
+   use sec_ham, only: h_reset, gauleg, legndr, def_rb_vals, xi, def_ncq
+   use rmx1_in, only: bug9, jpi_coupling
+   use error_out, only: error_check
+   use mpi_params, only: io_processor, taskid, numtasks
+!   use magma
+!   use mpi
+   include 'mpif.h'
+
+   integer, intent(in)      :: nc       ! # channels
+   integer, intent(in)      :: ncp
+   type(qno), intent(in)    :: sqno     ! scattering q#s
+   integer, intent(in)      :: spins(2) ! target spins
+   integer                  :: imesh    ! COARSE,FINE mesh - imesh=1,2
+   real(wp), allocatable    :: r(:)
+   real(wp)                 :: ral, rar
+   integer                  :: sect, st, status, nh, pp, qq
+   integer, save            :: nsect_fine, nsect_coarse, mesh_taskid, coarse_mesh_first_taskid
+   integer ( kind = 4 ) ierr
+
+   call sectors (nc) ! define sector boundaries
+   call gauleg  ! quadrature formula
+   call legndr  ! Legendre functions
+   call potl_pams (nc, nx, ncp)
+   call def_ncq (ncp) ! pass ethr offset in h_el
+
+! generate asymptotic potential coefficients
+   if (inc_lrp_prop) then
+   if (.not. jpi_coupling) then
+      call cfs (nc, tchl(ncp:), lchl(ncp:), sqno, spins)
+   else
+      call cfsj (nc, tchl(ncp:), lchl(ncp:), kschl(ncp:), sqno, spins)
+   end if
+   end if
+   allocate (r(nx), stat=status)
+   call error_check (status, 'pdiag: allocation r')
+   
+   nh = nc * nl
+   if ((bug9 > 0).and.io_processor) then
+      write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
+      write (fo,'(a,i4)') 'Partition Hamiltonian dimension = ', nh
+   end if
+   call flush(fo)
+
+    ! Allow tasks to cycle through both meshes ensuring best load-balancing
+    if(imesh==1) then ! fine mesh
+       mesh_taskid = taskid
+       nsect_fine = nsect
+    else ! coarse mesh
+
+      nsect_coarse = nsect
+      mesh_taskid = taskid - nsect_fine
+      coarse_mesh_first_taskid = MOD(nsect_fine,numtasks)
+      mesh_taskid = MOD((taskid - coarse_mesh_first_taskid + numtasks),numtasks)
+    end if
+
+!      if(numtasks<=nsect_fine) then ! numtasks < nsect_fine
+!        if(MOD(taskid,numtasks).ge.(nsect_fine-numtasks)) then
+!           mesh_taskid = mesh_taskid + numtasks
+!        else
+!           mesh_taskid = mesh_taskid + (2*numtasks)
+!        end if
+!      else ! numtasks > nsect_fine
+!        if(MOD(taskid,numtasks).lt.nsect_fine) mesh_taskid = mesh_taskid + numtasks 
+!      end if
+!    end if
+    
+!   do sect = startsect, endsect
+    do sect = 1, nsect
+      if(numtasks.gt.1) then
+        if (MOD(sect,numtasks)==MOD(mesh_taskid+1,numtasks)) then
+          if(imesh==1) then
+            write(fo,'(a,i0,a,i0,a,i0)')' FINE Region MPI Task ', taskid, ', Mesh Task ', mesh_taskid, ' calculating sector ', sect
+          else
+            write(fo,'(a,i0,a,i0,a,i0)')' COARSE Region MPI Task ', taskid, ', Mesh Task ', mesh_taskid, ' calculating sector ', sect
+          end if 
+          call flush(fo)
+        else
+          cycle
+        end if
+      end if
+      call start_sa_file (sqno, spins(1), sect, imesh)
+      ral = asect(sect)
+      rar = asect(sect+1)
+      r = 0.5_wp * ((rar - ral) * xi + rar + ral)
+      call def_rb_vals (ral, rar)
+      call potl (nc, nx, r)
+      call A_fill (nh) ! form sector H
+      call sh_diag (imesh, sect)      ! diagonalize sector H
+      call ampltd (ral, rar, nc, nl, nh)
+      call kill_az
+      call close_sa_file()
+   end do
+
+   call h_reset
+   call reset_sec
+   call reset_potl
+   if (inc_lrp_prop) call del_cf
+   deallocate (r, stat=status)
+   call error_check (status, 'pdiag: deallocation')
+   
+ end subroutine diag
+
+end module pdg_ctl

pfarm/cpu/src_mpi_omp_atom/potl.f90

+module sector_potl
+! generate sector potential
+! Time-stamp: "2003-09-19 11:52:40 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use potl_cofs, only: cf, lamd ! asymptotic potential coefficients
+  use rmx1_in, only: inc_lrp_prop, n_lambda
+  implicit none
+
+  integer, save                 :: ion    ! ionicity
+  integer, save                 :: lmx    ! # multipoles
+  integer, allocatable, save    :: lc(:)  ! channel orbital a.m.
+  integer, save                 :: nc
+  real(wp), allocatable, save   :: v(:,:,:)
+
+  private
+  public potl, potl_pams, reset_potl
+  public v, nc
+
+contains
+
+  subroutine potl_pams (nct, nx, ncp)
+! set potential parameters
+    use hfile_data, only: lchl
+!    use blacs, only: p_error
+    use scaling, only: get_ionicity
+    use error_out, only: error_check
+    integer, intent(in)   :: nct  ! # channels in partition
+    integer, intent(in)   :: nx   ! # radial points
+    integer, intent(in)   :: ncp  ! channel origin
+    integer               :: status
+
+    nc = nct
+    allocate (v(nx,nc,nc), lc(nc), stat=status)
+    call error_check (status, 'potl_pams: allocation')
+    lc(1:nc) = lchl(ncp:ncp+nc-1)
+    ion = get_ionicity()
+    lmx = n_lambda
+  end subroutine potl_pams
+
+  subroutine reset_potl
+    integer      :: status
+    deallocate (v, lc, stat=status)
+    if (status /= 0) then
+       write (fo,'(a)') 'reset_potl: deallocation error'
+       STOP
+    end if
+  end subroutine reset_potl
+
+  subroutine potl (nc, nx, r)
+! potl: potentials expanded in inverse powers of radial distance;
+    use rmx1_in, only: bug7, n_lambda
+    integer, intent(in)      :: nc              ! no. channels
+    integer, intent(in)      :: nx              ! no. basis fns
+    real(wp), intent(in)     :: r(nx)            ! radial points
+    real(wp)                 :: dg(nx)
+    real(wp)                 :: vt
+    integer                  :: i, j, k, locc, lbb, par,  lm
+
+    if (inc_lrp_prop) then  ! asymptotic potential included
+       lm = MIN(lmx, n_lambda)
+       v = 0.0_wp
+       do j = 1, nc
+          dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
+          locc = j*(j-1)/2
+          do i = 1, j
+             if (i == j) v(:,j,i) = dg
+             do k = 1, lm, 2
+                lbb = (k + 1) / 2
+                v(:,i,j) = v(:,i,j) + cf(locc+i,lbb) &
+                     / r**(lamd(locc+i,lbb)+1)
+             end do
+          end do
+       end do
+    else
+       v = 0.0_wp
+       do j = 1, nc
+          dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
+          v(:,j,j) = dg
+       end do
+    end if
+    if (bug7 > 0) then
+!       do i = 1, nc
+       do i = 1, 3
+!          do j = 1, nc
+          do j = 1, 3
+             write (fo,'(a,i4,i4)') ' Potential (v) matrix:, col, &
+                  &row ',i,j
+             write (fo,'(6f14.7)') v(:,i,j)
+          end do
+       end do
+    end if
+  end subroutine potl
+end module sector_potl

pfarm/cpu/src_mpi_omp_atom/potl_coefs.f90

+module potl_cofs
+!  Generate asymptotic coefficients from target multipoles
+! Time-stamp: "2003-09-23 09:42:38 cjn"
+! modified for jpi_coupling april to dec 2006 vmb
+  use precisn, only: wp
+  use io_units, only: fo
+  implicit none
+
+  private
+  public cfs, cfsj, cf, lamd, del_cf, rdtmom, rdjtmom
+
+  logical, save      :: gams = .false.
+  integer, parameter :: ngam = 500             ! n entries in gam table
+  real(wp), save     :: gam(ngam)              ! gam table
+
+  type tmblock    ! transition moment type definition
+     sequence
+     integer           :: spn         ! spin value
+     integer           :: ntm         ! # transition moments
+     integer, pointer  :: lmds(:)     ! multipole index
+     integer, pointer  :: lbls(:)     ! locator array
+     real(wp), pointer :: tms(:)     ! transition moments
+  end type tmblock
+
+  integer, save                    :: ntmb    ! # target spins
+  type(tmblock), allocatable, save :: tmb(:)   ! target transition moments
+
+  real(wp), allocatable, save  :: cf(:,:)
+  integer, allocatable, save   :: lamd(:,:)
+
+  integer, save                :: nt    ! # targets
+  integer, allocatable, save   :: tj(:), tl(:), tspn(:), tpty(:)
+  real(wp), allocatable, save  :: tm(:,:) 
+
+contains
+
+  subroutine cfs (nc, tchl, lchl, sqno, spins)
+! generate scaled asymptotic potential coefficients
+    use slp, only: qno, val_qno
+    use hfile_data, only: lamax
+!    use blacs, only: ctxt, p_error
+    use scaling, only: scale_cf
+    use error_out, only: error_check
+
+    integer, intent(in)    :: nc  ! no channels
+    integer, intent(in)    :: tchl(nc) ! target sequence of channel
+    integer, intent(in)    :: lchl(nc) ! orb. a.m. of channel
+    type(qno), intent(in)  :: sqno     ! scattering a. #s
+    integer, intent(in)    :: spins(2) ! target spins required
+    real(wp)               :: phz, z, drac
+    integer                :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
+    integer                :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
+    integer                :: nspnt,  qn(3), lrgl, nspn, npty, isp
+    integer                :: status, smx
+
+    smx = lamax
+    allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
+         stat=status)
+    call error_check (status, 'cfs: allocation error')
+    cf = 0.0_wp
+    lamd = 0
+
+    target_spins: do isp = 1, 2 ! 1 spin value or 2
+       nspnt = spins(isp)
+       if (nspnt == -999) cycle ! only one spin value
+       qn = val_qno(sqno)
+       lrgl = qn(1)   ! scattering orb am
+       nspn = qn(2)   ! scattering spin multiplicity
+       npty = qn(3)   ! scattering parity
+
+       call get_tm (nspnt, smx, nt)
+
+       phz = 1
+       if (MOD(lrgl+npty,2) == 1) phz = -1
+       do j = 1, nc
+          tarj = tchl(j)
+          if (tspn(tarj) /= nspnt) cycle
+          lj = lchl(j)
+          tjl = tl(tarj)
+          tjs = tspn(tarj)
+          tjp = tpty(tarj)
+          locc = j*(j-1)/2
+          do i = 1, j
+             tari = tchl(i)
+             if (tspn(tari) /= nspnt) cycle
+             li= lchl(i)
+             til = tl(tari)
+             tis = tspn(tari)
+             tip = tpty(tari)
+             if (tari > tarj) then ! ensure input target sequencing
+                tt = tarj
+                tarj = tari
+                tari = tt
+             end if
+             loct = tarj * (tarj - 1) / 2 + tari
+             lbx = MIN(li+lj,til+tjl,smx)    ! min lambda
+             lbn = MAX(ABS(li-lj),ABS(til-tjl),1) ! max lambda
+             do lb = lbn, lbx, 1
+                if (MOD(tjp+tip+lb,2) == 1) cycle
+                lbb = (lb + 1) / 2
+                if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
+                drac = dracah(2*til, 2*tjl, 2*li, 2*lj, 2*lb, 2*lrgl)
+                z = phz * drac * rme(li, lj, lb)
+                cf(locc+i,lbb) = tm(loct,lbb) * z
+                lamd(locc+i,lbb) = lb
+             end do
+          end do
+       end do
+       deallocate (tm, stat=status)
+       if (status /= 0) STOP
+    end do target_spins
+    call scale_cf (cf, lamd)
+  end subroutine cfs
+
+  subroutine cfsj (nc, tchl, lchl, kschl, sqno, spins)
+! generate scaled asymptotic potential coefficients
+    use slp, only: qno, val_qno
+    use hfile_data, only: lamax
+!    use blacs, only: ctxt, p_error
+    use scaling, only: scale_cf
+    use error_out, only: error_check
+    integer, intent(in)    :: nc  ! no channels
+    integer, intent(in)    :: tchl(nc) ! target sequence of channel
+    integer, intent(in)    :: lchl(nc) ! orb. a.m. of channel
+    integer, intent(in)    :: kschl(nc) ! k q. no
+    type(qno), intent(in)  :: sqno     ! scattering q. #s
+    integer, intent(in)    :: spins(2) ! k splits required
+    real(wp)               :: phz, z, drac
+    integer                :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
+    integer                :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
+    integer                :: ki, kj, kbar, tij, tjj
+    integer                :: kspnt,  qn(3), lrgl, nspn, npty, isp
+    integer                :: status, smx
+
+    smx = lamax
+    allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
+         stat=status)
+    call error_check (status, 'cfs: allocation error')
+    cf = 0.0_wp
+    lamd = 0
+
+    channel_k_splits: do isp = 1, 2 ! 1 k value or 2
+       kspnt = spins(isp) ! value of ksplit = 2k
+       if (kspnt == -999) cycle ! only one k value
+       qn = val_qno(sqno)
+       lrgl = qn(1)   ! scattering orb am = 2J
+       nspn = qn(2)   ! scattering spin multiplicity = 0
+       npty = qn(3)   ! scattering parity
+
+       call get_tm (0, smx, nt) ! target moments are not split
+       phz = 1
+       if (MOD(kspnt, 2) == 0) kbar=kspnt/2
+       if (MOD(kspnt, 2) == 1) kbar=(kspnt+1)/2! kbar remains real
+       if (MOD(npty+kbar, 2).eq.1) phz = -1
+       do j = 1, nc
+          kj = kschl(j)
+          if (kj /= kspnt) cycle
+          tarj = tchl(j)
+          lj = lchl(j)
+          tjj = tj(tarj)
+          tjl = tl(tarj)
+          tjs = tspn(tarj)
+          tjp = tpty(tarj)
+          locc = j*(j-1)/2
+          do i = 1, j
+             tari = tchl(i)
+             li = lchl(i)
+             ki = kschl(i)
+             if (ki /= kj) cycle
+             tij = tj(tari)
+             til = tl(tari)
+             tis = tspn(tari)
+             tip = tpty(tari)
+             if (tari > tarj) then ! ensure input target sequencing
+                tt = tarj
+                tarj = tari
+                tari = tt
+             end if
+             loct = tarj * (tarj - 1) / 2 + tari
+             lbx = MIN(li+lj,smx)    ! max lambda
+             lbn = MAX(ABS(li-lj),1) ! min lambda
+             do lb = lbn, lbx, 1
+                if (MOD(tjp+tip+lb,2) == 1) cycle
+                lbb = (lb + 1) / 2
+                if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
+                drac = dracah(tij, tjj, 2*li, 2*lj, 2*lb, kspnt)
+                z = phz * drac * rme(li, lj, lb)
+                cf(locc+i,lbb) = tm(loct,lbb) * z
+                lamd(locc+i,lbb) = lb
+             end do
+          end do
+       end do
+       deallocate (tm, stat=status)
+       if (status /= 0) STOP
+    end do channel_k_splits
+    call scale_cf (cf, lamd)
+  end subroutine cfsj
+
+  subroutine del_cf
+!    use blacs, only: p_error
+    use error_out, only: error_check
+    integer        :: status
+    deallocate (cf, lamd, stat=status)
+    call error_check (status, 'del_cf: deallocation error')
+  end subroutine del_cf
+
+  function dracah (i, j, k, l, m, n)
+! racah coefficients
+    real(wp)           :: dracah
+    integer, intent(in) :: i, j, k, l, m, n ! input momenta * 2
+    integer :: j1, j2, j3, j4, j5, j6, j7, numin, numax, icount, kk, ki
+
+    if (.not.gams) then
+       call gam_tab
+       gams = .true.
+    end if
+    j1 = i + j + m
+    j2=  k + l + m
+    j3 = i + k + n
+    j4 = j + l + n
+    if (MOD(j1,2) /= 0 .or. MOD(j2,2) /= 0 .or. MOD(j3,2) /= 0 .or. &
+         MOD(j4,2) /= 0) then
+       dracah = 0.0_wp
+       return
+    end if
+    if (2*MAX(i,j,m) > j1 .or. 2*MAX(k,l,m) > j2 .or. 2*MAX(i,k,n) > j3 &
+         .or. 2*MAX(j,l,n) > j4) then
+       dracah = 0.0_wp
+       return
+    end if
+    j1 = j1 / 2
+    j2 = j2 / 2
+    j3 = j3 / 2
+    j4 = j4 / 2
+    j5 = (i + j + k + l) / 2
+    j6 = (i + l + m + n) / 2
+    j7 = (j + k + m + n) / 2
+    numin = MAX(j1, j2, j3, j4) + 1
+    numax = MIN(j5, j6, j7) + 1
+    dracah = 1.0_wp
+    icount = 0
+    if (numin < numax) then
+       do kk = numin+1, numax
+          ki = numax - icount
+          dracah = 1.0_wp - (dracah * real((ki * (j5-ki+2) * (j6-ki+2) * &
+               (j7-ki+2)),wp) / real((ki-1-j1) * (ki-1-j2) * (ki-1-j3) * &
+               (ki-1-j4),wp))
+          icount = icount + 1
+       end do
+    end if
+    dracah = dracah * (-1.0_wp)**(j5+numin+1) *                     &
+         EXP((gam(numin+1) - gam(numin-j1) - gam(numin-j2) -         &
+         gam(numin-j3) - gam(numin-j4) - gam(j5+2-numin) -           &
+         gam(j6+2-numin) - gam(j7+2-numin)) + ((gam(j1+1-i) +        &
+         gam(j1+1-j) + gam(j1+1-m) - gam(j1+2) + gam(j2+1-k) +       &
+         gam(j2+1-l) + gam(j2+1-m) - gam(j2+2) + gam(j3+1-i) +       &
+         gam(j3+1-k) + gam(j3+1-n) - gam(j3+2) + gam(j4+1-j) +       &
+         gam(j4+1-l) + gam(j4+1-n) - gam(j4+2)) / 2.0_wp))
+  end function dracah
+
+  function rme (l, lp, k)
+! reduced matrix element (l//c(k)//lp)
+! Ref: Fano and Racah, Irreducible Tensorial Sets, Ch 14, P. 81
+    real(wp)           :: rme
+    integer, intent(in) :: l, lp           ! initial, final orbital a.m.
+    integer, intent(in) :: k               ! order of spherical tensor
+    integer             :: i2g, ig, i1, i2, i3, imax, ier
+    real(wp)           :: qusqrt
+
+    if (.not.gams) then
+       call gam_tab
+       gams = .true.
+    end if
+    if (k < ABS(l-lp) .or. k > l+lp) then
+       write (fo,'(a,/,a,2i3)') 'rme: a.m. &
+            & coupling error', 'rme zeroed for l, lp, k = ', l, lp, k
+       rme = 0.0_wp
+       return
+    end if
+    i2g = l + lp + k
+    ig = i2g / 2
+    if (MOD(i2g,2) /= 0) then
+       rme = 0.0_wp
+       return
+    end if
+    if (ig < 0) then
+       write (fo,'(a,/,a,2i3)') 'rme: &
+            & angle does not match', &
+            'rme zeroed for l, lp, k = ', l, lp, k
+       rme = 0.0_wp
+       return
+    elseif (ig == 0) then
+       rme = 1.0_wp
+       return
+    end if
+    i1 = ig - l
+    i2 = ig - lp
+    i3 = ig - k
+    imax = MAX(ig+1, i2g+2, 2*i3+1, 2*i2+1, 2*i1+1)
+    if (imax > ngam) then
+       write (fo,'(a,2i6)') 'rme: table overflow, imax, ngam = ', &
+            imax, ngam
+       STOP
+    end if
+    qusqrt = log(real(2*l+1,wp)) + log(real(2*lp+1,wp)) + gam(2*i1+1) +&
+         gam(2*i2+1) + gam(2*i3+1) - gam(i2g+2)
+    rme = EXP(0.5_wp * qusqrt + gam(ig+1) - gam(i1+1) - gam(i2+1) -     &
+         gam(i3+1))
+  end function rme
+
+  subroutine gam_tab
+! logs of factorials table
+    real(wp) :: x
+    integer   :: i
+
+    gam(1) = 1.0_wp
+    gam(2) = 1.0_wp
+    x = 2.0_wp
+    do i = 3, 25
+       gam(i) = gam(i-1) * x
+       x = x + 1.0_wp
+    end do
+    gam(1:25) = LOG(gam(1:25))
+    x = 25.0_wp
+    do i = 26, ngam
+       gam(i) = gam(i-1) + LOG(x)
+       x = x + 1.0_wp
+    end do
+  end subroutine gam_tab
+
+  subroutine rdtmom
+    ! read file of target moments
+    !    use blacs, only: ictxt, p_error, io_processor
+    use io_units, only: fo, itm
+    use xdr_files, only: xdr_io, open_xdr, close_xdr
+    use rmx1_in, only: bug4, xdr, fltm
+    use error_out, only: error_check
+    real(wp), allocatable           :: tmom(:,:)
+    integer, allocatable            :: ntmom(:)
+    integer                         :: ibuf(2)
+    integer                         :: lamind, isp, ijl, ist, jst
+    integer                         :: lamlo, lamup, lam, lamn, nijl
+    integer                         :: ltmom, status, ios, s, i, j
+    integer                         :: nsv, spn, min_lam, max_lam, ntm
+    integer                         :: min_spn, max_spn, itmom, lamax, ip
+
+    if (xdr) then    ! output file used XDR format
+       itmom = open_xdr (file=TRIM(fltm), action='read')
+       call xdr_io (itmom, ibuf, 2)
+       nt = ibuf(1)
+       lamax = ibuf(2)
+       allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
+       if (status /= 0) then
+          write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
+          STOP
+       end if
+       call xdr_io (itmom, tl, nt)
+       call xdr_io (itmom, tspn, nt)
+       call xdr_io (itmom, tpty, nt)
+       call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
+       min_spn = ibuf(1)
+       max_spn = ibuf(2)
+    else             ! fortran unformatted output file
+       itmom = itm
+       open (unit=itmom, file=TRIM(fltm), form='unformatted', &
+            status='old', action='read', iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'rdtmom: error opening TARMOM, &
+               &iostat = ', ios
+          stop
+       end if
+       read (itmom) nt, lamax
+       allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
+       call error_check (status, 'rdtmom: allocation error(1)')
+       read (itmom) tl(1:nt)
+       read (itmom) tspn(1:nt)
+       read (itmom) tpty(1:nt)
+       read (itmom) min_spn, max_spn
+    end if
+
+    lamind = (lamax + 1) / 2
+    ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
+    nsv = (max_spn - min_spn) / 2 + 1
+    allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
+    if (status /= 0) then
+       write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
+       STOP
+    end if
+    if (bug4 > 0) then ! echo data read from TARMOM
+       write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
+            ' max lambda =', lamax
+       write (fo,'(a)') 'Target L:'
+       write (fo,'(25i3)') tl(1:nt)
+       write (fo,'(a)') 'Target S:'
+       write (fo,'(25i3)') tspn(1:nt)
+       write (fo,'(a)') 'Target P:'
+       write (fo,'(25i3)') tpty(1:nt)
+       write (fo,'(a)') 'Multiplicity limits:'
+       write (fo,'(25i3)') min_spn, max_spn
+    end if
+
+    s = 0
+    spins: do spn = min_spn, max_spn, 2
+       s = s + 1
+       if (xdr) then
+          call xdr_io (itmom, ibuf, 2)
+          isp = ibuf(1)
+          ntm = ibuf(2)
+          ntmom(s) = ntm
+          call xdr_io (itmom, tmom(1:ntm,s), ntm)
+       else
+          read (itmom) isp, ntm
+          ntmom(s) = ntm
+          read (itmom) tmom(1:ntm,s)
+       end if
+       if (bug4 > 0) then
+          write (fo,'(a,i4)') 'Target spin = ', isp
+          write (fo,'(a,i6)') 'Number of moments = ', ntm
+          write (fo,'(5f14.6)') tmom(1:ntm,s)
+       end if
+    end do spins
+
+    if (xdr) then
+       call close_xdr (itmom)
+    else
+       close (unit=itmom, status='keep')
+    end if
+
+    ! unpack tmom:
+    ntmb = nsv
+    allocate (tmb(nsv), stat=status)
+    call error_check (status, 'rdtmom: tmb allocation error')
+    s = 0
+    do spn = min_spn, max_spn, 2
+       s = s + 1
+       ntm = ntmom(s)
+       allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
+            stat=status)
+       call error_check (status, 'rdtmom: allocation error(5)')
+       tmb(s)%spn = spn
+       tmb(s)%ntm = ntm
+       ijl = 0
+       do i = 1, nt
+          if (tspn(i) /= spn) cycle
+          ip = i * (i - 1) / 2
+          do j = 1, i
+             if (tspn(j) /= spn) cycle
+             min_lam = ABS(tl(i) - tl(j))
+             if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
+                  min_lam + 1
+             if (min_lam == 0) min_lam = 2
+             max_lam = MIN(tl(i)+tl(j), lamax)
+             if ((min_lam+1)/2 <= (max_lam+1)/2) then
+                do lam = min_lam, max_lam, 2
+                   ijl = ijl + 1
+                   tmb(s)%lbls(ijl) = ip + j
+                   tmb(s)%lmds(ijl) = lam
+                   tmb(s)%tms(ijl) = tmom(ijl,s)
+                end do
+             end if
+          end do
+       end do
+    end do
+    deallocate (tmom, ntmom, stat=status)
+    call error_check (status, 'rdtmom: deallocation error')
+  end subroutine rdtmom
+
+  subroutine rdjtmom
+    ! read file of target moments when jpi_coupling is true
+    !    use blacs, only: ictxt, p_error, io_processor
+    use io_units, only: fo, itm
+    use xdr_files, only: xdr_io, open_xdr, close_xdr
+    use rmx1_in, only: bug4, xdr, fltm
+    use error_out, only: error_check
+
+    real(wp), allocatable           :: tmom(:,:)
+    integer, allocatable            :: ntmom(:)
+    integer                         :: ibuf(2)
+    integer                         :: lamind, isp, ijl, ist, jst
+    integer                         :: lamlo, lamup, lam, lamn, nijl
+    integer                         :: ltmom, status, ios, s, i, j
+    integer                         :: nsv, spn, min_lam, max_lam, ntm
+    integer                         :: min_spn, max_spn, itmom, lamax, ip
+
+    if (xdr) then    ! output file used XDR format
+       itmom = open_xdr (file=TRIM(fltm), action='read')
+       call xdr_io (itmom, ibuf, 2)
+       nt = ibuf(1)
+       lamax = ibuf(2)
+       allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
+       if (status /= 0) then
+          write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
+          STOP
+       end if
+       call xdr_io (itmom, tj, nt)
+       call xdr_io (itmom, tl, nt)
+       call xdr_io (itmom, tspn, nt)
+       call xdr_io (itmom, tpty, nt)
+       call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
+       min_spn = ibuf(1)
+       max_spn = ibuf(2)
+    else             ! fortran unformatted output file
+       itmom = itm
+       open (unit=itmom, file=TRIM(fltm), form='unformatted', &
+            status='old', action='read', iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'rdjtmom: error opening TARMOM, &
+               &iostat = ', ios
+          stop
+       end if
+       read (itmom) nt, lamax
+       allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
+       call error_check (status, 'rdjtmom: allocation error(1)')
+       read (itmom) tj(1:nt)
+       read (itmom) tl(1:nt)
+       read (itmom) tspn(1:nt)
+       read (itmom) tpty(1:nt)
+       read (itmom) min_spn, max_spn !set zero in fine.f90
+    end if
+
+    lamind = (lamax + 1) / 2
+    ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
+    nsv = (max_spn - min_spn) / 2 + 1 !nsv=1 here
+    allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
+    if (status /= 0) then
+       write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
+       STOP
+    end if
+    if (bug4 > 0) then ! echo data read from TARMOM
+       write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
+            ' max lambda =', lamax
+       write (fo,'(a)') 'Target 2J:'
+       write (fo,'(25i3)') tj(1:nt)
+       write (fo,'(a)') 'Target L:'
+       write (fo,'(25i3)') tl(1:nt)
+       write (fo,'(a)') 'Target S:'
+       write (fo,'(25i3)') tspn(1:nt)
+       write (fo,'(a)') 'Target P:'
+       write (fo,'(25i3)') tpty(1:nt)
+       write (fo,'(a)') 'Multiplicity limits:'
+       write (fo,'(25i3)') min_spn, max_spn
+    end if
+
+    s = 0
+    spins: do spn = min_spn, max_spn, 2! once thro loop
+       s = s + 1
+       if (xdr) then
+          call xdr_io (itmom, ibuf, 2)
+          isp = ibuf(1)
+          ntm = ibuf(2)
+          ntmom(s) = ntm
+          call xdr_io (itmom, tmom(1:ntm,s), ntm)
+       else
+          read (itmom) isp, ntm
+          ntmom(s) = ntm
+          read (itmom) tmom(1:ntm,s)
+       end if
+       if (bug4 > 0) then
+          write (fo,'(a,i4)') 'Notional Target spin = ', isp
+          write (fo,'(a,i6)') 'Number of moments = ', ntm
+          write (fo,'(5f14.6)') tmom(1:ntm,s)
+       end if
+    end do spins
+
+    if (xdr) then
+       call close_xdr (itmom)
+    else
+       close (unit=itmom, status='keep')
+    end if
+
+    ! unpack tmom:
+    ntmb = nsv
+    allocate (tmb(nsv), stat=status)
+    call error_check (status, 'rdjtmom: tmb allocation error')
+    s = 0
+    do spn = min_spn, max_spn, 2
+       s = s + 1
+       ntm = ntmom(s)
+       allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
+            stat=status)
+       call error_check (status, 'rdtmom: allocation error(5)')
+       tmb(s)%spn = spn
+       tmb(s)%ntm = ntm
+       ijl = 0
+       do i = 1, nt
+          ip = i * (i - 1) / 2
+          do j = 1, i
+             !             min_lam = ABS(tl(i) - tl(j))
+             !             if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
+             !                  min_lam + 1
+             min_lam = MOD(tpty(i)+tpty(j),2)
+             if (min_lam == 0) min_lam = 2
+             !             max_lam = MIN(tl(i)+tl(j), lamax)
+             max_lam = lamax
+             if ((min_lam+1)/2 <= (max_lam+1)/2) then
+                do lam = min_lam, max_lam, 2
+                   ijl = ijl + 1
+                   tmb(s)%lbls(ijl) = ip + j
+                   tmb(s)%lmds(ijl) = lam
+                   tmb(s)%tms(ijl) = tmom(ijl,s)
+                end do
+             end if
+          end do
+       end do
+    end do
+    deallocate (tmom, ntmom, stat=status)
+    call error_check (status, 'rdjtmom: deallocation error')
+  end subroutine rdjtmom
+
+  subroutine get_tm (spn, smx, nt)
+! obtain expanded array of transition moments, tm for spin spn
+    use rmx1_in, only: bug5
+    integer, intent(in)       :: spn       ! target spin
+    integer, intent(in)       :: smx       ! max multipole
+    integer, intent(in)       :: nt        ! # targets
+    integer                   :: i, j, lp, ntm, l, ier, status, lb
+
+    do i = 1, ntmb
+       if (.not.associated(tmb(i)%tms)) exit
+       if (tmb(i)%spn /= spn) cycle ! not spin required
+       ntm = tmb(i)%ntm      ! # transition moments for spin value
+       allocate (tm(nt*(nt+1)/2,(smx+1)/2), stat=status)
+       if (status /= 0) STOP
+       tm = 0.0_wp
+       do j = 1, ntm
+          l = tmb(i)%lmds(j)
+          lb = tmb(i)%lbls(j)
+          tm(lb, (l+1)/2) = tmb(i)%tms(j)
+       end do
+       if (bug5 > 0) then
+          do j = 1, ntm
+             l = tmb(i)%lmds(j)
+             lb = tmb(i)%lbls(j)
+             write (fo,'(3i6,a,f14.6)') j, lb, l, ' tms = ', tmb(i)%tms(j)
+          end do
+       end if
+       return  ! required tm set has been located
+    end do
+    write (fo,'(a,i6)') 'get_tm: spin not found, spn = ', spn
+    STOP
+  end subroutine get_tm
+
+end module potl_cofs

pfarm/cpu/src_mpi_omp_atom/precisn.f90

+module precisn
+! define precision-related paramters
+! Time-stamp: "98/03/25 08:30:46 cjn"
+  implicit none
+  public
+  integer, parameter    :: ibyte = 4
+  integer, parameter    :: rbyte = 8
+  integer, parameter    :: zbyte = 16
+  integer, parameter    :: sp=selected_real_kind(6)
+  integer, parameter    :: wp=selected_real_kind(12)
+  integer, parameter    :: ep1=selected_real_kind(19)
+  integer, parameter    :: ip_long = selected_int_kind(18)
+  integer, parameter    :: ep = MAX(ep1, wp)
+!  real(wp), parameter   :: acc8 = epsilon(1.0_wp)
+!  real(wp), parameter   :: acc16 = epsilon(1.0_ep)
+!  real(wp), parameter   :: fpmax = huge(1.0_wp) * acc8
+!  real(wp), parameter   :: fpmin = tiny(1.0_wp) / acc8
+!  real(wp), save        :: fplmax
+!  real(wp), save        :: fplmin
+  real(wp), parameter   :: acc8 = 2.0e-16_wp
+  real(wp), parameter   :: acc16 = 3.0e-33_ep
+  real(wp), parameter   :: fpmax = 1.0e60_wp
+  real(wp), parameter   :: fpmin = 1.0e-60_wp
+  real(wp), parameter   :: fplmax = 140.0_wp
+  real(wp), parameter   :: fplmin = -140.0_wp
+end module precisn

pfarm/cpu/src_mpi_omp_atom/rd_Hfile.f90

+module hfile_data
+! read and store target and channel data from HAM H-file
+! Time-stamp: "2003-10-03 17:28:25 cjn"
+! reading of jpi_coupling data included april to dec 2006 vmb
+  use precisn, only: wp
+  use io_units, only: fo
+  use mpi_params, only: io_processor
+  implicit none
+
+  integer, save                :: iunit = 77
+
+  integer, save                :: nelc, nz, lrang2, lamax, ntarg
+  real(wp), save               :: rmatr, bbloch
+  integer, allocatable, save   :: jtarg(:), ltarg(:), starg(:), ptarg(:)
+  real(wp), allocatable, save  :: etarg(:)
+  integer, save                :: nc, nc1, nc2, ns
+  integer, allocatable, save   :: lchl(:), tchl(:), kschl(:)
+  real(wp), allocatable, save  :: ethr(:)
+
+  private
+  public readh1, readh2, cparm, degtar, asc_order
+  public readhj1, readhj2
+  public get_spins, get_split, split_chls, reorder_chls
+  public nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
+  public jtarg, ltarg, starg, ptarg, etarg
+  public lchl, tchl, kschl, nc, ns, ethr, nc1, nc2
+
+contains
+
+  subroutine cparm
+!   cparm : channel parameters eth, lch
+    real(wp)                :: e0      ! initial target energy
+    integer                  :: i, status
+
+    allocate (ethr(nc), stat=status)
+    e0 = MINVAL(etarg(1:ntarg))
+    do i = 1, nc
+       ethr(i) = etarg(tchl(i)) - e0
+    end do
+  end subroutine cparm
+
+  subroutine asc_order (en, nc, en_order)
+!  define pointer array eorder, giving values of en in ascending order
+    integer, intent(in)   :: nc           ! no scattering channels
+    real(wp), intent(in)  :: en(nc)       ! array to be ordered
+    integer, intent(out)  :: en_order(nc) ! order pointers to en
+    logical               :: mask(nc)     ! automatic mask array
+    integer :: i, j(1)
+
+    mask(1:nc) = .true.
+    do i = 1, nc
+       j = MINLOC(en,mask)
+       en_order(i) = j(1)
+       mask(j(1)) = .false.
+    end do
+  end subroutine asc_order
+
+  subroutine readh1
+! Read part of asymptotic file that is independent of lrgl, spin, parity
+    use xdr_files, only: open_xdr, xdr_io
+    use rmx1_in, only: filh, xdr, bug1, bug2
+    use scaling, only: set_charge, scale_radius, scale_etarg
+    integer                 :: ihbuf(5)
+    real(wp)                :: rhbuf(2), etgr
+    integer                 :: i, n, l, status
+    real(wp), allocatable   :: cfbut(:)
+    real(wp)                :: e0, scale
+
+
+    if(io_processor) write (fo,'(a)') 'Reading H file: ' // TRIM(filh)
+
+    if (xdr) then
+       iunit = open_xdr (file=TRIM(filh), action='read')
+       call xdr_io (iunit, ihbuf, 5)
+       nelc = ihbuf(1)
+       nz = ihbuf(2)
+       lrang2 = ihbuf(3)
+       lamax = ihbuf(4)
+       ntarg = ihbuf(5)
+       call xdr_io (iunit, rhbuf, 2)
+       rmatr = rhbuf(1)
+       bbloch = rhbuf(2)
+
+       allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), &
+            ptarg(ntarg), stat=status)
+       call xdr_io (iunit, etarg, ntarg)
+       call xdr_io (iunit, ltarg, ntarg)
+       call xdr_io (iunit, starg, ntarg)
+       call xdr_io (iunit, ptarg, ntarg)
+       allocate (cfbut(3), stat=status)
+       do l = 1, lrang2
+          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
+       end do
+    else          ! normal fortran unformatted output
+       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
+       allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), &
+            ptarg(ntarg), cfbut(3*lrang2),  stat=status)
+       read (iunit) etarg(1:ntarg)
+       read (iunit) ltarg(1:ntarg)
+       read (iunit) starg(1:ntarg)
+       read (iunit) ptarg(1:ntarg)
+       read (iunit) cfbut
+    end if
+
+    if ((bug2 > 0).and.io_processor) then
+       write (fo,'(a)') 'Subroutine readh1: sector matrix setup'
+       write (fo,'(a)') 'Data on the H file:'
+       write (fo,'(5i5,2f12.6)') nelc, nz, lrang2, lamax, ntarg, &
+            rmatr, bbloch
+       write (fo,'(5f14.6)') etarg
+       write (fo,'(5i6)') ltarg
+       write (fo,'(5i6)') starg
+       write (fo,'(5i6)') ptarg
+    end if
+
+! process input data:
+    call set_charge (nz, nelc)
+    call scale_radius (rmatr)
+    call scale_etarg (etarg)     ! etarg now in scaled Ryd
+    if ((bug1 > 0).and.io_processor) then
+      write (fo,'(/,a,/)') 'Target states'
+      write (fo,'(10x,a,5x,a,3x,a,8x,a,/,43x,a)') 'index', 'total l', &
+         '(2*s+1)', 'energy', 'scaled ryd'
+    end if
+    e0 = etarg(1)
+    do i = 1, ntarg
+       etgr = (etarg(i) - e0)           ! convert to Rydbergs
+       if ((bug1 > 0).and.io_processor) write (fo,'(3x,3i10,3x,f12.6)') i, ltarg(i), starg(i), etgr
+    end do
+    deallocate (cfbut, stat=status)
+  end subroutine readh1
+
+  subroutine readhj1
+! Read part of asymptotic file that is independent of lrgl, spin, parity
+! here lrgl corresponds to 2J and spin is zero
+    use xdr_files, only: open_xdr, xdr_io
+    use rmx1_in, only: filh, xdr, bug1, bug2
+    use scaling, only: set_charge, scale_radius, scale_etarg
+    integer                 :: ihbuf(5)
+    real(wp)                :: rhbuf(2), etgr
+    integer                 :: i, n, l, status
+    real(wp), allocatable   :: cfbut(:)
+    real(wp)                :: e0, scale
+
+
+    if(io_processor) write (fo,'(a)') 'Reading H file: ' // TRIM(filh)
+
+    if (xdr) then
+       iunit = open_xdr (file=TRIM(filh), action='read')
+       call xdr_io (iunit, ihbuf, 5)
+       nelc = ihbuf(1)
+       nz = ihbuf(2)
+       lrang2 = ihbuf(3)
+       lamax = ihbuf(4)
+       ntarg = ihbuf(5)
+       call xdr_io (iunit, rhbuf, 2)
+       rmatr = rhbuf(1)
+       bbloch = rhbuf(2)
+
+       allocate (etarg(ntarg), jtarg(ntarg), ltarg(ntarg), &
+            starg(ntarg), ptarg(ntarg), stat=status)
+       call xdr_io (iunit, etarg, ntarg)
+       call xdr_io (iunit, jtarg, ntarg)
+       call xdr_io (iunit, ltarg, ntarg)
+       call xdr_io (iunit, starg, ntarg)
+       call xdr_io (iunit, ptarg, ntarg)
+       allocate (cfbut(3), stat=status)
+       do l = 1, lrang2
+          call xdr_io (iunit, cfbut, 3)   ! skip Buttle data for rmx1
+       end do
+    else          ! normal fortran unformatted output
+       read (iunit) nelc, nz, lrang2, lamax, ntarg, rmatr, bbloch
+       allocate (etarg(ntarg), jtarg(ntarg), ltarg(ntarg), starg(ntarg), &
+            ptarg(ntarg), cfbut(3*lrang2),  stat=status)
+       read (iunit) etarg(1:ntarg)
+       read (iunit) jtarg(1:ntarg)
+       read (iunit) ltarg(1:ntarg)
+       read (iunit) starg(1:ntarg)
+       read (iunit) ptarg(1:ntarg)
+       read (iunit) cfbut
+    end if
+
+    if ((bug2 > 0).and.io_processor) then
+       write (fo,'(a)') 'Subroutine readhj1: sector matrix setup'
+       write (fo,'(a)') 'Data on the H file:'
+       write (fo,'(5i5,2f12.6)') nelc, nz, lrang2, lamax, ntarg, &
+            rmatr, bbloch
+       write (fo,'(5f14.6)') etarg
+       write (fo,'(5i6)') jtarg
+       write (fo,'(5i6)') ltarg
+       write (fo,'(5i6)') starg
+       write (fo,'(5i6)') ptarg
+    end if
+
+! process input data:
+    call set_charge (nz, nelc)
+    call scale_radius (rmatr)
+    call scale_etarg (etarg)     ! etarg now in scaled Ryd
+    if ((bug1 > 0).and.io_processor) then
+      write (fo,'(/,a,/)') 'Target states'
+      write (fo,'(10x,a,5x,a,3x,a,3x,a,8x,a,/,43x,a)') 'index', &
+        '  2j   ', 'total l', '(2*s+1)', 'energy', 'scaled ryd'
+    end if
+    e0 = etarg(1)
+    do i = 1, ntarg
+       etgr = (etarg(i) - e0)           ! convert to Rydbergs
+       if((bug1 > 0).and.io_processor) write (fo,'(3x,4i10,3x,f12.6)') i, jtarg(i), ltarg(i), starg(i), etgr
+    end do
+    deallocate (cfbut, stat=status)
+  end subroutine readhj1
+
+  subroutine get_spins (st1, st2)
+! find target spins which contribute to scattering SLp symmetry
+    integer, intent(out)   :: st1, st2 ! target spins for SLp case
+    integer                :: i
+
+    st2 = -999
+    st1 = starg(tchl(1))  ! 1st target spin
+    do i = 2, nc
+       if (starg(tchl(i)) == st1) cycle
+       st2 = starg(tchl(i))
+       exit
+    end do
+  end subroutine get_spins
+
+  subroutine split_chls_old (schl1, schl2)
+! split the channels into two possible spins
+    integer, intent(out)   :: schl1(nc), schl2(nc) ! channel seqs
+    integer                :: ispch(nc) ! new channel order
+    integer                :: st1, st2
+    integer                :: i, nts
+
+    st1 = starg (tchl(1))    ! first target spin
+    nts = 1
+    nc1 = 1
+    nc2 = 0
+    schl1(1) = 1
+    channels: do i = 2, nc
+       if (starg(tchl(i)) == st1) then ! first target channel
+          nc1 = nc1 + 1
+          schl1(nc1) = i
+       else
+          if (nts == 1) then   ! found a second target spin
+             nts = nts + 1
+             st2 = starg(tchl(i))
+             nc2 = 1
+             schl2(nc2) = i
+          else if (starg(tchl(i)) == st2) then ! second target channel
+             nc2 = nc2 + 1
+             schl2(nc2) = i
+          else
+             write (fo,'(a)') 'split_chls: Channel data inconsistent'
+             write (fo,'(a,2i6)') 'st1, st2, i, s = ', st1, st2, i, &
+                  starg(tchl(i))
+             stop
+          end if
+       end if
+    end do channels
+    if (nc /= nc1 + nc2) then
+       write (fo,'(a,3i4)') 'split_chls: error, nc1, nc2, nc = ', &
+            nc1, nc2, nc
+       stop
+    end if
+  end subroutine split_chls_old
+
+  subroutine get_split (st1, st2)
+! find splitting parameters which contribute to scattering symmetry
+! in SLP symmetries splitting parameters correspond to target spin
+! in JP symmetries splitting parameters correspond to channel q no k
+    integer, intent(out)   :: st1, st2 ! splitting parameters
+    integer                :: i
+
+    st2 = -999
+    st1 = kschl(1)  ! 1st splitting parameter
+    do i = 2, nc
+       if (kschl(i) == st1) cycle
+       st2 = kschl(i)
+       exit
+    end do
+  end subroutine get_split
+
+  subroutine split_chls (schl1, schl2)
+! this is changed but should be applicable to LS and J
+! split the channels into two possible spins
+    integer, intent(out)   :: schl1(nc), schl2(nc) ! channel seqs
+    integer                :: ispch(nc) ! new channel order
+    integer                :: st1, st2
+    integer                :: i, nts
+
+    st1 = kschl(1)    ! first ks value, either target spin or k
+    nts = 1
+    nc1 = 1
+    nc2 = 0
+    schl1(1) = 1
+    channels: do i = 2, nc
+       if (kschl(i) == st1) then ! first channel split
+          nc1 = nc1 + 1
+          schl1(nc1) = i
+       else
+          if (nts == 1) then   ! found a second ks value
+             nts = nts + 1
+             st2 = kschl(i)
+             nc2 = 1
+             schl2(nc2) = i
+          else if (kschl(i) == st2) then ! second channel split
+             nc2 = nc2 + 1
+             schl2(nc2) = i
+          else
+             write (fo,'(a)') 'split_chls: Channel data inconsistent'
+             write (fo,'(a,2i6)') 'st1, st2, i, ks = ', st1, st2, i, &
+                    kschl(i)
+             stop
+          end if
+       end if
+    end do channels
+    if (nc /= nc1 + nc2) then
+       write (fo,'(a,3i4)') 'split_chls: error, nc1, nc2, nc = ', &
+            nc1, nc2, nc
+       stop
+    end if
+    write (fo,*) ' nc1,nc2 ', nc1,nc2
+    write (fo,*) 'schl1 ', (schl1(i),i=1,nc1)
+    write (fo,*) 'schl2 ', (schl2(i),i=1,nc2)
+  end subroutine split_chls
+
+  subroutine reorder_chls
+! order channel arrays lchl, tchl, kschl according to split parameters
+    integer, allocatable      :: tmp_lchl(:), tmp_tchl(:), tmp_kschl(:)
+    integer, allocatable      :: schl1(:), schl2(:)
+    integer                   :: i, j, ip, jp, nc2, nct, status
+
+    allocate (schl1(nc), schl2(nc), stat=status)
+    call split_chls (schl1, schl2)
+    if (nc1 == nc) then ! only one target spin for this SLp case
+       deallocate (schl1, schl2, stat=status)
+       return
+    end if
+    nc2 = nc - nc1
+    allocate (tmp_lchl(nc), tmp_tchl(nc), tmp_kschl(nc), stat=status)
+    tmp_lchl = lchl
+    tmp_tchl = tchl
+    tmp_kschl = kschl
+    do i = 1, nc1
+       ip = schl1(i)
+       lchl(i) = tmp_lchl(ip)
+       tchl(i) = tmp_tchl(ip)
+       kschl(i) = tmp_kschl(ip)
+    end do
+    do i = 1, nc2
+       jp = nc1 + i
+       ip = schl2(i)
+       lchl(jp) = tmp_lchl(ip)
+       tchl(jp) = tmp_tchl(ip)
+       kschl(jp) = tmp_kschl(ip)
+    end do
+    deallocate (tmp_lchl, tmp_tchl, tmp_kschl, stat=status)
+  end subroutine reorder_chls
+
+  subroutine degtar (degeny)
+! make target states truly degenerate according to degeny (scaled Rydbergs)
+! etarg is in scaled Rydbergs and in energy order
+    real(wp), intent(in) :: degeny  ! degeneracy criterion
+    real(wp)             :: etdeg(ntarg) ! tmp modified target energies
+    integer               :: i, iloop, j, nd
+    real(wp)             :: eav, esum, etargi
+
+    i = 1
+    do iloop = 1, ntarg
+       etargi = etarg(i)
+       esum = etargi
+       nd = 1
+       do j = i+1, ntarg
+          if ((etarg(j)- etargi) > degeny) exit
+          esum = esum + etarg(j)
+          nd = nd + 1
+       end do
+       eav = esum / nd ! average of near degenerate states
+       do j = i, i+nd-1
+          etdeg(j) = eav
+       end do
+       i = i + nd
+       if (i > ntarg) exit
+    end do
+    etarg = etdeg
+  end subroutine degtar
+
+  subroutine readh2 (lrgl1, nspn1, npty1)
+! read H-file data for a particular scattering symmetry
+    use xdr_files, only: xdr_io
+    use rmx1_in, only: xdr, bug3
+    integer, intent(in)      :: lrgl1   ! scattering orb a.m.
+    integer, intent(in)      :: nspn1   ! scattering spin multiplicity
+    integer, intent(in)      :: npty1   ! scattering parity
+    integer, allocatable     :: nltarg(:)
+    real(wp), allocatable    :: wmat(:), eig(:)
+    integer                  :: loc, lrgl2, nspn2, npty2
+    integer                  :: ibuf(6)
+    integer                  :: i, j, k, it, nt, n, more2, nl, status
+    integer                  :: ncs
+    logical                  :: match
+
+    more2 = 1
+    match = .false.
+    nt = ntarg
+    symmetries: do while (more2 == 1) ! search for required SLp
+       if (xdr) then
+          call xdr_io (iunit, ibuf, 6)
+          lrgl2 = ibuf(1);   nspn2 = ibuf(2);  npty2 = ibuf(3)
+          nc = ibuf(4);  ns = ibuf(5);   more2 = ibuf(6)
+          ncs = MAX(nc, ns)
+          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
+                    eig(ns), wmat(ncs), stat=status)
+          call xdr_io (iunit, nltarg, nt)
+          call xdr_io (iunit, lchl, nc)
+          call xdr_io (iunit, eig, ns)  ! skip over eig, wmat
+          do j = 1, ns
+             call xdr_io (iunit, loc)   ! state number
+             call xdr_io (iunit, wmat(1:nc), nc)
+          end do
+       else
+          read (iunit) lrgl2, nspn2, npty2, nc, ns, more2
+          ncs = MAX(nc, ns)
+          allocate (lchl(nc), nltarg(nt), tchl(nc), kschl(nc), &
+                    eig(ns), wmat(ncs), stat=status)
+          read (iunit) nltarg
+          read (iunit) lchl
+          read (iunit) eig
+          do j = 1, ns
+             read (iunit) loc
+             read (iunit) wmat(1:nc)
+          end do
+       end if
+
+! form tchl array:
+       i = 0
+       do it  = 1, nt
+          nl = nltarg(it)
+          do n = 1, nl
+             i = i + 1
+             tchl(i) = it
+          end do
+       end do
+
+! form kschl array:  ! this array will be used for channel splitting
+       do i = 1, nc
+          kschl(i) = starg(tchl(i))
+       end do
+
+       if (lrgl1 == lrgl2 .and. nspn1 == nspn2 .and. npty1 == npty2) &
+            then
+          match = .true.
+          deallocate (eig, wmat, nltarg, stat=status)
+          exit
+       end if
+       deallocate (lchl, tchl, kschl, eig, wmat, nltarg, stat=status)
+       if (more2 == 2) exit
+       exit
+    end do symmetries
+
+    if (.NOT.match) then
+       write (fo,'(a,3i5)') 'readh2: unmatched symmetry, SLp = ', &
+            nspn1, lrgl1, npty1
+       stop
+    end if
+    if (bug3 /= 0) then
+       write (fo,'(a,i6,a,i6)') 'nc = ', nc, '    ns = ', ns
+       write (fo,'(/,a)') 'Channel data:'
+       write (fo,'(a)') '         target sequence   orbital a.m.'
+       do i = 1, nc
+          write (fo,'(i4,15x,i4,5x,i4)') i, tchl(i), lchl(i)
+       end do
+    end if
+  end subroutine readh2
+
+  subroutine readhj2 (lrgl1, nspn1, npty1)
+! read H-file data for a particular scattering symmetry
+! here lrgl corresponds to 2J and spin multiplicity is zero
+    use xdr_files, only: xdr_io
+    use rmx1_in, only: xdr, bug3
+    integer, intent(in)      :: lrgl1   ! scattering orb a.m.=2J or -2K
+    integer, intent(inout)   :: nspn1   ! 0 or weighting factor info for K-file
+    integer, intent(in)      :: npty1   ! scattering parity
+    integer, allocatable     :: nltarg(:)
+    real(wp), allocatable    :: wmat(:), eig(:)
+    integer                  :: loc, lrgl2, nspn2, npty2
+    integer                  :: ibuf(6)
+    integer                  :: i, j, k, it, nt, n, more2, nl, status
+    integer                  :: ncs
+    logical                  :: match
+
+    more2 = 1
+    match = .false.
+    nt = ntarg
+    symmetries: do while (more2 == 1) ! search for required SLp
+       if (xdr) then
+          call xdr_io (iunit, ibuf, 6)
+          lrgl2 = ibuf(1);   nspn2 = ibuf(2);  npty2 = ibuf(3)
+          nc = ibuf(4);  ns = ibuf(5);   more2 = ibuf(6)
+          ncs = MAX(nc, ns)
+          allocate (lchl(nc), kschl(nc), nltarg(nt), tchl(nc), eig(ns),& 
+                   wmat(ncs), stat=status)
+
+          call xdr_io (iunit, nltarg, nt)
+          call xdr_io (iunit, lchl, nc)
+          call xdr_io (iunit, kschl, nc)
+          call xdr_io (iunit, eig, ns)  ! skip over eig, wmat
+          do j = 1, ns
+             call xdr_io (iunit, loc)   ! state number
+             call xdr_io (iunit, wmat(1:nc), nc)
+          end do
+       else
+          read (iunit) lrgl2, nspn2, npty2, nc, ns, more2
+          ncs = MAX(nc, ns)
+          allocate (lchl(nc), kschl(nc), nltarg(nt), tchl(nc), eig(ns),& 
+                   wmat(ncs), stat=status)
+
+          read (iunit) nltarg
+          read (iunit) lchl
+          read (iunit) kschl
+          read (iunit) eig
+          do j = 1, ns
+             read (iunit) loc
+             read (iunit) wmat(1:nc)
+          end do
+       end if
+
+! form tchl array:
+       i = 0
+       do it  = 1, nt
+          nl = nltarg(it)
+          do n = 1, nl
+             i = i + 1
+             tchl(i) = it
+          end do
+       end do
+
+       if (lrgl1 == lrgl2 .and. npty1 == npty2) then !nspn1 not defined by namelist
+          match = .true.
+          deallocate (eig, wmat, nltarg, stat=status)
+          exit
+       end if
+       nspn1 = nspn2 !use the value on the h-file (weighting factor for K-files)
+       deallocate (lchl, kschl, tchl, eig, wmat, nltarg, stat=status)
+       if (more2 == 2) exit
+       exit
+    end do symmetries
+
+    if (.NOT.match) then
+       write (fo,'(a,3i5)') 'readhj2: unmatched symmetry, SLp = ', &
+            nspn1, lrgl1, npty1
+       stop
+    end if
+    if (bug3 /= 0) then
+       write (fo,'(a,i6,a,i6)') 'nc = ', nc, '    ns = ', ns
+       write (fo,'(/,a)') 'Channel data:'
+       write (fo,'(a)') '  target index       l          2k'
+       do i = 1, nc
+          write (fo,'(i4,6x,i4,5x,i4,5x,i4)') i, tchl(i), lchl(i),&
+                 kschl(i)
+       end do
+    end if
+  end subroutine readhj2
+
+end module hfile_data

pfarm/cpu/src_mpi_omp_atom/rmx1.f90

+program rmx1
+  !  Electron-atom R-matrix external region calculation - STAGE 1
+  ! Time-stamp: "2003-09-23 17:06:43 cjn"
+  ! modified for jpi_coupling april to dec 2006 vmb
+  use precisn, only: wp, ip_long
+  use io_units, only: fo
+  use slp, only: qno, def_qno
+  use rmx1_in, only: targ_low_sce, targ_high_sce, rd_nmlst, split_prop,&
+       degeny, lrgl1, nspn1, npty1, ne, bug1,       &
+       start_mesh, jpi_coupling
+  use hfile_data, only: readh1, degtar, readh2, get_spins, &
+       reorder_chls, cparm, nc, nc1, nc2, &
+       readhj1, readhj2, get_split
+  use energy_grid, only: energy_max
+  use potl_cofs, only: rdtmom, rdjtmom
+  use pdg_ctl, only: diag
+  use error_out, only: error_check
+  use mpi_params, only: taskid, numtasks, io_processor
+!  use mpi
+
+  !  use blacs, only: nprocs, iam, ictxt, io_processor, p, q, ctxt, myrow,&
+  !       mycol, p0, q0, set_cur_ctxt, p_error
+  implicit none
+  include 'mpif.h'
+  real(wp)                    :: e0
+  real(wp), allocatable, save :: ethr(:) ! threshold energies
+  integer, pointer            :: en_order(:) ! energy-order chl index
+  real(wp)                    :: rafin, rmatr ! prop. limits
+  integer                     :: smx = 1
+  type(qno)                   :: sqno
+  integer                     :: i, j, ion, ncp, ij, status
+  integer                     :: st1, st2
+  real(wp)                    :: t0, t1
+  integer                     :: nmesh, imesh
+  integer                     :: imycol, imyrow, ip, iq
+  integer                     :: BLACS_PNUM
+  integer(ip_long)           :: c0, c1, cr
+  integer, allocatable        :: map(:,:)
+  integer ( kind = 4 )        :: errcode, ierr
+
+  call MPI_INIT( ierr )
+  call MPI_COMM_RANK( MPI_COMM_WORLD, taskid, ierr )
+  call MPI_COMM_SIZE( MPI_COMM_WORLD, numtasks, ierr )
+  if(taskid==0) io_processor=.true.
+
+  call cpu_time (t0)
+  call system_clock (count=c0)
+
+  if (io_processor) then ! This is the i/o processor
+     write (fo,'(//,15x,a)') '============='
+     write (fo,'(15x,a)') 'Program RMX1'
+     write (fo,'(15x,a,//)') '============='
+     write (fo,'(a)') 'R-Matrix External Region: Sector Diagonalization'
+     write (fo,'(a,i6)') 'Number of MPI Tasks = ',numtasks 
+     call flush(fo)  
+  end if
+
+  call rd_nmlst    ! read input namelist
+
+  ! targ_high_sce > 0 signals scattering energies between targets targ_low_sce and
+  ! targ_high_sce. i.e. energies in the resonance region => nmesh = 1
+  ! ne(1) > 0 signals only a single step
+  if (targ_high_sce > 0 .or. ne(1) > 0) then
+     nmesh = 1
+  else
+     nmesh = 2
+  end if
+
+  if (.not. jpi_coupling) then
+     call readh1    !  read first record on H file
+  else
+     call readhj1 
+  end if
+
+  ! set near degenerate target states as strictly degenerate: iterate
+  degeny = 0.1_wp * degeny
+  call degtar (degeny)
+  degeny = 5.0_wp * degeny
+  call degtar (degeny)
+  degeny = 2.0_wp * degeny
+  call degtar (degeny)
+  call degtar (degeny)
+
+  ! define required SLp or Jp case:
+  sqno = def_qno(lrgl1, nspn1, npty1)
+  if(io_processor) then
+     write (fo,'(/a4,i3,a8,i3,a8,i2)')  'S = ', sqno%nspn,&
+          '    L = ', sqno%lrgl, '    P = ', sqno%npty
+  end if
+  if (.not. jpi_coupling) then
+     call readh2 (lrgl1, nspn1, npty1)
+  else
+     call readhj2 (lrgl1, nspn1, npty1)
+  end if
+  call get_split (st1, st2)
+
+  if (split_prop .and. st2 /= -999) then ! there is channel splitting
+     call reorder_chls
+     nc2 = nc - nc1
+     if (io_processor) then
+        write (fo,'(/a/)') 'Channel Splitting Used'
+        write (fo,'(a,i6)') 'Number of spin 1 Channels = ', nc1
+        write (fo,'(a,i6/)') 'Number of Spin 2 channels = ', nc2
+     end if
+  else   ! no channel_splitting
+     if (io_processor) then
+        write (fo,'(/a/)') 'No Channel Splitting Specified'
+        write (fo,'(a,i6)') 'Number of Channels = ', nc
+     end if
+     
+     nc1 = nc
+     split_prop = .false.
+  end if
+  call cparm   ! form ethr using reordered channel data
+  !  end if
+
+  !  call distribute_indata
+  if (.not. jpi_coupling ) then
+     call rdtmom ! read target transition moments and distribute
+     !     call distribute_hdata (st1, st2, lrgl1, nspn1, npty1, start_mesh, &
+     !       nmesh)
+  else
+     call rdjtmom 
+     !     call distribute_hjdata (st1, st2, lrgl1, nspn1, npty1, start_mesh, &
+     !       nmesh)
+  end if
+  sqno = def_qno(lrgl1, nspn1, npty1)
+
+  ! MPI/LAPACK/GPU version only supports nmesh<=2
+  if (nmesh.gt.2) then
+    if(io_processor) write(fo,'(a)') '**** Stopping run - MPI/LAPACK/GPU version does not support runs with nmesh > 2 ****'
+    call MPI_ABORT(MPI_COMM_WORLD,errcode,ierr)
+  end if
+  mesh_loop: do imesh = start_mesh, nmesh
+     call energy_max (imesh, nmesh)   ! find ebig
+     !        call distribute_edata ! broadcast ebig
+
+     ! initialize full or 1st split partition sector data
+     ncp = 1
+     call diag (nc1, ncp, (/st1,st2/), sqno, imesh)
+
+     if (nc1 /= nc) then ! initialize 2nd split partition data
+        if(io_processor) write(fo,'(a,2i0)') ' Initialize 2nd split partition data '
+        nc2 = nc - nc1
+        ncp = nc1 + 1
+        call diag (nc2, ncp, (/st2, -999/), sqno, imesh)
+     end if
+  end do mesh_loop
+
+  call cpu_time (t1)
+  write (fo,'(a,i0,a,f16.4,a)') 'Task ', taskid, ' ; End of RMX1 - Local CPU time     = ', t1 - t0, ' secs'
+  call system_clock (count=c1, count_rate=cr)
+  write (fo,'(a,i0,a,f16.4,a)') 'Task ', taskid, ' ; End of RMX1 - Local Elapsed time = ', REAL(c1-c0,wp) / REAL(cr,wp), ' secs'
+
+  call MPI_BARRIER(MPI_COMM_WORLD,ierr)
+  if (io_processor) then
+     call cpu_time (t1)
+     write (fo,'(a,f16.4,a)') 'End of RMX1 - Global CPU time     = ', t1 - t0, ' secs'
+     call system_clock (count=c1, count_rate=cr)
+     write (fo,'(a,f16.4,a)') 'End of RMX1 - Global Elapsed time = ', REAL(c1-c0,wp) / &
+          REAL(cr,wp), ' secs'
+  end if
+  call MPI_FINALIZE(ierr)
+
+end program rmx1