text:
c This module/program is part of DaliLite (c) L. Holm 1999
c
c
c----------------------------------------------------------------------
c
c	stack module
c		getnextbest: return best unique alignment or exception flag
c		clearstack: after overflow, make space by removing baddies
c		swapstack: chunk1 <-> chunk2
c		putchunk, putstack: pack, save chunk
c		getchunk, getstack: read, unpack chunk
c		declump: blank out alignment trace from stack
c		split: divide search space
c
c	variables
c		chunk: 1+[sum(is=1,ns)mi(seglist(is))+30]/31; est+candidates
c		link(): linked list of score estimates, points to stack
c			link(1,*) lower, link(2,*) higher score
c		bestest: points to beginning of linked list
c		stack: in chunks
c		top: starts from 1, points to estimate of chunk
c		ex(): energy table
c		start(iseg,jseg): pointer to ex(), index<0 means score is zero
c		ns: number of segments
c		mi: maximal number of residues in search space, including nul
c		ni,ci: work search space [packed]
c		ci0: initial search space [packed]
c
c----------------------------------------------------------------------
c
	subroutine split(est,ns,seglist,ex,top,stack,scorecutoff,ni,ci,mi,
     $		chunk,bestest,link,start,link_pointer,forward_pointer,
     $		backward_pointer,firstest,trans,lseqtl)
	implicit none
	include 'parsizes.for'
	integer re_estimate,start(maxseg,maxseg),trans(maxres0,maxseg)
	integer ns,seglist(maxseg),mi(maxseg),chunk,bestest,
     $ link(2,maxstack2)
	integer est,ex(exdim),scorecutoff,top,stack(maxstack)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
	logical lseqtl
c
c	split xseg according to ex(xseg,yseg,...)
c		find emax(xres)= max of ex(xseg,yseg,xres,...)
c		erange=smallest to largest emax(...)
c		lcand1(xres)=lcand(xseg,xres).and.emax(xres).le.midpoint
c		lcand2(xres)=lcand(xseg,xres).and.emax(xres).gt.midpoint
c		special case: smallest.eq.largest => assign every second to 1,2
c
	integer mini,maxi,xres,e,emax(maxres0),midpoint,
     $ ess(maxseg,maxseg)
	integer f,n1,ni(maxseg),ci(maxres0,maxseg),xr,yr
	integer iwhere,jw,ijstart
c
	integer emaxim,is,js,iseg,jseg,xseg,yseg,seed
	integer ni1,ni2,ci1(maxres0),ci2(maxres0),ci3(maxres0),
     $ ci4(maxres0)
	logical lclosed(maxseg)
	integer bestxres,ni3,ni4,ni5,ci5(maxres0),pr,ni6,ci6(maxres0)
c!!
	seed=1234567
c
	do is=1,ns
		iseg=seglist(is)
		lclosed(iseg)=(ni(iseg).le.1)
	end do
c
!	write(*,*) 'split',est,ns,top,scorecutoff
	call get_ess(ess,est,ns,seglist,ni,ci,ex,start,mi,trans,lseqtl)
	if(est.le.scorecutoff) return
c
	xseg=0
	yseg=0
	emaxim=-infinit
	do is=1,ns
	  iseg=seglist(is)
	  call getemax(ess(iseg,iseg),emaxim,lclosed,iseg,iseg,
     $ xseg,yseg,seed)
	  do js=1,is-1
	    jseg=seglist(js)
	    call getemax(ess(iseg,jseg),emaxim,lclosed,iseg,jseg,xseg,
     $ yseg,seed)
	  end do
	end do
	if((xseg.eq.0).and.(yseg.eq.0)) return
c
	if(start(xseg,yseg).lt.0) then
		maxi=0
		mini=0
		n1=0
	else
		maxi=-infinit
		mini=+infinit
		if(xseg.eq.yseg) then
			iwhere=start(xseg,xseg)
			do xr=1,ni(xseg)
				xres=ci(xr,xseg)
				emax(xres)=ex(iwhere+xres)
			end do
		else if(xseg.lt.yseg) then
			ijstart=start(xseg,yseg)
			jw=mi(yseg)
			do xr=1,ni(xseg)
				e=-infinit
				xres=ci(xr,xseg)
				iwhere=ijstart+(xres-1)*jw
				do yr=1,ni(yseg)
					e=max(e,ex(iwhere+ci(yr,yseg)))
				end do
				emax(xres)=e
			end do
		else if(xseg.gt.yseg) then
			do xr=1,ni(xseg)
				emax(ci(xr,xseg))=-infinit
			end do
			ijstart=start(yseg,xseg)
			jw=mi(xseg)
			do yr=1,ni(yseg)
				iwhere=ijstart+(ci(yr,yseg)-1)*jw
				do xr=1,ni(xseg)
					xres=ci(xr,xseg)
					e=emax(xres)
					emax(xres)=max(e,ex(iwhere+xres))
				end do
			end do
		end if
		n1=ni(xseg)
	end if
	bestxres=ci(1,xseg)
	do xr=1,n1
		xres=ci(xr,xseg)
		e=emax(xres)
		if(xr.gt.2) then
			if(e.gt.emax(bestxres)) bestxres=xres
		end if
		maxi=max(maxi,e)
		mini=min(mini,e)
	end do
c
c	split 90 % - 10 %
c
	midpoint=mini+(maxi-mini)*9/10
	ni1=0
	ni2=0
	if(mini.eq.maxi) then
		do xr=1,ni(xseg),2
			ni1=ni1+1
			ci1(ni1)=ci(xr,xseg)
		end do
		do xr=2,ni(xseg),2
			ni2=ni2+1
			ci2(ni2)=ci(xr,xseg)
		end do
	else
		do xr=1,ni(xseg)
			xres=ci(xr,xseg)
			if(emax(xres).le.midpoint) then
				ni2=ni2+1
				ci2(ni2)=ci(xr,xseg)
			else
				ni1=ni1+1
				ci1(ni1)=ci(xr,xseg)
			end if
		end do
	end if
!	if(ni1.le.0.or.ni2.le.0.or.ni1+ni2.ne.ni(xseg)) then
!	  write(*,*) 'mid:',midpoint,mini,maxi,xseg,yseg,ni(xseg),ni1,ni2
!	  stop
!	end if
c
c	stack better half
c
	if(lseqtl) then 	! split better half left&nul-top&right ci3,c4
		ni3=0
		ni4=0
		ni5=0
		ni6=0
		do xres=1,ni1
			pr=trans(ci1(xres),xseg)
			if(pr.eq.nul) then
				ni6=ni6+1
				ci6(ni6)=ci1(xres)
			else if(ci1(xres).lt.bestxres) then
				ni3=ni3+1
				ci3(ni3)=ci1(xres)
			else
				ni4=ni4+1
				ci4(ni4)=ci1(xres)
			end if
		end do
		if(ni6.gt.0) call copyandput(ni6,ci6,ni,ci,est,ns,
     $			seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $          	scorecutoff,top,stack,chunk,bestest,link,
     $			link_pointer,forward_pointer,backward_pointer,
     $			firstest)
		if(lskipflag) return
		if(ni5.gt.0) call copyandput(ni5,ci5,ni,ci,est,ns,
     $			seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $          	scorecutoff,top,stack,chunk,bestest,link,
     $			link_pointer,forward_pointer,backward_pointer,
     $			firstest)
		if(lskipflag) return
		if(ni3.gt.0) call copyandput(ni3,ci3,ni,ci,est,ns,
     $			seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $          	scorecutoff,top,stack,chunk,bestest,link,
     $			link_pointer,forward_pointer,backward_pointer,
     $			firstest)
		if(lskipflag) return
		if(ni4.gt.0) call copyandput(ni4,ci4,ni,ci,est,ns,
     $			seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $          	scorecutoff,top,stack,chunk,bestest,link,
     $			link_pointer,forward_pointer,backward_pointer,
     $			firstest)
		if(lskipflag) return
	else
		call copyandput(ni1,ci1,ni,ci,est,ns,
     $			seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $          	scorecutoff,top,stack,chunk,bestest,link,
     $			link_pointer,forward_pointer,backward_pointer,
     $			firstest)
 		if(lskipflag) return
	end if
c
c       stack lower half
c
	call copyandput(ni2,ci2,ni,ci,est,ns,
     $		seglist,ess,xseg,ex,start,mi,trans,lseqtl,
     $         	scorecutoff,top,stack,chunk,bestest,link,
     $		link_pointer,forward_pointer,backward_pointer,
     $		firstest)
	if(lskipflag) return
!       write(*,*) 'top',top

	return
	end
c
c----------------------------------------------------------------------
c
        subroutine copyandput(ni1,ci1,ni,ci,est,ns,seglist,ess,xseg,ex,
     $          start,mi,trans,lseqtl,scorecutoff,top,stack,chunk,
     $		bestest,link,link_pointer,forward_pointer,
     $          backward_pointer,firstest)
        implicit none
        include 'parsizes.for'
        integer re_estimate,start(maxseg,maxseg),trans(maxres0,maxseg)
        integer ns,seglist(maxseg),mi(maxseg),chunk,bestest
	integer link(2,maxstack2)
        integer est,ex(exdim),scorecutoff,top,stack(maxstack)
        integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $          backward_pointer(0:boxdim)
        real firstest
        logical lseqtl
        integer xres,ess(maxseg,maxseg),f,ni(maxseg),ci(maxres0,maxseg)
        integer xseg,ni1,ci1(maxres0)
c
        if(ni1.eq.0) return
        ni(xseg)=ni1
        do xres=1,ni1
                ci(xres,xseg)=ci1(xres)
        end do
!        call excludeimpossible(ni,ci,ns,seglist,trans,lseqtl)
        f=re_estimate(est,ns,seglist,ess,xseg,ex,ni,ci,start,
     $          mi,trans,lseqtl)
        if(f.gt.scorecutoff) then
                top=top+chunk
                call putrequest(f,ni,ci,ns,seglist,mi,top,stack,chunk,
     $                  bestest,link,link_pointer,
     $                  forward_pointer,backward_pointer,firstest)
		if(lskipflag) return
        end if

        return
        end
c
c----------------------------------------------------------------------
c
	subroutine excludeimpossible(ni,ci,ns,seglist,trans,lseqtl)
	implicit none
	include 'parsizes.for'
	integer ni(maxseg),ci(maxres0,maxseg),trans(maxres0,maxseg)
	integer ns,seglist(maxseg)
	logical lseqtl
c
	integer is,iseg,leftmost,ritemost,ir,irx,i,nx
	logical lhasnul(maxseg)
c
	if(.not.lseqtl) return
	nx=0
	do is=1,ns
		iseg=seglist(is)
		lhasnul(iseg)=(trans(ci(ni(iseg),iseg),iseg).eq.nul)
!		type *,'before',iseg,(trans(ci(ir,iseg),iseg),ir=1,ni(iseg))
	end do
	leftmost=-infinit
	do is=1,ns
	  iseg=seglist(is)
	  if(.not.lhasnul(iseg)) then
		irx=0
		do ir=1,ni(iseg)
			i=trans(ci(ir,iseg),iseg)
!			if(i.le.leftmost.and.i.ne.nul) irx=ir
			if(i.lt.leftmost.and.i.ne.nul) irx=ir
		end do
		nx=nx+irx
		! exclude 1...irx
		if(irx.gt.0) then
			i=0
			do ir=irx+1,ni(iseg)
				i=i+1
				ci(i,iseg)=ci(ir,iseg)
			end do
			ni(iseg)=i
		end if
		if(ni(iseg).gt.0) then
			if(.not.lhasnul(iseg)) then
				leftmost=trans(ci(1,iseg),iseg)
			else if (ni(iseg).gt.1) then
				leftmost=trans(ci(1,iseg),iseg)
			end if
		else
			leftmost=+infinit
			goto 99
		end if
		i=ci(1,iseg)
!		type *,'leftmost:',lhasnul(iseg),iseg,leftmost,i,
!     $			trans(i,iseg),irx
	  end if
	end do
	if(leftmost.eq.infinit) goto 99
10	ritemost=+infinit
	do is=ns,1,-1
		iseg=seglist(is)
		if(.not.lhasnul(iseg)) then
			irx=0
			do ir=ni(iseg),1,-1
				i=trans(ci(ir,iseg),iseg)
!				if(i.ge.ritemost.and.i.ne.nul) irx=ir
				if(i.gt.ritemost.and.i.ne.nul) irx=ir
			end do
			if(irx.gt.0) nx=nx+ni(iseg)-irx
			! exclude irx...ni(iseg) except nul
			if(irx.gt.0) then
				if(lhasnul(iseg)) then
					ci(irx,iseg)=ci(ni(iseg),iseg)
					ni(iseg)=irx
			  	else
					ni(iseg)=irx-1
				end if
			end if
			if(ni(iseg).gt.0) then
				if(.not.lhasnul(iseg)) then
				  ritemost=trans(ci(ni(iseg),iseg),iseg)
				else if(ni(iseg).gt.1) then
				  ritemost=trans(ci(ni(iseg)-1,iseg),iseg)
				end if
			else
				ritemost=-infinit
				goto 99
			end if
		end if
		i=ci(ni(iseg),iseg)
!		type *,'ritemost:',iseg,ritemost,i,trans(i,iseg),irx
	end do

99	continue
!	if(nx.gt.0) then
!	  do is=1,ns
!		iseg=seglist(is)
!		type *,'after',iseg,(trans(ci(ir,iseg),iseg),ir=1,ni(iseg))
!	  end do
!	  type *,'nx',nx
!	end if

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine getnextbest(scorecutoff,stack,top,chunk,ns,seglist,
     $		mi,bestest,ali,flag,link,ex,start,est,link_pointer,
     $		forward_pointer,backward_pointer,firstest,trans,lseqtl)
c
c	flags:	1 = empty stack
c		2 = top score below cutoff
c		4 = no candidates
c		5 = unique alignment, returned in ali()
c		3 = overflow
c
	implicit none
	include 'parsizes.for'
	integer scorecutoff,flag,stack(maxstack),top,chunk,ns,mi(maxseg)
	integer link(2,maxstack2),ali(maxseg),bestest,
     $ seglist(maxseg),ex(exdim)
	integer start(maxseg,maxseg),trans(maxres0,maxseg)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
	logical lseqtl
c
	integer is,est,nmin,nmax,ni(maxseg),ci(maxres0,maxseg),ir,iseg,i,m
c	write(*,*) 'getnextbest',scorecutoff,top,chunk,ns,bestest,firstest
	est=-infinit
	flag=0
c	check that bestest is!
c	call checkstack(stack,bestest,top,chunk,'getnext BUG ')
	do while(flag.eq.0)
c
c		check exceptions
c
		if(top.le.0) then
			flag=1
			est=-infinit
			return
		end if
		if(bestest.le.0) then
	write(*,*) scorecutoff,(stack(i),i=1,10),top,chunk,bestest
	write(*,*) (link(1,i),link(2,i),i=1,top)
	write(*,*) 'FATAL ERROR -- bestest out of range'
			lskipflag=.true.
			return
		end if
		if(stack(bestest).lt.scorecutoff) flag=2
		if(top+2*chunk.gt.maxstack) flag=3
c!!
		if(top.gt.maxstack-100) flag=3
		if(flag.ne.0) return
c
c		continue splitting / export unique alignment
c
		call toprequest(bestest,top,stack,chunk,link,mi,ns,seglist,
     $			ni,ci,est,link_pointer,forward_pointer,
     $			backward_pointer,firstest)
		if(lskipflag) return
c		write(*,*) 'top search space'
c		do is=1,ns
c			iseg=seglist(is)
c			write(*,*) is,iseg,mi(iseg),ni(iseg),
c     $				(ci(ir,iseg),ir=1,ni(iseg))
c		end do
		call getnminnmax(ni,ns,seglist,nmin,nmax)
		if(nmin.eq.0) flag=4
		if((nmin.eq.1).and.(nmax.eq.1)) flag=5
		if(flag.eq.0) call split(est,ns,seglist,ex,top,stack,
     $			scorecutoff,ni,ci,mi,chunk,bestest,link,start,
     $			link_pointer,forward_pointer,backward_pointer,firstest,
     $			trans,lseqtl)
		if(lskipflag) return
		if(flag.eq.5) then
			do is=1,ns
				ali(seglist(is))=ci(1,seglist(is))
			end do
		end if
c		write(*,*) 'flag is',flag
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine lpack(lwork,i,l)
	implicit none
	include 'parsizes.for'
	logical lwork(maxres0*maxseg),l(31)
	integer i,j

	do j=0,30
		l(j+1)=lwork(i+j)
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine putrequest(est,ni,ci,ns,seglist,mi,to,stack,
     $ chunk,bestest,link,link_pointer,forward_pointer,
     $ backward_pointer,firstest)
	implicit none
	include 'parsizes.for'
	integer ni(maxseg),ci(maxres0,maxseg),ns,mi(maxseg)
	integer seglist(maxseg)
	integer to,stack(maxstack),chunk,est,bestest,link(2,maxstack2)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer i,j,is,ir,irray((1+(maxres0+30)/31)*maxseg),s,iseg,n
	logical lwork(maxres0*maxseg),l(31),lw(maxres0)
c
c	write(*,*) 'putrequest',est,ns,to,chunk,bestest,firstest
c	call checkstack(stack,bestest,to,chunk,'put BUG 1')
	n=0
	do is=1,ns
		iseg=seglist(is)
		do ir=1,mi(iseg)
			lw(ir)=.false.
		end do
		do ir=1,ni(iseg)
			lw(ci(ir,iseg))=.true.
		end do
		do ir=1,mi(iseg)
			n=n+1
			lwork(n)=lw(ir)
		end do
	end do
	do i=n+1,min(n+30,maxres0*maxseg)
		lwork(i)=.false.
	end do
	j=1
	irray(1)=est
	do i=1,n,31
		call lpack(lwork,i,l)
		call bitpack(s,l,31)
		j=j+1
		irray(j)=s
c		write(*,*) 'l',i,j,s
	end do
c	write(*,*) 'put irray',(irray(i),i=1,j)
	if(to+chunk.gt.maxstack) then
	write(*,*)  'FATAL ERROR in putrequest'
		lskipflag=.true.
		return
	end if
	call putchunk(to,irray,stack,chunk)
	if(lskipflag) return
	call createlink(to,bestest,link,stack,link_pointer,
     $		forward_pointer,backward_pointer,firstest,to,chunk)
	if(lskipflag) return
c	call checkstack(stack,bestest,to,chunk,'put BUG 2')

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine declump(scorecutoff,stack,top,chunk,link,bestest,ns,
     $		seglist,mi,ali,ex,start,singletcutoff,link_pointer,
     $		forward_pointer,backward_pointer,firstest,trans,lseqtl)
	implicit none
	include 'parsizes.for'
	integer stack(maxstack),top,chunk,link(2,maxstack2),bestest
	integer scorecutoff,ns,mi(maxseg),ali(maxseg),seglist(maxseg)
	integer ex(exdim),start(maxseg,maxseg),singletcutoff
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim),trans(maxres0,maxseg)
	real firstest
	logical lseqtl
c
	integer i,irray((1+(maxres0+30)/31)*maxseg),est,ci(maxres0,maxseg)
	integer ni(maxseg),oldtop,is,ir,jr,ess(maxseg,maxseg),
     $ ires,iseg,nold
	logical lchange
	integer m
c
c
c	write(*,500) scorecutoff,top,chunk,bestest,ns,(ali(seglist(i)),i=1,ns)
	oldtop=top
	top=1-chunk
	bestest=0
	call init_link_pointer(link_pointer,forward_pointer,
     $ backward_pointer)
	do i=1,oldtop,chunk
		call getchunk(i,irray,stack,chunk)
		call translatechunk(irray,mi,ns,seglist,ni,ci,est,chunk)
c
c		only declump chunks below singletcutoff !
c
c????
		if(est.gt.singletcutoff) goto 20
c
c		blank out ali()
c
		lchange=.false.
		do is=1,ns
			iseg=seglist(is)
			nold=ni(iseg)
			m=ali(iseg)
			do ir=1,nold
			  ires=ci(ir,iseg)
			  if(ires.eq.m) then
			    if(trans(ires,iseg).ne.nul) then
c				write(*,*) 'remove',ires,iseg
				do jr=ir+1,ni(iseg)
					ci(jr-1,iseg)=ci(jr,iseg)
				end do
				ni(iseg)=ni(iseg)-1
c				no candidates left !
				if(ni(iseg).eq.0) then
					est=-infinit
					goto 29
				end if
				lchange=.true.
				goto 19
			    end if
			  end if
			end do
19			continue
		end do
c
c		re_estimate score if changes made
c
		if(lchange) call get_ess(ess,est,ns,seglist,ni,ci,ex,start,mi,
     $			trans,lseqtl)
20		if(est.gt.scorecutoff) then
			top=top+chunk
			call putrequest(est,ni,ci,ns,seglist,mi,
     $ top,stack,chunk,bestest,link,link_pointer,
     $ forward_pointer,backward_pointer,firstest)
			if(lskipflag) return
		else
	end if
29	end do
500	format('declump',4i12,85i4)
c   500 format('declump',4i12,<maxseg+5>i4)
	return
	end
c
c----------------------------------------------------------------------
c
	subroutine clearstack(cutoff,stack,top,chunk,link,bestest,
     $ link_pointer,forward_pointer,backward_pointer,firstest)
c
c	remove all chunks that have a score below cutoff
c
	implicit none
	include 'parsizes.for'
	integer cutoff,stack(maxstack),top,chunk,link(2,maxstack2),bestest
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer i,dn,oldtop,x,getk,e
c
c	compress chunks
c
c	write(*,*) 'clearstack in ',cutoff,top,chunk,bestest
	oldtop=top
	top=1-chunk
	do i=1,oldtop,chunk
		if((stack(i).gt.cutoff).and.(top.lt.i)) then
			top=top+chunk
			call swapstack(i,top,stack,chunk,link,bestest,
     $			 link_pointer,forward_pointer,backward_pointer,firstest)
			if(lskipflag) return
		end if
	end do
	call init_link_pointer(link_pointer,forward_pointer,
     $ backward_pointer)
	do i=1,top,chunk
		e=stack(i)
		x=getk(e,i,stack,link_pointer,forward_pointer,
     $			backward_pointer,firstest)
		if(lskipflag) return
	end do
c
c	curtail linked list
c
	do i=1,top,chunk
		dn=link(1,i)
		if((dn.gt.0).and.(stack(dn).le.cutoff)) then
			link(1,i)=0
c			write(*,*) 'clearstack out',cutoff,top,chunk,bestest
			return
		end if
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine init_link_pointer(link_pointer,forward_pointer,
     $		backward_pointer)
	implicit none
	include 'parsizes.for'
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
c
	integer i
c
	do i=0,boxdim
		link_pointer(i)=0
		forward_pointer(i)=boxdim
		backward_pointer(i)=0
	end do
c
	return
	end
c
c----------------------------------------------------------------------
c
	subroutine createlink_ok(i,bestest,link,stack,link_pointer,
     $		forward_pointer,backward_pointer,firstest,top,chunk)
c
c	link in new chunk 'i'
c
	implicit none
	include 'parsizes.for'
	integer i,bestest,link(2,maxstack2),stack(maxstack),top,chunk
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer j,k,est
c
	if(i.gt.maxstack2) then
	write(*,*)  'FATAL ERROR link overflow - increase maxstack2 '
		lskipflag=.true.
		return
	end if
	link(1,i)=0
	link(2,i)=0
	est=stack(i)
	j=bestest
	k=0
	do while(j.gt.0)
		if(stack(j).gt.est) k=j
		j=link(1,j)
	end do
	if(k.eq.0) then
		if(bestest.ne.0) then
			link(1,i)=bestest
			link(2,bestest)=i
		end if
		bestest=i
	else
		j=link(1,k)
		link(1,i)=j
		if(j.ne.0) link(2,j)=i
		link(1,k)=i
		link(2,i)=k
	end if

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine createlink(i,bestest,link,stack,link_pointer,
     $		forward_pointer,backward_pointer,firstest,top,chunk)
c
c	link in new chunk 'i'
c
	implicit none
	include 'parsizes.for'
	integer i,bestest,link(2,maxstack2),stack(maxstack)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim),top,chunk
	real firstest
c
	integer j,k,est,getk,x
c
c	write(*,*) 'createlink',i,bestest,firstest
	if(i.gt.maxstack2) then
	write(*,*)  'FATAL ERROR link overflow - increase maxstack2 '
		lskipflag=.true.
		return
	end if
	link(1,i)=0
	link(2,i)=0
	est=stack(i)
!	j=bestest
!	k=0
!	do while(j.gt.0)
!		if(stack(j).gt.est) k=j
!		j=link(1,j)
!	end do
!	x=k
	j=getk(est,i,stack,link_pointer,forward_pointer,
     $		backward_pointer,firstest)
	if(lskipflag) return
c	write(*,*) 'create',bestest,stack(bestest),j,stack(j),i,est
	k=0
	if(j.eq.0) j=bestest
	do while(j.gt.0)
		if(stack(j).gt.est) k=j
		j=link(1,j)
	end do
!	if(k.ne.x) write(*,*) 'BUG',x,stack(max(1,x)),bestest,
!     $		stack(max(1,bestest)),i,est,k,stack(max(1,k)),firstest
!	k=x
	if(k.eq.0) then
		if(bestest.ne.0) then
			link(1,i)=bestest
			link(2,bestest)=i
		end if
		bestest=i
c		call checkstack(stack,bestest,top,chunk,'create BUG 1')
	else
		j=link(1,k)
		link(1,i)=j
		if(j.ne.0) link(2,j)=i
		link(1,k)=i
		link(2,i)=k
c		call checkstack(stack,bestest,top,chunk,'create BUG 2')
	end if

	return
	end
c
c----------------------------------------------------------------------
c
	function getk(est,i,stack,link_pointer,forward_pointer,
     $		backward_pointer,firstest)
	implicit none
	include 'parsizes.for'
	integer getk,est,i,stack(maxstack)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer box,k,b,funbox
c
c	returns stack-pointer to a next larger estimate, 0 if new best
c
c	box=funbox(est,firstest)
	box=max(1,int(float(est)/firstest*float(boxdim-1)))
	if(box.gt.boxdim) then
	write(*,*)  'FATAL ERROR -- boxfun overflow'
		lskipflag=.true.
		return
	end if
	if(link_pointer(box).eq.0) then
		do b=box+1,forward_pointer(box)
			backward_pointer(b)=box
		end do
		do b=backward_pointer(box),box-1
			forward_pointer(b)=box
		end do
		link_pointer(box)=i
		box=forward_pointer(box)
	else
		if(est.ge.stack(link_pointer(box))) then
			link_pointer(box)=i
			box=forward_pointer(box)
		end if
	end if
	if(box.eq.boxdim) then
		k=0
	else
		k=link_pointer(box)
	end if

	getk=k
c	write(*,*) 'getk',est,i,k,box,(link_pointer(b),b=0,boxdim)

	return
	end
c
c----------------------------------------------------------------------
c
	function funbox(est,firstest)
	implicit none
	include 'parsizes.for'
	integer est,funbox
	real firstest,r
c
	r=float(est)/firstest
	funbox=max(1,int(r*r*(boxdim-1)))
c
	return
	end
c
c----------------------------------------------------------------------
c
	subroutine removebox(est,i,stack,link,link_pointer,
     $ forward_pointer,backward_pointer,firstest)
	implicit none
	include 'parsizes.for'
	integer est,i,link(2,maxstack2),stack(maxstack)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer b,box,k,bwp,fwp,funbox
c
c	box=funbox(est,firstest)
	box=max(1,int(float(est)/firstest*(boxdim-1)))
c	case: middle occupant of box -> no action
	if(link_pointer(box).ne.i) return
	k=link(1,i)
c	case: top is last in box, last in stack
	if(k.eq.0) then
		link_pointer(box)=0
		bwp=backward_pointer(box)
		fwp=forward_pointer(box)
		do b=box+1,fwp
			backward_pointer(b)=bwp
		end do
		do b=bwp,box-1
			forward_pointer(b)=fwp
		end do
		return
	end if
	if(k.gt.0) b=max(1,int(float(stack(k))/firstest*(boxdim-1)))
c	case: last in box
	if(b.ne.box) then
		link_pointer(box)=0
		bwp=backward_pointer(box)
		fwp=forward_pointer(box)
		do b=box+1,fwp
			backward_pointer(b)=bwp
		end do
		do b=bwp,box-1
			forward_pointer(b)=fwp
		end do
		return
	end if
c	case: remainders in box
	if(b.eq.box) then
		link_pointer(box)=k
		return
	end if

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine toprequest(bestest,top,stack,chunk,link,mi,ns,seglist,
     $		ni,ci,est,link_pointer,forward_pointer,backward_pointer,
     $		firstest)
	implicit none
	include 'parsizes.for'
	integer bestest,top,stack(maxstack),chunk,link(2,maxstack2),ns
	integer ni(maxseg),ci(maxres0,maxseg),est,mi(maxseg),
     $ seglist(maxseg)
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer irray((1+(maxres0+30)/31)*maxseg),s,a
c
c	move best chunk to top of stack
c	translate chunk contents to ni,ci
c	decrement stack counter
c
c	call checkstack(stack,bestest,top,chunk,'top BUG 1')
	a=bestest
	if(top.ne.bestest) call swapstack(a,top,stack,chunk,link,bestest,
     $		link_pointer,forward_pointer,backward_pointer,firstest)
	if(lskipflag) return
c	check that top is really maximum!
c	call checkstack(stack,bestest,top,chunk,'top BUG 2')
	call getchunk(top,irray,stack,chunk)
	call translatechunk(irray,mi,ns,seglist,ni,ci,est,chunk)
	s=stack(bestest)
	call removebox(s,bestest,stack,link,link_pointer,
     $		forward_pointer,backward_pointer,firstest)
c	call checkstack(stack,bestest,top,chunk,'top BUG 3')
	call removetoplink(top,chunk,link,bestest)
c	call checkstack(stack,bestest,top,chunk,'top BUG 4')
c	write(*,*) 'toprequest',bestest,top,chunk,ns,est
c	write(*,*) 'ni',(ni(seglist(is)),is=1,ns)

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine checkstack(stack,bestest,top,chunk,text)
	implicit none
	include 'parsizes.for'
	integer bestest,top,stack(maxstack),chunk
	character*40 text
c
	integer i,m
c
!	type *,'checkstack',bestest,top,chunk,text
	if(bestest.gt.0) then
	  m=stack(bestest)
	  do i=1,top-1,chunk
		if(stack(i).gt.m) then
			write(*,*) top,bestest,stack(top),stack(bestest)
			write(*,*) (stack(m),m=1,top,chunk)
			write(*,*) text
	write(*,*) ' FATAL ERROR in checkstack'
			lskipflag=.true.
			return
		end if
	  end do
	end if
c
	return
	end
c
c----------------------------------------------------------------------
c
	subroutine removetoplink(top,chunk,link,bestest)
	implicit none
	include 'parsizes.for'
	integer link(2,maxstack2),bestest,top,chunk
c
	bestest=link(1,bestest)
	if(bestest.gt.0) link(2,bestest)=0
	top=top-chunk

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine translatechunk(irray,mi,ns,seglist,ni,ci,est,chunk)
	implicit none
	include 'parsizes.for'
	integer irray((1+(maxres0+30)/31)*maxseg),mi(maxseg)
	integer ns,est,ni(maxseg)
	integer chunk,seglist(maxseg),ci(maxres0,maxseg)
c
	integer i,j,ires,is,iseg,k
	logical l(31)
c
	est=irray(1)
	do is=1,ns
		ni(seglist(is))=0
	end do
	is=0
	ires=0
	do i=2,chunk
		call bitunpack(irray(i),l,1073741823,31)
		do j=1,31
			if(ires.eq.0) then
				is=is+1
				if(is.gt.ns) return
				iseg=seglist(is)
			end if
			ires=ires+1
			if(l(j)) then
				ni(iseg)=ni(iseg)+1
				ci(ni(iseg),iseg)=ires
			end if
c			zero ires counter before next segment
			if(ires.eq.mi(iseg)) ires=0
		end do
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine swaplinks_new(a,b,link,bestest)
	implicit none
	include 'parsizes.for'
	integer a,b,link(2,maxstack2),bestest
c
	integer dn1,dn2,up1,up2
c
	dn1=link(1,a)
	up1=link(2,a)
	dn2=link(1,b)
	up2=link(2,b)

	if(up1.eq.b.or.up2.eq.a) then	! adjacent --a-b--
		link(1,a)=up2
		link(1,b)=a
		if(dn1.ne.0) link(1,dn1)=b
		link(2,a)=b
		link(2,b)=dn1
		if(up2.ne.0) link(2,up2)=a
	else			! non-adjacent --a-x-b--
		link(1,b)=dn1
		link(1,a)=dn2
		link(2,b)=up1
		link(2,a)=up2
		if(dn1.ne.0) link(2,dn1)=a
		if(dn2.ne.0) link(2,dn2)=b
		if(up1.ne.0) link(1,up1)=b
		if(up2.ne.0) link(1,up2)=a
	end if

	if(up1.eq.0) bestest=b
	if(up2.eq.0) bestest=a

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine swaplinks(a,b,link,bestest)
	implicit none
	include 'parsizes.for'
	integer bestest,a,b,link(2,maxstack2)
c
	integer dn1,dn2,up1,up2
c
	dn1=link(1,a)
	up1=link(2,a)
	if(dn1.eq.b) dn1=a
	dn2=link(1,b)
	up2=link(2,b)
	if(up2.eq.a) up2=b
	link(1,b)=dn1
	link(2,b)=up1
	link(1,a)=dn2
	link(2,a)=up2
	if(dn1.ne.0) link(2,dn1)=b
	if(up1.ne.0) link(1,up1)=b
	if(dn2.ne.0) link(2,dn2)=a
	if(up2.ne.0) link(1,up2)=a
	if(up1.eq.0) bestest=b
	if(up2.eq.0) bestest=a

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine swapstack(from,to,stack,chunk,link,bestest,
     $ link_pointer,forward_pointer,backward_pointer,firstest)
	implicit none
	include 'parsizes.for'
	integer from,to,stack(maxstack),chunk,link(2,maxstack2),bestest
	integer link_pointer(0:boxdim),forward_pointer(0:boxdim),
     $		backward_pointer(0:boxdim)
	real firstest
c
	integer i,a,b,e,getk
c
c	write(*,*) 'swap',from,stack(from),to,stack(to),firstest
	if(from.eq.to) return
	e=stack(from)
	call removebox(e,from,stack,link,link_pointer,
     $		forward_pointer,backward_pointer,firstest)
	e=stack(to)
	call removebox(e,to,stack,link,link_pointer,
     $		forward_pointer,backward_pointer,firstest)
	a=from
	b=to
	do i=1,chunk
		e=stack(b)
		stack(b)=stack(a)
		stack(a)=e
		a=a+1
		b=b+1
	end do
	call swaplinks(from,to,link,bestest)
	e=stack(from)
	i=getk(e,from,stack,link_pointer,
     $		forward_pointer,backward_pointer,firstest)
	if(lskipflag) return
	e=stack(to)
	i=getk(e,to,stack,link_pointer,
     $		forward_pointer,backward_pointer,firstest)
	if(lskipflag) return

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine putchunk(to,irray,stack,chunk)
	implicit none
	include 'parsizes.for'
	integer to,stack(maxstack),chunk,irray(chunk)
c
	integer i,a
c
	a=to
	if(a+chunk.gt.maxstack) then
		write(*,*) to,chunk
	write(*,*)  'FATAL ERROR stack overflow - increase maxstack'
		lskipflag=.true.
		return
	end if
	do i=1,chunk
		stack(a)=irray(i)
		a=a+1
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine getchunk(from,irray,stack,chunk)
	implicit none
	include 'parsizes.for'
	integer from,stack(maxstack),chunk,irray(chunk)
c
	integer i,a
c
	a=from
	do i=1,chunk
		irray(i)=stack(a)
		a=a+1
	end do

	return
	end
c
c----------------------------------------------------------------------
c
	subroutine getnminnmax(ni,ns,seglist,nmin,nmax)
	implicit none
	include 'parsizes.for'
	integer ns,ni(maxseg),nmin,nmax,seglist(maxseg)
c
	integer is,n
c
	nmin=infinit
	nmax=-infinit
	do is=1,ns
		n=ni(seglist(is))
		nmin=min(nmin,n)
		nmax=max(nmax,n)
	end do

	return
	end
c
c----------------------------------------------------------------------
c

max_stars_repo_name: gkrash/lcspromatch

max_issues_repo_name: gkrash/lcspromatch

max_forks_repo_name: gkrash/lcspromatch


! { dg-do compile }
! { dg-options "-fwhole-file" }
! Test the fixes for the first two problems in PR40011
!
! Contributed by Dominique d'Humieres <dominiq@lps.ens.fr>
!
! This function would not compile because -fwhole-file would
! try repeatedly to resolve the function because of the self
! reference.
RECURSIVE FUNCTION eval_args(q)  result (r)
  INTEGER NNODE
  PARAMETER (NNODE  = 10)
  TYPE NODE
    SEQUENCE
    INTEGER car
    INTEGER cdr
  END TYPE NODE
  TYPE(NODE) heap(NNODE)
  INTEGER r, q
  r = eval_args(heap(q)%cdr)
END FUNCTION eval_args

function test(n)
  real, dimension(2) :: test
  integer            :: n
  test = n
  return
end function test

program arr     ! The error was not picked up causing an ICE
  real, dimension(2) :: res
  res = test(2) ! { dg-error "Explicit interface required" }
  print *, res
end program


text:
ccccccccccccccccccccccccccccccccccccccccccc
c Interpolate over the kth element of z by by using its 4 valid neighbors.

      subroutine interp_4(k,vnk,x,y,z,ezk)
      implicit none
c Inputs...
      integer k 		!mark to be filled
      integer vnk(4)		!valid neighbors of k
      real*4 x(202),y(202),z(202)
c Output...
      real*4 ezk		!estimated value of z(k)

c Local variables...
      real*4 a(4,4),b(4)
      integer i,ki,ind

c Model the behavior of z in the local image region as a bilinear
c transfomation using its four neighbors as data points.

      do i=1,4
         ki = vnk(i)
         a(i,1) = 1
         a(i,2) = x(ki)
         a(i,3) = y(ki)
         a(i,4) = x(ki)*y(ki)
         b(i) = z(ki)
      enddo

      call simq(a,b,4,ind)
      if (ind.eq.0) then
         ezk = b(1) + b(2)*x(k) + b(3)*y(k)
     &		+ b(4)*x(k)*y(k)
      else
         call prnt(2,1,k,'***simq err at mark.')
         ezk = -999
      endif

      return
      end

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


program demo_esc
! read stdin and run it through M_escape::esc to display without color
use M_escape, only : esc, esc_mode
implicit none
character(len=1024) :: line
integer :: ios
   line=''
   call esc_mode(manner='plain')
   do
      read(*,'(a)',iostat=ios)line
      write(*,'(a)') esc(trim(line))
      if(ios.ne.0)exit
   enddo
   write(*,'(a)',advance='no') esc('<reset>')
end program demo_esc


text:
stk.v.12.0
WrittenBy    STK_v12.0.0

BEGIN Facility

    Name		 Cape_Canaveral

    BEGIN CentroidPosition

        CentralBody		 Earth
        DisplayCoords		 Geodetic
        EcfLatitude		  2.8500000000000004e+01
        EcfLongitude		 -8.0570000000000007e+01
        EcfAltitude		 -2.8000000000000000e+01
        HeightAboveGround		  0.0000000000000000e+00
        ComputeTrnMaskAsNeeded		 Off
        SaveTrnMaskDataInBinary		 Off
        LightingGrazingType		 GroundModel
        DisplayAltRef		 Ellipsoid
        UseTerrainInfo		 Off
        NumAzRaysInMask		 360
        TerrainNormalMode		 UseCbShape

    END CentroidPosition

    BEGIN Extensions

        BEGIN LaserCAT
        END LaserCAT

        BEGIN ExternData
        END ExternData

        BEGIN RFI
        END RFI

        BEGIN ADFFileData
        END ADFFileData

        BEGIN AccessConstraints
            LineOfSight IncludeIntervals

            UsePreferredMaxStep No
            PreferredMaxStep 360
        END AccessConstraints

        BEGIN ObjectCoverage
        END ObjectCoverage

        BEGIN Desc
            BEGIN LongText
Name:           Cape Canaveral Air Force Station
Country:        UnitedStates
Location:       Cape Canaveral, Florida
Status:         Active
Type:           LaunchSite
Alternate name: CCAFS (Cape Canaveral Air Force Station)
Alternate name: AFETR (Air Force Eastern Test Range)
Launch Pads:    36 pads
Notes:          Approximate location. Cape Canaveral Air Force Station is part of the US Air Force Space Command's 45th Space Wing (45 SW), headquartered at nearby Patrick Air Force Base.
Sources:                       http://en.wikipedia.org/wiki/Cape_Canaveral_Air_Force_Station
Last updated:   2009-12-02
            END LongText
        END Desc

        BEGIN Atmosphere
<?xml version = "1.0" standalone = "yes"?>
<SCOPE>
    <VAR name = "InheritAtmosAbsorptionModel">
        <BOOL>true</BOOL>
    </VAR>
    <VAR name = "AtmosAbsorptionModel">
        <VAR name = "Simple_Satcom">
            <SCOPE Class = "AtmosphericAbsorptionModel">
                <VAR name = "Version">
                    <STRING>"1.0.1 a"</STRING>
                </VAR>
                <VAR name = "ComponentName">
                    <STRING>"Simple_Satcom"</STRING>
                </VAR>
                <VAR name = "Type">
                    <STRING>"Simple Satcom"</STRING>
                </VAR>
                <VAR name = "SurfaceTemperature">
                    <QUANTITY Dimension = "Temperature" Unit = "K">
                        <REAL>293.15</REAL>
                    </QUANTITY>
                </VAR>
                <VAR name = "WaterVaporConcentration">
                    <QUANTITY Dimension = "Density" Unit = "g*m^-3">
                        <REAL>7.5</REAL>
                    </QUANTITY>
                </VAR>
            </SCOPE>
        </VAR>
    </VAR>
    <VAR name = "EnableLocalRainData">
        <BOOL>false</BOOL>
    </VAR>
    <VAR name = "LocalRainIsoHeight">
        <QUANTITY Dimension = "DistanceUnit" Unit = "m">
            <REAL>2000</REAL>
        </QUANTITY>
    </VAR>
    <VAR name = "LocalRainRate">
        <QUANTITY Dimension = "SlowRate" Unit = "mm*hr^-1">
            <REAL>1</REAL>
        </QUANTITY>
    </VAR>
    <VAR name = "LocalSurfaceTemp">
        <QUANTITY Dimension = "Temperature" Unit = "K">
            <REAL>293.15</REAL>
        </QUANTITY>
    </VAR>
</SCOPE>        END Atmosphere

        BEGIN RadarCrossSection
<?xml version = "1.0" standalone = "yes"?>
<SCOPE>
    <VAR name = "Inherit">
        <BOOL>true</BOOL>
    </VAR>
</SCOPE>        END RadarCrossSection

        BEGIN RadarClutter
<?xml version = "1.0" standalone = "yes"?>
<SCOPE>
    <VAR name = "Inherit">
        <BOOL>true</BOOL>
    </VAR>
</SCOPE>        END RadarClutter

        BEGIN Identification
        END Identification

        BEGIN Crdn
        END Crdn

        BEGIN Graphics

            BEGIN Attributes

                MarkerColor		 #00ff00
                LabelColor		 #00ff00
                LineStyle		 0
                MarkerStyle		 5
                FontStyle		 0

            END Attributes

            BEGIN Graphics

                Show		 On
                Inherit		 On
                IsDynamic		 Off
                ShowLabel		 On
                ShowAzElMask		 Off
                ShowAzElFill		 Off
                AzElFillStyle		 7
                AzElFillAltTranslucency		 0.5
                UseAzElColor		 Off
                AzElColor		 #ffffff
                MinDisplayAlt		 0
                MaxDisplayAlt		 10000000
                NumAzElMaskSteps		 1
                ShowAzElAtRangeMask		 Off
                ShowAzElAtRangeFill		 Off
                AzElFillRangeTranslucency		 0.5
                AzElAtRangeFillStyle		 7
                UseAzElAtRangeColor		 Off
                AzElAtRangeColor		 #ffffff
                MinDisplayRange		 0
                MaxDisplayRange		 1000000
                NumAzElAtRangeMaskSteps		 10

                BEGIN RangeContourData
                    Show		 Off
                    ShowRangeFill		 Off
                    RangeFillTranslucency		 0.5
                    LabelUnits		 4
                    NumDecimalDigits		 3

                END RangeContourData

            END Graphics

            BEGIN DisplayTimes
                DisplayType		 AlwaysOn
            END DisplayTimes
        END Graphics

        BEGIN ContourGfx
            ShowContours		 Off
        END ContourGfx

        BEGIN Contours
            ActiveContourType		 Radar Cross Section

            BEGIN ContourSet Radar Cross Section
                Altitude		 0
                ShowAtAltitude		 Off
                Projected		 On
                Relative		 On
                ShowLabels		 Off
                LineWidth		 1
                DecimalDigits		 1
                ColorRamp		 On
                ColorRampStartColor		 #ff0000
                ColorRampEndColor		 #0000ff
                BEGIN ContourDefinition
                    BEGIN CntrAntAzEl
                        CoordinateSystem		 0
                        BEGIN AzElPatternDef
                            SetResolutionTogether		 0
                            NumAzPoints		 361
                            AzimuthRes		 1
                            MinAzimuth		 -180
                            MaxAzimuth		 180
                            NumElPoints		 91
                            ElevationRes		 1
                            MinElevation		 0
                            MaxElevation		 90
                        END AzElPatternDef
                    END CntrAntAzEl
                    BEGIN RCSContour
                        Frequency		 2997924580
                        ComputeType		 0
                    END RCSContour
                END ContourDefinition
            END ContourSet
        END Contours

        BEGIN VO
        END VO

        BEGIN 3dVolume
            ActiveVolumeType		 Radar Cross Section

            BEGIN VolumeSet Radar Cross Section
                Scale		 100
                MinimumDisplayedRcs		 1
                Frequency		  1.4500000000000000e+10
                ShowAsWireframe		 0
                BEGIN AzElPatternDef
                    SetResolutionTogether		 0
                    NumAzPoints		 50
                    AzimuthRes		 7.346938775510203
                    MinAzimuth		 -180
                    MaxAzimuth		 180
                    NumElPoints		 50
                    ElevationRes		 3.673469387755102
                    MinElevation		 0
                    MaxElevation		 180
                END AzElPatternDef
                ColorMethod		 1
                MinToMaxStartColor		 #ff0000
                MinToMaxStopColor		 #0000ff
                RelativeToMaximum		 0
            END VolumeSet
            BEGIN VolumeGraphics
                ShowContours		 No
                ShowVolume		 No
            END VolumeGraphics
        END 3dVolume

    END Extensions

    BEGIN SubObjects

    END SubObjects

END Facility

max_stars_repo_name: ningwersen/STKCodeExamples

max_issues_repo_name: ndelago/STKCodeExamples

max_forks_repo_name: ndelago/STKCodeExamples


text:
    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    ParseCSDataBlock.f90
    !> \brief   Parse the data block section of a ChemSage data-file.
    !> \author  M.H.A. Piro
    !> \date    Mar. 4, 2018
    !> \sa      ParseCSDataFile.f90
    !> \sa      ParseCSDataBlockGibbs.f90
    !> \sa      ParseCSDataBlockQKTO.f90
    !> \sa      ParseCSDataBlockRKMP.f90
    !> \sa      ParseCSDataBlockSUBL.f90
    !> \todo    I may have to revisit this source file if SUBG phases should be treated as
    !!           sublattice phases; this is not yet clear. At this time, it seems unnecessary
    !!           so I excluded them.
    !
    !
    ! DISCLAIMER
    ! ==========
    !
    ! All of the programming herein is original unless otherwise specified and is completely
    ! independent of ChemApp and related products, including Solgas, Solgasmix, Fact, FactSage
    ! and ChemSage.
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer      Description of change
    !   ----            ----------      ---------------------
    !   10/06/2011      M.H.A. Piro     Original code
    !   12/20/2012      M.H.A. Piro     Check that the solution phase type is supported.
    !   03/04/2018      M.H.A. Piro     Added capability to handle SUBG phases.
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to parse the "data block" section of a ChemSage data-file.
    !! The "data block" section of this data-file contains the thermodynamic data of all species and
    !! phases in the system, including the stoichiometry coefficients of all species and coefficients
    !! of the standard Gibbs energy equations.
    !!
    !! First, all solution phases are parsed, then all pure condensed phases and finally the fictive phases
    !! are parsed.  Within the solution phase loop, the solution phase name and type are first parsed, then
    !! the Gibbs energy equations of the pure species, and finally the mixing parameters are parsed.  Since
    !! the excess mixing parameters are model dependent, individual subroutines are assigned for each solution
    !! phase model.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    ! INFO                      A scalar integer that indicates a successful exit or identifies an error.
    ! nSpeciesCS                Number of species in the system (combined solution species and pure
    !                            separate phases).
    ! nGibbsEqSpecies           Number of Gibbs energy equations for a particular species.
    ! iSpeciesAtomsCS           Integer matrix representing the number of atoms of a particular
    !                            elements in a species.
    ! iParticlesPerMoleCS       An integer vector containing the number of particles per mole of the
    !                            constituent species formula mass.  The default value is 1.
    ! cSolnPhaseNameCS          The name of a solution phase.
    ! cSolnPhaseTypeCS          The type of a solution phase.
    ! cSolnPhaseTypeSupport     A character array representing solution phase types that are supported.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine ParseCSDataBlock

    USE ModuleParseCS

    implicit none

    integer               :: i, j, k, iCounterGibbsEqn, nCountSublattice
    real(8)               :: dDummy
    real(8),dimension(15) :: dTempVec

    ! Intialize variables:
    iCounterGibbsEqn    = 0
    iParticlesPerMoleCS = 1
    nParamCS            = 0
    dTempVec            = 0D0
    nParamPhaseCS       = 0
    nCountSublattice    = 0
    dGibbsMagneticCS    = 0D0

    ! Loop through all solution phases:
    LOOP_SolnPhases: do i = 1, nSolnPhasesSysCS

        ! Entry 1: Read name of phase:
        read (1,*,IOSTAT = INFO) cSolnPhaseNameCS(i)

        if (INFO /= 0) then
            INFO = 1100 + i
            return
        end if

        ! Entry 2: Read model name (i.e., type of solution phase):
        read (1,*,IOSTAT = INFO) cSolnPhaseTypeCS(i)

        if (INFO /= 0) then
            INFO = 1200 + i
            return
        end if

        ! Check to make sure that this solution phase type is supported:
        INFO = 123
        LOOP_SolnTypeSupport: do j = 1, nSolnTypeSupport
            if (cSolnPhaseTypeCS(i) == cSolnPhaseTypeSupport(j)) then
                ! This solution phase type is supported:
                INFO = 0
                exit LOOP_SolnTypeSupport
            end if
        end do LOOP_SolnTypeSupport

        ! Report an error if this solution phase type is not supported and return:
        if (INFO == 123) then
            INFO = 17
            return
        end if

        ! Check if the solution phase contains magnetic ordering terms, and if so,
        ! read in the terms:
        if ((cSolnPhaseTypeCS(i) == 'RKMPM').OR.(cSolnPhaseTypeCS(i) == 'SUBLM')) then
            ! Magnetic ordering terms are present.

            j = nSpeciesPhaseCS(i-1) + 1
            read (1,*,IOSTAT = INFO) dGibbsMagneticCS(j,3:4)

            ! Record an error if necessesary:
            if (INFO /= 0) then
                INFO = 1100 + i
                return
            end if

        elseif (cSolnPhaseTypeCS(i) == 'SUBG') then
            nSROPhasesCS = nSROPhasesCS + 1

            ! ---> I'm not sure why, but SUBG phases appear to have some strange magnetic terms (?).
            read (1,*,IOSTAT = INFO) dDummy

            ! Read in two integers representing the number of species and the number of pairs:
            read (1,*,IOSTAT = INFO) nPairsSROCS(nSROPhasesCS,1:2)

        end if

        ! Loop through species in solution phase:
        LOOP_SpeciesInSolnPhase: do j = nSpeciesPhaseCS(i-1) + 1, nSpeciesPhaseCS(i)

            ! SUBG phases contain a certain number of species, which are necessarily less
            ! than the number of pair fractions. The # of species indicated in the
            ! header file actually represents the number of pairs. Therefore, there are
            ! fewer species listed than what has been allocated.
            if (cSolnPhaseTypeCS(i) == 'SUBG') then
                if (j >= nSpeciesPhaseCS(i-1) + 1 + nPairsSROCS(nSROPhasesCS,1) ) then
                    exit LOOP_SpeciesInSolnPhase
                end if
            end if

            ! Store the magnetic ordering terms for each solution:
            k = nSpeciesPhaseCS(i-1) + 1
            dGibbsMagneticCS(j,3:4) = dGibbsMagneticCS(k,3:4)

            ! Store the phase index corresponding to the current species:
            iPhaseCS(j) = i

            ! The following subroutine parses the Gibbs energy equations (entries 3-5):
            call ParseCSDataBlockGibbs(i,j,iCounterGibbsEqn)

            ! Entry 6: Definition of temperature and pressure dependence terms (I don't understand the reasoning):
            if (cSolnPhaseTypeCS(i) == 'QKTO') then
                dTempVec = 0D0
                read (1,*,IOSTAT = INFO) dTempVec(1:2)
                dTempVec(1:2) = INT(dTempVec(1:2)) - (/1, 1/)

                if (SUM(dTempVec(1:2)) /= 0) then
                    INFO = 1600 + i
                    return
                end if
            end if

            if (INFO /= 0) return

        end do LOOP_SpeciesInSolnPhase

        ! Check the type of solution phase to interpret mixing parameters:
        select case (cSolnPhaseTypeCS(i))

            ! Ideal mixture model
            case ('IDMX')
                ! Do nothing.

            ! Quasichemical Kohler-Toop model
            case ('QKTO')
                call ParseCSDataBlockQKTO(i)

            ! Redlich-Kister-Muggiano-Polynomial model
            case ('RKMP', 'RKMPM')

                call ParseCSDataBlockRKMP(i)

            ! Compound Energy Formalism (sublattice) model:
            case ('SUBL', 'SUBLM')

                ! Count the number of phases with a sublattice:
                nCountSublattice      = nCountSublattice + 1
                iPhaseSublatticeCS(i) = nCountSublattice

                call ParseCSDataBlockSUBL(i, nCountSublattice)

            ! Quadruplet quasichemical model:
            case ('SUBG')

                ! Parse the data-block section for SUBG phases:
                call ParseCSDataBlockSUBG(i)

            case default

                ! The solution phase type is not supported. Report an error.
                INFO = 17
                return

        end select ! End checking the type of solution phase.

        ! Return if an error has been recorded:
        if (INFO /= 0) return

        ! Record the index of the mixing parameter for this phase:
        nParamPhaseCS(i) = nParamCS

    end do LOOP_SolnPhases      ! Variable i

    ! Recount the number of charged phases:
    nChargedPhaseCS = 0
    do i = 1, nSolnPhasesSysCS
        if ((cSolnPhaseTypeCS(i) == 'SUBL').OR.(cSolnPhaseTypeCS(i) == 'SUBLM')) nChargedPhaseCS = nChargedPhaseCS + 1
    end do

    allocate(iParamPassCS(nParamCS))
    iParamPassCS = 0

    ! Begin parsing pure condensed phases:
    LOOP_PureConPhase: do j = nSpeciesPhaseCS(nSolnPhasesSysCS) + 1, nSpeciesCS
        ! The phase index of a pure separate phase is set to zero:
        iPhaseCS(j) = 0

        call ParseCSDataBlockGibbs(iPhaseCS(j),j,iCounterGibbsEqn)

        if (INFO /= 0) exit LOOP_PureConPhase

    end do LOOP_PureConPhase

    return

end subroutine ParseCSDataBlock

max_stars_repo_name: nrc-fuels/thermochimica

max_issues_repo_name: nrc-fuels/thermochimica

max_forks_repo_name: nrc-fuels/thermochimica


text:
      subroutine fluxplot
c******************************************************************************
c     This subroutine creates plots of flux curves
c******************************************************************************
      end

max_stars_repo_name: DanielAndreasen/pymoog

max_issues_repo_name: DanielAndreasen/pymoog

max_forks_repo_name: DanielAndreasen/pymoog


module m_fields

    implicit none

    ! TODO: what is this variable
    real*8, allocatable :: fields(:, :, :, :)
    real*8, allocatable :: pressure_field(:, :, :, :)

!================================================================================
contains
!================================================================================

    subroutine m_fields_init

        use m_io
        use m_parameters
        implicit none

        integer :: n

        n = 3 + n_scalars + n_les

        allocate (fields(nx + 2, ny, nz, n), stat=ierr)

        if (ierr .ne. 0) then
            write (out, *) "Cannot allocate fields, stopping."
            call my_exit(-1)
        end if
        fields = zip

        allocate (pressure_field(nx + 2, ny, nz, 3), stat=ierr)

        if (ierr .ne. 0) then
            write (out, *) "Cannot allocate fields, stopping."
            call my_exit(-1)
        end if
        pressure_field = zip

        write (out, "('Allocated ',i3,' fields.')") n
        call flush (out)

        return
    end subroutine m_fields_init

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

    subroutine m_fields_exit
        use m_io
        implicit none

        if (allocated(fields)) deallocate (fields)

        write (out, *) 'fields deallocated.'
        call flush (out)

        return
    end subroutine m_fields_exit

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

!--------------------------------------------------------------------------------
!  Subroutine that broadcasts the array "fields" from the hydro part to the
!  "stats" part of the code
!--------------------------------------------------------------------------------
    subroutine fields_to_stats

        use m_openmpi
        use m_parameters
        use m_io
        implicit none

        integer :: n_field, n_proc, k, ratio, n_scalars_bcast

        ! broadcasting time from the hydro root process to the whole world
!!$    write(out,*) "Broadcasting fields to stats part."
!!$    call flush(out)

        ! first send it to the root process of the stats part
        count = 1
        tag = 0
        if (iammaster) then
            if (task .eq. 'hydro') call MPI_SEND(TIME, count, MPI_REAL8, id_root_stats, tag, MPI_COMM_WORLD, mpi_err)
            if (task .eq. 'stats') call MPI_RECV(TIME, count, MPI_REAL8, id_root_hydro, tag, MPI_COMM_WORLD, mpi_status, mpi_err)
!!$       write(out,*) "Exchanged information between master processors."
!!$       call flush(out)
        end if
        ! then broadcast it over the "stats" communicator
        if (task .eq. 'stats') then
            call MPI_BCAST(TIME, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
            ! checking if we need to start advancing scalars
            if (.not. int_scalars .and. TIME .gt. TSCALAR) then
                int_scalars = .true.
                write (out, *) "Starting to move the scalars."
                call flush (out)
            end if
        end if
!!$    write(out,*) "Broadcasted to slave processors."
!!$    call flush(out)

        ! figuring out how many scalars to broadcast:
        ! 0     if we do not move scalars
        ! all   if we move scalars
        n_scalars_bcast = 0
        if (int_scalars) n_scalars_bcast = n_scalars
!!$    write(out,*) "Number of scalars to broadcast:", n_scalars_bcast
!!$    call flush(out)

        if (numprocs_hydro .ge. numprocs_stats) then

            ! using the code structure:
            ! since the size of the arrays is always 2^n, there is always 2^k slabs
            ! of a "hydro" array that correspond to q slab of the "stats" array

            ratio = numprocs_hydro/numprocs_stats

            select case (task)

            case ('hydro')

                ! sending the fields array to the corresponding process in "stats"
                do n_field = 1, 3 + n_scalars_bcast
                    count = (nx + 2)*ny*nz
                    id_to = numprocs_hydro + floor(real(myid_world)/real(ratio))
                    tag = myid_world*(3 + n_scalars_bcast) + n_field - 1

!!$             write(out,*) "Sending :", n_field, count, id_to, tag
!!$             call flush(out)

                    call MPI_ISEND(fields(1, 1, 1, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                    call MPI_WAIT(mpi_request, mpi_status, mpi_err)

!!$             write(out,*) "Sent."
!!$             call flush(out)

                end do
            case ('stats')

                ! receiving fields from hydro processors
                do n_proc = 0, ratio - 1
                    id_from = myid*ratio + n_proc
                    count = (nx + 2)*ny*nz/ratio
                    k = (nz/ratio)*n_proc + 1
                    do n_field = 1, 3 + n_scalars_bcast
                        tag = id_from*(3 + n_scalars_bcast) + n_field - 1

!!$                write(out,*) "Receiving :", n_field, count, id_from, tag
!!$                call flush(out)

                        call MPI_IRECV(fields(1, 1, k, n_field), count, MPI_REAL8, &
                                       id_from, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                        call MPI_WAIT(mpi_request, mpi_status, mpi_err)

!!$                write(out,*) "Received."
!!$                call flush(out)

                    end do
                end do

            end select

        else

            ! now doing the same for the case when more processors are involved in
            ! the "stat" part than in "hydro" part

            ratio = numprocs_stats/numprocs_hydro

            select case (task)

            case ('hydro')

                do n_proc = 0, ratio - 1
                    count = (nx + 2)*ny*nz/ratio
                    id_to = numprocs_hydro + myid_world*ratio + n_proc
                    k = (nz/ratio)*n_proc + 1
                    do n_field = 1, 3 + n_scalars_bcast
                        tag = id_to*(3 + n_scalars_bcast) + n_field - 1
                        call MPI_ISEND(fields(1, 1, k, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                        call MPI_WAIT(mpi_request, mpi_status, mpi_err)
                    end do
                end do

            case ('stats')

                do n_field = 1, 3 + n_scalars_bcast
                    count = (nx + 2)*ny*nz
                    id_from = floor(real(myid)/real(ratio))
                    tag = myid_world*(3 + n_scalars_bcast) + n_field - 1
                    call MPI_RECV(fields(1, 1, 1, n_field), count, MPI_REAL8, &
                                  id_from, tag, MPI_COMM_WORLD, mpi_status, mpi_err)
                    call MPI_WAIT(mpi_request, mpi_status, mpi_err)
                end do

            end select

        end if

!!$    write(out,*) "broadcasted fields to stats"
!!$    call flush(out)

        return
    end subroutine fields_to_stats

!================================================================================
!================================================================================
!--------------------------------------------------------------------------------
!  Subroutine that broadcasts the arrays that contain velocities in the x-space
!  to the "parts" part of the code, for tracking particles
!  The velocities are contained in the arrays wrk1...3
!  The whole ideology remains similar to the subroutine fields_to_stats, except
!  for the fact that the parts part of the code receives the velocities into
!  the "fields" array.
!--------------------------------------------------------------------------------
    subroutine fields_to_parts

        use m_openmpi
        use m_parameters
        use m_io
        use m_work
        implicit none

        integer :: n_field, n_proc, k, ratio, n_scalars_bcast

        ! if there are zero particles, return
        if (nptot .eq. 0) return

        ! broadcasting time and timestep (dt) from the hydro root process
        ! first send it to the root process of the "parts" part
        count = 1
        if (iammaster) then
            tag = 0
            if (task .eq. 'hydro') call MPI_SEND(TIME, count, MPI_REAL8, id_root_parts, tag, MPI_COMM_WORLD, mpi_err)
            if (task .eq. 'parts') call MPI_RECV(TIME, count, MPI_REAL8, id_root_hydro, tag, MPI_COMM_WORLD, mpi_status, mpi_err)
            tag = 1
            if (task .eq. 'hydro') call MPI_SEND(dt, count, MPI_REAL8, id_root_parts, tag, MPI_COMM_WORLD, mpi_err)
            if (task .eq. 'parts') call MPI_RECV(dt, count, MPI_REAL8, id_root_hydro, tag, MPI_COMM_WORLD, mpi_status, mpi_err)
        end if
        ! then broadcast them over the "parts" communicator
        if (task .eq. 'parts') call MPI_BCAST(TIME, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
        if (task .eq. 'parts') call MPI_BCAST(dt, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)

        ! if have not yet started moving particles, return
        if (TIME .lt. starttime_particles) return

        ! if this is the first timestep when we need to start moving particles,
        ! change the int_particle variable
        if (TIME .ge. starttime_particles .and. .not. int_particles) then
            write (out, *) 'Starting to move particles'
            call flush (out)
            int_particles = .true.
        end if

        ! figure out if we broadcast scalars (currently not)
        n_scalars_bcast = 0
        ! if (int_scalars) n_scalars_bcast = n_scalars

        if (numprocs_hydro .ge. numprocs_parts) then

            ! using the code structure:
            ! since the size of the arrays is always 2^n, there is always 2^k slabs
            ! of a "hydro" array that correspond to q slab of the "parts" array

            ratio = numprocs_hydro/numprocs_parts

            select case (task)

            case ('hydro')

                ! sending the fields array to the corresponding process in "stats"
                do n_field = 1, 3 + n_scalars_bcast
                    count = (nx + 2)*ny*nz
                    id_to = id_root_parts + floor(real(myid_world)/real(ratio))
                    tag = myid_world*(3 + n_scalars_bcast) + n_field - 1

                    ! if the paricles are advected by fully resolved velocity
                    ! ( that is, particles_filter_size=0) then send the fully resolved
                    ! velocity to the "parts" task
!   Else, if the particles are advected by locally averaged velofity, send
!   the velocities in the Fourier form.  They will be locally averaged and
!   processed by the "parts" part of the code
                    if (particles_filter_size .le. 0.d0) then
                        call MPI_ISEND(wrk(1, 1, 1, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                    else
                        call MPI_ISEND(fields(1, 1, 1, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
!!$
!!$                write(out,*) id_to, n_field, fields(:,1,1,n_field)
!!$                call flush(out)
!!$
                    end if
                    call MPI_WAIT(mpi_request, mpi_status, mpi_err)
                end do
            case ('parts')

                ! receiving fields from hydro processors
                do n_proc = 0, ratio - 1
                    id_from = myid*ratio + n_proc
                    count = (nx + 2)*ny*nz/ratio
                    k = (nz/ratio)*n_proc + 1
                    do n_field = 1, 3 + n_scalars_bcast
                        tag = id_from*(3 + n_scalars_bcast) + n_field - 1

                        call MPI_IRECV(fields(1, 1, k, n_field), count, MPI_REAL8, &
                                       id_from, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                        call MPI_WAIT(mpi_request, mpi_status, mpi_err)
!!$
!!$                write(out,*) 'rec',id_from, n_field, fields(:,1,k,n_field)
!!$                call flush(out)
!!$
                    end do
                end do

            end select

        else

            ! now doing the same for the case when more processors are involved in
            ! the "parts" part than in "hydro" part

            ratio = numprocs_parts/numprocs_hydro

            select case (task)

            case ('hydro')

                do n_proc = 0, ratio - 1
                    count = (nx + 2)*ny*nz/ratio
                    id_to = id_root_parts + myid_world*ratio + n_proc
                    k = (nz/ratio)*n_proc + 1
                    do n_field = 1, 3 + n_scalars_bcast
                        tag = id_to*(3 + n_scalars_bcast) + n_field - 1
                        if (particles_filter_size .le. 0.d0) then
                           call MPI_ISEND(wrk(1, 1, k, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                        else
                        call MPI_ISEND(fields(1, 1, k, n_field), count, MPI_REAL8, id_to, tag, MPI_COMM_WORLD, mpi_request, mpi_err)
                        end if
                        call MPI_WAIT(mpi_request, mpi_status, mpi_err)
                    end do
                end do

            case ('parts')

                do n_field = 1, 3 + n_scalars_bcast
                    count = (nx + 2)*ny*nz
                    id_from = floor(real(myid)/real(ratio))
                    tag = myid_world*(3 + n_scalars_bcast) + n_field - 1
                    call MPI_RECV(fields(1, 1, 1, n_field), count, MPI_REAL8, &
                                  id_from, tag, MPI_COMM_WORLD, mpi_status, mpi_err)
                    call MPI_WAIT(mpi_request, mpi_status, mpi_err)
                end do

            end select

        end if

        return
    end subroutine fields_to_parts

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

!!$
!!$  subroutine check_bcast
!!$
!!$    use m_parameters
!!$    use m_io
!!$    use m_work
!!$    implicit none
!!$
!!$    integer :: i,j,k,n
!!$    real*8 :: zyu
!!$
!!$    ! defining the fields to be cosines (testing)
!!$    zyu = 1.0
!!$    if (task.eq.'hydro') then
!!$       do n = 1,3
!!$          do k = 1,nz
!!$             do j = 1,ny
!!$                do i = 1,nx
!!$
!!$
!!$                   fields(i,j,k,n) = sin(dble(n*i)*dx)
!!$
!!$                end do
!!$             end do
!!$          end do
!!$       end do
!!$    end if
!!$
!!$
!!$    call m_fields_bcast
!!$
!!$
!!$    do n = 1,3
!!$
!!$       if (task.eq.'hydro') write(fname,"('hydro',i1,'.arr')") n
!!$       if (task.eq.'stats') write(fname,"('stats',i1,'.arr')") n
!!$
!!$       tmp4(:,:,:) = fields(1:nx,:,:,n)
!!$       call write_tmp4
!!$    end do
!!$
!!$    return
!!$  end subroutine check_bcast
!!$

end module m_fields


text:
      DOUBLE PRECISION FUNCTION PHO_GRSF2(X,S,Alp,Bet,A,B,Ga,Gb,Gc,Gd,
     &   Ge,Gep)
      IMPLICIT NONE
      DOUBLE PRECISION A , Alp , B , Bet , Ga , Gb , Gc , Gd , Ge ,
     &                 Gep , S , X
      SAVE

      PHO_GRSF2 = (S*X**A*(Ga+Gb*SQRT(X)+Gc*X**B)
     &            +S**Alp*EXP(-Ge+SQRT(Gep*S**Bet*LOG(1.D0/X))))
     &            *(1.D0-X)**Gd

      END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
    !-------------------------------------------------------------------
    ! Overloaded definitions for (+):
    !

    ELEMENTAL FUNCTION add_rpe (x) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var) :: z
        z%sbits = significand_bits(x)
        z = +(x%val)
    END FUNCTION add_rpe

    ELEMENTAL FUNCTION add_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val + y%val
    END FUNCTION add_rpe_rpe

    ELEMENTAL FUNCTION add_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val + y
    END FUNCTION add_rpe_integer

    ELEMENTAL FUNCTION add_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val + y
    END FUNCTION add_rpe_long

    ELEMENTAL FUNCTION add_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val + y
    END FUNCTION add_rpe_real

    ELEMENTAL FUNCTION add_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val + y
    END FUNCTION add_rpe_realalt

    ELEMENTAL FUNCTION add_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x + y%val
    END FUNCTION add_integer_rpe

    ELEMENTAL FUNCTION add_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x + y%val
    END FUNCTION add_long_rpe

    ELEMENTAL FUNCTION add_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x + y%val
    END FUNCTION add_real_rpe

    ELEMENTAL FUNCTION add_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x + y%val
    END FUNCTION add_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (-):
    !

    ELEMENTAL FUNCTION sub_rpe (x) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var) :: z
        z%sbits = significand_bits(x)
        z = -(x%val)
    END FUNCTION sub_rpe

    ELEMENTAL FUNCTION sub_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val - y%val
    END FUNCTION sub_rpe_rpe

    ELEMENTAL FUNCTION sub_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val - y
    END FUNCTION sub_rpe_integer

    ELEMENTAL FUNCTION sub_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val - y
    END FUNCTION sub_rpe_long

    ELEMENTAL FUNCTION sub_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val - y
    END FUNCTION sub_rpe_real

    ELEMENTAL FUNCTION sub_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val - y
    END FUNCTION sub_rpe_realalt

    ELEMENTAL FUNCTION sub_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x - y%val
    END FUNCTION sub_integer_rpe

    ELEMENTAL FUNCTION sub_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x - y%val
    END FUNCTION sub_long_rpe

    ELEMENTAL FUNCTION sub_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x - y%val
    END FUNCTION sub_real_rpe

    ELEMENTAL FUNCTION sub_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x - y%val
    END FUNCTION sub_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (*):
    !

    ELEMENTAL FUNCTION mul_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val * y%val
    END FUNCTION mul_rpe_rpe

    ELEMENTAL FUNCTION mul_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val * y
    END FUNCTION mul_rpe_integer

    ELEMENTAL FUNCTION mul_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val * y
    END FUNCTION mul_rpe_long

    ELEMENTAL FUNCTION mul_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val * y
    END FUNCTION mul_rpe_real

    ELEMENTAL FUNCTION mul_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val * y
    END FUNCTION mul_rpe_realalt

    ELEMENTAL FUNCTION mul_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x * y%val
    END FUNCTION mul_integer_rpe

    ELEMENTAL FUNCTION mul_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x * y%val
    END FUNCTION mul_long_rpe

    ELEMENTAL FUNCTION mul_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x * y%val
    END FUNCTION mul_real_rpe

    ELEMENTAL FUNCTION mul_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x * y%val
    END FUNCTION mul_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (/):
    !

    ELEMENTAL FUNCTION div_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val / y%val
    END FUNCTION div_rpe_rpe

    ELEMENTAL FUNCTION div_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val / y
    END FUNCTION div_rpe_integer

    ELEMENTAL FUNCTION div_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val / y
    END FUNCTION div_rpe_long

    ELEMENTAL FUNCTION div_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val / y
    END FUNCTION div_rpe_real

    ELEMENTAL FUNCTION div_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val / y
    END FUNCTION div_rpe_realalt

    ELEMENTAL FUNCTION div_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x / y%val
    END FUNCTION div_integer_rpe

    ELEMENTAL FUNCTION div_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x / y%val
    END FUNCTION div_long_rpe

    ELEMENTAL FUNCTION div_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x / y%val
    END FUNCTION div_real_rpe

    ELEMENTAL FUNCTION div_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x / y%val
    END FUNCTION div_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (.GE.):
    !

    ELEMENTAL FUNCTION ge_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GE. y%val
    END FUNCTION ge_rpe_rpe

    ELEMENTAL FUNCTION ge_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GE. y
    END FUNCTION ge_rpe_integer

    ELEMENTAL FUNCTION ge_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GE. y
    END FUNCTION ge_rpe_long

    ELEMENTAL FUNCTION ge_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GE. y
    END FUNCTION ge_rpe_real

    ELEMENTAL FUNCTION ge_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GE. y
    END FUNCTION ge_rpe_realalt

    ELEMENTAL FUNCTION ge_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GE. y%val
    END FUNCTION ge_integer_rpe

    ELEMENTAL FUNCTION ge_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GE. y%val
    END FUNCTION ge_long_rpe

    ELEMENTAL FUNCTION ge_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GE. y%val
    END FUNCTION ge_real_rpe

    ELEMENTAL FUNCTION ge_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GE. y%val
    END FUNCTION ge_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (.LE.):
    !

    ELEMENTAL FUNCTION le_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LE. y%val
    END FUNCTION le_rpe_rpe

    ELEMENTAL FUNCTION le_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LE. y
    END FUNCTION le_rpe_integer

    ELEMENTAL FUNCTION le_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LE. y
    END FUNCTION le_rpe_long

    ELEMENTAL FUNCTION le_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LE. y
    END FUNCTION le_rpe_real

    ELEMENTAL FUNCTION le_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LE. y
    END FUNCTION le_rpe_realalt

    ELEMENTAL FUNCTION le_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LE. y%val
    END FUNCTION le_integer_rpe

    ELEMENTAL FUNCTION le_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LE. y%val
    END FUNCTION le_long_rpe

    ELEMENTAL FUNCTION le_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LE. y%val
    END FUNCTION le_real_rpe

    ELEMENTAL FUNCTION le_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LE. y%val
    END FUNCTION le_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (.GT.):
    !

    ELEMENTAL FUNCTION gt_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GT. y%val
    END FUNCTION gt_rpe_rpe

    ELEMENTAL FUNCTION gt_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GT. y
    END FUNCTION gt_rpe_integer

    ELEMENTAL FUNCTION gt_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GT. y
    END FUNCTION gt_rpe_long

    ELEMENTAL FUNCTION gt_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GT. y
    END FUNCTION gt_rpe_real

    ELEMENTAL FUNCTION gt_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .GT. y
    END FUNCTION gt_rpe_realalt

    ELEMENTAL FUNCTION gt_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GT. y%val
    END FUNCTION gt_integer_rpe

    ELEMENTAL FUNCTION gt_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GT. y%val
    END FUNCTION gt_long_rpe

    ELEMENTAL FUNCTION gt_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GT. y%val
    END FUNCTION gt_real_rpe

    ELEMENTAL FUNCTION gt_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .GT. y%val
    END FUNCTION gt_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (.LT.):
    !

    ELEMENTAL FUNCTION lt_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LT. y%val
    END FUNCTION lt_rpe_rpe

    ELEMENTAL FUNCTION lt_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LT. y
    END FUNCTION lt_rpe_integer

    ELEMENTAL FUNCTION lt_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LT. y
    END FUNCTION lt_rpe_long

    ELEMENTAL FUNCTION lt_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LT. y
    END FUNCTION lt_rpe_real

    ELEMENTAL FUNCTION lt_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val .LT. y
    END FUNCTION lt_rpe_realalt

    ELEMENTAL FUNCTION lt_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LT. y%val
    END FUNCTION lt_integer_rpe

    ELEMENTAL FUNCTION lt_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LT. y%val
    END FUNCTION lt_long_rpe

    ELEMENTAL FUNCTION lt_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LT. y%val
    END FUNCTION lt_real_rpe

    ELEMENTAL FUNCTION lt_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x .LT. y%val
    END FUNCTION lt_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (==):
    !

    ELEMENTAL FUNCTION eq_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val == y%val
    END FUNCTION eq_rpe_rpe

    ELEMENTAL FUNCTION eq_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val == y
    END FUNCTION eq_rpe_integer

    ELEMENTAL FUNCTION eq_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val == y
    END FUNCTION eq_rpe_long

    ELEMENTAL FUNCTION eq_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val == y
    END FUNCTION eq_rpe_real

    ELEMENTAL FUNCTION eq_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val == y
    END FUNCTION eq_rpe_realalt

    ELEMENTAL FUNCTION eq_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x == y%val
    END FUNCTION eq_integer_rpe

    ELEMENTAL FUNCTION eq_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x == y%val
    END FUNCTION eq_long_rpe

    ELEMENTAL FUNCTION eq_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x == y%val
    END FUNCTION eq_real_rpe

    ELEMENTAL FUNCTION eq_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x == y%val
    END FUNCTION eq_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (/=):
    !

    ELEMENTAL FUNCTION ne_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val /= y%val
    END FUNCTION ne_rpe_rpe

    ELEMENTAL FUNCTION ne_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        LOGICAL :: z
        z = x%val /= y
    END FUNCTION ne_rpe_integer

    ELEMENTAL FUNCTION ne_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val /= y
    END FUNCTION ne_rpe_long

    ELEMENTAL FUNCTION ne_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val /= y
    END FUNCTION ne_rpe_real

    ELEMENTAL FUNCTION ne_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        LOGICAL :: z
        z = x%val /= y
    END FUNCTION ne_rpe_realalt

    ELEMENTAL FUNCTION ne_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x /= y%val
    END FUNCTION ne_integer_rpe

    ELEMENTAL FUNCTION ne_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x /= y%val
    END FUNCTION ne_long_rpe

    ELEMENTAL FUNCTION ne_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x /= y%val
    END FUNCTION ne_real_rpe

    ELEMENTAL FUNCTION ne_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        LOGICAL :: z
        z = x /= y%val
    END FUNCTION ne_realalt_rpe

    !-------------------------------------------------------------------
    ! Overloaded definitions for (**):
    !

    ELEMENTAL FUNCTION pow_rpe_rpe (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val ** y%val
    END FUNCTION pow_rpe_rpe

    ELEMENTAL FUNCTION pow_rpe_integer (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER, INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val ** y
    END FUNCTION pow_rpe_integer

    ELEMENTAL FUNCTION pow_rpe_long (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        INTEGER(KIND=8), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val ** y
    END FUNCTION pow_rpe_long

    ELEMENTAL FUNCTION pow_rpe_real (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val ** y
    END FUNCTION pow_rpe_real

    ELEMENTAL FUNCTION pow_rpe_realalt (x, y) RESULT (z)
        TYPE(rpe_var), INTENT(IN) :: x
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x%val ** y
    END FUNCTION pow_rpe_realalt

    ELEMENTAL FUNCTION pow_integer_rpe (x, y) RESULT (z)
        INTEGER, INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x ** y%val
    END FUNCTION pow_integer_rpe

    ELEMENTAL FUNCTION pow_long_rpe (x, y) RESULT (z)
        INTEGER(KIND=8), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x ** y%val
    END FUNCTION pow_long_rpe

    ELEMENTAL FUNCTION pow_real_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_REAL_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x ** y%val
    END FUNCTION pow_real_rpe

    ELEMENTAL FUNCTION pow_realalt_rpe (x, y) RESULT (z)
        REAL(KIND=RPE_ALTERNATE_KIND), INTENT(IN) :: x
        TYPE(rpe_var), INTENT(IN) :: y
        TYPE(rpe_var) :: z
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        z = x ** y%val
    END FUNCTION pow_realalt_rpe

max_stars_repo_name: sparonuz/rpe

max_issues_repo_name: sparonuz/rpe

max_forks_repo_name: sparonuz/rpe


text:
C ===============
      SUBROUTINE ECOLLAPSE(EC)
C
C Front end for ECOLLAPSE2
C
C John Kuszewski 11/12/98
C ===============
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'collapse.inc'
      INCLUDE 'heap.inc'
C i/o
      DOUBLE PRECISION EC
C local vbls
C begin
      CALL ECOLLAPSE2(EC, HEAP(CLPSSELPTR))
      RETURN
      END
C ===============
      SUBROUTINE ECOLLAPSE2(EC, SELECTEDATOMS)
C
C Energy term for a simple collapse-inducing force
C
C Just forces the radius of gyration to be equal to
C the given value
C
C Perhaps I'll someday change this to handle mass-weighting
C and selections of atoms instead of all atoms
C
C JJK 6/14/96
C ===============
      IMPLICIT NONE
C i/o
      DOUBLE PRECISION EC
      INTEGER SELECTEDATOMS(*)
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'numbers.inc'
      INCLUDE 'collapse.inc'
      INCLUDE 'deriv.inc'
C local vbls
      DOUBLE PRECISION CENX, CENY, CENZ, SUM, N, RGYR
      DOUBLE PRECISION DEDXI, DEDYI, DEDZI
      INTEGER COUNT
C begin
C
C first, calculate the centroid of the structure
C
      N = NSELECTED
      CENX = ZERO
      CENY = ZERO
      CENZ = ZERO
      DO COUNT = 1, NSELECTED
           CENX = CENX + X(SELECTEDATOMS(COUNT))
           CENY = CENY + Y(SELECTEDATOMS(COUNT))
           CENZ = CENZ + Z(SELECTEDATOMS(COUNT))
      END DO
      CENX = CENX / N
      CENY = CENY / N
      CENZ = CENZ / N
C
C Now calculate radius of gyration
C
C Rg = RMS (xi - xcm)
C    = Sqrt(Sum(ri - rcm)^2)) / Sqrt(N)
C
C where N is the number of items in the sum
C
      SUM = ZERO
      DO COUNT = 1, NSELECTED
           SUM = SUM + (X(SELECTEDATOMS(COUNT)) - CENX)**2 +
     &                 (Y(SELECTEDATOMS(COUNT)) - CENY)**2 +
     &                 (Z(SELECTEDATOMS(COUNT)) - CENZ)**2
      END DO
      RGYR = SQRT(SUM) / SQRT(N)
      EC = KCOLLAPSE * (RGYR - RTARGET)**2
C
C Now calculate forces
C
      DO COUNT = 1, NSELECTED
           DEDXI = 2 * KCOLLAPSE * (RGYR - RTARGET) *
     &          (X(SELECTEDATOMS(COUNT)) - CENX) /
     &          (SQRT(N) * SQRT(SUM))
           DEDYI = 2 * KCOLLAPSE * (RGYR - RTARGET) *
     &          (Y(SELECTEDATOMS(COUNT)) - CENY) /
     &          (SQRT(N) * SQRT(SUM))
           DEDZI = 2 * KCOLLAPSE * (RGYR - RTARGET) *
     &          (Z(SELECTEDATOMS(COUNT)) - CENZ) /
     &          (SQRT(N) * SQRT(SUM))
           DX(SELECTEDATOMS(COUNT)) = DX(SELECTEDATOMS(COUNT)) + DEDXI
           DY(SELECTEDATOMS(COUNT)) = DY(SELECTEDATOMS(COUNT)) + DEDYI
           DZ(SELECTEDATOMS(COUNT)) = DZ(SELECTEDATOMS(COUNT)) + DEDZI
      END DO
      RETURN
      END
C ================
      SUBROUTINE READCOLLAPSE
C
C Reads the input for collapse-energy stuff.
C
C JJK 6/14/96
C ================
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'collapse.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'funct.inc'
C local vbls
C begin
C
      CALL PUSEND('COLLAPSE>')
      DO WHILE (.NOT.DONE)
           CALL NEXTWD('COLLAPSE>')
           CALL MISCOM('COLLAPSE>',USED)
           IF (.NOT.USED) THEN
C
                IF (WD(1:4).EQ.'HELP') THEN
C-DOCUMENTATION-SOURCE-BEGIN
                WRITE(DUNIT,'(8X,A)')
     &'COLLapse <collapse-statement> END ',
     &'<COLLapse-statement>:== ',
     &'   FORCe <real> | sets the force constant (default 1.0)',
     &'   TARGet <real> | sets the target radius of gyration',
     &'      | (default 1.0)',
     &'   SELEction <sel> | selects the atoms which are to be',
     &'      | used to determine Rgyr and which get the forces',
     &'      | (default: no atoms selected) '
C-DOCUMENTATION-SOURCE-END
C
C set the force constant
C
                ELSE IF (WD(1:4).EQ.'FORC') THEN
                     CALL NEXTF('force constant =', KCOLLAPSE)
C
C set the target radius of gyration
C
                ELSE IF (WD(1:4).EQ.'TARG') THEN
                     CALL NEXTF('target radius of gyration', RTARGET)
C
C get the selection of atoms to include
C
                ELSE IF (WD(1:4).EQ.'SELE') THEN
                     IF (CLPSSELPTR.EQ.0) THEN
                          CLPSSELPTR=ALLHP(INTEG4(NATOM))
                     END IF
                     CALL READCOLLAPSE2(HEAP(CLPSSELPTR))
C
C check for END statement
C
                ELSE
                     CALL CHKEND('COUPLINGS>', DONE)
                END IF
           END IF
      END DO
      DONE = .FALSE.
      RETURN
      END
C ==============
      SUBROUTINE READCOLLAPSE2 (SEL)
C
C Reads the selection of atoms to include in the collapse term
C
C John Kuszewski 11/12/98
C ==============
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'collapse.inc'
C i/o
      INTEGER SEL(*)
C local vbls
C begin
      CALL SELCTA(SEL, NSELECTED, X, Y, Z,.TRUE.)
      CALL MAKIND(SEL, NATOM, NSELECTED)
      RETURN
      END
C ==============
      SUBROUTINE COLLAPSEINIT
C
C initializes stuff for the collapse energy term
C
C JJK 6/14/96
C ==============
      IMPLICIT NONE
C include files
      INCLUDE 'collapse.inc'
C begin
      KCOLLAPSE = 1.0D0
      RTARGET = 1.0D0
      NSELECTED = 0
      CLPSSELPTR = 0
      RETURN
      END


text:
module probdata_module

!     Combustor variables
      double precision, save ::  pamb, Tinit, uinit
      double precision, save :: Tinner, Touter
      double precision, save :: Rsplit
      double precision, save :: swrlang

      double precision :: uin, uout

      double precision :: rho_l, rhoe_l

!     These help specify which specific problem
      integer        , save ::  probtype,idir

      double precision, save :: split(3)

end module probdata_module

max_stars_repo_name: wjge/PeleC

max_issues_repo_name: marient/PeleC

max_forks_repo_name: marient/PeleC