text:
c
c Example analysis for "p p > t j [QCD]" process.
c Example analysis for "p p > t j $$ w+ w- [QCD]" process.
c Example analysis for "p p > w+ > t j  [QCD]" process.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_begin(nwgt,weights_info)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer nwgt
      character*(*) weights_info(*)
      integer j,kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      character*6 cc(2)
      data cc/'|T@NLO','|T@LO '/
      include 'dbook.inc'
      call inihist
      nwgt_analysis=nwgt
      if (nwgt_analysis*48.gt.nplots/4) then
         write (*,*) 'error in analysis_begin: '/
     &        /'too many histograms, increase NPLOTS to',
     &        nwgt_analysis*48*4
         stop 1
      endif
      do j=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*48+(j-1)*24
      call bookup(l+ 1,'t pt        '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+ 2,'t log pt    '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+ 3,'t y         '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+ 4,'t eta       '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      call bookup(l+ 5,'j1 pt       '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+ 6,'j1 log pt   '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+ 7,'j1 y        '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+ 8,'j1 eta      '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      call bookup(l+ 9,'j2 pt       '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+10,'j2 log pt   '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+11,'j2 y        '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+12,'j2 eta      '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      call bookup(l+13,'bj1 pt      '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+14,'bj1 log pt  '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+15,'bj1 y       '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+16,'bj1 eta     '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      call bookup(l+17,'bj2 pt      '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+18,'bj2 log pt  '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+19,'bj2 y       '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+20,'bj2 eta     '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      call bookup(l+21,'syst pt     '//cc(j)//weights_info(kk)
     &     ,5d0,0d0,200d0)
      call bookup(l+22,'syst log pt '//cc(j)//weights_info(kk)
     &     ,0.05d0,0d0,5d0)
      call bookup(l+23,'syst y      '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
      call bookup(l+24,'syst eta    '//cc(j)//weights_info(kk)
     &     ,0.25d0,-6d0,6d0)
c
      enddo
      enddo
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_end(xnorm)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      character*14 ytit
      double precision xnorm
      integer i
      integer kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      include 'dbook.inc'
      call open_topdrawer_file
      call mclear
      do i=1,NPLOTS
         call mopera(i,'+',i,i,xnorm,0.d0)
         call mfinal(i)
      enddo
      ytit='sigma per bin '
      do i=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*48+(i-1)*24
      call multitop(l+ 1,2,3,'t pt',' ','LOG')
      call multitop(l+ 2,2,3,'t log pt',' ','LOG')
      call multitop(l+ 3,2,3,'t y',' ','LOG')
      call multitop(l+ 4,2,3,'t eta',' ','LOG')
c
      call multitop(l+ 5,2,3,'j1 pt',' ','LOG')
      call multitop(l+ 6,2,3,'j1 log pt',' ','LOG')
      call multitop(l+ 7,2,3,'j1 y',' ','LOG')
      call multitop(l+ 8,2,3,'j1 eta',' ','LOG')
c
      call multitop(l+ 9,2,3,'j2 pt',' ','LOG')
      call multitop(l+10,2,3,'j2 log pt',' ','LOG')
      call multitop(l+11,2,3,'j2 y',' ','LOG')
      call multitop(l+12,2,3,'j2 eta',' ','LOG')
c
      call multitop(l+13,2,3,'bj1 pt',' ','LOG')
      call multitop(l+14,2,3,'bj1 log pt',' ','LOG')
      call multitop(l+15,2,3,'bj1 y',' ','LOG')
      call multitop(l+16,2,3,'bj1 eta',' ','LOG')
c
      call multitop(l+17,2,3,'bj2 pt',' ','LOG')
      call multitop(l+18,2,3,'bj2 log pt',' ','LOG')
      call multitop(l+19,2,3,'bj2 y',' ','LOG')
      call multitop(l+20,2,3,'bj2 eta',' ','LOG')
c
      call multitop(l+21,2,3,'syst pt',' ','LOG')
      call multitop(l+22,2,3,'syst log pt',' ','LOG')
      call multitop(l+23,2,3,'syst y',' ','LOG')
      call multitop(l+24,2,3,'syst eta',' ','LOG')
c
      enddo
      enddo
      call close_topdrawer_file
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      include 'nexternal.inc'
      include 'cuts.inc'
      integer istatus(nexternal)
      integer iPDG(nexternal)
      double precision p(0:4,nexternal)
      double precision wgts(*)
      integer ibody,mu,count_bj,count_j
      double precision wgt,var
      integer i,j,kk,l,nwgt_analysis,itop
      common/c_analysis/nwgt_analysis
      double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,pjet(0:3
     $     ,nexternal),ptj1,yj1,etaj1,ptbj1,ybj1,etabj1,ptj2,yj2,etaj2
     $     ,ptbj2,ybj2,etabj2,p_top(0:3),pttop,ytop,etatop,psyst(0:3)
     $     ,ptsyst,ysyst,etasyst,p_bjet(0:3,nexternal)
      integer nQCD,nb,is_b(nexternal),njet,jet(nexternal),nbjet
      logical is_b_jet(nexternal)
      double precision getpt,getrapidity,getpseudorap
      external getpt,getrapidity,getpseudorap
      double precision pi
      parameter (pi=3.14159265358979312d0)
      if (nexternal.ne.5) then

         write (*,*) 'error #1 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
         write (*,*) 'error #2 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
         write (*,*) 'error #3 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
         write (*,*) 'error #4 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(5)).le.5 .or. ipdg(5).eq.21)) then
         write (*,*) 'error #4 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif
      if (ipdg(3).ne.6) then
         write (*,*) 'error #5 in analysis_fill: '/
     &        /'only for process "p p > t j [QCD]"'
         stop 1
      endif

      nQCD=0
      nb=0
      do i=1,nexternal
         if(abs(ipdg(i)).le.5.or.ipdg(i).eq.21) then
            nQCD=nQCD+1
            do j=0,3
               pQCD(j,nQCD)=p(j,i)
            enddo
         endif
         if (abs(ipdg(i)).eq.5) then
            nb=nb+1
            is_b(nb)=nQCD
         endif
         if(abs(ipdg(i)).eq.6)itop=i
      enddo
      if(nb.ne.1)then
         write(*,*)'not one b',nb
         stop
      endif

c Define jet clustering parameters (from cuts.inc via the run_card.dat)
      palg=JETALGO              ! jet algorithm: 1.0=kt, 0.0=C/A, -1.0 = anti-kt
      rfj=JETRADIUS             ! the radius parameter
      sycut=PTJ                 ! minimum transverse momentum

c******************************************************************************
c     call FASTJET to get all the jets
c
c     INPUT:
c     input momenta:               pQCD(0:3,nexternal), energy is 0th component
c     number of input momenta:     nQCD
c     radius parameter:            rfj
c     minumum jet pt:              sycut
c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
c
c     OUTPUT:
c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
c     the number of jets (with pt > SYCUT):    njet
c     the jet for a given particle 'i':        jet(i),   note that this is
c     the particle in pQCD, which doesn't necessarily correspond to the particle
c     label in the process
c
      call amcatnlo_fastjetppgenkt(pQCD,nQCD,rfj,sycut,palg,pjet,njet,jet)
c
c******************************************************************************
      if(njet.gt.2)then
         write(*,*)'more than two jets',njet
         stop
      endif

C b-jet
      do i=1,njet
         is_b_jet(i)=.false.
         if(jet(is_b(nb)).eq.i)is_b_jet(i)=.true.
      enddo

      do i=0,3
         p_top(i)=p(i,itop)
      enddo
      pttop = getpt(p_top)
      ytop  = getrapidity(p_top)
      etatop= getpseudorap(p_top)

      count_j=0
      count_bj=0
      do i=1,njet
         if(.not.is_b_jet(i))then
            count_j=count_j+1
            if(count_j.eq.1)then
               ptj1 = getpt(pjet(0,i))
               yj1  = getrapidity(pjet(0,i))
               etaj1= getpseudorap(pjet(0,i))
               do mu=0,3
                  psyst(mu)=p_top(mu)+pjet(mu,i)
               enddo
               ptsyst = getpt(psyst)
               ysyst  = getrapidity(psyst)
               etasyst= getpseudorap(psyst)
            elseif(count_j.eq.2)then
               ptj2 = getpt(pjet(0,i))
               yj2  = getrapidity(pjet(0,i))
               etaj2= getpseudorap(pjet(0,i))
            endif
         elseif(is_b_jet(i))then
            count_bj=count_bj+1
            if (count_bj.eq.1) then
               ptbj1 = getpt(pjet(0,i))
               ybj1  = getrapidity(pjet(0,i))
               etabj1= getpseudorap(pjet(0,i))
            elseif (count_bj.eq.2) then
               ptbj2 = getpt(pjet(0,i))
               ybj2  = getrapidity(pjet(0,i))
               etabj2= getpseudorap(pjet(0,i))
            endif
         endif
      enddo
      nbjet=count_bj
c fill the histograms
      do i=1,2
         do kk=1,nwgt_analysis
            www=wgts(kk)
            l=(kk-1)*48+(i-1)*24
            if (ibody.ne.3 .and.i.eq.2) cycle
            call mfill(l+1,pttop,www)
            if(pttop.gt.0d0) call mfill(l+2,log10(pttop),www)
            call mfill(l+3,ytop,www)
            call mfill(l+4,etatop,www)
            if(njet.ge.1)then
               call mfill(l+5,ptj1,www)
               if (ptj1.gt.0d0) call mfill(l+6,log10(ptj1),www)
               call mfill(l+7,yj1,www)
               call mfill(l+8,etaj1,www)
               call mfill(l+21,ptsyst,www)
               if(ptsyst.gt.0d0) call mfill(l+22,log10(ptsyst),www)
               call mfill(l+23,ysyst,www)
               call mfill(l+24,etasyst,www)
            endif
            if(njet.ge.2)then
               call mfill(l+9,ptj2,www)
               if(ptj2.gt.0d0) call mfill(l+10,log10(ptj2),www)
               call mfill(l+11,yj2,www)
               call mfill(l+12,etaj2,www)
            endif
            if(nbjet.ge.1)then
               call mfill(l+13,ptbj1,www)
               if (ptbj1.gt.0d0) call mfill(l+14,log10(ptbj1),www)
               call mfill(l+15,ybj1,www)
               call mfill(l+16,etabj1,www)
            endif
            if(nbjet.ge.2)then
               call mfill(l+17,ptbj2,www)
               if (ptbj2.gt.0d0) call mfill(l+18,log10(ptbj2),www)
               call mfill(l+19,ybj2,www)
               call mfill(l+20,etabj2,www)
            endif
         enddo
      enddo

 999  return
      end

      function getrapidity(p)
      implicit none
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y,p(0:3)
      parameter (tiny=1.d-5)
c
      en=p(0)
      pl=p(3)
      xplus=en+pl
      xminus=en-pl
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
        if( (xplus/xminus).gt.tiny )then
          y=0.5d0*log( xplus/xminus )
        else
          y=sign(1.d0,pl)*1.d8
        endif
      else
        y=sign(1.d0,pl)*1.d8
      endif
      getrapidity=y
      return
      end

      function getpseudorap(p)
      implicit none
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th,p(0:3)
      parameter (tiny=1.d-5)
c
      en=p(0)
      ptx=p(1)
      pty=p(2)
      pl=p(3)
      pt=sqrt(ptx**2+pty**2)
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
        eta=sign(1.d0,pl)*1.d8
      else
        th=atan2(pt,pl)
        eta=-log(tan(th/2.d0))
      endif
      getpseudorap=eta
      return
      end

      function getpt(p)
      implicit none
      real*8 getpt,p(0:3)
      getpt=dsqrt(p(1)**2+p(2)**2)
      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
C     Scalar replacement in the presence of conditional control flow
C
C     Steve Carr and Ken Kennedy
C
C     Software - Practive and Experience, Vol. 24, No. 1, pp. 51-77, 1994
C
C     Sixth example, page 24: Livermore Loop 6 (or 11) cannot be handled
C     by PIPS because we do not exploit dependence arcs and because
C     regions fuse here two independent references

      subroutine scalarization36(b, w, n)
      real b(n,n), w(n)

      do i = 2, n
         do k = 1, i-1
            w(i) = w(i)+b(i,k)*w(i-k)
         enddo
      enddo

      print *, w

      end

ext: f


text:
subroutine MakeGravGradGridDH(cilm, lmax, gm, r0, a, f, vxx, vyy, vzz, vxy, &
                              vxz, vyz, n, sampling, lmax_calc, extend, &
                              exitstatus)
!------------------------------------------------------------------------------
!
!   Given the 4pi normalized gravitational spherical harmonic coefficients
!   CILM, this subroutine will compute 2D Driscol and Healy sampled grids of
!   the six components of the gravity "gradient" tensor in a local
!   north-oriented reference frame:
!
!       (Vxx,   Vxy,    Vxz)
!       (Vyx,   Vyy,    Vyz)
!       (Vzx,   Vzy,    Vzz)
!
!   where X points NORTH, Y points WEST, and Z points UPWARD. The gravitational
!   potential is defined as
!
!       V = GM/r Sum_{l=0}^LMAX (R0/r)^l Sum_{m=-l}^l C_{lm} Y_{lm},
!
!   Laplace's equation implies that Vxx + Vyy + Vzz = 0, and the gravity tensor
!   is symmetric. The components are calculated according to eq. 1 in
!   Petrovskaya and Vershkov (2006, J. Geod, 80, 117-127), which is based on
!   eq. 3.28 in Reed (1973, Ohio State Univ., Dept. Geod. Sci., Rep. 201,
!   Columbus, OH). Note that Reed's equations are in terms of latitude, and
!   that the Y axis points East:
!
!       Vzz = Vrr
!       Vxx = 1/r Vr + 1/r^2 Vtt
!       Vyy = 1/r Vr + 1/r^2 /tan(t) Vt + 1/r^2 /sin(t)^2 Vpp
!       Vxy = 1/r^2 /sin(t) Vtp - cos(t)/sin(t)^2 /r^2 Vp
!       Vxz = 1/r^2 Vt - 1/r Vrt
!       Vyz = 1/r^2 /sin(t) Vp - 1/r /sin(t) Vrp
!
!   where r, t, p stand for radius, theta, and phi, and subscripts on V denote
!   partial derivatives.
!
!   The output grid are in units of 1/s and are cacluated on a flattened
!   ellipsoid with semi-major axis A and flattening F. To obtain units of
!   Eotvos (10^-9 s^-1), multiply by 10^9.
!
!   When SAMPLING = 1, the output grids contain N samples in latitude from
!   90 to -90 + interval and N samples in longitude from 0 to 360-2*interval,
!   where N=2*(LMAX+1) and interval=180/N. When SAMPLING = 2, the grids are
!   equally spaced in degrees latitude and longitude with dimension (N x 2N).
!   If the optional parameter EXTEND is set to 1, the output grids will contain
!   an extra column corresponding to 360 E and an extra row corresponding to
!   90 S, which increases each of the dimensions of the grid by one.
!
!   Calling Parameters
!
!       IN
!           cilm        The 4-pi normalized gravitational potential spherical
!                       harmonic coefficients.
!           lmax        The maximum spherical harmonic degree of the function,
!                       used to determine the number of samples N.
!           GM          Product of the gravitatonal constant and the planet's
!                       mass.
!           r0          Reference radius of potential coefficients.
!           a           The semimajor axis of the flattened ellipsoid.
!           f           Flattening of the planet.
!
!       IN, OPTIONAL
!           sampling    (1) Grid is N latitudes by N longitudes (default).
!                       (2) Grid is N by 2N. The higher frequencies resulting
!                       from this oversampling in longitude are discarded, and
!                       hence not aliased into lower frequencies.
!           lmax_calc   The maximum spherical harmonic degree to evaluate
!                       the coefficients up to.
!           extend      If 1, return a grid that contains an additional column
!                       and row corresponding to 360 E longitude and 90 S
!                       latitude, respectively.
!
!       OUT
!           Vxx         x-x component of the gravity gradient tensor.
!           Vyy         y-y component of the gravity gradient tensor.
!           Vzz         z-z component of the gravity gradient tensor.
!           Vxy         x-y component of the gravity gradient tensor.
!           Vxz         x-z component of the gravity gradient tensor.
!           Vyz         y-z component of the gravity gradient tensor.
!           N           Number of samples in latitude. Number of samples in
!                       longitude is N when sampling is 1 (default), and is
!                       2N when sampling is 2.
!
!       OPTIONAL (OUT)
!           exitstatus  If present, instead of executing a STOP when an error
!                       is encountered, the variable exitstatus will be
!                       returned describing the error.
!                       0 = No errors;
!                       1 = Improper dimensions of input array;
!                       2 = Improper bounds for input variable;
!                       3 = Error allocating memory;
!                       4 = File IO error.
!
!   Notes:
!       1.  If lmax is greater than the the maximum spherical harmonic
!           degree of the input coefficients, then the coefficients will be
!           zero padded.
!       2.  Latitude is geocentric latitude.
!
!   Copyright (c) 2005-2019, SHTOOLS
!   All rights reserved.
!
!------------------------------------------------------------------------------
    use FFTW3
    use ftypes
    use, intrinsic :: iso_c_binding

    implicit none

    real(dp), intent(in) :: cilm(:,:,:), gm, r0, a, f
    real(dp), intent(out) :: vxx(:,:), vyy(:,:), vzz(:,:), vxy(:,:), vxz(:,:),&
                             vyz(:,:)
    integer, intent(in) :: lmax
    integer, intent(out) :: n
    integer, intent(in), optional :: sampling, lmax_calc, extend
    integer, intent(out), optional :: exitstatus
    integer :: l, m, i, l1, m1, lmax_comp, i_eq, i_s, astat(4), nlong, &
               nlat_out, nlong_out, extend_grid
    real(dp) :: grid(4*lmax+4), pi, theta, scalef, rescalem, u, p, dpl, dpl2, &
                dpl2s, pmm, sint, pm1, pm2, z, tempr, r_ex, lat, &
                prefactor(lmax), coefr0, coefrs0, coeft0, coefts0, coefp0, &
                coefps0, coefrr0, coefrrs0, coefrt0, coefrts0, coefrp0, &
                coefrps0, coeftp0, coeftps0, coefpp0, coefpps0, coeftt0, &
                coeftts0
    complex(dp) :: coef(2*lmax+3), coefr(2*lmax+3), coefrs(2*lmax+3), &
                   coeft(2*lmax+3), coefts(2*lmax+3), coefp(2*lmax+3), &
                   coefps(2*lmax+3), tempc, coefrr(2*lmax+3), &
                   coefrrs(2*lmax+3), coefrt(2*lmax+3), coefrts(2*lmax+3), &
                   coefrp(2*lmax+3), coefrps(2*lmax+3), coeftp(2*lmax+3), &
                   coeftps(2*lmax+3), coefpp(2*lmax+3), coefpps(2*lmax+3), &
                   coeftt(2*lmax+3), coeftts(2*lmax+3)
    type(C_PTR) :: plan
    real(dp), save, allocatable :: ff1(:,:), ff2(:,:), sqr(:)
    integer(int1), save, allocatable :: fsymsign(:,:)
    integer, save :: lmax_old = 0

!$OMP   threadprivate(ff1, ff2, sqr, fsymsign, lmax_old)

    if (present(exitstatus)) exitstatus = 0

    n = 2 * lmax + 2

    if (present(sampling)) then
        if (sampling == 1) then
            nlong = n
        else if (sampling == 2) then
            nlong = 2 * n
        else
            print*, "Error --- MakeGravGradGridDH"
            print*, "Optional parameter SAMPLING must be 1 (N by N) " // &
                    "or 2 (N by 2N)."
            print*, "Input value is ", sampling
            if (present(exitstatus)) then
                exitstatus = 2
                return
            else
                stop
            end if
        end if
    else
        nlong = n
    end if

    if (present(extend)) then
        if (extend == 0) then
            extend_grid = 0
            nlat_out = n
            nlong_out = nlong
        else if (extend == 1) then
            extend_grid = 1
            nlat_out = n + 1
            nlong_out = nlong + 1
        else
            print*, "Error --- MakeGravGradGridDH"
            print*, "Optional parameter EXTEND must be 0 or 1."
            print*, "Input value is ", extend
            if (present(exitstatus)) then
                exitstatus = 2
                return
            else
                stop
            end if
        end if
    else
        extend_grid = 0
        nlat_out = n
        nlong_out = nlong
    end if

    if (size(cilm(:,1,1)) < 2) then
        print*, "Error --- MakeGravGradGridDH"
        print*, "CILM must be dimensioned as (2, *, *)."
        print*, "Input dimension is ", size(cilm(:,1,1)), size(cilm(1,:,1)), &
                size(cilm(1,1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    end if

    if (size(vxx(:,1)) < nlat_out .or. size(vxx(1,:)) < nlong_out .or. &
            size(vyy(:,1)) < nlat_out .or. size(vyy(1,:)) < nlong_out .or. &
            size(vzz(:,1)) < nlat_out .or. size(vzz(1,:)) < nlong_out .or. &
            size(vxy(:,1)) < nlat_out .or. size(vxy(1,:)) < nlong_out .or. &
            size(vxz(:,1)) < nlat_out .or. size(vxz(1,:)) < nlong_out .or. &
            size(vyz(:,1)) < nlat_out .or. size(vyz(1,:)) < nlong_out) then
        print*, "Error --- MakeGravGradGridDH"
        print*, "VXX, VYY, VZZ, VXY, VXZ, and VYZ must be dimensioned " // &
                "as: ", nlat_out, nlong_out
        print*, "Input dimensions are ", size(vxx(:,1)), size(vxx(1,:)), &
            size(vyy(:,1)), size(vyy(1,:)), size(vzz(:,1)), size(vzz(1,:)), &
            size(vxy(:,1)), size(vxy(1,:)), size(vxz(:,1)), size(vxz(1,:)), &
            size(vyz(:,1)), size(vyz(1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    end if

    if (cilm(1,1,1) /= 1.0_dp .and. f /= 0.0_dp) then
        print*, "Warning --- MakeGravGradGridDH"
        print*, "The degree-0 term of the spherical harmonic coefficients is "
        print*, "not equal to 1. The variation in gravity resulting from "
        print*, "variations in radius of the flattened ellipsoid will not be "
        print*, "taken into account."
        print*, "C00 = ", cilm(1,1,1)
        print*, "F = ", f
    end if

    pi = acos(-1.0_dp)

    scalef = 1.0e-280_dp

    if (present(lmax_calc)) then
        if (lmax_calc > lmax) then
            print*, "Error --- MakeGravGradGridDH"
            print*, "LMAX_CALC must be less than or equal to LMAX."
            print*, "LMAX = ", lmax
            print*, "LMAX_CALC = ", lmax_calc
            if (present(exitstatus)) then
                exitstatus = 2
                return
            else
                stop
            end if

        else
            lmax_comp = min(lmax, size(cilm(1,1,:))-1, size(cilm(1,:,1))-1, &
                            lmax_calc)

        end if

    else
        lmax_comp = min(lmax, size(cilm(1,1,:))-1, size(cilm(1,:,1))-1)

    end if

    !--------------------------------------------------------------------------
    !
    !   Calculate recursion constants used in computing Legendre functions.
    !
    !--------------------------------------------------------------------------
    if (lmax_comp /= lmax_old) then

        if (allocated (sqr)) deallocate (sqr)
        if (allocated (ff1)) deallocate (ff1)
        if (allocated (ff2)) deallocate (ff2)
        if (allocated (fsymsign)) deallocate (fsymsign)

        allocate (sqr(2 * lmax_comp + 1), stat=astat(1))
        allocate (ff1(lmax_comp+1,lmax_comp+1), stat=astat(2))
        allocate (ff2(lmax_comp+1,lmax_comp+1), stat=astat(3))
        allocate (fsymsign(lmax_comp+1,lmax_comp+1), stat=astat(4))

        if (sum(astat(1:4)) /= 0) then
            print*, "Error --- MakeGravGradGridDH"
            print*, "Problem allocating arrays SQR, FF1, FF2, or FSYMSIGN", &
                    astat(1), astat(2), astat(3), astat(4)
            if (present(exitstatus)) then
                exitstatus = 3
                return
            else
                stop
            end if

        end if

        !----------------------------------------------------------------------
        !
        !   Calculate signs used for symmetry of Legendre functions about
        !   equator. For the first derivative in theta, these signs are
        !   reversed.
        !
        !----------------------------------------------------------------------
        do l = 0, lmax_comp, 1
            do m = 0, l, 1
                if (mod(l-m,2) == 0) then
                    fsymsign(l+1,m+1) = 1

                else
                    fsymsign(l+1,m+1) = -1

                end if

            end do

        end do

        !----------------------------------------------------------------------
        !
        !   Precompute square roots of integers that are used several times.
        !
        !----------------------------------------------------------------------
        do l=1, 2*lmax_comp+1
            sqr(l) = sqrt(dble(l))
        end do

        !----------------------------------------------------------------------
        !
        !   Precompute multiplicative factors used in recursion relationships
        !       P(l,m) = x*f1(l,m)*P(l-1,m) - P(l-2,m)*f2(l,m)
        !       k = l*(l+1)/2 + m + 1
        !   Note that prefactors are not used for the case when m=l as a
        !   different recursion is used. Furthermore, for m=l-1, Plmbar(l-2,m)
        !   is assumed to be zero.
        !
        !----------------------------------------------------------------------
        if (lmax_comp /= 0) then
            ff1(2,1) = sqr(3)
            ff2(2,1) = 0.0_dp

        end if

        do l = 2, lmax_comp, 1
            ff1(l+1,1) = sqr(2*l-1) * sqr(2*l+1) / dble(l)
            ff2(l+1,1) = dble(l-1) * sqr(2*l+1) / sqr(2*l-3) / dble(l)

            do m = 1, l-2, 1
                ff1(l+1,m+1) = sqr(2*l+1) * sqr(2*l-1) / sqr(l+m) / sqr(l-m)
                ff2(l+1,m+1) = sqr(2*l+1) * sqr(l-m-1) * sqr(l+m-1) &
                               / sqr(2*l-3) / sqr(l+m) / sqr(l-m)
            end do

            ff1(l+1,l) = sqr(2*l+1) * sqr(2*l-1) / sqr(l+m) / sqr(l-m)
            ff2(l+1,l) = 0.0_dp

        end do

        lmax_old = lmax_comp

    end if

    !--------------------------------------------------------------------------
    !
    !   Create generic plan for grid.
    !
    !--------------------------------------------------------------------------
    plan = fftw_plan_dft_c2r_1d(nlong, coef(1:nlong/2+1), grid(1:nlong), &
                                FFTW_MEASURE)

    !--------------------------------------------------------------------------
    !
    !   Determine Clms one l at a time by intergrating over latitude.
    !
    !--------------------------------------------------------------------------
    i_eq = n / 2 + 1  ! Index correspondong to zero latitude

    ! First do equator

    r_ex = a
    theta = pi / 2.0_dp
    z = 0.0_dp
    u = 1.0_dp
    lat = 0.0_dp

    coefr(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefr0 = 0.0_dp

    coefrr(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefrr0 = 0.0_dp

    coeft(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coeft0 = 0.0_dp

    coeftp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coeftp0 = 0.0_dp

    coefp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefp0 = 0.0_dp

    coefrt(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefrt0 = 0.0_dp

    coefrp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefrp0 = 0.0_dp

    coefpp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coefpp0 = 0.0_dp

    coeftt(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
    coeftt0 = 0.0_dp

    pm2 = 1.0_dp

    tempr = -cilm(1,1,1) * pm2 ! l = 0
    coefr0 = coefr0 + tempr

    tempr = 2 * cilm(1,1,1) * pm2  ! l = 0
    coefrr0 = coefrr0 + tempr

    ! derivative in theta and phi of l=0 term is 0. No need to calculate this

    if (lmax_comp /= 0) then    ! l = 1
        prefactor(1) = r0 / r_ex

        do l = 2, lmax_comp, 1
            prefactor(l) = prefactor(l-1) * r0 / r_ex
        end do

        pm1 = 0.0_dp

        dpl = ff1(2,1)
        tempr = cilm(1,2,1) * dpl * prefactor(1)
        coeft0 = coeft0 + tempr

        tempr = cilm(1,2,1) * dpl * (-2) * prefactor(1)  ! -2 = (l+1) prefactor
        coefrt0 = coefrt0 + tempr

    end if

    do l = 2, lmax_comp, 1
        l1 = l + 1
        p = - ff2(l1,1) * pm2
        tempr = cilm(1,l1,1) * p * (-l1) * prefactor(l)
        coefr0 = coefr0 + tempr

        tempr = cilm(1,l1,1) * p * (l1) * (l1+1) * prefactor(l)
        coefrr0 = coefrr0 + tempr

        dpl = l * ( sqr(2*l+1) / sqr(2*l-1) * pm1 )
        tempr = cilm(1,l1,1) * dpl * prefactor(l)
        coeft0 = coeft0 + tempr

        tempr = cilm(1,l1,1) * dpl * (-l1) * prefactor(l)
        coefrt0 = coefrt0 + tempr

        dpl2 = -l*l1 * p
        tempr = cilm(1,l1,1) * dpl2 * prefactor(l)
        coeftt0 = coeftt0 + tempr

        pm2 = pm1
        pm1 = p

    end do

    pmm = sqr(2) * scalef

    rescalem = 1.0_dp / scalef

    do m = 1, lmax_comp-1, 1

        m1 = m + 1

        pmm = pmm * sqr(2*m+1) / sqr(2*m)
        pm2 = pmm

        tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 * (-m-1) &
                * prefactor(m)    ! (m,m)
        coefr(m1) = coefr(m1) + tempc

        tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 * (m+1) &
                * (m+2) * prefactor(m) ! (m,m)
        coefrr(m1) = coefrr(m1) + tempc

        tempc = cmplx(cilm(2,m1,m1), cilm(1,m1,m1), dp) * pm2 &
                * prefactor(m) * m ! (m,m)
        coefp(m1) = coefp(m1) + tempc

        tempc = - cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 &
                * prefactor(m) * m**2 ! (m,m)
        coefpp(m1) = coefpp(m1) + tempc

        tempc = cmplx(cilm(2,m1,m1), cilm(1,m1,m1), dp) * pm2 *(-m-1) &
                * prefactor(m) * m ! (m,m)
        coefrp(m1) = coefrp(m1) + tempc

        dpl2 = -(m*m1 -(m**2)) * pm2
        tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * dpl2 &
                * prefactor(m) ! (m,m)
        coeftt(m1) = coeftt(m1) + tempc

        pm1 = 0.0_dp

        dpl = sqr(2*m+3) * pmm
        tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * dpl &
                * prefactor(m+1)    ! (m+1,m)
        coeft(m1) = coeft(m1) + tempc

        tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * dpl * (-m-2) &
                * prefactor(m+1)   ! (m+1,m)
        coefrt(m1) = coefrt(m1) + tempc

        tempc = cmplx(cilm(2,m1+1,m1), cilm(1,m1+1,m1), dp) * dpl &
                * prefactor(m+1) * m  ! (m+1,m)
        coeftp(m1) = coeftp(m1) + tempc

        do l = m + 2, lmax_comp, 1
            l1 = l + 1
            p = - ff2(l1,m1) * pm2
            tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p * (-l1) &
                    * prefactor(l)
            coefr(m1) = coefr(m1) + tempc

            tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p * (l1) &
                    * (l1+1) * prefactor(l)
            coefrr(m1) = coefrr(m1) + tempc

            tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * p &
                    * prefactor(l) * m
            coefp(m1) = coefp(m1) + tempc

            tempc = - cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p &
                    * prefactor(l) * m**2
            coefpp(m1) = coefpp(m1) + tempc

            tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * p * (-l1) &
                    * prefactor(l) * m
            coefrp(m1) = coefrp(m1) + tempc

            dpl = sqr(2*l+1) * sqr(l-m) * sqr(l+m) / sqr(2*l-1) * pm1
            tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * dpl &
                    * prefactor(l)
            coeft(m1) = coeft(m1) + tempc

            tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * dpl * (-l1) &
                    * prefactor(l)
            coefrt(m1) = coefrt(m1) + tempc

            tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * dpl * &
                    prefactor(l) * m
            coeftp(m1) = coeftp(m1) + tempc

            dpl2 = - (l * l1 -(m**2) / u**2) * p
            tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * dpl2 &
                    * prefactor(l)
            coeftt(m1) = coeftt(m1) + tempc

            pm2 = pm1
            pm1 = p

        end do

        coefr(m1) = coefr(m1) * rescalem
        coefrr(m1) = coefrr(m1) * rescalem
        coeft(m1) = coeft(m1) * rescalem
        coeftt(m1) = coeftt(m1) * rescalem
        coefrt(m1) = coefrt(m1) * rescalem
        coefp(m1) = coefp(m1) * rescalem
        coefpp(m1) = coefpp(m1) * rescalem
        coeftp(m1) = coeftp(m1) * rescalem
        coefrp(m1) = coefrp(m1) * rescalem

    end do

    if (lmax_comp /= 0) then

        pmm = pmm * sqr(2*lmax_comp+1) / sqr(2*lmax_comp) * rescalem
        tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                        - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                        * (-lmax_comp-1) * prefactor(lmax_comp)
        coefr(lmax_comp+1) = coefr(lmax_comp+1) + tempc

        tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                        - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                        * (lmax_comp+1) * (lmax_comp+2) * prefactor(lmax_comp)
        coefrr(lmax_comp+1) = coefrr(lmax_comp+1) + tempc

        tempc = cmplx(cilm(2,lmax_comp+1,lmax_comp+1), &
                        cilm(1,lmax_comp+1,lmax_comp+1), dp) * pmm &
                        * prefactor(lmax_comp) * lmax_comp
        coefp(lmax_comp+1) = coefp(lmax_comp+1) + tempc

        tempc = - cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                        - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                        * prefactor(lmax_comp) * lmax_comp**2
        coefpp(lmax_comp+1) = coefpp(lmax_comp+1) + tempc

        tempc = cmplx(cilm(2,lmax_comp+1,lmax_comp+1), &
                        cilm(1,lmax_comp+1,lmax_comp+1), dp) * pmm &
                        * (-lmax_comp-1) * prefactor(lmax_comp) * lmax_comp
        coefrp(lmax_comp+1) = coefrp(lmax_comp+1) + tempc

        dpl2 = -(lmax_comp*(lmax_comp+1)-(lmax_comp**2)) * pmm
        tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                        - cilm(2,lmax_comp+1,lmax_comp+1), dp) * dpl2 &
                        * prefactor(lmax_comp)
        coeftt(lmax_comp+1) = coeftt(lmax_comp+1) + tempc

    end if

    coefr0 = coefr0 * gm / r_ex**2
    coefr(2:lmax+1) = coefr(2:lmax+1) * gm / r_ex**2

    coefrt0 = - coefrt0 * gm / r_ex**2
    coefrt(2:lmax+1) = - coefrt(2:lmax+1) * gm / r_ex**2

    coefrr0 = coefrr0 * gm / r_ex**3
    coefrr(2:lmax+1) = coefrr(2:lmax+1) * gm / r_ex**3

    coeft0 = - coeft0 * gm / r_ex
    coeft(2:lmax+1) = - coeft(2:lmax+1) * gm / r_ex

    coeftt0 = coeftt0 * gm / r_ex
    coeftt(2:lmax+1) = coeftt(2:lmax+1) * gm / r_ex

    coeftp0 = - coeftp0 * gm / r_ex
    coeftp(2:lmax+1) = - coeftp(2:lmax+1) * gm / r_ex

    coefp0 = coefp0 * gm / r_ex
    coefp(2:lmax+1) = coefp(2:lmax+1) * gm / r_ex

    coefpp0 = coefpp0 * gm / r_ex
    coefpp(2:lmax+1) = coefpp(2:lmax+1) * gm / r_ex

    coefrp0 = coefrp0 * gm / r_ex**2
    coefrp(2:lmax+1) = coefrp(2:lmax+1) * gm / r_ex**2

    ! Vzz = Vrr
    coef(1) = cmplx(coefrr0, 0.0_dp, dp)
    coef(2:lmax+1) = coefrr(2:lmax+1) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vzz(i_eq,1:nlong) = grid(1:nlong)

    ! Vxx = 1/r Vr + 1/r^2 Vtt
    coef(1) = cmplx(coefr0/r_ex + coeftt0/r_ex**2, 0.0_dp, dp)
    coef(2:lmax+1) = (coefr(2:lmax+1)/r_ex + coeftt(2:lmax+1)/r_ex**2 ) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vxx(i_eq,1:nlong) = grid(1:nlong)

    ! Vyy = 1/r Vr + 1/r^2 /tan(t) Vt + 1/r^2 /sin(t)^2 Vpp
    coef(1) = cmplx(coefr0/r_ex + coefpp0/(r_ex**2), 0.0_dp, dp)
    coef(2:lmax+1) = (coefr(2:lmax+1)/r_ex &
                     + coefpp(2:lmax+1)/(r_ex**2) ) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vyy(i_eq,1:nlong) = grid(1:nlong)

    ! Vxy = 1/r^2 /sin(t) Vtp - cos(t)/sin(t)^2 /r^2 Vp
    coef(1) = cmplx(coeftp0/r_ex**2, 0.0_dp, dp)
    coef(2:lmax+1) = (coeftp(2:lmax+1)/r_ex**2 ) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vxy(i_eq,1:nlong) = grid(1:nlong)

    ! Vxz = 1/r^2 Vt - 1/r Vrt
    coef(1) = cmplx(coeft0/r_ex**2 - coefrt0/r_ex, 0.0_dp, dp)
    coef(2:lmax+1) = (coeft(2:lmax+1) / r_ex**2 &
                    - coefrt(2:lmax+1)/r_ex ) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vxz(i_eq,1:nlong) = grid(1:nlong)

    ! Vyz = 1/r^2 /sin(t) Vp - 1/r /sin(t) Vrp
    coef(1) = cmplx(coefp0/r_ex**2 - coefrp0/r_ex, 0.0_dp, dp)
    coef(2:lmax+1) = (coefp(2:lmax+1) / r_ex**2 &
                    - coefrp(2:lmax+1)/r_ex ) / 2.0_dp

    if (present(sampling)) then
        if (sampling == 2) then
            coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
        end if
    end if

    call fftw_execute_dft_c2r(plan, coef, grid)
    vyz(i_eq,1:nlong) = grid(1:nlong)

    do i = 1, i_eq - 1, 1

        i_s = 2 * i_eq - i

        theta = pi * dble(i-1) / dble(n)
        z = cos(theta)
        u = sqrt( (1.0_dp-z) * (1.0_dp+z) )
        sint = sin(theta)
        lat = pi / 2.0_dp - theta

        if (i == 1) then      ! Reference ellipsoid radius
            r_ex = a * (1.0_dp - f)

        else
            r_ex = (1.0_dp + tan(lat)**2) / &
                   (1.0_dp  + tan(lat)**2 / (1.0_dp - f)**2 )
            r_ex = a * sqrt(r_ex)

        end if

        coefr(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefr0 = 0.0_dp
        coefrs(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrs0 = 0.0_dp

        coefrr(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrr0 = 0.0_dp
        coefrrs(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrrs0 = 0.0_dp

        coeft(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coeft0 = 0.0_dp
        coefts(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefts0 = 0.0_dp

        coeftp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coeftp0 = 0.0_dp
        coeftps(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coeftps0 = 0.0_dp

        coefp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefp0 = 0.0_dp
        coefps(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefps0 = 0.0_dp

        coefrt(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrt0 = 0.0_dp
        coefrts(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrts0 = 0.0_dp

        coefrp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrp0 = 0.0_dp
        coefrps(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefrps0 = 0.0_dp

        coefpp(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefpp0 = 0.0_dp
        coefpps(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coefpps0 = 0.0_dp

        coeftt(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coeftt0 = 0.0_dp
        coeftts(1:lmax+2) = cmplx(0.0_dp, 0.0_dp, dp)
        coeftts0 = 0.0_dp

        pm2 = 1.0_dp

        tempr = -cilm(1,1,1) * pm2 ! l = 0
        coefr0 = coefr0 + tempr
        coefrs0 = coefrs0 + tempr   ! fsymsign is always 1 for l=m=0

        tempr = 2 * cilm(1,1,1) * pm2  ! l = 0
        coefrr0 = coefrr0 + tempr
        coefrrs0 = coefrrs0 + tempr

        ! derivative in theta and phi of l=0 term is 0, so no need to
        ! calculate this

        if (lmax_comp /= 0) then    ! l = 1
            prefactor(1) = r0 / r_ex

            do l = 2, lmax_comp, 1
                prefactor(l) = prefactor(l-1) * r0 / r_ex
            end do

            pm1 = ff1(2,1) * z

            ! -2 = (l+1) prefactor
            tempr = cilm(1,2,1) * pm1 * (-2) * prefactor(1)
            coefr0 = coefr0 + tempr
            coefrs0 = coefrs0 - tempr   ! fsymsign = -1

            ! 6 = (l+1)*(l+2) prefactor
            tempr = cilm(1,2,1) * pm1 * (6) * prefactor(1)
            coefrr0 = coefrr0 + tempr
            coefrrs0 = coefrrs0 - tempr     ! fsymsign = -1

            ! dpl is the first derivative with respect to Z
            dpl = ff1(2,1)
            tempr = cilm(1,2,1) * dpl * prefactor(1)
            coeft0 = coeft0 + tempr
            coefts0 = coefts0 + tempr   ! reverse fsymsign

            ! -2 = (l+1) prefactor
            tempr = cilm(1,2,1) * dpl * (-2) * prefactor(1)
            coefrt0 = coefrt0 + tempr
            coefrts0 = coefrts0 + tempr     ! reverse fsymsign

            ! dpl2 is the second derivative with respect to THETA. Must
            ! multiply dpl by -sin(theta) in the recurrence relation
            ! note that the first term is symmetric about the equator according
            ! the fsymsign, whereas the second term differs
            ! by -1.
            dpl2 = -2 * pm1 + z * dpl
            dpl2s = 2 * pm1 + z * dpl
            tempr = cilm(1,2,1) * prefactor(1)
            coeftt0 = coeftt0 + tempr * dpl2
            coeftts0 = coeftts0 + tempr * dpl2s

        end if

        do l = 2, lmax_comp, 1
            l1 = l + 1
            p = ff1(l1,1) * z * pm1 - ff2(l1,1) * pm2
            tempr = cilm(1,l1,1) * p * (-l1) * prefactor(l)
            coefr0 = coefr0 + tempr
            coefrs0 = coefrs0 + tempr * fsymsign(l1,1)

            tempr = cilm(1,l1,1) * p * (l1) * (l1+1) * prefactor(l)
            coefrr0 = coefrr0 + tempr
            coefrrs0 = coefrrs0 + tempr * fsymsign(l1,1)

            dpl = l * ( sqr(2*l+1) / sqr(2*l-1) * pm1 - z * p ) / u**2
            tempr = cilm(1,l1,1) * dpl * prefactor(l)
            coeft0 = coeft0 + tempr
            coefts0 = coefts0 - tempr * fsymsign(l1,1)  ! reverse fsymsign

            tempr = cilm(1,l1,1) * dpl * (-l1) * prefactor(l)
            coefrt0 = coefrt0 + tempr
            coefrts0 = coefrts0 - tempr * fsymsign(l1,1)

            dpl2 = -l * l1 * p + z * dpl
            dpl2s = -l * l1 * p * fsymsign(l1,1) - z * dpl * fsymsign(l1,1)
            tempr = cilm(1,l1,1) * prefactor(l)
            coeftt0 = coeftt0 + tempr * dpl2
            coeftts0 = coeftts0 + tempr * dpl2s

            pm2 = pm1
            pm1 = p

        end do

        pmm = sqr(2) * scalef

        rescalem = 1.0_dp / scalef

        do m = 1, lmax_comp-1, 1

            m1 = m + 1
            rescalem = rescalem * u

            pmm = pmm * sqr(2*m+1) / sqr(2*m)
            pm2 = pmm

            tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 * (-m-1) &
                    * prefactor(m)    ! (m,m)
            coefr(m1) = coefr(m1) + tempc
            coefrs(m1) = coefrs(m1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 * (m+1) &
                    * (m+2) * prefactor(m) ! (m,m)
            coefrr(m1) = coefrr(m1) + tempc
            coefrrs(m1) = coefrrs(m1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(2,m1,m1), cilm(1,m1,m1), dp) * pm2 &
                    * prefactor(m) * m ! (m,m)
            coefp(m1) = coefp(m1) + tempc
            coefps(m1) = coefps(m1) + tempc ! fsymsign = 1

            tempc = - cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * pm2 &
                    * prefactor(m) * m**2 ! (m,m)
            coefpp(m1) = coefpp(m1) + tempc
            coefpps(m1) = coefpps(m1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(2,m1,m1), cilm(1,m1,m1), dp) * pm2 *(-m-1) &
                    * prefactor(m) * m ! (m,m)
            coefrp(m1) = coefrp(m1) + tempc
            coefrps(m1) = coefrps(m1) + tempc ! fsymsign = 1

            dpl = -m * z * pm2 / u**2
            tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * dpl &
                    * prefactor(m)  ! (m,m)
            coeft(m1) = coeft(m1) + tempc
            coefts(m1) = coefts(m1) - tempc ! reverse fsymsign

            tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) * dpl * (-m-1) &
                    * prefactor(m) ! (m,m)
            coefrt(m1) = coefrt(m1) + tempc
            coefrts(m1) = coefrts(m1) - tempc ! reverse fsymsign

            tempc = cmplx(cilm(2,m1,m1), cilm(1,m1,m1), dp) * dpl &
                    * prefactor(m) * m ! (m,m)
            coeftp(m1) = coeftp(m1) + tempc
            coeftps(m1) = coeftps(m1) - tempc ! reverse fsymsign

            dpl2 = -(m*m1 - (m**2)/u**2) * pm2 + z * dpl
            dpl2s = -(m*m1 - (m**2)/u**2) * pm2 - z * dpl
            tempc = cmplx(cilm(1,m1,m1), - cilm(2,m1,m1), dp) &
                    * prefactor(m)   ! (m,m)
            coeftt(m1) = coeftt(m1) + tempc * dpl2
            coeftts(m1) = coeftts(m1) + tempc * dpl2s

            pm1 = z * ff1(m1+1,m1) * pm2
            tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * pm1 &
                    * (-m-2) * prefactor(m+1)  ! (m+1,m)
            coefr(m1) = coefr(m1) + tempc
            coefrs(m1) = coefrs(m1) - tempc ! fsymsign = -1

            tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * pm1 &
                    * (m+2) * (m+3) * prefactor(m+1)   ! (m+1,m)
            coefrr(m1) = coefrr(m1) + tempc
            coefrrs(m1) = coefrrs(m1) - tempc ! fsymsign = -1

            tempc = cmplx(cilm(2,m1+1,m1), cilm(1,m1+1,m1), dp) * pm1 &
                    * prefactor(m+1) * m ! (m+1,m)
            coefp(m1) = coefp(m1) + tempc
            coefps(m1) = coefps(m1) - tempc ! fsymsign = -1

            tempc = - cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * pm1 &
                    * prefactor(m+1) * m**2   ! (m+1,m)
            coefpp(m1) = coefpp(m1) + tempc
            coefpps(m1) = coefpps(m1) - tempc ! fsymsign = -1

            tempc = cmplx(cilm(2,m1+1,m1), cilm(1,m1+1,m1), dp) * pm1 &
                    * (-m-2) * prefactor(m+1) * m    ! (m+1,m)
            coefrp(m1) = coefrp(m1) + tempc
            coefrps(m1) = coefrps(m1) - tempc ! fsymsign = -1

            dpl = (sqr(2*m+3) * pmm - z * (m+1) * pm1) / u**2
            tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * dpl &
                    * prefactor(m+1)    ! (m+1,m)
            coeft(m1) = coeft(m1) + tempc
            coefts(m1) = coefts(m1) + tempc ! reverse fsymsign

            tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) * dpl &
                    * (-m-2) * prefactor(m+1)   ! (m+1,m)
            coefrt(m1) = coefrt(m1) + tempc
            coefrts(m1) = coefrts(m1) + tempc   ! reverse fsymsign

            tempc = cmplx(cilm(2,m1+1,m1), cilm(1,m1+1,m1), dp) * dpl &
                    * prefactor(m+1) * m  ! (m+1,m)
            coeftp(m1) = coeftp(m1) + tempc
            coeftps(m1) = coeftps(m1) + tempc ! reverse fsymsign

            dpl2 = -(m1*(m1+1) - (m**2)/u**2) * pm1 + z * dpl
            dpl2s = (m1*(m1+1) - (m**2)/u**2) * pm1 + z * dpl
            tempc = cmplx(cilm(1,m1+1,m1), - cilm(2,m1+1,m1), dp) &
                    * prefactor(m+1) ! (m+1,m)
            coeftt(m1) = coeftt(m1) + tempc  * dpl2
            coeftts(m1) = coeftts(m1) + tempc * dpl2s

            do l=m+2, lmax_comp, 1
                l1 = l + 1
                p = z * ff1(l1,m1) * pm1 - ff2(l1,m1) * pm2
                tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p * (-l1) &
                        * prefactor(l)
                coefr(m1) = coefr(m1) + tempc
                coefrs(m1) = coefrs(m1) + tempc * fsymsign(l1,m1)

                tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p * (l1) &
                        * (l1+1) * prefactor(l)
                coefrr(m1) = coefrr(m1) + tempc
                coefrrs(m1) = coefrrs(m1) + tempc * fsymsign(l1,m1)

                tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * p &
                        * prefactor(l) * m
                coefp(m1) = coefp(m1) + tempc
                coefps(m1) = coefps(m1) + tempc * fsymsign(l1,m1)

                tempc = - cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * p &
                        * prefactor(l) * m**2
                coefpp(m1) = coefpp(m1) + tempc
                coefpps(m1) = coefpps(m1) + tempc * fsymsign(l1,m1)

                tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * p * (-l1) &
                        * prefactor(l) * m
                coefrp(m1) = coefrp(m1) + tempc
                coefrps(m1) = coefrps(m1) + tempc * fsymsign(l1,m1)

                dpl = ( sqr(2*l+1) * sqr(l-m) * sqr(l+m) / sqr(2*l-1) &
                       * pm1 - z * l * p) / u**2
                tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * dpl &
                        * prefactor(l)
                coeft(m1) = coeft(m1) + tempc
                ! reverse fsymsign
                coefts(m1) = coefts(m1) - tempc * fsymsign(l1,m1)

                tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) * dpl &
                        * (-l1) * prefactor(l)
                coefrt(m1) = coefrt(m1) + tempc
                ! reverse fsymsign
                coefrts(m1) = coefrts(m1) - tempc * fsymsign(l1,m1)

                tempc = cmplx(cilm(2,l1,m1), cilm(1,l1,m1), dp) * dpl &
                        * prefactor(l) * m
                coeftp(m1) = coeftp(m1) + tempc
                ! reverse fsymsign
                coeftps(m1) = coeftps(m1) - tempc * fsymsign(l1,m1)

                dpl2 = -(l * l1 -(m**2)/u**2) * p + z * dpl
                dpl2s = -(l * l1 -(m**2)/u**2) * p * fsymsign(l1,m1) &
                       - z * dpl * fsymsign(l1,m1)
                tempc = cmplx(cilm(1,l1,m1), - cilm(2,l1,m1), dp) &
                        * prefactor(l)
                coeftt(m1) = coeftt(m1) + tempc * dpl2
                coeftts(m1) = coeftts(m1) +  tempc * dpl2s

                pm2 = pm1
                pm1 = p

            end do

            coefr(m1) = coefr(m1) * rescalem
            coefrs(m1) = coefrs(m1) * rescalem

            coefrr(m1) = coefrr(m1) * rescalem
            coefrrs(m1) = coefrrs(m1) * rescalem

            coeft(m1) = coeft(m1) * rescalem
            coefts(m1) = coefts(m1) * rescalem

            coeftt(m1) = coeftt(m1) * rescalem
            coeftts(m1) = coeftts(m1) * rescalem

            coefrt(m1) = coefrt(m1) * rescalem
            coefrts(m1) = coefrts(m1) * rescalem

            coefp(m1) = coefp(m1) * rescalem
            coefps(m1) = coefps(m1) * rescalem

            coefpp(m1) = coefpp(m1) * rescalem
            coefpps(m1) = coefpps(m1) * rescalem

            coeftp(m1) = coeftp(m1) * rescalem
            coeftps(m1) = coeftps(m1) * rescalem

            coefrp(m1) = coefrp(m1) * rescalem
            coefrps(m1) = coefrps(m1) * rescalem

        end do

        if (lmax_comp /= 0) then

            rescalem = rescalem * u

            pmm = pmm * sqr(2*lmax_comp+1) / sqr(2*lmax_comp) * rescalem
            tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                            * (-lmax_comp-1) * prefactor(lmax_comp)
            coefr(lmax_comp+1) = coefr(lmax_comp+1) + tempc
            coefrs(lmax_comp+1) = coefrs(lmax_comp+1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                            * (lmax_comp+1) * (lmax_comp+2) &
                            * prefactor(lmax_comp)
            coefrr(lmax_comp+1) = coefrr(lmax_comp+1) + tempc
            coefrrs(lmax_comp+1) = coefrrs(lmax_comp+1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(2,lmax_comp+1,lmax_comp+1), &
                            cilm(1,lmax_comp+1,lmax_comp+1), dp) * pmm &
                            * prefactor(lmax_comp) * lmax_comp
            coefp(lmax_comp+1) = coefp(lmax_comp+1) + tempc
            coefps(lmax_comp+1) = coefps(lmax_comp+1) + tempc   ! fsymsign = 1

            tempc = - cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) * pmm &
                            * prefactor(lmax_comp) * lmax_comp**2
            coefpp(lmax_comp+1) = coefpp(lmax_comp+1) + tempc
            coefpps(lmax_comp+1) = coefpps(lmax_comp+1) + tempc ! fsymsign = 1

            tempc = cmplx(cilm(2,lmax_comp+1,lmax_comp+1), &
                            cilm(1,lmax_comp+1,lmax_comp+1), dp) * pmm &
                            * (-lmax_comp-1) * prefactor(lmax_comp) * lmax_comp
            coefrp(lmax_comp+1) = coefrp(lmax_comp+1) + tempc
            coefrps(lmax_comp+1) = coefrps(lmax_comp+1) + tempc ! fsymsign = 1

            dpl = -lmax_comp * z * pmm / u**2
            tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) * dpl &
                            * prefactor(lmax_comp)
            coeft(lmax_comp+1) = coeft(lmax_comp+1) + tempc
            ! reverse fsymsign
            coefts(lmax_comp+1) = coefts(lmax_comp+1) - tempc

            tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) * dpl &
                            * (-lmax_comp-1) * prefactor(lmax_comp)
            coefrt(lmax_comp+1) = coefrt(lmax_comp+1) + tempc
            ! reverse fsymsign
            coefrts(lmax_comp+1) = coefrts(lmax_comp+1) - tempc

            tempc = cmplx(cilm(2,lmax_comp+1,lmax_comp+1), &
                            cilm(1,lmax_comp+1,lmax_comp+1), dp) * dpl &
                            * prefactor(lmax_comp) * lmax_comp
            coeftp(lmax_comp+1) = coeftp(lmax_comp+1) + tempc
            coeftps(lmax_comp+1) = coeftps(lmax_comp+1) - tempc

            dpl2 = -(lmax_comp*(lmax_comp+1)-(lmax_comp**2)/u**2) * pmm &
                   + z * dpl
            dpl2s = -(lmax_comp*(lmax_comp+1)-(lmax_comp**2)/u**2) * pmm &
                    - z * dpl
            tempc = cmplx(cilm(1,lmax_comp+1,lmax_comp+1), &
                            - cilm(2,lmax_comp+1,lmax_comp+1), dp) &
                            * prefactor(lmax_comp)
            coeftt(lmax_comp+1) = coeftt(lmax_comp+1) + tempc * dpl2
            coeftts(lmax_comp+1) = coeftts(lmax_comp+1) + tempc * dpl2s

        end if

        ! Note that the first angular derivatives are with repsect to z,
        ! but that the second is with respect to theta.

        coefr0 = coefr0 * gm / r_ex**2
        coefr(2:lmax+1) = coefr(2:lmax+1) * gm / r_ex**2

        coefrs0 = coefrs0 * gm / r_ex**2
        coefrs(2:lmax+1) = coefrs(2:lmax+1) * gm / r_ex**2

        coefrt0 = -sint * coefrt0 * gm / r_ex**2
        coefrt(2:lmax+1) = -sint * coefrt(2:lmax+1) * gm / r_ex**2

        coefrts0 = -sint * coefrts0 * gm / r_ex**2
        coefrts(2:lmax+1) = -sint * coefrts(2:lmax+1) * gm / r_ex**2

        coefrr0 = coefrr0 * gm / r_ex**3
        coefrr(2:lmax+1) = coefrr(2:lmax+1) * gm / r_ex**3

        coefrrs0 = coefrrs0 * gm / r_ex**3
        coefrrs(2:lmax+1) = coefrrs(2:lmax+1) * gm / r_ex**3

        coeft0 = -sint * coeft0 * gm / r_ex
        coeft(2:lmax+1) = -sint * coeft(2:lmax+1) * gm / r_ex

        coefts0 = -sint * coefts0 * gm / r_ex
        coefts(2:lmax+1) = -sint * coefts(2:lmax+1) * gm / r_ex

        coeftt0 = coeftt0 * gm / r_ex
        coeftt(2:lmax+1) = coeftt(2:lmax+1) * gm / r_ex

        coeftts0 = coeftts0 * gm / r_ex
        coeftts(2:lmax+1) = coeftts(2:lmax+1) * gm / r_ex

        coeftp0 = -sint * coeftp0 * gm / r_ex
        coeftp(2:lmax+1) = -sint * coeftp(2:lmax+1) * gm / r_ex

        coeftps0 = -sint * coeftps0 * gm / r_ex
        coeftps(2:lmax+1) = -sint * coeftps(2:lmax+1) * gm / r_ex

        coefp0 = coefp0 * gm / r_ex
        coefp(2:lmax+1) = coefp(2:lmax+1) * gm / r_ex

        coefps0 = coefps0 * gm / r_ex
        coefps(2:lmax+1) = coefps(2:lmax+1) * gm / r_ex

        coefpp0 = coefpp0 * gm / r_ex
        coefpp(2:lmax+1) = coefpp(2:lmax+1) * gm / r_ex

        coefpps0 = coefpps0 * gm / r_ex
        coefpps(2:lmax+1) = coefpps(2:lmax+1) * gm / r_ex

        coefrp0 = coefrp0 * gm / r_ex**2
        coefrp(2:lmax+1) = coefrp(2:lmax+1) * gm / r_ex**2

        coefrps0 = coefrps0 * gm / r_ex**2
        coefrps(2:lmax+1) = coefrps(2:lmax+1) * gm / r_ex**2

        ! Vzz = Vrr
        coef(1) = cmplx(coefrr0, 0.0_dp, dp)
        coef(2:lmax+1) = coefrr(2:lmax+1) / 2.0_dp

        if (present(sampling)) then
            if (sampling == 2) then
                coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
            end if

        end if

        call fftw_execute_dft_c2r(plan, coef, grid)
        vzz(i,1:nlong) = grid(1:nlong)

        if (i == 1) then
            ! These derivatives are undefined at the pole
            vxx(1,1:nlong) = 0.0_dp
            vyy(1,1:nlong) = 0.0_dp
            vxy(1,1:nlong) = 0.0_dp
            vxz(1,1:nlong) = 0.0_dp
            vyz(1,1:nlong) = 0.0_dp

        else
            ! Vxx = 1/r Vr + 1/r^2 Vtt
            coef(1) = cmplx(coefr0/r_ex + coeftt0/r_ex**2, 0.0_dp, dp)
            coef(2:lmax+1) = (coefr(2:lmax+1)/r_ex &
                             + coeftt(2:lmax+1)/r_ex**2 ) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vxx(i,1:nlong) = grid(1:nlong)

            ! Vyy = 1/r Vr + 1/r^2 /tan(t) Vt + 1/r^2 /sin(t)^2 Vpp
            coef(1) = cmplx(coefr0/r_ex + coeft0/(r_ex**2)/tan(theta) &
                            + coefpp0/(r_ex**2)/u**2, 0.0_dp, dp)
            coef(2:lmax+1) = (coefr(2:lmax+1)/r_ex &
                            + coeft(2:lmax+1)/(r_ex**2)/tan(theta) + &
                            coefpp(2:lmax+1)/(r_ex**2)/u**2 ) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vyy(i,1:nlong) = grid(1:nlong)

            ! Vxy = 1/r^2 /sin(t) Vtp - cos(t)/sin(t)^2 /r^2 Vp
            coef(1) = cmplx(coeftp0/sint/r_ex**2 &
                            - coefp0/(r_ex**2)*z/u**2, 0.0_dp, dp)
            coef(2:lmax+1) = (coeftp(2:lmax+1)/sint/r_ex**2 &
                            - coefp(2:lmax+1)/(r_ex**2)*z/u**2 ) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vxy(i,1:nlong) = grid(1:nlong)

            ! Vxz = 1/r^2 Vt - 1/r Vrt
            coef(1) = cmplx(coeft0/r_ex**2 - coefrt0/r_ex, 0.0_dp, dp)
            coef(2:lmax+1) = (coeft(2:lmax+1)/r_ex**2 &
                            - coefrt(2:lmax+1)/r_ex ) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vxz(i,1:nlong) = grid(1:nlong)

            ! Vyz = 1/r^2 /sin(t) Vp - 1/r /sin(t) Vrp
            coef(1) = cmplx(coefp0/sint/r_ex**2 - coefrp0/sint/r_ex, 0.0_dp, dp)
            coef(2:lmax+1) = (coefp(2:lmax+1)/sint/r_ex**2 &
                            - coefrp(2:lmax+1)/sint/r_ex ) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vyz(i,1:nlong) = grid(1:nlong)

        end if

        ! don't compute value for south pole when extend = 0.
        if (.not. (i == 1 .and. extend_grid == 0) ) then
            ! Vzz = Vrr
            coef(1) = cmplx(coefrrs0, 0.0_dp, dp)
            coef(2:lmax+1) = coefrrs(2:lmax+1) / 2.0_dp

            if (present(sampling)) then
                if (sampling == 2) then
                    coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                end if
            end if

            call fftw_execute_dft_c2r(plan, coef, grid)
            vzz(i_s,1:nlong) = grid(1:nlong)

            if (i == 1) then
                ! These derivatives are undefined at the pole
                vxx(1,1:nlong) = 0.0_dp
                vyy(1,1:nlong) = 0.0_dp
                vxy(1,1:nlong) = 0.0_dp
                vxz(1,1:nlong) = 0.0_dp
                vyz(1,1:nlong) = 0.0_dp

            else
                ! Vxx = 1/r Vr + 1/r^2 Vtt
                coef(1) = cmplx(coefrs0/r_ex + coeftts0/r_ex**2, 0.0_dp, dp)
                coef(2:lmax+1) = (coefrs(2:lmax+1)/r_ex &
                                  + coeftts(2:lmax+1)/r_ex**2 ) / 2.0_dp

                if (present(sampling)) then
                    if (sampling == 2) then
                        coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                    end if
                end if

                call fftw_execute_dft_c2r(plan, coef, grid)
                vxx(i_s,1:nlong) = grid(1:nlong)

                ! Vyy = 1/r Vr + 1/r^2 /tan(t) Vt + 1/r^2 /sin(t)^2 Vpp
                coef(1) = cmplx(coefrs0/r_ex + coefts0/(r_ex**2)/tan(theta) &
                                + coefpps0/(r_ex**2)/u**2, 0.0_dp, dp)
                coef(2:lmax+1) = (coefrs(2:lmax+1)/r_ex &
                                  + coefts(2:lmax+1)/(r_ex**2)/tan(theta) + &
                                  coefpps(2:lmax+1)/(r_ex**2)/u**2 ) / 2.0_dp

                if (present(sampling)) then
                    if (sampling == 2) then
                        coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                    end if
                end if

                call fftw_execute_dft_c2r(plan, coef, grid)
                vyy(i_s,1:nlong) = grid(1:nlong)

                ! Vxy = 1/r^2 /sin(t) Vtp - cos(t)/sin(t)^2 /r^2 Vp
                coef(1) = cmplx(coeftps0/sint/r_ex**2 &
                                - coefps0/(r_ex**2)*z/u**2, 0.0_dp, dp)
                coef(2:lmax+1) = (coeftps(2:lmax+1)/sint/r_ex**2 &
                                  - coefps(2:lmax+1)/(r_ex**2)*z/u**2 ) &
                                  / 2.0_dp

                if (present(sampling)) then
                    if (sampling == 2) then
                        coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                    end if
                end if

                call fftw_execute_dft_c2r(plan, coef, grid)
                vxy(i_s,1:nlong) = grid(1:nlong)

                ! Vxz = 1/r^2 Vt - 1/r Vrt
                coef(1) = cmplx(coefts0/r_ex**2 - coefrts0/r_ex, 0.0_dp, dp)
                coef(2:lmax+1) = (coefts(2:lmax+1)/r_ex**2 &
                                  - coefrts(2:lmax+1)/r_ex ) / 2.0_dp

                if (present(sampling)) then
                    if (sampling == 2) then
                        coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                    end if
                end if

                call fftw_execute_dft_c2r(plan, coef, grid)
                vxz(i_s,1:nlong) = grid(1:nlong)

                ! Vyz = 1/r^2 /sin(t) Vp - 1/r /sin(t) Vrp
                coef(1) = cmplx(coefps0/sint/r_ex**2 - coefrps0/sint/r_ex, &
                                0.0_dp, dp)
                coef(2:lmax+1) = (coefps(2:lmax+1)/sint/r_ex**2 &
                                  - coefrps(2:lmax+1)/sint/r_ex ) / 2.0_dp

                if (present(sampling)) then
                    if (sampling == 2) then
                        coef(lmax+2:2*lmax+3) = cmplx(0.0_dp, 0.0_dp, dp)
                    end if
                end if

                call fftw_execute_dft_c2r(plan, coef, grid)
                vyz(i_s,1:nlong) = grid(1:nlong)

            end if

        end if

    end do

    if (extend_grid == 1) then
        vxx(1:nlat_out, nlong_out) = vxx(1:nlat_out, 1)
        vyy(1:nlat_out, nlong_out) = vyy(1:nlat_out, 1)
        vzz(1:nlat_out, nlong_out) = vzz(1:nlat_out, 1)
        vxy(1:nlat_out, nlong_out) = vxy(1:nlat_out, 1)
        vxz(1:nlat_out, nlong_out) = vxz(1:nlat_out, 1)
        vyz(1:nlat_out, nlong_out) = vyz(1:nlat_out, 1)
    end if

    call fftw_destroy_plan(plan)

end subroutine MakeGravGradGridDH

max_stars_repo_name: megies/SHTOOLS

max_issues_repo_name: megies/SHTOOLS

max_forks_repo_name: megies/SHTOOLS


text:
      DOUBLE PRECISION   FUNCTION PDLATRA( N, A, IA, JA, DESCA )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      INTEGER            IA, JA, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      DOUBLE PRECISION   A( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLATRA computes the trace of an N-by-N distributed matrix sub( A )
*  denoting A( IA:IA+N-1, JA:JA+N-1 ). The result is left on every
*  process of the grid.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on i.e the
*          order of the distributed submatrix sub( A ).  N >= 0.
*
*  A       (local input) DOUBLE PRECISION pointer into the local memory
*          to an array of dimension ( LLD_A, LOCc(JA+N-1) ). This array
*          contains the local pieces of the distributed matrix the trace
*          is to be computed.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  ====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            ICURCOL, ICURROW, II, IOFFA, J, JB, JJ, JN,
     $                   LDA, LL, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   TRACE
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DGSUM2D, INFOG2L
*     ..
*     .. External Functions ..
      INTEGER            ICEIL
      EXTERNAL           ICEIL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
*
      TRACE = ZERO
      IF( N.EQ.0 ) THEN
         PDLATRA = TRACE
         RETURN
      END IF
*
      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II, JJ,
     $              ICURROW, ICURCOL )
*
      JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
      JB = JN-JA+1
      LDA = DESCA( LLD_ )
      IOFFA = II + ( JJ - 1 ) * LDA
*
*     Handle first diagonal block separately
*
      IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
         DO 10 LL = IOFFA, IOFFA + (JB-1)*(LDA+1), LDA+1
            TRACE = TRACE + A( LL )
   10    CONTINUE
      END IF
      IF( MYROW.EQ.ICURROW )
     $   IOFFA = IOFFA + JB
      IF( MYCOL.EQ.ICURCOL )
     $   IOFFA = IOFFA + JB*LDA
      ICURROW = MOD( ICURROW+1, NPROW )
      ICURCOL = MOD( ICURCOL+1, NPCOL )
*
*     Loop over the remaining block of columns
*
      DO 30 J = JN+1, JA+N-1, DESCA( NB_ )
         JB = MIN( JA+N-J, DESCA( NB_ ) )
*
         IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
            DO 20 LL = IOFFA, IOFFA + (JB-1)*(LDA+1), LDA+1
               TRACE = TRACE + A( LL )
   20       CONTINUE
         END IF
         IF( MYROW.EQ.ICURROW )
     $      IOFFA = IOFFA + JB
         IF( MYCOL.EQ.ICURCOL )
     $      IOFFA = IOFFA + JB*LDA
         ICURROW = MOD( ICURROW+1, NPROW )
         ICURCOL = MOD( ICURCOL+1, NPCOL )
   30 CONTINUE
*
      CALL DGSUM2D( DESCA( CTXT_ ), 'All', ' ', 1, 1, TRACE, 1, -1,
     $              MYCOL )
*
      PDLATRA = TRACE
*
      RETURN
*
*     End of PDLATRA
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
MODULE MODULE_BOX

USE MODULE_PBC
USE MODULE_COM

CONTAINS

  !! This function is not used by gro, xtc and trr

  !! box(1,1) = v1_x, box(1,2) = v1_y, box(1,3) = v1_z
  !! box(2,1) = v2_x, box(2,2) = v2_y, box(2,3) = v2_z
  !! box(3,1) = v3_x, box(3,2) = v3_y, box(3,3) = v3_z

  !! Para gromacs:
  !!
  !! box(1,1) = v1_x, box(1,2) = 0,    box(1,3) = 0
  !! box(2,1) = v2_x, box(2,2) = v2_y, box(2,3) = 0
  !! box(3,1) = v3_x, box(3,2) = v3_y, box(3,3) = v3_z
  !!
  !! cell(1,1)= v1_x, cell(1,2)= alpha cell(1,3)= beta
  !! cell(2,1)= 0,    cell(2,2)= v2_y, cell(2,3)= gamma
  !! cell(3,1)= 0,    cell(3,2)= 0,    cell(3,3)= v3_z

  SUBROUTINE BOX2INVBOX (box,invbox,nframes)

    IMPLICIT NONE
    INTEGER,INTENT(IN)::nframes
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(IN)::box
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(OUT)::invbox

    invbox=0.0d0

    invbox(:,1,1)=1.0d0/box(:,1,1)
    invbox(:,2,2)=1.0d0/box(:,2,2)
    invbox(:,3,3)=1.0d0/box(:,3,3)

    invbox(:,2,1)=-box(:,2,1)/(box(:,1,1)*box(:,2,2))
    invbox(:,3,1)=(box(:,2,1)*box(:,3,2)-box(:,3,1)*box(:,2,2))/(box(:,1,1)*box(:,2,2)*box(:,3,3))
    invbox(:,3,2)=-box(:,3,2)/(box(:,2,2)*box(:,3,3))

  END SUBROUTINE BOX2INVBOX

  SUBROUTINE CELL2BOX (cell,box,volume,ortho,nframes)

    IMPLICIT NONE
    INTEGER,INTENT(IN)::nframes
    INTEGER::ii
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(INOUT)::cell
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(OUT)::box
    DOUBLE PRECISION,DIMENSION(nframes),INTENT(OUT)::volume
    INTEGER,INTENT(OUT)::ortho
    DOUBLE PRECISION::alpha,beta,gamm,x,y,z
    DOUBLE PRECISION,PARAMETER::fact_pi_div_180=1.745329251994329547437168059786927E-0002

    DO ii=1,nframes

      ortho=0
      box=0.0d0
      alpha=cell(ii,1,2)
      beta=cell(ii,1,3)
      gamm=cell(ii,2,3)
      x=cell(ii,1,1)
      y=cell(ii,2,2)
      z=cell(ii,3,3)

      box(ii,1,1)=x
      IF ((alpha==90.0d0).and.(beta==90.0d0).and.(gamm==90.0d0)) THEN
         box(ii,2,2)=y
         box(ii,3,3)=z
         volume(ii)=x*y*z
         ortho=1
      ELSE IF ((alpha==0.0d0).and.(beta==0.0d0).and.(gamm==0.0d0)) THEN
         cell(ii,1,2)=90.0d0
         cell(ii,1,3)=90.0d0
         cell(ii,2,3)=90.0d0
         box(ii,2,2)=y
         box(ii,3,3)=z
         volume(ii)=x*y*z
         ortho=1
      ELSE
         alpha=alpha/fact_pi_div_180
         beta=beta/fact_pi_div_180
         gamm=gamm/fact_pi_div_180
         box(ii,2,1)=y*cos(gamm)
         box(ii,2,2)=y*sin(gamm)  ! sqrt(y**2-box(2,2)**2)
         box(ii,3,1)=z*cos(beta)
         box(ii,3,2)=z*(cos(alpha)-cos(beta)*cos(gamm))/sin(gamm)
         box(ii,3,3)=sqrt(z*z-box(ii,3,1)**2-box(ii,3,2)**2)
         volume(ii)=box(ii,1,1)*box(ii,2,2)*box(ii,3,3)
      END IF

    END DO
  END SUBROUTINE CELL2BOX

  !! This function is used by gro, xtc and trr
  SUBROUTINE BOX2CELL (box,cell,volume,ortho,nframes)

    IMPLICIT NONE
    INTEGER,INTENT(IN)::nframes
    INTEGER::ii
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(OUT)::cell
    DOUBLE PRECISION,DIMENSION(nframes,3,3),INTENT(IN)::box
    DOUBLE PRECISION,DIMENSION(nframes),INTENT(OUT)::volume
    INTEGER,INTENT(OUT)::ortho
    DOUBLE PRECISION::x,y,z
    DOUBLE PRECISION,PARAMETER::fact_pi_div_180=1.745329251994329547437168059786927E-0002

    ortho=0
    cell=0.0d0
    volume=0.0d0

    DO ii=1,nframes

      IF ((box(ii,2,1)==0.0d0).and.(box(ii,3,1)==0.0d0).and.(box(ii,3,2)==0.0d0)) THEN
         volume(ii)=box(ii,1,1)*box(ii,2,2)*box(ii,3,3)
         cell(ii,1,1)=box(ii,1,1)
         cell(ii,2,2)=box(ii,2,2)
         cell(ii,3,3)=box(ii,3,3)
         cell(ii,1,2)=90.0d0
         cell(ii,1,3)=90.0d0
         cell(ii,2,3)=90.0d0
         ortho=1
      ELSE
         volume(ii)=box(ii,1,1)*box(ii,2,2)*box(ii,3,3)
         x=box(ii,1,1)
         y=sqrt(box(ii,2,1)**2+box(ii,2,2)**2)
         z=sqrt(box(ii,3,1)**2+box(ii,3,2)**2+box(ii,3,3)**2)
         cell(ii,1,1)=x
         cell(ii,2,2)=y
         cell(ii,3,3)=z
         cell(ii,2,3)=(acos(dot_product(box(ii,2,:),box(ii,1,:))/(x*y)))/fact_pi_div_180 ! gamma
         cell(ii,1,3)=(acos(dot_product(box(ii,3,:),box(ii,1,:))/(x*z)))/fact_pi_div_180 ! beta
         cell(ii,1,2)=(acos(dot_product(box(ii,3,:),box(ii,2,:))/(x*y)))/fact_pi_div_180 ! alpha
      END IF

    END DO

  END SUBROUTINE BOX2CELL

  SUBROUTINE WRAP (coors,box,ortho,natoms,nframes)

    INTEGER,INTENT(IN)::ortho,natoms,nframes
    double precision,DIMENSION(nframes,3,3),INTENT(IN)::box
    double precision,dimension(nframes,natoms,3),intent(INOUT)::coors

    DOUBLE PRECISION::Lx,Ly,Lz,x,y,z

    INTEGER::ii,jj,nn

    DO jj=1,nframes

      Lx=box(jj,1,1)
      Ly=box(jj,2,2)
      Lz=box(jj,3,3)

      IF (ortho==1) THEN

         DO ii=1,natoms
            x=coors(jj,ii,1)
            y=coors(jj,ii,2)
            z=coors(jj,ii,3)
            IF (x<0.0d0) THEN
               nn=CEILING(abs(x)/Lx)
               x=x+nn*Lx
            ELSE IF (x>=Lx) THEN
               nn=INT(x/Lx)
               x=x-nn*Lx
            END IF
            IF (y<0.0d0) THEN
               nn=CEILING(abs(y)/Ly)
               y=y+nn*Ly
            ELSE IF (y>=Ly) THEN
               nn=INT(y/Ly)
               y=y-nn*Ly
            END IF
            IF (z<0.0d0) THEN
               nn=CEILING(abs(z)/Lz)
               z=z+nn*Lz
            ELSE IF (z>=Lz) THEN
               nn=INT(z/Lz)
               z=z-nn*Lz
            END IF
            coors(jj,ii,:)=(/x,y,z/)
         END DO

      ELSE

         print*,'WRAP function not implemented for not orthorhombic unit cells.'

      END IF

    END DO

  END SUBROUTINE WRAP

  SUBROUTINE UNWRAP (coors, atoms_indices_serialized, atoms_series_starts, bonds_indices_serialized, &
          & bonds_series_starts, frame_indices, box, inv, ortho, n_frames, n_atoms, n_atoms_indices, &
          & n_atoms_series, n_bonds_indices, n_bonds_series, n_frames_indices)

    IMPLICIT NONE

    INTEGER,INTENT(IN):: n_frames, n_atoms, n_atoms_indices, n_atoms_series
    INTEGER,INTENT(IN):: n_bonds_indices, n_bonds_series, n_frames_indices, ortho
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms,3),INTENT(INOUT)::coors
    INTEGER,DIMENSION(n_atoms_indices),INTENT(IN)::atoms_indices_serialized ! Molecules in this case
    INTEGER,DIMENSION(n_atoms_series+1),INTENT(IN)::atoms_series_starts
    INTEGER,DIMENSION(n_bonds_indices),INTENT(IN)::bonds_indices_serialized
    INTEGER,DIMENSION(n_bonds_series+1),INTENT(IN)::bonds_series_starts
    INTEGER,DIMENSION(n_frames_indices),INTENT(IN)::frame_indices
    DOUBLE PRECISION,DIMENSION(n_frames,3,3),INTENT(IN)::box,inv

    INTEGER:: ii, jj, kk, ll
    INTEGER:: molecule_start, molecule_end, molecule_natoms
    INTEGER:: n_left, n_storage
    INTEGER:: atom_id1, atom_id2
    LOGICAL,DIMENSION(n_atoms)::aux_filter
    DOUBLE PRECISION,DIMENSION(3)::vect_aux
    INTEGER,DIMENSION(:),ALLOCATABLE::left,storage
    INTEGER::frame_index

    aux_filter(:)=.FALSE.

    DO ii=1,n_atoms_series

        molecule_start=atoms_series_starts(ii)+1
        molecule_end=atoms_series_starts(ii+1)
        molecule_natoms=molecule_end-molecule_start+1

        n_left=1
        ALLOCATE(left(n_left),storage(molecule_natoms))
        jj=atoms_indices_serialized(molecule_start)+1
        left(1)=jj
        aux_filter(jj)=.TRUE.

        DO WHILE (n_left>0)
            n_storage=0
            DO jj=1,n_left
                atom_id1 = left(jj)
                DO kk=bonds_series_starts(atom_id1)+1,bonds_series_starts(atom_id1+1)
                    atom_id2=bonds_indices_serialized(kk)+1
                    IF (aux_filter(atom_id2).eqv..FALSE.) THEN
                        DO ll=1,n_frames_indices
                            frame_index=frame_indices(ll)+1
                            vect_aux(:)=coors(frame_index,atom_id2,:)-coors(frame_index,atom_id1,:)
                            CALL PBC(vect_aux,box(frame_index,:,:),inv(frame_index,:,:),ortho)
                            coors(frame_index,atom_id2,:)=coors(frame_index,atom_id1,:)+vect_aux(:)
                        END DO
                        aux_filter(atom_id2)=.TRUE.
                        n_storage=n_storage+1
                        storage(n_storage)=atom_id2
                    END IF
                END DO
            END DO
            n_left=n_storage
            IF (n_left>0) THEN
                DEALLOCATE(left)
                ALLOCATE(left(n_left))
                left(:)=storage(1:n_left)
            END IF
        END DO

        DEALLOCATE(left,storage)

    END DO

  END SUBROUTINE UNWRAP

  SUBROUTINE MINIMUM_IMAGE_CONVENTION(coors, reference_coors, atoms_indices_serialized, atoms_series_starts, frame_indices, &
          & box, inv, ortho, n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frames_indices)

    IMPLICIT NONE

    INTEGER,INTENT(IN):: n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frames_indices, ortho

    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms,3),INTENT(INOUT)::coors
    DOUBLE PRECISION,DIMENSION(n_frames,1,3),INTENT(INOUT)::reference_coors
    INTEGER,DIMENSION(n_atoms_indices),INTENT(IN)::atoms_indices_serialized ! Molecules in this case
    INTEGER,DIMENSION(n_atoms_series+1),INTENT(IN)::atoms_series_starts
    INTEGER,DIMENSION(n_frames_indices),INTENT(IN)::frame_indices
    DOUBLE PRECISION,DIMENSION(n_frames,3,3),INTENT(IN)::box,inv

    INTEGER::ii,jj,ll
    INTEGER::atom_index, frame_index
    DOUBLE PRECISION,DIMENSION(n_frames,1,3)::com,vect_aux
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms_series,3)::centers_molecules

    centers_molecules = GEOMETRICAL_CENTER(coors, atoms_indices_serialized, atoms_series_starts, frame_indices, &
        & n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frames_indices)

    DO ii=1, n_atoms_series
        com(:,1,:) = centers_molecules(:,ii,:)
        vect_aux = com - reference_coors
        CALL PBC_TENSOR(vect_aux, box, inv, ortho, n_frames, n_atoms)
        DO jj=atoms_series_starts(ii)+1, atoms_series_starts(ii+1)
            atom_index = atoms_indices_serialized(jj)+1
            DO ll=1,n_frames_indices
                frame_index=frame_indices(ll)+1
                coors(frame_index, atom_index, :) = coors(frame_index, atom_index, :) - com(frame_index, 1, :) +&
                & vect_aux(frame_index, 1, :) + reference_coors(frame_index, 1, :)
            END DO
        END DO
    END DO

  END SUBROUTINE MINIMUM_IMAGE_CONVENTION

END MODULE MODULE_BOX

max_stars_repo_name: LMMV/MolModMT

max_issues_repo_name: LMMV/MolModMT

max_forks_repo_name: LMMV/MolModMT