!!!!!!!!!!!!!!!!!!!!!?????????????????????????
!!!!!!!!!!!!!!!!!!! one line change here just  ????????
!!!!!!!!!!!!!!!!!for using subroutine outgrads1 ????
!!!!!!!!!!!!!!!!!!!!!?????????????????????????
subroutine antest_maps0(mype,theta0f,z0f)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    antest_maps0
!   prgmmr:
!
! abstract: this routine creates output maps of background error correlations
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2010-03-30  zhu    - use nvars from control_vectors
!   2010-12-05  pondeca - change variable names: "tv" is now "t" , and "st" is "sf"
!   2013-10-19  todling - metguess now holds background
!   2016-04-19  pondeca - updated to use bundle_work as argument to ansmoothrf
!
!   input argument list:
!    mype
!    theta0f
!    z0f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: kvar_start,kvar_end,var_names,pf2aP1,indices
  use gridmod, only: nsig,nsig1o,nlon,nlat,istart,jstart,lat2,lon2
  use constants, only: zero_single,zero,one,rd_over_cp,r100
  use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
  use guess_grids, only: ntguessig,ges_prsl
  use fgrid2agrid_mod, only: fgrid2agrid
  use control_vectors, only: nvars
  use gsi_bundlemod, only: gsi_bundlecreate
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetvar
  use gsi_bundlemod, only: gsi_bundleputvar
  use gsi_bundlemod, only: gsi_bundledestroy
  use gsi_bundlemod, only: gsi_grid
  use gsi_bundlemod, only: gsi_gridcreate
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use mpeu_util, only: die
  implicit none

  integer(i_kind),intent(in   ) :: mype
  type(gsi_bundle):: bundle_work
  type(gsi_grid) :: grid
  character(2) :: names2dwork(1),names3dwork(4)

  real(r_single) ,intent(in   ) :: theta0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: z0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)

  character(len=*),parameter::myname='antest_maps0'
  real(r_kind),dimension(lat2,lon2,nsig):: twork,qwork,stwork,vpwork
  real(r_kind),dimension(lat2,lon2):: pwork

  real(r_kind) tempf(nlat,nlon),tempc(pf2aP1%nlatf,pf2aP1%nlonf)
  real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat)

  character(80) ref_plotcor
  character(80) var_plotcor
  integer(i_kind) i_plotcor,j_plotcor,k_plotcor

  real(r_kind),dimension(:,:  ),pointer::ges_z_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()

  real(r_kind)h00,h000
  integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot
  integer(i_kind) it,mm1,ier,istatus
  integer(i_kind) iloc,jloc,kloc

  integer(i_kind):: iref

  integer(i_kind):: nsig_aux

  names2dwork(1)='ps'
  names3dwork(1)='sf'
  names3dwork(2)='vp'
  names3dwork(3)='t'
  names3dwork(4)='q'

!*********************************************************************
!          variable names expected for var_plotcor are
!
!    sf  -- stream function
!    vp  -- velocity potential
!    ps  -- surface pressure
!    t   -- virtual temperature
!    q   -- specific humidity
!    note: add more as desired
!*********************************************************************
! Make choices here!
  i_plotcor=nlat/2
  j_plotcor=nlon/2
  k_plotcor=min(nsig,25)
  iloc=i_plotcor-istart(mype+1)+2
  jloc=j_plotcor-jstart(mype+1)+2
  kloc=k_plotcor
  var_plotcor='sf'
!Note: Must call this subroutine from anprewgt_reg.f90
!Make sure statement has been uncommented!
! End of choice section
!*********************************************************************

  call gsi_gridcreate(grid,lat2,lon2,nsig)
  ref_plotcor='theta'
  ier=0
  it=ntguessig
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
  ier=ier+istatus
  if(ier/=0) call die(myname,'missing fields, ier= ', ier)

  do 200 iref=1,5

     if (iref==1) var_plotcor='sf'
     if (iref==2) var_plotcor='vp'
     if (iref==3) var_plotcor='ps'
     if (iref==4) var_plotcor='t'
     if (iref==5) var_plotcor='q'

     ref_plotcor=var_plotcor
     lunin=1
     if(mype==0) then
        open(lunin,file="cormaps_"//trim(var_plotcor),form='unformatted')
        rewind lunin
     end if
     ivar_plot=0
     do ivar=1,nvars
        if(trim(var_names(ivar))==trim(var_plotcor)) then
           ivar_plot=ivar
           exit
        end if
     end do
     if(ivar_plot==0) then
        write(0,*)' in antest_maps0, variable ',trim(var_plotcor),'  not found.  program stops'
        call mpi_finalize(ierror)
        stop
     end if

     twork=zero
     pwork=zero
     qwork=zero
     stwork=zero
     vpwork=zero

     if(ivar_plot==1) stwork(iloc,jloc,kloc)=one
     if(ivar_plot==2) vpwork(iloc,jloc,kloc)=one
     if(ivar_plot==3) pwork(iloc,jloc)=one
     if(ivar_plot==4) twork(iloc,jloc,kloc)=one
     if(ivar_plot==5) qwork(iloc,jloc,kloc)=one

     call gsi_bundlecreate(bundle_work,grid,'bundle work',ier, &
              names2d=names2dwork,names3d=names3dwork,bundle_kind=r_kind)

     call gsi_bundleputvar(bundle_work,'ps',pwork,ier)
     call gsi_bundleputvar(bundle_work,'sf',stwork,ier)
     call gsi_bundleputvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundleputvar(bundle_work,'t' ,twork,ier)
     call gsi_bundleputvar(bundle_work,'q' ,qwork,ier)

     call ansmoothrf(bundle_work)

     call gsi_bundlegetvar(bundle_work,'ps',pwork,ier)
     call gsi_bundlegetvar(bundle_work,'sf',stwork,ier)
     call gsi_bundlegetvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundlegetvar(bundle_work,'t' ,twork,ier)
     call gsi_bundlegetvar(bundle_work,'q' ,qwork,ier)

     call gsi_bundledestroy(bundle_work,ier)

     if(mype==0) write(0,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), &
            '  i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

     if(ivar_plot==1) h00=stwork(iloc,jloc,kloc)
     if(ivar_plot==2) h00=vpwork(iloc,jloc,kloc)
     if(ivar_plot==3) h00=pwork(iloc,jloc)
     if(ivar_plot==4) h00=twork(iloc,jloc,kloc)
     if(ivar_plot==5) h00=qwork(iloc,jloc,kloc)

     call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror)

     stwork=stwork/h000
     vpwork=vpwork/h000
     pwork=pwork/h000
     twork=twork/h000
     qwork=qwork/h000

!     output original pot temp  (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=ges_tv_it(i,j,k)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!             output "smoothed pot temp"

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=theta0f(i,j,k-indices%kps+1)
              end do
           end do
           call fgrid2agrid(pf2aP1,tempc,tempf)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     if(ivar_plot==3) then
        nsig_aux=1
      else
        nsig_aux=nsig
     endif

     mm1=mype+1
     do k=1,nsig_aux
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              if(ivar_plot==1) outwork(jglob,iglob)=stwork(i,j,k)
              if(ivar_plot==2) outwork(jglob,iglob)=vpwork(i,j,k)
              if(ivar_plot==3) outwork(jglob,iglob)=pwork(i,j)
              if(ivar_plot==4) outwork(jglob,iglob)=twork(i,j,k)
              if(ivar_plot==5) outwork(jglob,iglob)=qwork(i,j,k)
           end do
        end do
!             very slow way to move field from local processor to processor 0

        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!     output original terrain (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     outwork=zero_single
     do j=2,lon2-1
        jglob=jstart(mm1)-2+j
        do i=2,lat2-1
           iglob=istart(mm1)-2+i
           outwork(jglob,iglob)=ges_z_it(i,j)
        end do
     end do
     call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
     if(mype==0) write(lunin) outwork0

!             output "smoothed terrain"

     if(mype==0)PRINT*,'IN ANPREWGT_REG,KPS,KPE=',indices%KPS,indices%KPE
     do k=1, 1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=z0f(i,j,k)!  theta0f(i,j,k-indices%kps+1)
              end do
           end do
           call fgrid2agrid(pf2aP1,tempc,tempf)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     close(lunin)
!    if(mype>-1000) then
!       call mpi_finalize(i)
!       stop
!    end if

200 continue

end subroutine antest_maps0
!-------------------------------------------------------------------------------------

subroutine antest_maps0_subdomain_option(mype,theta0f,z0f)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    antest_maps0_subdomain_option
!   prgmmr:
!
! abstract: this routine creates output maps of background error correlations
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2013-01-15  parrish - ansmoothrf_reg_subdomain_option now uses bundle input.
!                             Add changes to update from old to new version.
!   2013-10-19  todling - metguess now holds background
!
!   input argument list:
!    mype
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use constants, only: max_varname_length
  use anberror, only: kvar_start,kvar_end,var_names,levs_jdvar,indices,pf2aP1
  use gridmod, only: nsig,nlon,nlat,istart,jstart,lat2,lon2
  use constants, only: zero_single,zero,one,rd_over_cp,r100
  use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
  use guess_grids, only: ntguessig,ges_prsl
  use control_vectors, only: nvars
  use gsi_bundlemod, only: gsi_bundlecreate
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetvar
  use gsi_bundlemod, only: gsi_bundleputvar
  use gsi_bundlemod, only: gsi_bundledestroy
  use gsi_bundlemod, only: gsi_grid
  use gsi_bundlemod, only: gsi_gridcreate
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use mpeu_util, only: die

  implicit none

  integer(i_kind),intent(in   ) :: mype
  real(r_single) ,intent(in   ) :: theta0f(lat2,lon2,nsig)
  real(r_single) ,intent(in   ) :: z0f(lat2,lon2,nsig)

  character(len=*),parameter::myname='antest_maps0_subdomain_option'
  real(r_kind),dimension(lat2,lon2,nsig):: twork,qwork,stwork,vpwork
  real(r_kind),dimension(lat2,lon2):: pwork
  real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat)
  character(80) ref_plotcor
  character(80) var_plotcor
  character(80) plotname
  integer(i_kind) i_plotcor,j_plotcor,k_plotcor
  integer(i_kind) iloc,jloc,kloc,ier,istatus

  real(r_kind),dimension(:,:  ),pointer::ges_z_it =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_it=>NULL()

  real(r_kind) h00,h000
  integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot
  integer(i_kind) it,mm1
  integer(i_kind) lvar
  integer(i_kind):: ips,ipe,jps,jpe,kps,kpe
  integer(i_kind):: nlatf,nlonf
  type(gsi_bundle):: bundle_work
  type(gsi_grid) :: grid
  character(max_varname_length) :: names2dwork(1),names3dwork(4)

  names2dwork(1)='ps'
  names3dwork(1)='sf'
  names3dwork(2)='vp'
  names3dwork(3)='t'
  names3dwork(4)='q'
  ips=indices%ips; ipe=indices%ipe
  jps=indices%jps; jpe=indices%jpe
  kps=indices%kps; kpe=indices%kpe
  nlatf=pf2aP1%nlatf
  nlonf=pf2aP1%nlonf

!*********************************************************************
!          variable names expected for var_plotcor are
!
!    sf  -- stream function
!    vp  -- velocity potential
!    ps  -- surface pressure
!    t   -- virtual temperature
!    q   -- specific humidity
!    oz  -- ozone
!    sst -- sea surface temperature
!    stl -- skin temp over land
!    sti -- skin temp over ice
!    cw  -- cloud water
!*********************************************************************
! Make choices here!
  i_plotcor=460 !440!  430
  j_plotcor=265!  250
  k_plotcor=1
  iloc=i_plotcor-istart(mype+1)+2
  jloc=j_plotcor-jstart(mype+1)+2
!  var_plotcor='sf'
!Note: Must call this subroutine from anprewgt_reg.f90
!Make sure statement has been uncommented!
! End of choice section
!*********************************************************************

     call gsi_gridcreate(grid,lat2,lon2,nsig)
  ref_plotcor='theta'
  ier=0
  it=ntguessig
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
  ier=ier+istatus
  if(ier/=0) call die(myname,'missing fields, ier= ', ier)

  do 200 lvar=1,5
     if (lvar==1)  var_plotcor='sf'
     if (lvar==2)  var_plotcor='vp'
     if (lvar==3)  var_plotcor='ps'
     if (lvar==4)  var_plotcor='t'
     if (lvar==5)  var_plotcor='q'
     lunin=1
     if(mype==0) then
!       open(lunin,file='cormaps_'//trim(var_plotcor),form='unformatted')
!       rewind lunin
     end if
     ivar_plot=0
     do ivar=1,nvars
        if(trim(var_names(ivar))==trim(var_plotcor)) then
           ivar_plot=ivar
           exit
        end if
     end do
     if(ivar_plot==0) then
        write(6,*)' in antest_maps0, variable ',trim(var_plotcor),'  not found.  program stops'
        call mpi_finalize(ierror)
        stop
     end if
     kloc=k_plotcor
     twork=zero
     pwork=zero
     qwork=zero
     stwork=zero
     vpwork=zero
     if(i_plotcor>=ips.and.i_plotcor<=ipe.and. &
        j_plotcor>=jps.and.j_plotcor<=jpe) then
        if(ivar_plot==1)     stwork(iloc,jloc,kloc)=one
        if(ivar_plot==2) vpwork(iloc,jloc,kloc)=one
        if(ivar_plot==3) pwork(iloc,jloc)=one
        if(ivar_plot==4) twork(iloc,jloc,kloc)=one
        if(ivar_plot==5) qwork(iloc,jloc,kloc)=one
     end if

     call gsi_bundlecreate(bundle_work,grid,'bundle work',ier, &
              names2d=names2dwork,names3d=names3dwork,bundle_kind=r_kind)
     call gsi_bundleputvar(bundle_work,'ps',pwork,ier)
     call gsi_bundleputvar(bundle_work,'sf',stwork,ier)
     call gsi_bundleputvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundleputvar(bundle_work,'t' ,twork,ier)
     call gsi_bundleputvar(bundle_work,'q' ,qwork,ier)

     call ansmoothrf_reg_subdomain_option(bundle_work)
     call gsi_bundlegetvar(bundle_work,'ps',pwork,ier)
     call gsi_bundlegetvar(bundle_work,'sf',stwork,ier)
     call gsi_bundlegetvar(bundle_work,'vp',vpwork,ier)
     call gsi_bundlegetvar(bundle_work,'t' ,twork,ier)
     call gsi_bundlegetvar(bundle_work,'q' ,qwork,ier)
     call gsi_bundledestroy(bundle_work,ier)
 !   if(mype==0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, &
 !                nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1
     if(mype==0) write(6,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), &
            '  i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

    !---------------------in case we haven't normalized, divide by value of correlation point
     h00=zero
     if(i_plotcor>=ips.and.i_plotcor<=ipe.and. &
        j_plotcor>=jps.and.j_plotcor<=jpe) then
        if(ivar_plot==1)     h00=stwork(iloc,jloc,kloc)
        if(ivar_plot==2) h00=vpwork(iloc,jloc,kloc)
        if(ivar_plot==3) h00=pwork(iloc,jloc)
        if(ivar_plot==4) h00=twork(iloc,jloc,kloc)
        if(ivar_plot==5) h00=qwork(iloc,jloc,kloc)
     end if
     call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
     stwork=stwork/h000
     vpwork=vpwork/h000
     pwork=pwork/h000
     twork=twork/h000
     qwork=qwork/h000

!     output original pot temp  (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=ges_tv_it(i,j,k)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'theta')
     end do

!             output "smoothed pot temp"

     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=theta0f(i,j,k)
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'sm_theta')
     end do

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        kloc=levs_jdvar(k)
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              if(ivar_plot==1)     outwork(jglob,iglob)=stwork(i,j,kloc)
              if(ivar_plot==2) outwork(jglob,iglob)=vpwork(i,j,kloc)
              if(ivar_plot==3) outwork(jglob,iglob)=pwork(i,j)
              if(ivar_plot==4) outwork(jglob,iglob)=twork(i,j,kloc)
              if(ivar_plot==5) outwork(jglob,iglob)=qwork(i,j,kloc)
           end do
        end do
!             very slow way to move field from local processor to processor 0

        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) then
           write(plotname,'("sub_",a)')trim(var_plotcor)
           call outgrads1(outwork0,nlon,nlat,plotname)
        end if
     end do

!     output original terrain (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     outwork=zero_single
     do j=2,lon2-1
        jglob=jstart(mm1)-2+j
        do i=2,lat2-1
           iglob=istart(mm1)-2+i
           outwork(jglob,iglob)=ges_z_it(i,j)
        end do
     end do
     call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
     !if(mype==0) write(lunin) outwork0
     if(mype==0) call outgrads1(outwork0,nlon,nlat,'z')

!             output "smoothed terrain"

     if(mype == 0)PRINT*,'IN ANPREWGT_REG_subdomain_option,KPS,KPE=',KPS,KPE

     do k=1,1
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=z0f(i,j,k)
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        !if(mype==0) write(lunin) outwork0
        if(mype==0) call outgrads1(outwork0,nlon,nlat,'sm_z')
     end do

     !close(lunin)
200 continue

end subroutine antest_maps0_subdomain_option

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

subroutine outgrads1(f,nx,ny,label)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    outgrads1
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-18  lueken - added subprogram doc block
!   2012-12-11  parrish - assign np a value.
!
!   input argument list:
!    label
!    nx,ny
!    f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  use kinds, only: i_kind,r_single
  implicit none

  character(*)   ,intent(in   ) :: label
  integer(i_kind),intent(in   ) :: nx,ny
  real(r_single) ,intent(in   ) :: f(nx,ny)

  integer(i_kind) i,l,next,last,np,ntime,ioutdat,ioutcor,koutmax
  real(r_single) rlonmap0,undef,dlonmap,pinc,startp,rlatmap0,dlatmap
  character(80) dsdes,dsdat
  character(80) datdes(1000)
  character(1) blank
  data blank/' '/
  data undef/-9.99e33_r_single/

  ioutcor=10
  ioutdat=11
  np=1

  write(dsdes,'(a,".ctl")')trim(label)
  write(dsdat,'(a,".grd")')trim(label)
  open(unit=ioutcor,file=dsdes,form='formatted')
  open(unit=ioutdat,file=dsdat,form='unformatted')
  ntime=1
  rlonmap0=1._r_single
  dlonmap=1._r_single
  rlatmap0=1._r_single
  dlatmap=1._r_single
  startp=1._r_single
  pinc=1._r_single
  koutmax=1
  do i=1,1000
     write(datdes(i),'(80a1)')(blank,l=1,80)
  end do
  write(datdes(1),'("DSET ",a)')trim(dsdat)
  write(datdes(2),'("options big_endian sequential")')
  write(datdes(3),'("TITLE ",a)')trim(label)
  write(datdes(4),'("UNDEF ",e11.2)')undef
  write(datdes(5),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nx,rlonmap0,dlonmap
  write(datdes(6),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')ny,rlatmap0,dlatmap
  next=7
  write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc
  next=next+1
  write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax
  next=next+1
  write(datdes(next),'("VARS 1")')
  next=next+1
  write(datdes(next),'("f   ",i5," 99 f   ")')np
  next=next+1
  write(datdes(next),'("ENDVARS")')
  last=next
  write(ioutcor,'(a80)')(datdes(i),i=1,last)
  close(ioutcor)

  write(ioutdat) f
  close(ioutdat)

end subroutine outgrads1


text: Address(Balsam Place) is a residential Culdesacs culdesac in East Davis that branches off of Temple Drive. It is north of Slide Hill Park.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!=======================================================================
      SUBROUTINE MULTBC(N, K, M, C, TEMP, B)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!=======================================================================
!       called by: DVDRVR
!
!       Multiplies B(N,K)*C(K,M) and stores it in B(N,M)
!       Used for collapsing the expanding basis to current estimates,
!       when basis becomes too large, or for returning the results back
!
!       Subroutines called
!       DINIT, DGEMV, DCOPY
!-----------------------------------------------------------------------
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:22  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE dgemv_I
      USE dcopy_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: N
      INTEGER  :: K
      INTEGER  :: M
      REAL(DOUBLE)  :: C(K*M)
      REAL(DOUBLE)  :: TEMP(M)
      REAL(DOUBLE)  :: B(N*K)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IROW
!-----------------------------------------------
!-----------------------------------------------------------------------
      DO IROW = 1, N
!              CALL DINIT(M,0.d0,TEMP,1)
         CALL DGEMV ('Transp', K, M, 1.D0, C, K, B(IROW), N, 0.D0, TEMP, 1)
         CALL DCOPY (M, TEMP, 1, B(IROW), N)
      END DO

      RETURN
      END SUBROUTINE MULTBC


text:
C    @(#)srtbrnch.f	20.6 11/11/97
      	subroutine srtbrnch
c
      	include 'ipfinc/parametr.inc'

	include 'ipfinc/blank.inc'
      	include 'ipfinc/branch.inc'
      	include 'ipfinc/bus.inc'
      	include 'ipfinc/prt.inc'
      	include 'ipfinc/qsdup.inc'

        common /scratch/ nbr, array(2,100)
        integer array

        common /is_batch / is_batch

        integer p, pnxt, q, oldp, komparbr, sect, ipack_4
        character id * 1, xbuf * 120

        external komparbr, swap_br

        if (ltot .lt. 1) then
           write (errbuf(1), 100)
  100      format('No branch data in system.')
           if (is_batch .eq. 0) then
              call prterx ('E',1)
              return
           else
              call prterx ('F',1)
              call erexit()
           endif
      	else
           do i = 1, ntot
              nbr = 0
              p = kbsdta(16,i)
              do while (p .gt. 0)
                 nbr = nbr + 1
                 q = brnch_ptr(p)
                 array(1,nbr) = p
                 k2 = ky(p)
                 id = brid(p)
                 sect = brsect(p)
                 if (sect .gt. 0 .and. q .lt. 0) then
                    sect = 10 - sect
                 endif
                 array(2,nbr) = ipack_4 (inp2alf(k2), ichar(id),
     &                                   sect, brtype(p))
                 p = brnch_nxt(p)
              end do
              if (nbr .eq. 0) then
                 if (bus(i) .ne. '~~~~~~~~') then
                   write (errbuf(1), 110) bus(i), base(i)
  110              format('No branches for bus ', a8, f6.1)
                   call prterx ('E',1)
                 endif
              else if (nbr .gt. 1) then
                 dupsw = .false.
                 call shellsrt (1, nbr, komparbr, swap_br)
                 if (dupsw) then
                    j = 1
                    do while (j .lt. nbr)
                       if (array(2,j) .eq. array(2,j+1)) then
                          p = array(1,j)
                          pnxt = array(1,j+1)
                          if (brtype(p) .ne. 1) then
                             write (errbuf(1),122)
  122                        format(' Duplicate branch records.')
                             errbuf(2) = ' '
                             call bcdbrn(p,xbuf)
                             write (errbuf(3),160) xbuf(1:80)
  160                        format(13x,'(',a80,')')
                             call bcdbrn(pnxt,xbuf)
                             write (errbuf(4),160) xbuf(1:80)
                             call prterx ('W',4)
                          endif
                          do k = j+1, nbr-1
                             array(1,k) = array(1,k+1)
                             array(2,k) = array(2,k+1)
                          enddo
                          nbr = nbr - 1
                       else
                          j = j + 1
                       endif
                    enddo
                 endif
                 oldp = array(1,1)
                 kbsdta(16,i) = oldp
                 do j = 2, nbr
                    p = array(1,j)
                    brnch_nxt(oldp) = p
                    oldp = p
                 end do
                 brnch_nxt(oldp) = 0
              endif
           enddo
        endif
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
Orchid Rockrose (Cistus x purpureus)  is a fragrant hybrid of two Mediterranean Cistus species (Cistus ladanifer x Cistus creticus).  This hybrid blooms locally in weather late spring and continues through summer. The rich papery pink blooms are about threeinches across and beautifully marked with a deep magenta eye at the base of each petal contrasted by a dense cluster of deep yellow stamens in the center. The buds sit quietly for several days, but as each flower opens for the day it flattens out before gently dropping each petal in the afternoon heat or curling in slightly by evening.

This relatively fast growing evergreen shrub clothed in narrow dark green foliage grows up to 5 feet tall and roots along its stem  perfect for hillsides. Some well placed specimens can be found spilling onto walkways at the eastern end of the Arboretum, near Solano Park. Ideally suited to Davis, these hardy plants grow in all types of soil, prefer full sun and are drought tolerant. The dry, woody perfume of the rockrose leaves comes from its parent plants, which have been used for centuries in perfumery. Because of its tendency to layer along its stems, Orchid Rockrose is easily propagated by layering or from cuttings. Though it roots easily, this particular hybrid is selected for its compact habit, so its easily tamed for backyard settings.

For a listing of other plants found in Davis, visit our Town Flora.


text:
      DOUBLE PRECISION FUNCTION DK211(I,A,B,X)
      DOUBLE PRECISION F6211, A, B, X, XX, C1, C2, AAJ, C3
      DIMENSION  X(1)
      XX = X(I)
      IF ( (B * XX) ** 2 - A ** 2 ) 100,1,200
    1 CONTINUE
      IF (A .NE. B * XX) GO TO 50
      F6211 = 0.5D0 * (DLOG (DABS(2.0D0 * B * XX)) ) **2
      DK211 = F6211
      RETURN
   50 CONTINUE
      F6211 = 0.0D0
      DK211 = F6211
      RETURN
  100 CONTINUE
      F6211 = DLOG(DABS(A))* DLOG(DABS(XX))
      C1 =-B * XX / A
      C2 = 1.0D0
      J = 0
  110 J = J + 1
      AAJ = J
      C2 = C2 * C1
      C3 = C2 / (AAJ ** 2)
      F6211 = F6211 - C3
      IF(DABS(C3) .GT. 0.1D-5)   GO TO 110
      DK211 = F6211
      RETURN
  200 CONTINUE
      F6211 = (DLOG(DABS(B* XX)) ** 2) / 2.0D0
      C1 =-A / (B * XX)
      C2 = 1.0D0
      J = 0
  210 J = J + 1
      AAJ = J
      C2 = C2 * C1
      C3 = C2 / (AAJ ** 2)
      F6211 = F6211 + C3
      IF(DABS(C3) .GT. 0.1D-5)   GO TO 210
      DK211 = F6211
      RETURN
      END


MODULE interface_tools
  use cosmosis_modules
  IMPLICIT none

  integer, parameter :: dl=8
  type ini_settings
     real(dl):: zmin,zmax,kmin,kmax,dz,dk
     integer ::nk_steps ,nz_steps
  end type ini_settings

  type pk_settings
     real(dl), dimension(:), allocatable :: redshifts
     real(dl), dimension(:), allocatable :: kh
     real(dl), dimension(:,:), allocatable :: matpower
     real(dl), dimension(:,:), allocatable :: ddmat
     integer num_z, num_k
  end type pk_settings

contains

  function load_matter_power(block, PK) result(status)
    use cosmosis_modules
    integer(cosmosis_block) :: block
    integer(cosmosis_status) :: status
    type(pk_settings) :: PK
    real(dl), allocatable, dimension(:) :: k, z
    real(dl), allocatable, dimension(:,:) :: P

    !Get the data columns from the fits data
    status = 0

    !Load k, z, P
    status = datablock_get_double_grid(block, matter_power_lin_section, &
         "K_H", k, "Z", z, "P_K", P)

    if (status .ne. 0) then
       write(*,*) "Could not find K_H, Z, or P_K in block"
       return
    endif

    !Fill in data structure
    PK%num_k = size(k)
    PK%num_z = size(z)
    call allocate_matterpower(PK)
    PK%kh = k
    PK%redshifts = z
    PK%matpower = P

    !Clean up
    deallocate(k, z, P)

  end function load_matter_power

    DOUBLE PRECISION FUNCTION dz_growth(z,xa,ya,ya2) ! growth function fitting
    IMPLICIT none
    DOUBLE PRECISION, intent(IN) :: z
    DOUBLE PRECISION, allocatable, dimension(:) , intent(IN) :: xa,ya,ya2
    DOUBLE PRECISION :: a,b,hh,x, ypval, yppval
    INTEGER :: jlo,dim
    dim=0
    dim=size(ya, dim=1)
    x=z+0.00000000000000000001d0
    call spline_cubic_val( dim, xa, ya, ya2, z, dz_growth, ypval, yppval)
    return
  END FUNCTION dz_growth

  subroutine allocate_matterpower(PK)
    type(pk_settings) :: PK
    allocate(PK%redshifts(PK%num_z))
    allocate(PK%kh(PK%num_k))
    allocate(PK%matpower(PK%num_k,PK%num_z))
    allocate(PK%ddmat(PK%num_k,PK%num_z))
  end subroutine allocate_matterpower

  subroutine deallocate_matterpower(PK)
    type(pk_settings) :: PK
    deallocate(PK%redshifts)
    deallocate(PK%kh)
    deallocate(PK%matpower)
    deallocate(PK%ddmat)
  end subroutine deallocate_matterpower

END MODULE interface_tools


text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
program MT2rp
!===============================================================================
!  MT2rp outputs the radiation pattern for a moment tensor, including the source
!  polarisation.

   use moment_tensor

   implicit none

   integer, parameter :: rs = 8
   real(rs) :: M(6), azi, inc, P, SV, SH, j
   logical :: dirs_on_stdin = .true., MTs_on_stdin = .true.
   integer :: i, iostat

   ! Get command line parameters and options
   call get_args

   ! Loop over input on stdin
   if (dirs_on_stdin .or. MTs_on_stdin) then
      iostat = 0
      do while (iostat == 0)
         if (dirs_on_stdin .and. .not.MTs_on_stdin) then
            read(*,*,iostat=iostat) azi, inc
         else if (.not.dirs_on_stdin .and. MTs_on_stdin) then
            read(*,*,iostat=iostat) (M(i), i=1,6)
         else
            read(*,*,iostat=iostat) azi, inc, (M(i), i=1,6)
         endif

         ! Check EOF or bad input
         if (iostat < 0) exit
         if (iostat > 0) then
            write(0,'(a)') 'MT2rp: Error: Problem reading values from stdin'
            stop
         endif

         call mt_radiation_pattern(M, azi, inc, P, SV, SH, j)
         call write_output
      enddo

   ! If no lines to process on stdin, do that one
   else
      call mt_radiation_pattern(M, azi, inc, P, SV, SH, j)
      call write_output
   endif

contains
   !============================================================================
   subroutine usage()
   !============================================================================
      write(0,'(a)') &
         'Usage: MT2rp (options) < [Mrr Mtt Mpp Mrt Mrp Mtp] [azi] [inc]', &
         'Output:', &
         '   azi, inc, P, SV, SH, j', &
         'Directions (degrees):', &
         '   inc : Angle from down towards upwards', &
         '   azi : Azimuth measured from north towards east', &
         '   j   : S wave polarisation looking along ray, measured clokwise from upwards', &
         'Options:', &
         '   -d [azi inc] : Read directions from command line; MTs are on stdin', &
         '                  (can be used with -m)', &
         '   -m [Mrr Mtt Mpp Mrt Mrp Mtp] : Read MT from command line; read dirs', &
         '                  from stdin (can be used with -d)'
      stop
   end subroutine usage
   !----------------------------------------------------------------------------

   !============================================================================
   subroutine write_output()
   !============================================================================
      write(*,'(2(f7.2,1x),3(e12.6,1x),f7.2)') azi, inc, P, SV, SH, j
   end subroutine write_output
   !----------------------------------------------------------------------------

   !============================================================================
   subroutine get_args()
   !============================================================================
      implicit none
      integer :: iarg, narg
      character(len=250) :: arg

      narg = command_argument_count()
      iarg = 1
      do while (iarg <= narg)
         call get_command_argument(iarg, arg)
         select case(arg)
            case('-d')
               call get_command_argument(iarg+1, arg)
               read(arg,*) azi
               call get_command_argument(iarg+2, arg)
               read(arg,*) inc
               dirs_on_stdin = .false.
               iarg = iarg + 3
            case('-h')
               call usage
            case('-m')
               call get_mt_from_cmd_line(iarg+1)
               MTs_on_stdin = .false.
               iarg = iarg + 7
            case default
               write(0,'(a)') 'MT2rp: Error: Unrecognised option "'//trim(arg)//'"'
               stop
         end select
      enddo

   end subroutine get_args
   !----------------------------------------------------------------------------

   !===============================================================================
   subroutine get_mt_from_cmd_line(iarg)
   !===============================================================================
      implicit none
      integer, intent(in) :: iarg
      character(len=250) :: arg

      do i=iarg, iarg+5
         call get_command_argument(i, arg)
         read(arg,*,iostat=iostat) M(i-iarg+1)
         if (iostat /= 0) then
            write(0,'(a,i0.0,a)') 'MT2rp: Error: Cannot get value ', i-iarg+1, &
               ' of MT from command line argument "'//trim(arg)//'"'
            stop
         endif
      enddo

   end subroutine get_mt_from_cmd_line
   !-------------------------------------------------------------------------------

end program MT2rp
!-------------------------------------------------------------------------------

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
module mod_cond_triangular
  use mod_triangular
  use mod_prec
  use mod_error_id
  use mod_utility
  implicit none

  private

  public :: lower_left_nullvec, z_lower_left_nullvec, d_lower_left_nullvec, &
       lower_right_nullvec, z_lower_right_nullvec, d_lower_right_nullvec, &
       lower_min_sv, d_lower_min_sv, z_lower_min_sv, &
       lower_max_sv, d_lower_max_sv, z_lower_max_sv, &
       upper_min_sv, d_upper_min_sv, z_upper_min_sv, &
       upper_max_sv, d_upper_max_sv, z_upper_max_sv

  interface lower_left_nullvec
     module procedure z_lower_left_nullvec, d_lower_left_nullvec
  end interface lower_left_nullvec

  interface lower_right_nullvec
     module procedure z_lower_right_nullvec, d_lower_right_nullvec
  end interface lower_right_nullvec

  interface lower_min_sv
     module procedure  d_lower_min_sv, z_lower_min_sv
  end interface lower_min_sv

  interface lower_max_sv
     module procedure  d_lower_max_sv, z_lower_max_sv
  end interface lower_max_sv

  interface upper_min_sv
     module procedure  d_upper_min_sv, z_upper_min_sv
  end interface upper_min_sv

  interface upper_max_sv
     module procedure  d_upper_max_sv, z_upper_max_sv
  end interface upper_max_sv

  integer(kind=int32), parameter :: default_maxit=20
  real(kind=dp), parameter :: default_tol=4*eps, default_tolres=4*eps

contains

  real(kind=dp) function d_upper_max_sv(r,u1,v1,res,tolres0, &
       maxit0,error) result(sigmau)
    real(kind=dp), dimension(:,:), intent(in) :: r
    real(kind=dp), dimension(:), intent(out) :: u1,v1,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit
    real(kind=dp) :: tolres, sigmav, sigmau0, maxr
    type(routine_info), parameter :: info=info_d_upper_max_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(r,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(r,2)) then
       call set_error(2, info, error); return
    else if (n /= size(u1) .or. n /= size(v1) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxr=maxabs(r)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxr*default_tolres,tolres0)

    u1=0.0_dp
    v1=0.0_dp
    if (maxr==0.0_dp) then
       u1(1)=1.0_dp; v1(1)=1.0_dp
       sigmau=0.0_dp
       res=0.0_dp
       call pop_id(error)
       return
    else
       ! Get initial vectors.
       call upper_tr_left_multiply_linpack(r,u1,v1)
       sigmau=norm2(v1)
       v1=v1/sigmau
       sigmav=norm2(u1) ! not really sigmav.
       sigmau=sigmau/sigmav
       u1=u1/sigmav
       k=1
       do while (k < maxit)
          ! Update u1
          res=u1; sigmau0=sigmau
          call upper_left_multiply(r,v1,u1)
          sigmav=norm2(u1)
          u1=u1/sigmav
          ! Update v1
          call upper_tr_left_multiply(r,u1,v1)
          sigmau=norm2(v1)
          v1=v1/sigmau
          res = u1*sigmav-res*sigmau0
          if (norm2(res) < tolres) then
             call pop_id(error)
             return
          end if
          k=k+1
       end do
       if (tolres==0.0_dp) then
          ! maxits is OK if tolerance is zero.
          call pop_id(error)
       else
          call set_error(4, info, error); return
       end if
    end if
  end function d_upper_max_sv

  real(kind=dp) function z_upper_max_sv(r,u1,v1,res,tolres0, &
       maxit0,error) result(sigmau)
    complex(kind=dp), dimension(:,:), intent(in) :: r
    complex(kind=dp), dimension(:), intent(out) :: u1,v1,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit
    real(kind=dp) :: tolres, sigmav, sigmau0, maxr
    type(routine_info), parameter :: info=info_z_upper_max_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(r,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(r,2)) then
       call set_error(2, info, error); return
    else if (n /= size(u1) .or. n /= size(v1) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxr=maxabs(r)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxr*default_tolres,tolres0)

    u1=(0.0_dp,0.0_dp)
    v1=(0.0_dp,0.0_dp)
    if (maxr==0.0_dp) then
       u1(1)=(1.0_dp,0.0_dp); v1(1)=(1.0_dp,0.0_dp)
       sigmau=0.0_dp
       res=(0.0_dp,0.0_dp)
       call pop_id(error)
       return
    else
       ! Get initial vectors.
       call upper_tr_left_multiply_linpack(r,u1,v1)
       sigmau=norm2(v1)
       v1=v1/sigmau
       sigmav=norm2(u1) ! not really sigmav.
       sigmau=sigmau/sigmav
       u1=u1/sigmav
       k=1
       do while (k < maxit)
          ! Update u1
          res=u1; sigmau0=sigmau
          call upper_left_multiply(r,v1,u1)
          sigmav=norm2(u1)
          u1=u1/sigmav
          ! Update v1
          call upper_tr_left_multiply(r,u1,v1)
          sigmau=norm2(v1)
          v1=v1/sigmau
          res = u1*sigmav-res*sigmau0
          if (norm2(res) < tolres) then
             call pop_id(error)
             return
          end if
          k=k+1
       end do
       if (tolres==0.0_dp) then
          ! maxits is OK if tolerance is zero.
          call pop_id(error)
       else
          call set_error(4, info, error); return
       end if
    end if
  end function z_upper_max_sv

  real(kind=dp) function d_upper_min_sv(r,un,vn,res,tolres0, &
       maxit0,error) result(sigmau)
    real(kind=dp), dimension(:,:), intent(in) :: r
    real(kind=dp), dimension(:), intent(out) :: un,vn,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit, p, q
    real(kind=dp) :: tolres, sigmav, maxr
    type(routine_info), parameter :: info=info_d_upper_min_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(r,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(r,2)) then
       call set_error(2, info, error); return
    else if (n /= size(un) .or. n /= size(vn) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxr=maxabs(r)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxr*default_tolres,tolres0)

    un=0.0_dp
    vn=0.0_dp
    if (maxr==0.0_dp) then
       un(1)=1.0_dp; vn(1)=1.0_dp
       sigmau=0.0_dp
       res=0.0_dp
       call pop_id(error)
       return
    else
       ! treat very small diagonal elements as hard zeros.
       p=first_zero_diagonal(r, maxr*eps*eps)
       q=reverse_first_zero_diagonal(r, maxr*eps*eps)
       if (p<=n .or. q >= 1) then
          vn(p)=1.0_dp
          un(q)=1.0_dp
          if (p > 1) call upper_left_invert(r(1:p-1,1:p-1),vn(1:p-1),-r(1:p-1,p))
          if (q < n) call upper_right_invert(un(q+1:n),r(q+1:n,q+1:n),-r(q,q+1:n))
          ! un=conjg(un)
          un=un/norm2(un); vn=vn/norm2(vn)
          sigmau=0.0_dp
          ! Residual is R*vn-sigma*un or just R*vn
          res=0.0_dp
          call upper_left_multiply(r(1:p,1:p),vn(1:p),res(1:p))
          call pop_id(error); return
       else
          ! Get initial vectors.
          call upper_tr_left_invert_linpack(r,un,vn)
          sigmau=1/norm2(un)
          un=sigmau*un
          sigmau=sigmau*norm2(vn)
          k=1
          do while (k < maxit)
             ! Update v
             call upper_left_invert(r,vn,un)
             sigmav=1/norm2(vn)
             vn=sigmav*vn
             ! Update un
             res=un
             call upper_tr_left_invert(r,un,vn)
             sigmau=1/norm2(un)
             un=sigmau*un
             res=res*sigmav - un*sigmau
             if (norm2(res) < tolres) then
                call pop_id(error)
                return
             end if
             k=k+1
          end do
          if (tolres==0.0_dp) then
             ! maxits is OK if tolerance is zero.
             call pop_id(error)
          else
             call set_error(4, info, error); return
          end if
       end if
    end if
  end function d_upper_min_sv

  real(kind=dp) function z_upper_min_sv(r,un,vn,res,tolres0, &
       maxit0,error) result(sigmau)
    complex(kind=dp), dimension(:,:), intent(in) :: r
    complex(kind=dp), dimension(:), intent(out) :: un,vn,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit, p, q
    real(kind=dp) :: tolres, sigmav, maxr
    type(routine_info), parameter :: info=info_z_upper_min_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(r,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(r,2)) then
       call set_error(2, info, error); return
    else if (n /= size(un) .or. n /= size(vn) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxr=maxabs(r)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxr*default_tolres,tolres0)

    un=(0.0_dp,0.0_dp)
    vn=(0.0_dp,0.0_dp)
    if (maxr==0.0_dp) then
       un(1)=(1.0_dp,0.0_dp); vn(1)=(1.0_dp,0.0_dp)
       sigmau=0.0_dp
       res=(0.0_dp,0.0_dp)
       call pop_id(error)
       return
    else
       ! treat very small diagonal elements as hard zeros.
       p=first_zero_diagonal(r, maxr*eps*eps)
       q=reverse_first_zero_diagonal(r, maxr*eps*eps)
       if (p<=n .or. q >= 1) then
          vn(p)=(1.0_dp,0.0_dp)
          un(q)=(1.0_dp,0.0_dp)
          if (p > 1) call upper_left_invert(r(1:p-1,1:p-1),vn(1:p-1),-r(1:p-1,p))
          if (q < n) call upper_right_invert(un(q+1:n),r(q+1:n,q+1:n),-r(q,q+1:n))
          un=conjg(un)
          un=un/norm2(un); vn=vn/norm2(vn)
          sigmau=0.0_dp
          ! Residual is R*vn-sigma*un or just R*vn
          res=(0.0_dp,0.0_dp)
          call upper_left_multiply(r(1:p,1:p),vn(1:p),res(1:p))
          call pop_id(error); return
       else
          ! Get initial vectors.
          call upper_tr_left_invert_linpack(r,un,vn)
          sigmau=1/norm2(un)
          un=sigmau*un
          sigmau=sigmau*norm2(vn)
          k=1
          do while (k < maxit)
             ! Update v
             call upper_left_invert(r,vn,un)
             sigmav=1/norm2(vn)
             vn=sigmav*vn
             ! Update un
             res=un
             call upper_tr_left_invert(r,un,vn)
             sigmau=1/norm2(un)
             un=sigmau*un
             res=res*sigmav - un*sigmau
             if (norm2(res) < tolres) then
                call pop_id(error)
                return
             end if
             k=k+1
          end do
          if (tolres==0.0_dp) then
             ! maxits is OK if tolerance is zero.
             call pop_id(error)
          else
             call set_error(4, info, error); return
          end if
       end if
    end if
  end function z_upper_min_sv

  real(kind=dp) function d_lower_max_sv(l,u1,v1,res,tolres0, &
       maxit0,error) result(sigmau)
    real(kind=dp), dimension(:,:), intent(in) :: l
    real(kind=dp), dimension(:), intent(out) :: u1,v1,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit
    real(kind=dp) :: tolres, sigmav, sigmau0, maxl
    type(routine_info), parameter :: info=info_d_lower_max_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(l,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(l,2)) then
       call set_error(2, info, error); return
    else if (n /= size(u1) .or. n /= size(v1) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxl=maxabs(l)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxl*default_tolres,tolres0)

    u1=0.0_dp
    v1=0.0_dp
    if (maxl==0.0_dp) then
       u1(1)=1.0_dp; v1(1)=1.0_dp
       sigmau=0.0_dp
       res=0.0_dp
       call pop_id(error)
       return
    else
       ! Get initial vectors.
       call lower_tr_left_multiply_linpack(l,u1,v1)
       sigmau=norm2(v1)
       v1=v1/sigmau
       sigmav=norm2(u1) ! not really sigmav.
       sigmau=sigmau/sigmav
       u1=u1/sigmav
       k=1
       do while (k < maxit)
          ! Update u1
          res=u1; sigmau0=sigmau
          call lower_left_multiply(l,v1,u1)
          sigmav=norm2(u1)
          u1=u1/sigmav
          ! Update v1
          call lower_tr_left_multiply(l,u1,v1)
          sigmau=norm2(v1)
          v1=v1/sigmau
          res = u1*sigmav-res*sigmau0
          if (norm2(res) < tolres) then
             call pop_id(error)
             return
          end if
          k=k+1
       end do
       if (tolres==0.0_dp) then
          ! maxits is OK if tolerance is zero.
          call pop_id(error)
       else
          call set_error(4, info, error); return
       end if
    end if
  end function d_lower_max_sv

  real(kind=dp) function z_lower_max_sv(l,u1,v1,res,tolres0, &
       maxit0,error) result(sigmau)
    complex(kind=dp), dimension(:,:), intent(in) :: l
    complex(kind=dp), dimension(:), intent(out) :: u1,v1,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit
    real(kind=dp) :: tolres, sigmav, sigmau0, maxl
    type(routine_info), parameter :: info=info_z_lower_max_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(l,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(l,2)) then
       call set_error(2, info, error); return
    else if (n /= size(u1) .or. n /= size(v1) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxl=maxabs(l)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxl*default_tolres,tolres0)

    u1=(0.0_dp,0.0_dp)
    v1=(0.0_dp,0.0_dp)
    if (maxl==0.0_dp) then
       u1(1)=(1.0_dp,0.0_dp); v1(1)=(1.0_dp,0.0_dp)
       sigmau=0.0_dp
       res=(0.0_dp,0.0_dp)
       call pop_id(error)
       return
    else
       ! Get initial vectors.
       call lower_tr_left_multiply_linpack(l,u1,v1)
       sigmau=norm2(v1)
       v1=v1/sigmau
       sigmav=norm2(u1) ! not really sigmav.
       sigmau=sigmau/sigmav
       u1=u1/sigmav
       k=1
       do while (k < maxit)
          ! Update u1
          res=u1; sigmau0=sigmau
          call lower_left_multiply(l,v1,u1)
          sigmav=norm2(u1)
          u1=u1/sigmav
          ! Update v1
          call lower_tr_left_multiply(l,u1,v1)
          sigmau=norm2(v1)
          v1=v1/sigmau
          res = u1*sigmav-res*sigmau0
          if (norm2(res) < tolres) then
             call pop_id(error)
             return
          end if
          k=k+1
       end do
       if (tolres==0.0_dp) then
          ! maxits is OK if tolerance is zero.
          call pop_id(error)
       else
          call set_error(4, info, error); return
       end if
    end if
  end function z_lower_max_sv

  real(kind=dp) function d_lower_min_sv(l,un,vn,res,tolres0, &
       maxit0,error) result(sigmau)
    real(kind=dp), dimension(:,:), intent(in) :: l
    real(kind=dp), dimension(:), intent(out) :: un,vn,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit, p, q
    real(kind=dp) :: tolres, sigmav, maxl
    type(routine_info), parameter :: info=info_d_lower_min_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(l,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(l,2)) then
       call set_error(2, info, error); return
    else if (n /= size(un) .or. n /= size(vn) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxl=maxabs(l)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxl*default_tolres,tolres0)

    un=0.0_dp
    vn=0.0_dp
    if (maxl==0.0_dp) then
       un(1)=1.0_dp; vn(1)=1.0_dp
       sigmau=0.0_dp
       res=0.0_dp
       call pop_id(error)
       return
    else
       ! treat very small diagonal elements as hard zeros.
       p=first_zero_diagonal(l, maxl*eps*eps)
       q=reverse_first_zero_diagonal(l, maxl*eps*eps)
       if (p<=n .or. q >= 1) then
          un(p)=1.0_dp
          vn(q)=1.0_dp
          if (p > 1) call lower_right_invert(un(1:p-1),l(1:p-1,1:p-1),-l(p,1:p-1))
          if (q < n) call lower_left_invert(l(q+1:n, q+1:n), vn(q+1:n), -l(q+1:n,q))
          un=un/norm2(un); vn=vn/norm2(vn)
          sigmau=0.0_dp
          ! Residual is L*vn-sigma*un or just L*vn
          res=0.0_dp
          call lower_left_multiply(l(q:n,q:n),vn(q:n),res(q:n))
          call pop_id(error); return
       else
          ! Get initial vectors.
          call lower_tr_left_invert_linpack(l,un,vn)
          sigmau=1/norm2(un)
          un=sigmau*un
          sigmau=sigmau*norm2(vn)
          k=1
          do while (k < maxit)
             ! Update v
             call lower_left_invert(l,vn,un)
             sigmav=1/norm2(vn)
             vn=sigmav*vn
             ! Update un
             res=un
             call lower_tr_left_invert(l,un,vn)
             sigmau=1/norm2(un)
             un=sigmau*un
             res=res*sigmav - un*sigmau
             if (norm2(res) < tolres) then
                call pop_id(error)
                return
             end if
             k=k+1
          end do
          if (tolres==0.0_dp) then
             ! maxits is OK if tolerance is zero.
             call pop_id(error)
          else
             call set_error(4, info, error); return
          end if
       end if
    end if
  end function d_lower_min_sv

  real(kind=dp) function z_lower_min_sv(l,un,vn,res,tolres0, &
       maxit0,error) result(sigmau)
    complex(kind=dp), dimension(:,:), intent(in) :: l
    complex(kind=dp), dimension(:), intent(out) :: un,vn,res
    real(kind=dp), intent(in), optional :: tolres0
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in), optional :: maxit0

    integer(kind=int32) :: n, k, maxit, p, q
    real(kind=dp) :: tolres, sigmav, maxl
    type(routine_info), parameter :: info=info_z_lower_min_sv
    !
    sigmau=0.0_dp
    if (failure(error)) return
    call push_id(info, error)

    n=size(l,1)
    if (n<1) then
       call set_error(1, info, error); return
    else if (n /= size(l,2)) then
       call set_error(2, info, error); return
    else if (n /= size(un) .or. n /= size(vn) .or. n /= size(res)) then
       call set_error(3, info, error); return
    end if

    maxl=maxabs(l)
    maxit=equals_option(default_maxit,maxit0)
    tolres=equals_option(maxl*default_tolres,tolres0)

    un=(0.0_dp,0.0_dp)
    vn=(0.0_dp,0.0_dp)
    if (maxl==0.0_dp) then
       un(1)=(1.0_dp,0.0_dp); vn(1)=(1.0_dp,0.0_dp)
       sigmau=0.0_dp
       res=(0.0_dp,0.0_dp)
       call pop_id(error)
       return
    else
       ! treat very small diagonal elements as hard zeros.
       p=first_zero_diagonal(l, maxl*eps*eps)
       q=reverse_first_zero_diagonal(l, maxl*eps*eps)
       if (p<=n .or. q >= 1) then
          un(p)=(1.0_dp,0.0_dp)
          vn(q)=(1.0_dp,0.0_dp)
          if (p > 1) call lower_right_invert(un(1:p-1),l(1:p-1,1:p-1),-l(p,1:p-1))
          if (q < n) call lower_left_invert(l(q+1:n, q+1:n), vn(q+1:n), -l(q+1:n,q))
          un=conjg(un)
          un=un/norm2(un); vn=vn/norm2(vn)
          sigmau=0.0_dp
          ! Residual is L*v-sigma*u or just L*v
          res=(0.0_dp,0.0_dp)
          call lower_left_multiply(l(q:n,q:n),vn(q:n),res(q:n))
          call pop_id(error); return
       else
          ! Get initial vectors.
          call lower_tr_left_invert_linpack(l,un,vn)
          sigmau=1/norm2(un)
          un=sigmau*un
          sigmau=sigmau*norm2(vn)
          k=1
          do while (k < maxit)
             ! Update vn
             call lower_left_invert(l,vn,un)
             sigmav=1/norm2(vn)
             vn=sigmav*vn
             ! Update un
             res=un
             call lower_tr_left_invert(l,un,vn)
             sigmau=1/norm2(un)
             un=sigmau*un
             res=res*sigmav - un*sigmau
             if (norm2(res) < tolres) then
                call pop_id(error)
                return
             end if
             k=k+1
          end do
          if (tolres==0.0_dp) then
             ! maxits shouldn't raise an error if the tolerance is zero.
             call pop_id(error)
          else
             call set_error(4, info, error); return
          end if
       end if
    end if
  end function z_lower_min_sv

  ! null vector of a lower triangular matrix.
  ! error: 1 no null vector within tolerance
  ! p returns a value such that x(p+1:n)=0
  ! If tol=0.0_dp, just do maxits iterations.
  subroutine d_lower_left_nullvec(x,l,tol,maxit, p, error)
    real(kind=dp), dimension(:,:), intent(in) :: l
    real(kind=dp), dimension(:), intent(out) :: x
    real(kind=dp), intent(in) :: tol
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in) :: maxit
    integer(kind=int32), intent(out) :: p

    integer(kind=int32) :: n, k
    real(kind=dp) :: nrmx
    real(kind=dp), dimension(size(x)) :: y
    type(routine_info), parameter :: info=info_d_lower_left_nullvec
    !
    if (failure(error)) return
    call push_id(info, error)

    n=size(l,1)
    x=0.0_dp
    y=0.0_dp
    p=first_zero_diagonal(l, tol)
    If (p==1) then
       x(1)=1.0_dp
       call pop_id(error); return
    else if (p <= n) then
       ! if a small enough diagonal element is found, solve
       ! for a null vector directly.
       x(p)=1.0_dp
       call lower_right_invert(x(1:p-1),l(1:p-1, 1:p-1), -l(p,1:p-1))
       call pop_id(error); return
    else
       ! Otherwise use the Linpack condition estimator approach.
       p=n
       call lower_right_invert_linpack(x,l,y)
       nrmx=norm2(x)
       x=x/nrmx
       nrmx=nrmx/norm2(y)
       k=1
       do while (k < maxit)
          if (1/nrmx < tol) then
             call pop_id(error)
             return
          end if
          call lower_tr_right_invert(y,l,x)
          y=y/norm2(y)
          call lower_right_invert(x,l,y)
          nrmx=norm2(x)
          x=x/nrmx
          k=k+1
       end do
       if (tol==0.0_dp) then
          call pop_id(error)
       else
          call set_error(1, info, error); return
       end if
    end if
  end subroutine d_lower_left_nullvec

  ! null vector of a lower triangular matrix.
  subroutine z_lower_left_nullvec(x,l,tol,maxit,p, error)
    complex(kind=dp), dimension(:,:), intent(in) :: l
    complex(kind=dp), dimension(:), intent(out) :: x
    real(kind=dp), intent(in) :: tol
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in) :: maxit
    integer(kind=int32), intent(out) :: p
    !
    integer(kind=int32) :: n, k
    real(kind=dp) :: nrmx
    complex(kind=dp), dimension(size(x)) :: y
    type(routine_info), parameter :: info=info_z_lower_left_nullvec
    !
    n=size(l,1)
    if (failure(error)) return
    call push_id(info, error)

    x=(0.0_dp, 0.0_dp)
    y=(0.0_dp, 0.0_dp)
    p=first_zero_diagonal(l, tol)
    if (p==1) then
       x(p)=(1.0_dp, 0.0_dp)
       call pop_id(error); return
    else if (p <= n) then
       ! if a small enough diagonal element is found, solve
       ! for a null vector directly.
       x(p)=(1.0_dp, 0.0_dp)
       call lower_right_invert(x(1:p-1),l(1:p-1, 1:p-1), -l(p,1:p-1))
       call pop_id(error); return
    else
       ! Otherwise use the Linpack condition estimator.
       p=n
       call lower_right_invert_linpack(x,l,y)
       nrmx=norm2(x)
       x=x/nrmx
       nrmx=nrmx/norm2(y)
       k=1
       do while (k < maxit)
          if (1/nrmx < tol) then
             call pop_id(error)
             return
          end if
          call lower_tr_right_invert(y,l,x)
          y=y/norm2(y)
          call lower_right_invert(x,l,y)
          nrmx=norm2(x)
          x=x/nrmx
          k=k+1
       end do
       if (tol==0.0_dp) then
          call pop_id(error)
       else
          call set_error(1, info, error); return
       end if
    end if
  end subroutine z_lower_left_nullvec

  ! right null vector of a lower triangular matrix.
  ! Returns p such that x(1:p-1)=0
  subroutine d_lower_right_nullvec(x,l,tol,maxit, p, error)
    real(kind=dp), dimension(:,:), intent(in) :: l
    real(kind=dp), dimension(:), intent(out) :: x
    real(kind=dp), intent(in) :: tol
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in) :: maxit
    integer(kind=int32), intent(out) :: p
    !
    integer(kind=int32) :: n, k
    real(kind=dp) :: nrmx
    real(kind=dp), dimension(size(x)) :: y
    type(routine_info), parameter :: info=info_d_lower_right_nullvec
    !
    n=size(l,1)
    if (failure(error)) return
    call push_id(info, error)
    x=0.0_dp
    y=0.0_dp
    p=reverse_first_zero_diagonal(l, tol)
    if (p==n) then
       x(n)=1.0_dp
       call pop_id(error); return
    else if (p >= 1) then
       ! if a small enough diagonal element is found, solve
       ! for a null vector directly.
       x(p)=1.0_dp
       call lower_left_invert(l(p+1:n,p+1:n),x(p+1:n),-l(p+1:n,p))
       call pop_id(error); return
    else
       ! Otherwise use the Linpack condition estimator approach.
       p=1
       call lower_left_invert_linpack(l,x,y)
       nrmx=norm2(x)
       x=x/nrmx
       nrmx=nrmx/norm2(y)
       k=1
       do while (k < maxit)
          if (1/nrmx < tol) then
             call pop_id(error); return
          end if
          call lower_tr_left_invert(l,y,x)
          y=y/norm2(y)
          call lower_left_invert(l,x,y)
          nrmx=norm2(x)
          x=x/nrmx
          k=k+1
       end do
       if (tol==0.0_dp) then
          call pop_id(error)
       else
          call set_error(1, info, error); return
       end if
    end if
  end subroutine d_lower_right_nullvec

  ! right null vector of a lower triangular matrix.
  subroutine z_lower_right_nullvec(x, l, tol, maxit, p, error)
    complex(kind=dp), dimension(:,:), intent(in) :: l
    complex(kind=dp), dimension(:), intent(out) :: x
    real(kind=dp), intent(in) :: tol
    type(error_info), intent(inout), optional :: error
    integer(kind=int32), intent(in) :: maxit
    integer(kind=int32), intent(out) :: p
    !
    integer(kind=int32) :: n, k
    real(kind=dp) :: nrmx
    complex(kind=dp), dimension(size(x)) :: y
    type(routine_info), parameter :: info=info_z_lower_right_nullvec
    !
    n=size(l,1)
    if (failure(error)) return
    call push_id(info, error)
    x=(0.0_dp, 0.0_dp)
    y=(0.0_dp, 0.0_dp)
    p=reverse_first_zero_diagonal(l, tol)
    if (p==n) then
       x(n)=(1.0_dp, 0.0_dp)
       call pop_id(error); return
    else if (p >= 1) then
       ! if a small enough diagonal element is found, solve
       ! for a null vector directly.
       x(p)=(1.0_dp, 0.0_dp)
       call lower_left_invert(l(p+1:n,p+1:n),x(p+1:n),-l(p+1:n,p))
       call pop_id(error); return
    else
       ! Otherwise use the Linpack condition estimator approach.
       p=1
       call lower_left_invert_linpack(l,x,y)
       nrmx=norm2(x)
       x=x/nrmx
       nrmx=nrmx/norm2(y)
       k=1
       do while (k < maxit)
          if (1/nrmx < tol) then
             call pop_id(error)
             return
          end if
          call lower_tr_left_invert(l,y,x)
          y=y/norm2(y)
          call lower_left_invert(l,x,y)
          nrmx=norm2(x)
          x=x/nrmx
          k=k+1
       end do
       if (tol==0.0_dp) then
          call pop_id(error)
       else
          call set_error(1, info, error); return
       end if
    end if
  end subroutine z_lower_right_nullvec

end module mod_cond_triangular

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


text:
c------------------------------------------------------------------
      subroutine rent(atm,res,rnum,chn,nent)
c
c enter character strings res,atm into hash table for radii
c by assigning it a number nent
c
	include 'qdiffpar4.h'
c
c check to see if there is room
c
	if(irtot.eq.nrlist) then
	  write(6,*)' radii list full- increase nrlist'
	  stop
	end if

	n = ichash(atm,res,rnum,chn)
	if(irnumb(n).ne.0) then
c
c slot filled
c run down linked list
c
9000	continue
	if(irlink(n).eq.0)goto 9001
            n = irlink(n)
	goto 9000
9001	continue
c
c  search for empty slot
c
         new = 1
9002	continue
	if(irnumb(new).eq.0)goto 9003
            new = new + 1
	goto 9002
9003	continue
c
c found one- add link
c
         irlink(n) = new
         n = new
      end if

c
c slot empty
c fill slot
c
      irnumb(n) = nent
      irlink(n) = 0
	irtot = irtot + 1
      return
      end

max_stars_repo_name: caixiuhong/Stable-MCCE

max_issues_repo_name: caixiuhong/Stable-MCCE

max_forks_repo_name: caixiuhong/Stable-MCCE


text:
      SUBROUTINE PUTDSS ( GROUP,LOCAT,PARAM,INTL,TSTYPE,DES,CFSCT,
     . CIRBLK,ZTYPE,UNITS,FACTOR,NSIG,IRND,IFIRST,NXTDAT,DATA,DQ,
     . JDATE,MTIME,IREV,LFNOUT,IOLVL,IFLTAB, NTOTAL )
C
C     PUTDSS ASSEMBLES THE PARTICULAR DATA NEEDED TO WRITE TO DSS
C     CALLS THE APPROPRIATE DSS TIME-SERIES ENTRY ROUTINE TO
C     POST THE DATA TO THE DSS FILE.
C
C     EXTERNAL VARIABLES
      CHARACTER GROUP*(*),LOCAT*(*),TSTYPE*(*),DES*(*),CIRBLK*(*),IREV*1
      CHARACTER DQ(*)*(*),PARAM*(*),ZTYPE*(*),UNITS*(*), CFSCT*(*)
      INTEGER  NXTDAT(*),MTIME(*)
C     INTEGER JDATE(*)                                                  H
      INTEGER*4 JDATE(*), JULS, JULE, INTL                              M
      INTEGER*4 NPERS, NOPERS, inctim, idummy, NVALS
      REAL DATA(*)
C
C     LOCAL VARIABLES
      CHARACTER CBUFF*80
      PARAMETER ( MAXBLK=3000 )
      INTEGER IHEAD(40)
      CHARACTER CPATH*80,CD*32,CE*32
      DIMENSION BUFFER(2*MAXBLK+100),ZDATA(MAXBLK),ZTIME(MAXBLK)
      LOGICAL LPOPT,DEBUG
C
C
C ****************************************************************************
C INITIALIZATION
C ****************************************************************************
C
      DEBUG = LPOPT('D')
C
C     COMMON PATHNAME PARTS
      CALL LFLNB( GROUP,1,LEN(GROUP),IBEG,NA )
      IF ( NA.EQ.0 ) THEN
         NA = 1
      ELSE IF (IBEG.GT.1) THEN
         CBUFF = GROUP
         GROUP = ' '
         GROUP(1:NA) = CBUFF(IBEG:IBEG+NA-1)
      ELSE
      ENDIF
      CALL LFLNB( LOCAT,1,LEN(LOCAT),IBEG,NB )
      IF ( NB.EQ.0 ) THEN
         NB = 1
      ELSE IF (IBEG.GT.1) THEN
         CBUFF = LOCAT
         LOCAT = ' '
         LOCAT(1:NB) = CBUFF(IBEG:IBEG+NB-1)
      ELSE
      ENDIF
      CALL LFLNB( PARAM,1,LEN(PARAM),IBEG,NC )
      IF ( NC.EQ.0 ) THEN
         NC = 1
      ELSE IF (IBEG.GT.1) THEN
         CBUFF = PARAM
          PARAM =  ' '
         PARAM(1:NC) = CBUFF(IBEG:IBEG+NC-1)
      ELSE
      ENDIF
      CALL LFLNB( DES,1,LEN(DES),IBEG,NF )
      IF ( NF.EQ.0 ) THEN
         NF = 1
      ELSE IF (IBEG.GT.1) THEN
         CBUFF = DES
         DES = ' '
         DES(1:NF) = CBUFF(IBEG:IBEG+NF-1)
      ELSE
      ENDIF
C
      IF (DES(1:1).EQ.'*') THEN
         CALL SHSAVK('G',0,KYR,KMO,KDA,KHR,KMN,KSE)
         KMN = 0
         WRITE (DES, 9) CFSCT,KYR,KMO,KDA,KHR,KMN
 9       FORMAT(A,': ',I4.4, '.'I2.2,'.',I2.2,'-',2I2.2)
         CALL CHRLNB(DES, NF)
      ENDIF
C
C     REGULAR INTERVAL TIME-SERIES DATA
      IF ( TSTYPE.EQ.'RG' ) THEN
         J = IFIRST
   10    CONTINUE
         ncount=0
         DO 20 I=1,MAXBLK
            ZDATA(I) = -901.
   20    CONTINUE
         ISTIME = MTIME(J)
         JULS   = JDATE(J)
         IF (DEBUG) WRITE ( LFNOUT,* ) 'JULS =',JULS,'  ISTIME =',ISTIME
   30    CONTINUE
         IF (DEBUG) WRITE ( LFNOUT,* ) 'J =',J,'  JDATE =',JDATE(J),
     1   '   MTIME = ',MTIME(J)
c        write (lfnout,*) 'putdss - intl:     ', intl
         NPERS = NOPERS ( INTL,0,JULS,ISTIME,JDATE(J),MTIME(J) )
         if(debug)write (lfnout,*) 'npers: ', npers
         IF ( NPERS.LT.0 ) THEN
            WRITE ( LFNOUT,*) '*** ERROR *** DATA OUT OF ORDER'
            CBUFF = '    OBSERVATION ON           AT      HRS '
            CALL JULDAT (JDATE(J), 104, CBUFF(20:28), N)
            MI = M2IHM ( MTIME(J), CBUFF(33:36) )
            WRITE ( LFNOUT,* ) CBUFF
            GO TO 999
         ENDIF
         NVALS = 1 + NPERS
         ncount=ncount+1
         if(nvals.lt.ncount) nvals=ncount
         ncount=nvals
         IF ( DEBUG ) WRITE ( LFNOUT,* ) 'NVALS =',NVALS
         IF ( NVALS.LE.MAXBLK ) THEN
            IF ((NSIG.NE.-1).AND.(IRND.NE.-999)) THEN
               CALL ROUND (LFNOUT,DATA(J),1,NSIG,IRND)
            ENDIF
            IF ( ABS(DATA(J)+9999.).GT.1.0E-5 ) THEN
               ZDATA(NVALS) = DATA(J) * FACTOR
            ELSE
            ENDIF
            IETIME = MTIME(J)
            JULE   = JDATE(J)
            J = NXTDAT(J)
            NDATA = NVALS
         ENDIF
         IF ( NVALS.GE.MAXBLK.OR.J.EQ.0 ) THEN
            IF (DEBUG)WRITE (LFNOUT,* ) 'NDATA=',NDATA,'  JULE =',JULE,
     1                 '  IETIME=',IETIME
            IF ( DEBUG ) WRITE (LFNOUT,'(8F10.1)') (ZDATA(I),I=1,NDATA)
            NDATA = MAXBLK
            IF ( IREV.EQ.'R'.OR.LPOPT('R') ) THEN
               ISTAT = 0
            ELSE
               ISTAT = -1
            ENDIF
c
            idum = 0
            dum1 = 0.
            dum2 = 0.
c
c     write (*,*) ' '
c     write (*,*) 'juls:   ', juls
c     write (*,*) 'jule:   ', jule
c     write (*,*) 'istime: ', istime
c     write (*,*) 'ietime: ', ietime
c     write (*,*) 'intl:   ', intl
c     write (*,*) 'ndata:  ', ndata
c
c     do 47 ic=1,3
c        write (*,*) 'zdata(',ic,'): ', zdata(ic)
c  47 continue
c
c     call zset ( 'MLVL', ' ', 9 )
c
            IDUM = 0
      iflag=0
      if(j.eq.0) then
        idummy=inctim(intl,iflag,ncount-1,juls,istime,jule,ietime)
      endif
      if(debug) then
      write(lfnout,*)'j ncount nvals ndata',j,ncount,nvals,ndata
      write(lfnout,*)'juls,istime,jule,ietime',juls,istime,jule,ietime
      endif
            CALL ZPTDTS ( IFLTAB, GROUP, NA, LOCAT, NB, PARAM, NC,
     .                    DES, NF, JULS, ISTIME, JULE, IETIME, INTL,
     .                    DUM1, IDUM, DUM2, ZDATA, NDATA, UNITS,
     .                    ZTYPE, ISTAT )
c
            IF ( ISTAT.GE.10 ) THEN
               WRITE ( LFNOUT,* ) '*** COULD NOT WRITE TO DSS ',
     1            ' DSS ERROR CODE =', ISTAT
            ELSE IF ( ISTAT.GT.0.AND.IOLVL.EQ.1 ) THEN
               WRITE ( LFNOUT,* ) '*** WARNING -- ',
     1         ' DSS ERROR CODE =', ISTAT
            ENDIF
            NTOTAL =NTOTAL + NDATA
         ELSE
            GO TO 30
         ENDIF
         IF ( J.GT.0 ) GO TO 10
C
C     IRREGULAR INTERVAL DATA
C
      ELSE IF ( TSTYPE.EQ.'IR' ) THEN
         NVALS = 1
         J = IFIRST
         BDATE = FLOAT(JDATE(J))
         JULS  = JDATE(J)
            IF (DEBUG) WRITE (LFNOUT,* ) 'JULS =', JULS, '  BDATE = ',
     1     BDATE
   50    CONTINUE
         IF ((NSIG.NE.-1).AND.(IRND.NE.-999)) THEN
            CALL ROUND(LFNOUT,DATA(J),1,NSIG,IRND)
         ENDIF
         IF ( DEBUG ) WRITE (LFNOUT,* ) J,DATA(J)
         IF ( ABS(DATA(J)+9999.).GT.1.0E-5 ) THEN
            ZDATA(NVALS) = DATA(J) * FACTOR
         ELSE
            ZDATA(NVALS) = -901.
         ENDIF
         ZTIME(NVALS)=FLOAT(JDATE(J)-JULS)+FLOAT(MTIME(J))*.0006944444
C        IF ( DEBUG ) WRITE ( LFNOUT,* ) N,ZDATA(N),ZTIME(N)
         IF ( NVALS.LT.MAXBLK.AND.NXTDAT(J).GT.0 ) THEN
            J = NXTDAT(J)
            NVALS = NVALS + 1
            GO TO 50
         ELSE
         ENDIF
            IF ( IREV.EQ.'R'.OR.LPOPT('R') ) THEN
               INFLAG = 1
            ELSE
               INFLAG = 0
            ENDIF
            IF ( DEBUG ) WRITE (LFNOUT,* ) 'INFLAG =',INFLAG
            CALL JULDAT (JULS, 104, CD, N)
            CD(1:2)='01'
            IF ( DEBUG ) WRITE ( LFNOUT,'(1X,20A4)') CD
            CE = CIRBLK
            CALL LFLNB( CIRBLK,1,10,IE,NE )
            IF ( DEBUG ) WRITE ( LFNOUT,'(1X,20A4)') CE
         CPATH=' '
         CALL ZFPN ( GROUP,NA,LOCAT,NB,PARAM,NC,CD,9,CE,NE,DES,NF,
     1               CPATH,NPATH )
         IF ( DEBUG ) WRITE (LFNOUT,'(1X,20A4)') CPATH
         IF ( DEBUG ) WRITE ( LFNOUT,'(1X,4(F10.5,F10.2 ))') ( ZTIME(I),
     1      ZDATA(I),I = 1,NVALS )
         NDATA = MAXBLK
         NBUFF = 2*MAXBLK+100
         IT = 0
         CALL ZPIRTS ( IFLTAB,CPATH,NPATH,BUFFER,NBUFF,IHEAD,IT,ZTIME,
     1   ZDATA,NVALS,BDATE,UNITS,ZTYPE,INFLAG,ISTAT )
         IF ( ISTAT.GE.10 ) THEN
            WRITE ( LFNOUT,* ) '*** COULD NOT WRITE TO DSS ',
     1         ' DSS ERROR CODE =', ISTAT
         ELSE IF ( ISTAT.GT.0.AND.IOLVL.EQ.1 ) THEN
            WRITE ( LFNOUT,* ) '*** WARNING -- ',
     1         ' DSS ERROR CODE =', ISTAT
         ELSE
         NTOTAL = NTOTAL + NVALS
         ENDIF
            IF ( NXTDAT(J).GT.0 ) THEN
               NVALS = 1
               J = NXTDAT(J)
               GO TO 50
            ELSE
            ENDIF
      ELSE
      ENDIF
C
  999 CONTINUE
      RETURN
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
! This code is from QEDMOD with minor modifications.
!
!   V.M. Shabaev, I.I. Tupitsyn, V.A. Yerokhin, QEDMOD: Fortran program for calculating
!   the model Lamb-shift operator, Computer Physics Communications, Volume 189, 2015,
!   Pages 175-181, ISSN 0010-4655, http://dx.doi.org/10.1016/j.cpc.2014.12.002.
!
! For a description of the changes made, see the README.
!
! Original code is distributed under the CPC non-profit use licence
! (http://cpc.cs.qub.ac.uk/licence/licence.html). Redistributed here as part of
! GRASP2K with permission from the authors.
!
! Originally from: f/potgen/wk.f
!
c       =================================================
        subroutine SHAB_wk(maxii,r,v_wk)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c       This subroutine computes Wichmann-Kroll point nuclear
c       vacuum-polarization potential by means of
C       the analytical approximation formulas from
c       A G Fainshtein, N L Manakov and A A Nekipelov,
c       J.Phys.E: At.Mol.Opt.Phys. v.23 (1990) 559-569.
c
c       Input parameters:
c       r(i)    .... Radial grid
c       maxii   .... the size of the r array
c
c       Output parameters:
c       wk(i)   .... Wichmann-Kroll potential multiplied by r(i).
c
c       Common block parameters:
c       cl     ..... Speed of light
c       z      ..... Nuclear charge
c       ii     ..... Number of grid points
c       h, al, bt ..... Parameter of the radial grid
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        implicit real*8 (a-h,o-z)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        common /SHAB_const/cl,z
        real*8 r(maxii),v_wk(maxii)
        integer ii
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c       Calculation Wichmann-Kroll potential times on r
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        call SHAB_wk_init
c
        ii=maxii-20
	Do i=1,ii
          x=r(i)*cl
          v_wk(i)=SHAB_u_wk(x)*r(i)
	End do
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        imax=ii
        do i=ii,1,-1
	  If (dabs(v_wk(i)/r(i)).gt.1.D-40) go to 200
          imax=i
          v_wk(i)=0.d0
        enddo
c       - - - - - - - - - - - - - - - - - - - - - - - - -
 200    v_wk(ii+3)=imax
        v_wk(ii+4)=0.d0
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c        open (unit=15,file='wk.dat',status='unknown')
c        write(15,'(f22.14,16x,a)') z,'z'
c        write(15,'(i22,16x,a)') ii,'ii'
c        write(15,'(d22.14,16x,a)')  h,'h'
c        write(15,'(d22.14,16x,a)') al,'al'
c        write(15,'(d22.14,16x,a)') bt,'bt'
c        write(15,*)
c        write(15,'(17x,a,23x,a,19x)') 'r','wk'
c        do i=1,ii
c          write (15,'(i6,3d24.14)') i,r(i),v_wk(i)
c        enddo
c        write(15,*)
c        write(15,'(d22.14,16x,a)') v_wk(ii+3),'wk(ii+3)'
c        write(15,'(d22.14,16x,a)') v_wk(ii+4),'wk(ii+4)'
c        close(unit=15)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        return
        end
c       =================================================
        subroutine SHAB_wk_init
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        implicit real*8 (a-h,o-z)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c       A G Fainshtein, N L Manakov and A A Nekipelov,
c       J.Phys.E: At.Mol.Opt.Phys. v.23 (1990) 559-569.
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c       Initialization
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        real*8 f,p,q,a
        integer s
        common /SHAB_wk_coeff/f(15,5),p(0:35),q(0:15),
     &  a(4,0:7,4),s(4,4)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        do m=1,15
        do n=1,5
          f(m,n) = 0.d0
        enddo
        enddo
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        p(0)=  2.094 001 968 517 539 d-02
        p(1)= -6.328 682 056 509 706 d-02
        p(2)=  1.656 675 047 469 226 d-01
        p(3)=  6.218 254 052 794 603 d-02
        p(4)=  8.561 450 880 149 488 d-01
        p(5)= -1.324 099 530 525 646 d+00
        p(6)=  1.022 894 990 055 742 d-01
        p(7)=  1.818 708 880 809 046 d-01
        p(8)=  5.002 561 255 243 687 d-04
        p(9)=  2.167 778 174 657 628 d-01
        p(10)= 5.453 835 857 080 859 d-01
        p(11)= 3.137 663 571 113 079 d-01
        p(12)=-1.619 774 880 547 481 d-01
        p(13)=-6.749 903 602 516 955 d-02
        p(14)=-2.935 144 138 820 853 d-04
        p(15)= 9.755 800 734 714 044 d-02
        p(16)=-1.994 016 822 313 688 d-01
        p(17)=-4.777 623 414 220 549 d-02
        p(18)= 7.503 472 520 308 676 d-03
        p(19)= 3.440 742 164 594 281 d-05

        p(20)= 2.829 420 9807 d-03
        p(21)= 1.419 551 8237 d-04
        p(22)= 1.348 983 3978 d-05
        p(23)= 2.013 400 8556 d-07
        p(24)= 8.602 918 4536 d-06
        p(25)=-1.787 016 5724 d-05
        p(26)= 2.326 960 7856 d-05
        p(27)=-1.779 846 3760 d-05
        p(28)= 8.727 474 4996 d-06
        p(29)=-2.923 560 6110 d-06
        p(30)= 6.892 936 0802 d-07
        p(31)=-1.149 169 9126 d-07
        p(32)= 1.330 283 5807 d-08
        p(33)=-1.019 098 8753 d-09
        p(34)= 4.650 804 1566 d-11
        p(35)=-9.578 592 4366 d-13
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        q(0)  = 2.09400223235d-02
        q(1)  = 1.42425723499d-02
        q(2)  = 0.76419630941d-02
        q(3)  = 0.36124746165d-02
        q(4)  = 0.15118743043d-02
        q(5)  = 0.05247168246d-02
        q(6)  = 0.01797758276d-02
        q(7)  =-0.04761316284d-02
        q(8)  = 0.19713924639d-02
        q(9)  =-0.62421312791d-02
        q(10) = 1.29259217859d-02
        q(11) =-1.89507086023d-02
        q(12) = 1.89089237003d-02
        q(13) =-1.24611625668d-02
        q(14) = 0.48793792307d-02
        q(15) =-0.08858438257d-02
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        f(5,1) = 0.3812769d+00
        f(5,2) = 0.0266795d+00
        f(5,3) = 0.0000860d+00

        f(6,1) = 0.4946161d+00
        f(6,2) = 0.0458093d+00
        f(6,3) = 0.0002918d+00
        f(6,4) = 6.d-08

        f(7,1) = 0.6150374d+00
        f(7,2) = 0.0707397d+00
        f(7,3) = 0.0007216d+00
        f(7,4) = 5.d-07

        f(8,1) = 0.7412290d+00
        f(8,2) = 0.1014082d+00
        f(8,3) = 0.0014705d+00
        f(8,4) = 23.0d-07
        f(8,5) = 1.0d-10

        f(9,1) =0.8721898d+00
        f(9,2) =0.1376368d+00
        f(9,3) =0.0026288d+00
        f(9,4) =69.d-07
        f(9,5) =2.0d-09

        f(10,1) =1.0071432d+00
        f(10,2) =0.1791838d+00
        f(10,3) =0.0042776d+00
        f(10,4) =167.0d-07
        f(10,5) =8.0d-09

        f(11,1) = 1.1454776d+00
        f(11,2) = 0.2557770d+00
        f(11,3) = 0.0064866d+00
        f(11,4) = 346.0d-07
        f(11,5) = 3.0d-08

        f(12,1) = 1.2867042d+00
        f(12,2) = 0.2771342d+00
        f(12,3) = 0.0093134d+00
        f(12,4) = 643.0d-07
        f(12,5) = 9.0d-08

        f(13,1) = 1.4304276d+00
        f(13,1) = 0.3329753d+00
        f(13,1) = 0.0128043d+00
        f(13,1) = 1097.0d-07
        f(13,1) = 2.0d-07

        f(14,1) = 1.5763244d+00
        f(14,1) = 0.3930296d+00
        f(14,1) = 0.0169953d+00
        f(14,1) = 1750.0d-07
        f(14,1) = 5.0d-07

        f(15,1) = 1.7241274d+00
        f(15,1) = 0.4570401d+00
        f(15,1) = 0.0219132d+00
        f(15,1) = 2644.0d-07
        f(15,1) = 9.0d-07
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        do i=1,4
        do m=0,7
        do n=1,4
          a(i,m,n)=0.d0
        enddo
        enddo
        enddo
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        a(1,0,1) = 2.9675d+00
        a(1,0,2) = 7.8894d-01
        a(1,0,3) =-3.0201d-04
        a(1,0,4) = 4.9725d-04

        a(1,1,1) =-4.4288d+00
        a(1,1,2) =-2.9946d+01
        a(1,1,3) =-9.0441d-01
        a(1,1,4) = 2.4826d+00

        a(1,2,1) = 5.9339d01
        a(1,2,2) = 8.4881d+02
        a(1,2,3) = 5.8388d+01
        a(1,2,4) =-1.0745d+02

        a(1,3,1) =-5.6820d+02
        a(1,3,2) =-1.4888d+04
        a(1,3,3) =-1.5899d+03
        a(1,3,4) = 2.7128d+03

        a(1,4,1) = 4.6131d+03
        a(1,4,2) = 1.6003d+05
        a(1,4,3) = 2.2880d+04
        a(1,4,4) =-3.8249d+04

        a(1,5,1) =-2.6282d+04
        a(1,5,2) =-1.0259d+06
        a(1,5,3) =-1.8036d+05
        a(1,5,4) = 3.0100d+05

        a(1,6,1) = 8.9448d+04
        a(1,6,2) = 3.5980d+06
        a(1,6,3) = 7.3716d+05
        a(1,6,4) =-1.2368d+06

        a(1,7,1) =-1.3412d+05
        a(1,7,2) =-5.3067d+06
        a(1,7,3) =-1.2228d+06
        a(1,7,4) = 2.0678d+06
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        a(2,0,1) =-6.64673d-01
        a(2,0,2) = 4.68647d-05
        a(2,0,3) = 2.56839d-02

        a(2,1,1) = 5.70540d-01
        a(2,1,2) =-9.83836d-04
        a(2,1,3) =-3.21021d-03

        a(2,2,1) =-9.43955d-01
        a(2,2,2) = 8.28746d-03
        a(2,2,3) = 1.01163d-02

        a(2,3,1) = 1.18405d+00
        a(2,3,2) = 5.94372d-02
        a(2,3,3) =-1.62093d-02

        a(2,4,1) =-8.91689d-01
        a(2,4,2) =-5.31431d-02
        a(2,4,3) = 1.17762d-02

        a(2,5,1) = 3.64707d-01
        a(2,5,2) = 7.52653d-03
        a(2,5,3) =-4.69814d-03

        a(2,6,1) =-6.17902d-02
        a(2,6,2) =-4.54081d-04
        a(2,6,3) = 7.76469d-04
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        a(3,0,1) = 1.93743d-02
        a(3,0,2) = 1.64381d-01
        a(3,0,3) = 2.53250d-02

        a(3,1,1) =-6.43399d-01
        a(3,1,2) =-3.05149d-01
        a(3,1,3) = 3.56651d-03

        a(3,2,1) = 2.56565d-01
        a(3,2,2) = 2.31560d-01
        a(3,2,3) =-8.14980d-03

        a(3,3,1) =-1.09235d-01
        a(3,3,2) =-9.49948d-02
        a(3,3,3) = 6.21168d-03

        a(3,4,1) = 4.14226d-02
        a(3,4,2) = 2.22450d-02
        a(3,4,3) =-3.29527d-03

        a(3,5,1) =-7.28744d-03
        a(3,5,2) =-2.80168d-03
        a(3,5,3) = 7.60885d-04

        a(3,6,1) = 4.54914d-04
        a(3,6,2) = 1.47382d-04
        a(3,6,3) =-6.03993d-05
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        a(4,0,1) =-2.4734758d+09
        a(4,0,2) =-1.5115262d+04

        a(4,1,1) = 3.2270784d+09
        a(4,1,2) = 1.4945586d+04

        a(4,2,1) =-1.7207915d+09
        a(4,2,2) =-5.8361389d+03

        a(4,3,1) = 4.7635740d+08
        a(4,3,2) = 1.1476935d+03

        a(4,4,1) =-7.1414759d+07
        a(4,4,2) =-1.2098290d+02

        a(4,5,1) = 5.4305859d+06
        a(4,5,2) = 6.5439613d+00

        a(4,6,1) =-1.6479928d+05
        a(4,6,2) =-1.4296421d-01
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        s(1,1) = 1
        s(1,2) = 1
        s(1,3) = 0
        s(1,4) = 0

        s(2,1) = 0
        s(2,2) =-3
        s(2,3) = 0
        s(2,4) = 0

        s(3,1) =-1
        s(3,2) = 2
        s(3,3) = 0
        s(3,4) = 0

        s(4,1) =-13
        s(4,2) =-6
        s(4,3) = 0
        s(4,4) = 0
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        return
        end
c       =================================================
        real*8 function SHAB_u_wk(x)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        implicit real*8 (a-h,o-z)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
c       A G Fainshtein, N L Manakov and A A Nekipelov,
c       J.Phys.E: At.Mol.Opt.Phys. v.23 (1990) 559-569.
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        common /SHAB_const/cl,z
        real*8 break (5)
        parameter (pi=3.1415926535 8979323846d0)
        real*8 f,p,q,a
        integer s
        common /SHAB_wk_coeff/f(15,5),p(0:35),q(0:15),
     &  a(4,0:7,4),s(4,4)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        alz=z/cl
        dlam=1.d0-dsqrt(1.d0-alz)*dsqrt(1.d0+alz)
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        if (x.lt.10.d0) then
          i=1
          if (x.gt.4.d0) then
            i=4
          else
            if (x.gt.1.5d0) then
              i=3
            else
              if (x.gt.0.15d0) i=2
            endif
          endif
c
          v3=0.d0
          if (x.le.4.25d0) then
            do n=0,8
              v3=v3+p(n)*x**(n-1)
            enddo
            do n=0,5
              v3=v3+p(n+9)*x**(n+2)*dlog(x)
            enddo
            do n=0,4
              v3=v3+p(n+15)*x**(n+3)*dlog(x)**2
            enddo
          else
            do n=0,15
              v3=v3+p(n+20)*(1.d1/x)**(2*n)
            enddo
            v3=v3/x**5
          endif
          v3=v3/cl*alz**3
c
          SHAB_u_wk=v3
c
          if (x.lt.0.15d0) then
            q5=0.d0
            do n=1,15
              q5=q5+q(n)*dlam**n
            enddo
            q5=q5*alz**3/cl
            SHAB_u_wk=SHAB_u_wk+q5/x
          endif
c
          sum=1.d0
          do n=1,4
          do m=0,7
            sum=sum+a(i,m,n)*(dlam**n)*x**(m+s(i,n))
          enddo
          enddo
          SHAB_u_wk=SHAB_u_wk/sum
c
        else
c
          SHAB_u_wk=2.d0/225.d0/x**5+59.d0/1323.d0/x**7+(1977.d0+
     &    20.d0*alz*alz)/4725.d0/x**9+(1258680.d0+34144.d0*alz*alz)/
     &    190575/x**11+(1960420032.d0+93618070.d0*alz*alz+
     &    96096.d0*alz**4)/12297285.d0/x**13
          SHAB_u_wk=SHAB_u_wk*alz**3
c
          prod=576.d0/x**13
          do m=5,15
            prod=prod*m*m/(x*x)
            sum=0.d0
            do n=1,5
              sum=sum+f(m,n)*alz**(2*n+1)
            enddo
            SHAB_u_wk=SHAB_u_wk+sum*prod
          enddo
          SHAB_u_wk=SHAB_u_wk/(pi*cl)
        endif
c       - - - - - - - - - - - - - - - - - - - - - - - - -
        return
        end

max_stars_repo_name: compas/grasp-rci-qed

max_issues_repo_name: compas/grasp-rci-qed

max_forks_repo_name: compas/grasp-rci-qed


text:
      PROGRAM PRIMES2
C
C     PROGRAM TO PRINT THE FIRST N PRIMES
C
      READ(*,*)N
      I=1
      J=0
 10   I=I+1
      IS = SQRT(FLOAT(I))
      DO 20 K=2,IS
      IF (I/K*K .EQ.I) GOTO 10
 20   CONTINUE
      WRITE(*,*)I
      J=J+1
      IF (J.LT.N) GOTO 10
      END

max_stars_repo_name: felx/detritus

max_issues_repo_name: felx/detritus

max_forks_repo_name: felx/detritus


text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  December 27, 1994
C   Authors:        Joseph W.H. Liu
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C****************     ETREE ..... ELIMINATION TREE     *****************
C***********************************************************************
C***********************************************************************
C
C   WRITTEN BY JOSEPH LIU (JUL 17, 1985)
C
C   PURPOSE:
C       TO DETERMINE THE ELIMINATION TREE FROM A GIVEN ORDERING AND
C       THE ADJACENCY STRUCTURE.  THE PARENT VECTOR IS RETURNED.
C
C   INPUT PARAMETERS:
C       NEQNS           -   NUMBER OF EQUATIONS.
C       (XADJ,ADJNCY)   -   THE ADJACENCY STRUCTURE.
C       (PERM,INVP)     -   PERMUTATION AND INVERSE PERMUTATION VECTORS
C
C   OUTPUT PARAMETERS:
C       PARENT          -   THE PARENT VECTOR OF THE ELIMINATION TREE.
C
C   WORKING PARAMETERS:
C       ANCSTR          -   THE ANCESTOR VECTOR.
C
C***********************************************************************
C
      SUBROUTINE  ETREE (   NEQNS , XADJ  , ADJNCY, PERM  , INVP  ,
     &                      PARENT, ANCSTR                          )
C
C***********************************************************************
C
        INTEGER*4           ADJNCY(*)     , ANCSTR(*)     ,
     &                      INVP(*)       , PARENT(*)     ,
     &                      PERM(*)
C
        INTEGER*4           NEQNS
        INTEGER*4           XADJ(*)
C
C***********************************************************************
C
        INTEGER*4           I     , J     , JSTOP , JSTRT , NBR   ,
     &                      NEXT  , NODE
C
C***********************************************************************
C
        IF  ( NEQNS .LE. 0 )  RETURN
C
        DO  400  I = 1, NEQNS
            PARENT(I) = 0
            ANCSTR(I) = 0
            NODE = PERM(I)
C
            JSTRT = XADJ(NODE)
            JSTOP = XADJ(NODE+1) - 1
            IF  ( JSTRT .LE. JSTOP )  THEN
                DO  300  J = JSTRT, JSTOP
                    NBR = ADJNCY(J)
                    NBR = INVP(NBR)
                    IF  ( NBR .LT. I )  THEN
C                       -------------------------------------------
C                       FOR EACH NBR, FIND THE ROOT OF ITS CURRENT
C                       ELIMINATION TREE.  PERFORM PATH COMPRESSION
C                       AS THE SUBTREE IS TRAVERSED.
C                       -------------------------------------------
  100                   CONTINUE
                            IF  ( ANCSTR(NBR) .EQ. I )  GO TO 300
                            IF  ( ANCSTR(NBR) .GT. 0 )  THEN
                                NEXT = ANCSTR(NBR)
                                ANCSTR(NBR) = I
                                NBR = NEXT
                                GO TO 100
                            ENDIF
C                       --------------------------------------------
C                       NOW, NBR IS THE ROOT OF THE SUBTREE.  MAKE I
C                       THE PARENT NODE OF THIS ROOT.
C                       --------------------------------------------
                        PARENT(NBR) = I
                        ANCSTR(NBR) = I
                    ENDIF
  300           CONTINUE
            ENDIF
  400   CONTINUE
C
        RETURN
      END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv