module scaling
! Time-stamp: "2003-09-18 07:51:08 cjn"
  use precisn, only: wp
  use io_units, only: fo
  implicit none

  private
  public set_charge, charge, get_charge, scale_etarg
  public scale_cf, scale_radius, get_ionicity, scale_cf_farm, scale_cf_farm_pnz

  real(wp), parameter         :: atry = 2.0_wp
  real(wp), save              :: iz_res = -999
  logical, save               :: charged_case
  integer, save               :: ion  ! ionicity

contains

  subroutine set_charge (nz, nelc)
! define overall Coulomb charge on target
    integer, intent(in)      :: nz   ! nuclear charge
    integer, intent(in)      :: nelc ! # target electrons

    iz_res = nz - nelc
    if (iz_res == 0) then
       charged_case = .false.
       ion = 0
    else
       charged_case = .true.
       if (iz_res > 0) then
          ion = 1
       else
          ion = iz_res
       end if
    end if
  end subroutine set_charge

  function get_ionicity () 
    integer :: get_ionicity
    get_ionicity = ion
  end function get_ionicity

  subroutine print_charge
    write (fo,'(a,i6)') 'Residual charge (iz_res) = ', iz_res
  end subroutine print_charge

  subroutine get_charge (ires_charge)
    integer, intent(out)   :: ires_charge

    if (iz_res == -999) then ! not initialized - warn
       write(fo,'(a)') 'get_charge: attempt to access uninitialized &
            &charge'
    else
      ires_charge = iz_res
   end if
  end subroutine get_charge

  subroutine charge (charge_log)
    logical, intent(inout)  :: charge_log

    charge_log = charged_case
  end subroutine charge

  subroutine scale_radius (rmatr)
    real(wp), intent(inout)   :: rmatr

    if (charged_case) rmatr = REAL(iz_res) * rmatr
  end subroutine scale_radius

  subroutine scale_etarg (etarg)
    real(wp),intent(inout)    :: etarg(:)
    real(wp)                  :: fac

    if (charged_case) then
       fac = atry / REAL(iz_res*iz_res,wp)
       etarg = fac * etarg
    else
       fac = atry
       etarg = etarg * fac
    end if
  end subroutine scale_etarg

  subroutine scale_cf (cf,lamd)
    real(wp), intent(inout)    :: cf(:,:)
    integer, intent(in)        :: lamd(:,:)
    real(wp)                   :: flion
    integer                    :: i,j,iplim,jplim

    iplim = size(cf,dim=1)
    jplim = size(cf,dim=2)
    if (charged_case) then
       flion = REAL(iz_res,wp)
       do j = 1, jplim
          do i = 1, iplim
             cf(i,j) = atry * cf(i,j) * flion**(lamd(i,j)-1)
          end do
       end do
    else
       cf = cf * atry
    end if
  end subroutine scale_cf

  subroutine scale_cf_farm(cf_farm)
    real(wp), intent(inout)       :: cf_farm(:,:,:)
    real(wp)                   :: fac,flion
    integer                    :: lamax, lam

! note that the cf_farm are aussumed to include the facotr of two.
    if(charged_case) then
       lamax = size(cf_farm, dim=3)
       flion = REAL(iz_res,wp)
       fac = 1.0_wp
       do lam = 2, lamax
          fac=fac*flion
          cf_farm(:,:,lam) = fac * cf_farm(:,:,lam)
       end do
    end if

  end subroutine scale_cf_farm

  subroutine scale_cf_farm_pnz(cf_nz)
    real(wp), intent(inout)       :: cf_nz(:,:)
    real(wp)                   :: fac,flion
    integer                    :: lamax, lam

! note that the cf_farm are aussumed to include the facotr of two.
    if(charged_case) then
       lamax = size(cf_nz, dim=2)
       flion = REAL(iz_res,wp)
       fac = 1.0_wp
       do lam = 2, lamax
          fac=fac*flion
          cf_nz(:,lam) = fac * cf_nz(:,lam)
       end do
    end if

  end subroutine scale_cf_farm_pnz

end module scaling
