--- /dev/null
+ module io_base
+!-----------------------------------------------------------------------
+ use names
+ use io_units
+ implicit none
+!-----------------------------------------------------------------------------
+! Max. number of AA residues
+ integer,parameter :: maxres=6000!1200
+! Appr. max. number of interaction sites
+ integer,parameter :: maxres2=2*maxres
+! parameter (maxres6=6*maxres)
+! parameter (mmaxres2=(maxres2*(maxres2+1)/2))
+!-----------------------------------------------------------------------------
+! Max. number of S-S bridges
+! integer,parameter :: maxss=20
+!-----------------------------------------------------------------------------
+! Max. number of derivatives of virtual-bond and side-chain vectors in theta
+! or phi.
+! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! readrtns_CSA.F
+!-----------------------------------------------------------------------------
+ subroutine read_bridge
+! Read information about disulfide bridges.
+ use geometry_data, only: nres
+ use energy_data
+ use control_data, only:out1file
+ use MPI_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+! include 'COMMON.DBASE'
+! include 'COMMON.THREAD'
+! include 'COMMON.TIME1'
+! include 'COMMON.SETUP'
+!el local variables
+ integer :: i,j,ierror
+
+! Read bridging residues.
+ read (inp,*) ns
+ if (ns.gt.0) then
+ allocate(iss(ns))
+ read (inp,*) (iss(i),i=1,ns)
+ endif
+
+! print *,'ns=',ns
+ if(me.eq.king.or..not.out1file) &
+ write (iout,*) 'ns=',ns
+ if (ns.gt.0) &
+ write(iout,*) ' iss:',(iss(i),i=1,ns)
+! Check whether the specified bridging residues are cystines.
+ do i=1,ns
+ if (itype(iss(i)).ne.1) then
+ if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
+ 'Do you REALLY think that the residue ',&
+ restyp(itype(iss(i))),i,&
+ ' can form a disulfide bridge?!!!'
+ write (*,'(2a,i3,a)') &
+ 'Do you REALLY think that the residue ',&
+ restyp(itype(iss(i))),i,&
+ ' can form a disulfide bridge?!!!'
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,ierror)
+ stop
+#endif
+ endif
+ enddo
+! Read preformed bridges.
+ if (ns.gt.0) then
+ read (inp,*) nss
+ if (nss.gt.0) then
+ if(.not.allocated(ihpb)) allocate(ihpb(nss))
+ if(.not.allocated(jhpb)) allocate(jhpb(nss))
+ read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
+
+ if(.not.allocated(dhpb)) allocate(dhpb(nss))
+ if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
+ if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
+ if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
+! el Initialize the bridge array
+ do i=1,nss
+ dhpb(i)=0.0D0
+ enddo
+ endif
+!--------------------
+ if(fg_rank.eq.0) &
+ write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+ if (nss.gt.0) then
+ write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+ nhpb=nss
+! Check if the residues involved in bridges are in the specified list of
+! bridging residues.
+ do i=1,nss
+ do j=1,i-1
+ if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
+ .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
+ write (iout,'(a,i3,a)') 'Disulfide pair',i,&
+ ' contains residues present in other pairs.'
+ write (*,'(a,i3,a)') 'Disulfide pair',i,&
+ ' contains residues present in other pairs.'
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,ierror)
+ stop
+#endif
+ endif
+ enddo
+ do j=1,ns
+ if (ihpb(i).eq.iss(j)) goto 10
+ enddo
+ write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+ 10 continue
+ do j=1,ns
+ if (jhpb(i).eq.iss(j)) goto 20
+ enddo
+ write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+ 20 continue
+ dhpb(i)=d0cm
+ forcon(i)=akcm
+ enddo
+ do i=1,nss
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ enddo
+ endif
+ endif
+! write(iout,*) "end read_bridge"
+ return
+ end subroutine read_bridge
+!-----------------------------------------------------------------------------
+ subroutine read_x(kanal,*)
+
+ use geometry_data
+ use energy_data
+ use geometry, only:int_from_cart1
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! Read coordinates from input
+!
+!el local variables
+ integer :: l,k,j,i,kanal
+
+ read(kanal,'(8f10.5)',end=10,err=10) &
+ ((c(l,k),l=1,3),k=1,nres),&
+ ((c(l,k+nres),l=1,3),k=nnt,nct)
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ call int_from_cart1(.false.)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ endif
+ enddo
+
+ return
+ 10 return 1
+ end subroutine read_x
+!-----------------------------------------------------------------------------
+ subroutine read_threadbase
+
+ use geometry_data
+ use energy_data
+ use compare_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+! include 'COMMON.DBASE'
+! include 'COMMON.THREAD'
+! include 'COMMON.TIME1'
+
+!el local variables
+ integer :: k,j,i
+
+! Read pattern database for threading.
+ read (icbase,*) nseq
+ allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
+ allocate(nres_base(3,nseq)) !(3,maxseq)
+ allocate(str_nam(nseq)) !(maxseq)
+ do i=1,nseq
+ read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
+ nres_base(2,i),nres_base(3,i)
+ read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
+ nres_base(1,i))
+! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
+! & nres_base(2,i),nres_base(3,i)
+! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
+! & nres_base(1,i))
+ enddo
+ close (icbase)
+ if (weidis.eq.0.0D0) weidis=0.1D0
+ do i=nnt,nct
+ do j=i+2,nct
+ nhpb=nhpb+1
+ ihpb(nhpb)=i
+ jhpb(nhpb)=j
+ forcon(nhpb)=weidis
+ enddo
+ enddo
+ read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
+ write (iout,'(a,i5)') 'nexcl: ',nexcl
+ write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
+ return
+ end subroutine read_threadbase
+!-----------------------------------------------------------------------------
+#ifdef WHAM_RUN
+!el subroutine read_angles(kanal,iscor,energ,iprot,*)
+ subroutine read_angles(kanal,*)
+
+ use geometry_data
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'COMMON.INTERACT'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+ character(len=80) :: lineh
+ integer :: iscor,iprot,ic
+ real(kind=8) :: energ
+#else
+ subroutine read_angles(kanal,*)
+
+ use geometry_data
+ ! use energy
+ ! use control
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTROL'
+#endif
+! Read angles from input
+!
+!el local variables
+ integer :: i,kanal
+#ifdef WHAM_RUN
+ read(kanal,'(a80)',end=10,err=10) lineh
+ read(lineh(:5),*,err=8) ic
+ read(lineh(6:),*,err=8) energ
+ goto 9
+ 8 ic=1
+ print *,'error, assuming e=1d10',lineh
+ energ=1d10
+ nss=0
+ 9 continue
+ read(lineh(18:),*,end=10,err=10) nss
+ IF (NSS.LT.9) THEN
+ read (lineh(20:),*,end=10,err=10) &
+ (IHPB(I),JHPB(I),I=1,NSS),iscor
+ ELSE
+ read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
+ read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
+ I=9,NSS),iscor
+ ENDIF
+! print *,"energy",energ," iscor",iscor
+#endif
+ read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
+ read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
+ read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
+ read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
+
+ do i=1,nres
+! 9/7/01 avoid 180 deg valence angle
+ if (theta(i).gt.179.99d0) theta(i)=179.99d0
+!
+ theta(i)=deg2rad*theta(i)
+ phi(i)=deg2rad*phi(i)
+ alph(i)=deg2rad*alph(i)
+ omeg(i)=deg2rad*omeg(i)
+ enddo
+ return
+ 10 return 1
+ end subroutine read_angles
+!-----------------------------------------------------------------------------
+ subroutine reada(rekord,lancuch,wartosc,default)
+
+! implicit none
+ character*(*) :: rekord,lancuch
+ real(kind=8) :: wartosc,default
+ integer :: iread !,ilen
+!el external ilen
+ iread=index(rekord,lancuch)
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,err=10,end=10) wartosc
+ return
+ 10 wartosc=default
+ return
+ end subroutine reada
+!-----------------------------------------------------------------------------
+ subroutine readi(rekord,lancuch,wartosc,default)
+
+! implicit none
+ character*(*) :: rekord,lancuch
+ integer :: wartosc,default
+ integer :: iread !,ilen
+!el external ilen
+ iread=index(rekord,lancuch)
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,err=10,end=10) wartosc
+ return
+ 10 wartosc=default
+ return
+ end subroutine readi
+!-----------------------------------------------------------------------------
+ subroutine multreadi(rekord,lancuch,tablica,dim,default)
+
+! implicit none
+ integer :: dim,i
+ integer :: tablica(dim),default
+ character*(*) :: rekord,lancuch
+ character(len=80) :: aux
+ integer :: iread !,ilen
+!el external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end subroutine multreadi
+!-----------------------------------------------------------------------------
+ subroutine multreada(rekord,lancuch,tablica,dim,default)
+
+! implicit none
+ integer :: dim,i
+ real(kind=8) :: tablica(dim),default
+ character*(*) :: rekord,lancuch
+ character(len=80) :: aux
+ integer :: iread !,ilen
+!el external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end subroutine multreada
+!-----------------------------------------------------------------------------
+ subroutine card_concat(card,to_upper)
+
+! dla UNRESA to_upper jest zawsze .true.
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+ character(*) :: card
+ character(len=80) :: karta !,ucase
+ logical :: to_upper
+!el external ilen
+ read (inp,'(a)') karta
+ if (to_upper) karta=ucase(karta)
+ card=' '
+ do while (karta(80:80).eq.'&')
+ card=card(:ilen(card)+1)//karta(:79)
+ read (inp,'(a)') karta
+ if (to_upper) karta=ucase(karta)
+ enddo
+ card=card(:ilen(card)+1)//karta
+ return
+ end subroutine card_concat
+!-----------------------------------------------------------------------------
+ subroutine read_dist_constr
+ use MPI_data
+ ! use control
+ use geometry, only: dist
+ use geometry_data
+ use control_data
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.SBRIDGE'
+ integer,dimension(2,100) :: ifrag_,ipair_
+ real(kind=8),dimension(100) :: wfrag_,wpair_
+ character(len=640) :: controlcard
+
+!el local variables
+ integer :: i,k,j,ddjk,ii,jj,itemp
+ integer :: nfrag_,npair_,ndist_
+ real(kind=8) :: dist_cut
+
+! write (iout,*) "Calling read_dist_constr"
+! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+! call flush(iout)
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+! write (iout,*) "IFRAG"
+! do i=1,nfrag_
+! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+! enddo
+! write (iout,*) "IPAIR"
+! do i=1,npair_
+! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+! enddo
+ if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
+ if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
+ if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
+ if(.not.allocated(forcon)) allocate(forcon(maxdim))
+
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
+ ifrag_(2,i)=nstart_sup+nsup-1
+! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+! write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king) &
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king) &
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king) &
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+ nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ endif
+ enddo
+ call flush(iout)
+ return
+ end subroutine read_dist_constr
+!-----------------------------------------------------------------------------
+#ifdef WINIFL
+ subroutine flush(iu)
+ return
+ end subroutine flush
+#endif
+#ifdef AIX
+ subroutine flush(iu)
+ call flush_(iu)
+ return
+ end subroutine flush
+#endif
+!-----------------------------------------------------------------------------
+ subroutine copy_to_tmp(source)
+
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ character*(*) :: source
+ character(len=256) :: tmpfile
+! integer ilen
+!el external ilen
+ logical :: ex
+ tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
+ inquire(file=tmpfile,exist=ex)
+ if (ex) then
+ write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
+ " to temporary directory..."
+ write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
+ call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
+ endif
+ return
+ end subroutine copy_to_tmp
+!-----------------------------------------------------------------------------
+ subroutine move_from_tmp(source)
+
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ character*(*) :: source
+! integer ilen
+!el external ilen
+ write (*,*) "Moving ",source(:ilen(source)),&
+ " from temporary directory to working directory"
+ write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
+ call system("/bin/mv "//source(:ilen(source))//" "//curdir)
+ return
+ end subroutine move_from_tmp
+!-----------------------------------------------------------------------------
+! misc.f
+!-----------------------------------------------------------------------------
+! $Date: 1994/10/12 17:24:21 $
+! $Revision: 2.5 $
+
+ logical function find_arg(ipos,line,errflag)
+
+ integer, parameter :: maxlen = 80
+ character(len=80) :: line
+ character(len=1) :: empty=' ',equal='='
+ logical :: errflag
+ integer :: ipos
+! This function returns .TRUE., if an argument follows keyword keywd; if so
+! IPOS will point to the first non-blank character of the argument. Returns
+! .FALSE., if no argument follows the keyword; in this case IPOS points
+! to the first non-blank character of the next keyword.
+
+ do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
+ ipos=ipos+1
+ enddo
+ errflag=.false.
+ if (line(ipos:ipos).eq.equal) then
+ find_arg=.true.
+ ipos=ipos+1
+ do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
+ ipos=ipos+1
+ enddo
+ if (ipos.gt.maxlen) errflag=.true.
+ else
+ find_arg=.false.
+ endif
+
+ return
+ end function find_arg
+!-----------------------------------------------------------------------------
+ logical function find_group(iunit,jout,key1)
+
+ character*(*) :: key1
+ character(len=80) :: karta !,ucase
+ integer :: iunit,jout
+ integer :: ll !,ilen
+!EL external ilen
+!EL logical lcom
+ rewind (iunit)
+ karta=' '
+ ll=ilen(key1)
+ do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
+ read (iunit,'(a)',end=10) karta
+ enddo
+ write (jout,'(2a)') '> ',karta(1:78)
+ find_group=.true.
+ return
+ 10 find_group=.false.
+ return
+ end function find_group
+!-----------------------------------------------------------------------------
+ logical function iblnk(charc)
+ character(len=1) :: charc
+ integer :: n
+ n = ichar(charc)
+ iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
+ return
+ end function iblnk
+!-----------------------------------------------------------------------------
+ integer function ilen(string)
+ character*(*) :: string
+!EL logical :: iblnk
+
+ ilen = len(string)
+1 if ( ilen .gt. 0 ) then
+ if ( iblnk( string(ilen:ilen) ) ) then
+ ilen = ilen - 1
+ goto 1
+ endif
+ endif
+ return
+ end function ilen
+!-----------------------------------------------------------------------------
+ integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
+ integer :: nkey,i,ikey,narg
+ character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
+! character(len=16) :: ucase
+
+ do i=1,narg
+ if (ucase(keywd).eq.keywdset(i,ikey)) then
+! Match found
+ in_keywd_set=i
+ return
+ endif
+ enddo
+! No match to the allowed set of keywords if this point is reached.
+ in_keywd_set=0
+ return
+ end function in_keywd_set
+!-----------------------------------------------------------------------------
+ character function lcase(string)
+ integer :: i, k, idiff
+ character*(*) :: string
+ character(len=1) :: c
+ character(len=40) :: chtmp
+!
+ i = len(lcase)
+ k = len(string)
+ if (i .lt. k) then
+ k = i
+ if (string(k+1:) .ne. ' ') then
+ chtmp = string
+ endif
+ endif
+ idiff = ichar('a') - ichar('A')
+ lcase = string
+ do 99 i = 1, k
+ c = string(i:i)
+ if (lge(c,'A') .and. lle(c,'Z')) then
+ lcase(i:i) = char(ichar(c) + idiff)
+ endif
+ 99 continue
+ return
+ end function lcase
+!-----------------------------------------------------------------------------
+ logical function lcom(ipos,karta)
+ character(len=80) :: karta
+ character :: koment(2) = (/'!','#'/)
+ integer :: ipos,i
+
+ lcom=.false.
+ do i=1,2
+ if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
+ enddo
+ return
+ end function lcom
+!-----------------------------------------------------------------------------
+ logical function lower_case(ch)
+ character*(*) :: ch
+ lower_case=(ch.ge.'a' .and. ch.le.'z')
+ return
+ end function lower_case
+!-----------------------------------------------------------------------------
+ subroutine mykey(line,keywd,ipos,blankline,errflag)
+
+! This subroutine seeks a non-empty substring keywd in the string LINE.
+! The substring begins with the first character different from blank and
+! "=" encountered right to the pointer IPOS (inclusively) and terminates
+! at the character left to the first blank or "=". When the subroutine is
+! exited, the pointer IPOS is moved to the position of the terminator in LINE.
+! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
+! only separators or the maximum length of the data line (80) has been reached.
+! The logical variable ERRFLAG is set at .TRUE. if the string
+! consists only from a "=".
+ integer, parameter :: maxlen=80
+ character(len=1) :: empty=' ',equal='=',comma=','
+ character*(*) :: keywd
+ character(len=80) :: line
+ logical :: blankline,errflag !EL,lcom
+ integer :: ipos,istart,iend
+ errflag=.false.
+ do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
+ ipos=ipos+1
+ enddo
+ if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
+! At this point the rest of the input line turned out to contain only blanks
+! or to be commented out.
+ blankline=.true.
+ return
+ endif
+ blankline=.false.
+ istart=ipos
+! Checks whether the current char is a separator.
+ do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
+ .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
+ ipos=ipos+1
+ enddo
+ iend=ipos-1
+! Error flag set to .true., if the length of the keyword was found less than 1.
+ if (iend.lt.istart) then
+ errflag=.true.
+ return
+ endif
+ keywd=line(istart:iend)
+ return
+ end subroutine mykey
+!-----------------------------------------------------------------------------
+ subroutine numstr(inum,numm)
+ character(len=10) :: huj='0123456789'
+ character*(*) :: numm
+ integer :: inumm,inum,inum1,inum2
+ inumm=inum
+ inum1=inumm/10
+ inum2=inumm-10*inum1
+ inumm=inum1
+ numm(3:3)=huj(inum2+1:inum2+1)
+ inum1=inumm/10
+ inum2=inumm-10*inum1
+ inumm=inum1
+ numm(2:2)=huj(inum2+1:inum2+1)
+ inum1=inumm/10
+ inum2=inumm-10*inum1
+ inumm=inum1
+ numm(1:1)=huj(inum2+1:inum2+1)
+ return
+ end subroutine numstr
+!-----------------------------------------------------------------------------
+ function ucase(string)
+ integer :: i, k, idiff
+ character(*) :: string
+ character(len=len(string)) :: ucase
+ character(len=1) :: c
+ character(len=40) :: chtmp
+!
+ i = len(ucase)
+ k = len(string)
+ if (i .lt. k) then
+ k = i
+ if (string(k+1:) .ne. ' ') then
+ chtmp = string
+ endif
+ endif
+ idiff = ichar('a') - ichar('A')
+ ucase = string
+ do 99 i = 1, k
+ c = string(i:i)
+ if (lge(c,'a') .and. lle(c,'z')) then
+ ucase(i:i) = char(ichar(c) - idiff)
+ endif
+ 99 continue
+ return
+ end function ucase
+!-----------------------------------------------------------------------------
+! geomout.F
+!-----------------------------------------------------------------------------
+ subroutine pdbout(etot,tytul,iunit)
+
+ use geometry_data, only: c,nres
+ use energy_data
+ ! use control
+ use compare_data
+ use MD_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.HEADER'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.DISTFIT'
+! include 'COMMON.MD'
+!el character(len=50) :: tytul
+ character*(*) :: tytul
+ character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
+ integer,dimension(nres) :: ica !(maxres)
+
+!el local variables
+ integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
+ real(kind=8) :: etot
+ integer :: nres2
+ nres2=2*nres
+
+ if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
+
+ write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
+!model write (iunit,'(a5,i6)') 'MODEL',1
+ if (nhfrag.gt.0) then
+ do j=1,nhfrag
+ iti=itype(hfrag(1,j))
+ itj=itype(hfrag(2,j))
+ if (j.lt.10) then
+ write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
+ 'HELIX',j,'H',j,&
+ restyp(iti),hfrag(1,j)-1,&
+ restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+ else
+ write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
+ 'HELIX',j,'H',j,&
+ restyp(iti),hfrag(1,j)-1,&
+ restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+ endif
+ enddo
+ endif
+
+ if (nbfrag.gt.0) then
+
+ do j=1,nbfrag
+
+ iti=itype(bfrag(1,j))
+ itj=itype(bfrag(2,j)-1)
+
+ write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
+ 'SHEET',1,'B',j,2,&
+ restyp(iti),bfrag(1,j)-1,&
+ restyp(itj),bfrag(2,j)-2,0
+
+ if (bfrag(3,j).gt.bfrag(4,j)) then
+
+ itk=itype(bfrag(3,j))
+ itl=itype(bfrag(4,j)+1)
+
+ write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
+ 'SHEET',2,'B',j,2,&
+ restyp(itl),bfrag(4,j),&
+ restyp(itk),bfrag(3,j)-1,-1,&
+ "N",restyp(itk),bfrag(3,j)-1,&
+ "O",restyp(iti),bfrag(1,j)-1
+
+ else
+
+ itk=itype(bfrag(3,j))
+ itl=itype(bfrag(4,j)-1)
+
+
+ write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
+ 'SHEET',2,'B',j,2,&
+ restyp(itk),bfrag(3,j)-1,&
+ restyp(itl),bfrag(4,j)-2,1,&
+ "N",restyp(itk),bfrag(3,j)-1,&
+ "O",restyp(iti),bfrag(1,j)-1
+
+
+
+ endif
+
+ enddo
+ endif
+
+ if (nss.gt.0) then
+ do i=1,nss
+ if (dyn_ss) then
+ write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
+ 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
+ 'CYS',jdssb(i)-nnt+1
+ else
+ write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
+ 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
+ 'CYS',jhpb(i)-nnt+1-nres
+ endif
+ enddo
+ endif
+
+ iatom=0
+ ichain=1
+ ires=0
+ do i=nnt,nct
+ iti=itype(i)
+ if (iti.eq.ntyp1) then
+ ichain=ichain+1
+ ires=0
+ write (iunit,'(a)') 'TER'
+ else
+ ires=ires+1
+ iatom=iatom+1
+ ica(i)=iatom
+ write (iunit,10) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,i),j=1,3),vtot(i)
+ if (iti.ne.10) then
+ iatom=iatom+1
+ write (iunit,20) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,nres+i),j=1,3),&
+ vtot(i+nres)
+ endif
+ endif
+ enddo
+ write (iunit,'(a)') 'TER'
+ do i=nnt,nct-1
+ if (itype(i).eq.ntyp1) cycle
+ if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+ write (iunit,30) ica(i),ica(i+1)
+ else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ write (iunit,30) ica(i),ica(i+1),ica(i)+1
+ else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+ write (iunit,30) ica(i),ica(i)+1
+ endif
+ enddo
+ if (itype(nct).ne.10) then
+ write (iunit,30) ica(nct),ica(nct)+1
+ endif
+ do i=1,nss
+ if (dyn_ss) then
+ write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+ else
+ write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+ endif
+ enddo
+ write (iunit,'(a6)') 'ENDMDL'
+ 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+ 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+ 30 FORMAT ('CONECT',8I5)
+ return
+ end subroutine pdbout
+!-----------------------------------------------------------------------------
+ subroutine MOL2out(etot,tytul)
+! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
+! format.
+ use geometry_data, only: c
+ use energy_data
+ ! use control
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.HEADER'
+! include 'COMMON.SBRIDGE'
+ character(len=32) :: tytul,fd
+ character(len=3) :: zahl
+ character(len=6) :: res_num,pom !,ucase
+
+!el local variables
+ integer :: i,j
+ real(kind=8) :: etot
+
+#ifdef AIX
+ call fdate_(fd)
+#elif (defined CRAY)
+ call date(fd)
+#else
+ call fdate(fd)
+#endif
+ write (imol2,'(a)') '#'
+ write (imol2,'(a)') &
+ '# Creating user name: unres'
+ write (imol2,'(2a)') '# Creation time: ',&
+ fd
+ write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
+ write (imol2,'(a)') tytul
+ write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
+ write (imol2,'(a)') 'SMALL'
+ write (imol2,'(a)') 'USER_CHARGES'
+ write (imol2,'(a)') '\@<TRIPOS>ATOM'
+ do i=nnt,nct
+ write (zahl,'(i3)') i
+ pom=ucase(restyp(itype(i)))
+ res_num = pom(:3)//zahl(2:)
+ write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
+ enddo
+ write (imol2,'(a)') '\@<TRIPOS>BOND'
+ do i=nnt,nct-1
+ write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
+ enddo
+ do i=1,nss
+ write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
+ enddo
+ write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
+ do i=nnt,nct
+ write (zahl,'(i3)') i
+ pom = ucase(restyp(itype(i)))
+ res_num = pom(:3)//zahl(2:)
+ write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
+ enddo
+ 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
+ 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
+ return
+ end subroutine MOL2out
+!-----------------------------------------------------------------------------
+ subroutine intout
+
+ use geometry_data
+ use energy_data, only: itype
+ ! use control
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.TORSION'
+!el local variables
+ integer :: i,iti
+
+ write (iout,'(/a)') 'Geometry of the virtual chain.'
+ write (iout,'(7a)') ' Res ',' d',' Theta',&
+ ' Phi',' Dsc',' Alpha',' Omega'
+ do i=1,nres
+ iti=itype(i)
+ write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
+ rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
+ rad2deg*omeg(i)
+ enddo
+ return
+ end subroutine intout
+!-----------------------------------------------------------------------------
+#ifdef CLUSTER
+ subroutine briefout(it,ener,free)!,plik)
+#else
+ subroutine briefout(it,ener)
+#endif
+ use geometry_data
+ use energy_data
+ ! use control
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.SBRIDGE'
+! print '(a,i5)',intname,igeom
+!el local variables
+ integer :: i,it
+ real(kind=8) :: ener,free
+! character(len=80) :: plik
+! integer :: iii
+
+#if defined(AIX) || defined(PGI)
+ open (igeom,file=intname,position='append')
+#else
+ open (igeom,file=intname,access='append')
+#endif
+#ifdef WHAM_RUN
+! iii=igeom
+ igeom=iout
+#endif
+ IF (NSS.LE.9) THEN
+#ifdef CLUSTER
+ WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+ ELSE
+ WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
+#else
+ WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+ ELSE
+ WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
+#endif
+ WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
+ ENDIF
+! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+ WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+ WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
+! if (nvar.gt.nphi+ntheta) then
+ write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
+ write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
+! endif
+ close(igeom)
+ 180 format (I5,F12.3,I2,9(1X,2I3))
+ 190 format (3X,11(1X,2I3))
+ 200 format (8F10.4)
+ return
+ end subroutine briefout
+!-----------------------------------------------------------------------------
+#ifdef WINIFL
+ subroutine fdate(fd)
+ character(len=32) :: fd
+ write(fd,'(32x)')
+ return
+ end subroutine fdate
+#endif
+!-----------------------------------------------------------------------------
+#ifdef WHAM_RUN
+ real(kind=8) function gyrate(jcon)
+#else
+ real(kind=8) function gyrate()
+#endif
+
+ use geometry_data, only: c
+ ! use geometry
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CHAIN'
+ real(kind=8) :: rg
+ real(kind=8),dimension(3) :: cen
+!el local variables
+ integer :: i,j,jcon
+
+ do j=1,3
+ cen(j)=0.0d0
+ enddo
+
+ do i=nnt,nct
+ do j=1,3
+ cen(j)=cen(j)+c(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cen(j)=cen(j)/dble(nct-nnt+1)
+ enddo
+ rg = 0.0d0
+ do i = nnt, nct
+ do j=1,3
+ rg = rg + (c(j,i)-cen(j))**2
+ enddo
+ end do
+#ifdef WHAM_RUN
+ gyrate = dsqrt(rg/dble(nct-nnt+1))
+#else
+ gyrate = sqrt(rg/dble(nct-nnt+1))
+#endif
+ return
+ end function gyrate
+#ifdef WHAM_RUN
+!-----------------------------------------------------------------------------
+! readrtns.F WHAM
+ subroutine reads(rekord,lancuch,wartosc,default)
+! implicit none
+ character*(*) :: rekord,lancuch,wartosc,default
+ character(len=80) :: aux
+ integer :: lenlan,lenrec,iread,ireade
+!el external ilen
+!el logical iblnk
+!el external iblnk
+ lenlan=ilen(lancuch)
+ lenrec=ilen(rekord)
+ iread=index(rekord,lancuch(:lenlan)//"=")
+! print *,"rekord",rekord," lancuch",lancuch
+! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+lenlan+1
+! print *,"iread",iread
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
+ iread=iread+1
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ enddo
+! print *,"iread",iread
+ if (iread.gt.lenrec) then
+ wartosc=default
+ return
+ endif
+ ireade=iread+1
+! print *,"ireade",ireade
+ do while (ireade.lt.lenrec .and. &
+ .not.iblnk(rekord(ireade:ireade)))
+ ireade=ireade+1
+ enddo
+ wartosc=rekord(iread:ireade)
+ return
+ end subroutine reads
+#endif
+!-----------------------------------------------------------------------------
+! permut.F
+!-----------------------------------------------------------------------------
+ subroutine permut(isym)
+
+ use geometry_data, only: tabperm
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+
+ integer :: n,isym
+! logical nextp
+!el external nextp
+ integer,dimension(isym) :: a
+! parameter(n=symetr)
+!el local variables
+ integer :: kkk,i
+
+ n=isym
+ if (n.eq.1) then
+ tabperm(1,1)=1
+ return
+ endif
+ kkk=0
+ do i=1,n
+ a(i)=i
+ enddo
+ 10 print *,(a(i),i=1,n)
+ kkk=kkk+1
+ do i=1,n
+ tabperm(kkk,i)=a(i)
+! write (iout,*) "tututu", kkk
+ enddo
+ if(nextp(n,a)) go to 10
+ return
+ end subroutine permut
+!-----------------------------------------------------------------------------
+ logical function nextp(n,a)
+
+ integer :: n,i,j,k,t
+! logical :: nextp
+ integer,dimension(n) :: a
+ i=n-1
+ 10 if(a(i).lt.a(i+1)) go to 20
+ i=i-1
+ if(i.eq.0) go to 20
+ go to 10
+ 20 j=i+1
+ k=n
+ 30 t=a(j)
+ a(j)=a(k)
+ a(k)=t
+ j=j+1
+ k=k-1
+ if(j.lt.k) go to 30
+ j=i
+ if(j.ne.0) go to 40
+ nextp=.false.
+ return
+ 40 j=j+1
+ if(a(j).lt.a(i)) go to 40
+ t=a(i)
+ a(i)=a(j)
+ a(j)=t
+ nextp=.true.
+ return
+ end function nextp
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module io_base