text:
c This module/program is part of DaliLite (c) L. Holm 1999
c
c
c-------------------------------------------------------------------------------
c
	subroutine mclean(imax0,ngene,genepool,gblocksize,blocksize,
     $		nres1,nres2,d1,d2,bestscore,lsave,initfragali,bestali1,
     $		itrimx,ikillx,bestresali1,score,killcutx,ind,lkill,lfrozen,
     $		outfragali,reflen,nclonex,expectscore)
	implicit none
	include 'gagasizes.for'
	integer imax,ngene,blocksize,gblocksize,nres1,nres2,itrimx,reflen
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real totscore,expectscore
	logical lsave
	integer*2 genepool(5,maxpair),initfragali(maxres,popsize),
     $		bestali1(maxres),outfragali(maxres)
	integer nclonex
c
c	input
c		imax0 -- max # Monte Carlo steps
c		itrimx -- trim every itrimx'th step
c		gblocksize -- used by genepool
c		blocksize -- blocksize used by alignment
c			N.B.: chops up genepool to genes of length blocksize
c		initfragali -- set of blocksize-genes to start from
c		ltop -- TRUE means sequential constraint
c		lsave -- TRUE means use rescore,dscore; FALSE means calculate all
c		ikillx -- call killer every ikillx'th step
c		bestresali1 -- saved alignments for killer
c		score -- scores of saved alignments
c		killcutx -- max allowed sequence identity for killer
c		ind -- identifies alignment for killer
c		lfrozen -- for killer
cc		reflen -- reference lenali for killer (in last refinement)
c		nclonex -- # clones to test in initfragali
c		expectscore -- return after 1st cycle if totscore is below
c	output
c		outfragali -- optimized set of genes
c		bestscore -- score of initfragali
c		bestali1 -- per residue alignment
c		lkill -- TRUE if alignment is converged (return if so)
c		lfrozen -- TRUE if alignment is converged
c
	integer cnt(maxres),prev1(0:maxres+1),next1(0:maxres+1)
	integer nextres(0:maxres+1),next2(0:maxres+1)
	integer*2 ali1(maxres),ali2(maxres),fragali1(maxres),fragali2(maxres)
 	logical lset(maxpair),ltest(maxres,maxres)
	integer ncand, candij(2,maxpair*4)
	integer ntest,testij(2,maxpair),test(maxres,-maxres:maxres)
	real rescore(maxres,maxres),dscore(maxres,maxres)
	integer ikillx,ind,imax0
	real score(popsize),killcutx
	integer*2 bestresali1(maxres,popsize)
	logical lkill,lfrozen(popsize)
c
c	local bookkeeping
c		cnt(i) -- multiplicity of genes at residue i
c		ali1(i),ali2(j) -- per residue alignment
c		fragali1(i),fragali2(j) -- per gene alignment
c		prev1(i),next1(i) -- lookup with respect to fragali1(i)
c		next2(j) -- lookup with respect to fragali2(j)
c		nextres(i) -- lookup with respect to ali1(i)
c		lset(ix) -- TRUE if gene is in fragali1
c		ltest(i,j) -- TRUE if residue pair has been appended to cand-list
c		ncand -- length of candidate (residue pairs i,j) list
c		candij(2,ix) -- pointer from candidate to residue pair
c		rescore(i,j) -- worth of residue pair in current alignment
c		dscore(i,j) -- how much residue pair would add to current ali
c		ntest -- length of active gene list
c		testij(2,it) -- pointer from gene list to first residue pair
c		test(i,j) -- pointer from genes first residue pair to gene list
c
c	make sure to delete incompatible residues before adding new ones !
c
c
c	functions
c
	real getp
	real ran
c
c	locals
c
	integer istep,nacc,ndel,ncand0,lenali
	logical lfr
	real alfa,delta,bestscore
	real dx,p
	integer i,j,k,ix
	integer nrem,rem(2,maxres),nnew,new(2,maxres),ngenedel,genedel(maxres)
	integer di(maxpair),ve(maxpair),q,imp,iclone
	real impscore
c
	integer jnk1(dim1)
	common /junk1/jnk1
	equivalence (jnk1,test)
	integer jnk23(dim23)
	common /junk23/jnk23
	equivalence (jnk23(1),rescore)
	equivalence (jnk23(maxres*maxres+1),dscore)
	integer jnk5(dim5)
	common /junk5/jnk5
	equivalence (jnk5(1),candij)
	equivalence (jnk5(8*maxpair+1),testij)
	equivalence (jnk5(10*maxpair+1),di)
	equivalence (jnk5(11*maxpair+1),ve)
	equivalence (jnk5(12*maxpair+1),lset)
	integer jnk6(dim6)
	common /junk6/jnk6
	equivalence (jnk6,ltest)
c
	imax=imax0
	if(lverb)write(*,*)'This is Monte Carlo',imax,ngene,blocksize,
     $		gblocksize,ltop,lsave
c
c	initialize
c
	bestscore=0.0
	call initialize(cnt,ali1,ali2,fragali1,fragali2,prev1,next1,
     $		ltest,ncand,rescore,dscore,nres1,nres2,ntest,nextres,test,
     $		next2)
c
c	load initali
c
	lenali=0
	totscore=0.0
	do i=1,nres1
		j=initfragali(i,1)
		if(j.ne.0) then
			call appendcand(i,j,blocksize,
     $				lset,ltest,dscore,ncand,
     $				candij,ali1,ali2,nres1,d1,d2,
     $				nextres,lsave,ntest,test,testij)
c			call printcandlist(ncand,candij,dscore,rescore)
			call testaddition(i,j,lsave,blocksize,prev1,next1,
     $				ali1,cnt,nres1,d1,d2,rescore,dscore,
     $				dx,nrem,rem,nnew,new,ngenedel,genedel,
     $				nres2,fragali1,fragali2,nextres,next2)

			call changeconfig(ngenedel,genedel,nrem,rem,i,j,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)

		end if
	end do
	if(lverb)write(*,*) 'initial totscore',totscore
c!!!	call gettotscore(ali1,d1,d2,nres1,totscore)
c	write(*,500) (initfragali(i,1),i=1,nres1)
c	write(*,500) (ali1(i),i=1,nres1)
c
c	Monte Carlo step
c
       do iclone=1,nclonex
	lfr=.false.
	ncand0=0
	alfa=50.0
	delta=0.0
	imp=0
	istep=0
	impscore=0.0
	do while(istep.lt.imax)
		istep=istep+1
c
c		kill !
c
		if(mod(istep,ikillx).eq.1) then
			call killer(bestresali1,ali1,bestscore/100,nres1,
     $				blocksize,score,killcutx,lfrozen,lkill,
     $				ind,reflen)
			if(lkill) then
				if(lverb)write(*,*) ' killing converged ',ind
				return
			end if
		end if
c
c		heatshock if score improvement is stalled
c
		if(istep-imp.eq.1) then
			alfa=50.0
			delta=0.0
		else if(istep-imp.eq.2) then
			alfa=1.0/sqrt(max(10000.0,totscore))
			delta=alfa/1000
		end if
		ndel=0
		nacc=0
		alfa=alfa+delta
c
c		pull more genes as long as find more candidates
c
		if(.not.lfr) then
			call pullgene(gblocksize,blocksize,ngene,genepool,
     $				lset,ltest,dscore,ncand,candij,fragali1,next1,
     $				ali1,ali2,nres1,nres2,d1,d2,nextres,
     $				lsave,ntest,test,testij)
c			call printcandlist(ncand,candij,dscore,rescore)
			if(ncand.eq.ncand0) then
				lfr=.true.
			else
				imax=imax+1
			end if
			ncand0=ncand
		end if
c
c		adding phase  WARM
c
c
c		test in randomized order
c
		do ix=1,ntest
			di(ix)=ix
			ve(ix)=nint(100.0*ran(seed))
		end do
		call j_index_Qsort(1,ntest,maxpair,di,ve)
		do q=1,ntest
			ix=di(q)
			i=testij(1,ix)
			j=testij(2,ix)
			if(fragali1(i).ne.j) then
				if(fragali1(i).eq.j) then
					write(*,*) 'WARNING: duplicate gene',i,j
					goto 19
     				end if
				call testaddition(i,j,lsave,blocksize,
     $					prev1,next1,ali1,cnt,nres1,d1,d2,
     $					rescore,dscore,dx,nrem,rem,nnew,new,
     $					ngenedel,genedel,nres2,
     $					fragali1,fragali2,nextres,next2)
				p=getp(alfa*dx)
				if(ran(seed).lt.p) then
c				  write(*,*)'accepted',i,j,dx

				  call changeconfig(ngenedel,genedel,
     $		nrem,rem,i,j,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)
				  if(totscore.gt.impscore) then
					imp=istep
					impscore=totscore
				  end if
				end if
19			end if
		end do
c!!!		call gettotscore(ali1,d1,d2,nres1,totscore)
c
c		trimming phase COLD
c
		if(mod(istep,itrimx).eq.0) then
		 i=0
		 do while(next1(i).le.nres1)
		  i=next1(i)
		  j=fragali1(i)
c		  Trim from ends only
		  if((cnt(i).eq.1).or.(cnt(i+blocksize-1).eq.1))then
			ngenedel=1
			genedel(1)=i
			call testdeletion(ngenedel,genedel,blocksize,cnt,
     $				rescore,ali1,nextres,nres1,d1,d2,lsave,dx,
     $				nrem,rem)
c			p=getp(alfa*dx)
c			if(ran(seed).lt.p) then
c			only accept improvements !
			if(dx.gt.0.0) then
c				write(*,*)'trimming',i,j,dx

				  call changeconfig(ngenedel,genedel,
     $		nrem,rem,0,0,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)
				  if(totscore.gt.impscore) then
					imp=istep
					impscore=totscore
				  end if

			end if
		  end if
		 end do
		end if
		if(lverb)write(*,510) istep,imp,nacc,ndel,lenali,ncand,
     $			totscore,bestscore,alfa
c!!!		call gettotscore(ali1,d1,d2,nres1,totscore)
c		if(lsave)call print(ncand,candij,dscore,rescore,ali1,
c     $			next1,d1,d2,nres1)
c
c		compare to expected score
c
		if(totscore.lt.expectscore) return
	end do
c
c	mutate fragali1 to initfragali(iclone+1)
c
	if(iclone.lt.nclonex) then
		if(lverb)write(*,*) ' clone',iclone
		if(lverb)write(*,500) (fragali1(i),i=1,nres1)
		if(lverb)write(*,500) (initfragali(i,iclone+1),i=1,nres1)
c
c		additions, replacements
c
		do i=1,nres1
		  j=initfragali(i,iclone+1)
		  k=fragali1(i)
		  if((j.ne.0).and.(j.ne.k)) then
				call testaddition(i,j,lsave,blocksize,
     $					prev1,next1,ali1,cnt,nres1,d1,d2,
     $					rescore,dscore,dx,nrem,rem,nnew,new,
     $					ngenedel,genedel,nres2,
     $					fragali1,fragali2,nextres,next2)
c				write(*,*) 'add',i,j
c				write(*,500) (rem(1,ix),rem(2,ix),ix=1,nrem)
c				write(*,500) (new(1,ix),new(2,ix),ix=1,nnew)

				call changeconfig(ngenedel,genedel,
     $		nrem,rem,i,j,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)

		  end if
		end do
c!!!		call gettotscore(ali1,d1,d2,nres1,totscore)
c
c		deletions
c
		do i=1,nres1
		  j=initfragali(i,iclone+1)
		  k=fragali1(i)
		  if((k.ne.0).and.(k.ne.j)) then
			ngenedel=1
			genedel(1)=i
			call testdeletion(ngenedel,genedel,blocksize,cnt,
     $				rescore,ali1,nextres,nres1,d1,d2,lsave,dx,
     $				nrem,rem)
c			write(*,*) 'delete',i,k
c			write(*,500) (rem(1,ix),rem(2,ix),ix=1,nrem)

			call changeconfig(ngenedel,genedel,
     $		nrem,rem,0,0,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)

		  end if
		end do
c!!!		call gettotscore(ali1,d1,d2,nres1,totscore)
	end if
	imax=imax0
       end do
c
500	format(20i4)
510	format('step',6i6,2f10.2,f10.5)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine changeconfig(ngenedel,genedel,nrem,rem,i,j,
     $		nnew,new,fragali1,fragali2,ali1,ali2,next1,prev1,next2,
     $		nextres,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,ndel,nacc,test,blocksize,lsave,totscore,bestscore,
     $		nres1,outfragali,bestali1,dx)
c
c	input:
c		ngenedel,genedel -- genes to delete
c		i,j -- gene to add
c		nrem,rem -- residues to delete
c		nnew,new -- residues to add
c		istep -- for imp
c		dx -- change in totscore
c
c	output:
c		imp -- updated if score improved
c		totscore -- updated
c		outfragali -- bestever fragali1 updated
c		bestali1 -- bestever ali1 updated
c
	implicit none
	include 'gagasizes.for'
	integer ngenedel,genedel(maxres),nrem,rem(2,maxres),i,j,nnew
	integer new(2,maxres)
	integer*2 fragali1(maxres),fragali2(maxres),ali1(maxres),ali2(maxres)
	integer next1(0:maxres+1),prev1(0:maxres+1),next2(0:maxres+1)
	integer nextres(0:maxres+1),candij(2,maxpair*4),ncand,lenali
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real rescore(maxres,maxres)
	real dscore(maxres,maxres),totscore,bestscore,dx
	integer cnt(maxres),ndel,nacc,test(maxres,-maxres:maxres),blocksize
	logical lset(maxpair),lsave
	integer nres1
	integer*2 outfragali(maxres),bestali1(maxres)
c
	integer k
c
	call dodeletions(ngenedel,genedel,nrem,rem,fragali1,fragali2,ali1,ali2,
     $		next1,prev1,lset,candij,ncand,d1,d2,rescore,dscore,
     $		lenali,cnt,ndel,lsave, test,nextres,blocksize,next2)
	if(i.ne.0) call doadditions(i,j,nnew,new,fragali1,fragali2,ali1,ali2,
     $		next1,prev1,lset,candij,ncand,d1,d2,rescore,dscore,lenali,
     $		cnt,blocksize,nacc,lsave,test,nextres,next2,nres1)
	totscore=totscore+dx
	if(ldebug) then
		write(*,*) 'totscore now',totscore,nacc,ndel
		write(*,500) (ali1(k),k=1,nres1)
		write(*,500) (fragali1(k),k=1,nres1)
		write(*,500) (cnt(k),k=1,nres1)
		write(*,500) (nextres(k),k=1,nres1)
c		if(lsave)call print(ncand,candij,dscore,
c     $			rescore,ali1,next1,d1,d2,nres1)
		write(*,500) (rem(1,k),rem(2,k),k=1,nrem)
		write(*,500) (new(1,k),new(2,k),k=1,nnew)
		call gettotscore(ali1,d1,d2,nres1,totscore)
		call checkcnt(cnt,fragali1,nres1,blocksize)
	end if
	if(totscore.gt.bestscore) then
		call saveali(nres1,fragali1,outfragali,totscore,bestscore)
		call saveali(nres1,ali1,bestali1,totscore,bestscore)
	end if
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine print(ncand,candij,dscore,rescore,
     $		ali1,next1,d1,d2,nres1)
	implicit none
	include 'gagasizes.for'
	integer ncand,candij(2,maxpair*4),next1(0:maxres+1),nres1
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real dscore(maxres,maxres)
	real rescore(maxres,maxres)
	integer*2 ali1(maxres)
c
	real getrescore
c
	integer k,ix,jx
	real x,a,b
c
	write(*,*) 'This is print',ncand
	do k=1,ncand
		ix=candij(1,k)
		jx=candij(2,k)
		x=getrescore(ix,jx,d1,d2,ali1,next1,nres1)
		a=dscore(ix,jx)
		b=rescore(ix,jx)
		if(abs(a+b-x).gt.0.1) write(*,*)'dscore',ix,jx,a,b,x,jx-ali1(ix)
	end do
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine dodeletions(ngenedel,genedel,nrem,rem,fragali1,fragali2,
     $		ali1,ali2,next1,prev1,lset,candij,ncand,d1,d2,
     $		rescore,dscore,lenali,cnt,ndel,lsave,test,nextres,blocksize,
     $		next2)
	implicit none
	include 'gagasizes.for'
	integer ngenedel,genedel(maxres),nrem,rem(2,maxres)
	integer next2(0:maxres+1)
	integer*2 fragali1(maxres),fragali2(maxres),ali1(maxres),ali2(maxres)
	integer next1(0:maxres+1),prev1(0:maxres+1),cnt(maxres),blocksize
	logical lset(maxpair)
	integer candij(2,maxpair*4),ncand
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real rescore(maxres,maxres),dscore(maxres,maxres)
	integer lenali,ndel,test(maxres,-maxres:maxres),nextres(0:maxres+1)
	logical lsave
c
	real scorefun
c
	integer irem,i,j,j1,k,l,ix,jx,in,ip
	real s,dx
c
	if(ldebug)then
		write(*,*) 'This is dodeletions',ngenedel
		write(*,500) (genedel(i),i=1,ngenedel)
		write(*,500) (rem(1,i),rem(2,i),i=1,nrem)
		write(*,500) (ali1(i),i=1,38)
		write(*,500) (nextres(i),i=0,39)
	end if
c
c	delete genes
c
	do irem=1,ngenedel
		i=genedel(irem)
		j=fragali1(i)
		fragali1(i)=0
		fragali2(abs(j))=0
		lset(test(i,j))=.false.
		do l=0,blocksize-1
			cnt(i+l)=cnt(i+l)-1
			if(cnt(i+l).lt.0) write(*,*)
     $				' WARNING: subzero cnt ' ,i,j,i+l,j+l
		end do
c
c		update prev1,next
c
		ip=prev1(i)
		in=next1(i)
		do ix=i,in
			prev1(ix)=ip
		end do
		do ix=ip,i
			next1(ix)=in
		end do
		j1=abs(j)
		in=next2(j1)
		ix=j1
		do while((ix.gt.0).and.(next2(ix-1).eq.j1))
			ix=ix-1
			next2(ix)=in
		end do
	end do
c
c	delete residues
c
	do irem=1,nrem
		i=rem(1,irem)
		j=rem(2,irem)
c
c		update nextres
c
		in=nextres(i)
		ix=i
c		write(*,*) 'nextres',in,i,nextres(ix-1)
		do while((ix.gt.0).and.(nextres(ix-1).eq.i))
			ix=ix-1
			nextres(ix)=in
c			if(ldebug)write(*,*) 'set nextres',ix,' to ',in
		end do
		if(lsave) then
c
c			update dscore(x,y), rescore(x,y)
c
			do k=1,ncand
				ix=candij(1,k)
				jx=candij(2,k)
				dx=0.0
				if((i.ne.ix).and.(j.ne.jx)) then
					s=-scorefun(d1(ix,i),d2(jx,j))
					dx=s+s
				end if
				if(ali1(ix).eq.jx) then
					rescore(ix,jx)=rescore(ix,jx)+dx
				else
					dscore(ix,jx)=dscore(ix,jx)+dx
				end if
			end do
			dscore(i,j)=rescore(i,j)
			rescore(i,j)=0.0
		end if
c
c		delete residue pair, update lenali
c
		ali1(i)=0
		ali2(j)=0
		lenali=lenali-1
		ndel=ndel+1
	end do
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine doadditions(i,j,nnew,new,fragali1,fragali2,
     $		ali1,ali2,next1,prev1,lset,candij,ncand,d1,d2,
     $		rescore,dscore,lenali,cnt,blocksize,nacc,lsave,test,nextres,
     $		next2,nres1)
	implicit none
	include 'gagasizes.for'
	integer nnew,new(2,maxres),i,j,nres1
	integer*2 fragali1(maxres),fragali2(maxres),ali1(maxres),ali2(maxres)
	integer next1(0:maxres+1),prev1(0:maxres+1),cnt(maxres)
	logical lset(maxpair)
	integer candij(2,maxpair*4),ncand,test(maxres,-maxres:maxres)
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real rescore(maxres,maxres),dscore(maxres,maxres)
	integer lenali,blocksize,nacc,nextres(0:maxres+1),next2(0:maxres+1)
	logical lsave
c
	real scorefun
c
	integer inew,k,l,i1,j1,ix,jx,in,ip,iold,jold
	real s,dx
c
	if(ldebug)then
		write(*,*) 'This is doadditions',i,j,nnew,test(i,j)
		write(*,500) (new(1,k),new(2,k),k=1,nnew)
		write(*,500) (ali1(k),k=1,38)
		write(*,500) (nextres(k),k=0,39)
	end if
c
c	add gene
c
	fragali1(i)=j
	if(j.lt.0) fragali2(-j)=-i
	if(j.gt.0) fragali2(j)=i
	lset(test(i,j))=.true.
	do l=0,blocksize-1
		cnt(i+l)=cnt(i+l)+1
		if(cnt(i+l).gt.blocksize) write(*,*)
     $			' WARNING: excess cnt',i+l,j+l
	end do
c
c	update prev1,next
c
	ip=prev1(i)
	in=next1(i)
	do ix=ip,i-1
		next1(ix)=i
	end do
	do ix=i+1,in
		prev1(ix)=i
	end do
	j1=abs(j)
	in=next2(j1)
	ix=j1
	do while((ix.gt.0).and.(next2(ix-1).eq.in))
		ix=ix-1
		next2(ix)=j1
	end do
c
c	add residues
c
	do inew=1,nnew
		i1=new(1,inew)
		j1=new(2,inew)
c
c		update nextres
c
		in=nextres(i1)
		ix=i1
c		write(*,*) 'nextres',in,i,nextres(ix-1)
		do while((ix.gt.0).and.(nextres(ix-1).eq.in))
			ix=ix-1
			nextres(ix)=i1
c			if(ldebug) write(*,*) 'set nextres',ix,' to ',i1
		end do
		if(lsave) then
c
c			update dscore(x,y), rescore(x,y)
c
			do k=1,ncand
				ix=candij(1,k)
				jx=candij(2,k)
				dx=0.0
				if((i1.ne.ix).and.(j1.ne.jx)) then
					s=scorefun(d1(ix,i1),d2(jx,j1))
					dx=s+s
				end if
				if(ali1(ix).eq.jx) then
					rescore(ix,jx)=rescore(ix,jx)+dx
				else
					dscore(ix,jx)=dscore(ix,jx)+dx
				end if
				if(ldebug) write(*,510) k,i1,j1,ix,jx,
     $					dscore(ix,jx),rescore(ix,jx),dx
			end do
			rescore(i1,j1)=dscore(i1,j1)
			dscore(i1,j1)=0.0
		end if
c
c		add residue pair, update lenali
c
		iold=ali2(j1)
		jold=ali1(i1)
		if((iold.ne.0).or.(jold.ne.0)) then
			write(*,*) ' WARNING: undeleted ',i1,j1,iold,jold
			write(*,500) (ali1(k),k=1,nres1)
			write(*,500) (fragali1(k),k=1,nres1)
		end if
		if(iold.ne.0) ali1(iold)=0
		if(jold.ne.0) ali2(jold)=0
		ali1(i1)=j1
		ali2(j1)=i1
		lenali=lenali+1
		nacc=nacc+1
	end do
c
500	format(20i4)
510	format('update',5i5,3f10.3)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine testaddition(i,j,lsave,blocksize,prev1,next1,ali1,cnt,
     $		nres1,d1,d2,rescore,dscore,dx,nrem,rem,nnew,new,ngenedel,
     $		genedel,nres2,fragali1,fragali2,nextres,next2)
	implicit none
	include 'gagasizes.for'
	integer  ngenedel,genedel(maxres)
	integer i,j,nrem,rem(2,maxres),nnew,new(2,maxres),blocksize
	integer prev1(0:maxres+1),next1(0:maxres+1),cnt(maxres),nres1,nres2
	integer*2 ali1(maxres),fragali1(maxres),fragali2(maxres)
	logical lsave
	real dx
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real rescore(maxres,maxres),dscore(maxres,maxres)
	integer nextres(0:maxres+1),next2(0:maxres+1)
c
	real getrescore
	real scorefun
c
c	in:
c		i,j -- gene to be added
c		ltop -- TRUE means topology constraint
c	out:
c		dx -- change in score on addition
c	    	nrem,rem -- list of residues to delete
c		nnew,new -- list of residues to add
c		ngenedel,genedel -- list of genes to delete
c
c		N.B.: nrem,nnew are initialized to zero
c
	integer k,l,inew,i1,i2
	real s
c
	if(ldebug) then
		write(*,*) 'This is testaddition',i,j,ltop,lsave
		write(*,500) (ali1(k),k=1,nres1)
	end if
c
c	initialize
c
	ngenedel=0
	nrem=0
	nnew=0
c
c	optional topology check in range <i-blocksize+1, > i+blocksize-1
c		both in i- and j-directions
c
	if(ltop)call checktopo(i,j,blocksize,prev1,next1,fragali1,ngenedel,
     $		genedel,nres1)
c
c	which genes overlap with gene (i,j) ?
c	check range i-blocksize+1 ... i+blocksize-1
c                   j-blocksize+1 ... j+blocksize-1
c
	call geneoverlap(i,j,fragali1,fragali2,next1,blocksize,
     $		nres1,nres2,ngenedel,genedel,next2)
c
c	which residues go with genedel ?
c
	call gowithgenes(i,j,ngenedel,genedel,ali1,cnt,blocksize,
     $		nrem,rem,.true.)
c
c	new residues
c
	do l=0,blocksize-1
		if(ali1(i+l).ne.abs(j+l)) then
			nnew=nnew+1
			new(1,nnew)=i+l
			new(2,nnew)=abs(j+l)
		end if
	end do
c
c	if(ldebug) then
c		write(*,*) 'genedel:',ngenedel
c		write(*,500) (genedel(k),k=1,ngenedel)
c		write(*,*) 'rem:',nrem
c		write(*,500) (rem(1,k),rem(2,k),k=1,nrem)
c		write(*,*) 'new:',nnew
c		write(*,500) (new(1,k),new(2,k),k=1,nnew)
c	end if
c
c	calculate dx
c
c	add new
c
	dx=0.0
	do inew=1,nnew
	  i1=new(1,inew)
	  i2=new(2,inew)
	  if(i1.le.maxres.and.i2.le.maxres) then
		if(lsave) then
			dx=dx+dscore(i1,i2)
                        if(ldebug) then
                          write(*,*) 'add inew',inew,i1,i2,dscore(i1,i2),dx
                        end if
		else
			dx=dx+getrescore(i1,i2,d1,d2,ali1,nextres,nres1)
		end if
c
c		add new x new
c
		do k=1,inew-1
			s=scorefun(d1(i1,new(1,k)),d2(i2,new(2,k)))
			dx=dx+s+s
		        if(ldebug) then
                         write(*,*) 'add new*new ',inew,k,i1,i2,
     $	new(1,k),new(2,k),s+s,dx,d1(i1,new(1,k)),d2(i2,new(2,k))
                        end if
		end do
c
c		subtract new x rem
c
		do k=1,nrem
			if((i1.ne.rem(1,k)).and.(i2.ne.rem(2,k))) then
				s=scorefun(d1(i1,rem(1,k)),d2(i2,rem(2,k)))
				dx=dx-s-s
                                if(ldebug) then
                           write(*,*) 'subtract new*old ',inew,k,i1,i2,
     $	new(1,k),new(2,k),-s-s,dx,d1(i1,rem(1,k)),d2(i2,rem(2,k))
                                end if
			end if
		end do
	  end if
	end do
c
c	subtract rem x rem
c
	do k=1,nrem
		if(lsave) then
			dx=dx-rescore(rem(1,k),rem(2,k))
		else
			dx=dx-getrescore(rem(1,k),rem(2,k),d1,d2,
     $				ali1,nextres,nres1)
		end if
		do l=1,k-1
			s=scorefun(d1(rem(1,k),rem(1,l)),
     $				d2(rem(2,k),rem(2,l)))
			dx=dx+s+s
		end do
	end do
c
	if(ldebug)write(*,*) 'dx: ',dx
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine geneoverlap(i,j,fragali1,fragali2,next1,blocksize,
     $		nres1,nres2,ngenedel,genedel,next2)
	implicit none
	include 'gagasizes.for'
	integer i,j,ngenedel,genedel(maxres),blocksize,nres1,nres2,
     $		next1(0:maxres),next2(0:maxres+1)
	integer*2 fragali1(maxres),fragali2(maxres)
c
	integer k,l,j1
c
	if(ldebug)then
		write(*,*) 'This is geneoverlap',i,j
		write(*,500) (next1(k),k=1,nres1)
		write(*,500) (next2(k),k=1,nres2)
	end if
c
c	check i direction
c
	k=next1(max(0,i-blocksize))
	do while(k.le.min(nres1,i+blocksize-1))
		l=fragali1(k)
		if(k-l.ne.i-j) then
			ngenedel=ngenedel+1
			genedel(ngenedel)=k
		end if
		k=next1(k)
	end do
c
c	check j direction
c
	if(j.gt.0) then
c		parallel - parallel
		l=next2(max(0,j-blocksize))
		do while(l.le.min(nres2,j+blocksize-1))
			k=fragali2(l)
			if((k.gt.0).and.(k-l.ne.i-j).and.
     $				((k.ge.i+blocksize).or.(k.le.i-blocksize)))then
				ngenedel=ngenedel+1
				genedel(ngenedel)=k
			end if
			l=next2(l)
		end do
c		parallel - antiparallel: always clash
		l=next2(max(0,j-1))
		do while(l.le.min(nres2,j+blocksize+blocksize-1))
			k=fragali2(l)
			if((k.lt.0).and.
     $				((-k.ge.i+blocksize).or.(-k.le.i-blocksize)))then
				ngenedel=ngenedel+1
				genedel(ngenedel)=-k
			end if
			l=next2(l)
		end do
	else if(j.lt.0) then
		j1=-j
c		antiparallel - parallel: always clash
		l=next2(max(0,j1-blocksize-blocksize+1))
		do while(l.le.j1)
			k=fragali2(l)
			if((k.gt.0).and.
     $				((k.ge.i+blocksize).or.(k.le.i-blocksize)))then
				ngenedel=ngenedel+1
				genedel(ngenedel)=k
			end if
			l=next2(l)
		end do
c		antiparallel - antiparallel
		l=next2(max(0,j1-blocksize))
		do while(l.le.min(nres2,j1+blocksize-1))
			k=fragali2(l)
			if((k.lt.0).and.(k-l.ne.i-j).and.
     $				((-k.ge.i+blocksize).or.(-k.le.i-blocksize)))then
				ngenedel=ngenedel+1
				genedel(ngenedel)=-k
			end if
			l=next2(l)
		end do
	else
		stop 'geneoverlap'
	end if
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	function getrescore(i1,i2,d1,d2,ali1,nextres,nres1)
c
c	calculates rescore(i1,i2) or dscore(i1,i2) on the fly
c
	implicit none
	include 'gagasizes.for'
	integer i1,i2,nextres(0:maxres+1),nres1
	integer*2 ali1(maxres)
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	real getrescore
c
	real scorefun
c
	real dx,s
	integer x,k
        integer*2 zero
c
        zero=0
	dx=scorefun(zero,zero)
	k=0
	do while(nextres(k).le.nres1)
		k=nextres(k)
c	do k=1,nres1
c	  if(ali1(k).ne.0) then
		x=ali1(k)
		if((k.ne.i1).and.(x.ne.i2)) then
			s=scorefun(d1(i1,k),d2(i2,x))
			dx=dx+s+s
		end if
c	  end if
	end do
	getrescore=dx
	if(ldebug) write(*,*) 'This was getrescore',i1,i2,dx
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine checktopo(i,j,blocksize,prev1,next1,fragali1,ngenedel,
     $		genedel,nres1)
c
c	topology check in range <i-blocksize+1, > i+blocksize-1
c				<j-blocksize+1, > j+blocksize-1
c
	implicit none
	include 'gagasizes.for'
	integer i,j,blocksize,ngenedel,genedel(maxres)
	integer prev1(0:maxres+1),next1(0:maxres+1),nres1
	integer*2 fragali1(maxres)
c
	integer k,l
c
	if(ldebug)then
		write(*,*) 'This is checktopo'
		write(*,500) ngenedel,(fragali1(k),k=1,nres1)
		write(*,500) (next1(k),k=0,nres1+1)
		write(*,500) (prev1(k),k=0,nres1+1)
	end if
	k=min(nres1+1,i+blocksize-1)
	do while(next1(k).le.nres1)
		k=next1(k)
		l=fragali1(k)
		if(l.le.j-blocksize) then
			ngenedel=ngenedel+1
			genedel(ngenedel)=k
		else
			goto 19
		end if
	end do
19	k=max(0,i-blocksize+1)
	do while(prev1(k).gt.0)
		k=prev1(k)
		l=fragali1(k)
		if(l.ge.j+blocksize) then
			ngenedel=ngenedel+1
			genedel(ngenedel)=k
		else
			goto 29
		end if
	end do
29	continue
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine testdeletion(ngenedel,genedel,
     $		blocksize,cnt,rescore,ali1,nextres,nres1,
     $		d1,d2,lsave,dx,nrem,rem)
c
c	in:
c		ngenedel,genedel -- list of genes to be deleted
c	out:
c		dx -- change in score on deletion
c	    	nrem,rem -- list of residues to delete
c
	implicit none
	include 'gagasizes.for'
	integer blocksize,cnt(maxres),nextres(0:maxres+1),
     $		nrem,rem(2,maxres),nres1,ngenedel,genedel(maxres)
	integer*2 ali1(maxres)
	real rescore(maxres,maxres),dx
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	logical lsave
c
	real scorefun
	real getrescore
c
	integer i,j,i1,i2,j1,j2,l,k,ix
	real s
c
c	write(*,*) 'This is testdeletion',i,j
c
c	which residues go with gene (i,j) ?
c
	do ix=1,ngenedel
	  i=genedel(ix)
	  j=ali1(i)
	  call gowithgenes(i,j,1,genedel,ali1,cnt,blocksize,nrem,rem,.false.)
c
c	  calculate dx
c
	  dx=0.0
	  do l=1,nrem
		i1=rem(1,l)
		i2=rem(2,l)
c
c		subtract rescore
c
		if(lsave) then
			dx=dx-rescore(i1,i2)
		else
			dx=dx-getrescore(i1,i2,d1,d2,ali1,nextres,nres1)
		end if
c
c		add rem x rem
c
		do k=1,l-1
			j1=rem(1,k)
			j2=rem(2,k)
			s=scorefun(d1(i1,j1),d2(i2,j2))
			dx=dx+s+s
		end do
	  end do
	end do
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine gowithgenes(i1,j1,ngenedel,genedel,ali1,cnt,blocksize,
     $		nrem,rem,laddition)
c
c	in:	ngenedel, genedel -- list of genes to delete
c		laddition -- TRUE if adding, FALSE if deleting genes
c		i1,j1 -- gene to append (don't remove identical residue pairs)
c	out:	nrem rem -- list of residues to delete (x,y always positive)
c	N.B.:	nrem is initialized to zero !
c
	implicit none
	include 'gagasizes.for'
	integer blocksize,cnt(maxres),nrem,rem(2,maxres),i1,j1
	integer ngenedel,genedel(maxres)
	integer*2 ali1(maxres)
	logical laddition
c
	integer tmp(maxres),l,k,i,m
c
c	copy cnt() to tmp()
c
	do i=1,ngenedel
		do l=0,blocksize-1
			k=abs(genedel(i)+l)
			tmp(k)=cnt(k)
		end do
	end do
c
c	subtract counts from tmp(), remove residue if empty
c
	nrem=0
	do i=1,ngenedel
		do l=0,blocksize-1
			k=abs(genedel(i)+l)
			tmp(k)=tmp(k)-1
			if(tmp(k).eq.0) then
			  if(laddition) then
c				testaddition
				m=k-i1
				if((m.ge.0).and.(m.le.blocksize-1).and.
     $					(ali1(k).eq.abs(j1+m))) goto 19
				nrem=nrem+1
				rem(1,nrem)=k
				rem(2,nrem)=ali1(k)
			  else
c				testdeletion
				nrem=nrem+1
				rem(1,nrem)=k
				rem(2,nrem)=ali1(k)
			  end if
19			end if
		end do
	end do
c
	if(ldebug) then
		write(*,*) 'This was gowithgenes',ngenedel,nrem
		write(*,500) (genedel(i),i=1,ngenedel)
		write(*,500) (rem(1,i),rem(2,i),i=1,nrem)
	end if
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine appendcand(i,j,blocksize,
     $		lset,ltest,dscore,ncand,candij,
     $		ali1,ali2,nres1,d1,d2,nextres,lsave,
     $		ntest,test,testij)
c
c	in: gene (i,j), where j may be negative
c	out: residue pair candidates (x,y), where y is always positive !
c
	implicit none
	include 'gagasizes.for'
	integer blocksize,nres1,i,j
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	logical lsave
	integer*2 ali1(maxres),ali2(maxres)
	logical lset(maxpair),ltest(maxres,maxres)
	integer ncand, candij(2,maxpair*4)
	integer ntest,testij(2,maxpair),test(maxres,-maxres:maxres),
     $		nextres(0:maxres+1)
	real dscore(maxres,maxres)
c
	real addscore
c
	integer l,j1
c
	if(ldebug) write(*,*) 'This is appendcand',i,j,ncand,lsave,ntest
c
c	add to gene-candidate list
c
	if(test(i,j).ne.0) then
		if(lverb)write(*,*) 'skip ',i,j
		return
	end if
	ntest=ntest+1
	testij(1,ntest)=i
	testij(2,ntest)=j
	test(i,j)=ntest
	lset(ntest)=.false.
	if(lsave) then
c
c	  add to residue-candidate lists
c
	  do l=0,blocksize-1
	   j1=abs(j+l)
	   if(.not.ltest(i+l,j1)) then
		ncand=ncand+1
		candij(1,ncand)=i+l
		candij(2,ncand)=j1
		ltest(i+l,j1)=.true.
c
c		calculate dscore(i+l,j1)
c
		if(ali1(i+l).ne.j1) dscore(i+l,j1)
     $			  =addscore(i+l,j1,ali1,ali2,nextres,nres1,d1,d2)
	   end if
	  end do
	end if
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	function addscore(i,j,ali1,ali2,nextres,nres1,d1,d2)
c
c	returns the marginal score of i,j relative to ali1, excluding iold,jold
c	j must be positive
c
	implicit none
	include 'gagasizes.for'
	real addscore
	integer i,j,nres1
	integer*2 ali1(maxres),ali2(maxres)
	integer nextres(0:maxres+1)
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
c
	real scorefun
c
	integer*2 a,b,zero
	real s,dx
	integer k,iold,jold
c
	if(j.lt.0) write(*,*) 'This is addscore',i,j
c
        zero=0
	dx=scorefun(zero,zero)
	iold=ali2(j)
	jold=ali1(i)
	k=0
	do while(nextres(k).le.nres1)
		k=nextres(k)
		if((k.ne.iold).and.(ali1(k).ne.jold)) then
			a=d1(i,k)
			b=d2(j,ali1(k))
			s=scorefun(a,b)
			dx=dx+s+s
			if(ldebug)write(*,*) 'pair',i,j,k,ali1(k),s+s,dx
		end if
	end do
	addscore=dx

	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine pullgene(gblocksize,blocksize,ngene,genepool,
     $		lset,ltest,dscore,ncand,candij,fragali1,next1,
     $		ali1,ali2,nres1,nres2,d1,d2,nextres,lsave,
     $		ntest,test,testij)
	implicit none
	include 'gagasizes.for'
	integer ngene,blocksize,gblocksize,nres1,nres2
	integer*2 d1(maxres,maxres),d2(maxres,maxres)
	logical lsave
	integer*2 genepool(5,maxpair)
	integer*2 ali1(maxres),ali2(maxres),fragali1(maxres)
	logical lset(maxpair),ltest(maxres,maxres)
	integer ncand, candij(2,maxpair*4)
	real dscore(maxres,maxres)
	integer ntest,testij(2,maxpair),test(maxres,-maxres:maxres),
     $		nextres(0:maxres+1),next1(0:maxres+1)
c
	logical lhalfin
c
	integer igene,i,j,l,p,q
c
	if(ldebug)then
		write(*,*) 'This is pullgene',ncand,ntest,nres1,nres2,ngene
		write(*,500) (ali1(i),i=1,nres1)
		write(*,500) (nextres(i),i=1,nres1)
	end if
c
	do igene=1,ngene
		do l=0,gblocksize-blocksize
			i=genepool(1,igene)+l
			j=genepool(3,igene)+l
			p=genepool(2,igene)+l
			q=genepool(4,igene)+l
c
c			gene may be antiparallel, residue pairs may not !
c
c			append one-half to cand-list if other-half is
c				at least partly aligned in ali1
c
			if(test(p,q).eq.0) then
			  if(lhalfin(i,j,fragali1,blocksize,next1,nres1))
     $				call appendcand(p,q,blocksize,
     $				lset,ltest,dscore,ncand,candij,ali1,ali2,
     $				nres1,d1,d2,nextres,lsave,
     $				ntest,test,testij)
			end if
			if(test(i,j).eq.0) then
			  if(lhalfin(p,q,fragali1,blocksize,next1,nres1))
     $				call appendcand(i,j,blocksize,
     $				lset,ltest,dscore,ncand,candij,ali1,ali2,
     $				nres1,d1,d2,nextres,lsave,
     $				ntest,test,testij)
			end if
		end do
	end do
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	function lhalfin(i,j,fragali1,blocksize,next1,nres1)
	implicit none
	include 'gagasizes.for'
	integer i,j,blocksize,next1(0:maxres+1),nres1
	integer*2 fragali1(maxres)
	logical lhalfin
c
	integer k,l
c
	lhalfin=.false.
	k=next1(max(0,i-blocksize))
	do while((.not.lhalfin).and.(k.le.i+blocksize-1).and.(k.le.nres1))
		l=fragali1(k)
		lhalfin=(k-l.eq.i-j)
		k=next1(k)
	end do
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine printcandlist(ncand,candij,dscore,rescore)
	implicit none
	include 'gagasizes.for'
	integer ncand,candij(2,maxpair*4)
	real dscore(maxres,maxres),rescore(maxres,maxres)
c
	integer i,p,q
c
c	print candlist !
c
	write(*,*) 'This is the candidate list',ncand
	do i=1,ncand
		p=candij(1,i)
		q=candij(2,i)
		write(*,500) i,p,q,dscore(p,q),rescore(p,q)
	end do
c
500	format('candidate',3i5,2f10.3,2x,l5)
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine initialize(cnt,ali1,ali2,fragali1,fragali2,
     $		prev1,next1,ltest,ncand,rescore,dscore,
     $		nres1,nres2,ntest,nextres,test,next2)
	implicit none
	include 'gagasizes.for'
	integer cnt(maxres),prev1(0:maxres+1),next1(0:maxres+1)
	integer nres1,nres2
	integer*2 ali1(maxres),ali2(maxres),fragali1(maxres),fragali2(maxres)
	logical ltest(maxres,maxres)
	integer ncand, ntest,nextres(0:maxres+1),next2(0:maxres+1)
	real rescore(maxres,maxres),dscore(maxres,maxres)
	integer test(maxres,-maxres:maxres)
c
	integer i,j
c
	call initali1(ali1,nres1)
	call initali1(ali2,nres2)
	call initali1(fragali1,nres1)
	call initali1(fragali2,nres2)
	do i=1,nres1
		cnt(i)=0
		do j=1,nres2
			ltest(i,j)=.false.
			rescore(i,j)=0.0
			dscore(i,j)=0.0
		end do
	end do
	do i=1,nres1
		do j=-nres2,nres2
			test(i,j)=0
		end do
	end do
	do i=0,nres1+1
		prev1(i)=0
		next1(i)=nres1+1
		nextres(i)=nres1+1
	end do
	do i=0,nres2+1
		next2(i)=nres2+1
	end do
	ncand=0
	ntest=0
c
	return
	end
c
c-------------------------------------------------------------------------------
c
	subroutine checkcnt(cnt,fragali1,nres1,blocksize)
	implicit none
	include 'gagasizes.for'
	integer cnt(maxres),nres1,blocksize
	integer*2 fragali1(maxres)
c
	integer i,l,tmp(maxres)
	logical lp
c
	lp=.false.
	do i=1,nres1
		tmp(i)=0
	end do
	do i=1,nres1
		if(fragali1(i).ne.0) then
			do l=0,blocksize-1
				tmp(i+l)=tmp(i+l)+1
			end do
		end if
	end do
	do i=1,nres1
		if(tmp(i).ne.cnt(i)) then
			write(*,*) 'cnt error',i,cnt(i),tmp(i)
			lp=.true.
		end if
	end do
	if(lp) then
		write(*,500) (fragali1(i),i=1,nres1)
		write(*,500) (cnt(i),i=1,nres1)
	end if
c
500	format(20i4)
c
	return
	end
c
c-------------------------------------------------------------------------------
c

max_stars_repo_name: gkrash/lcspromatch

max_issues_repo_name: gkrash/lcspromatch

max_forks_repo_name: gkrash/lcspromatch


text:
    module solver_preparation

    use custom_precision

    implicit none

    private

    integer                                     :: irow, icol, i

    interface serialise_grid_nodes
    procedure serialise_grid_nodes_sp
    procedure serialise_grid_nodes_dp
    end interface

    interface map_vector_components_to_grid
    procedure map_vector_components_to_grid_sp
    procedure map_vector_components_to_grid_dp
    end interface

    interface fill_coefficient_matrix_poisson
    procedure fill_coefficient_matrix_poisson_sp
    procedure fill_coefficient_matrix_poisson_dp
    end interface

    interface fill_coefficient_matrix_poisson_dirichlet_left
    procedure fill_coefficient_matrix_poisson_dirichlet_left_sp
    procedure fill_coefficient_matrix_poisson_dirichlet_left_dp
    end interface

    interface fill_coefficient_matrix_poisson_neumann_right
    procedure fill_coefficient_matrix_poisson_neumann_right_sp
    procedure fill_coefficient_matrix_poisson_neumann_right_dp
    end interface

    interface fill_coefficient_matrix_poisson_dirichlet_bottom
    procedure fill_coefficient_matrix_poisson_dirichlet_bottom_sp
    procedure fill_coefficient_matrix_poisson_dirichlet_bottom_dp
    end interface

    interface fill_coefficient_matrix_poisson_dirichlet_top
    procedure fill_coefficient_matrix_poisson_dirichlet_top_sp
    procedure fill_coefficient_matrix_poisson_dirichlet_top_dp
    end interface

    public get_linear_index
    public serialise_grid_nodes
    public map_vector_components_to_grid
    public fill_coefficient_matrix_poisson
    public fill_coefficient_matrix_poisson_dirichlet_left
    public fill_coefficient_matrix_poisson_neumann_right
    public fill_coefficient_matrix_poisson_dirichlet_bottom
    public fill_coefficient_matrix_poisson_dirichlet_top

    contains

    !#####################################
    pure function get_row_id(i, stride) result(irow)
    integer, intent(in)         :: i, stride
    integer                     :: irow

    irow = mod(i, stride)
    end function get_row_id

    !#####################################
    pure function get_column_id(i, stride) result(icol)
    integer, intent(in)         :: i, stride
    integer                     :: icol

    icol = 1 + i / stride
    end function get_column_id

    !#####################################
    pure function get_linear_index(irow, icol, stride) result(id)
    integer, intent(in)         :: irow, icol, stride
    integer                     :: id

    id = (icol - 1) * stride + irow
    end function get_linear_index

    !#####################################
    subroutine serialise_grid_nodes_sp(grid, vector, nrows, ncols)
    real(SP), allocatable, intent(in)           :: grid(:,:)
    real(SP), allocatable, intent(inout)        :: vector(:)
    integer, intent(in)                         :: nrows, ncols

    i = 1
    do icol = 1, ncols
        do irow = 1, nrows
            vector(i) = grid(irow, icol)
            i = i + 1
        end do
    end do
    end subroutine serialise_grid_nodes_sp

    !#####################################
    subroutine serialise_grid_nodes_dp(grid, vector, nrows, ncols)
    real(DP), allocatable, intent(in)           :: grid(:,:)
    real(DP), allocatable, intent(inout)        :: vector(:)
    integer, intent(in)                         :: nrows, ncols

    i = 1
    do icol = 1, ncols
        do irow = 1, nrows
            vector(i) = grid(irow, icol)
            i = i + 1
        end do
    end do
    end subroutine serialise_grid_nodes_dp

    !#####################################
    subroutine map_vector_components_to_grid_sp(vector, n, grid, nrows, ncols)
    real(SP), allocatable, intent(in)           :: vector(:)
    real(SP), allocatable, intent(inout)        :: grid(:,:)
    integer, intent(in)                         :: n, nrows, ncols

    do i = 1, n
        irow = get_row_id(i, nrows)
        icol = get_column_id(i, nrows)
        grid(irow, icol) = vector(i)
    end do
    end subroutine map_vector_components_to_grid_sp

    !#####################################
    subroutine map_vector_components_to_grid_dp(vector, n, grid, nrows, ncols)
    real(DP), allocatable, intent(in)           :: vector(:)
    real(DP), allocatable, intent(inout)        :: grid(:,:)
    integer, intent(in)                         :: n, nrows, ncols

    do i = 1, n
        irow = get_row_id(i, nrows)
        icol = get_column_id(i, nrows)
        grid(irow, icol) = vector(i)
    end do
    end subroutine map_vector_components_to_grid_dp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_sp(mat, nrows, ncols, dx, dy)
    real(SP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols
    real(SP), intent(in)                        :: dx, dy

    real(SP)                                    :: dx2, dy2
    integer                                     :: k, n

    dx2 = dx * dx
    dy2 = dy * dy
    mat(:,:) = 0.
    n = nrows * ncols

    do icol = 2, ncols - 1
        do irow = 2, nrows - 1
            k = get_linear_index(irow, icol, nrows)
            mat(k, k) = -((2. / dx2) + (2. / dy2))

            mat(k, k - 1) = 1. / dy2
            mat(k, k + 1) = 1. / dy2

            mat(k, k - nrows) = 1. / dx2
            mat(k, k + nrows) = 1. / dx2
        end do
    end do
    end subroutine fill_coefficient_matrix_poisson_sp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dp(mat, nrows, ncols, dx, dy)
    real(DP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols
    real(DP), intent(in)                        :: dx, dy

    real(DP)                                    :: dx2, dy2
    integer                                     :: k, n

    dx2 = dx * dx
    dy2 = dy * dy
    mat(:,:) = 0.d0
    n = nrows * ncols

    do icol = 2, ncols - 1
        do irow = 2, nrows - 1
            k = get_linear_index(irow, icol, nrows)
            mat(k, k) = -((2.d0 / dx2) + (2.d0 / dy2))

            mat(k, k - 1) = 1.d0 / dy2
            mat(k, k + 1) = 1.d0 / dy2

            mat(k, k - nrows) = 1.d0 / dx2
            mat(k, k + nrows) = 1.d0 / dx2
        end do
    end do
    end subroutine fill_coefficient_matrix_poisson_dp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_left_sp(mat, nrows, ncols)
    real(SP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    icol = 1
    do irow = 1, nrows
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_left_sp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_left_dp(mat, nrows, ncols)
    real(DP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    icol = 1
    do irow = 1, nrows
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.d0
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_left_dp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_neumann_right_sp(mat, nrows, ncols, dx, dy)
    real(SP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols
    real(SP), intent(in)                        :: dx, dy

    integer                                     :: k
    real(SP)                                    :: dx2, dy2

    dx2 = dx * dx
    dy2 = dy *dy

    icol = ncols
    do irow = 2, nrows - 1
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = -(7. / (2. * dx2) + (2. / dy2))
        mat(k, k - 1) = 1. / dy2
        mat(k, k + 1) = 1. / dy2
        mat(k, k - nrows) = 4. / dx2
        mat(k, k + nrows) = -1. / (2. * dx2)
    end do
    end subroutine fill_coefficient_matrix_poisson_neumann_right_sp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_neumann_right_dp(mat, nrows, ncols, dx, dy)
    real(DP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols
    real(DP), intent(in)                        :: dx, dy

    integer                                     :: k
    real(DP)                                    :: dx2, dy2

    dx2 = dx * dx
    dy2 = dy *dy

    icol = ncols
    do irow = 2, nrows - 1
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = -(7.d0 / (2.d0 * dx2) + (2.d0 / dy2))
        mat(k, k - 1) = 1.d0 / dy2
        mat(k, k + 1) = 1.d0 / dy2
        mat(k, k - nrows) = 4.d0 / dx2
        mat(k, k + nrows) = -1.d0 / (2.d0 * dx2)
    end do
    end subroutine fill_coefficient_matrix_poisson_neumann_right_dp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_bottom_sp(mat, nrows, ncols)
    real(SP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    irow = 1
    do icol = 1, ncols
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_bottom_sp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_bottom_dp(mat, nrows, ncols)
    real(DP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    irow = 1
    do icol = 1, ncols
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.d0
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_bottom_dp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_top_sp(mat, nrows, ncols)
    real(SP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    irow = nrows
    do icol = 1, ncols
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_top_sp

    !#####################################
    subroutine fill_coefficient_matrix_poisson_dirichlet_top_dp(mat, nrows, ncols)
    real(DP), allocatable, intent(inout)        :: mat(:,:)
    integer, intent(in)                         :: nrows, ncols

    integer                                     :: k

    irow = nrows
    do icol = 1, ncols
        k = get_linear_index(irow, icol, nrows)
        mat(k, k) = 1.d0
    end do
    end subroutine fill_coefficient_matrix_poisson_dirichlet_top_dp

    end module solver_preparation

max_stars_repo_name: aliakatas/LinearSolver

max_issues_repo_name: aliakatas/LinearSolver

max_forks_repo_name: aliakatas/LinearSolver


text:
The original Silo Café & Pub has been replaced by the new Gunrock Pub.

This was (and still is) the only place on campus that serves alcohol every day during the week. The Pub was run by Sodexho, as is the new Gunrock Pub.

Draft Beer:
  Sierra Nevada Pale Ale, Boont Amber Ale, Fat Tire: $4.25 pint, $10.00 pitcher
  Guinness: $4.50 pint, $12.00 pitcher

Bottled Beer:
  Domestic: $3.25
  Imported: $4.25

Wine:
  Glass: $3.75  $4.25
  Bottle: $14.00  $18.00

Food:
  Entrees: $7.25  $10.00

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE PHO_SASVMD(Iset,Kf,X,Q2,P2,Alam,Xpga,Vxpga)
      IMPLICIT NONE
      REAL aem , aem2pi , Alam , alam3 , alam5 , P2 , p2div , p2eff ,
     &     pmb , pmc , Q2 , q2div , q2eff , s , s2 , s3 , s4 , sbt ,
     &     sch , sll
      REAL Vxpga , X , x1 , xbot , xchm , xglu , xl , Xpga , xsea ,
     &     xsea0 , xval
      INTEGER Iset , Kf , kfa , kfl , nfp , nfq
C...Purpose: to evaluate the VMD parton distributions of a photon,
C...evolved homogeneously from an initial scale P2 to Q2.
C...Does not include dipole suppression factor.
C...ISET is parton distribution set, see above;
C...additionally ISET=0 is used for the evolution of an anomalous photon
C...which branched at a scale P2 and then evolved homogeneously to Q2.
C...ALAM is the 4-flavour Lambda, which is automatically converted
C...to 3- and 5-flavour equivalents as needed.
      SAVE
      DIMENSION Xpga(-6:6) , Vxpga(-6:6)

C  input/output channels
      INCLUDE 'inc/poinou'

      DATA pmc/1.3/ , pmb/4.6/ , aem/0.007297/ , aem2pi/0.0011614/

C...Reset output.
      DO kfl = -6 , 6
         Xpga(kfl) = 0.
         Vxpga(kfl) = 0.
      END DO
      kfa = ABS(Kf)

C...Calculate Lambda; protect against unphysical Q2 and P2 input.
      alam3 = Alam*(pmc/Alam)**(2./27.)
      alam5 = Alam*(Alam/pmb)**(2./23.)
      p2eff = MAX(P2,1.2*alam3**2)
      IF ( kfa.EQ.4 ) p2eff = MAX(p2eff,pmc**2)
      IF ( kfa.EQ.5 ) p2eff = MAX(p2eff,pmb**2)
      q2eff = MAX(Q2,p2eff)

C...Find number of flavours at lower and upper scale.
      nfp = 4
      IF ( p2eff.LT.pmc**2 ) nfp = 3
      IF ( p2eff.GT.pmb**2 ) nfp = 5
      nfq = 4
      IF ( q2eff.LT.pmc**2 ) nfq = 3
      IF ( q2eff.GT.pmb**2 ) nfq = 5

C...Find s as sum of 3-, 4- and 5-flavour parts.
      s = 0.
      IF ( nfp.EQ.3 ) THEN
         q2div = pmc**2
         IF ( nfq.EQ.3 ) q2div = q2eff
         s = s + (6./27.)*LOG(LOG(q2div/alam3**2)/LOG(p2eff/alam3**2))
      END IF
      IF ( nfp.LE.4 .AND. nfq.GE.4 ) THEN
         p2div = p2eff
         IF ( nfp.EQ.3 ) p2div = pmc**2
         q2div = q2eff
         IF ( nfq.EQ.5 ) q2div = pmb**2
         s = s + (6./25.)*LOG(LOG(q2div/Alam**2)/LOG(p2div/Alam**2))
      END IF
      IF ( nfq.EQ.5 ) THEN
         p2div = pmb**2
         IF ( nfp.EQ.5 ) p2div = p2eff
         s = s + (6./23.)*LOG(LOG(q2eff/alam5**2)/LOG(p2div/alam5**2))
      END IF

C...Calculate frequent combinations of x and s.
      x1 = 1. - X
      xl = -LOG(X)
      s2 = s**2
      s3 = s**3
      s4 = s**4

C...Evaluate homogeneous anomalous parton distributions below or
C...above threshold.
      IF ( Iset.EQ.0 ) THEN
         IF ( Q2.LE.P2 .OR. (kfa.EQ.4 .AND. Q2.LT.pmc**2) .OR.
     &        (kfa.EQ.5 .AND. Q2.LT.pmb**2) ) THEN
            xval = X*1.5*(X**2+x1**2)
            xglu = 0.
            xsea = 0.
         ELSE
            xval = (1.5/(1.-0.197*s+4.33*s2)*X**2+(1.5+2.10*s)
     &             /(1.+3.29*s)*x1**2+5.23*s/(1.+1.17*s+19.9*s3)*X*x1)
     &             *X**(1./(1.+1.5*s))*(1.-X**2)**(2.667*s)
            xglu = 4.*s/(1.+4.76*s+15.2*s2+29.3*s4)
     &             *X**(-2.03*s/(1.+2.44*s))*(x1*xl)**(1.333*s)
     &             *((4.*X**2+7.*X+4.)*x1/3.-2.*X*(1.+X)*xl)
            xsea = s2/(1.+4.54*s+8.19*s2+8.05*s3)
     &             *X**(-1.54*s/(1.+1.29*s))*x1**(2.667*s)
     &             *((8.-73.*X+62.*X**2)*x1/9.+(3.-8.*X**2/3.)
     &             *X*xl+(2.*X-1.)*X*xl**2)
         END IF

C...Evaluate set 1D parton distributions below or above threshold.
      ELSE IF ( Iset.EQ.1 ) THEN
         IF ( Q2.LE.P2 .OR. (kfa.EQ.4 .AND. Q2.LT.pmc**2) .OR.
     &        (kfa.EQ.5 .AND. Q2.LT.pmb**2) ) THEN
            xval = 1.294*X**0.80*x1**0.76
            xglu = 1.273*X**0.40*x1**1.76
            xsea = 0.100*x1**3.76
         ELSE
            xval = 1.294/(1.+0.252*s+3.079*s2)*X**(0.80-0.13*s)
     &             *x1**(0.76+0.667*s)*xl**(2.*s)
            xglu = 7.90*s/(1.+5.50*s)*EXP(-5.16*s)
     &             *X**(-1.90*s/(1.+3.60*s))*x1**1.30*xl**(0.50+3.*s)
     &             + 1.273*EXP(-10.*s)*X**0.40*x1**(1.76+3.*s)
            xsea = (0.1-0.397*s2+1.121*s3)/(1.+5.61*s2+5.26*s3)
     &             *X**(-7.32*s2/(1.+10.3*s2))
     &             *x1**((3.76+15.*s+12.*s2)/(1.+4.*s))
            xsea0 = 0.100*x1**3.76
         END IF

C...Evaluate set 1M parton distributions below or above threshold.
      ELSE IF ( Iset.EQ.2 ) THEN
         IF ( Q2.LE.P2 .OR. (kfa.EQ.4 .AND. Q2.LT.pmc**2) .OR.
     &        (kfa.EQ.5 .AND. Q2.LT.pmb**2) ) THEN
            xval = 0.8477*X**0.51*x1**1.37
            xglu = 3.42*X**0.255*x1**2.37
            xsea = 0.
         ELSE
            xval = 0.8477/(1.+1.37*s+2.18*s2+3.73*s3)*X**(0.51+0.21*s)
     &             *x1**1.37*xl**(2.667*s)
            xglu = 24.*s/(1.+9.6*s+0.92*s2+14.34*s3)*EXP(-5.94*s)
     &             *X**((-0.013-1.80*s)/(1.+3.14*s))*x1**(2.37+0.4*s)
     &             *xl**(0.32+3.6*s) + 3.42*EXP(-12.*s)
     &             *X**0.255*x1**(2.37+3.*s)
            xsea = 0.842*s/(1.+21.3*s-33.2*s2+229.*s3)
     &             *X**((0.13-2.90*s)/(1.+5.44*s))*x1**(3.45+0.5*s)
     &             *xl**(2.8*s)
            xsea0 = 0.
         END IF

C...Evaluate set 2D parton distributions below or above threshold.
      ELSE IF ( Iset.EQ.3 ) THEN
         IF ( Q2.LE.P2 .OR. (kfa.EQ.4 .AND. Q2.LT.pmc**2) .OR.
     &        (kfa.EQ.5 .AND. Q2.LT.pmb**2) ) THEN
            xval = X**0.46*x1**0.64 + 0.76*X
            xglu = 1.925*x1**2
            xsea = 0.242*x1**4
         ELSE
            xval = (1.+0.186*s)/(1.-0.209*s+1.495*s2)*X**(0.46+0.25*s)
     &             *x1**((0.64+0.14*s+5.*s2)/(1.+s))*xl**(1.9*s)
     &             + (0.76+0.4*s)*X*x1**(2.667*s)
            xglu = (1.925+5.55*s+147.*s2)/(1.-3.59*s+3.32*s2)
     &             *EXP(-18.67*s)
     &             *X**((-5.81*s-5.34*s2)/(1.+29.*s-4.26*s2))
     &             *x1**((2.-5.9*s)/(1.+1.7*s))*xl**(9.3*s/(1.+1.7*s))
            xsea = (0.242-0.252*s+1.19*s2)/(1.-0.607*s+21.95*s2)
     &             *X**(-12.1*s2/(1.+2.62*s+16.7*s2))*x1**4*xl**s
            xsea0 = 0.242*x1**4
         END IF

C...Evaluate set 2M parton distributions below or above threshold.
      ELSE IF ( Iset.EQ.4 ) THEN
         IF ( Q2.LE.P2 .OR. (kfa.EQ.4 .AND. Q2.LT.pmc**2) .OR.
     &        (kfa.EQ.5 .AND. Q2.LT.pmb**2) ) THEN
            xval = 1.168*X**0.50*x1**2.60 + 0.965*X
            xglu = 1.808*x1**2
            xsea = 0.209*x1**4
         ELSE
            xval = (1.168+1.771*s+29.35*s2)*EXP(-5.776*s)
     &             *X**((0.5+0.208*s)/(1.-0.794*s+1.516*s2))
     &             *x1**((2.6+7.6*s)/(1.+5.*s))*xl**(5.15*s/(1.+2.*s))
     &             + (0.965+22.35*s)/(1.+18.4*s)*X*x1**(2.667*s)
            xglu = (1.808+29.9*s)/(1.+26.4*s)*EXP(-5.28*s)
     &             *X**((-5.35*s-10.11*s2)/(1.+31.71*s))
     &             *x1**((2.-7.3*s+4.*s2)/(1.+2.5*s))
     &             *xl**(10.9*s/(1.+2.5*s))
            xsea = (0.209+0.644*s2)/(1.+0.319*s+17.6*s2)
     &             *X**((-0.373*s-7.71*s2)/(1.+0.815*s+11.0*s2))
     &             *x1**(4.+s)*xl**(0.45*s)
            xsea0 = 0.209*x1**4
         END IF
      END IF

C...Threshold factors for c and b sea.
      sll = LOG(LOG(q2eff/Alam**2)/LOG(p2eff/Alam**2))
      xchm = 0.
      IF ( Q2.GT.pmc**2 .AND. Q2.GT.1.001*p2eff ) THEN
         sch = MAX(0.,LOG(LOG(pmc**2/Alam**2)/LOG(p2eff/Alam**2)))
         IF ( Iset.EQ.0 ) THEN
            xchm = xsea*(1.-(sch/sll)**2)
         ELSE
            xchm = MAX(0.,xsea-xsea0*x1**(2.667*s))*(1.-sch/sll)
         END IF
      END IF
      xbot = 0.
      IF ( Q2.GT.pmb**2 .AND. Q2.GT.1.001*p2eff ) THEN
         sbt = MAX(0.,LOG(LOG(pmb**2/Alam**2)/LOG(p2eff/Alam**2)))
         IF ( Iset.EQ.0 ) THEN
            xbot = xsea*(1.-(sbt/sll)**2)
         ELSE
            xbot = MAX(0.,xsea-xsea0*x1**(2.667*s))*(1.-sbt/sll)
         END IF
      END IF

C...Fill parton distributions.
      Xpga(0) = xglu
      Xpga(1) = xsea
      Xpga(2) = xsea
      Xpga(3) = xsea
      Xpga(4) = xchm
      Xpga(5) = xbot
      Xpga(kfa) = Xpga(kfa) + xval
      DO kfl = 1 , 5
         Xpga(-kfl) = Xpga(kfl)
      END DO
      Vxpga(kfa) = xval
      Vxpga(-kfa) = xval

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
! =========================================================================
elemental function gsw_pt0_from_t_ice (t, p)
! =========================================================================
!
!  Calculates potential temperature of ice Ih with a reference pressure of
!  0 dbar, from in-situ temperature, t.
!
!  t   =  in-situ temperature  (ITS-90)                           [ deg C ]
!  p   =  sea pressure                                             [ dbar ]
!         ( i.e. absolute pressure - 10.1325 dbar )
!
!  pt0_ice  =  potential temperature of ice Ih with reference pressure of
!              zero dbar (ITS-90)                                 [ deg C ]
!--------------------------------------------------------------------------

use gsw_mod_teos10_constants, only : gsw_t0

use gsw_mod_toolbox, only : gsw_gibbs_ice_part_t, gsw_gibbs_ice_pt0
use gsw_mod_toolbox, only : gsw_gibbs_ice_pt0_pt0

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: t, p

real (r8) :: gsw_pt0_from_t_ice

integer :: number_of_iterations
real (r8) :: dentropy, dentropy_dt, pt0_ice
real (r8) :: pt0_ice_old, ptm_ice, true_entropy

! This is the starting polynomial for pt0 of ice Ih.
real (r8), parameter :: s1 = -2.256611570832386e-4_r8
real (r8), parameter :: s2 = -6.045305921314694e-7_r8
real (r8), parameter :: s3 =  5.546699019612661e-9_r8
real (r8), parameter :: s4 =  1.795030639186685e-11_r8
real (r8), parameter :: s5 =  1.292346094030742e-9_r8

real (r8), parameter :: p1 = -2.259745637898635e-4_r8
real (r8), parameter :: p2 =  1.486236778150360e-9_r8
real (r8), parameter :: p3 =  6.257869607978536e-12_r8
real (r8), parameter :: p4 = -5.253795281359302e-7_r8
real (r8), parameter :: p5 =  6.752596995671330e-9_r8
real (r8), parameter :: p6 =  2.082992190070936e-11_r8

real (r8), parameter :: q1 = -5.849191185294459e-15_r8
real (r8), parameter :: q2 =  9.330347971181604e-11_r8
real (r8), parameter :: q3 =  3.415888886921213e-13_r8
real (r8), parameter :: q4 =  1.064901553161811e-12_r8
real (r8), parameter :: q5 = -1.454060359158787e-10_r8
real (r8), parameter :: q6 = -5.323461372791532e-13_r8

true_entropy = -gsw_gibbs_ice_part_t(t,p)

if (t.lt.-45.0_r8 .and. t.gt.-273.0_r8) then

    pt0_ice = t + p*(p1 + p*(p2 + p3*t) + t*(p4 + t*(p5 + p6*t)))

    if (pt0_ice.lt.-gsw_t0) pt0_ice = -gsw_t0

    ! we add 0.05d0 to the initial estimate of pt0_ice at
    ! temps less than -273 to ensure that it is never less than -273.15.
    if (pt0_ice.lt.-273.0_r8) pt0_ice = pt0_ice + 0.05_r8

    dentropy_dt = -gsw_gibbs_ice_pt0_pt0(pt0_ice)

    do number_of_iterations = 1, 3
        pt0_ice_old = pt0_ice
        dentropy = -gsw_gibbs_ice_pt0(pt0_ice_old) - true_entropy
        pt0_ice = pt0_ice_old - dentropy/dentropy_dt
        ptm_ice = 0.5_r8*(pt0_ice + pt0_ice_old)
        dentropy_dt = -gsw_gibbs_ice_pt0_pt0(ptm_ice)
        pt0_ice = pt0_ice_old - dentropy/dentropy_dt
    end do

else

    pt0_ice = t + p*(s1 + t*(s2 + t*(s3 + t*s4)) + s5*p)
    dentropy_dt = -gsw_gibbs_ice_pt0_pt0(pt0_ice)

    pt0_ice_old = pt0_ice
    dentropy = -gsw_gibbs_ice_pt0(pt0_ice_old) - true_entropy

    pt0_ice = pt0_ice_old - dentropy/dentropy_dt
    ptm_ice = 0.5_r8*(pt0_ice + pt0_ice_old)
    dentropy_dt = -gsw_gibbs_ice_pt0_pt0(ptm_ice)
    pt0_ice = pt0_ice_old - dentropy/dentropy_dt

end if

if (pt0_ice.lt.-273.0_r8) then

    pt0_ice = t + p*(q1 + p*(q2 + q3*t) + t*(q4 + t*(q5 + q6*t)))

    ! add 0.01d0 to the initial estimate of pt_ice used in the derivative to
    ! ensure that it is never less than -273.15d0 because the derivative
    ! approaches zero at absolute zero.
    dentropy_dt = -gsw_gibbs_ice_pt0_pt0(pt0_ice+0.01_r8)

    do number_of_iterations = 1, 3
        pt0_ice_old = pt0_ice
        dentropy = -gsw_gibbs_ice_pt0(pt0_ice_old) - true_entropy
        pt0_ice = pt0_ice_old - dentropy/dentropy_dt
        ptm_ice = 0.5_r8*(pt0_ice + pt0_ice_old)
        ! add 0.01d0 to the estimate of ptm_ice for temperatures less than
        ! -273 to ensure that they are never less than -273.15d0 because
        ! the derivative approaches zero at absolute zero and the addition
        ! of 0.01d0 degrees c ensures that when we divide by the derivatve
        ! in the modified newton routine the function does not blow up.
        ptm_ice = ptm_ice + 0.01_r8
        dentropy_dt = -gsw_gibbs_ice_pt0_pt0(ptm_ice)
        pt0_ice = pt0_ice_old - dentropy/dentropy_dt
    end do

end if

! For temperatures less than -273.1 degsC the maximum error is less than
! 2x10^-7 degsC. For temperatures between -273.1 and 273 the maximum error
! is less than 8x10^-8 degsC, and for temperatures greater than -273 degsC the
! maximum error is 1.5x10^-12 degsC.   These errors are over the whole
! ocean depths with p varying between 0 and 10,000 dbar, while the in-situ
! temperature varied independently between -273.15 and +2 degsC.

gsw_pt0_from_t_ice = pt0_ice

return
end function

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

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


!> \file mo_functions.f90

!> \brief Special functions such as gamma function.

!> \details This module provides special functions such as the gamma function.

!> \authors Matthias Cuntz
!> \date May 2014
MODULE mo_functions

  ! Provide special functions.

  ! Written  Matthias Cuntz, May 2014
  ! Modified Matthias Cuntz, Dec 2017 - morris
  !          Stephan Thober, Aug 2020 - added beta function

  ! License
  ! -------
  ! This file is part of the JAMS Fortran package, distributed under the MIT License.
  !
  ! Copyright (c) 2014-2017 Matthias Cuntz - mc (at) macu (dot) de
  !
  ! Permission is hereby granted, free of charge, to any person obtaining a copy
  ! of this software and associated documentation files (the "Software"), to deal
  ! in the Software without restriction, including without limitation the rights
  ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  ! copies of the Software, and to permit persons to whom the Software is
  ! furnished to do so, subject to the following conditions:
  !
  ! The above copyright notice and this permission notice shall be included in all
  ! copies or substantial portions of the Software.
  !
  ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  ! SOFTWARE.

  USE mo_kind, ONLY: i4, i8, sp, dp

  IMPLICIT NONE

  PUBLIC :: factln    ! ln(n!)
  PUBLIC :: factorial ! n!
  PUBLIC :: gammln    ! ln(gamma)
  PUBLIC :: gamm      ! gamm
  PUBLIC :: morris    ! elementary effects test function after Morris (1991)
  PUBLIC :: beta      ! beta function

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         factln
  !
  !     PURPOSE
  !         Logarithm of factorial.
  !
  !>        \brief \f$ \ln(n!) \f$.
  !
  !>        \details Elemental function of the logarithm of the factorial:
  !>        \f[ \ln(n!) \f]
  !>        with
  !>        \f[ n! = 1 \dot 2 \dot ... \dot n \f]
  !>
  !>         Uses n! = \f$ \Gamma(n+1) \f$.
  !
  !     INTENT(IN)
  !>        \param[in] "integer(i4/i8) :: n"        input
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(dp) :: factln — \f$ \ln(n!) \f$
  !
  !     RESTRICTIONS
  !>       \note return is always double precision.
  !
  !     EXAMPLE
  !         fac = exp(factln(5))
  !         fac -> 120
  !         -> see also example in test directory
  !
  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date May 2014
  INTERFACE factln
     MODULE PROCEDURE factln_i4, factln_i8
  END INTERFACE factln

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         factorial
  !
  !     PURPOSE
  !         Factorial.
  !
  !>        \brief \f$ n! \f$.
  !
  !>        \details Elemental function of the factorial:
  !>        \f[ n! = 1 \dot 2 \dot ... \dot n \f]
  !>
  !>         Uses n! = \f$ \Gamma(n+1) \f$.
  !
  !     INTENT(IN)
  !>        \param[in] "integer(i4/i8) :: n"        input
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(dp) :: factln — \f$ \ln(n!) \f$
  !
  !     RESTRICTIONS
  !>       \note return is always double precision.
  !
  !     EXAMPLE
  !         fac = factorial(5)
  !         fac -> 120
  !         -> see also example in test directory
  !
  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date Feb 2013
  ! ------------------------------------------------------------------
  INTERFACE factorial
     MODULE PROCEDURE factorial_i4, factorial_i8
  END INTERFACE factorial

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         gamm
  !
  !     PURPOSE
  !         The gamma function.
  !
  !>        \brief \f$ \Gamma(z) \f$.
  !
  !>        \details Elemental function of the gamma function:
  !>        \f[ \Gamma = \Int_0^\Inf t^{z-1} e^{-t} dt \f]
  !>
  !>         Uses Lanczos-type approximation to ln(gamma) for z > 0.
  !>         The function calculates gamma from abs(z).
  !>         Accuracy: About 14 significant digits except for small regions in the vicinity of 1 and 2.
  !
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: z"        input
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: gamm — \f$ \Gamma(z) \f$
  !
  !     RESTRICTIONS
  !>       \note z must be >0. gammln is calculated from abs(z).
  !
  !     EXAMPLE
  !         vec = (/ 1., 2, 3., -9., 5., 6. /)
  !         gl  = gamm(vec)
  !         -> see also example in test directory
  !
  !     LITERATURE
  !         Lanczos, C. ''A precision approximation of the gamma function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
  !
  !     HISTORY
  !>        \author Alan Miller, 1 Creswick Street, Brighton, Vic. 3187, Australia, e-mail: amiller @ bigpond.net.au
  !>        \date Oct 1996
  !         Modified, Matthias Cuntz, May 2014
  INTERFACE gamm
     MODULE PROCEDURE gamm_sp, gamm_dp
  END INTERFACE gamm

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         gammln
  !
  !     PURPOSE
  !         Logarithm of the gamma function.
  !
  !>        \brief \f$ \ln(\Gamma(z)) \f$.
  !
  !>        \details Elemental function of the logarithm of the gamma function:
  !>        \f[ \Gamma = \Int_0^\Inf t^{z-1} e^{-t} dt \f]
  !>
  !>         Uses Lanczos-type approximation to ln(gamma) for z > 0.
  !>         The function calculates ln(gamma) from abs(z).
  !>         Accuracy: About 14 significant digits except for small regions in the vicinity of 1 and 2.
  !
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: z"        input
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: gammln — \f$ \ln\Gamma(z) \f$
  !
  !     RESTRICTIONS
  !>       \note z must be >0. gammln is calculated from abs(z).
  !
  !     EXAMPLE
  !         vec = (/ 1., 2, 3., -9., 5., 6. /)
  !         gl  = gammln(vec)
  !         -> see also example in test directory
  !
  !     LITERATURE
  !         Lanczos, C. ''A precision approximation of the gamma function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
  !
  !     HISTORY
  !>        \author Alan Miller, 1 Creswick Street, Brighton, Vic. 3187, Australia, e-mail: amiller @ bigpond.net.au
  !>        \date Oct 1996
  !         Modified, Matthias Cuntz, May 2014
  INTERFACE gammln
     MODULE PROCEDURE gammln_sp, gammln_dp
  END INTERFACE gammln

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         morris
  !
  !     PURPOSE
  !>        \details 20-dimension test function for elementary effects after Morris (1991).
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp), dimension(20[,npoints]) :: x"        number of points
  !>        \param[in] "real(sp/dp)                          :: beta0"    parameters
  !>        \param[in] "real(sp/dp), dimension(20)           :: beta1"    parameters
  !>        \param[in] "real(sp/dp), dimension(20,20)        :: beta2"    parameters
  !>        \param[in] "real(sp/dp), dimension(20,20,20)     :: beta3"    parameters
  !>        \param[in] "real(sp/dp), dimension(20,20,20,20)  :: beta4"    parameters
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp)[, dimension(npoints)] :: morris — result of Morris function
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         -> see example in test directory
  !
  !     LITERATURE
  !         Morris (1991), Factorial sampling plans for preliminary computational experiments,
  !         Technometrics 33, 161-174.
  !
  !     HISTORY
  !>        \author Matthias Cuntz, Juliane Mai, Mar 2015 in Python
  !>        \date Mar 2015
  !         Modified, Matthias Cuntz, Dec 2017 - ported to Fortran
  INTERFACE morris
     MODULE PROCEDURE morris_0d_dp, morris_1d_dp, &
          morris_0d_sp, morris_1d_sp
  END INTERFACE morris

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         beta
  !
  !     PURPOSE
  !>        \details calculates beta function
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp)                          :: alpha"    parameters
  !>        \param[in] "real(sp/dp)                          :: beta"     parameters
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp)[, dimension(npoints)] :: beta — result of beta function
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         -> see example in test directory
  !
  !     LITERATURE
  !         https://en.wikipedia.org/wiki/Beta_function
  !
  !     HISTORY
  !>        \author Stephan Thober
  !>        \date Aug 2020
  !         Modified
  INTERFACE beta
     MODULE PROCEDURE beta_sp, beta_dp
  END INTERFACE beta

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

  PRIVATE

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

CONTAINS

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

  ELEMENTAL PURE FUNCTION factln_i4(n)

    IMPLICIT NONE

    INTEGER(i4), INTENT(IN) :: n
    REAL(dp)                :: factln_i4

    factln_i4 = gammln(real(n,dp)+1.0_dp)

  end function factln_i4

  ELEMENTAL PURE FUNCTION factln_i8(n)

    IMPLICIT NONE

    INTEGER(i8), INTENT(IN) :: n
    REAL(dp)                :: factln_i8

    factln_i8 = gammln(real(n,dp)+1.0_dp)

  end function factln_i8

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

  elemental pure function factorial_i4(n)
    ! Returns the factorial n!
    IMPLICIT NONE

    INTEGER(i4), INTENT(IN) :: n
    INTEGER(i4)             :: factorial_i4

    factorial_i4 = nint(exp(factln(n)),i4)

  end function factorial_i4

  elemental pure function factorial_i8(n)
    ! Returns the factorial n!
    IMPLICIT NONE

    INTEGER(i8), INTENT(IN) :: n
    INTEGER(i8)             :: factorial_i8

    factorial_i8 = nint(exp(factln(n)),i8)

  end function factorial_i8

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

  elemental pure function gamm_dp(z)

    IMPLICIT NONE

    REAL(dp), INTENT(IN) :: z
    REAL(dp)             :: gamm_dp

    gamm_dp = exp(gammln(z))

  end function gamm_dp

  elemental pure function gamm_sp(z)

    IMPLICIT NONE

    REAL(sp), INTENT(IN) :: z
    REAL(sp)             :: gamm_sp

    gamm_sp = exp(gammln(z))

  end function gamm_sp

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

  ELEMENTAL PURE FUNCTION gammln_dp(z)
    !  Uses Lanczos-type approximation to ln(gamma) for z > 0.
    !  Reference:
    !       Lanczos, C. ''A precision approximation of the gamma
    !               function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
    !  Accuracy: About 14 significant digits except for small regions
    !            in the vicinity of 1 and 2.
    !  Programmer: Alan Miller
    !              1 Creswick Street, Brighton, Vic. 3187, Australia
    !  e-mail: amiller @ bigpond.net.au
    !  Latest revision - 14 October 1996

    IMPLICIT NONE

    REAL(dp), INTENT(IN) :: z
    REAL(dp)             :: gammln_dp

    ! Local variables
    REAL(dp), parameter :: a(9) = (/ 0.9999999999995183_dp, 676.5203681218835_dp, &
         -1259.139216722289_dp, 771.3234287757674_dp, &
         -176.6150291498386_dp, 12.50734324009056_dp, &
         -0.1385710331296526_dp, 0.9934937113930748E-05_dp, &
         0.1659470187408462E-06_dp /)
    REAL(dp), parameter :: zero = 0.0_dp,   &
         one = 1.0_dp, &
         lnsqrt2pi =  0.9189385332046727_dp, &
         half = 0.5_dp, &
         sixpt5 = 6.5_dp, &
         seven = 7.0_dp
    REAL(dp)  :: tmp
    INTEGER   :: j

    ! if (z <= 0.0_dp) stop 'Error gammln_dp: z <= 0'
    gammln_dp = zero
    tmp = abs(z) + seven
    DO j = 9, 2, -1
       gammln_dp = gammln_dp + a(j)/tmp
       tmp = tmp - one
    END DO
    gammln_dp = gammln_dp + a(1)
    gammln_dp = LOG(gammln_dp) + lnsqrt2pi - (abs(z) + sixpt5) + (abs(z) - half)*LOG(abs(z) + sixpt5)

  END FUNCTION gammln_dp

  ELEMENTAL PURE FUNCTION gammln_sp(z)
    !  Uses Lanczos-type approximation to ln(gamma) for z > 0.
    !  Reference:
    !       Lanczos, C. ''A precision approximation of the gamma
    !               function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
    !  Accuracy: About 14 significant digits except for small regions
    !            in the vicinity of 1 and 2.
    !  Programmer: Alan Miller
    !              1 Creswick Street, Brighton, Vic. 3187, Australia
    !  e-mail: amiller @ bigpond.net.au
    !  Latest revision - 14 October 1996

    IMPLICIT NONE

    REAL(sp), INTENT(IN) :: z
    REAL(sp)             :: gammln_sp

    ! Local variables
    REAL(sp), parameter :: a(9) = (/ 0.9999999999995183_sp, 676.5203681218835_sp, &
         -1259.139216722289_sp, 771.3234287757674_sp, &
         -176.6150291498386_sp, 12.50734324009056_sp, &
         -0.1385710331296526_sp, 0.9934937113930748E-05_sp, &
         0.1659470187408462E-06_sp /)
    REAL(sp), parameter :: zero = 0.0_sp,   &
         one = 1.0_sp, &
         lnsqrt2pi =  0.9189385332046727_sp, &
         half = 0.5_sp, &
         sixpt5 = 6.5_sp, &
         seven = 7.0_sp
    REAL(sp)  :: tmp
    INTEGER   :: j

    ! if (z <= 0.0_sp) stop 'Error gammln_sp: z <= 0'
    gammln_sp = zero
    tmp = abs(z) + seven
    DO j = 9, 2, -1
       gammln_sp = gammln_sp + a(j)/tmp
       tmp = tmp - one
    END DO
    gammln_sp = gammln_sp + a(1)
    gammln_sp = LOG(gammln_sp) + lnsqrt2pi - (abs(z) + sixpt5) + (abs(z) - half)*LOG(abs(z) + sixpt5)

  END FUNCTION gammln_sp

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

  FUNCTION morris_0d_dp(x, beta0, beta1, beta2, beta3, beta4)

    use mo_kind, only: dp, i4

    implicit none

    real(dp), dimension(20),          intent(in) :: x
    real(dp),                         intent(in) :: beta0
    real(dp), dimension(20),          intent(in) :: beta1
    real(dp), dimension(20,20),       intent(in) :: beta2
    real(dp), dimension(20,20,20),    intent(in) :: beta3
    real(dp), dimension(20,20,20,20), intent(in) :: beta4
    real(dp)                                     :: morris_0d_dp

    ! local
    real(dp),    dimension(20) :: om
    integer(i4), dimension(3) :: ii
    integer(i4) :: i, j, l, s, nn

    om = 2._dp*(x - 0.5_dp)
    ii = (/2, 4, 6/) + 1 ! +1 Python -> Fortran
    om(ii) = 2._dp*(1.1_dp*x(ii)/(x(ii)+0.1_dp) - 0.5_dp)

    nn = size(x)

    morris_0d_dp = beta0
    do i=1, nn
       morris_0d_dp = morris_0d_dp + beta1(i)*om(i)
       do j=i+1, nn
          morris_0d_dp = morris_0d_dp + beta2(i,j)*om(i)*om(j)
          do l=j+1, nn
             morris_0d_dp = morris_0d_dp + beta3(i,j,l)*om(i)*om(j)*om(l)
             do s=l+1, nn
                morris_0d_dp = morris_0d_dp + beta4(i,j,l,s)*om(i)*om(j)*om(l)*om(s)
             end do
          end do
       end do
    end do

  END FUNCTION morris_0d_dp

  FUNCTION morris_1d_dp(x, beta0, beta1, beta2, beta3, beta4)

    use mo_kind, only: dp, i4

    implicit none

    real(dp), dimension(:,:),         intent(in) :: x
    real(dp),                         intent(in) :: beta0
    real(dp), dimension(20),          intent(in) :: beta1
    real(dp), dimension(20,20),       intent(in) :: beta2
    real(dp), dimension(20,20,20),    intent(in) :: beta3
    real(dp), dimension(20,20,20,20), intent(in) :: beta4
    real(dp), dimension(size(x,2))               :: morris_1d_dp

    ! local
    real(dp),    dimension(20,size(x,2)) :: om
    integer(i4), dimension(3)            :: ii
    integer(i4) :: i, j, l, s, nn

    om = 2._dp*(x - 0.5_dp)
    ii = (/2, 4, 6/)
    om(ii,:) = 2._dp*(1.1_dp*x(ii,:)/(x(ii,:)+0.1_dp) - 0.5_dp)

    nn = size(x)

    morris_1d_dp(:) = beta0
    do i=1, nn
       morris_1d_dp(:) = morris_1d_dp(:) + beta1(i)*om(i,:)
       do j=i+1, nn
          morris_1d_dp(:) = morris_1d_dp(:) + beta2(i,j)*om(i,:)*om(j,:)
          do l=j+1, nn
             morris_1d_dp(:) = morris_1d_dp(:) + beta3(i,j,l)*om(i,:)*om(j,:)*om(l,:)
             do s=l+1, nn
                morris_1d_dp(:) = morris_1d_dp(:) + beta4(i,j,l,s)*om(i,:)*om(j,:)*om(l,:)*om(s,:)
             end do
          end do
       end do
    end do

  END FUNCTION morris_1d_dp

  FUNCTION morris_0d_sp(x, beta0, beta1, beta2, beta3, beta4)

    use mo_kind, only: sp, i4

    implicit none

    real(sp), dimension(20),          intent(in) :: x
    real(sp),                         intent(in) :: beta0
    real(sp), dimension(20),          intent(in) :: beta1
    real(sp), dimension(20,20),       intent(in) :: beta2
    real(sp), dimension(20,20,20),    intent(in) :: beta3
    real(sp), dimension(20,20,20,20), intent(in) :: beta4
    real(sp)                                     :: morris_0d_sp

    ! local
    real(sp),    dimension(20) :: om
    integer(i4), dimension(3) :: ii
    integer(i4) :: i, j, l, s, nn

    om = 2._sp*(x - 0.5_sp)
    ii = (/2, 4, 6/)
    om(ii) = 2._sp*(1.1_sp*x(ii)/(x(ii)+0.1_sp) - 0.5_sp)

    nn = size(x)

    morris_0d_sp = beta0
    do i=1, nn
       morris_0d_sp = morris_0d_sp + beta1(i)*om(i)
       do j=i+1, nn
          morris_0d_sp = morris_0d_sp + beta2(i,j)*om(i)*om(j)
          do l=j+1, nn
             morris_0d_sp = morris_0d_sp + beta3(i,j,l)*om(i)*om(j)*om(l)
             do s=l+1, nn
                morris_0d_sp = morris_0d_sp + beta4(i,j,l,s)*om(i)*om(j)*om(l)*om(s)
             end do
          end do
       end do
    end do

  END FUNCTION morris_0d_sp

  FUNCTION morris_1d_sp(x, beta0, beta1, beta2, beta3, beta4)

    use mo_kind, only: sp, i4

    implicit none

    real(sp), dimension(:,:),         intent(in) :: x
    real(sp),                         intent(in) :: beta0
    real(sp), dimension(20),          intent(in) :: beta1
    real(sp), dimension(20,20),       intent(in) :: beta2
    real(sp), dimension(20,20,20),    intent(in) :: beta3
    real(sp), dimension(20,20,20,20), intent(in) :: beta4
    real(sp), dimension(size(x,2))               :: morris_1d_sp

    ! local
    real(sp),    dimension(20,size(x,2)) :: om
    integer(i4), dimension(3)            :: ii
    integer(i4) :: i, j, l, s, nn

    om = 2._sp*(x - 0.5_sp)
    ii = (/2, 4, 6/)
    om(ii,:) = 2._sp*(1.1_sp*x(ii,:)/(x(ii,:)+0.1_sp) - 0.5_sp)

    nn = size(x)

    morris_1d_sp(:) = beta0
    do i=1, nn
       morris_1d_sp(:) = morris_1d_sp(:) + beta1(i)*om(i,:)
       do j=i+1, nn
          morris_1d_sp(:) = morris_1d_sp(:) + beta2(i,j)*om(i,:)*om(j,:)
          do l=j+1, nn
             morris_1d_sp(:) = morris_1d_sp(:) + beta3(i,j,l)*om(i,:)*om(j,:)*om(l,:)
             do s=l+1, nn
                morris_1d_sp(:) = morris_1d_sp(:) + beta4(i,j,l,s)*om(i,:)*om(j,:)*om(l,:)*om(s,:)
             end do
          end do
       end do
    end do

  END FUNCTION morris_1d_sp

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

  elemental pure function beta_sp(alpha, beta)

    use mo_kind, only: sp

    real(sp), intent(in) :: alpha
    real(sp), intent(in) :: beta
    real(sp)             :: beta_sp

    beta_sp = (gamm(alpha) * gamm(beta)) / gamm(alpha + beta)

  end function beta_sp

  elemental pure function beta_dp(alpha, beta)

    use mo_kind, only: dp

    real(dp), intent(in) :: alpha
    real(dp), intent(in) :: beta
    real(dp)             :: beta_dp

    beta_dp = (gamm(alpha) * gamm(beta)) / gamm(alpha + beta)

  end function beta_dp

END MODULE mo_functions


text:
Name: Mistress Althaea

Office: Trained Monkey

Personality: Definately bizarre.  If normalcy was legally ordained in this country she would be locked away forever.  Fortunately she looks normal enough, it is only when she starts pulling out the stencils and vaginas that people start to worry.

Contact: Its best to keep back 20 feet from this one at all times.

http://people.tribe.net/riotgrrl/

20060609 18:42:58 nbsp I liked your voter guide, but i think you rasterized the type at a low resolution. What distro are you running, btw? Users/ArlenAbraham

20060609 19:53:01 nbsp I liked the voter guide too. Would it be possible to use embedded fonts? Arlen, if it helps you any, this is what pdfinfo tells me about that document: http://www.bayworld.net/~nick/PERBvoterguide_1.pdf.txt. Users/NickSchmalenberger

20060628 22:07:25 nbsp I am running Ubuntu 6.06, WinXP and OS X on my dell box, thank you very much for the opportunity to brag about it. Users/MistressAlthaea

20060628 22:08:53 nbsp And the copy you got was probably a lores jpg version b/c i was lazy when creating the smaller pdf.  If you want to teach me better ways, I am all ears, oh wise ones. Users/MistressAlthaea

20060628 22:13:23 nbsp My comments on the Ubuntu 6.06 release: http://people.tribe.net/riotgrrl/blog/07621b59222a439fb6e449ed9fbafcef Users/MistressAlthaea

20060628 23:00:17 nbsp Kubuntu + EasyUbuntu, myself.  Ive got a server on Ubuntu Server, and if all goes well, by 2007 well be migrating everything over to Ubuntu Server from SUSE.  Nice little distro. Users/JabberWokky

20061103 17:44:54 nbsp I really like the voters guide, but you do not appear to be using outline fonts except in pages 24. Do you think you could fix this, it makes most of the document somewhat hard to read. Thanks. Users/NickSchmalenberger

20061103 17:50:08 nbsp If you need advice about InDesign, I know the editor of The Flatlander, and she uses InDesign. She might be able to give you some advice, and she would probably like the guide too. I can put you in contact with her if you want. Users/NickSchmalenberger

20090113 17:13:35 nbsp It was nice seeing you for winterness! More monkeys later? Users/BeckWehrle


PROGRAM stats_3
!
!  Purpose:
!    To calculate mean and the standard deviation of an input
!    data set, where each input value can be positive, negative,
!    or zero.
!
!  Record of revisions:
!      Date       Programmer          Description of change
!      ====       ==========          =====================
!    11/13/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Data dictionary: declare variable types, definitions, & units
INTEGER :: i         ! Loop index
INTEGER :: n = 0     ! The number of input samples.
REAL :: std_dev      ! The standard deviation of the input samples.
REAL :: sum_x = 0.   ! The sum of the input values.
REAL :: sum_x2 = 0.  ! The sum of the squares of the input values.
REAL :: x = 0.       ! An input data value.
REAL :: x_bar        ! The average of the input samples.

! Get the number of points to input.
WRITE (*,*) 'Enter number of points: '
READ  (*,*) n

! Check to see if we have enough input data.
IF ( n < 2 ) THEN ! Insufficient data

   WRITE (*,*) 'At least 2 values must be entered.'

ELSE ! we will have enough data, so let's get it.

   ! Loop to read input values.
   DO i = 1, n

      ! Read values
      WRITE (*,*) 'Enter number: '
      READ  (*,*) x
      WRITE (*,*) 'The number is ', x

      ! Accumulate sums.
      sum_x  = sum_x + x
      sum_x2 = sum_x2 + x**2

   END DO

   ! Now calculate statistics.
   x_bar = sum_x / real(n)
   std_dev = SQRT((real(n)*sum_x2 - sum_x**2) / (real(n)*real(n-1)))

   ! Tell user.
   WRITE (*,*) 'The mean of this data set is:', x_bar
   WRITE (*,*) 'The standard deviation is:   ', std_dev
   WRITE (*,*) 'The number of data points is:', n

END IF

END PROGRAM stats_3


text:
program textposition
    use kernel32
    implicit none
    integer(HANDLE) :: hConsole
    integer(BOOL) :: q

    hConsole = GetStdHandle(STD_OUTPUT_HANDLE)
    q = SetConsoleCursorPosition(hConsole, T_COORD(3, 6))
    q = WriteConsole(hConsole, loc("Hello"), 5, NULL, NULL)
end program

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData