hump.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	hump (const, t, tp, nm, lp, rmu, thd, ixy, nmu)

Function/Subroutine Documentation

hump()

subroutine hump	(	real*8	const,
		real*8, dimension(ntf,nstk)	t,
		real*8, dimension(nlyr,nstk)	tp,
			nm,
		real*8	lp,
		real*8, dimension(51)	rmu,
		real*8, dimension(nt !*********************************************************************** ! program to compute pbar for each mup when mu=mup and delfi=0 ! izia=0 ! write(6,*)'entering hump' pi=dacos(-1.0d0 ) conv =pi/180.0d0 p(:,:)=0.0d0 dp=0.10d0 dt=0.10d0 fi1=lp/2 ni1=fi1/dp+1 dp=fi1/dfloat(ni1) do k=1,ni1 fk=dfloat(k-1)*dp if(ixy.ne.0) then phi=(180-fi1+fk+dp/2.0d0)*conv else phi=(fk+dp/2.0d0)*conv endif sinth(k) = dsin(phi) costh(k) = dcos(phi) enddo deldeg=dt do kl=1,nm/2 deg=rmu(kl+1)-rmu(kl) nodeg=deg/deldeg deldeg=deg/dfloat(nodeg) nodp1=nodeg+1 xx=deg/2.0d0 yy=xx-deldeg/2.0d0 thetd=(rmu(kl)+rmu(kl+1))/2.0d0 amup = dcos(thetd*conv) do i=1,nodp1 fi=dfloat(i-1)*deldeg cthr(i)=(thetd-xx+fi)*conv enddo do i=1,nodeg fi=dfloat(i-1)*deldeg thr(i)=(thetd-yy+fi)*conv am(i) = dcos(thr(i)) dmu(i)=dcos(cthr(i))-dcos(cthr(i+1)) enddo ddmu =dcos(rmu(kl)*conv)-dcos(rmu(kl+1)*conv) amups=amup**2 do i=1,4 k=(kl-1)*4+i do j=1,4 tp(k,j)=0.0d0 enddo enddo do i=1,nodeg if(ixy.ne.0) then amu=-am(i) else amu=am(i) endif p(:,:)=0.0d0 amusq=amu**2 amumu=amu*amup anunu = dsqrt((1.0d0-amups)*(1.0d0-amusq)) do j=1,ni1 siin=sinth(j) coosn=costh(j) copsi=anunu+amumu*coosn sifi=copsi*coosn sfisq=siin**2 cfisq=coosn**2 sisq=copsi**2 a=amup*coosn b=amup*copsi c=amu*coosn d=amu*copsi e=sfisq*(amusq+amups) g=amumu*sfisq h=sfisq*(amusq-amups) x=amumu+anunu*coosn yz= dacos(x)/conv yzz=10.0d0*yz m=idint(yzz)+1 if(m.eq.0)m=1 ya=yz-thd(m) dthd=thd(m)-thd(m+1) call xntpln(yz,thd(m),thd(m+1),t(m,1),t(m+1,1),t1) call xntpln(yz,thd(m),thd(m+1),t(m,2),t(m+1,2),t2) call xntpln(yz,thd(m),thd(m+1),t(m,3),t(m+1,3),t3) call xntpln(yz,thd(m),thd(m+1),t(m,4),t(m+1,4),t4) p(3,3)=p(3,3)+(sifi-g)*(t1+t2)+(cfisq+sisq-e)*t3 p(2,2)=p(2,2)+sisq*t1+cfisq*t2+2.0d0 *sifi * t3 p(1,2)=p(1,2)+sfisq*(amups*t1+amusq*t2+2.0d0 * amumu * p(2,1)=p(2,1)+sfisq*(amups*t2+amusq*t1+2.0d0 * amumu * p(1,1)=p(1,1)+cfisq*t1+sisq*t2+2.0d0 * sifi *t3 p(4,3)=p(4,3)+(sisq-cfisq-h)*t4 p(3,4)=p(3,4)+(cfisq-sisq-h)*t4 p(4,4)=p(4,4)+(sifi+g)*(t1+t2)+(cfisq+sisq+e)*t3 enddo do k=1,4 kk=(kl-1)*4+k do l=1,4 ttpp1=const*p(k,1)*dmu(i)*dp/fi1 ttpp2=const*p(k,2)*dmu(i)*dp/fi1 ttpp=0.5d0*(ttpp1+ttpp2) tp(kk,l)=p(k,l)*dmu(i)*dp/fi1+tp(kk,l) km=kk-(kl-1)*4 enddo enddo enddo do k=1,4 kk=(kl-1)*4+k do l=1,4 tp(kk,l)=tp(kk,l)/ddmu enddo enddo id=(kl-1)*4 ppbar=const*(tp(id+1,1)+tp(id+1,2)+tp(id+2,1)+tp(id+2,2))	thd,
			ixy,
			nmu
	)

Definition at line 2 of file hump.f.