text:
      subroutine tqente(trans,cod,beam,al,bz,irad)
      use mathfun, only: sqrtl
      use element_drift_common, only:tsoldz
      use temw,only:tmulbs
      implicit none
      integer*4 ,intent(in):: irad
      real*8 ,intent(in):: al,bz
      real*8 ,intent(inout):: trans(6,12),cod(6),beam(42)
      real*8 trans1(6,6),trans2(6,6),pr,a,pz,r,pxi,pyi,bzh,
     $     z0,f,g,fx,fy,dpz,rpz
      real*8, parameter ::ptmin=0.999999d0
      pr=1.d0+cod(6)
      if(bz .eq. 0.d0)then
        pxi=cod(2)
        pyi=cod(4)
        a=pxi**2+pyi**2
        dpz=-a/pr/(1.d0+sqrtl(1.d0-a/pr**2))
        pz=pr+dpz
        r=-dpz/pz/pr*al
        rpz=pz/pr
        cod(1)=cod(1)+r*pxi
        cod(3)=cod(3)+r*pyi
        cod(5)=cod(5)-(2.d0+rpz)*a/2.d0/(pr+pz)*r
        f=al/pz**3
        fx=f*pxi
        fy=f*pyi
        g=dpz*(1.d0+rpz+rpz**2)
        if(irad .gt. 6)then
          call tinitr(trans2)
        endif
        trans2(1,2)=r+pxi*fx
        trans2(1,4)=pyi*fx
        trans2(1,6)=g*fx
        trans2(3,2)=pxi*fy
        trans2(3,4)=r+pyi*fy
        trans2(3,6)=g*fy
        trans2(5,2)=trans2(1,6)
        trans2(5,4)=trans2(3,6)
        trans2(5,6)=-a*g/pr*f
        trans(1:5:2,:)=trans(1:5:2,:)+
     $       matmul(trans2(1:5:2,2:6:2),trans(2:6:2,:))
      else
        bzh=bz*.5d0
c cod are canonical!
        pxi=cod(2)+bzh*cod(3)
        pyi=cod(4)-bzh*cod(1)
        a=pxi**2+pyi**2
        pz=pr*sqrtl(1.d0-a/pr**2)
        r=a/pr/(pr+pz)*al
        call tinitr(trans2)
        trans2(2,3)= bzh
        trans2(4,1)=-bzh
        trans2(5,2)= pxi/pr/pz*al
        trans2(5,3)= pxi*bzh/pr/pz*al
        trans2(5,4)= pyi/pr/pz*al
        trans2(5,1)=-pyi*bzh/pr/pz*al
        trans2(5,5)=0.d0
        trans2(5,6)=-a/pr**2/pz*al
        cod(2)=pxi
        cod(4)=pyi
        z0=cod(5)
        call tsoldz(trans1,cod,r,0.d0,0.d0,bz,.false.)
        trans2=matmul(trans1,trans2)
c        call tmultr(trans2,trans1,6)
        cod(5)=z0-(a/pr/(pr+pz))**2*(2.d0*pr+pz)/pz/2.d0*al
        trans2(5,5)=1.d0
        trans2(5,6)=a**2*(pr**2+pr*pz+pz**2)/(pr*pz)**3/(pr+pz)*al
        cod(2)=cod(2)-bzh*cod(3)
        cod(4)=cod(4)+bzh*cod(1)
        trans2(2,1)=trans2(2,1)-bzh*trans2(3,1)
        trans2(4,1)=trans2(4,1)+bzh*trans2(1,1)
        trans2(2,2)=trans2(2,2)-bzh*trans2(3,2)
        trans2(4,2)=trans2(4,2)+bzh*trans2(1,2)
        trans2(2,3)=trans2(2,3)-bzh*trans2(3,3)
        trans2(4,3)=trans2(4,3)+bzh*trans2(1,3)
        trans2(2,4)=trans2(2,4)-bzh*trans2(3,4)
        trans2(4,4)=trans2(4,4)+bzh*trans2(1,4)
        trans2(2,6)=trans2(2,6)-bzh*trans2(3,6)
        trans2(4,6)=trans2(4,6)+bzh*trans2(1,6)
        trans2(5,1)=
     $       -trans2(1,1)*trans2(2,6)
     $       +trans2(2,1)*trans2(1,6)
     $       -trans2(3,1)*trans2(4,6)
     $       +trans2(4,1)*trans2(3,6)
        trans2(5,2)=
     $       -trans2(1,2)*trans2(2,6)
     $       +trans2(2,2)*trans2(1,6)
     $       -trans2(3,2)*trans2(4,6)
     $       +trans2(4,2)*trans2(3,6)
        trans2(5,3)=
     $       -trans2(1,3)*trans2(2,6)
     $       +trans2(2,3)*trans2(1,6)
     $       -trans2(3,3)*trans2(4,6)
     $       +trans2(4,3)*trans2(3,6)
        trans2(5,4)=
     $       -trans2(1,4)*trans2(2,6)
     $       +trans2(2,4)*trans2(1,6)
     $       -trans2(3,4)*trans2(4,6)
     $       +trans2(4,4)*trans2(3,6)
        call tmultr5(trans,trans2,irad)
c      write(*,'(a/,6(1p6g11.4/))')
c     $     'tqente-2 ',((trans(i,j),j=1,6),i=1,6)
c      write(*,*)'with irad: ',irad
      endif
      if(irad .gt. 6)then
        call tmulbs(beam,trans2,.true.)
      endif
      return
      end

max_stars_repo_name: KatsOide/SAD

max_issues_repo_name: KatsOide/SAD

max_forks_repo_name: KatsOide/SAD


text:
      PROGRAM xdfridr
C     driver for routine dfridr
      REAL dfridr,dx,err,func,h,x
      EXTERNAL func
1     write(*,*) 'input x,h '
      read(*,*,END=999) x,h
      dx=dfridr(func,x,h,err)
      write(*,*) 'DFRIDR=',dx,1./cos(x)**2,err
      goto 1
999   END
      REAL FUNCTION func(x)
      REAL x
      func=tan(x)
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xdfridr.for


text:
      subroutine qqb_z2jetx(p,msq,mqq,msqx,msqx_cs)
      implicit none
c---Matrix element squared averaged over initial colors and spins
c     q(-p1)+qbar(-p2) --> Z +g(p5) +g(p6)
c                          |
c                          --> l(p3)+a(p4)
c
c--all momenta incoming
      include 'constants.f'
      include 'masses.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'zcouple.f'
      include 'ewcharge.f'
      include 'sprods_com.f'
      include 'zprods_com.f'
      include 'flags.f'
      include 'lc.f'
      include 'part.f'
      integer i,j,k,f,pq,pl,nquark,nup,ndo,
     .   j1,j2,j3,icol
      double precision msq(-nf:nf,-nf:nf),p(mxpart,4),fac,faclo,
     .   qqbZgg2_cs(0:2,2,2),qbqZgg2_cs(0:2,2,2),
     .   qgZqg2_cs(0:2,2,2),gqZqg2_cs(0:2,2,2),
     .   qgZgq2_cs(0:2,2,2),gqZgq2_cs(0:2,2,2),
     .   qbgZqbg2_cs(0:2,2,2),gqbZqbg2_cs(0:2,2,2),
     .   qbgZgqb2_cs(0:2,2,2),gqbZgqb2_cs(0:2,2,2),
     .   ggZqbq2_cs(0:2,2,2),ggZqqb2_cs(0:2,2,2),ggtemp(0:2),
     .   qqbZgg2(2,2),
     .   qgZqg2(2,2),
c     .   gqbZgqb2(2,2),
     .   gqZgq2(2,2),
     .   qgZgq2(2,2),
c     .   gqbZqbg2(2,2),
     .   gqZqg2(2,2),
     .   ggZqbq2(2,2)
c     .   qbgZgqb2(2,2),
c     .   ggZqqb2(2,2),
c     .   qbgZqbg2(2,2),
c     .   qbqZgg2(2,2)
      double precision tup,tdo

      double complex qRb_a(2,2,2),qRb_b(2,2,2)
      double complex qqb_a(2,2,2),qqb_b(2,2,2),prop

      double complex qbq_a(2,2,2),qbq_b(2,2,2)
      double complex qbR_a(2,2,2),qbR_b(2,2,2)

      double complex qq_a(2,2,2),qq_b(2,2,2)
      double complex qR_a(2,2,2),qR_b(2,2,2)

      double complex qbRb_a(2,2,2),qbRb_b(2,2,2)
      double complex qbqb_a(2,2,2),qbqb_b(2,2,2)

      double complex qRb_ax(2,2,2),qRb_bx(2,2,2)
      double complex qqb_ax(2,2,2),qqb_bx(2,2,2)

      double complex qbR_ax(2,2,2),qbR_bx(2,2,2)
      double complex qbq_ax(2,2,2),qbq_bx(2,2,2)

      double precision mqq(0:2,fn:nf,fn:nf),msq0,msq1,msq2
      double precision msqx(0:2,-nf:nf,-nf:nf,-nf:nf,-nf:nf)
      double precision msqx_cs(0:2,-nf:nf,-nf:nf)

      integer rcolourchoice

      logical rGflag

      integer,parameter::swap(2)=(/2,1/),swap1(0:2)=(/0,2,1/)

c--- if we're calculating the REAL or VIRT matrix elements, we
c--- need all the colour structures, but want to preserve
c--- the actual value of colourchoice
      if ((part .eq. 'real') .or. (part .eq. 'virt')) then
        rcolourchoice=colourchoice
        colourchoice=0
      endif
c--- if we're calculating the REAL matrix elements with Qflag=TRUE,
c    the subtraction terms involve the (Gflag=TRUE) matrix elements
      if ((part .eq. 'real') .and. (Qflag .eqv. .true.)) then
        rGflag=Gflag
        Gflag=.true.
      endif

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo

      call spinoru(6,p,za,zb)

      prop=s(3,4)/dcmplx((s(3,4)-zmass**2),zmass*zwidth)

c--- calculate 2-quark, 2-gluon amplitudes
      if (Gflag) then
        call z2jetsq(1,2,3,4,5,6,za,zb,qqbZgg2)
        call storecsz(qqbZgg2_cs)
        call z2jetsq(1,5,3,4,2,6,za,zb,qgZqg2)
        call storecsz(qgZqg2_cs)
        call z2jetsq(2,5,3,4,1,6,za,zb,gqZqg2)
        call storecsz(gqZqg2_cs)
        call z2jetsq(1,6,3,4,2,5,za,zb,qgZgq2)
        call storecsz(qgZgq2_cs)
        call z2jetsq(2,6,3,4,1,5,za,zb,gqZgq2)
        call storecsz(gqZgq2_cs)
C --NB this is the matrix element for gg->Z qb(5) q(6)
        call z2jetsq(5,6,3,4,1,2,za,zb,ggZqbq2)
        call storecsz(ggZqbq2_cs)

        do j=1,2
        do k=1,2
c        qbqZgg2(j,k)=qqbZgg2(swap(j),k)
c        qbgZqbg2(j,k)=qgZqg2(swap(j),k)
c        gqbZqbg2(j,k)=gqZqg2(swap(j),k)
c        qbgZgqb2(j,k)=qgZgq2(swap(j),k)
c        gqbZgqb2(j,k)=gqZgq2(swap(j),k)
c        ggZqqb2(j,k)=ggZqbq2(swap(j),k)
        do i=0,2
        qbqZgg2_cs(i,j,k)=qqbZgg2_cs(swap1(i),swap(j),k)
        qbgZqbg2_cs(i,j,k)=qgZqg2_cs(swap1(i),swap(j),k)
        gqbZqbg2_cs(i,j,k)=gqZqg2_cs(swap1(i),swap(j),k)
        qbgZgqb2_cs(i,j,k)=qgZgq2_cs(swap1(i),swap(j),k)
        gqbZgqb2_cs(i,j,k)=gqZgq2_cs(swap1(i),swap(j),k)
        ggZqqb2_cs(i,j,k)=ggZqbq2_cs(swap1(i),swap(j),k)
        enddo
        enddo
        enddo

c        call z2jetsq(2,1,3,4,5,6,za,zb,qbqZgg2)
c        call storecsz(qbqZgg2_cs)
c        call z2jetsq(5,1,3,4,2,6,za,zb,qbgZqbg2)
c        call storecsz(qbgZqbg2_cs)
c        call z2jetsq(5,2,3,4,1,6,za,zb,gqbZqbg2)
c        call storecsz(gqbZqbg2_cs)
c        call z2jetsq(6,1,3,4,2,5,za,zb,qbgZgqb2)
c        call storecsz(qbgZgqb2_cs)
c        call z2jetsq(6,2,3,4,1,5,za,zb,gqbZgqb2)
c        call storecsz(gqbZgqb2_cs)
C --NB this is the matrix element for gg->Z q(5) qb(6)
c        call z2jetsq(6,5,3,4,1,2,za,zb,ggZqqb2)
c        call storecsz(ggZqqb2_cs)

        fac=v*xn/four*(esq*gsq)**2
        do pq=1,2
        do pl=1,2
        do i=0,2
          qqbZgg2_cs(i,pq,pl) = half*aveqq*fac*qqbZgg2_cs(i,pq,pl)
          qbqZgg2_cs(i,pq,pl) = half*aveqq*fac*qbqZgg2_cs(i,pq,pl)
          gqZqg2_cs(i,pq,pl)  = aveqg*fac*gqZqg2_cs(i,pq,pl)
          qgZqg2_cs(i,pq,pl)  = aveqg*fac*qgZqg2_cs(i,pq,pl)
          gqbZqbg2_cs(i,pq,pl)= aveqg*fac*gqbZqbg2_cs(i,pq,pl)
          qbgZqbg2_cs(i,pq,pl)= aveqg*fac*qbgZqbg2_cs(i,pq,pl)

          gqZgq2_cs(i,pq,pl)  = aveqg*fac*gqZgq2_cs(i,pq,pl)
          qgZgq2_cs(i,pq,pl)  = aveqg*fac*qgZgq2_cs(i,pq,pl)
          gqbZgqb2_cs(i,pq,pl)= aveqg*fac*gqbZgqb2_cs(i,pq,pl)
          qbgZgqb2_cs(i,pq,pl)= aveqg*fac*qbgZgqb2_cs(i,pq,pl)

          ggZqbq2_cs(i,pq,pl) = avegg*fac*ggZqbq2_cs(i,pq,pl)
          ggZqqb2_cs(i,pq,pl) = avegg*fac*ggZqqb2_cs(i,pq,pl)
       enddo

        qqbZgg2(pq,pl) = qqbZgg2_cs(1,pq,pl)+qqbZgg2_cs(2,pq,pl)
     .                  +qqbZgg2_cs(0,pq,pl)
        gqZqg2(pq,pl)  = gqZqg2_cs(1,pq,pl) +gqZqg2_cs(2,pq,pl)
     .                  +gqZqg2_cs(0,pq,pl)
        qgZqg2(pq,pl)  = qgZqg2_cs(1,pq,pl)  +qgZqg2_cs(2,pq,pl)
     .                  +qgZqg2_cs(0,pq,pl)
        gqZgq2(pq,pl)  = gqZgq2_cs(1,pq,pl) +gqZgq2_cs(2,pq,pl)
     .                  +gqZgq2_cs(0,pq,pl)
        qgZgq2(pq,pl)  = qgZgq2_cs(1,pq,pl)  +qgZgq2_cs(2,pq,pl)
     .                  +qgZgq2_cs(0,pq,pl)
        ggZqbq2(pq,pl) = ggZqbq2_cs(1,pq,pl) +ggZqbq2_cs(2,pq,pl)
     .                  +ggZqbq2_cs(0,pq,pl)
c        gqbZqbg2(pq,pl)= gqbZqbg2_cs(1,pq,pl)+gqbZqbg2_cs(2,pq,pl)
c     .                  +gqbZqbg2_cs(0,pq,pl)
c        qbqZgg2(pq,pl) = qbqZgg2_cs(1,pq,pl)+qbqZgg2_cs(2,pq,pl)
c     .                  +qbqZgg2_cs(0,pq,pl)
c        qbgZqbg2(pq,pl)= qbgZqbg2_cs(1,pq,pl)+qbgZqbg2_cs(2,pq,pl)
c     .                  +qbgZqbg2_cs(0,pq,pl)
c        gqbZgqb2(pq,pl)= gqbZgqb2_cs(1,pq,pl)+gqbZgqb2_cs(2,pq,pl)
c     .                  +gqbZgqb2_cs(0,pq,pl)
c        qbgZgqb2(pq,pl)= qbgZgqb2_cs(1,pq,pl)+qbgZgqb2_cs(2,pq,pl)
c     .                  +qbgZgqb2_cs(0,pq,pl)
c        ggZqqb2(pq,pl) = ggZqqb2_cs(1,pq,pl) +ggZqqb2_cs(2,pq,pl)
c     .                  +ggZqqb2_cs(0,pq,pl)
        enddo
        enddo
      endif

      if (Qflag) then
      call spinoru(6,p,za,zb)

c--- qRb->qRb
      call ampqqb_qqb(1,5,2,6,qRb_a,qRb_b)
c--- qR->qR
c instead of calling ampqqb_qqb(1,5,6,2,qR_a,qR_b)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qR_a(j1,j2,j3)=+qRb_a(j1,swap(j2),j3)
      qR_b(j1,j2,j3)=-qRb_b(j1,swap(j2),j3)
      enddo
      enddo
      enddo
c--- qbR->qbR
      call ampqqb_qqb(6,1,5,2,qbR_a,qbR_b)

c--- qbRb->qbRb
c instead of calling ampqqb_qqb(5,1,2,6,qbRb_a,qbRb_b)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qbRb_a(j1,j2,j3)=-qRb_a(swap(j1),j2,j3)
      qbRb_b(j1,j2,j3)=+qRb_b(swap(j1),j2,j3)
      enddo
      enddo
      enddo

c--- qqb->qqb
      call ampqqb_qqb(1,2,5,6,qqb_a,qqb_b)
c--- qbq->qqb
c instead of calling ampqqb_qqb(2,1,5,6,qbq_a,qbq_b)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qbq_a(j1,j2,j3)=-qqb_a(swap(j1),j2,j3)
      qbq_b(j1,j2,j3)=+qqb_b(swap(j1),j2,j3)
      enddo
      enddo
      enddo

c--- qq->qq
      call ampqqb_qqb(1,6,5,2,qq_a,qq_b)
c--- qbqb->qbqb
c instead of calling ampqqb_qqb(6,1,2,5,qbqb_a,qbqb_b)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qbqb_a(j1,j2,j3)=-qq_a(swap(j1),swap(j2),j3)
      qbqb_b(j1,j2,j3)=-qq_b(swap(j1),swap(j2),j3)
      enddo
      enddo
      enddo

c--- four extra amplitudes added for 4-label matrix elements

c--- qRb->Rbq
c instead of calling ampqqb_qqb(1,6,2,5,qRb_ax,qRb_bx)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qRb_ax(j1,j2,j3)=-qbqb_a(swap(j1),j2,j3)
      qRb_bx(j1,j2,j3)=+qbqb_b(swap(j1),j2,j3)
      enddo
      enddo
      enddo

c--- qqb->qbq
c instead of calling ampqqb_qqb(1,2,6,5,qqb_ax,qqb_bx)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qqb_ax(j1,j2,j3)=+qqb_a(j1,swap(j2),j3)
      qqb_bx(j1,j2,j3)=-qqb_b(j1,swap(j2),j3)
      enddo
      enddo
      enddo

c--- qbR->qbR
c instead of calling ampqqb_qqb(5,1,6,2,qbR_ax,qbR_bx)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qbR_ax(j1,j2,j3)=-qRb_a(swap(j1),swap(j2),j3)
      qbR_bx(j1,j2,j3)=-qRb_b(swap(j1),swap(j2),j3)
      enddo
      enddo
      enddo

c--- qbq->qbq
c instead of calling ampqqb_qqb(2,1,6,5,qbq_ax,qbq_bx)
      do j1=1,2
      do j2=1,2
      do j3=1,2
      qbq_ax(j1,j2,j3)=+qbq_a(j1,swap(j2),j3)
      qbq_bx(j1,j2,j3)=-qbq_b(j1,swap(j2),j3)
      enddo
      enddo
      enddo

      faclo=4d0*V*gsq**2*esq**2*aveqq
      endif

      if (Gflag) then
      do j=-nf,nf
      do k=-nf,nf

      do icol=0,2
        msqx_cs(icol,j,k)=0d0
      enddo

      if( j .ne. 0 .and. k .ne. 0 .and. j .ne. -k) goto 19

      if     ((j .eq. 0) .and. (k .eq. 0)) then

           do icol=0,2
           ggtemp(icol)=0d0
           do nquark=1,nf
           msqx(icol,j,k,-nquark,nquark)=
     .      +abs(Q(nquark)*q1+L(nquark)*l1*prop)**2*ggZqbq2_cs(icol,1,1)
     .      +abs(Q(nquark)*q1+R(nquark)*r1*prop)**2*ggZqbq2_cs(icol,2,2)
     .      +abs(Q(nquark)*q1+L(nquark)*r1*prop)**2*ggZqbq2_cs(icol,1,2)
     .      +abs(Q(nquark)*q1+R(nquark)*l1*prop)**2*ggZqbq2_cs(icol,2,1)
           msqx(icol,j,k,nquark,-nquark)=
     .      +abs(Q(nquark)*q1+L(nquark)*l1*prop)**2*ggZqqb2_cs(icol,1,1)
     .      +abs(Q(nquark)*q1+R(nquark)*r1*prop)**2*ggZqqb2_cs(icol,2,2)
     .      +abs(Q(nquark)*q1+L(nquark)*r1*prop)**2*ggZqqb2_cs(icol,1,2)
     .      +abs(Q(nquark)*q1+R(nquark)*l1*prop)**2*ggZqqb2_cs(icol,2,1)
           ggtemp(icol)=ggtemp(icol)+msqx(icol,j,k,-nquark,nquark)
           enddo
           msqx_cs(icol,j,k)=ggtemp(icol)

           enddo
      elseif ((j .gt. 0) .and. (k .lt. 0)) then
          do icol=0,2
             msqx_cs(icol,j,k)=
     .       +abs(Q(j)*q1+L(j)*l1*prop)**2*qqbZgg2_cs(icol,1,1)
     .       +abs(Q(j)*q1+R(j)*r1*prop)**2*qqbZgg2_cs(icol,2,2)
     .       +abs(Q(j)*q1+L(j)*r1*prop)**2*qqbZgg2_cs(icol,1,2)
     .       +abs(Q(j)*q1+R(j)*l1*prop)**2*qqbZgg2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo

c---Statistical factor already included above
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
          do icol=0,2
             msqx_cs(icol,j,k)=
     .       +abs(Q(k)*q1+L(k)*l1*prop)**2*qbqZgg2_cs(icol,1,1)
     .       +abs(Q(k)*q1+R(k)*r1*prop)**2*qbqZgg2_cs(icol,2,2)
     .       +abs(Q(k)*q1+L(k)*r1*prop)**2*qbqZgg2_cs(icol,1,2)
     .       +abs(Q(k)*q1+R(k)*l1*prop)**2*qbqZgg2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo
      elseif ((j .gt. 0) .and. (k .eq. 0)) then
          do icol=0,2
c--- Madloop check: full Z in Born (no photon)
c             msqx_cs(icol,j,k)=
c     .       +abs(Q(j)*q1*0+L(j)*l1*prop)**2*qgZqg2_cs(icol,1,1)
c     .       +abs(Q(j)*q1*0+R(j)*r1*prop)**2*qgZqg2_cs(icol,2,2)
c     .       +abs(Q(j)*q1*0+L(j)*r1*prop)**2*qgZqg2_cs(icol,1,2)
c     .       +abs(Q(j)*q1*0+R(j)*l1*prop)**2*qgZqg2_cs(icol,2,1)
c--- Madloop check: axial coupling of Z only
c             msqx_cs(icol,j,k)=
c     .       +abs(Q(j)*q1*0+0.25d0/sin2w**2*prop)**2*qgZqg2_cs(icol,1,1)
c     .       +abs(Q(j)*q1*0+0.25d0/sin2w**2*prop)**2*qgZqg2_cs(icol,2,2)
c     .       +abs(Q(j)*q1*0+0.25d0/sin2w**2*prop)**2*qgZqg2_cs(icol,1,2)
c     .       +abs(Q(j)*q1*0+0.25d0/sin2w**2*prop)**2*qgZqg2_cs(icol,2,1)

c--- normal case
             msqx_cs(icol,j,k)=
     .       +abs(Q(j)*q1+L(j)*l1*prop)**2*qgZqg2_cs(icol,1,1)
     .       +abs(Q(j)*q1+R(j)*r1*prop)**2*qgZqg2_cs(icol,2,2)
     .       +abs(Q(j)*q1+L(j)*r1*prop)**2*qgZqg2_cs(icol,1,2)
     .       +abs(Q(j)*q1+R(j)*l1*prop)**2*qgZqg2_cs(icol,2,1)
             msqx(icol,j,k,j,k)=msqx_cs(icol,j,k)
             msqx(icol,j,k,k,j)=
     .       +abs(Q(j)*q1+L(j)*l1*prop)**2*qgZgq2_cs(icol,1,1)
     .       +abs(Q(j)*q1+R(j)*r1*prop)**2*qgZgq2_cs(icol,2,2)
     .       +abs(Q(j)*q1+L(j)*r1*prop)**2*qgZgq2_cs(icol,1,2)
     .       +abs(Q(j)*q1+R(j)*l1*prop)**2*qgZgq2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo
      elseif ((j .lt. 0) .and. (k .eq. 0)) then
          do icol=0,2
             msqx_cs(icol,j,k)=
     .       +abs(Q(-j)*q1+L(-j)*l1*prop)**2*qbgZqbg2_cs(icol,1,1)
     .       +abs(Q(-j)*q1+R(-j)*r1*prop)**2*qbgZqbg2_cs(icol,2,2)
     .       +abs(Q(-j)*q1+L(-j)*r1*prop)**2*qbgZqbg2_cs(icol,1,2)
     .       +abs(Q(-j)*q1+R(-j)*l1*prop)**2*qbgZqbg2_cs(icol,2,1)
             msqx(icol,j,k,j,k)=msqx_cs(icol,j,k)
             msqx(icol,j,k,k,j)=
     .       +abs(Q(-j)*q1+L(-j)*l1*prop)**2*qbgZgqb2_cs(icol,1,1)
     .       +abs(Q(-j)*q1+R(-j)*r1*prop)**2*qbgZgqb2_cs(icol,2,2)
     .       +abs(Q(-j)*q1+L(-j)*r1*prop)**2*qbgZgqb2_cs(icol,1,2)
     .       +abs(Q(-j)*q1+R(-j)*l1*prop)**2*qbgZgqb2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo
      elseif ((j .eq. 0) .and. (k .gt. 0)) then
          do icol=0,2
             msqx_cs(icol,j,k)=
     .       +abs(Q(k)*q1+L(k)*l1*prop)**2*gqZqg2_cs(icol,1,1)
     .       +abs(Q(k)*q1+R(k)*r1*prop)**2*gqZqg2_cs(icol,2,2)
     .       +abs(Q(k)*q1+L(k)*r1*prop)**2*gqZqg2_cs(icol,1,2)
     .       +abs(Q(k)*q1+R(k)*l1*prop)**2*gqZqg2_cs(icol,2,1)
             msqx(icol,j,k,k,j)=msqx_cs(icol,j,k)
             msqx(icol,j,k,j,k)=
     .       +abs(Q(k)*q1+L(k)*l1*prop)**2*gqZgq2_cs(icol,1,1)
     .       +abs(Q(k)*q1+R(k)*r1*prop)**2*gqZgq2_cs(icol,2,2)
     .       +abs(Q(k)*q1+L(k)*r1*prop)**2*gqZgq2_cs(icol,1,2)
     .       +abs(Q(k)*q1+R(k)*l1*prop)**2*gqZgq2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo
      elseif ((j .eq. 0) .and. (k .lt. 0)) then
          do icol=0,2
             msqx_cs(icol,j,k)=
     .       +abs(Q(-k)*q1+L(-k)*l1*prop)**2*gqbZqbg2_cs(icol,1,1)
     .       +abs(Q(-k)*q1+R(-k)*r1*prop)**2*gqbZqbg2_cs(icol,2,2)
     .       +abs(Q(-k)*q1+L(-k)*r1*prop)**2*gqbZqbg2_cs(icol,1,2)
     .       +abs(Q(-k)*q1+R(-k)*l1*prop)**2*gqbZqbg2_cs(icol,2,1)
             msqx(icol,j,k,k,j)=msqx_cs(icol,j,k)
             msqx(icol,j,k,j,k)=
     .       +abs(Q(-k)*q1+L(-k)*l1*prop)**2*gqbZgqb2_cs(icol,1,1)
     .       +abs(Q(-k)*q1+R(-k)*r1*prop)**2*gqbZgqb2_cs(icol,2,2)
     .       +abs(Q(-k)*q1+L(-k)*r1*prop)**2*gqbZgqb2_cs(icol,1,2)
     .       +abs(Q(-k)*q1+R(-k)*l1*prop)**2*gqbZgqb2_cs(icol,2,1)
             msqx_cs(icol,j,k)=msqx_cs(icol,j,k)
          enddo
      endif
      msq(j,k)=msqx_cs(0,j,k)+msqx_cs(1,j,k)+msqx_cs(2,j,k)
   19 continue
      enddo
      enddo
      endif

      if (Qflag) then

        do j=-nf,nf
        do k=-nf,nf

        do i=0,2
        mqq(i,j,k)=0d0
        enddo

        if ((j .gt. 0) .and. (k .gt. 0)) then
c----QQ case
          call msq_qq(j,k,qR_a,qR_b,qq_a,qq_b,prop,msq0,msq1,msq2)

          msqx(0,j,k,j,k)=faclo*msq0
          msqx(1,j,k,j,k)=faclo*msq1
          msqx(2,j,k,j,k)=faclo*msq2

          mqq(0,j,k)=msqx(0,j,k,j,k)
          mqq(1,j,k)=msqx(1,j,k,j,k)
          mqq(2,j,k)=msqx(2,j,k,j,k)

          if (j .ne. k) then
            call msq_qq(j,k,qq_a,qq_b,qR_a,qR_b,prop,msq0,msq1,msq2)
            msqx(0,j,k,k,j)=faclo*msq0
            msqx(1,j,k,k,j)=faclo*msq1
            msqx(2,j,k,k,j)=faclo*msq2
         endif

          elseif ((j .lt. 0) .and. (k .lt. 0)) then
c----QbQb case
          call msq_qq(-j,-k,qbRb_a,qbRb_b,qbqb_a,qbqb_b,prop,
     .                msq0,msq1,msq2)

          msqx(0,j,k,j,k)=faclo*msq0
          msqx(1,j,k,j,k)=faclo*msq1
          msqx(2,j,k,j,k)=faclo*msq2

          mqq(0,j,k)=msqx(0,j,k,j,k)
          mqq(1,j,k)=msqx(1,j,k,j,k)
          mqq(2,j,k)=msqx(2,j,k,j,k)

          if (j .ne. k) then
            call msq_qq(-j,-k,qbqb_a,qbqb_b,qbRb_a,qbRb_b,prop,
     .                  msq0,msq1,msq2)

            msqx(0,j,k,k,j)=faclo*msq0
            msqx(1,j,k,k,j)=faclo*msq1
            msqx(2,j,k,k,j)=faclo*msq2
          endif

        elseif ((j .gt. 0) .and. (k .lt. 0)) then
C---Q-Qb case
          call msq_qqb(j,k,qRb_a,qRb_b,qqb_a,qqb_b,prop,
     .                msq0,msq1,msq2,tup,tdo)

          msqx(0,j,k,j,k)=faclo*msq0
          msqx(1,j,k,j,k)=faclo*msq1
          msqx(2,j,k,j,k)=faclo*msq2

          mqq(0,j,k)=msqx(0,j,k,j,k)
          mqq(1,j,k)=msqx(1,j,k,j,k)
          mqq(2,j,k)=msqx(2,j,k,j,k)

          if (j .eq. -k) then
            do f=1,5,2
              if (f .ne. j) then
                msqx(0,j,k,f,-f)=zip
                msqx(1,j,k,f,-f)=faclo*tdo
                msqx(2,j,k,f,-f)=zip
              endif
            enddo
            do f=2,4,2
              if (f .ne. j) then
                msqx(0,j,k,f,-f)=zip
                msqx(1,j,k,f,-f)=faclo*tup
                msqx(2,j,k,f,-f)=zip
              endif
           enddo
            if ((j.eq.1).or.(j.eq.3).or.(j.eq.5)) then
              nup=2
              ndo=nf-3
            else
              nup=1
              ndo=nf-2
            endif
           mqq(1,j,k)=mqq(1,j,k)+faclo*(dfloat(nup)*tup+dfloat(ndo)*tdo)
         endif

         call msq_qqb(j,k,qRb_ax,qRb_bx,qqb_ax,qqb_bx,prop,
     .             msq0,msq1,msq2,tup,tdo)

         msqx(0,j,k,k,j)=faclo*msq0
         msqx(1,j,k,k,j)=faclo*msq1
         msqx(2,j,k,k,j)=faclo*msq2

          if (j .eq. -k) then
            do f=1,5,2
              if (f .ne. j) then
                msqx(0,j,k,-f,f)=zip
                msqx(1,j,k,-f,f)=faclo*tdo
                msqx(2,j,k,-f,f)=zip
              endif
            enddo
            do f=2,4,2
              if (f .ne. j) then
                msqx(0,j,k,-f,f)=zip
                msqx(1,j,k,-f,f)=faclo*tup
                msqx(2,j,k,-f,f)=zip
              endif
            enddo
         endif

        elseif ((j .lt. 0) .and. (k .gt. 0)) then
C---Qb-Q case
          call msq_qqb(-j,-k,qbR_a,qbR_b,qbq_a,qbq_b,prop,
     .                msq0,msq1,msq2,tup,tdo)

          msqx(0,j,k,k,j)=faclo*msq0
          msqx(1,j,k,k,j)=faclo*msq1
          msqx(2,j,k,k,j)=faclo*msq2

          mqq(0,j,k)=msqx(0,j,k,k,j)
          mqq(1,j,k)=msqx(1,j,k,k,j)
          mqq(2,j,k)=msqx(2,j,k,k,j)

          if (j .eq. -k) then
            do f=1,5,2
              if (f .ne. k) then
                msqx(0,j,k,f,-f)=zip
                msqx(1,j,k,f,-f)=faclo*tdo
                msqx(2,j,k,f,-f)=zip
              endif
            enddo
            do f=2,4,2
              if (f .ne. k) then
                msqx(0,j,k,f,-f)=zip
                msqx(1,j,k,f,-f)=faclo*tup
                msqx(2,j,k,f,-f)=zip
              endif
           enddo
            if ((k.eq.1).or.(k.eq.3).or.(k.eq.5)) then
              nup=2
              ndo=nf-3
            else
              nup=1
              ndo=nf-2
            endif
           mqq(1,j,k)=mqq(1,j,k)+faclo*(dfloat(nup)*tup+dfloat(ndo)*tdo)
         endif

         call msq_qqb(-j,-k,qbR_ax,qbR_bx,qbq_ax,qbq_bx,prop,
     .             msq0,msq1,msq2,tup,tdo)

         msqx(0,j,k,j,k)=faclo*msq0
         msqx(1,j,k,j,k)=faclo*msq1
         msqx(2,j,k,j,k)=faclo*msq2

          if (j .eq. -k) then
            do f=1,5,2
              if (f .ne. k) then
                msqx(0,j,k,-f,f)=zip
                msqx(1,j,k,-f,f)=faclo*tdo
                msqx(2,j,k,-f,f)=zip
              endif
            enddo
            do f=2,4,2
              if (f .ne. k) then
                msqx(0,j,k,-f,f)=zip
                msqx(1,j,k,-f,f)=faclo*tup
                msqx(2,j,k,-f,f)=zip
              endif
            enddo
         endif

      endif

      msq(j,k)=msq(j,k)+mqq(0,j,k)+mqq(1,j,k)+mqq(2,j,k)

      enddo
      enddo

      endif

  999 continue
c--- restore proper colourchoice if necessary
      if ((part .eq. 'real') .or. (part .eq. 'virt')) then
        colourchoice=rcolourchoice
      endif
c--- restore proper parton sub-process selection, if necessary
      if ((part .eq. 'real') .and. (Qflag .eqv. .true.)) then
        Gflag=rGflag
      endif

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
PROGRAM driver
  use srHNF_brute
  use symmetry, only: get_lattice_pointGroup
  use vector_matrix_utilities
  use num_types
  use utilities
  use io_utils
  use numerical_utilities

  implicit none

  integer :: i, sfound, nfound, status, j
  integer, allocatable :: srHNFs(:,:,:), isrHNFs(:,:,:), temp_storage(:,:,:)
  real(dp), pointer :: pg(:,:,:)
  real(dp), allocatable :: gens(:,:,:)

  integer nMin, nMax ! Numbers of various things
  integer k ! Number of colors/label types (i.e., binary, ternary, etc.)
  integer LatDim ! 2D or 3D parent lattice?
  integer nD ! Number of sites in the basis (i.e., number of points in the multilattice)
  real(dp), pointer :: d(:,:) => null()
  character(1) :: latTyp
  real(dp)  eps, tstart, tend
  real(dp) :: parLV(3,3), invLV(3,3), abgens(3,3)

  character(80) title, fname
  logical fullLab,concCheck
  integer, pointer :: cRange(:,:)
  integer, pointer :: label(:,:)
  integer, pointer :: digit(:)
  integer, pointer :: equivalencies(:)

  if (iargc()>=1) then
     call getarg(1,fname)
  else
     fname = "struct_enum.in"
  endif
  call read_input(title,LatDim,parLV,nD,d,k,equivalencies,nMin,nMax,eps&
       &,fullLab,label,digit,fname,cRange,concCheck) ! Read in parent lattice vectors, etc.

  ! lat = reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
  ! lat = reshape((/0.5,0.5,0.0,0.5,0.0,0.5,0.0,0.5,0.5/),(/3,3/))
  ! lat = reshape((/0.5,0.5,-0.5,0.5,-0.5,0.5,-0.5,0.5,0.5/),(/3,3/))
  ! vol = 1000

  call cpu_time(tstart)
  ! Need to replace get_lattice_pointGroup with get space group
  call get_lattice_pointGroup(parLV,pg)
  call find_pg_gens(pg,gens,1e-10_dp)

  ! print *, "pg size", size(pg,3)
  ! print *, "lat vecs", parLV

  do i=1, size(gens,3)
     ! print *, "cbas gen", i
     ! do j=1,3
     !    print *, gens(j,:,i)
     ! end do
     print *, "latbas gen", i
     call matrix_inverse(parLV,invLV)
     abgens = matmul(invLV,matmul(gens(:,:,i),parLV))
     do j=1,3
        print *, abgens(j,:)
     end do
  end do

end PROGRAM driver

max_stars_repo_name: glwhart/autoGR

max_issues_repo_name: glwhart/autoGR

max_forks_repo_name: glwhart/autoGR


text:
submodule (potential_mumps) potential2d

use grid, only: gridflag
use calculus, only : grad2D1_curv_alt, grad2D3, grad2D3_curv_periodic
use PDEelliptic, only: elliptic2D_polarization,elliptic2D_polarization_periodic,elliptic2D_cart

implicit none (type, external)

contains

module procedure potential2D_polarization
!! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS
!! INCLUDES FULL OF POLARIZATION CURRENT, INCLUDING CONVECTIVE
!! TERMS.  VELOCITIES SHOULD BE TRIMMED (WITHOUT GHOST CELLS).
!! THIS VERSION OF THE *INTEGRATED* POTENTIAL SOLVER OBVIATES
!! ALL OTHERS SINCE A PURELY ELECTRSTATIC FORM CAN BE RECOVERED
!! BY ZEROING OUT THE INERTIAL CAPACITANCE.
!!
!! THIS FORM IS INTENDED TO  WORK WITH CURVILINEAR MESHES.
!! NOTE THAT THE FULL GRID VARIABLES (X%DX3ALL, ETC.) MUST
!! BE USED HERE!!!

real(wp), dimension(1:size(SigP2,1),1:size(SigP2,2)) :: gradSigH2,gradSigH3
integer :: u
integer :: lx2,lx3

lx2=x%lx2all    !use full grid sizes
lx3=x%lx3all

!gradSigH2=grad2D1(SigH,x,1,lx2)   !x2 is now 1st index and x3 is second...  This one appears to be the problem.  This issue here is that grad2D1 automatically uses x%dx1 as the differential element...
gradSigH2=grad2D1_curv_alt(SigH,x,1,lx2)
!! note the alt since we need to use dx2 as differential...  Tricky bug/feature
gradSigH3=grad2D3(SigH,x,1,lx3)
!! awkward way of handling this special case derivative which uses x3 as the differential to operate on a 2D array.

potential2D_polarization=elliptic2D_polarization(srcterm,SigP2,SigP3,SigH,gradSigH2,gradSigH3,Cm,v2,v3, &
                                  Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,x%dx1,x%dx1i,x%dx2all,x%dx2iall, &
                                  x%dx3all,x%dx3iall,Phi0,perflag,it)

end procedure potential2D_polarization

module procedure potential2D_polarization_periodic

!! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS
!! INCLUDES FULL OF POLARIZATION CURRENT, INCLUDING CONVECTIVE
!! TERMS.  VELOCITIES SHOULD BE TRIMMED (WITHOUT GHOST CELLS).
!! THIS VERSION OF THE *INTEGRATED* POTENTIAL SOLVER OBVIATES
!! ALL OTHERS SINCE A PURELY ELECTRSTATIC FORM CAN BE RECOVERED
!! BY ZEROING OUT THE INERTIAL CAPACITANCE.
!!
!! THIS FORM IS INTENDED TO  WORK WITH CARTESIAN MESHES ONLY.
!! NOTE THAT THE FULL GRID VARIABLES (X%DX3ALL, ETC.) MUST
!! BE USED HERE!!!
!!
!! THIS FUNCTION WORKS ON A PERIODIC MESH BY USING A CIRCULANT MATRIX

real(wp), dimension(1:size(SigP,1),1:size(SigP,2)+1) :: gradSigH2,gradSigH3

integer :: lx2,lx3

lx2=x%lx2all    !use full grid sizes
lx3=x%lx3all

!ZZZ - THESE NEED TO BE CHANGED INTO CIRCULAR/PERIODIC DERIVATIVES FOR THE X3 DIRECTION
gradSigH2=grad2D1_curv_alt(SigH,x,1,lx2)   !note the alt since we need to use dx2 as differential...  Tricky bug/feature
gradSigH3=grad2D3_curv_periodic(SigH,x,1,lx3)    !circular difference

potential2D_polarization_periodic=elliptic2D_polarization_periodic(srcterm,SigP,SigH,gradSigH2,gradSigH3, &
                          Cm,v2,v3,Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,x%dx1,x%dx1i,x%dx2all,x%dx2iall,x%dx3all, &
                          x%dx3iall,Phi0,perflag,it)

end procedure potential2D_polarization_periodic

module procedure potential2D_fieldresolved
!! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 2D USING MUMPS
!! ASSUME THAT WE ARE RESOLVING THE POTENTIAL ALONG THE FIELD
!! LINE AND THAT IT VARIES IN X1 AND X3 (X2 IS NOMINALL JUST
!! ONE ELEMENT.  LEFT AND RIGHT BOUNDARIES (IN X3) ARE ASSUMED
!! TO USE DIRICHLET BOUNARY CONDITIONS, WHILE THE (ALTITUDE)
!! TOP CAN BE NEUMANN OR DIRICHLET.  BOTTOM (ALTITUDE)
!! IS ALWAYS ASSUMED TO BE DIRICHLET.

real(wp), dimension(size(Vmaxx1,1),size(Vmaxx1,2)) :: Vmaxx1alt      !in case we need to do some transformations to convert current into potential

integer :: lx1,lx2all,lx3all

lx1=size(sig0,1)
lx2all=size(sig0,2)
lx3all=size(sig0,3)

if (flagdirich==0) then      !convert current into potential normal derivative
  Vmaxx1alt=-1*Vmaxx1*x%h1all(lx1,1:lx2all,1:lx3all)/sig0(lx1,1:lx2all,1:lx3all)
!  print*, 'Transforming boundary from current...'
!  print*, shape(Vmaxx1),shape(x%h1all),shape(sig0)
else                         !Dirichlet boundary conditions, don't change
  Vmaxx1alt=Vmaxx1
end if

potential2D_fieldresolved=elliptic2D_cart(srcterm,sig0,sigP,Vminx1,Vmaxx1alt,Vminx3,Vmaxx3, &
                     x%dx1,x%dx1i,x%dx3all,x%dx3iall,flagdirich,perflag,gridflag,it)

end procedure potential2D_fieldresolved

end submodule potential2d

max_stars_repo_name: jklenzing/gemini

max_issues_repo_name: jklenzing/gemini

max_forks_repo_name: jklenzing/gemini


text:
      SUBROUTINE PHO_SETPDF(Idpdg,Ityp,Ipar,Iset,Iext,Ipaval,Mode)
C***************************************************************
C
C     assigns  PDF numbers to particles
C
C     input:      IDPDG    PDG number of particle
C                 ITYP     particle type
C                 IPAR     PDF paramertization
C                 ISET     number of set
C                 IEXT     library number for PDF calculation
C                 IPAVAL   (only output)
C                          1 PDF with valence quarks
C                          0 PDF without valence quarks
C                 MODE     -1   add entry to table
C                           1   read from table
C                           2   map particle IPAR to existing entry
C                           2   output of table
C
C***************************************************************
      IMPLICIT NONE
      INTEGER i , idcmp , Idpdg , ientry , Iext , ii , Ipar , Ipaval ,
     &        ipdfs , Iset , Ityp , ityp1 , Mode
      SAVE

C  input/output channels
      INCLUDE 'inc/poinou'
C  event debugging information
      INCLUDE 'inc/podebg'
C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
      INCLUDE 'inc/pohdfl'

      DIMENSION ipdfs(5,50)
      DATA ientry/0/

      IF ( Mode.EQ.1 ) THEN
         i = 1
         IF ( Idpdg.EQ.81 ) THEN
            idcmp = IDEqp(1)
            Ipaval = IHFls(1)
         ELSE IF ( Idpdg.EQ.82 ) THEN
            idcmp = IDEqp(2)
            Ipaval = IHFls(2)
         ELSE
            idcmp = Idpdg
            Ipaval = 1
         END IF
 50      IF ( idcmp.EQ.ipdfs(1,i) ) THEN
            Ityp = ipdfs(2,i)
            Ipar = ipdfs(3,i)
            Iset = ipdfs(4,i)
            Iext = ipdfs(5,i)
            IF ( LPRi.GT.4 .AND. IDEb(80).GE.15 )
     &            WRITE (LO,'(1X,A,I7,5X,3I4)')
     &            'PHO_SETPDF: ID,IPAR,ISET,IEXT' , idcmp , Ipar ,
     &           Iset , Iext
            RETURN
         END IF
         i = i + 1
         IF ( i.GT.ientry ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,I7)')
     &            'PHO_SETPDF: no PDF assigned to ' , idcmp
            CALL PHO_ABORT
         END IF
         GOTO 50

      ELSE IF ( Mode.EQ.2 ) THEN
         ii = 0
         DO i = 1 , ientry
            ii = ii + 1
            IF ( Idpdg.EQ.ipdfs(1,i) ) GOTO 100
         END DO

 100     IF ( ii.LE.ientry ) THEN
            ientry = ientry + 1
            ipdfs(1,ientry) = Ipar
            ipdfs(2,ientry) = ipdfs(2,ii)
            ipdfs(3,ientry) = ipdfs(3,ii)
            ipdfs(4,ientry) = ipdfs(4,ii)
            ipdfs(5,ientry) = ipdfs(5,ii)
            IF ( (LPRi.GT.4) .AND. (IDEb(80).GT.0) )
     &            WRITE (LO,'(1X,A,I7,A,I7)') 'PHO_SETPDF: Particle ' ,
     &           Ipar , ' mapped to ' , Idpdg
         ELSE
            WRITE (LO,'(1X,A,I7,A,I7)') 'PHO_SETPDF: PDF of particle' ,
     &             Ipar , 'unknown.'
            CALL PHO_ABORT
         END IF

      ELSE IF ( Mode.EQ.-1 ) THEN
         DO i = 1 , ientry
            IF ( Idpdg.EQ.ipdfs(1,i) ) THEN
               IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,5I6)')
     &               'PHO_SETPDF: overwrite old particle PDF' , Idpdg ,
     &              ipdfs(2,i) , ipdfs(3,i) , ipdfs(4,i) , ipdfs(5,i)
               GOTO 150
            END IF
         END DO
         i = ientry + 1
         IF ( i.GT.50 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,/1x,6I6)')
     &            'PHO_SETPDF:ERROR: no space left in IPDFS:' , i ,
     &           Idpdg , ipdfs(2,i) , ipdfs(3,i) , ipdfs(4,i) ,
     &           ipdfs(5,i)
            STOP
         END IF
         ientry = i
 150     ipdfs(1,i) = Idpdg
         IF ( Idpdg.EQ.990 ) THEN
            ityp1 = 20
         ELSE IF ( Idpdg.EQ.22 ) THEN
            ityp1 = 3
         ELSE IF ( ABS(Idpdg).LT.1000 ) THEN
            ityp1 = 2
         ELSE
            ityp1 = 1
         END IF
         ipdfs(2,i) = ityp1
         ipdfs(3,i) = Ipar
         ipdfs(4,i) = Iset
         ipdfs(5,i) = Iext
      ELSE IF ( Mode.EQ.-2 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A)')
     &         'PHO_SETPDF: PDFs assigned by user:'
         DO i = 1 , ientry
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,I4,A,I7,A,4I5)') i ,
     &           '  particle:' , ipdfs(1,i) , '   PDF-set  ' ,
     &           ipdfs(2,i) , ipdfs(3,i) , ipdfs(4,i) , ipdfs(5,i)
         END DO
      ELSE
         IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,I5)')
     &         'PHO_SETPDF:ERROR: invalid mode ' , Mode
      END IF
      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


module utility_deallocate

    use define_indexes
    use define_globals
    use define_settings, only: setting

    implicit none

!%-------------------------------------------------------------------------
!% Description:
!%   This has all the array and variables deallocation.
!%   All the top-level storage should be deallocated in this module
!%-------------------------------------------------------------------------

    private

    integer           :: deallocation_status
    character(len=99) ::              emsg

    public :: util_deallocate_network_data
    public :: util_deallocate_partitioning_arrays
    public :: util_deallocate_check

contains
!%
!%==========================================================================
!% PUBLIC
!%==========================================================================
!%
    subroutine util_deallocate_network_data()
        !%------------------------------------------------------------------
        !% Description:
        !%   deallocate all the data relate to defined network, including:
        !%   linkX, nodeX, elemX, faceX
        !%-------------------------------------------------------------------
        !% Declarations:
            character(64) :: subroutine_name = 'util_deallocate_network_data'
        !%-------------------------------------------------------------------
        !% Preliminaries
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        !print *, 'call util_deallocate_linknode() ',this_image()
        call util_deallocate_linknode()

        !print *, 'call util_deallocate_elemX_faceX() ', this_image()
        call util_deallocate_elemX_faceX()

        !print *, 'call util_deallocate_elem_boundary_ghost() ', this_image()
        call util_deallocate_elem_boundary_ghost()

        !print *, 'call util_deallocate_columns() ', this_image()
        call util_deallocate_columns()

        !print *, 'call util_deallocate_bc() ', this_image()
        call util_deallocate_bc()

        !%-------------------------------------------------------------------
        !% closing
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_network_data
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_partitioning_arrays()
        !%-------------------------------------------------------------------
        !% Description
        !%   deallocates the partitioning arrays storage used for storing partitioned
        !%   network information
        !%-------------------------------------------------------------------
        !% Declarations
            character(64) :: subroutine_name = 'util_deallocate_partitioning_arrays'
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        deallocate(adjacent_links, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg,'adjacent_links')

        deallocate(elem_per_image, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elem_per_image')

        deallocate(image_full, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'image_full')

        !%-------------------------------------------------------------------
        !% closing
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_partitioning_arrays
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_check(deallocation_status, emsg, locationstring)
        !%-------------------------------------------------------------------
        !% Description:
        !%   Checks deallocation status and stops if there is an error
        !%-------------------------------------------------------------------
        !% Declaratoins
            integer,            intent(in   ) :: deallocation_status
            character(len=*),   intent(in   ) :: emsg
            character(len=*),   intent(in   ) :: locationstring !% unique identifier of location
            character(64):: subroutine_name = 'util_deallocate_check'
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%utility) &
            write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

            if (deallocation_status > 0) then
                print *, trim(emsg)
                print *, 'variable trying to deallocate = ',trim(locationstring)
                stop
            end if

        !%-------------------------------------------------------------------
            if (setting%Debug%File%utility) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_check
!%
!%==========================================================================
!% PRIVATE
!%==========================================================================
!%
    subroutine util_deallocate_linknode()
        !%-------------------------------------------------------------------
        !% Description:
        !%   deallocates the link and node storage used for the coarse representation
        !%   of the network connectivity
        !%
        !% Method:
        !%   Release the memory used by the tables node%I, link%I, node%R, link%R, node%YN,
        !%   link%YN.
        !%   These are defined in globals.f08), and allocated in allocate_storage.f08
        !%   Every time memory is deallocated, the utility_check_deallocation functionality
        !%   (from utility.f08) is used to determine wheter or not there was an error during
        !%   the deallocation.
        !%-------------------------------------------------------------------
        !% Declarations
            integer       :: ii
            character(64) :: subroutine_name = 'util_deallocate_linknode'
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        deallocate(node%I, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'node%I')

        deallocate(link%I, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'link%I')

        deallocate(node%R, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'node%R')

        deallocate(node%YN, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'node%YN')

        if (allocated(node%P%have_flowBC)) then
            deallocate(node%P%have_flowBC, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check(deallocation_status, emsg, 'node%P%Have_flowBC')
        end if

        if (allocated(node%P%have_headBC)) then
            deallocate(node%P%have_headBC, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check(deallocation_status, emsg, 'node%P%have_headBC')
        end if

        !% Deallocate link/node names
        ! do ii = 1, size(link%Names)
        !     deallocate(link%Names(ii)%str, stat=deallocation_status, errmsg=emsg)
        !     call util_deallocate_check(deallocation_status, emsg, 'link%Names(ii)%str')
        ! end do
        do ii = 1, size(node%Names)
            deallocate(node%Names(ii)%str, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check(deallocation_status, emsg, 'node%Names(ii)%str')
        end do

        deallocate(node%Names, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'node%Names')

        ! deallocate(link%Names, stat=deallocation_status, errmsg=emsg)
        ! call util_deallocate_check(deallocation_status, emsg, 'link%Names')

        !% brh20211217 --
        !% the followin commented asobsolete approach to storing output nodes and links
        !deallocate(node_output_idx,stat=deallocation_status,errmsg=emsg)
        !call util_deallocate_check(deallocation_status, emsg, 'node_output_idx')

        !deallocate(link_output_idx,stat=deallocation_status,errmsg=emsg)
        !call util_deallocate_check(deallocation_status, emsg, 'link_output_idx')

        ! if (allocated(node%P%have_output)) then
        !     deallocate(node%P%have_output,stat=deallocation_status,errmsg=emsg)
        !     call util_deallocate_check(deallocation_status, emsg, 'node%P%have_output')
        ! end if

        ! if (allocated(link%P%have_output)) then
        !     deallocate(link%P%have_output,stat=deallocation_status,errmsg=emsg)
        !     call util_deallocate_check(deallocation_status, emsg, 'link%P%have_output')
        ! end if

        !%-------------------------------------------------------------------
        !% closing
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_linknode
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_elemX_faceX ()
        !%-------------------------------------------------------------------
        ! Description:
        !   Simply deallocate the elemX and faceX thatwe assigned across all employed images
        !%-------------------------------------------------------------------
        !% Declarations:
            character(64) :: subroutine_name = 'util_deallocate_elemX_faceX'
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        !==== elem deallocation ====
        !print *, 'elemR'
        deallocate(elemR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemR')

        !print *, 'elemI'
        deallocate(elemI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemI')

        !print *, 'elemYN'
        deallocate(elemYN, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemYN')

        !print *, 'elemP'
        deallocate(elemP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemP')

        !print *, 'elemPGetm'
        deallocate(elemPGetm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemPGetm')

        !print *, 'elemPGac'
        deallocate(elemPGac, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemPGac')

        !print *, 'elemPGalltm'
        deallocate(elemPGalltm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemPGalltm')

        !print *, 'elemSI'
        deallocate(elemSI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemSI')

        !print *, 'elemSR'
        deallocate(elemSR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'elemSR')

        !==== face deallocation ====
        !print *, 'faceR'
        deallocate(faceR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'faceR')

        !print *, 'faceI'
        deallocate(faceI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'faceI')

        !print *, 'faceYN'
        deallocate(faceYN, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'faceYN')

        !print *, 'faceP'
        deallocate(faceP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'faceP')

        !%-------------------------------------------------------------------
        !% closing
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_elemX_faceX
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_elem_boundary_ghost ()
        !%-------------------------------------------------------------------
        ! Description:
        !   Simply deallocate the elemB%R and elemGR that we assigned across all employed images
        !%-------------------------------------------------------------------
        !% Declarations:
            character(64) :: subroutine_name = 'util_deallocate_elem_boundary_ghost'
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        !==== elemB deallocation ====
        if (num_images() > 1) then
            deallocate(elemB%R, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'elemB%R')

            deallocate(elemB%I, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'elemB%I')

            deallocate(elemGR, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check(deallocation_status, emsg, 'elemGR')
        end if
        !%-------------------------------------------------------------------
        !% closing
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_elem_boundary_ghost
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_columns()
        !%-------------------------------------------------------------------
        !% Description:
        !%   All variables are stored in col_elemX(:) arrays, release the memory here
        !%   in all images
        !%-------------------------------------------------------------------
            character(64) :: subroutine_name = 'util_deallocate_columns'
        !%-------------------------------------------------------------------
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        !==== col_elemI====
        deallocate(col_elemI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemI')
        !==== col_elemP====
        deallocate(npack_elemP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'npack_elemP')

        deallocate(col_elemP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemP')
        !==== col_elemPG====
        deallocate(npack_elemPGalltm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'npack_elemPGalltm')

        deallocate(col_elemPGalltm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemPGalltm')

        deallocate(npack_elemPGac, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'npack_elemPGac')

        deallocate(col_elemPGac, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemPGac')

        deallocate(npack_elemPGetm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'npack_elemPGetm')

        deallocate(col_elemPGetm, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemPGetm')
        !==== col_elemR ====
        deallocate(col_elemR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemR')
        !==== col_elemSI ====
        deallocate(col_elemSI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemSI')
        !==== col_elemSR ====
        deallocate(col_elemSR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemSR')
        !==== col_elemSGR ====
        deallocate(col_elemSGR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemSGR')
        !==== col_elemWDI ====
        deallocate(col_elemWDI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemWDI')
        !==== col_elemWDR ====
        deallocate(col_elemWDR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemWDR')
        !==== col_elemYN ====
        deallocate(col_elemYN, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_elemYN')
        !==== col_faceI ====
        deallocate(col_faceI, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_faceI')
        !==== col_faceM ====
       ! deallocate(col_faceM, stat=deallocation_status, errmsg=emsg)
       ! call util_deallocate_check(deallocation_status, emsg, 'col_faceM')
        !==== col_faceP ====
        deallocate(npack_faceP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'npack_faceP')

        deallocate(col_faceP, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_faceP')
        !==== col_faceR ====
        deallocate(col_faceR, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_faceR')
        !==== col_faceYN ====
        deallocate(col_faceYN, stat=deallocation_status, errmsg=emsg)
        call util_deallocate_check(deallocation_status, emsg, 'col_faceYN')

        !%-------------------------------------------------------------------
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_columns
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine util_deallocate_bc()
        !%-------------------------------------------------------------------
            character(64) :: subroutine_name = 'util_deallocate_bc'
        !%-------------------------------------------------------------------
            !if (crashYN) return
            if (setting%Debug%File%utility_deallocate) &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------

        if (N_flowBC > 0) then
            deallocate(BC%flowI, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%flowI')

            deallocate(BC%flowR_timeseries, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%flowR_timeseries')

            deallocate(BC%flowIdx, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%flowIdx')

            deallocate(BC%flowRI, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%flowRI')

            if (allocated(BC%P%BCup)) then
                deallocate(BC%P%BCup, stat=deallocation_status, errmsg=emsg)
                call util_deallocate_check (deallocation_status, emsg, 'BC%P%BCup')
            end if

            if (allocated(BC%P%BClat)) then
                deallocate(BC%P%BClat, stat=deallocation_status, errmsg=emsg)
                call util_deallocate_check (deallocation_status, emsg, 'BC%P%BClat')
            end if
        end if

        if (N_headBC > 0) then
            deallocate(BC%headI, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%headI')

            deallocate(BC%headR_timeseries, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%headR_timeseries')

            deallocate(BC%headIdx, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%headIdx')

            deallocate(BC%headRI, stat=deallocation_status, errmsg=emsg)
            call util_deallocate_check (deallocation_status, emsg, 'BC%headRI')

            if (allocated(BC%P%BCdn)) then
                deallocate(BC%P%BCdn, stat=deallocation_status, errmsg=emsg)
                call util_deallocate_check (deallocation_status, emsg, 'BC%P%BCdn')
            end if
        end if
        !%-------------------------------------------------------------------
        if (setting%Debug%File%utility_deallocate) &
        write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine util_deallocate_bc
!%
!%==========================================================================
!% END MODULE
!%==========================================================================
!%
end module utility_deallocate


text:
!-----------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in    !
!  continuous development by various groups and is based on information !
!  from these groups: Federal Government employees, contractors working !
!  within a United States Government contract, and non-Federal sources  !
!  including research institutions.  These groups give the Government   !
!  permission to use, prepare derivative works of, and distribute copies!
!  of their work in the CMAQ system to the public and to permit others  !
!  to do so.  The United States Environmental Protection Agency         !
!  therefore grants similar permission to use the CMAQ system software, !
!  but users are requested to provide copies of derivative works or     !
!  products designed to operate in the CMAQ system to the United States !
!  Government without restrictions as to use by others.  Software       !
!  that is used with the CMAQ system but distributed under the GNU      !
!  General Public License or the GNU Lesser General Public License is   !
!  subject to their copyright restrictions.                             !
!-----------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/yoj/arc/JPROC/src/driver/jproc_table/o3scal.f,v 1.5 2011/10/29 01:03:53 sjr Exp $

C what(1) key, module and SID; SCCS file; date and time of last delta:
C @(#)o3scal.F	1.1 /project/mod3/JPROC/src/driver/jproc_table/SCCS/s.o3scal.F 23 May 1997 12:44:20

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       SUBROUTINE O3SCAL ( O3, HO3, XLAT, DOBNEW )

C*********************************************************************
C
C  adjustment of O3 profiles to a user-selected dobson value.
C   select value of DOBNEW in main program
C   if don't want to use, don't call this subroutine
C
C*********************************************************************

      IMPLICIT NONE

      INCLUDE 'JVALPARMS.EXT'    ! jproc parameters

C...........ARGUMENTS and their descriptions

      REAL         O3( MXLEV )        ! ozone profile
      REAL         HO3                ! ozone scale height
      REAL         XLAT               ! latitudinal band
      REAL         DOBNEW             ! new dobson unit

C...........LOCAL VARIABLES and their descriptions:

      LOGICAL, SAVE :: FIRSTTIME = .TRUE.

      INTEGER      I                  ! level index

      REAL         DOBSREF            ! reference dobson unit
      REAL         SCALE_FACTOR

C*********************************************************************
C     begin body of subroutine O3SCAL2

      IF ( FIRSTTIME ) THEN
        FIRSTTIME = .FALSE.
        WRITE ( 6, 2001 )
      END IF

C...convert O3 value into DU

      DOBSREF = O3( MXLEV ) * 1.0E5 * HO3
      DO I = 1, MXLEV
        DOBSREF = DOBSREF + O3( I ) * 1.0E5
      END DO
      DOBSREF = DOBSREF / 2.687E16

C...apply scale factor to O3 profile if a non-zero TOC value provided

      IF ( DOBNEW .GT. 0.0 ) THEN
        SCALE_FACTOR = DOBNEW / DOBSREF
        WRITE( 6, 2003 ) XLAT, DOBSREF, DOBNEW, SCALE_FACTOR
        DO I = 1, MXLEV
          O3( I ) = O3( I ) * SCALE_FACTOR
        END DO
      ELSE
        WRITE( 6, 2005 ) XLAT, DOBSREF
      END IF

C...formats

2001  FORMAT( 1X, 'Rescaling O3 profile based on TOC data' )
2003  FORMAT( 7X, 'Latitude:', 1X, F5.1, ', oldDOBS=', F9.4, ', newDOBS=',
     &        F9.4, ', scale factor=', F7.4 )
2005  FORMAT( 7X, 'Latitude:', 1X, F5.1, ', DOBS=', F9.4, ', not rescaled' )

      RETURN
      END

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
    ! This module does post-processing of .data files. For example to importance sample new data,
    ! correct approximate theory (eg. add in lensing), or to compute missing theory (e.g. matter power).

    module ImportanceSampling
    use settings
    use GeneralTypes
    use CalcLike
    use ParamPointSet
    use IO
    implicit none
    private

    Type, extends(TTheoryLikelihoodUser) :: TImportanceSampler
        logical  redo_likelihoods, redo_theory
        real(mcp) :: redo_skip = 100
        character(LEN=:), allocatable :: redo_datafile, redo_outroot,redo_like_name
        real(mcp) :: redo_likeoffset = 0
        real(mcp) :: redo_temperature = 1
        integer :: redo_thin = 1
        logical :: redo_change_like_only = .false.
        !This last one is for comparing goodness of fit
        !After importance sampling, you can recompute the likelihoods without the new data, but
        !keeping the weights from the importance sampling, and thereby asses whether the mean
        !likelihood wrt the original distribution of the parameter space after importance sampling
        !is similar to that after, in which case the datasets intersect in a region of high likelihood

        logical :: redo_nochange = .false.
        !for adding likelihood derived parameters without affecting original distribution

        logical :: redo_add = .false.
        !if just want to add new datasets rather than re-computing the entire likelihood

        logical :: redo_from_text = .false.
        !Redo from text files if .data files not available

        logical :: redo_no_new_data  = .false. !true to make no new .data files to save space

        logical ::  redo_output_txt_theory = .false. !output theory data as text for each point (e.g. for rainbow plots)
        character(LEN=:), allocatable :: redo_output_txt_root

        logical :: redo_auto_likescale = .true. !If big difference is log-likelihood, automatically rescale to O(1) weights
        real(mcp) :: redo_max_logLike_diff = 10
        integer :: redo_auto_likescale_count = 5 !number to check weights before rescaling
    contains
    procedure :: ReadParams => TImportanceSampler_ReadParams
    procedure :: Init => TImportanceSampler_Init
    procedure :: ImportanceSample => TImportanceSampler_ImportanceSample
    end Type TImportanceSampler

    !not supported any more    logical :: txt_theory = .false. !True to put P_k in output chains

    public TImportanceSampler
    contains

    subroutine TImportanceSampler_ReadParams(this, Ini)
    class(TImportanceSampler) :: this
    class(TSettingIni) :: Ini

    this%redo_likelihoods = Ini%Read_Logical('redo_likelihoods')
    this%redo_theory = Ini%read_Logical('redo_theory')
    this%redo_datafile = Ini%Read_String('redo_datafile')
    this%redo_outroot = Ini%Read_String('redo_outroot')
    this%redo_like_name = Ini%Read_String('redo_like_name')

    call Ini%Read('redo_likeoffset',this%redo_likeoffset)
    call Ini%Read('redo_temp',this%redo_temperature)
    call Ini%Read('redo_change_like_only',this%redo_change_like_only)
    call Ini%Read('redo_nochange',this%redo_nochange)
    call Ini%Read('redo_add',this%redo_add)
    call Ini%Read('redo_from_text',this%redo_from_text)
    call Ini%Read('redo_no_new_data',this%redo_no_new_data)
    call Ini%Read('redo_skip',this%redo_skip)
    call Ini%Read('redo_thin',this%redo_thin,min=1)
    call Ini%Read('redo_output_txt_theory',this%redo_output_txt_theory)
    call Ini%Read('redo_auto_likescale',this%redo_auto_likescale)
    call Ini%Read('redo_max_logLike_diff',this%redo_max_logLike_diff)
    call Ini%Read('redo_auto_likescale_count',this%redo_auto_likescale_count)

    if (this%redo_from_text .and. (this%redo_add .or. this%redo_like_name/='')) &
        call Mpistop('redo_add and/or redo_like_name require .data files, not from text')

    if (this%redo_outroot == '') then
        this%redo_outroot =  File%ExtractPath(baseroot)//'post_'//File%ExtractName(baseroot)
    end if
    if (this%redo_output_txt_theory) this%redo_output_txt_root = Ini%Read_String('redo_output_txt_dir')

    end subroutine TImportanceSampler_ReadParams

    subroutine TImportanceSampler_Init(this, LikeCalculator)
    class(TImportanceSampler) :: this
    class(TLikeCalculator), target:: LikeCalculator

    select type (LikeCalculator)
    class is (TTheoryLikeCalculator)
        this%LikeCalculator => LikeCalculator
        class default
        call MpiStop('Importance sampling requires TTheoryLikeCalculator')
    end select

    end subroutine TImportanceSampler_Init

    subroutine TImportanceSampler_ImportanceSample(this,InputFile)
    class(TImportanceSampler) :: this
    character(LEN=*), intent(INOUT):: InputFile
    real(mcp) truelike,mult,like, like_diff
    real(mcp) weight_min, weight_max, mult_sum, mult_ratio, mult_max,weight
    real(mcp) max_like, max_truelike
    integer error,num, debug
    character (LEN=:), allocatable :: post_root, data_point_txt_root
    integer i
    Type (ParamSet), allocatable :: Params
    logical :: has_likes(DataLikelihoods%Count)
    class(TDataLikelihood), pointer :: DataLike
    logical :: first, has_chain = .true.
    integer(File_size_int)  last_file_loc
    integer :: at_beginning, num_used, redo_loop
    integer thin_acc
    class(TFileStream), pointer :: InF
    Type(TTextFile), target :: InChain
    Type(TBinaryFile), target :: OutData, InData
#ifdef MPI
    integer ierror
    real(mcp), allocatable :: like_diffs(:)
#endif

    flush_write = .false.
    debug = 0

    this%LikeCalculator%Temperature = this%redo_temperature

    if (Feedback>0 .and. this%redo_change_like_only) &
        write (*,*) 'Warning: only changing likelihoods not weights'

    if (Feedback>0 .and. this%redo_nochange) &
        write (*,*) 'Warning: calcalating new likelihoods but not using them'

    if (this%redo_datafile /= '') InputFile = this%redo_datafile

    if (this%redo_from_text) then
        call InChain%Open(trim(InputFile)//'.txt')
        if (this%redo_skip < 1) then
            this%redo_skip = nint(InChain%Lines()*this%redo_skip)
        end if
        InF => InChain
        if (.not. this%redo_theory) write (*,*) '**You probably want to set redo_theory**'
    else
        if (File%Exists(trim(InputFile)//'.data')) then
            call InData%Open(trim(InputFile)//'.data')
            InF => InData
        else
            write(*,*) 'Chain .data files does not exist: ', MpiRank+1
            has_chain =.false.
        end if
    end if

    post_root = this%redo_outroot

    if (MpiRank==0 .and.BaseParams%NameMapping%nnames/=0) then
        call BaseParams%OutputParamNames(post_root,params_used, add_derived=.true.)
        call BaseParams%OutputParamRanges(post_root)
    end if

    if (has_chain) then
        if (instance /= 0) post_root = numcat(post_root//'_',instance)

        if (Feedback > 0) then
            if (this%redo_from_text) then
                write (*,*) 'reading from: ' //  trim(InputFile)//'.txt'
            else
                write (*,*) 'reading from: ' //  trim(InputFile)//'.data'
            end if
            write (*,*) 'writing to: ' // trim(post_root)//'.*'
        end if

        if (this%redo_output_txt_theory) then
            if (this%redo_output_txt_root =='') then
                this%redo_output_txt_root  =  post_root
            else
                this%redo_output_txt_root =  File%Join(this%redo_output_txt_root, File%ExtractName(post_root))
            end if
            write (*,*) 'Writing text file data to ' // this%redo_output_txt_root
        end if

        write (*,*) 'Using temperature: ', this%LikeCalculator%Temperature

        redo_loop= 1
        do
            call InF%Rewind()
            call ChainOutFile%CreateFile(trim(post_root)//'.txt')
            if (.not. this%redo_no_new_data) call OutData%CreateFile(trim(post_root)//'.data')

            num = 0
            num_used = 0
            weight_min= 1e30_mcp
            weight_max = -1e30_mcp
            mult_sum = 0
            mult_ratio = 0
            mult_max = -1e30_mcp
            max_like = logZero
            max_truelike = logZero
            at_beginning=0
            first = .true.

            allocate(Params)
            call this%LikeCalculator%Config%NewTheory(Params%Theory)

            do
                if (this%redo_from_text) then
                    error = 0
                    Params%P(:num_params)= BaseParams%center
                    if (.not. IO_ReadChainRow(InChain, mult, like, Params%P, params_used)) exit
                    num=num+1
                    if (this%redo_skip>=1 .and. num<=this%redo_skip) cycle
                else
                    call Params%ReadModel(InF,has_likes, mult,like, error)
                    num=num+1
                    if (first .and. this%redo_like_name/='') then
                        first=.false.
                        do i=1, DataLikelihoods%Count
                            DataLike => DataLikelihoods%Item(i)
                            if (DataLike%name==this%redo_like_name) then
                                if (.not. has_likes(i)) &
                                    call MpiStop('does not currently have like named:'//trim(this%redo_like_name))
                                has_likes(i)=.true.
                                if (any(.not. has_likes)) call MpiStop('not all other likelihoods exist already')
                                has_likes(i)=.false.
                                this%redo_add =.true.
                                exit
                            end if
                        end do
                    end if
                    if (this%redo_skip>0.d0 .and. this%redo_skip<1) then
                        at_beginning=at_beginning+1
                        if (at_beginning==1) then
                            last_file_loc = InF%Position()
                            cycle
                        elseif (at_beginning==2) then
                            this%redo_skip = InF%Size()/(InF%Position() -last_file_loc) * this%redo_skip
                            if (Feedback > 0) print *,'skipping ',nint(this%redo_skip), ' models'
                        end if
                    else if (num<=this%redo_skip) then
                        cycle
                    end if
                end if

                if (this%redo_thin>1) then
                    if (abs(nint(mult) - mult) > 1e-4) &
                        call MpiStop('redo_thin can only be used with chains with integer weights')
                    thin_acc = thin_acc + nint(mult)
                    if (thin_acc >= this%redo_thin) then
                        mult = thin_acc / this%redo_thin
                        thin_acc = mod(thin_acc, this%redo_thin)
                    else
                        cycle
                    end if
                end if

                if (error ==1) then
                    if (num==0) call MpiStop('Error reading data file.')
                    exit
                end if

                if (this%redo_likelihoods .or. this%redo_add) then
                    !Check for new prior before calculating anything
                    if (this%LikeCalculator%CheckPriorCuts(Params)==logZero) then
                        if (Feedback >1) write(*,*) 'Model outside new prior bounds: skipped'
                        cycle
                    end if
                end if

                if (this%redo_theory) then
                    call this%LikeCalculator%GetTheoryForImportance(Params, error)
                else
                    error = 0
                end if

                if (error ==0) then
                    if (this%redo_likelihoods .or. this%redo_add) then
                        call this%LikeCalculator%UpdateTheoryForLikelihoods(Params)
                        if (this%redo_add) then
                            truelike = this%LikeCalculator%GetLogLikePost(Params, .not. has_likes)
                        else
                            truelike = this%LikeCalculator%GetLogLikePost(Params)
                        end if
                        if (truelike == logZero) then
                            weight = 0
                        else
                            weight = exp(like-truelike+this%redo_likeoffset)
                        end if

                        if (.not. this%redo_change_like_only .and. .not. this%redo_nochange)  mult = mult*weight
                    else
                        truelike = like
                        weight = 1
                    end if

                    if (this%redo_nochange) truelike = like

                    max_like = min(max_like,like)
                    max_truelike = min(max_truelike,truelike)

                    num_used=num_used+1

                    mult_ratio = mult_ratio + weight
                    mult_sum = mult_sum + mult

                    if (this%redo_auto_likescale .and. redo_loop==1 .and. num_used == this%redo_auto_likescale_count &
                        & .and. .not. this%redo_change_like_only .and. .not. this%redo_nochange) then
                    !Check log likelihoods scaled to give sensible weights. Rescale must be constant between chains
                    if (max_truelike /= logZero) then
                        like_diff = max_truelike - max_like
                    else
                        like_diff = logZero
                    end if
#ifdef MPI
                    allocate(like_diffs(MPIchains))
                    call MPI_Allgather(like_diff, 1, MPI_real_mcp, like_diffs, 1,  MPI_real_mcp, MPI_COMM_WORLD, ierror)
                    if (all(like_diffs==logZero)) then
                        like_diff=0
                    else
                        like_diff = sum(like_diffs, mask = like_diffs/=logZero) / count(like_diffs/=logZero)
                    end if
#endif
                    if (abs(like_diff) > this%redo_max_logLike_diff) then
                        this%redo_likeoffset = like_diff
                        if (Feedback > 0 .and. MpiRank==0) write(*,*) 'Re-starting with redo_likeoffset = ',this%redo_likeoffset
                        redo_loop=2
                        exit
                    end if
                    end if

                    if (mult /= 0) then
                        if (this%redo_output_txt_theory) then
                            data_point_txt_root = this%redo_output_txt_root // '_'//IntToStr(num)
                            call this%LikeCalculator%WriteParamPointTextData(data_point_txt_root, Params)
                            call this%LikeCalculator%WriteParamsHumanText(data_point_txt_root//'.pars', Params, truelike, weight)
                        end if
                        call Params%WriteParams(this%LikeCalculator%Config,mult,truelike)
                        if (.not. this%redo_no_new_data) call Params%WriteModel(OutData, truelike,mult)
                    else
                        if (Feedback >1 ) write (*,*) 'Zero weight: new like = ', truelike
                    end if

                    if (Feedback > 1) write (*,*) num, ' mult= ', real(mult), ' weight = ', real(weight)
                    weight_max = max(weight,weight_max)
                    weight_min = min(weight,weight_min)
                    mult_max = max(mult_max,mult)
                end if

            end do

            call OutData%Close()
            call ChainOutFile%Close()
            deallocate(Params)
            if (redo_loop/=2) exit
            redo_loop=3
        end do
        call InF%Close()

        if (Feedback>0) then
            write(*,*) 'finished. Processed ',num_used,' models'
            write (*,*) 'max weight= ',weight_max, ' min weight = ',weight_min
            write (*,*) 'mean mult  = ', mult_sum/num_used
            write (*,*) 'mean importance weight (approx evidence ratio) = ',mult_ratio/num_used
            write (*,*) 'effective number of samples =',mult_sum/mult_max
            write (*,*) 'Best redo_likeoffset = ',max_truelike - max_like
        end if

        if ((mult_ratio < 1e-6 .or. mult_ratio > 1e8) .and. .not.this%redo_change_like_only) then
            write (*,*) 'WARNING: use redo_likeoffset to rescale likelihoods'
        end if

    end if

    end subroutine TImportanceSampler_ImportanceSample

    end module ImportanceSampling

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD


text:
!
! Copyright (c) 2017, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! Sourced allocate of array of derived types.
program P
  implicit none
  integer, parameter :: SIZE1 = 5, SIZE2 = 3
  type A
    integer :: a
  end type A
  type B
    integer, allocatable :: bb
    type(A), allocatable :: b
  end type B

  call test_1_dim()
  call test_2_dim()

contains

  subroutine test_1_dim()
    integer :: i, n
    integer :: expected(SIZE1), actual(SIZE1)
    type(B), allocatable :: x(:), y(:)
    allocate(x(SIZE1))
    do i = 1, SIZE1
      allocate(x(i)%b)
      x(i)%b%a = i
    end do

    allocate(y, source=x)

    n = 0
    do i = 1, SIZE1
      n = n + 1
      expected(n) = x(i)%b%a
      x(i)%b%a = 0  ! should not change y
      actual(n) = y(i)%b%a
    end do
    call check(actual, expected, SIZE1)
  end subroutine

  subroutine test_2_dim()
    integer :: i, j, n
    integer :: expected(SIZE1*SIZE2), actual(SIZE1*SIZE2)
    type(B), allocatable :: x(:,:), y(:,:)

    allocate(x(SIZE1,SIZE2))
    do i = 1, SIZE1
      do j = 1, SIZE2
        allocate(x(i,j)%b)
        x(i,j)%b%a = 10*i + j
      end do
    end do

    allocate(y, source=x)

    n = 0
    do i = 1, SIZE1
      do j = 1, SIZE2
        n = n + 1
        expected(n) = x(i,j)%b%a
        x(i,j)%b%a = 0  ! should not change y
        actual(n) = y(i,j)%b%a
      end do
    end do
    call check(actual, expected, SIZE1*SIZE2)
  end subroutine

end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
module Mod_nsm_Output
   use MPI
   use Mod_iofile
   use Mod_int2str
   use Mod_PointerSetter
   use Mod_nsm_BaseBouope
   implicit none
   private
   public SetPointersOutput

   type, extends(PointerSetter) :: SPOutput
contains
      procedure :: SpecificSet => SpecificSetOutput
   end type
   type(SPOutput) :: SetPointersOutput

   real(rp), allocatable, dimension(:,:) :: force, momen
   integer(ip) :: ierr

contains

   subroutine SpecificSetOutput(d)
      implicit none
      class(SPOutput) :: d
      if (a%kfl_outfm==1) then
         if(a%kfl_openforce) call ConcatenateProcedures(ProcHook_Initializations,OpenFiles)
         call ConcatenateProcedures(ProcHook_Initializations,AllocForces)
         call ConcatenateProcedures(ProcHook_Initializations,ToZero)
         call ConcatenateProcedures(ProcHook_InGaussElmatsAssembly,ForceGauss)
         call ConcatenateProcedures(ProcHook_PostLoop,ComputeForce)
         call ConcatenateProcedures(ProcHook_PostLoop,OutputForce)
         call ConcatenateProcedures(ProcHook_Finalizations,DeallocForces)
      end if
   end subroutine

   subroutine AllocForces
      implicit none

      call a%Mesh%GetNbody(nbody)
      call a%Memor%alloc(ndime,nbody,force,'force','nsi_forces')
      call a%Memor%alloc(3,nbody,momen,'momen','nsi_forces')
   end subroutine

   subroutine DeallocForces
      implicit none

      call a%Mesh%GetNbody(nbody)
      call a%Memor%dealloc(ndime,nbody,force,'force','nsi_forces')
      call a%Memor%dealloc(3,nbody,momen,'momen','nsi_forces')
   end subroutine

   subroutine ToZero
      implicit none

      a%force=0.0_rp
      a%momen=0.0_rp
      force=0.0_rp
      momen=0.0_rp
   end subroutine

   subroutine OpenFiles
      implicit none
      character(150) :: fil_force

      call a%Mesh%GetNbody(nbody)
      if (a%MPIrank == a%MPIroot) then
         do ibody=1,nbody
            fil_force = trim(a%OutputFolder)//'/'//adjustl(trim(a%namda))//adjustl(trim(int2str(a%MPIrank)))//'.'//adjustl(trim(a%exmod))//adjustl(trim(int2str(ibody)))//'.frc'
            call iofile(zero,a%lun_force(ibody),fil_force,adjustl(trim(a%exmod))//'FORCES & MOMENTS')
            if (ndime==2) write (a%lun_force(ibody),12)
            if (ndime==3) write (a%lun_force(ibody),13)
         end do
         a%kfl_openforce = .false.
      endif

     !Formats
     12   format( &
         '$        Forces and Moments over bodies          ',/,          &
         '$    Time',15x,'Fx',15x,'Fy',15x,'Mz')
     13   format( &
         '$          Forces and Moments over bodies            ',/,      &
         '$ Time',15x,'Fx',15x,'Fy',15x,'Fz',15x,'Mx',15x,'My'15x,'Mz')
   end subroutine

   subroutine ForceGauss
      implicit none

      if (ibody>0) then
         !Forces
         do idime=1,ndime
            force(idime,ibody)=force(idime,ibody) + tract(idime)*dsurf*weightfactor
         end do

         !Moments
         if (ndime==2) then
            momen(3,ibody) = momen(3,ibody) &
               + tract(2)*(gpcod(1)-a%origm(1))*dsurf*weightfactor &
               - tract(1)*(gpcod(2)-a%origm(2))*dsurf*weightfactor
         elseif (ndime==3) then
            momen(1,ibody) = momen(1,ibody) &
               + tract(3)*(gpcod(2)-a%origm(2))*dsurf*weightfactor &
               - tract(2)*(gpcod(3)-a%origm(3))*dsurf*weightfactor
            momen(2,ibody) = momen(1,ibody) &
               + tract(1)*(gpcod(3)-a%origm(3))*dsurf*weightfactor &
               - tract(3)*(gpcod(1)-a%origm(1))*dsurf*weightfactor
            momen(3,ibody) = momen(3,ibody) &
               + tract(2)*(gpcod(1)-a%origm(1))*dsurf*weightfactor &
               - tract(1)*(gpcod(2)-a%origm(2))*dsurf*weightfactor
         end if
      end if
   end subroutine

   subroutine ComputeForce
      implicit none
      if (nbody>0) then
        !call MPI_ALLREDUCE(force,a%force,ndime*nbody,MPI_REAL8,MPI_SUM,a%MPIroot,a%MPIcomm,ierr)
        !call MPI_ALLREDUCE(momen,a%momen,3*nbody,MPI_REAL8,MPI_SUM,a%MPIroot,a%MPIcomm,ierr)
        call MPI_REDUCE(force,a%force,ndime*nbody,MPI_REAL8,MPI_SUM,a%MPIroot,a%MPIcomm,ierr)
        call MPI_REDUCE(momen,a%momen,3*nbody,MPI_REAL8,MPI_SUM,a%MPIroot,a%MPIcomm,ierr)
      end if
   end subroutine

   subroutine OutputForce
      implicit none

      call a%Mesh%GetNbody(nbody)
      if (a%MPIrank == a%MPIroot) then
         if (ndime==2) then
            do ibody = 1,nbody
               write(a%lun_force(ibody),20) a%ctime,&
                     -a%force(1,ibody)/a%adimf(1),  &
                     -a%force(2,ibody)/a%adimf(2),  &
                     -a%momen(3,ibody)/a%adimm(3)
            end do
         else if(ndime==3) then
            do ibody = 1,nbody
               write(a%lun_force(ibody),30) a%ctime,&
                     -a%force(1,ibody)/a%adimf(1),  &
                     -a%force(2,ibody)/a%adimf(2),  &
                     -a%force(3,ibody)/a%adimf(3),  &
                     -a%momen(1,ibody)/a%adimm(1),  &
                     -a%momen(2,ibody)/a%adimm(2),  &
                     -a%momen(3,ibody)/a%adimm(3)
            end do
         end if
      end if

      if (a%kfl_flush == 1) then
         do ibody = 1,nbody
            call flush(a%lun_force(ibody))
         end do
      endif

      !Formats
      20   format(1x,e12.6,2x,3(2x,e15.8))
      30   format(1x,e12.6,2x,6(2x,e15.8))
   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      SUBROUTINE JULDAT ( JUL, ISTYLE, CDATE, NDATE)
C
C     Converts an HEC style julian date (days since Dec 31, 1899),
C     Into a character date of various styles, as shown below
C
C     Input:
C        JUL:  The julian date, in days since Dec 31, 1899
C     Output:
C        CDATE:  A character variable to contain the date (should be
C                long enough to hold the date
C        NDATE:  The number of characters in CDATE
C
C
C  ISTYLE  Form   ISTYLE   Form      ISTYLE   Form      ISTYLE   Form
C      LC, 4 CH YR       LC, 2 CH YR     UC, 4 CH YR      UC, 2 CH YR
C  0:  June 2, 1985  10:  June 2, 85  100:  JUNE 2, 1985  110:  JUNE 2, 85
C  1:  Jun 2, 1985   11:  Jun 2, 85   101:  JUN 2, 1985   111:  JUN 2, 85
C  2:  2 June 1985   12:  2 June 85   102:  2 JUNE 1985   112:  2 JUNE 85
C  3:  June 1985     13:  June 85     103:  JUNE 1985     113:  JUNE 85
C  4:  02Jun1985     14:  02Jun85     104:  02JUN1985     114:  02JUN85
C  5:  2Jun1985      15:  2Jun85      105:  2JUN1985      115:  2JUN85
C  6:  Jun1985       16:  Jun85       106:  JUN1985       116:  JUN85
C  7:  02 Jun 1985   17:  02 Jun 85   107:  02 JUN 1985   117:  02 JUN 85
C  8:  2 Jun 1985    18:  2 Jun 85    108:  2 JUN 1985    118:  2 JUN 85
C  9:  Jun 1985      19:  Jun 85      109:  JUN 1985      119:  JUN 85
C
C     ISTYLE=-1:  CDATE = 6/2/85       ISTYLE=-11:  CDATE = 06/02/85
C     ISTYLE=-2:  CDATE = 6-2-85       ISTYLE=-12:  CDATE = 06-02-85
C
C     If ISTYLE is zero, it defaults to style 1.
C
C
      CHARACTER CDATE*(*)
      INTEGER*4 JUL                                                     MLu
C
      I = JLIYMD ( JUL, IYR, IMON, IDAY)
      IF (I.LT.0) THEN
         call strcpy(CDATE, ' ')
         NDATE = 1
         RETURN
      ENDIF
      CALL YMDDAT ( IYR, IMON, IDAY, ISTYLE, CDATE, NDATE, IERR)
	IF (IERR.LT.0) THEN
         call strcpy(CDATE, ' ')
         NDATE = 1
      ENDIF
C
      RETURN
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


PROGRAM main

  USE mo_kind, ONLY: dp, sp
  USE mo_fit,  ONLY: fitfun, fpoly_dp, fpoly_sp, polyfit, svdfit
  USE mo_fit,  ONLY: polyval
  use mo_ansi_colors, only: color, c_red, c_green

  IMPLICIT NONE

  LOGICAL                                :: isgood

  integer, parameter                     :: degree = 2
  integer                                :: i

  real(dp), dimension(11)                :: dx = (/ (i,i=0,10) /)
  real(dp), dimension(11)                :: dy = (/ 1,   6,  17,  34, 57,  86, 121, 162, 209, 262, 321 /)
  real(dp), dimension(11)                :: dsig
  real(dp), dimension(degree+1)          :: da, dw
  real(dp), dimension(degree+1,degree+1) :: dv
  real(dp)                               :: dchisq
  real(dp), dimension(2)                 :: dla
  real(dp)                               :: xdp = 3
  real(dp), dimension(4)                 :: xdpvec = (/ 0, 1, 2, 3 /)
  real(dp), dimension(3)                 :: polydp = (/ 1, 2, 1 /)

  real(sp), dimension(11)                :: sx = (/ (i,i=0,10) /)
  real(sp), dimension(11)                :: sy = (/ 1,   6,  17,  34, 57,  86, 121, 162, 209, 262, 321 /)
  real(sp), dimension(11)                :: ssig
  real(sp), dimension(degree+1)          :: sa, sw
  real(sp), dimension(degree+1,degree+1) :: sv
  real(sp)                               :: schisq
  real(sp), dimension(2)                 :: sla
  real(sp)                               :: xsp = 3
  real(sp), dimension(4)                 :: xspvec = (/ 0, 1, 2, 3 /)
  real(sp), dimension(3)                 :: polysp = (/ 1, 2, 1 /)

  Write(*,*) ''
  Write(*,*) 'Test mo_fit.f90'
  Write(*,*) ''
  Write(*,*) 'Output should read 1, 2, 3'

  ! Double precision
  Write(*,*) ''
  isgood = .true.
  da = polyfit(dx, dy, degree)
  write(*,*) da
  isgood = isgood .and. all((nint(1000._dp*da)-(/1000, 2000, 3000/)) == 0)
  dsig = 0.1
  call fitfun(dx, dy, dsig, da, fpoly_dp)
  write(*,*) da
  isgood = isgood .and. all((nint(1000._dp*da)-(/1000, 2000, 3000/)) == 0)
  call svdfit(dx, dy, dsig, da, dv, dw, dchisq, fpoly_dp)
  write(*,*) da
  isgood = isgood .and. all((nint(1000._dp*da)-(/1000, 2000, 3000/)) == 0)
  Write(*,*) ''
  Write(*,*) 'Output should read -44, 32'
  Write(*,*) ''
  dla = polyfit(dx, dy, 1)
  write(*,*) dla
  Write(*,*) ''
  Write(*,*) 'Output should read 16 & 1 4 9 16'
  Write(*,*) ''
  write(*,*) polyval( polydp, xdp)
  write(*,*) polyval( polydp, xdpvec)
  isgood = isgood .and. all((nint(1000._dp*polyval( polydp, xdpvec))-(/1000, 4000, 9000, 16000/)) == 0)

  if (isgood) then
     write(*,*) 'mo_fit double precision ', color('o.k.', c_green)
  else
     write(*,*) 'mo_fit double precision ', color('failed!', c_red)
  endif

  ! Single precision
  Write(*,*) ''
  isgood = .true.
  sa = polyfit(sx, sy, degree)
  write(*,*) sa
  isgood = isgood .and. all((nint(1000._sp*sa)-(/1000, 2000, 3000/)) == 0)
  ssig = 0.1
  call fitfun(sx, sy, ssig, sa, fpoly_sp)
  write(*,*) sa
  isgood = isgood .and. all((nint(1000._sp*sa)-(/1000, 2000, 3000/)) == 0)
  call svdfit(sx, sy, ssig, sa, sv, sw, schisq, fpoly_sp)
  write(*,*) sa
  isgood = isgood .and. all((nint(1000._sp*sa)-(/1000, 2000, 3000/)) == 0)
  Write(*,*) ''
  Write(*,*) 'Output should read -44, 32'
  Write(*,*) ''
  sla = polyfit(sx, sy, 1)
  write(*,*) sla
  Write(*,*) ''
  Write(*,*) 'Output should read 16 & 1 4 9 16'
  Write(*,*) ''
  write(*,*) polyval( polysp, xsp)
  write(*,*) polyval( polysp, xspvec)
  isgood = isgood .and. all((nint(1000._sp*polyval( polysp, xspvec))-(/1000, 4000, 9000, 16000/)) == 0)

  if (isgood) then
     write(*,*) 'mo_fit single precision ', color('o.k.', c_green)
  else
     write(*,*) 'mo_fit single precision ', color('failed!', c_red)
  endif

END PROGRAM main


text:
      PROGRAM main
C     THIS SUBROUTINE READS THE RUNSCRIPT AND ENTERS FILENAMES INTO DATA ARRAYS
C     FOR USE IN 'READS' AND 'READQ'. WHEN FINISHED THIS SUBROUTINE CALLS
C     'SOIL' WHICH IS THE MAIN SUBROUTINE FROM WHICH ALL OTHERS ARE CALLED
C
      include "parameters.h"
      include "filec.h"
      include "files.h"
      include "blkc.h"
      DIMENSION NA(250),ND(250)
      CHARACTER*16 DATA(30),DATAC(30,250,250),DATAP(JP,JY,JX)
     2,DATAM(JP,JY,JX),DATAX(JP),DATAY(JP),DATAZ(JP,JY,JX)
     3,OUTS(10),OUTP(10),OUTFILS(10,JY,JX),OUTFILP(10,JP,JY,JX)
      CHARACTER*3 CHOICE(102,20)
      CHARACTER*8 CDATE
      CHARACTER*80 BUF,PREFIX
      CALL GETCWD(BUF)
C
C     IDENTIFY OPERATING SYSTEM: DOS OR UNIX
C
      IF((.NOT.(BUF(1:1).EQ.'/'.OR.BUF(1:1).EQ.'~'))
     2.AND.BUF(2:2).EQ.':')THEN
      CALL GETARG(1,BUF)
      OPEN(5,FILE=BUF,STATUS='OLD')
      CALL GETARG(2,BUF)
      CALL CHDIR(BUF)
      PREFIX='.\\'
C     make output directory
      outdir=trim(buf)//'\\outputs\\'
      ELSE
      PREFIX='./'
C     make output directory
      outdir=trim(buf)//'/outputs/'
      ENDIF
      call system('mkdir -p '//trim(outdir))
C
C     READ INPUT FILES
C
10    FORMAT(A16)
      IGO=0
C
C     NUMBER OF COLUMNS AND ROWS
C
      READ(5,*)NHW,NVN,NHE,NVS
C
C     SITE FILE
C
      READ(5,10)DATA(1)
C
C     TOPOGRAPHY FILE
C
      READ(5,10)DATA(2)
C
C     READ THE NUMBER OF TIMES THE SCENARIOS IN THE MODEL RUN
C     ARE TO BE EXECUTED
C
100   READ(5,*)NAX,NDX
      IF(NAX.EQ.0.AND.NDX.EQ.0)GO TO 1000
C
C     NUMBER OF SCENES IN THE NEXT SCENARIO OF THE MODEL RUN
C     AND THE NUMBER OF TIMES THIS SCENARIO IS TO BE EXECUTED
C
      DO 105 NEX=1,NAX
      READ(5,*)NAY,NDY
      NA(NEX)=NAY
      ND(NEX)=NDY
C
C     FOR EACH SCENE IN THIS SCENARIO:
C
      DO 110 NE=1,NA(NEX)
C
C     WEATHER FILE
C
      READ(5,10,END=1000)DATAC(3,NE,NEX)
C
C     WEATHER OPTIONS
C
      READ(5,10)DATAC(4,NE,NEX)
C
C     LAND MANAGEMENT FILE
C
      READ(5,10)DATAC(9,NE,NEX)
C
C     PLANT MANAGEMENT FILE
C
      READ(5,10)DATAC(10,NE,NEX)
C
C     OUTPUT DATA CONTROL
C
      DO 115 N=21,30
      READ(5,10)DATAC(N,NE,NEX)
115   CONTINUE
110   CONTINUE
105   CONTINUE
C
C     RUN THIS SCRIPT
C
      DO 120 NTX=1,NDX
      DO 120 NEX=1,NAX
      DO 120 NT=1,ND(NEX)
      DO 120 NE=1,NA(NEX)
      CALL SOIL(NA,ND,NT,NE,NAX,NDX,NTX,NEX,NHW,NHE,NVN,NVS)
      IGO=IGO+1
120   CONTINUE
C
C     SCRIPT COMPLETED, START NEXT SCRIPT
C
      GO TO 100
1000  STOP
      END

max_stars_repo_name: ThemosTsikas/ECOSYS

max_issues_repo_name: ThemosTsikas/ECOSYS

max_forks_repo_name: ThemosTsikas/ECOSYS