text:
*###[ ffdxc0:
	subroutine ffdxc0(cs3,ipi12,isoort,clogi,ilogi,xpi,dpipj,piDpj,
     +		xqi,dqiqj,qiDqj,sdel2,del2s,etalam,etami,delpsi,alph,
     +		ddel2s,ldel2s,npoin,ier)
***#[*comment:***********************************************************
*									*
*	Calculates the difference of two threepoint functions		*
*	C(3,...a) - C(4,...b)						*
*	For this we not only calculate the roots of the three-point	*
*	function y,z(1-4,3-4,1-3) but also the combinations		*
*									*
*		yzzy = y(,4,)*z(,3,) - z(,4,)*y(,3,)			*
*	and								*
*		yyzz = y(,4,) - z(,4,) - y(,3,) + z(,3,)		*
*									*
*	This is done explicitly for most special cases, so a lot of	*
*	lines of code result.  This may be shortened with a smart use	*
*	of indices, however, it is readable now.			*
*									*
*	Input:	xpi(6,3:4)	(real)	transformed mi,pi squared in Ci	*
*		dpipj(6,6,3:4)	(real)	xpi(i)-xpi(j)			*
*		piDpj(6,6,3:4)	(real)	pi(i).pi(j)			*
*		xqi(10,10)	(real)	transformed mi,pi squared in D	*
*		dqiqj(10,10)	(real)	xqi(i)-xqi(j)			*
*		qiDqj(10,10)	(real)	qi(i).qi(j)			*
*		sdel2		(real)	sqrt(delta_{p_1 p_2}^{p_1 p_2})	*
*		del2s(3,3:4)	(real)	delta_{p_i s_i}^{p_i s_i}	*
*		etalam(3:4)	(real)	delta_{s_1 s_2 s_3}^{s_1 s_2 s_3}
*					  /delta_{p_1 p_2}^{p_1 p_2}	*
*		etami(6,3:4)	(real)	m_i^2 - etalam			*
*		ddel2s(2:3)	(real)	del2s(i,3) - del2s(i,4)		*
*		alph(3)		(real)	alph(1)=alpha, alph(3)=1-alpha	*
*		ldel2s	 (logical)	indicates yes/no limit del2s->0	*
*									*
*	Output: cs3	 (complex)(160)	C0(3)-C0(4), not yet summed.	*
*		ipi12	 (integer)(6)	factors pi^2/12, not yet summed	*
*		slam	 (complex)	lambda(p1,p2,p3).		*
*		isoort	 (integer)(16)	indication of he method used	*
*		clogi	 (complex)(6)	log(-dyz(2,1,i)/dyz(2,2,i))	*
*		ilogi	 (integer)(6)	factors i*pi in this		*
*		ier	 (integer)	0=ok, 1=inaccurate, 2=error	*
*									*
*	Calls:	...							*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments:
*
	integer ipi12(16),isoort(16),ilogi(6),npoin,ier
	logical ldel2s
	DOUBLE COMPLEX cs3(160),clogi(6)
	DOUBLE PRECISION xqi(10),dqiqj(10,10),qiDqj(10,10),
     +		xpi(6,3:4),dpipj(6,6,3:4),piDpj(6,6,3:4),
     +		sdel2,del2s(3,3:4),etalam(3:4),etami(6,3:4),alph(3),
     +		ddel2s(2:3),delpsi(3,3:4)
*
*	local variables:
*
	integer i,j,k,l,ip,ier0,ii,ifirst,ieri(12),idone(6)
	logical lcompl
	DOUBLE COMPLEX c,csom,chck,cs(5),csdeli(3,3:4),csdel2,
     +		cy(4,3:4,3),cz(4,3:4,3),cdyz(2,2,3:4,3),cd2yzz(3:4,3),
     +		cpi(6,3:4),cpiDpj(6,6,3:4),cdyzzy(4,3),cdyyzz(2,3)
	DOUBLE PRECISION sdel2i(3,3:4),s(5),som,smax,absc,dfflo1,xhck,
     +		rloss,y(4,3:4,3),z(4,3:4,3),dyz(2,2,3:4,3),d2yzz(3:4,3),
     +		dy2z(4,3:4,3),dyzzy(4,3),dsdel2,xmax
	DOUBLE COMPLEX zxfflg,zfflog,zfflo1
*
*	common blocks:
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ check input:
	if ( ltest ) then
*	    call ffxhck(xpi(1,3),dpipj(1,1,3),6,ier)
*	    call ffxhck(xpi(1,4),dpipj(1,1,4),6,ier)
	    call ffxhck(xqi,dqiqj,10,ier)
	endif
*  #] check input:
*  #[ get y,z-roots:
	lcompl = .FALSE.
	if ( lwrite ) print '(a)','  ##[ get roots:'
	do 20 k=3,4
	do 10 i=1,3
*
*	get roots (y,z) and flag what to do: 0=nothing, 1=normal,
*	-1=complex
*
	ip = i+3
*	first get the roots
	if ( del2s(i,k) .le. 0 ) then
*	    real case
	    sdel2i(i,k) = sqrt(-del2s(i,k))
*	    then handle the special case Si = 0
	    if ( xpi(ip,k) .eq. 0 ) then
		if ( i .eq. 1 .and. alph(3) .eq. 0 .or.
     +		     i .eq. 3 .and. alph(1) .eq. 0 ) then
		    isoort(2*i-1+8*(k-3)) = 0
		    isoort(2*i+8*(k-3)) = 0
		    goto 10
		endif
	    endif
	    call ffxxyz(y(1,k,i),z(1,k,i),dyz(1,1,k,i),d2yzz(k,i),
     +		dy2z(1,k,i),i,sdel2,sdel2i(i,k),etalam(k),etami(1,k),
     +		delpsi(i,k),xpi(1,k),dpipj(1,1,k),piDpj(1,1,k),
     +		isoort(2*i-1+8*(k-3)),ldel2s,6,ier)
	else
*	    complex case
	    sdel2i(i,k) = sqrt(del2s(i,k))
	    csdeli(i,k) = DCMPLX(x0,sdel2i(i,k))
	    lcompl = .TRUE.
	    call ffcxyz(cy(1,k,i),cz(1,k,i),cdyz(1,1,k,i),cd2yzz(k,i),i,
     +		sdel2,sdel2i(i,k),etalam(k),etami(1,k),delpsi(i,k),xpi(
     +		1,k),piDpj(1,1,k),isoort(2*i-1+8*(k-3)),ldel2s,6,ier)
	endif
   10	continue
   20	continue
*  #] get y,z-roots:
*  #[ convert to complex if necessary:
	do 60 i=2,3
	    l = 2*i-1
	    if ( isoort(l).gt.0 .and. isoort(l+8).lt.0 ) then
		k = 3
*		we get  -5, -105 if they have equal roots, isoort=+2
*			-6, -106 if they have unequal roots, isoort=+1
		if ( .not.ldel2s ) then
		    isoort(l) = isoort(l)-7
		    isoort(l+1) = isoort(l+1)-7
		else
		    isoort(l) = isoort(l)-207
		    isoort(l+1) = isoort(l+1)-207
		endif
	    elseif ( isoort(l).lt.0 .and. isoort(l+8).gt.0 ) then
		k = 4
		if ( .not.ldel2s ) then
		    isoort(l+8) = isoort(l+8)-7
		    isoort(l+9) = isoort(l+9)-7
		else
		    isoort(l+8) = isoort(l+8)-207
		    isoort(l+9) = isoort(l+9)-207
		endif
	    else
		k = 0
	    endif
	    if ( k .ne. 0 ) then
		if ( lwrite ) print *,'ffdxc0: converting i,k=',i,k,
     +			' to complex'
		do 30 j=1,4
		    cy(j,k,i) = y(j,k,i)
		    cz(j,k,i) = z(j,k,i)
   30		continue
		do 50 j=1,2
		    do 40 l=1,2
			cdyz(l,j,k,i) = dyz(l,j,k,i)
   40		    continue
   50		continue
		cd2yzz(k,i) = d2yzz(k,i)
		csdeli(i,k) = sdel2i(i,k)
	    endif
   60	continue
*  #] convert to complex if necessary:
*  #[ get differences:
*
*	the only important differences are y4z3-z3y4 and (1-y4)(1-z3)-
*	(1-y3)(1-z4)
*
	do 100 i=1,12
	    ieri(i) = 0
  100	continue
*  #[	vertices (1):
	som = qiDqj(7,2)/sdel2
	if ( isoort(1) .ge. 0 ) then
*	    Note that the isoorts are equal as the vertex is equal.
*
*	    flag if we have a cancellation
*
	    if ( abs(som) .lt. xloss ) then
		isoort(1) = isoort(1) + 10
		isoort(9) = isoort(9) + 10
	    endif
	    do 110 k=1,4
		dyzzy(k,1) = som*z(k,3,1)
		if ( k .gt. 2 ) dyzzy(k,1) = -dyzzy(k,1)
  110	    continue
	else
	    if ( abs(som) .lt. xloss ) then
		isoort(1) = isoort(1) - 10
		isoort(9) = isoort(9) - 10
	    endif
	    do 120 k=1,4
		cdyzzy(k,1) = DBLE(som)*cz(k,3,1)
		if ( k .gt. 2 ) cdyzzy(k,1) = -cdyzzy(k,1)
  120	    continue
	    cdyyzz(1,1) = som
	    cdyyzz(2,1) = som
	    if ( lwrite ) then
		print *,'cdyyzz(11) =',cy(2,4,1)-cy(2,3,1),
     +			absc(cy(2,4,1))
		print *,'cdyyzz(11)+=',cdyyzz(1,1)
	    endif
	endif
*  #]	vertices (1):
*  #[	vertices (2):
	if ( isoort(3) .ge. 0 ) then
*  #[	    real case: (note that this implies isoort(11)>0)
	    ifirst = 0
	    do 150 j=1,2
	    do 140 k=1,2
		ii = 2*(j-1) + k
		dyzzy(ii,2) = y(2*j,4,2)*z(ii,3,2)-y(2*j,3,2)*z(ii,4,2)
		xmax = abs(y(2*j,4,2)*z(ii,3,2))
		if ( abs(dyzzy(ii,2)) .ge. xmax ) goto 140
		isoort(3) = isoort(3) + 10
		isoort(11) = isoort(11) + 10
 1000		format(a,i1,a,g22.14,g12.4)
		if ( lwrite ) print 1000,'dyzzy(',ii,'2)  = ',
     +			dyzzy(ii,2),xmax
		if ( ldel2s ) then
		    print *,'ffdxc0: not ready for del2s=0, real case'
		    goto 130
		endif
		if ( ifirst .le. 0 ) then
		    if ( ddel2s(2) .eq. 0 ) then
			dsdel2 = 0
		    else
			dsdel2 = ddel2s(2)/(sdel2i(2,3)+sdel2i(2,4))
		    endif
		endif
		if ( ifirst .le. 1 ) then
		    if ( j .eq. 1 ) then
			s(1) = xqi(6)*qiDqj(7,4)*qiDqj(5,4)/sdel2
			s(2) = -qiDqj(7,4)*sdel2i(2,3)
			s(3) = +qiDqj(6,4)*dsdel2
		    else
			s(1) = xqi(6)*qiDqj(7,2)*qiDqj(5,2)/sdel2
			s(2) = -qiDqj(7,2)*sdel2i(2,3)
			s(3) = +qiDqj(6,2)*dsdel2
		    endif
		endif
		if ( ifirst .le. 0 ) then
		    ifirst = 2
		    s(4) = -qiDqj(5,10)*qiDqj(7,4)*sdel2i(2,3)/sdel2
		    s(5) = delpsi(2,3)*dsdel2/sdel2
		endif
		if ( k .eq. 1 ) then
		    som = s(1) + s(2) + s(3) + s(4) + s(5)
		else
		    som = s(1) - s(2) - s(3) - s(4) - s(5)
		endif
		smax = max(abs(s(1)),abs(s(2)),abs(s(3)),abs(s(4)),
     +			abs(s(5)))/xqi(6)**2
		if ( lwrite ) then
		    print 1000,'dyzzy(',ii,'2)+ = ',som/xqi(6)**2,smax
*		    print *,(s(i)/xqi(6)**2,i=1,5)
		endif
		if ( smax .lt. xmax ) then
		    dyzzy(ii,2) = som/xqi(6)**2
		    xmax = smax
		endif
  130		continue
		if ( lwarn .and. abs(dyzzy(ii,2)) .lt. xloss*xmax ) then
		    call ffwarn(140,ieri(2*k+j-2),dyzzy(ii,2),xmax)
		endif
  140	    continue
	    ifirst = ifirst - 1
  150	    continue
*  #]	    real case:
	else
*  #[	    complex case:
	    ifirst = 0
	    do 180 j=1,2
	    do 170 k=1,2
		ii = 2*(j-1) + k
		cdyzzy(ii,2) = cy(2*j,4,2)*cz(ii,3,2)-cy(2*j,3,2)*
     +			cz(ii,4,2)
		xmax = absc(cy(2*j,4,2)*cz(ii,3,2))
		if ( absc(cdyzzy(ii,2)) .ge. xmax ) goto 170
		isoort(3) = isoort(3) - 10
		isoort(11) = isoort(11) - 10
 1002		format(a,i1,a,2g22.14,g12.4)
		if ( lwrite ) print 1002,'cdyzzy(',ii,'2)  =',
     +						cdyzzy(ii,2),xmax
		if ( ldel2s ) then
		    ip = 3
		else
		    ip = 6
		endif
		if ( mod(isoort(3),10).ne.0 .or. mod(isoort(11),10).ne.0
     +								) then
*
*		    one of the roots is really real
*
		    if ( ifirst .le. 0 ) then
			csdel2=DBLE(ddel2s(2))/(csdeli(2,3)+csdeli(2,4))
		    endif
		    if ( ifirst .le. 1 ) then
			if ( j .eq. 1 .neqv. ldel2s ) then
			    if ( .not.ldel2s ) then
				cs(1)=xqi(6)*qiDqj(7,4)*qiDqj(5,4)/sdel2
				cs(2) = -DBLE(qiDqj(7,4))*csdeli(2,3)
				cs(3) = +DBLE(qiDqj(6,4))*csdel2
			    else
				cs(1)=-xqi(3)*qiDqj(5,10)*qiDqj(7,2)/
     +								sdel2
				cs(2) = -DBLE(qiDqj(7,2))*csdeli(2,3)
				cs(3) = -DBLE(qiDqj(6,3))*csdel2
			    endif
			else
			    cs(1) = xqi(ip)*qiDqj(7,2)*qiDqj(5,2)/sdel2
			    cs(2) = -DBLE(qiDqj(7,2))*csdeli(2,3)
			    cs(3) = +DBLE(qiDqj(ip,2))*csdel2
			endif
		    endif
		    if ( ifirst .le. 0 ) then
			ifirst = 2
			if ( .not.ldel2s ) then
			    cs(4) = -DBLE(qiDqj(5,10)*qiDqj(7,4)/sdel2)*
     +				csdeli(2,3)
			else
			    cs(4) = -DBLE(qiDqj(5,3)*qiDqj(7,2)/sdel2)*
     +				csdeli(2,3)
			endif
			cs(5) = DBLE(delpsi(2,3)/sdel2)*csdel2
		    endif
		else
*
*		    both roots are complex
*
		    if ( ifirst .eq. 0 ) then
			dsdel2 = -ddel2s(2)/(sdel2i(2,3)+sdel2i(2,4))
			csdel2 = DCMPLX(x0,dsdel2)
		    endif
		    if ( ifirst .le. 1 ) then
			if ( j .eq. 1 .neqv. ldel2s ) then
			    if ( .not.ldel2s ) then
				cs(1)=xqi(6)*qiDqj(7,4)*qiDqj(5,4)/sdel2
				cs(2)=-DCMPLX(x0,qiDqj(7,4)*sdel2i(2,3))
				cs(3)=+DCMPLX(x0,qiDqj(6,3)*dsdel2)
			    else
				cs(1)=-xqi(3)*qiDqj(5,10)*qiDqj(7,2)/
     +								sdel2
				cs(2)=-DCMPLX(x0,qiDqj(7,2)*sdel2i(2,3))
				cs(3)=-DCMPLX(x0,qiDqj(6,3)*dsdel2)
			    endif
			else
			    cs(1) = xqi(ip)*qiDqj(7,2)*qiDqj(5,2)/sdel2
			    cs(2) = -DCMPLX(x0,qiDqj(7,2)*sdel2i(2,3))
			    cs(3) = +DCMPLX(x0,qiDqj(ip,2)*dsdel2)
			endif
		    endif
		    if ( ifirst .eq. 0 ) then
			ifirst = 2
			if ( .not.ldel2s ) then
			    cs(4) = -DCMPLX(x0,qiDqj(5,10)*qiDqj(7,4)*
     +						sdel2i(2,3)/sdel2)
			else
			    cs(4) = -DCMPLX(x0,qiDqj(5,3)*qiDqj(7,2)*
     +						sdel2i(2,3)/sdel2)
			endif
			cs(5) = DCMPLX(x0,delpsi(2,3)*dsdel2/sdel2)
		    endif
		endif
		if ( k .eq. 1 ) then
		    csom = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
		else
		    csom = cs(1) - cs(2) - cs(3) - cs(4) - cs(5)
		endif
		smax = max(absc(cs(1)),absc(cs(2)),absc(cs(3)),
     +			absc(cs(4)),absc(cs(5)))/xqi(ip)**2
		if ( lwrite ) then
		    print 1002,'cdyzzy(',ii,'2)+ =',csom/DBLE(xqi(ip))**
     +			2,smax
***		    print *,(cs(i)/DBLE(xqi(ip))**2,i=1,5)
		endif
		if ( smax .lt. xmax ) then
		    cdyzzy(ii,2) = csom/DBLE(xqi(ip))**2
		    xmax = smax
		endif
		if ( lwarn .and. absc(cdyzzy(ii,2)).lt.xloss*xmax ) then
		  call ffwarn(140,ieri(2*k+j-2),absc(cdyzzy(ii,2)),xmax)
		endif
  170	    continue
*
*	    get cdyyzz
*
	    if ( ldel2s ) then
		cdyyzz(j,2) = cdyz(2,j,4,2) - cdyz(2,j,3,2)
		xmax = absc(cdyz(2,j,4,2))
		if ( absc(cdyyzz(j,2)) .ge. xloss*xmax ) goto 175
		if ( lwrite ) print 1002,'cdyyzz(',j,'2) =',cdyyzz(j,2),
     +			xmax
		if ( ifirst .le. 0 ) then
		    if ( mod(isoort( 3),10).ne.0 .or.
     +			 mod(isoort(11),10).ne.0 ) then
			csdel2=DBLE(ddel2s(2))/(csdeli(2,3)+csdeli(2,4))
		    else
			dsdel2 = -ddel2s(2)/(sdel2i(2,3)+sdel2i(2,4))
			csdel2 = DCMPLX(x0,dsdel2)
		    endif
		endif
		cs(2) = csdel2/DBLE(xqi(3))
		cs(1) = qiDqj(5,3)*qiDqj(7,2)/(sdel2*xqi(3))
		if ( j .eq. 1 ) then
		    csom = cs(1) + cs(2)
		else
		    csom = cs(1) - cs(2)
		endif
		smax = absc(cs(1))
		if ( lwrite ) print 1002,'cdyyzz(',j,'2)+=',csom,smax
		if ( smax .lt. xmax ) then
		    cdyyzz(j,2) = csom
		    xmax = smax
		endif
		if ( lwarn .and. absc(cdyyzz(j,2)).lt.xloss*xmax ) then
		    call ffwarn(147,ieri(7+j),absc(cdyyzz(j,2)),xmax)
		endif
	    endif
*
*	    bookkeeping
*
  175	    continue
	    ifirst = ifirst - 1
  180	    continue
*  #]	    complex case:
	endif
*  #]	vertices (2):
*  #[	vertices (3):
	if ( isoort(5) .ge. 0 ) then
*  #[	    real case: (note that this implies isoort(15)>0)
	    ifirst = 0
	    do 210 j=1,2
	    do 200 k=1,2
		ii = 2*(j-1) + k
		dyzzy(ii,3) = y(2*j,4,3)*z(ii,3,3)-y(2*j,3,3)*z(ii,4,3)
		xmax = abs(y(2*j,4,3)*z(ii,3,3))
		if ( abs(dyzzy(ii,3)) .ge. xmax ) goto 200
		isoort(5) = isoort(5) + 10
		isoort(13) = isoort(13) + 10
		if ( lwrite ) print 1000,'dyzzy(',ii,'3)  = ',
     +			dyzzy(ii,3),xmax
		if ( ldel2s ) then
		    print *,'ffdxc0: not ready for del2s=0, real case'
		    goto 190
		endif
		if ( ifirst .le. 0 ) then
		    if ( ddel2s(2) .eq. 0 ) then
			dsdel2 = 0
		    else
			dsdel2 = ddel2s(3)/(sdel2i(3,3)+sdel2i(3,4))
		    endif
		endif
		if ( ifirst .le. 1 ) then
		    if ( j .eq. 1 ) then
			s(1) = xqi(8)*qiDqj(7,1)*qiDqj(5,1)/sdel2
			s(2) = +qiDqj(7,1)*sdel2i(3,3)
			s(3) = +qiDqj(9,1)*dsdel2
		    else
			s(1) = xqi(8)*qiDqj(7,4)*qiDqj(5,4)/sdel2
			s(2) = +qiDqj(7,4)*sdel2i(3,3)
			s(3) = +qiDqj(9,4)*dsdel2
		    endif
		endif
		if ( ifirst .le. 0 ) then
		    ifirst = 2
		    s(4) = -qiDqj(5,9)*qiDqj(7,1)*sdel2i(3,3)/sdel2
		    s(5) = delpsi(3,3)*dsdel2/sdel2
		endif
		if ( k .eq. 1 ) then
		    som = s(1) + s(2) + s(3) + s(4) + s(5)
		else
		    som = s(1) - s(2) - s(3) - s(4) - s(5)
		endif
		smax = max(abs(s(1)),abs(s(2)),abs(s(3)),abs(s(4)),
     +			abs(s(5)))/xqi(8)**2
		if ( lwrite ) then
		    print 1000,'dyzzy(',ii,'3)+ = ',som/xqi(8)**2,smax
***		    print *,(s(i)/xqi(8)**2,i=1,5)
		endif
		if ( smax .lt. xmax ) then
		    dyzzy(ii,3) = som/xqi(8)**2
		    xmax = smax
		endif
  190		continue
		if ( lwarn .and. abs(dyzzy(ii,3)) .lt. xloss*xmax ) then
		    call ffwarn(140,ieri(2*k+j+2),dyzzy(ii,3),xmax)
		endif
  200	    continue
	    ifirst = ifirst - 1
  210	    continue
*  #]	    real case:
	else
*  #[	    complex case:
	    ifirst = 0
	    do 240 j=1,2
	    do 230 k=1,2
		ii = 2*(j-1) + k
		cdyzzy(ii,3) = cy(2*j,4,3)*cz(ii,3,3)-cy(2*j,3,3)*
     +			cz(ii,4,3)
		xmax = absc(cy(2*j,4,3)*cz(ii,3,3))
		if ( absc(cdyzzy(ii,3)) .ge. xmax ) goto 230
		isoort(5) = isoort(5) - 10
		isoort(13) = isoort(13) - 10
		if ( lwrite ) print 1002,'cdyzzy(',ii,'3)  =',
     +						cdyzzy(ii,3),xmax
		if ( ldel2s ) then
		    ip = 3
		else
		    ip = 8
		endif
		if ( mod(isoort(3),10).ne.0 .or. mod(isoort(11),10).ne.0
     +								) then
*
*		    one of the roots is really real
*
		    if ( ifirst .le. 0 ) then
			csdel2=DBLE(ddel2s(3))/(csdeli(3,3)+csdeli(3,4))
		    endif
		    if ( ifirst .le. 1 ) then
			if ( j .eq. 1 ) then
			    cs(1) = xqi(ip)*qiDqj(7,1)*qiDqj(5,1)/sdel2
			    cs(2) = +DBLE(qiDqj(7,1))*csdeli(3,3)
			    if ( .not.ldel2s ) then
				cs(3) = +DBLE(qiDqj(9,1))*csdel2
			    else
				cs(3) = +DBLE(qiDqj(3,1))*csdel2
			    endif
			else
			    if ( .not.ldel2s ) then
				cs(1) = xqi(ip)*qiDqj(7,4)*qiDqj(5,4)/
     +								sdel2
				cs(2) = DBLE(qiDqj(7,4))*csdeli(3,3)
			    else
				cs(1) = xqi(ip)*qiDqj(7,1)*qiDqj(5,9)/
     +								sdel2
				cs(2) = DBLE(qiDqj(7,1))*csdeli(3,3)
			    endif
			    cs(3) = +DBLE(qiDqj(9,3))*csdel2
			endif
			if ( ldel2s ) cs(3) = -cs(3)
		    endif
		    if ( ifirst .le. 0 ) then
			ifirst = 2
			if ( .not.ldel2s ) then
			    cs(4) = -DBLE(qiDqj(5,9)*qiDqj(7,1)/sdel2)*
     +				csdeli(3,3)
			else
			    cs(4) = DBLE(qiDqj(5,4)*qiDqj(7,1)/sdel2)*
     +				csdeli(3,3)
			endif
			cs(5) = DBLE(delpsi(3,3)/sdel2)*csdel2
		    endif
		else
*
*		    both roots are complex
*
		    if ( ifirst .eq. 0 ) then
			dsdel2 = -ddel2s(3)/(sdel2i(3,3)+sdel2i(3,4))
			csdel2 = DCMPLX(x0,dsdel2)
		    endif
		    if ( ifirst .le. 1 ) then
			if ( j .eq. 1 ) then
			    cs(1) = xqi(ip)*qiDqj(7,1)*qiDqj(5,1)/sdel2
			    cs(2) = +DCMPLX(x0,qiDqj(7,1)*sdel2i(3,3))
			    if ( .not.ldel2s ) then
				cs(3) = +DCMPLX(x0,qiDqj(9,1)*dsdel2)
			    else
				cs(3) = +DCMPLX(x0,qiDqj(3,1)*dsdel2)
			    endif
			else
			    if ( .not.ldel2s ) then
				cs(1) = xqi(ip)*qiDqj(7,4)*qiDqj(5,4)/
     +								sdel2
				cs(2) =DCMPLX(x0,qiDqj(7,4)*sdel2i(3,3))
			    else
				cs(1) = xqi(ip)*qiDqj(7,1)*qiDqj(5,9)/
     +								sdel2
				cs(2) =DCMPLX(x0,qiDqj(7,1)*sdel2i(3,3))
			    endif
			    cs(3) = +DCMPLX(x0,qiDqj(9,3)*dsdel2)
			endif
			if ( ldel2s ) cs(3) = -cs(3)
		    endif
		    if ( ifirst .le. 0 ) then
			ifirst = 2
			if ( .not.ldel2s ) then
			    cs(4) = -DCMPLX(x0,qiDqj(5,9)*qiDqj(7,1)*
     +						sdel2i(3,3)/sdel2)
			else
			    cs(4) = DCMPLX(x0,qiDqj(5,4)*qiDqj(7,1)*
     +						sdel2i(3,3)/sdel2)
			endif
			cs(5) = DCMPLX(x0,delpsi(3,3)*dsdel2/sdel2)
		    endif
		endif
		if ( k .eq. 1 ) then
		    csom = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
		else
		    csom = cs(1) - cs(2) - cs(3) - cs(4) - cs(5)
		endif
		smax =max(absc(cs(1)),absc(cs(2)),absc(cs(3)),
     +			absc(cs(4)),absc(cs(5)))/xqi(ip)**2
		if ( lwrite ) then
		    print 1002,'cdyzzy(',ii,'3)+ =',csom/DBLE(xqi(ip))**
     +			2,smax
***		    print *,(cs(i)/DBLE(xqi(ip))**2,i=1,5)
		endif
		if ( smax .lt. xmax ) then
		    cdyzzy(ii,3) = csom/DBLE(xqi(ip))**2
		    xmax = smax
		endif
		if ( lwarn .and. absc(cdyzzy(ii,3)).lt.xloss*xmax ) then
		  call ffwarn(140,ieri(2*k+j+2),absc(cdyzzy(ii,3)),xmax)
		endif
  230	    continue
*
*	    get cdyyzz
*
	    if ( ldel2s ) then
		cdyyzz(j,3) = cdyz(2,j,4,3) - cdyz(2,j,3,3)
		xmax = absc(cdyz(2,j,4,3))
		if ( absc(cdyyzz(j,3)) .ge. xloss*xmax ) goto 235
		if ( lwrite ) print 1002,'cdyyzz(',j,'3) =',cdyyzz(j,3),
     +			xmax
		if ( ifirst .le. 0 ) then
		    if ( mod(isoort( 5),10).ne.0 .or.
     +			 mod(isoort(13),10).ne.0 ) then
			csdel2=DBLE(ddel2s(3))/(csdeli(3,3)+csdeli(3,4))
		    else
			dsdel2 = -ddel2s(3)/(sdel2i(3,3)+sdel2i(3,4))
			csdel2 = DCMPLX(x0,dsdel2)
		    endif
		endif
		cs(2) = -csdel2/DBLE(xqi(3))
		cs(1) = qiDqj(5,3)*qiDqj(7,1)/(sdel2*xqi(3))
		if ( j .eq. 1 ) then
		    csom = cs(1) + cs(2)
		else
		    csom = cs(1) - cs(2)
		endif
		smax = absc(cs(1))
		if ( lwrite ) print 1002,'cdyyzz(',j,'3)+=',csom,smax
		if ( smax .lt. xmax ) then
		    cdyyzz(j,3) = csom
		    xmax = smax
		endif
		if ( lwarn .and. absc(cdyyzz(j,3)).lt.xloss*xmax ) then
		    call ffwarn(147,ieri(9+j),absc(cdyyzz(j,3)),xmax)
		endif
	    endif
*
*	    bookkeeping
*
  235	    continue
	    ifirst = ifirst - 1
  240	    continue
*  #]	    complex case:
	endif
*  #]	vertices (3):
	ier0 = 0
	do 250 i = 1,12
	   ier0 = max(ier0,ieri(i))
  250	continue
	ier = ier + ier0
*  #] get differences:
*  #[ check differences:
	if ( ltest ) then
	    rloss = xloss**2*DBLE(10)**(-mod(ier,50))
	    do 300 i=1,3
		if ( isoort(2*i-1) .ne. isoort(2*i+7) ) goto 300
		do 290 j=1,2
		    if ( isoort(2*i-1) .ge. 0 ) then
			xhck = dyzzy(j,i) - y(2,4,i)*z(j,3,i)
     +				      + z(j,4,i)*y(2,3,i)
			if ( rloss*abs(xhck) .gt. precx*max(
     +				abs(y(2,4,i)*z(j,3,i)),
     +				abs(z(j,4,i)*y(2,3,i))) ) print *,
     +			'ffdxc0: error: ','dyzzy(',j,i,') <> terms, ',
     +				dyzzy(j,i),y(2,4,i)*z(j,3,i),
     +				z(j,4,i)*y(2,3,i),xhck
			xhck = dyzzy(j+2,i) - y(4,4,i)*z(j+2,3,i)
     +				       + z(j+2,4,i)*y(4,3,i)
			if ( rloss*abs(xhck) .gt. precx*max(
     +				abs(y(4,4,i)*z(j+2,3,i)),
     +				abs(z(j+2,4,i)*y(4,3,i))) ) print *,
     +			'ffdxc0: error: ','dyzzy(',j+2,i,') <> terms, ',
     +				dyzzy(j+2,i),y(4,4,i)*z(j+2,3,i),
     +				z(j+2,4,i)*y(4,3,i),xhck
		    else
			chck = cdyzzy(j,i) - cy(2,4,i)*cz(j,3,i)
     +				      + cz(j,4,i)*cy(2,3,i)
			if ( rloss*absc(chck) .gt. precc*max(
     +				abs(cy(2,4,i)*cz(j,3,i)),
     +				abs(cz(j,4,i)*cy(2,3,i))) ) print *,
     +			'ffdxc0: error: ','cdyzzy(',j,i,') <> terms, ',
     +				cdyzzy(j,i),cy(2,4,i)*cz(j,3,i),
     +				cz(j,4,i)*cy(2,3,i),chck
			chck = cdyzzy(j+2,i) - cy(4,4,i)*cz(j+2,3,i)
     +				       + cz(j+2,4,i)*cy(4,3,i)
			if ( rloss*absc(chck) .gt. precc*max(
     +				abs(cy(4,4,i)*cz(j+2,3,i)),
     +				abs(cz(j+2,4,i)*cy(4,3,i))) ) print *,
     +			'ffdxc0: error: ','cdyzzy(',j+2,i,') <> terms,',
     +				cdyzzy(j+2,i),cy(4,4,i)*cz(j+2,3,i),
     +				cz(j+2,4,i)*cy(4,3,i),chck
		    endif
  290		continue
  300	    continue
	endif
*  #] check differences:
*  #[ write output:
	if ( lwrite ) then
	    print *,'ffdxc0: found roots:'
	    do 320 k=3,4
	    do 310 i=1,3
		print *,'  k = ',i
		if ( isoort(2*i+8*(k-3)) .gt. 0 ) then
		    print *,'  ym,ym1 = ',y(1,k,i),y(3,k,i),
     +			' (not used)'
		    print *,'  yp,yp1 = ',y(2,k,i),y(4,k,i)
		    print *,'  zm,zm1 = ',z(1,k,i),z(3,k,i)
		    print *,'  zp,zp1 = ',z(2,k,i),z(4,k,i)
		elseif ( isoort(2*i+8*(k-3)) .eq. 0 ) then
		    if ( isoort(2*i-1+8*(k-3)) .eq. 0 ) then
			print *,'  no roots, all is zero'
		    else
			print *,'  yp,yp1 = ',y(2,k,i),y(4,k,i)
			print *,'  zp,zp1 = ',z(2,k,i),z(4,k,i)
		    endif
		else
		    print *,'  cym,cym1 = ',cy(1,k,i),cy(3,k,i),
     +			'(not used)'
		    print *,'  cyp,cyp1 = ',cy(2,k,i),cy(4,k,i)
		    print *,'  czm,czm1 = ',cz(1,k,i),cz(3,k,i)
		    print *,'  czp,czp1 = ',cz(2,k,i),cz(4,k,i)
		endif
  310	    continue
  320	    continue
	endif
	if ( lwrite ) print '(a)','  ##] get roots:'
*  #] write output:
*  #[ logarithms for 4point function:
*
*	Not yet made stable ...
*
	if ( npoin .eq. 4 ) then
	if ( lwrite ) print '(a)','  ##[ logarithms for Ai<0:'
	do 420 i = 1,3
	do 410 k = 3,4
	    ii = i+3*(k-3)
	    if ( ilogi(ii) .ne. -999 ) then
		idone(ii) = 0
		goto 410
	    endif
	    l = 2*i+8*(k-3)-1
	    if ((isoort(l).gt.0 .or. mod(isoort(l),10).le.-5) .and.
     +		(isoort(l+1).ge.0 .or. mod(isoort(l+1),10).le.-5)) then
*  #[		real case:
*
*		the real case (isoort=-5,-6: really real but complex for ffdcs)
*
		s(1) = -dyz(2,1,k,i)/dyz(2,2,k,i)
		if ( lwrite ) then
*		    fantasize imag part, but suppress error message
		    clogi(ii) = zxfflg(s(1),1,x1,ier0)
		    print *,'clogi  = ',clogi(ii)
		endif
		if ( abs(s(1)-1) .lt. xloss ) then
		    clogi(ii) = dfflo1(d2yzz(k,i)/dyz(2,2,k,i),ier)
		    ilogi(ii) = 0
		else
		    if ( abs(s(1)+1) .lt. xloss ) then
			clogi(ii) = dfflo1(-2*sdel2i(i,k)/(xpi(i+3,k)*
     +				dyz(2,2,k,i)),ier)
		    else
			clogi(ii) = zxfflg(abs(s(1)),0,x0,ier)
		    endif
		    if ( dyz(2,2,k,i).gt.0 .and. dyz(2,1,k,i).gt.0 )
     +								then
			ilogi(ii) = -1
		    elseif ( dyz(2,1,k,i).lt.0 .and. dyz(2,2,k,i).lt.0)
     +								then
			ilogi(ii) = +1
		    else
			ilogi(ii) = 0
		    endif
*		    in case del2s=0 and i=3 we pick up a minus sign, I think
		    if ( ldel2s .and. i .eq. 3 ) ilogi(ii) = -ilogi(ii)
		endif
		if ( lwrite ) print *,'clogi+ = ',clogi(ii)+
     +			DCMPLX(x0,pi)*ilogi(ii)
		idone(ii) = 1
*  #]		real case:
	    elseif ( isoort(l) .lt. 0 ) then
*  #[		complex case:
*		for stability split the unit circle up in 4*pi/2
*		(this may have to be improved to 8*pi/4...)
*
		ier0 = 0
		if ( lwrite ) then
		    if ( abs(DBLE(cdyz(2,1,k,i))) .lt. xalog2 .or.
     +			abs(DIMAG(cdyz(2,2,k,i))) .lt. xalog2 ) then
			csom = -DCMPLX(DBLE(cdyz(2,1,k,i))/xalog2,DIMAG(
     +				cdyz(2,1,k,i))/xalog2) /DCMPLX(DBLE(cdyz
     +				(2,2,k,i))/xalog2,DIMAG(cdyz(2,2,k,i))/
     +				xalog2)
		    else
			csom = -cdyz(2,1,k,i)/cdyz(2,2,k,i)
		    endif
		    clogi(ii)=zfflog(csom,0,c0,ier0)
		    print *,'isoort = ',isoort(2*i-1)
		    print *,'cdyz(2,1) = ',cdyz(2,1,k,i)
		    print *,'cdyz(2,2) = ',cdyz(2,2,k,i)
		    print *,'clogi = ',clogi(ii)
		endif
		if ( DBLE(cdyz(2,1,k,i)) .gt. abs(DIMAG(cdyz(2,1,k,i))))
     +								then
		    som =2*atan2(DIMAG(cdyz(2,1,k,i)),DBLE(
     +							cdyz(2,1,k,i)))
		    clogi(ii) = DCMPLX(x0,som)
		    if ( DIMAG(cdyz(2,1,k,i)) .gt. 0 ) then
			ilogi(ii) = -1
		    else
			ilogi(ii) = +1
		    endif

		elseif ( DBLE(cdyz(2,1,k,i)) .lt.
     +			-abs(DIMAG(cdyz(2,1,k,i))) ) then
		    if ( DIMAG(cdyz(2,1,k,i)) .eq. 0 ) then
			call fferr(82,ier)
			print *,'isoort = ',isoort(l),isoort(l+1)
		    endif
		    som = 2*atan2(-DIMAG(cdyz(2,1,k,i)),-DBLE(
     +							cdyz(2,1,k,i)))
		    clogi(ii) = DCMPLX(x0,som)
		    if ( DIMAG(cdyz(2,1,k,i)) .gt. 0 ) then
			ilogi(ii) = +1
		    else
			ilogi(ii) = -1
		    endif
		else
		    s(1) = -DBLE(cdyz(2,1,k,i))
		    s(2) = DIMAG(cdyz(2,1,k,i))
		    som = 2*atan2(s(1),s(2))
		    clogi(ii) = DCMPLX(x0,som)
		    ilogi(ii) = 0
		endif
		if ( lwrite ) print *,'clogi+= ',clogi(ii)+
     +			DCMPLX(x0,pi)*ilogi(ii)
		idone(ii) = 1
*  #]		complex case:
	    endif
*	    Note that we generate an error if isoort(l)=0
	    if ( lwrite ) then
		print *,'ffdxc0: ',ii,': ',clogi(ii),' + ',ilogi(ii),
     +			'*i*pi'
	endif
  410	continue
	if ( idone(ii) .ne. 0 .and. idone(ii-3) .ne. 0 .and.
     +	     absc(clogi(ii)-clogi(ii-3)).lt.xloss*absc(clogi(ii)) .and.
     +	     ilogi(ii).eq.ilogi(ii-3) ) then
*  #[	    subtract more smartly:
	    if ( isoort(l).gt.0 .and. isoort(l+1).ge.0 ) then
		if ( lwrite ) print *,'ffdxc0: extra logs not ready ',
     +			'in the real case'
		goto 420
	    else
		cs(1) = cdyzzy(1,i)
		cs(2) = cdyzzy(2,i)
		if ( i .eq. 1 ) then
		    cs(3) = 0
		else
		    if ( lwrite ) print *,'ffdxc0: extra logs not ',
     +			'ready for i <>1'
		    goto 420
		endif
		csom = cs(1) - cs(2) + cs(3)
		xmax = max(absc(cs(1)),absc(cs(2)),absc(cs(3)))
*		change this to "no warning and quit" later
		if ( lwarn .and. absc(csom) .lt. xmax ) then
		    goto 420
***		    call ffwarn(148,ier,absc(csom),xmax)
		endif
		if ( lwrite ) print *,'som was : ',clogi(ii-3)-clogi(ii)
		c = csom/(cdyz(2,2,3,i)*cdyz(2,1,4,i))
		c = zfflo1(c,ier)
		if ( lwrite ) print *,'som is  : ',c
*
*		the log is never much bigger than 1, so demand at least
*		accuracy to 0.1; this will catch all i*pi errors
*
		if ( abs(clogi(ii-3)-clogi(ii)-c).gt.0.1 ) then
		    print *,'ffdxc0: error in smart logs: ',clogi(ii-3)-
     +		    	clogi(ii),c,' not used'
		    goto 420
		endif
		clogi(ii-3) = c
		clogi(ii) = 0
	    endif
*  #]	    subtract more smartly:
	endif
  420	continue
*	An algorithm to obtain the sum of two small logarithms more
*	accurately has been put in ffcc0p, not yet here
	if ( lwrite ) print '(a)','  ##] logarithms for Ai<0:'
	endif
*  #] logarithms for 4point function:
*  #[ real case integrals:
	if ( .not. lcompl ) then
*	    normal case
	    do 510 i=1,3
		if ( lwrite ) print '(a,i1,a)','  ##[ dxs nr ',i,':'
		j = 2*i-1
		if ( isoort(j) .eq. 0 ) then
		    if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j),isoort(j+1)
		    endif
		    if ( isoort(j+8) .eq. 0 ) then
			if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j+8),isoort(j+9)
			endif
		    else
			call ffcxs3(cs3(20*i+61),ipi12(j+8),y(1,4,i),
     +			   z(1,4,i),dyz(1,1,4,i),d2yzz(4,i),dy2z(1,4,i),
     +			   xpi(1,4),piDpj(1,1,4),i,6,isoort(j+8),ier)
		    endif
		elseif ( isoort(j+8) .eq. 0 ) then
		    if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j),isoort(j+1)
		    endif
		    call ffcxs3(cs3(20*i-19),ipi12(j),y(1,3,i),
     +			   z(1,3,i),dyz(1,1,3,i),d2yzz(3,i),dy2z(1,3,i),
     +			   xpi(1,3),piDpj(1,1,3),i,6,isoort(j),ier)
		else
		    call ffdcxs(cs3(20*i-19),ipi12(j),y(1,3,i),z(1,3,i),
     +			dyz(1,1,3,i),d2yzz(3,i),dy2z(1,3,i),dyzzy(1,i),
     +			xpi,piDpj,i,6,isoort(j),ier)
		endif
		if ( lwrite ) print '(a,i1,a)','  ##] dxs nr ',i,':'
  510	    continue
	    isoort(7) = 0
	    isoort(8) = 0
*  #] real case integrals:
*  #[ complex case integrals:
	else
*	    convert xpi
	    do 540 k=3,4
*not	    cetami(1,k) = etami(1,k)
*used	    cetami(3,k) = etami(3,k)
	    do 530 i=1,6
		cpi(i,k) = xpi(i,k)
		do 520 j=1,6
		    cpiDpj(j,i,k) = piDpj(j,i,k)
  520		continue
  530	    continue
  540	    continue
	    do 550 i=1,3
		if ( lwrite ) print '(a,i1,a)','  ##[ dcs nr ',i,':'
		j = 2*i-1
		if ( isoort(j) .eq. 0 ) then
		    if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j),isoort(j+1)
		    endif
		    if ( isoort(j+8) .eq. 0 ) then
			if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j+8),isoort(j+9)
			endif
		    else
			call ffcxs3(cs3(20*i+61),ipi12(j+8),y(1,4,i),
     +			   z(1,4,i),dyz(1,1,4,i),d2yzz(4,i),dy2z(1,4,i),
     +			   xpi(1,4),piDpj(1,1,4),i,6,isoort(j+8),ier)
		    endif
		elseif ( isoort(j+8) .eq. 0 ) then
		    if ( lwrite ) then
			print *,'ffdxc0: xk=0, ma=mb/Si-0 -> S3 = 0'
			print *,'isoort:',isoort(j),isoort(j+1)
		    endif
		    call ffcxs3(cs3(20*i-19),ipi12(j),y(1,3,i),
     +			   z(1,3,i),dyz(1,1,3,i),d2yzz(3,i),dy2z(1,3,i),
     +			   xpi(1,3),piDpj(1,1,3),i,6,isoort(j),ier)
		elseif ( isoort(j) .gt. 0 ) then
		    if ( isoort(j+8) .gt. 0 ) then
			call ffdcxs(cs3(20*i-19),ipi12(j),y(1,3,i),
     +			   z(1,3,i),dyz(1,1,3,i),d2yzz(3,i),dy2z(1,3,i),
     +			   dyzzy(1,i),xpi,piDpj,i,6,isoort(j),ier)
		    else
			print *,'ffdxc0: error: should not occur!'
			call ffcxs3(cs3(20*i-19),ipi12(j),y(1,3,i),
     +			   z(1,3,i),dyz(1,1,3,i),d2yzz(3,i),dy2z(1,3,i),
     +			   xpi(1,3),piDpj(1,1,3),i,6,isoort(j),ier)
			call ffcs3(cs3(20*i+61),ipi12(j+8),cy(1,4,i),
     +			    cz(1,4,i),cdyz(1,1,4,i),cd2yzz(4,i),
     +			    cpi(1,4),cpiDpj(1,1,4),i,6,isoort(j+8),ier)
		    endif
		else
		    if ( isoort(j+8) .lt. 0 ) then
			call ffdcs(cs3(20*i-19),ipi12(j),cy(1,3,i),
     +			    cz(1,3,i),cdyz(1,1,3,i),cd2yzz(3,i),
     +			    cdyzzy(1,i),cdyyzz(1,i),cpi,cpiDpj,
     +			    i,6,isoort(j),ier)
		    else
			print *,'ffdxc0: error: should not occur!'
			call ffcs3(cs3(20*i-19),ipi12(j),cy(1,3,i),
     +			    cz(1,3,i),cdyz(1,1,3,i),cd2yzz(3,i),
     +			    cpi(1,3),cpiDpj(1,1,3),i,6,isoort(j),ier)
			call ffcxs3(cs3(20*i+61),ipi12(j+8),y(1,4,i),
     +			   z(1,4,i),dyz(1,1,4,i),d2yzz(4,i),dy2z(1,4,i),
     +			   xpi(1,4),piDpj(1,1,4),i,6,isoort(j+8),ier)
		    endif
		endif
		if ( lwrite ) print '(a,i1,a)','  ##] dcs nr ',i,':'
  550	    continue
	    isoort(7) = 0
	    isoort(8) = 0
	endif
	return
*  #] complex case integrals:
*###] ffdxc0:
	end


Module nescgle
Use math_constants
Use omp_lib
Use scgle
Implicit None

contains

!Subroutine that calculates the alpha vector used for the calculation of
!the S(k,u(t)) structure
Subroutine Calc_alphaSku(k_v,SkF_v,D0_val,alpha_v)
	implicit none
	real * 8, dimension(:), intent(in) :: k_v, SkF_v
	real * 8, dimension(:), intent(out) :: alpha_v
	real * 8, intent (in) :: D0_val
	alpha_v = 2.d0 * (k_v ** 2) * D0_val / SkF_v
End Subroutine

!Subroutine that calculates the structure factor of the system at a given
!u material time through an instantaneous isochoric quench
Function sku_quench(Sk0_v,SkF_v,alpha_v,u_val) result(sku)
	implicit none
	real * 8, dimension(:), intent(in) :: Sk0_v, SkF_v, alpha_v
	real * 8, intent(in) :: u_val
	real * 8, dimension(size(Sk0_v)) :: Sku
	integer :: i1
	If (u_val > 0.d0) Then
		sku = SkF_v + ( Sk0_v - SkF_v ) * exp( - alpha_v * u_val )
	Else
		sku= Sk0_v
	End If
	return
End Function

!Subroutine that calculates the value of the material time u in which a system arrest
!calculated up to a max value of u and gamma when an instantaneous isochoric quench
!is applied to the system
Function calc_ua_isoch_inst_quench(Sk0,SkF,alpha,lam,k,kw,rho,d) result(ua_val)
	implicit none
	real * 8, dimension(:), intent(in) :: Sk0, SkF, alpha, lam, k, kw
	real * 8, dimension(size(Sk0)) :: sku
	real * 8, intent(in) :: rho
	real * 8 :: gam_val
	integer, intent(in) :: d
	real * 8 :: ua_val
	logical :: convergence
	real * 8 :: u_low, u_up, gam_test_u, gam_test_l, error, u_test, gam_test
	integer :: i1
	Real * 8, parameter :: gam_max_val = 1.d20, tol = 1.d-8
	Real * 8, parameter :: ua_max_val = 1.d20, ua_min_val = 1.d-20
	u_low = 0.d0
	u_up = 1.d-4
	sku = sku_quench( Sk0,SkF,alpha,u_up)
	gam_test_u = gamma(d,rho,k,kw,sku,lam)
	sku = sku_quench(Sk0,SkF,alpha,u_low)
	gam_test_l = gamma(d,rho,k,kw,sku,lam)
	If ( gam_test_u > gam_max_val .OR. gam_test_u /= gam_test_u ) Then
		convergence = .false.
		do while (convergence .eqv. .false.)
			u_low = u_up
			u_up  = u_up * 2.d0
			sku = sku_quench( Sk0,SkF,alpha,u_up)
			gam_test = gamma(d,rho,k,kw,sku,lam)
			If ( gam_test < gam_max_val .or. u_up > ua_max_val )	convergence = .true.
		end do
	else if (gam_test_l < gam_max_val ) then
		u_up = u_low
		return
	end if
	ua_val  = u_up
	gam_val = gam_test
	error =  abs(u_up - u_low) / u_up
	if ( error > tol ) convergence = .false.
	if ( u_up > ua_max_val ) convergence = .true.
	do while (convergence .eqv. .false.)
		u_test = (u_up + u_low) / 2.d0
		sku = sku_quench( Sk0,SkF,alpha,u_test)
		gam_test = gamma(d,rho,k,kw,sku,lam)
		if ( gam_test < gam_max_val ) then
			u_up = u_test
			gam_val = gam_test
		else
			u_low = u_test
		end if
		error =  abs(u_up - u_low) / u_up
		if ( error < tol ) convergence = .true.
	end do
	ua_val = u_up
	print*, ua_val, gam_val, error
	return
End Function

!Subroutine that search in an isochore the arrest temperature and then

Real * 8 Function gamma_u(Sk0_v,SkF_v,alpha_v,u_val,lam_v,k_v,kw_v,phi_val,d)
	implicit none
	real * 8, dimension(:), intent(in) :: Sk0_v, SkF_v, alpha_v, lam_v, k_v, kw_v
	real * 8, intent(in) :: u_val, phi_val
	integer, intent(in) :: d
	real * 8, dimension(size(Sk0_v)) :: Sku_v
	Real * 8 :: rho
	rho = phi_val * 6.d0 / pi
	!Call Calc_Sk_u(Sk0_v,SkF_v,alpha_v,Sku_v,u_val)
	gamma_u = gamma(d,rho,k_v,kw_v,Sku_v,lam_v)
	return
End Function

!This subroutine is for the writting of the structure factor as well as the
!weights of the quadrature and the information file
Subroutine su_si_sf_writting(k,sku,ski,skf,kw,i_unit)
	Implicit None
	Real * 8, dimension(:), intent(in) :: k, sku, ski, skf, kw
	Integer :: i1
	Integer, intent(in) :: i_unit
	!k S(k) kw WRITTING
	Write (i_unit,*) "# //	k //	S(k;u)	//	Si(k)	//	Sf(k)	//	Weights		//"
	Do i1=1, size(k)
		Write(i_unit,*) k(i1),sku(i1),ski(i1),skf(i1),kw(i1)
	End Do
End Subroutine

!Subroutine diagram_phi_T(phi,Ti,Tfl,system,k,kw,lamk,)
!	Implicit None
	!!! Volume/Area fraction and # density !!!
!	Real * 8 :: phi, rho
	!!! Quench Initial Temperature // Quench Final Temperature !!!
!	Real * 8 :: Ti, Tf
	!!! Structure factors !!!
!	Real * 8, dimension(:) :: sku, ski, skf
	!!! Wave vectors !!!
!	Real * 8, dimension(:) :: k, kw
	!!! Radial distribution functions !!!
!	Real * 8, dimension(:) :: gru, gri, grf
	!!! Radial vectors !!!
!	Real * 8, dimension(:) :: r, rw
	!!! Differential Radial vectors !!!
!	Real * 8, dimension(:) :: dr, drw
	!!! Integer for the wave-vector parameter !!!
!	Integer :: ikp
	!!! System and system approximation !!!
!	Character ( len = 4 ) :: sys, sys_approx
	!!! System parameters !!!
!	Real * 8, dimension(:) :: s_params
	!!! Writting Folder !!!
!	Character ( len = 100 ) :: folder
	!!! Writting Wave vectors !!!
!	Real * 8, dimension(:) :: k_wr, kw_wr
	!!! Writting Structure factors !!!
!	Real * 8, dimension(:) :: sku_wr, ski_wr, skf_wr
	!!! Writting File Names !!!
!	Character(len=100) :: gam_file_name, S_file_name, g_file_name, S_inf_file_name
!	Character(len=100) :: gam_inf_file_name
	!!! Dummy Variables !!!
!	Real * 8 :: ua
!	Integer :: i1,i2,i3
!	Do i1 = 1, 4000
!		phi = i1 * 1.d-4
!		rho = phi * 6.d0 / pi !3D
!		Ti = 4.d0
!		swlambda = 1.5d0
!		s_params(1) = phi
!		s_params(2) = Ti
!		Do i2 = 1, knp
!			s_params(ikp) = kv(i2)
!			SkI(i2) = structure("s0",sys,sys_approx,s_params)
!		End Do
!		Do i2 = 1, knpgr
!			s_params(ikp) = kvgr(i2)
!			skigr(i2) = structure("s0",sys,sys_approx,s_params)
!		End Do
!		i3 = 800
!		ua_val = 0.d0
!		call file_name_constructor(sys,s_params,folder,"gamma",id,".dat",S_file_name)
!		call file_name_constructor(sys,s_params,folder,"info_ua",id,".inf",inf_file_name)
!		open (unit=104, file=trim(S_file_name), status="Replace")
!		open (unit=105, file=trim(inf_file_name), status="Replace")
!		call k_info_writting(kv,k_quadp,k_approx,105)
!		call lambda_info_writting(kc,105)
!		write(104,*) "#phi,   Tf,   gamma,    ua"
!		close(105)
!		Do while (ua_val < 1.d20)
!			T_low = i3 * 1.d-4
!			s_params(2) = T_low
!			Print*, phi, T_low
!			call file_name_constructor(sys,s_params,folder,"s_ua",id,".dat",S_file_name)
!			call file_name_constructor(sys,s_params,folder,"g_ua",id,".dat",g_file_name)
!			call file_name_constructor(sys,s_params,folder,"info_ua",id,".inf",inf_file_name)
!			Print*, S_file_name
!			open (unit=101, file=trim(S_file_name), status="Replace")
!			open (unit=102, file=trim(g_file_name), status="Replace")
!			open (unit=103, file=trim(inf_file_name), status="Replace")
!			Do i2 = 1, knp
!				s_params(ikp) = kv(i2)
!				SkF(i2) = structure("s0",sys,sys_approx,s_params)
!			End Do
!			Print*, "Skf Done"
!			call Calc_alphaSku(kv,SkF,D0,alphaku)
!			Print*, "alphaku Done"
!			ua_val = calc_ua_isoch_inst_quench(SkI,SkF,alphaku,lamkv,kv,kwv,rho,d)
!			skv = sku_quench(SkI,SkF,alphaku,ua_val)
!			if (ua_val < 1.d20) Then
!			 gam_val = gamma(d,rho,kv,kwv,skv,lamkv)
!			 ua_writting = ua_val
!			 Ta = T_low
!			 Print*, "gamma=",gam_val, "ua=",ua_val
!			 Write(104,*) phi, T_low, gam_val, ua_val
!			End If
!			Do i2 = 1, knpgr
!				s_params(ikp) = kvgr(i2)
!				skfgr(i2) = structure("s0",sys,sys_approx,s_params)
!			End Do
!			call Calc_alphaSku(kvgr,skfgr,D0,alphakgr)
!			skvgru = sku_quench(skigr,skfgr,alphakgr,ua_val)
!			call calc_gr_dependent_functions(rho,kvgr,kwvgr,skvgru,rv,drv,drw,grv,coor_n_rv,d)
!			call s_writting(kvgr,skvgru,kwvgr,101)
!			call g_writting(rv,grv,coor_n_rv,102)
!			call k_info_writting(kvgr,kgr_quadp,k_approx,103)
!			call r_info_writting(rv,r_quadp,r_approx,103)
!			call dr_info_writting(drv,drw,dr_quadp,dr_approx,103)
!			flush(101)
!			flush(102)
!			flush(103)
!			Close(101)
!			Close(102)
!			Close(103)
!			i3=i3+1
!		End Do
!		Write(201,*) phi, Ta, gam_val, ua_writting
!		flush(201)
!	End Do
!	Close(201)
!End Subroutine

End Module


text:
!
!     (c) 2020-2021 Sourcery, Inc.
!     This software was developed for the U.S. Nuclear Regulatory Commission (US NRC) under contract # 31310020D0006:
!     "Technical Assistance in Support of NRC Nuclear Regulatory Research for Materials, Waste, and Reactor Programs"
!
module output_data_m
  use input_data_m, only : input_data_t
  use oracle_interface, only : oracle_t
  use random_samples_m, only: random_samples_t
  implicit none

  private

  type, public, extends(oracle_t) :: output_data_t
    private
    type(input_data_t) input_data_
    real ,allocatable :: R_Tndt_(:)
    real ,allocatable :: K_hist_(:), Chemistry_content_(:,:), Chemistry_factor_(:)
    real ,allocatable :: CPI_(:)
    real ,allocatable :: CPI_avg_(:)
  contains
    procedure :: assign
    generic :: assignment(=) => assign
    procedure :: norm
    procedure :: subtract
    procedure :: write_formatted
    procedure :: R_Tndt
    procedure :: K_hist
    procedure :: Chemistry_content
    procedure :: Chemistry_factor
    procedure :: CPI
    procedure :: CPI_avg
    procedure :: nsim
    procedure :: ntime
  end type

  interface output_data_t

    pure module function default_constructor() result(new_output_data_t)
      implicit none
      type(output_data_t) new_output_data_t
    end function

    module function whole_shebang(input_data, random_samples) result(new_output_data)
      implicit none
      type(input_data_t), intent(in) :: input_data
      type(random_samples_t), intent(in) :: random_samples(:)
      type(output_data_t) :: new_output_data
    end function

    pure module function new_output_data(input_data, R_Tndt, K_hist, Chemistry_content, Chemistry_factor, CPI, CPI_avg)
      implicit none
      type(input_data_t), intent(in) :: input_data
      real, intent(in) :: R_Tndt(:)
      real, intent(in) :: K_hist(:), Chemistry_content(:,:), Chemistry_factor(:)
      real, intent(in) :: CPI(:)
      real, intent(in) :: CPI_avg(:)
      type(output_data_t) new_output_data
    end function

  end interface

  interface

    module subroutine assign(self, rhs)
      implicit none
      class(output_data_t), intent(inout) :: self
      class(oracle_t), intent(in) :: rhs
    end subroutine

    module function subtract(self, rhs) result(difference)
      !! result has components corresponding to subtracting rhs's components fron self object's components
      implicit none
      class(output_data_t), intent(in) :: self
      class(oracle_t), intent(in) :: rhs
      class(oracle_t), allocatable :: difference
    end function

    pure module function norm(self) result(norm_of_self)
      !! result is a norm of the array formed by concatenating the real components of self object
      implicit none
      class(output_data_t), intent(in) :: self
      real norm_of_self
    end function

    module subroutine write_formatted(self, unit, iotype, v_list, iostat, iomsg)
      implicit none
      class(output_data_t), intent(in) :: self
      integer, intent(in) :: unit
      character(*), intent(in) :: iotype
      integer, intent(in) :: v_list(:)
      integer, intent(out) :: iostat
      character(*), intent(inout) :: iomsg
    end subroutine

    module function R_Tndt(self) result(my_R_Tndt)
      implicit none
      class(output_data_t), intent(in) :: self
      real, allocatable :: my_R_Tndt(:)
    end function

    module function K_hist(self) result(my_K_hist)
      implicit none
      class(output_data_t), intent(in) :: self
      real, allocatable :: my_K_hist(:)
    end function

    module function Chemistry_content(self) result(my_Chemistry_content)
      implicit none
      class(output_data_t), intent(in) :: self
      real, allocatable :: my_Chemistry_content(:,:)
    end function

    module function Chemistry_factor(self) result(my_Chemistry_factor)
      implicit none
      class(output_data_t), intent(in) :: self
       real, allocatable :: my_Chemistry_factor(:)
    end function

    module function CPI(self) result(my_CPI)
      implicit none
      class(output_data_t), intent(in) :: self
      real ,allocatable :: my_CPI(:)
    end function

    module function CPI_avg(self) result(my_CPI_avg)
      implicit none
      class(output_data_t), intent(in) :: self
      real ,allocatable :: my_CPI_avg(:)
    end function

    module function nsim(self) result(my_nsim)
      implicit none
      class(output_data_t), intent(in) :: self
      integer my_nsim
    end function

    module function ntime(self) result(my_ntime)
      implicit none
      class(output_data_t), intent(in) :: self
      integer my_ntime
    end function

  end interface

end module output_data_m

max_stars_repo_name: sourceryinstitute/miniFAVOR

max_issues_repo_name: sourceryinstitute/miniFAVOR

max_forks_repo_name: sourceryinstitute/miniFAVOR


text:
Describe Users/William42 here.

20100921 16:28:08 nbsp WilliamWelcome to the Wiki! Sorry you had a bad experience with Acadian...not sure if posting daily is necessary as they have a clear long history on that page of treating people badly. Have you considered filing your complaints in the areas listed early on that page: http://www.yolocounty.org/index.aspx?recordid40&page632   http://daviswiki.org/Better_Business_Bureau ???? Users/PeterBoulay

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	PROGRAM NAMSND
C************************************************************************
C* PROGRAM NAMSND							*
C*									*
C* This programs decodes profile and surface data from NMC models and	*
C* creates both a surface and an upper-air GEMPAK sounding dataset.	*
C* The parameter packing files must be provided by the user.		*
C**									*
C* Log:									*
C* K. Brill/NMC		12/93						*
C* K. Brill/NMC		 7/94	Added SNBPRM				*
C* K. Brill/NMC		10/94	Changes for auxiliary sfc file		*
C* K. Brill/NMC		 1/95	Skip station for iret=2 in NAMRDD	*
C* K. Brill/NMC		 8/95	Add ASCII output option			*
C* K. Brill/EMC		 4/96	Add sfcfc2 for parms in aux file	*
C* K. Brill/EMC		 7/98	Changes for J. Woollen BUFR decoder	*
C* K. Brill/EMC		 9/98	Close BUFR file only if open		*
C* D. Kidwell/NCEP	11/98	snefil -> snbufr                        *
C* D. Kidwell/NCEP	12/98	SNMODL -> NAMSND, SNM... -> NAM...      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER	snbufr*132, snoutf*132, sfoutf*132, snprmf*132,
     +			sfprmf*132, timstn*32
	CHARACTER	dattim*20, dtmlst*20, stid*8, sidlst*8,
     +			jensid*3
	CHARACTER*4	parmsn (MMPARM), parmsf (MMPARM),
     +			prmsf2 (MMPARM), snbprm (MMPARM)
	REAL		sndata (LLMXLV*MMPARM), sfdata (2*MMPARM)
	REAL		sfdat2 (MMPARM)
	REAL		sncfac (MMPARM), sfcfac (MMPARM),
     +			sfcfc2 (MMPARM)
	CHARACTER*8	tblde (MMPARM)
	INTEGER		lunsf (2)
	LOGICAL		degree, done, proces, clsbfr
C------------------------------------------------------------------------
C*	Initialize TAE.
C
	CALL IP_INIT  ( respnd, ier )
	IF  ( ier .ne. 0 )  THEN
	    iret = -1
	    CALL ER_WMSG  ( 'NAMSND', iret, ' ', ier )
	    CALL SS_EXIT
	END IF
	CALL IP_IDNT  ( 'NAMSND', ier )
	done = .false.
	clsbfr = .false.
C
C*	Loop through program until user exits.
C
	DO WHILE  ( .not. done )
C
C*	    Get input.
C
	    CALL NAMINP  ( snbufr, snoutf, sfoutf, snprmf, sfprmf,
     +			   timstn, iret )
	    IF  ( iret .ne. 0 )  THEN
		CALL ER_WMSG  ( 'NAMSND', iret, ' ', ier )
		CALL SS_EXIT
	    END IF
	    proces = .true.
C
C*	    Open the model sounding file, process the parameter
C*	    files and open the output files.  Also read in the
C*	    conversion factors for the parameters.
C
	    CALL NAMOPN ( snbufr, snoutf, sfoutf, snprmf, sfprmf,
     +			  timstn, lunin, lunsn, lunsf,
     +			  jenlun, jenstn, jensid,
     +			  nprmsn, nprmsf, sncfac, sfcfac, parmsn,
     +			  parmsf, snbprm, npbfr, nprms2, prmsf2, sfcfc2,
     +			  tblde, nde, degree, clsbfr, iret )
	    IF  ( iret .ne. 0 )  THEN
		CALL ER_WMSG  ( 'NAMSND', iret, snbufr, ier )
		proces = .false.
	    END IF
C
C*	    Read the input file; write to the output file.
C
	    dtmlst = ' '
	    sidlst = ' '
	    isnlst = 0
	    istbyt = 1
	    DO WHILE ( proces )
		CALL NAMRDD ( lunin, tblde, nde, degree,
     +			      sncfac, sfcfac, nprmsn,
     +			      nprmsf, parmsn, parmsf, snbprm, npbfr,
     +			      nprms2, prmsf2, sfcfc2,
     +			      istbyt, dattim, stid, istnm,
     +			      slat, slon, selv, sndata, sfdata,
     +			      sfdat2, iclass, nz, ifhr, iret )
		IF ( iret .lt. 0 .and. iret .ne. -27 ) THEN
		    proces = .false.
		    CALL ER_WMSG ( 'NAMSND', iret, ' ', ier )
		ELSE IF ( iret .eq. 1 ) THEN
		    proces = .false.
		    CALL ER_WMSG ( 'NAMSND', iret, ' ', ier )
		ELSE IF ( iret .eq. -27 ) THEN
C
C*		    Just skip this station.
C
		    CALL ER_WMSG ( 'NAMSND', iret, 'WARNING', ier )
		END IF
C*
		IF ( iret .eq. 0 .and. jenlun .ne. 0
     +			.and. istnm .eq. jenstn ) THEN
C
C*		    Add to the ASCII file.
C
		    CALL NAMJEN ( jenlun, nprmsn, nprmsf, parmsn,
     +				  parmsf, dattim, jensid, jenstn,
     +				  slat, slon, selv, sndata, sfdata,
     +				  nz, ifhr, iret )
		END IF
C*
		IF ( iret .eq. 0 .and. lunsn .ne. 0 ) THEN
C
C*		    Set station and time for writing.
C
		    CALL NAMSST ( lunsn, lunsf, stid, istnm,
     +			          slat, slon, selv, dattim,
     +				  dtmlst, sidlst, isnlst, ierx )
		    IF ( ierx .eq. 0 .and. lunsn .ne. 0 ) THEN
     			CALL SN_WDAT ( lunsn, ifhr, nz, sndata,
     +     			       ier )
			IF ( ier .ne. 0 )
     +			    CALL ER_WMSG ( 'SN', ier, ' ', ierr )
		    END IF
		    IF ( ierx .eq. 0 .and. lunsf (1) .ne. 0 ) THEN
			CALL SF_WDAT ( lunsf (1), ifhr, sfdata, ier )
		        IF ( ier .ne. 0 )
     +			    CALL ER_WMSG ( 'SF', ier, ' ', ierr )
		    END IF
		    IF ( ierx .eq. 0 .and. lunsf (2) .ne. 0 ) THEN
			CALL SF_WDAT ( lunsf (2), ifhr, sfdat2, ier )
		        IF ( ier .ne. 0 )
     +			    CALL ER_WMSG ( 'SF', ier, ' ', ierr )
		    END IF
		END IF
	    END DO
C
C*	    Close files.
C
	    IF ( clsbfr ) CALL JB_CLOS  ( lunin, ier )
	    CALL SN_CLOS  ( lunsn, ier )
	    CALL SF_CLOS  ( lunsf(1), ier )
	    CALL SF_CLOS  ( lunsf(2), ier )
C
C*	    Call dynamic tutor.
C
	    CALL IP_DYNM  ( done, ier )
C*
	END DO
C
C*	Exit.
C
	CALL IP_EXIT  ( ier )
C*
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do compile }
! { dg-options "-fcoarray=lib" }

program p
   call s
contains
   subroutine s
      real :: x[*] = 1
      block
      end block
      x = 2
   end
end


text:
program MHD2D

!A 'CLASSICAL HYDRODYNAMIC' SOLVER (WAVES PROPAGATED THRU SOURCE TERMS)
!FOR THE IDEAL MHD EQUATIONS IN 2 DIMENSIONS.  M. ZETTERGREN, T. SYMONS, K. CHAN
!
! LOOP LEVEL PARALLELIZATION HAS BEEN IMPLEMENTED USING OPENMP

use phys_consts
use numerical
use MHD

implicit none

integer, parameter :: outunit=42,method=1
!integer, parameter :: npts1=1001, npts2=1001
!integer, parameter :: npts1=501, npts2=501
integer, parameter :: npts1=251, npts2=251
real(8), parameter :: tcfl=0.75d0

!real(8), parameter :: stride1=2d3,stride2=2d3
!real(8), parameter :: stride1=4d3,stride2=4d3
real(8), parameter :: stride1=8d3,stride2=8d3

real(8), dimension(npts1) :: dx1i
real(8), dimension(-1:npts1+2) :: x1
real(8), dimension(npts1+1) :: x1i
real(8), dimension(0:npts1+2) :: dx1

real(8), dimension(npts2) :: dx2i
real(8), dimension(-1:npts2+2) :: x2
real(8), dimension(npts2+1) :: x2i
real(8), dimension(0:npts2+2) :: dx2

integer :: it,lx1,lx2,ix1,ix2

real(8), dimension(-1:npts1+2,-1:npts2+2) :: rhom,rhou1,rhou2,rhou3,B1,B2,B3,rhoeps,u1,u2,u3

real(8), dimension(npts1,npts2) :: du1full,du2full,vcon1,vcon2
real(8), dimension(npts1,npts2) :: p,vA,vsnd,vmatter

real(8) :: t=0d0,dt=1d-6,dtout=5d0,tdur=120d0
real(8) :: dt1,dt2
real(8) :: toutnext

real(8) :: tstart,tfin                           !for assessing performance of code blocks

real(8) :: Ti=8000d0                             !temperature used to initialize internal energy
real(8) :: vsw=15d3*5d0                          !initial condition solar wind speed

!GRID CONSTRUCTION
write(*,*) 'Creating grid of size:  ',npts1,npts2
x1=[ (ix1*stride1, ix1=-2,npts1+1) ]
lx1=size(x1)-4                             !exclude ghost cells in count
x1=x1-0.5d0*(x1(1)+x1(lx1))    !center the grid
dx1=x1(0:lx1+2)-x1(-1:lx1+1)                !backward diffs
x1i(1:lx1+1)=0.5d0*(x1(0:lx1)+x1(1:lx1+1))    !interface data
dx1i=x1i(2:lx1+1)-x1i(1:lx1)

x2=[ (ix2*stride2, ix2=-2,npts2+1) ]
lx2=size(x2)-4                             !exclude ghost cells in count
x2=x2-0.5d0*(x2(1)+x2(lx2))    !center
dx2=x2(0:lx2+2)-x2(-1:lx2+1)                !backward diffs
x2i(1:lx2+1)=0.5d0*(x2(0:lx2)+x2(1:lx2+1))    !interface data
dx2i=x2i(2:lx1+1)-x2i(1:lx2)

!PREP OUTPUT FILE (BINARY STREAM), WRITE GRID SIZES AND GRID
write(*,*) 'Opening output file...'
open(outunit,file='MHD2D.dat',status='replace',access='stream')
write(outunit) lx1
write(outunit) x1
write(outunit) lx2
write(outunit) x2

!INITIAL CONDITIONS
write(*,*) 'Setting initial conditions...'
rhom=2d9*amu
rhou1=0d0
rhou2=2d9*amu*vsw
!rhou2=0d0
rhou3=0d0
B1=0d0
B2=0d0
B3=0d0
rhoeps=rhom*kb*Ti/amu/(gammap-1)

!SET VELOCITIES
u1=rhou1/rhom
u2=rhou2/rhom
u3=rhou3/rhom

!MAIN LOOP
toutnext=0d0
do while (t<tdur)
  call cpu_time(tstart)

  !VELOCITY DIFFERENCES FOR CALCULATING ARTIFICAL VISCOSITY
  call visc_vel_diffs(u1,u2,du1full,du2full)

  !TIME STEP DETERMINATION (A BIT UGLY...)
  !$omp parallel do
  do ix2=1,lx2
    vA(:,ix2)=sqrt((B1(1:lx1,ix2)**2+B2(1:lx1,ix2)**2+B3(1:lx1,ix2)**2)/mu0/rhom(1:lx1,ix2))
    p(:,ix2)=rhoeps(1:lx1,ix2)*(gammap-1d0)
    vsnd(:,ix2)=sqrt(gammap*p(:,ix2)/rhom(1:lx1,ix2))
    vmatter(:,ix2)=sqrt(u1(1:lx1,ix2)**2+u2(1:lx1,ix2)**2)
    vcon1(:,ix2)=abs(0.5d0*xi**2*min(du1full(:,ix2),0d0))
    vcon2(:,ix2)=abs(0.5d0*xi**2*min(du2full(:,ix2),0d0))
  end do

  dt1=tcfl*minval(dx1i/(maxval(vA,2)+ &
                      maxval(vsnd,2)+ &
                      maxval(vcon1,2)+ &
                      maxval(abs(u1(1:lx1,1:lx2)),2)))
  dt2=tcfl*minval(dx2i/(maxval(vA,1)+ &
                      maxval(vsnd,1)+ &
                      maxval(vcon2,1)+ &
                      maxval(abs(u2(1:lx1,1:lx2)),1)))
  dt=min(dt1,dt2)
  t=t+dt                                        !time after this step
  write(*,*) 'Time:  ',t,'  Time step:  ',dt
  write(*,*) '  Alfven speed:  ',maxval(pack(vA,.true.))/1d3
  write(*,*) '  Sound speed:  ',maxval(pack(vsnd,.true.))/1d3
  write(*,*) '  Matter speed:  ',maxval(pack(vmatter,.true.))/1d3
  write(*,*) '    viscous-1 speed:  ',maxval(pack(vcon1,.true.))/1d3, &
             ' viscous-2 speed:  ',maxval(pack(vcon2,.true.))/1d3

  !FILL GHOST CELLS
  call boundaries(t,x1,x2,rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3)

  !PICK A REGION IN WHICH TO INTRODUCE AN OBJECT INTO THE FLOW
  call flow_obs(t,x1,x2,rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3)

  !ADVECTION SUBSTEP
  call MHD_advection(dt,dx1,dx1i,dx2,dx2i,rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3)

  !SOURCES SUBSTEP
  call MHD_sources(dt,dx1,dx2,du1full,du2full,rhom,rhou1,rhou2,rhou3,rhoeps, &
                   B1,B2,B3,u1,u2,u3)

  !OUTPUT (RAW BINARY STREAM)
  if (t>toutnext) then
    write(*,*) 'Initiating output for time step:  ',t
    write(outunit) t
    write(outunit) rhom/amu
    write(outunit) u1
    write(outunit) u2
    write(outunit) u3
    write(outunit) B1
    write(outunit) B2
    write(outunit) B3
    write(outunit) rhoeps*(gammap-1)
    toutnext=toutnext+dtout
  end if

  call cpu_time(tfin)
  write(*,*) 'Time step completed in cpu time of:  ',tfin-tstart
end do

close(outunit)

contains

  subroutine boundaries(t,x1,x2,rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3)

    !------------------------------------------------------------
    !-------FILL THE GHOST CELLS TO ENFORCE BOUNDARY CONDITIONS
    !------------------------------------------------------------

    real(8), intent(in) :: t
    real(8), dimension(-1:), intent(in) :: x1
    real(8), dimension(-1:), intent(in) :: x2
    real(8), dimension(-1:,-1:), intent(inout) :: rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3

    real(8) :: vsw

    !SIZES WITHOUT GHOST CELLS
    lx1=size(rhom,1)-4; lx2=size(rhom,2)-4;

    !SET THE SOLAR WIND SPEED AT BOUNDARY (COULD BE TIME-DEPENDENT)
    vsw=15d3*5d0

    !BOUNDARY CONDITIONS (THESE ARE FREE-FLOW)
    rhom(-1:0,1:lx2)=spread(rhom(1,1:lx2),1,2)              !bottom
    rhom(lx1+1:lx1+2,1:lx2)=spread(rhom(1,1:lx2),1,2)       !top
    rhom(1:lx1,-1:0)=2d9*amu                                !left
    rhom(1:lx1,lx2+1:lx2+2)=spread(rhom(1:lx1,lx2),2,2)     !right

    rhou1(-1:0,1:lx2)=spread(rhou1(1,1:lx2),1,2)
    rhou1(lx1+1:lx1+2,1:lx2)=spread(rhou1(lx1,1:lx2),1,2)
    rhou1(1:lx1,-1:0)=spread(rhou1(1:lx1,1),2,2)
    rhou1(1:lx1,lx2+1:lx2+2)=spread(rhou1(1:lx1,lx2),2,2)

    rhou2(-1:0,1:lx2)=spread(rhou2(1,1:lx2),1,2)
    rhou2(lx1+1:lx1+2,1:lx2)=spread(rhou2(1,1:lx2),1,2)
    rhou2(1:lx1,-1:0)=2d9*amu*vsw
    rhou2(1:lx1,lx2+1:lx2+2)=spread(rhou2(1:lx1,lx2),2,2)

    rhou3(-1:0,1:lx2)=spread(rhou3(1,1:lx2),1,2)
    rhou3(lx1+1:lx1+2,1:lx2)=spread(rhou3(lx1,1:lx2),1,2)
    rhou3(1:lx1,-1:0)=spread(rhou3(1:lx1,1),2,2)
    rhou3(1:lx1,lx2+1:lx2+2)=spread(rhou3(1:lx1,lx2),2,2)

    B1(-1:0,1:lx2)=spread(B1(1,1:lx2),1,2)
    B1(lx1+1:lx1+2,1:lx2)=spread(B1(lx1,1:lx2),1,2)
    B1(1:lx1,-1:0)=spread(B1(1:lx1,1),2,2)
!    B1(1:lx1,-1:0)=1d-9
    B1(1:lx1,lx2+1:lx2+2)=spread(B1(1:lx1,lx2),2,2)

    B2(-1:0,1:lx2)=spread(B2(1,1:lx2),1,2)
    B2(lx1+1:lx1+2,1:lx2)=spread(B2(lx1,1:lx2),1,2)
    B2(1:lx1,-1:0)=spread(B2(1:lx1,1),2,2)
    B2(1:lx1,lx2+1:lx2+2)=spread(B2(1:lx1,lx2),2,2)

    B3(-1:0,1:lx2)=spread(B3(1,1:lx2),1,2)
    B3(lx1+1:lx1+2,1:lx2)=spread(B3(lx1,1:lx2),1,2)
!    B3(1:lx1,-1:0)=spread(B3(1:lx1,1),2,2)
    B3(1:lx1,-1:0)=1d-9    !works well
    B3(1:lx1,lx2+1:lx2+2)=spread(B3(1:lx1,lx2),2,2)

    rhoeps(-1:0,1:lx2)=spread(rhoeps(1,1:lx2),1,2)
    rhoeps(lx1+1:lx1+2,1:lx2)=spread(rhoeps(lx1,1:lx2),1,2)
    rhoeps(1:lx1,-1:0)=spread(rhoeps(1:lx1,1),2,2)
    rhoeps(1:lx1,lx2+1:lx2+2)=spread(rhoeps(1:lx1,lx2),2,2)

  end subroutine boundaries

  subroutine flow_obs(t,x1,x2,rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3)

    !------------------------------------------------------------
    !-------CREATES AN OBSTRUCTION TO FLOW ON THE GRID.  HAVE
    !-------THIS SUBPROGRAM DO NOTHING FOR A `FREE' INTERIOR.
    !------------------------------------------------------------

    real(8), intent(in) :: t
    real(8), dimension(-1:), intent(in) :: x1
    real(8), dimension(-1:), intent(in) :: x2
    real(8), dimension(-1:,-1:), intent(inout) :: rhom,rhou1,rhou2,rhou3,rhoeps,B1,B2,B3

    real(8) :: x1ctr,x2ctr,dist,objrad               !for solid object embedded in flow

    !SIZES WITHOUT GHOST CELLS
    lx1=size(rhom,1)-4; lx2=size(rhom,2)-4;

    !PICK A CIRCULAR REGION IN WHICH TO ZERO OUT VELOC.
    x1ctr=0d0
    x2ctr=0d0
    objrad=25d0*stride1/2d0
    do ix2=1,lx2
      do ix1=1,lx1
        dist=sqrt(x1(ix1)**2+x2(ix2)**2)
        if (dist<=objrad) then
          rhou1(ix1,ix2)=0d0
          rhou2(ix1,ix2)=0d0
          rhom(ix1,ix2)=2d9*amu
          rhoeps(ix1,ix2)=rhom(ix1,ix2)*kb*Ti/amu/(gammap-1)
          B1(ix1,ix2)=0d0
          B2(ix1,ix2)=0d0
          B3(ix1,ix2)=0d0
        end if
      end do
    end do

  end subroutine flow_obs

end program MHD2D

max_stars_repo_path: archive/MHD2D copy.f90

max_issues_repo_path: archive/MHD2D copy.f90

max_forks_repo_path: archive/MHD2D copy.f90


text:
      program lenddo

      real x(100), y(100)

      do i = 1, 100
         if(y(j).lt.0.) go to 1
         x(i) = 0.
 1    enddo

      do 2 i = 1, 99
         x(i) = 1.
 2    enddo

      do i = 1, 98
         x(i) = 2.
      enddo

c      do i = 1, 100
c         x(i) = 3.
c 4    continue

      do 5 j = 1, 99
         do 5 i = 1, 99
            x(i) = 1.+real(j)
 5    enddo

      end


text:
c
c program markibis
c
      include 'VICMAIN_FOR'
      subroutine main44

      parameter (ntable=200000)
c ntable = max # tiepoints
      integer*4 status,inmark,outmark,inibis,outibis
      real*4 buf(4,ntable)
      character*80 msg
      logical flag

c determine the type of the input file
      call xvunit(i,'INP',1,status,' ')
      call ibis_file_open(i,j,'READ',0,0,' ',' ',status)
      if(status.eq.1)then
        mode=1 ! ibis input
      else
        mode=2 ! mark input
      endif
      call ibis_file_close(j,' ',status)
      if(mode.eq.1) goto 200

c*********************************************************************
c load the tiepoint table from the mark file.
      call xvunit(inmark,'INP',1,status,' ')
      call xvopen(inmark,status,'IO_ACT','AS','U_FORMAT','REAL'
     +            ,'OPEN_ACT','AS',' ')
      call xvget(inmark,status,'NL',nl,'NS',ns,' ')
      npts=nl*ns/4
      if(npts.gt.ntable)then
        call xvmessage('Too many input points',' ')
        call abend
      endif
      k=1
      do i=1,nl
         call xvread(inmark,buf(1,k),status,'LINE',i,' ')
         k=k+ns/4
      enddo
      write(msg,*)nl,' mark tiepoint records read'
      call xvmessage(msg,' ')

c remove zeroes
      j=0
10    j=j+1
11      flag=.false.
        if(buf(1,j).eq.0.0) flag=.true.
        if(buf(2,j).eq.0.0) flag=.true.
        if(buf(3,j).eq.0.0) flag=.true.
        if(buf(4,j).eq.0.0) flag=.true.
        if(flag)then
          if(j.eq.npts)then ! last point is bad
            npts=npts-1
            goto 12
          endif
          do i=1,npts-j     ! move everything up one
            buf(1,j+i-1)=buf(1,j+i)
            buf(2,j+i-1)=buf(2,j+i)
            buf(3,j+i-1)=buf(3,j+i)
            buf(4,j+i-1)=buf(4,j+i)
          enddo
          npts=npts-1
          goto 11
        endif
      if(j.lt.npts) goto 10
12    continue
      write(msg,*)npts,' good tiepoints located'
      call xvmessage(msg,' ')

c switch the tiepoint locations so IN is columns 3,4 and OUT is 1,2
c which is the ibis tiepoint file convention.
      do i=1,npts
        tl=buf(1,i)
        ts=buf(2,i)
        buf(1,i)=buf(3,i)
        buf(2,i)=buf(4,i)
        buf(3,i)=tl
        buf(4,i)=ts
      enddo

c write the ibis file
      call xvunit(outibis,'OUT',1,status,' ')
      call iwrite_tiepoints(outibis,0,0,npts,buf,4)
      call xvmessage('Ibis file written',' ')
      return

c**********************************************************************
c read the ibis file
200   continue
      call xvunit(inibis,'INP',1,status,' ')
      i=0
      call iread_tiepoints(inibis,i,npts,ntable,buf,4)
      write(msg,*)npts,' tiepoints located'
      call xvmessage(msg,' ')

c fill the end of the buffer with zeroes
      nl=npts/32
      if(nl*32.lt.npts) nl=nl+1
      if(nl*32.gt.npts)then
        do i=npts+1,nl*32
          buf(1,i)=0.0
          buf(2,i)=0.0
          buf(3,i)=0.0
          buf(4,i)=0.0
        enddo
      endif

c switch the tiepoint locations so IN is columns 1,2 and OUT is 3,4
c which is the mark convention.
      do i=1,npts
        tl=buf(1,i)
        ts=buf(2,i)
        buf(1,i)=buf(3,i)
        buf(2,i)=buf(4,i)
        buf(3,i)=tl
        buf(4,i)=ts
      enddo

c write the mark file.
      call xvunit(outmark,'OUT',1,status,' ')
      call xvopen(outmark,status,'IO_ACT','AS','U_FORMAT','REAL',
     +            'U_NL',nl,'U_NS',128,'U_NB',1,
     +            'OP','WRITE','O_FORMAT','REAL','OPEN_ACT','AS',' ')
      do i=1,nl
         call xvwrit(outmark,buf(1,(i-1)*32+1),status,'LINE',i,' ')
      enddo
      call xvclose(outmark,status,'CLOS_ACT','FREE',' ')
      call xvmessage('Mark file written',' ')

      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


! shows how to properly specify character length and character array

character*(5) :: cb5  ! obsolete, don't use
character(5)  :: c5   ! good to use
character(*),parameter  :: ca(2)='hello'

if (.not.all([len(cb5)==5,len(c5)==5,len(ca)==5,size(ca)==2])) error stop
print *, ca

end program


text:
! RUN: bbc -fopenmp -emit-fir %s -o - | \
! RUN: FileCheck %s --check-prefix=FIRDialect
! RUN: bbc -fopenmp %s -o - | fir-opt --fir-to-llvm-ir | \
! RUN: FileCheck %s --check-prefix=LLVMDialect

! This test checks the lowering of atomic read

!FIRDialect: func @_QQmain() {
!FIRDialect: %[[VAR_A:.*]] = fir.address_of(@_QFEa) : !fir.ref<!fir.char<1>>
!FIRDialect: %[[VAR_B:.*]] = fir.address_of(@_QFEb) : !fir.ref<!fir.char<1>>
!FIRDialect: %[[VAR_C:.*]] = fir.alloca !fir.logical<4> {bindc_name = "c", uniq_name = "_QFEc"}
!FIRDialect: %[[VAR_D:.*]] = fir.alloca !fir.logical<4> {bindc_name = "d", uniq_name = "_QFEd"}
!FIRDialect: %[[VAR_E:.*]] = fir.address_of(@_QFEe) : !fir.ref<!fir.char<1,8>>
!FIRDialect: %[[VAR_F:.*]] = fir.address_of(@_QFEf) : !fir.ref<!fir.char<1,8>>
!FIRDialect: %[[VAR_G:.*]] = fir.alloca f32 {bindc_name = "g", uniq_name = "_QFEg"}
!FIRDialect: %[[VAR_H:.*]] = fir.alloca f32 {bindc_name = "h", uniq_name = "_QFEh"}
!FIRDialect: %[[VAR_X:.*]] = fir.alloca i32 {bindc_name = "x", uniq_name = "_QFEx"}
!FIRDialect: %[[VAR_Y:.*]] = fir.alloca i32 {bindc_name = "y", uniq_name = "_QFEy"}
!FIRDialect: omp.atomic.read %[[VAR_X]] = %[[VAR_Y]] memory_order(acquire)  hint(uncontended) : !fir.ref<i32>
!FIRDialect: omp.atomic.read %[[VAR_A]] = %[[VAR_B]] memory_order(relaxed) hint(none)  : !fir.ref<!fir.char<1>>
!FIRDialect: omp.atomic.read %[[VAR_C]] = %[[VAR_D]] memory_order(seq_cst)  hint(contended) : !fir.ref<!fir.logical<4>>
!FIRDialect: omp.atomic.read %[[VAR_E]] = %[[VAR_F]] hint(speculative) : !fir.ref<!fir.char<1,8>>
!FIRDialect: omp.atomic.read %[[VAR_G]] = %[[VAR_H]] hint(nonspeculative) : !fir.ref<f32>
!FIRDialect: omp.atomic.read %[[VAR_G]] = %[[VAR_H]] : !fir.ref<f32>
!FIRDialect: return
!FIRDialect: }

!LLVMDialect: llvm.func @_QQmain() {
!LLVMDialect: %[[LLVM_VAR_A:.*]] = llvm.mlir.addressof @_QFEa : !llvm.ptr<array<1 x i8>>
!LLVMDialect: %[[LLVM_VAR_B:.*]] = llvm.mlir.addressof @_QFEb : !llvm.ptr<array<1 x i8>>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect:  %[[LLVM_VAR_C:.*]] = llvm.alloca {{.*}} x i32 {bindc_name = "c", in_type = !fir.logical<4>, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEc"} : (i64) -> !llvm.ptr<i32>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect: %[[LLVM_VAR_D:.*]] = llvm.alloca {{.*}} x i32 {bindc_name = "d", in_type = !fir.logical<4>, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEd"} : (i64) -> !llvm.ptr<i32>
!LLVMDialect: %[[LLVM_VAR_E:.*]] = llvm.mlir.addressof @_QFEe : !llvm.ptr<array<8 x i8>>
!LLVMDialect: %[[LLVM_VAR_F:.*]] = llvm.mlir.addressof @_QFEf : !llvm.ptr<array<8 x i8>>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect: %[[LLVM_VAR_G:.*]] = llvm.alloca {{.*}} x f32 {bindc_name = "g", in_type = f32, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEg"} : (i64) -> !llvm.ptr<f32>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect: %[[LLVM_VAR_H:.*]] = llvm.alloca {{.*}} x f32 {bindc_name = "h", in_type = f32, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEh"} : (i64) -> !llvm.ptr<f32>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect: %[[LLVM_VAR_X:.*]] = llvm.alloca {{.*}} x i32 {bindc_name = "x", in_type = i32, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEx"} : (i64) -> !llvm.ptr<i32>
!LLVMDialect: {{.*}} = llvm.mlir.constant(1 : i64) : i64
!LLVMDialect: %[[LLVM_VAR_Y:.*]] = llvm.alloca {{.*}} x i32 {bindc_name = "y", in_type = i32, operand_segment_sizes = dense<0> : vector<2xi32>, uniq_name = "_QFEy"} : (i64) -> !llvm.ptr<i32>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_X]] = %[[LLVM_VAR_Y]] memory_order(acquire)  hint(uncontended) : !llvm.ptr<i32>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_A]] =  %[[LLVM_VAR_B]] memory_order(relaxed) hint(none) : !llvm.ptr<array<1 x i8>>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_C]] = %[[LLVM_VAR_D]] memory_order(seq_cst)  hint(contended) : !llvm.ptr<i32>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_E]] = %[[LLVM_VAR_F]] hint(speculative) : !llvm.ptr<array<8 x i8>>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_G]] = %[[LLVM_VAR_H]] hint(nonspeculative) : !llvm.ptr<f32>
!LLVMDialect: omp.atomic.read %[[LLVM_VAR_G]] = %[[LLVM_VAR_H]] : !llvm.ptr<f32>
!LLVMDialect: llvm.return
!LLVMDialect: }

program OmpAtomic

    use omp_lib
    integer :: x, y
    character :: a, b
    logical :: c, d
    character(8) :: e, f
    real g, h
    !$omp atomic acquire read hint(omp_sync_hint_uncontended)
       x = y
    !$omp atomic relaxed read hint(omp_sync_hint_none)
       a = b
    !$omp atomic read seq_cst hint(omp_sync_hint_contended)
       c = d
    !$omp atomic read hint(omp_sync_hint_speculative)
       e = f
    !$omp atomic read hint(omp_sync_hint_nonspeculative)
       g = h
    !$omp atomic read
       g = h
end program OmpAtomic

max_stars_repo_name: ornata/llvm-project

max_issues_repo_name: ornata/llvm-project

max_forks_repo_name: ornata/llvm-project


! Copyright 2005-2017 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!  Description: how to get values using keys.
!
!
program get
use grib_api
implicit none

  integer                                         ::  ifile
  integer                                         ::  iret
  integer                                         ::  igrib
  integer                                         ::  i
  real(kind=8)                                    ::  latitudeOfFirstPointInDegrees
  real(kind=8)                                    ::  longitudeOfFirstPointInDegrees
  real(kind=8)                                    ::  latitudeOfLastPointInDegrees
  real(kind=8)                                    ::  longitudeOfLastPointInDegrees
  real(kind=8)                                    ::  jDirectionIncrementInDegrees
  real(kind=8)                                    ::  iDirectionIncrementInDegrees
  integer(kind = 4)                               ::  numberOfPointsAlongAParallel
  integer(kind = 4)                               ::  numberOfPointsAlongAMeridian
  real(kind=8), dimension(:), allocatable         ::  values
  integer(kind = 4)                               ::  numberOfValues
  real(kind=8)                                    ::  average
  character(len = 256)                            ::  error

  call grib_open_file(ifile, &
       '../../data/reduced_latlon_surface.grib1','r')

!     a new grib message is loaded from file
!     igrib is the grib id to be used in subsequent calls
  call  grib_new_from_file(ifile,igrib)

!     get as a integer
  call grib_get(igrib,'numberOfPointsAlongAParallel', &
                               numberOfPointsAlongAParallel)
  write(*,*) 'numberOfPointsAlongAParallel=', &
              numberOfPointsAlongAParallel

!     get as a integer
  call grib_get(igrib,'numberOfPointsAlongAMeridian', &
                                       numberOfPointsAlongAMeridian)
  write(*,*) 'numberOfPointsAlongAMeridian=', &
              numberOfPointsAlongAMeridian

!     get as a real8
call grib_get(igrib, &
                      'latitudeOfFirstGridPointInDegrees', &
                       latitudeOfFirstPointInDegrees)
write(*,*) 'latitudeOfFirstGridPointInDegrees=', &
            latitudeOfFirstPointInDegrees

!     get as a real8
  call grib_get(igrib, &
               'longitudeOfFirstGridPointInDegrees', &
                longitudeOfFirstPointInDegrees)
  write(*,*) 'longitudeOfFirstGridPointInDegrees=', &
     longitudeOfFirstPointInDegrees

!     get as a real8
  call grib_get(igrib, &
     'latitudeOfLastGridPointInDegrees', &
     latitudeOfLastPointInDegrees)
  write(*,*) 'latitudeOfLastGridPointInDegrees=', &
     latitudeOfLastPointInDegrees

!     get as a real8
  call grib_get(igrib, &
     'longitudeOfLastGridPointInDegrees', &
      longitudeOfLastPointInDegrees)
  write(*,*) 'longitudeOfLastGridPointInDegrees=', &
              longitudeOfLastPointInDegrees

!     get as a real8
  call grib_get(igrib, &
                   'jDirectionIncrementInDegrees', &
                    jDirectionIncrementInDegrees)
  write(*,*) 'jDirectionIncrementInDegrees=', &
              jDirectionIncrementInDegrees

!     get as a real8
  call grib_get(igrib, &
       'iDirectionIncrementInDegrees', &
        iDirectionIncrementInDegrees)
  write(*,*) 'iDirectionIncrementInDegrees=', &
              iDirectionIncrementInDegrees

!     get the size of the values array
  call grib_get_size(igrib,'values',numberOfValues)
  write(*,*) 'numberOfValues=',numberOfValues

  allocate(values(2*numberOfValues), stat=iret)
!     get data values
  print*, size(values)
  call grib_get(igrib,'values',values)

  average = 0
  do i=1,numberOfValues
        average = average + values(i);
  enddo

  average =average / numberOfValues

  write(*,*)'There are ',numberOfValues, &
            ' average is ',average

  call grib_release(igrib)

  call grib_close_file(ifile)

  deallocate(values)
end program get


text:
	subroutine crgass(atm,res,rnum,chn,crg)
c
c actual charge assignment from the hash tables...
	include "qdiffpar4.h"
c
	character*3 sres
	character*1 schn
	character*4 snum
	character*6 satm
c
c	character*1 chn
c	character*3 res
c	character*4 rnum
c	character*6 atm
c
c
c assignchargeto atom, searching for decreasingly specific specification
c kas,29jn 89- add extra searches with wild card:
c search order:  atom, res, num, chain
c                 *              x      walter, modified to account for atom * wildcard
c				x     x    x   x
c		 		x     x    x
c				x     x        x
c				x     x
c				x          x   x
c				x          x
c				x              x
c				x
c
c note if no charge record found, is assumed to be 0.0
c
c b+++++++++++++++++++++wildcard for atom
      satm= '*     '
	sres = '   '
	snum = '    '
      call cfind(satm,sres,snum,chn,ifind,n)
      if(ifind.eq.0) then
c e++++++++++++++++++++++++++++++++++++++
	  call cfind(atm,res,rnum,chn,ifind,n)
	  if(ifind.eq.0) then
	    schn = chn
	    chn = ' '
	    call cfind(atm,res,rnum,chn,ifind,n)
	    if(ifind.eq.0) then
		chn = schn
		snum = rnum
		rnum = '    '
	      call cfind(atm,res,rnum,chn,ifind,n)
	      if(ifind.eq.0) then
	        schn = chn
	        chn = ' '
	        call cfind(atm,res,rnum,chn,ifind,n)
		  if(ifind.eq.0) then
		    chn = schn
		    rnum = snum
		    sres = res
		    res = '   '
	          call cfind(atm,res,rnum,chn,ifind,n)
	          if(ifind.eq.0) then
	            schn = chn
	            chn = ' '
	            call cfind(atm,res,rnum,chn,ifind,n)
	            if(ifind.eq.0) then
		        chn = schn
		        snum = rnum
		        rnum = '    '
	              call cfind(atm,res,rnum,chn,ifind,n)
	              if(ifind.eq.0) then
	                schn = chn
	                chn = ' '
	                call cfind(atm,res,rnum,chn,ifind,n)
			  endif
			endif
		    endif
		  end if
		end if
	    end if
	  end if
	end if
c
	if(ifind.ne.0) then
	crg=chrgvt(n)
	else
	crg=0.0
	end if
c
	return
	end

max_stars_repo_name: caixiuhong/Stable-MCCE

max_issues_repo_name: caixiuhong/Stable-MCCE

max_forks_repo_name: caixiuhong/Stable-MCCE


text:
module mpi_f08_types_10_75

   use, intrinsic :: ISO_C_BINDING

   integer(C_INT), parameter :: MPI_STATUS_SIZE = 5

 ! type, BIND(C) :: MPI_Comm
   type, PUBLIC, BIND(C) :: MPI_Comm
      integer(C_INT) :: MPI_VAL
   end type MPI_Comm

 ! This is the problem statement.
   type(MPI_Comm), protected, bind(C, name="ompi_f08_mpi_comm_world") :: MPI_COMM_WORLD
 ! type(MPI_Comm), protected, bind(C) :: MPI_COMM_WORLD
 ! type(MPI_Comm), public, bind(C, name="ompi_f08_mpi_comm_world") :: MPI_COMM_WORLD

end module mpi_f08_types_10_75

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!=======================================================================
C  MODULE FertType_mod
C  06/15/2014 CHP Written
!=======================================================================

      MODULE FertType_mod
!     Contains data definitions for flexible fertilizer inputs
      USE ModuleDefs

!     Data construct for control variables
      TYPE FertPropType
        REAL NO3_N_pct, NH4_N_pct, UREA_N_pct
        REAL UIEFF, UIDUR, NIEFF, NIDUR, NRL50, NSIGK
        CHARACTER*6 FertN, FertP, FertK, FertC, FertMg, FertS
        CHARACTER*3 NRFNC
        CHARACTER*35 FertName
        INTEGER Check
      END TYPE FertPropType

      INTEGER, PARAMETER :: NFertTypes = 999
      Type (FertPropType), dimension(NFertTypes), protected :: FertFile
      INTEGER Number_of_Fertilizers

      TYPE SlowReleaseNType
        INTEGER StartYRDOY
        INTEGER KMAX
        REAL N0, NRL50, KN, AN, BN, NO3_frac, NH4_frac, UREA_frac
        REAL FMIXEFF, PROF(NL)
        REAL CumRelToday, CumRelYesterday
        LOGICAL ACTIVE
      END TYPE SlowReleaseNType

      INTEGER, PARAMETER :: NSlowRelN = 5
      Type (SlowReleaseNType), dimension(NSlowRelN) :: SlowRelN
      INTEGER NActiveSR

      TYPE NIType
        REAL NIEFF
        INTEGER NIEND, NILYR
      END TYPE NIType
      TYPE (NIType) NIData

      TYPE UIType
        REAL UIEFF
        INTEGER UIEND, UILYR
      END TYPE UIType
      TYPE (UIType) UIData

!=======================================================================
      CONTAINS

!=======================================================================
!  FertTypeRead
!  Reads fertilizer characteristics from FERTCH???.SDA file.
!-----------------------------------------------------------------------
!  REVISION       HISTORY
!  2019-06-04 CHP Written
!=======================================================================

      SUBROUTINE FertTypeRead(CONTROL)
!-------------------------------------------------------------------
      IMPLICIT NONE
      SAVE

      TYPE (ControlType), INTENT(IN) :: CONTROL

      LOGICAL FEXIST
      INTEGER LUNF, LNUM, LINE, FOUND, ISECT, I
      INTEGER ERR, Ftype
      INTEGER PFLAG

      CHARACTER*8, PARAMETER :: ERRKEY = 'FertRead'
      CHARACTER*150 FILEF, CHARTEST

      CHARACTER*12 NAMEF
      CHARACTER*80 PATHSD

!-----------------------------------------------------------------------
!     Read from FERTCH???.SDA file
      FILEF = 'FERCH' // ModelVerTxt // '.SDA'
      LNUM = 0
      INQUIRE (FILE = FILEF, EXIST = FEXIST)

      IF (.NOT. FEXIST) THEN
        CALL PATH('STD',CONTROL%DSSATP,PATHSD,PFLAG,NAMEF)
        FILEF = TRIM(PATHSD) // FILEF
      ENDIF

      INQUIRE (FILE = FILEF, EXIST = FEXIST)
      IF (.NOT. FEXIST) THEN
        CALL Fert_ERROR(29, ERRKEY, FILEF, LNUM, LUNF)
        RETURN
      ENDIF

      CALL GETLUN('FILEF', LUNF)
      OPEN (LUNF,FILE=FILEF,STATUS='OLD',IOSTAT=ERR)
      IF (ERR > 0) THEN
        CALL Fert_ERROR(10, ERRKEY, FILEF, LNUM, LUNF)
        RETURN
      ENDIF

!     Set check variable to zero
      FertFile % Check = 0

!     Read characteristics for each fertilizer type
      REWIND (LUNF)
      CALL FIND2(LUNF,'@CDE',LINE,FOUND); LNUM = LNUM + LINE
      IF (FOUND == 1) THEN
        I = 0
        DO WHILE (.TRUE.)
          CALL IGNORE(LUNF,LNUM,ISECT,CHARTEST)
          IF (ISECT == 1) THEN
            READ(CHARTEST,'(2X,I3)',IOSTAT=ERR) Ftype
            IF (ERR .GT. 0) CALL Fert_ERROR(64, ERRKEY, FILEF,LNUM,LUNF)
            IF (Ftype .GT. NFertTypes)
     &        CALL Fert_ERROR(42, ERRKEY, FILEF, LNUM, LUNF)
            READ(CHARTEST,'(6X,A35,3X,8F6.0,3X,A3,F6.0,6A6)',IOSTAT=ERR)
     &        FertFile(Ftype) % FertName,
     &        FertFile(Ftype) % NO3_N_pct,
     &        FertFile(Ftype) % NH4_N_pct,
     &        FertFile(Ftype) % UREA_N_pct,
     &        FertFile(Ftype) % UIEFF,
     &        FertFile(Ftype) % UIDUR,
     &        FertFile(Ftype) % NIEFF,
     &        FertFile(Ftype) % NIDUR,
     &        FertFile(Ftype) % NRL50,
     &        FertFile(Ftype) % NRFNC,
     &        FertFile(Ftype) % NSIGK,
     &        FertFile(Ftype) % FertN,
     &        FertFile(Ftype) % FertP,
     &        FertFile(Ftype) % FertK,
     &        FertFile(Ftype) % FertC,
     &        FertFile(Ftype) % FertMg,
     &        FertFile(Ftype) % FertS
            IF (ERR .GT. 0) CALL Fert_ERROR(64, ERRKEY, FILEF,LNUM,LUNF)
            I = I + 1
            FertFile(Ftype) % Check = 1
          ELSE
            Number_of_Fertilizers = I
            EXIT
          ENDIF
        ENDDO

      ELSE
        CALL Fert_ERROR(64, ERRKEY, FILEF, LNUM, LUNF)
        RETURN
      ENDIF

      RETURN
      END SUBROUTINE FertTypeRead
!======================================================================

!=======================================================================
!  FertProps
!  Sends fertilizer properties to FertPlace as needed.
!-----------------------------------------------------------------------
!  REVISION       HISTORY
!  2019-06-04 CHP Written
!=======================================================================

      SUBROUTINE FertProps()
!-------------------------------------------------------------------
      IMPLICIT NONE
      SAVE
!-----------------------------------------------------------------------

      RETURN
      END SUBROUTINE FertProps
!======================================================================

!=======================================================================
!  FERTILIZERTYPE, Subroutine for fertilizer placement module.
!  Determines fertilizer type
!-----------------------------------------------------------------------
!  Revision history
!  09/26/2003 AJG Separated this code from FPLACE_C into a new
!                 subroutine, restructured it and added P.
!  12/17/2003 AJG Renamed all the CH_ and CHEM_ variables to P_ variables.
!  05/06/2004 AJG Removed a few errors and the line in which FERTYPE 17
!                 was set back to 12.
!  01/20/2005 CHP Moved to generic fertilizer placement module, removed
!                 "_C" from name and removed fertilizer distribution.
!
!-----------------------------------------------------------------------
!  Called : FPLACE_C
!  Calls  : DISTRIB_SOILPi_C
!=======================================================================

      SUBROUTINE FERTILIZERTYPE (ISWITCH,
     &    ANFER, APFER, AKFER, FERTYPE, FERTYPE_CDE,      !Input
     &    HASN, HASP, HASK, HASUI, HASNI, HASCR)          !Output

!     ------------------------------------------------------------------
      IMPLICIT  NONE

      LOGICAL HASN, HASP, HASK  ! N, P, K
      LOGICAL HASUI             ! Urease inhibitor
      LOGICAL HASNI             ! Nitrification inhibitor
      LOGICAL HASCR             ! Controlled release
      INTEGER FERTYPE
      REAL ANFER, APFER, AKFER
      CHARACTER*5  FERTYPE_CDE
      CHARACTER*6, PARAMETER :: ERRKEY = 'FPLACE'
      CHARACTER*35 FERTYPE_TEXT
      CHARACTER*78 MSG(5)

      Type (SwitchType) ISWITCH

      HASN = .FALSE.
      HASP = .FALSE.
      HASK = .FALSE.
      HASUI= .FALSE.
      HASNI= .FALSE.
      HASCR= .FALSE.

!     ------------------------------------------------------------------
!     Does this fertilizer have N?
      CALL HASE_check(FertFile(Fertype) % FertN, HASN)
      IF (.NOT. HASN .AND. ANFER > 1.E-6 .AND. ISWITCH%ISWNIT == 'Y')
     &          THEN
        CALL READ_DETAIL(5, 35, FERTYPE_CDE, '*Ferti', FERTYPE_TEXT)
        MSG(1) = "Invalid fertilizer data in experiment file."
        WRITE(MSG(2),'(A,A,A,A)')
     &    "Fertilizer type: ",FERTYPE_CDE, " ",FERTYPE_TEXT
        WRITE(MSG(3),'(A,F8.1,A)')"Amount N specified: ",ANFER," kg/ha"
        MSG(4)="Check DETAIL.CDE file for valid fertilizer types."
        MSG(5)="No N applied on this date."
        CALL WARNING(5, ERRKEY, MSG)
      ENDIF

!     Does this fertilizer have P?
      CALL HASE_check(FertFile(Fertype) % FertP, HASP)
      IF (.NOT. HASP .AND. APFER > 1.E-6 .AND. ISWITCH % ISWPHO == 'Y')
     &          THEN
        CALL READ_DETAIL(5, 35, FERTYPE_CDE, '*Ferti', FERTYPE_TEXT)
        MSG(1) = "Invalid fertilizer data in experiment file."
        WRITE(MSG(2),'(A,A,A,A)')
     &    "Fertilizer type: ",FERTYPE_CDE, " ",FERTYPE_TEXT
        WRITE(MSG(3),'(A,F8.1,A)')"Amount P specified: ",APFER," kg/ha"
        MSG(4)="Check DETAIL.CDE file for valid fertilizer types."
        MSG(5)="No P applied on this date."
        CALL WARNING(5, ERRKEY, MSG)
      ENDIF

!     Does this fertilizer have K?
      CALL HASE_check(FertFile(Fertype) % FertK, HASK)
      IF (.NOT. HASK .AND. AKFER > 1.E-6 .AND. ISWITCH % ISWPOT == 'Y')
     &          THEN
        CALL READ_DETAIL(5, 35, FERTYPE_CDE, '*Ferti', FERTYPE_TEXT)
        MSG(1) = "Invalid fertilizer data in experiment file."
        WRITE(MSG(2),'(A,A,A,A)')
     &    "Fertilizer type: ",FERTYPE_CDE, " ",FERTYPE_TEXT
        WRITE(MSG(3),'(A,F8.1,A)')"Amount K specified: ",AKFER," kg/ha"
        MSG(4)="Check DETAIL.CDE file for valid fertilizer types."
        MSG(5)="No K applied on this date."
        CALL WARNING(5, ERRKEY, MSG)
      ENDIF

!     Does this fertilizer have urease inhibitor?
      IF (FertFile(Fertype) % UIEFF .GT. 1.E-6 .AND.
     &    FertFile(Fertype) % UIDUR .GT. 1.E-6) THEN
        HASUI = .TRUE.
      ELSE
        HASUI = .FALSE.
      ENDIF

!     Does this fertilizer have nitrification inhibitor?
      IF (FertFile(Fertype) % NIEFF .GT. 1.E-6 .AND.
     &    FertFile(Fertype) % NIDUR .GT. 1.E-6) THEN
        HASNI = .TRUE.
      ELSE
        HASNI = .FALSE.
      ENDIF

!     Does this fertilizer have controlled or slow release?
      IF (FertFile(Fertype) % NRL50 .GT. 1.E-6) THEN
        HASCR = .TRUE.

!       Set the sigma k value based on curve type
        IF (FertFile(Fertype) % NRFNC .EQ. "LIN") THEN
          FertFile(Fertype) % NSIGK = 0.2
        ELSEIF (FertFile(Fertype) % NRFNC .EQ. "STP") THEN
          FertFile(Fertype) % NSIGK = 2.0
        ELSE
          IF (FertFile(Fertype) % NSIGK .LT. 0.2) THEN
            FertFile(Fertype) % NSIGK = 0.2
          ENDIF
          IF (FertFile(Fertype) % NSIGK .GT. 2.0) THEN
            FertFile(Fertype) % NSIGK = 2.0
          ENDIF
        ENDIF
      ELSE
        HASCR = .FALSE.
      ENDIF

      RETURN
      END SUBROUTINE FERTILIZERTYPE

!=======================================================================
!     Check if element (N, P, or K) is supposed to be found in fertilizer
      Subroutine HASE_check(FertE, HasE)
      Implicit None

      CHARACTER*6 FertE
      LOGICAL HasE
      INTEGER ERR
      REAL Fert_content

      IF (FertE(4:6) .EQ. 'var' .OR. FertE(4:6) .EQ. 'VAR') THEN
        HASE = .TRUE.
        RETURN
      ENDIF

      READ (FertE,'(F6.0)', IOSTAT=ERR) Fert_content

      IF (ERR > 0) THEN
        HASE = .FALSE.
        RETURN
      ENDIF

      IF (Fert_content .GT. 1.E-6) THEN
        HASE = .TRUE.
        RETURN
      ENDIF

      End Subroutine HASE_check

!=======================================================================
      SUBROUTINE Fert_ERROR(IERROR, ERRKEY, FILEF, LNUM, LUNF)
!     Error handling routine for fertilizer characteristics input

      IMPLICIT NONE

      CHARACTER* 6 ERRKEY
      CHARACTER*78 MSG(4)
      CHARACTER*150 FILEF
      INTEGER IERROR, LNUM, LUNF

      CLOSE (LUNF)

      MSG(2) = FILEF(1:78)
      WRITE(MSG(3),'(A,I4)') "Line number ", LNUM

      SELECT CASE(IERROR)

      CASE (10)
        MSG(1) ='Error opening fertilizer characteristics file.'
        WRITE(MSG(4),'(A)') "Program will stop."
        CALL WARNING(4,ERRKEY,MSG)
        CALL ERROR(ERRKEY,IERROR,FILEF,LNUM)

      CASE (29)
        MSG(1) = 'Fertilizer characteristics file does not exist.'
        WRITE(MSG(4),'(A)') "Program will stop."
        CALL WARNING(4,ERRKEY,MSG)
        CALL ERROR(ERRKEY,IERROR,FILEF,LNUM)

      CASE (42)
        MSG(1) ='Invalid fertilizer type in file.'
        WRITE(MSG(4),'(A)') "Program will stop."
        CALL WARNING(4,ERRKEY,MSG)
        CALL ERROR(ERRKEY,IERROR,FILEF,LNUM)

      CASE (64)
        MSG(1) ='Invalid format in fertilizer characteristics file.'
        WRITE(MSG(4),'(A)') "Program will stop."
        CALL WARNING(4,ERRKEY,MSG)
        CALL ERROR(ERRKEY,IERROR,FILEF,LNUM)

      END SELECT

      RETURN
      END SUBROUTINE Fert_ERROR
!======================================================================
      END MODULE FertType_mod
!======================================================================

ext: for

max_stars_repo_name: ACJingQi/dssat-csm-os

max_issues_repo_name: ACJingQi/dssat-csm-os

max_forks_repo_name: ACJingQi/dssat-csm-os


text:
  CALL MatrixColumnNorm(this,column)
  norm = MAXVAL(column)

max_stars_repo_name: Kokookster/NTPoly

max_issues_repo_name: omni-compiler/NTPoly

max_forks_repo_name: omni-compiler/NTPoly