hump.f

Go to the documentation of this file.

      subroutine hump(const,t,tp,nm,lp,rmu,thd,ixy,nmu)                       
c                                                                       
c***********************************************************************
      implicit real*8(a-h,o-z)                                          
      parameter(nw=250,nth=46,nsz=45,nph=91,nsg=10,nstk=4,ntrx=16,npti=8,
     1       nlyr=200,ntf=1801,nwnd=8,nrhum=10,nmd=3)
      real*8 yy,xx,const,lp
      real*8 t(ntf,nstk),am(nlyr),costh(nlyr),sinth(nlyr),p(nstk,nstk)
      real*8 tp(nlyr,nstk),cthr(nlyr),thr(nlyr),dmu(nlyr),rmu(51),thd(ntf)
c***********************************************************************
c    program to compute pbar for each mup when mu=mup and delfi=0       
c                                                                       
      izia=0
c     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 *t3 )
                  p(2,1)=p(2,1)+sfisq*(amups*t2+amusq*t1+2.0d0 * amumu *t3)
                  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))*0.5d0
      enddo
      if(izia.eq.1)stop
c***********************************************************************
c***********************************************************************
      return                                                            
      end                                                               
c***********************************************************************