text:
*******************************************
*******************************************

!           GRID-1.00

*******************************************
*******************************************

      subroutine grid

      include 'param-1.00.inc'
      include 'com-1.00.inc'

!     s grid parameters: gridding for parallel dynamics

      real s(nz,nf)
      real qs(nz,nf)
      real xs(nz,nf),ys(nz,nf),zs(nz,nf)
      real brs(nz,nf)

!     p grid parameters: gridding for perpendicular dynamics

      real qp(nzp1,nfp1),pp(nzp1,nfp1)
      real xp(nzp1,nfp1,2),yp(nzp1,nfp1,2),zp(nzp1,nfp1,2)
      real brp(nzp1,nfp1),blatp(nzp1,nfp1),blonp(nzp1,nfp1)
      real grp(nzp1,nfp1),glatp(nzp1,nfp1),glonp(nzp1,nfp1)
      real altp(nzp1,nfp1)

      do j = 1,nf
        k    = nf + 1 - j
        dx   = float(j-1)/float(nf-1)
        x2   = dx * sinh ( gamp )
        x3   = alog ( x2 + sqrt ( x2**2 + 1. ) ) / gamp
        grad_inp(k) = rmin + ( rmax - rmin ) * ( 1. - x3 )
      enddo

      call igrf_sub(glon_in)

! MS: Incorporate eccentric dipole fit to IGRF.  The fit was done by
! taken swaths in geographic longitude and fitting the best eccentric
! dipole.  We want the simulation plane to go through the point
! specified in the namelist, but the fit is done in terms of
! geographic longitude at the equator.  So iterate: Take a guess at
! an initial geographic longitude and find the eccentric dipole
! parameters. Find the magnetic longitude (use gtob_norm) of the
! namelist point and a point at the equator.  Are the two magnetic
! longitudes the same?  If not, adjust the guess for the geographic
! longitude and try again.  Rinse, repeat, wipe hands on pants.

! MS: Note that since the magnetic field is tilted (and SAMI models
! one magnetic longitude) points far away from the equator will not be
! correct.  They are at different geographic longitudes so their IGRF
! parameters will be different (and hence the g to b to g conversions
! will be different). I'm not yet sure how significant a problem this
! is.

       phitemp = 0.

       do j = 1,20
         call igrf_sub(phitemp)
         call gtob ( temp1,temp2,temp3,grad_in+re,glat_in,glon_in )
         call gtob ( brad,blon0,blat,grad_in+re,0.,phitemp )
         phitemp = phitemp + (temp2-blon0) + 360.
         phitemp = mod(phitemp,360.)
       enddo

! MS: Check for convergence

       if (abs(temp2-blon0) .gt. 1.e-2) then
         print *,'Failure to converge to initial magnetic latitude.'
         stop
       endif

      call gtob (  brad,blon0,blat
     .            ,grad_in+re
     .            ,glat_in
     .            ,glon_in          )

      pvalue0  = brad / re / cos(blat*po180) ** 2
      rgmin    = altmin + re

      call btog ( pvalue0*re,blon0,0.
     .           ,grad_max
     .           ,glat_max
     .           ,glon_max             )

!     s grid

      do j = 1,nf

        i    = 0
        delp = .01 / re
        delg = grad_inp(j)

        do while ( abs(delg) .gt. .01 )
          if ( grad_inp(j) + re .lt. grad_max .and. j .eq. 1 ) then
                pvalue = pvalue0 - i * delp
          else
                pvalue = pvalue0 + i * delp
          endif
          call btog ( pvalue*re,blon0,0.
     .               ,grad
     .               ,glat
     .               ,glon                 )
          delg = grad_inp(j) + re - grad
          i = i + 1
        enddo

        pvalue0 = pvalue

!        print *, ' j,grad_inp(j),pvalue',j,grad_inp(j),pvalue

!       north minimum: find q value at northern most point

        blata   = 0.
        blatb   = 80.
        blatc   = 0.5 * ( blata + blatb )
        rba     = pvalue * re * cos ( blata * po180 ) ** 2
        rbb     = pvalue * re * cos ( blatb * po180 ) ** 2
        rbc     = pvalue * re * cos ( blatc * po180 ) ** 2
        call  btog ( rba,blon0,blata,rga,glat,glon )
        call  btog ( rbb,blon0,blatb,rgb,glat,glon )
        call  btog ( rbc,blon0,blatc,rgc,glat,glon )
        delrg   = .01
        delrc   = abs ( rgc - rgmin )
        qminn   = ( re / rbc ) ** 2 * sin ( blatc * po180 )

        iminn   = 0

        do while ( delrc .gt. delrg .and. iminn .lt. 20)
          iminn   = iminn + 1
          if ( rgc .lt. rgmin ) blatb = blatc
          if ( rgc .gt. rgmin ) blata = blatc
          blatc   = 0.5 * ( blata + blatb )
          rbc     = pvalue * re * cos ( blatc * po180 ) ** 2
          call  btog ( rbc,blon0,blatc,rgc,glat,glon )
          delrc   = abs ( rgc - rgmin )
          qminn   = ( re / rbc ) ** 2 * sin ( blatc * po180 )
        enddo

!       south minimum: find q value at southern most point

        blata   = 0.
        blatb   = -80.
        blatc   = 0.5 * ( blata + blatb )
        rba     = pvalue * re * cos ( blata * po180 ) ** 2
        rbb     = pvalue * re * cos ( blatb * po180 ) ** 2
        rbc     = pvalue * re * cos ( blatc * po180 ) ** 2
        call  btog ( rba,blon0,blata,rga,glat,glon )
        call  btog ( rbb,blon0,blatb,rgb,glat,glon )
        call  btog ( rbc,blon0,blatc,rgc,glat,glon )

        delrg   = .01
        delrc   = abs ( rgc - rgmin )
        qmins   = ( re / rbc ) ** 2 * sin ( blatc * po180 )

        imins   = 0

        do while ( delrc .gt. delrg )
          imins   = imins + 1
          if ( rgc .lt. rgmin ) blatb = blatc
          if ( rgc .gt. rgmin ) blata = blatc
          blatc   = 0.5 * ( blata + blatb )
          rbc     = pvalue * re * cos ( blatc * po180 ) ** 2
          call  btog ( rbc,blon0,blatc,rgc,glat,glon )
          delrc   = abs ( rgc - rgmin )
          qmins   = ( re / rbc ) ** 2 * sin ( blatc * po180 )
          rbms    = rbc
        enddo

!       determine distribution of points (q,s) along field line
!       see red book p. 249 (millward et al.)

        r    = rbms

        nzh   = (nz-1)/2
        delqs = qmins
        delqn = qminn

!       firs half of field line

        do i = 1,nzh+1

          x    = -1. + ( float(i) - 1 ) / float(nzh)
          xx   = x * sinh ( gams * qmins )
          qtmp = -alog ( xx + sqrt ( xx**2 + 1. ) ) / gams

!         x    = -1. + 2. * ( float(i) - 1 ) / float(nz-1)
!         xx   = x/2. * ( sinh(gams*qminn) * (x+1) +
!     .                   sinh(gams*qmins) * (x-1) )
!         qtmp = alog ( xx + sqrt ( xx**2 + 1. ) ) / gams

          r = re*qp_solve(qtmp,pvalue)

          qs(i,j)       = qtmp
          ps(i,j)       = pvalue
          brs(i,j)      = r
          br_norm       = r / re
          blats(i,j)  = asin ( qs(i,j) * br_norm ** 2 ) * rtod
          call btog ( brs(i,j),blon0,blats(i,j),
     .                grs(i,j),glats(i,j),glons(i,j) )

!         we need to calculate some quantities which depend on the usual
!         angular definition of theta

          theta = acos ( qs(i,j) * (brs(i,j)/re) ** 2 )  ! theta is in radians

!         get the dimensionless magnetic field: bm = b/b0
!         and sine / cosine of dip angle

          bms(i,j) = sqrt ( ( 2. * cos(theta) ) ** 2
     .                           + sin(theta)   ** 2 ) / br_norm ** 3

!          sindips(i,j) = 2.* cos (theta) / ( bms(i,j) * br_norm **3 )
!          cosdips(i,j) =     sin (theta) / ( bms(i,j) * br_norm **3 )

          call vector( brs(i,j),blon0,blats(i,j),
     .                 arg(i,j),athg(i,j),aphig(i,j) )

!         we define the dimensional grid points given by s = q * re

          s(i,j)  = qs(i,j)  * re
          alts(i,j) = grs(i,j) - re   ! altitude above earth

!         grid in cartesian coordinates

          xs(i,j)    = brs(i,j) * cos( blats(i,j)*po180 ) *
     .                            sin( blon0     *po180 )
          ys(i,j)    = brs(i,j) * sin( blats(i,j)*po180 )
          zs(i,j)    = brs(i,j) * cos( blats(i,j)*po180 ) *
     .                            cos( blon0     *po180 )

        enddo

!       second half of field line

        do i = nz,nzh+2,-1

          x    = ( float(i) - ( float(nzh) + 1.) ) / float(nzh)
          xx   = x * sinh ( gams * qminn )
          qtmp = alog ( xx + sqrt ( xx**2 + 1. ) ) / gams

!         x    = -1. + 2. * ( float(i) - 1 ) / float(nz-1)
!         xx   = x/2. * ( sinh(gams*qminn) * (x+1) +
!     .                   sinh(gams*qmins) * (x-1) )
!         qtmp = alog ( xx + sqrt ( xx**2 + 1. ) ) / gams

          r = re*qp_solve(qtmp,pvalue)

          qs(i,j)       = qtmp
          ps(i,j)       = pvalue
          brs(i,j)      = r
          br_norm       = r / re
          blats(i,j)  = asin ( qs(i,j) * br_norm ** 2 ) * rtod
          call btog ( brs(i,j),blon0,blats(i,j),
     .                grs(i,j),glats(i,j),glons(i,j) )

!         we need to calculate some quantities which depend on the usual
!         angular definition of theta

          theta = acos ( qs(i,j) * (brs(i,j)/re) ** 2 )  ! theta is in radians

!         get the dimensionless magnetic field: bm = b/b0
!         and sine / cosine of dip angle

          bms(i,j) = sqrt ( ( 2. * cos(theta) ) ** 2
     .                           + sin(theta)   ** 2 ) / br_norm ** 3

!          sindips(i,j) = 2.* cos (theta) / ( bms(i,j) * br_norm **3 )
!          cosdips(i,j) =     sin (theta) / ( bms(i,j) * br_norm **3 )

          call vector( brs(i,j),blon0,blats(i,j),
     .                 arg(i,j),athg(i,j),aphig(i,j) )

!         we define the dimensional grid points given by s = q * re

          s(i,j)  = qs(i,j)  * re
          alts(i,j) = grs(i,j) - re   ! altitude above earth

!         grid in cartesian coordinates

          xs(i,j)    = brs(i,j) * cos( blats(i,j)*po180 ) *
     .                            sin( blon0     *po180 )
          ys(i,j)    = brs(i,j) * sin( blats(i,j)*po180 )
          zs(i,j)    = brs(i,j) * cos( blats(i,j)*po180 ) *
     .                            cos( blon0     *po180 )

        enddo

      enddo ! j loop (number of field lines)

!     following distances are in cm

      do j = 1,nf
        do i = 2,nz-1
          ds(i,j)   = ( s(i,j)   - s(i-1,j) ) * 1.e5
          d2s(i,j)  = ( s(i+1,j) - s(i-1,j) ) * 1.e5
          d22s(i,j) = .5 * d2s(i,j)
        enddo
        ds(1,j)    = ds(2,j)
        ds(nz,j)   = ds(nz-1,j)
        d2s(1,j)   = d2s(2,j)
        d2s(nz,j)  = d2s(nz-1,j)
        d22s(1,j)  = d22s(2,j)
        d22s(nz,j) = d22s(nz-1,j)
      enddo

!     calculate dels: grid length along field line using cartesian geometry
!     and convert to cm from km

      do j = 1,nf
        do i = 1,nz-1
          x2  = ( xs(i+1,j) - xs(i,j) ) ** 2
          y2  = ( ys(i+1,j) - ys(i,j) ) ** 2
          zz2 = ( zs(i+1,j) - zs(i,j) ) ** 2
          dels(i,j) = sqrt ( x2 + y2 + zz2 ) * 1.e5
        enddo
        dels(nz,j) = dels(nz-1,j)
      enddo

!     now do gridding for p mesh

!     calculate new set of q's  (qp)

      do j = 1,nf
        do i = 1,nz-1
          qp(i+1,j) = 0.5 * ( qs(i,j) + qs(i+1,j) )
        enddo
        qp(1,j)    = qs(1,j)  + 0.5 * ( qs(1,j) - qs(2,j) )
        qp(nzp1,j) = qs(nz,j) + 0.5 * ( qs(nz,j) - qs(nz-1,j) )
      enddo
      do i = 1,nzp1
        qp(i,nfp1) = qp(i,nf) + ( qp(i,nf) - qp(i,nf-1) )
      enddo

      do i = 1,nz
        do j = 1,nf-1
          pp(i,j+1) = 0.5 * ( ps(i,j) + ps(i,j+1) )
        enddo
        pp(i,1)    = ps(i,1)  + 0.5 * ( ps(i,1)  - ps(i,2)    )
        pp(i,nfp1) = ps(i,nf) + 0.5 * ( ps(i,nf) - ps(i,nf-1) )
      enddo
      do j = 1,nfp1
        pp(nzp1,j) = pp(nz,j) + 0.5 * ( pp(nz,j) - pp(nz-1,j) )
      enddo

!     now calculate the grid for nfp1 to 1

      do j = nfp1,1,-1
        do i = 1,nzp1
          pvalue   = pp(i,j)
          if ( j .eq. nfp1 .and. i .ne. nzp1 ) then
            r    = brs(i,nf)
          elseif ( j .eq. nfp1 .and. i .eq. nzp1 ) then
            r    = brs(nz,nf)
          else
            r    = brp(i,j+1)
          endif
          qtmp = qp(i,j)
          r    = re*qp_solve(qtmp,pvalue)

          brp(i,j)    = r
          br_norm     = r / re
          blatp(i,j)  = asin ( qp(i,j) * br_norm ** 2 ) * rtod
          blonp(i,j)  = blon0

          call btog ( brp(i,j),blonp(i,j),blatp(i,j),
     .                grp(i,j),glatp(i,j),glonp(i,j) )

          altp(i,j)   = grp(i,j) - re

          delphi =  5.
          blon1  =  blon0 - delphi
          if ( blon1 .ge. 360. ) blon1 = blon1 - 360.
          blon2  =  blon0 + delphi
          if ( blon2 .ge. 360. ) blon2 = blon2 - 360.

!         grid in cartesian coordinates

          xp(i,j,1)   = brp(i,j) * cos( blatp(i,j) * po180 ) *
     .                             sin( blon1      * po180 )
          yp(i,j,1)   = brp(i,j) * sin( blatp(i,j) * po180 )
          zp(i,j,1)   = brp(i,j) * cos( blatp(i,j) * po180 ) *
     .                             cos( blon1      * po180 )

          xp(i,j,2)   = brp(i,j) * cos( blatp(i,j) * po180 ) *
     .                             sin( blon2      * po180 )
          yp(i,j,2)   = brp(i,j) * sin( blatp(i,j) * po180 )
          zp(i,j,2)   = brp(i,j) * cos( blatp(i,j) * po180 ) *
     .                             cos( blon2      * po180 )
        enddo
      enddo

!     calculate geometric parameters: cell area and volume

!     vol is cell volume

      call volume ( xp,yp,zp )

!     areap is cell face area in j-direction (p)
!     areas is cell face area in i-direction (s)

      call area ( xp,yp,zp )

!     xdels is line distance in i-direction (s)
!     xdelp is line distance in j-direction (p)
!     xdelh is line distance in k-direction (phi)

      call line ( xp,yp,zp )

!     normal: calculates normal to cell face in s-direction
!     vnormal: calcuates e x b direction

      call  normal ( xp,yp,zp )
      call vnormal ( qs,blon0,brs )

      end

*******************************************
*******************************************

!            gtob

*******************************************
*******************************************

      subroutine gtob(brad,blond,blatd,grad,glatd,glond)

!     conversion from geographic to
!     offset centered dipole coordinates
!     brad: radius in the offset dipole system
!     grad: radius in the geocentric system
!     (g joyce june 1998)

      include 'param-1.00.inc'
      include 'com-1.00.inc'

!     coordinates of dipole in km relative to center of earth

!      x0 = -392.
!      y0 =  258.
!      z0 =  179.

!     convert to radians

      glat = glatd * po180
      glon = glond * po180

!     rotate in longitude

      xg = grad * cos ( glat ) * cos ( glon - plon )
      yg = grad * cos ( glat ) * sin ( glon - plon )
      zg = grad * sin ( glat )

!     rotate in latitude

      xmm = xg * sin ( plat ) - zg * cos ( plat )
      ymm = yg
      zmm = xg * cos ( plat ) + zg * sin ( plat )

!     calculate offset in geographic polar (p. 248 redbook)

      d = sqrt(x0**2 + y0**2 + z0**2)
      thd = acos( z0/d )
      phd = atan2( y0,x0 )

!     change to colatitude

      cotn = pie/2. - plat
      phin = plon

!     calculate offset in cartesian rotated system

      thp = acos(cos(cotn)*cos(thd)
     .     + sin(cotn)*sin(thd)*cos(phd-phin))
      cosphp = ( cos(cotn)*cos(thp) - cos(thd) )
     .         / (sin(cotn) * sin(thp))
      sinphp = sin(thd) * sin(phd-phin) / sin(thp)

      x0p = d*sin(thp)*cosphp
      y0p = d*sin(thp)*sinphp
      z0p = d*cos(thp)

!     geomagnetic latitude and longitude in degrees

      xm = xmm - x0p
      ym = ymm - y0p
      zm = zmm - z0p

!     magnetic lat and long converted to degrees

      brad  = sqrt ( xm ** 2 + ym ** 2 + zm ** 2 )
      blatd = asin ( zm / brad ) * rtod
      blond = ( atan2 ( ym/brad,xm/brad ) ) * rtod
      if (blond .lt. -0.01) blond = blond+360.

      return
      end

*******************************************
*******************************************

!            btog

*******************************************
*******************************************

      subroutine btog ( brad,blond,blatd,grad,glatd,glond )

!     conversion from centered dipole to geographic coordinates
!     plat,plon =  coords of north magnetic pole (in param2d.9.14e.inc)
!     the geomagnetic longitude is measured east from the
!     geomagnetic prime meridian at 291 degrees east geographic.
!     brad: radius in the offset dipole frame
!     grad: radius in the geocentric frame
!     (g joyce june 1998)

      include 'param-1.00.inc'
      include 'com-1.00.inc'

      real brad,blond,blatd,grad,glatd,glond

!     coordinates of dipole in km relatvie to center of earth

!      x0 = -392.
!      y0 =  258.
!      z0 =  179.

!     convert magnetic lat and long to radians

      blonr = blond * po180
      blatr = blatd * po180

!     position of point in geomagnetic coords
!     first get xm ym zm in the eccentric dipole system

      xm = brad *cos ( blatr ) * cos ( blonr )
      ym = brad *cos ( blatr ) * sin ( blonr )
      zm = brad *sin ( blatr )

!     calculate offset in geographic polar (p. 248 redbook)

      d = sqrt(x0**2 + y0**2 + z0**2)
      thd = acos( z0/d )
      phd = atan2( y0,x0 )

!     change to colatitude

      cotn = pie/2. - plat
      phin = plon

!     calculate offset in cartesian rotated system

      thp = acos(cos(cotn)*cos(thd)
     .      + sin(cotn)*sin(thd)*cos(phd-phin))
      cosphp = ( cos(cotn)*cos(thp) - cos(thd) )
     .         / (sin(cotn) * sin(thp))
      sinphp = sin(thd) * sin(phd-phin) / sin(thp)

      x0p = d*sin(thp)*cosphp
      y0p = d*sin(thp)*sinphp
      z0p = d*cos(thp)

!     next shift to the tilted dipole

      xmm = xm + x0p
      ymm = ym + y0p
      zmm = zm + z0p

!     r is invariant under the rotations of the tilted dipole

      grad = sqrt ( xmm ** 2 + ymm ** 2 + zmm ** 2 )

!     rotate coords in north-south direction

      xg =  xmm * sin ( plat ) + zmm * cos ( plat )
      yg =  ymm
      zg = -xmm * cos ( plat ) + zmm * sin ( plat )

!     geographic latitude and longitude converted to degrees

      glatd = asin ( zg / grad ) * rtod
      glond = ( plon + atan2 ( yg/grad,xg/grad ) ) * rtod

      if (glond .ge. 360) glond = glond-  360.

      return
      end

*******************************************
*******************************************

!             volume

*******************************************
*******************************************

        subroutine volume(x,y,z)

!       calculate cell volume
!       break each cell into
!       twelve tetrahedrons and use the formula:
!           V = (1/6) a . ( b x c )
!       where
!           a: vector from A to B
!           b: vector from A to C
!           c: vector from A to D
!       and node A is the 'midpoint' coordinate

        include 'param-1.00.inc'
        include 'com-1.00.inc'

        real x(nzp1,nfp1,2),y(nzp1,nfp1,2),z(nzp1,nfp1,2)
        real voli(nz,nf),volj(nz,nf),volk(nz,nf)

!       volume from sidei

        do k = 1,1
          do j = 1,nf
            do i = 1,nz

              xmid  = 0.5 * ( x(i,j,k) + x(i+1,j+1,k+1) )
              ymid  = 0.5 * ( y(i,j,k) + y(i+1,j+1,k+1) )
              zmid  = 0.5 * ( z(i,j,k) + z(i+1,j+1,k+1) )

              ax1 = x(i,j,k) - xmid
              ay1 = y(i,j,k) - ymid
              az1 = z(i,j,k) - zmid

              bx1 = x(i,j+1,k) - xmid
              by1 = y(i,j+1,k) - ymid
              bz1 = z(i,j+1,k) - zmid

              cx1 = x(i,j,k+1) - xmid
              cy1 = y(i,j,k+1) - ymid
              cz1 = z(i,j,k+1) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v1 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i,j+1,k+1) - xmid
              ay2 = y(i,j+1,k+1) - ymid
              az2 = z(i,j+1,k+1) - zmid

              v2 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              ax1 = x(i+1,j,k) - xmid
              ay1 = y(i+1,j,k) - ymid
              az1 = z(i+1,j,k) - zmid

              bx1 = x(i+1,j+1,k) - xmid
              by1 = y(i+1,j+1,k) - ymid
              bz1 = z(i+1,j+1,k) - zmid

              cx1 = x(i+1,j,k+1) - xmid
              cy1 = y(i+1,j,k+1) - ymid
              cz1 = z(i+1,j,k+1) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v3 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i+1,j+1,k+1) - xmid
              ay2 = y(i+1,j+1,k+1) - ymid
              az2 = z(i+1,j+1,k+1) - zmid

              v4 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              voli(i,j) = v1 + v2 + v3 + v4

            enddo
          enddo
        enddo

!       volume from sidej

        do k = 1,1
          do j = 1,nf
            do i = 1,nz

              xmid = 0.5 * ( x(i,j,k) + x(i+1,j+1,k+1) )
              ymid = 0.5 * ( y(i,j,k) + y(i+1,j+1,k+1) )
              zmid = 0.5 * ( z(i,j,k) + z(i+1,j+1,k+1) )

              ax1 = x(i,j,k) - xmid
              ay1 = y(i,j,k) - ymid
              az1 = z(i,j,k) - zmid

              bx1 = x(i+1,j,k) - xmid
              by1 = y(i+1,j,k) - ymid
              bz1 = z(i+1,j,k) - zmid

              cx1 = x(i,j,k+1) - xmid
              cy1 = y(i,j,k+1) - ymid
              cz1 = z(i,j,k+1) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v1 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i+1,j,k+1) - xmid
              ay2 = y(i+1,j,k+1) - ymid
              az2 = z(i+1,j,k+1) - zmid

              v2 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              ax1 = x(i,j+1,k) - xmid
              ay1 = y(i,j+1,k) - ymid
              az1 = z(i,j+1,k) - zmid

              bx1 = x(i+1,j+1,k) - xmid
              by1 = y(i+1,j+1,k) - ymid
              bz1 = z(i+1,j+1,k) - zmid

              cx1 = x(i,j+1,k+1) - xmid
              cy1 = y(i,j+1,k+1) - ymid
              cz1 = z(i,j+1,k+1) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v3 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i+1,j+1,k+1) - xmid
              ay2 = y(i+1,j+1,k+1) - ymid
              az2 = z(i+1,j+1,k+1) - zmid

              v4 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              volj(i,j) = v1 + v2 + v3 + v4

            enddo
          enddo
        enddo

!       volume from sidek

        do k = 1,1
          do j = 1,nf
            do i = 1,nz

              xmid = 0.5 * ( x(i,j,k) + x(i+1,j+1,k+1) )
              ymid = 0.5 * ( y(i,j,k) + y(i+1,j+1,k+1) )
              zmid = 0.5 * ( z(i,j,k) + z(i+1,j+1,k+1) )

              ax1 = x(i,j,k) - xmid
              ay1 = y(i,j,k) - ymid
              az1 = z(i,j,k) - zmid

              bx1 = x(i+1,j,k) - xmid
              by1 = y(i+1,j,k) - ymid
              bz1 = z(i+1,j,k) - zmid

              cx1 = x(i,j+1,k) - xmid
              cy1 = y(i,j+1,k) - ymid
              cz1 = z(i,j+1,k) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v1 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i+1,j+1,k) - xmid
              ay2 = y(i+1,j+1,k) - ymid
              az2 = z(i+1,j+1,k) - zmid

              v2 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              ax1 = x(i,j,k+1) - xmid
              ay1 = y(i,j,k+1) - ymid
              az1 = z(i,j,k+1) - zmid

              bx1 = x(i+1,j,k+1) - xmid
              by1 = y(i+1,j,k+1) - ymid
              bz1 = z(i+1,j,k+1) - zmid

              cx1 = x(i,j+1,k+1) - xmid
              cy1 = y(i,j+1,k+1) - ymid
              cz1 = z(i,j+1,k+1) - zmid

              dx1 =    by1 * cz1 - bz1 * cy1
              dy1 = -( bx1 * cz1 - bz1 * cx1 )
              dz1 =    bx1 * cy1 - by1 * cx1

              v3 = 0.166667 *
     .             abs(( ax1 * dx1 + ay1 * dy1 + az1 * dz1 ))

              ax2 = x(i+1,j+1,k+1) - xmid
              ay2 = y(i+1,j+1,k+1) - ymid
              az2 = z(i+1,j+1,k+1) - zmid

              v4 = 0.166667 *
     .             abs(( ax2 * dx1 + ay2 * dy1 + az2 * dz1 ))

              volk(i,j) = v1 + v2 + v3 + v4

            enddo
          enddo
        enddo

        do j = 1,nf
          do i = 1,nz
            vol(i,j) = 1.e15 *
     .                 ( voli(i,j) + volj(i,j) + volk(i,j) )

          enddo
        enddo

        return
        end

*******************************************
*******************************************

!            area

*******************************************
*******************************************

        subroutine area(x,y,z)

!       calculate areas of cell sides
!       break each quadrilateral side into
!       two triangles and use the formula:
!           A = (1/2)|a x b|
!       where
!           a: vector from A to B
!           b: vector from A to C

        include 'param-1.00.inc'
        include 'com-1.00.inc'

        real x(nzp1,nfp1,2),y(nzp1,nfp1,2),z(nzp1,nfp1,2)

!       sidei (s-direction)

        do k = 1,1
          do j = 1,nf
            do i = 1,nzp1

              ax1 = x(i,j+1,k) - x(i,j,k)
              ay1 = y(i,j+1,k) - y(i,j,k)
              az1 = z(i,j+1,k) - z(i,j,k)

              bx1 = x(i,j,k+1) - x(i,j,k)
              by1 = y(i,j,k+1) - y(i,j,k)
              bz1 = z(i,j,k+1) - z(i,j,k)

              cx1 =    ay1 * bz1 - az1 * by1
              cy1 = -( ax1 * bz1 - az1 * bx1 )
              cz1 =    ax1 * by1 - ay1 * bx1

              a1 = 0.5 * sqrt ( cx1*cx1 + cy1*cy1 + cz1*cz1 )

              ax2 = x(i,j+1,k) - x(i,j+1,k+1)
              ay2 = y(i,j+1,k) - y(i,j+1,k+1)
              az2 = z(i,j+1,k) - z(i,j+1,k+1)

              bx2 = x(i,j,k+1) - x(i,j+1,k+1)
              by2 = y(i,j,k+1) - y(i,j+1,k+1)
              bz2 = z(i,j,k+1) - z(i,j+1,k+1)

              cx2 =    ay2 * bz2 - az2 * by2
              cy2 = -( ax2 * bz2 - az2 * bx2 )
              cz2 =    ax2 * by2 - ay2 * bx2

              a2 = 0.5 * sqrt ( cx2*cx2 + cy2*cy2 + cz2*cz2 )

              areas(i,j) = ( a1 + a2 ) * 1.e10

            enddo
          enddo
        enddo

!       sidej (p-direction)

        do k = 1,1
          do j = 1,nfp1
            do i = 1,nz

              ax1 = x(i+1,j,k) - x(i,j,k)
              ay1 = y(i+1,j,k) - y(i,j,k)
              az1 = z(i+1,j,k) - z(i,j,k)

              bx1 = x(i,j,k+1) - x(i,j,k)
              by1 = y(i,j,k+1) - y(i,j,k)
              bz1 = z(i,j,k+1) - z(i,j,k)

              cx1 =    ay1 * bz1 - az1 * by1
              cy1 = -( ax1 * bz1 - az1 * bx1 )
              cz1 =    ax1 * by1 - ay1 * bx1

              a1 = 0.5 * sqrt ( cx1*cx1 + cy1*cy1 + cz1*cz1 )

              ax2 = x(i+1,j,k) - x(i+1,j,k+1)
              ay2 = y(i+1,j,k) - y(i+1,j,k+1)
              az2 = z(i+1,j,k) - z(i+1,j,k+1)

              bx2 = x(i,j,k+1) - x(i+1,j,k+1)
              by2 = y(i,j,k+1) - y(i+1,j,k+1)
              bz2 = z(i,j,k+1) - z(i+1,j,k+1)

              cx2 =    ay2 * bz2 - az2 * by2
              cy2 = -( ax2 * bz2 - az2 * bx2 )
              cz2 =    ax2 * by2 - ay2 * bx2

              a2 = 0.5 * sqrt ( cx2*cx2 + cy2*cy2 + cz2*cz2 )

              areap(i,j) = ( a1 + a2 ) * 1.e10

            enddo
          enddo
        enddo

        return
        end

*******************************************
*******************************************

!            line

*******************************************
*******************************************

        subroutine line(x,y,z)

!       calculate length of cell sides

        include 'param-1.00.inc'
        include 'com-1.00.inc'

        real x(nzp1,nfp1,2),y(nzp1,nfp1,2),z(nzp1,nfp1,2)
        real xdelh(nzp1,nfp1)

!       xdels (s-direction)

        do k = 1,2
          do j = 1,nfp1
            do i = 1,nz

              ax1 = x(i+1,j,k) - x(i,j,k)
              ay1 = y(i+1,j,k) - y(i,j,k)
              az1 = z(i+1,j,k) - z(i,j,k)

              xdels(i,j,k) = sqrt ( ax1*ax1 + ay1*ay1 + az1*az1 )
     .                       * 1.e5

            enddo
          enddo
        enddo

!       xdelp (p-direction)

        do k = 1,2
          do j = 1,nf
            do i = 1,nzp1

              ax1 = x(i,j+1,k) - x(i,j,k)
              ay1 = y(i,j+1,k) - y(i,j,k)
              az1 = z(i,j+1,k) - z(i,j,k)

              xdelp(i,j,k) = sqrt ( ax1*ax1 + ay1*ay1 + az1*az1 )
     .                       * 1.e5

            enddo
          enddo
        enddo

!       xdelh (phi-direction)

          do j = 1,nfp1
            do i = 1,nzp1

              ax1 = x(i,j,2) - x(i,j,1)
              ay1 = y(i,j,2) - y(i,j,1)
              az1 = z(i,j,2) - z(i,j,1)

              xdelh(i,j) = sqrt ( ax1*ax1 + ay1*ay1 + az1*az1 )
     .                     * 1.e5

            enddo
          enddo

        return
        end

*******************************************
*******************************************

!            normal

*******************************************
*******************************************

        subroutine normal ( x,y,z )

!       calculate unit normal direction to cell face
!       normal: c = a x b / |a x b|

        include 'param-1.00.inc'
        include 'com-1.00.inc'

        real x(nzp1,nfp1,2),y(nzp1,nfp1,2),z(nzp1,nfp1,2)

!       norms (normal to cell face in s-direction)

        do j = 1,nf
          do i = 1,nzp1

            ax1 = x(i,j+1,2) - x(i,j,1)
            ay1 = y(i,j+1,2) - y(i,j,1)
            az1 = z(i,j+1,2) - z(i,j,1)

            bx1 = x(i,j,2) - x(i,j+1,1)
            by1 = y(i,j,2) - y(i,j+1,1)
            bz1 = z(i,j,2) - z(i,j+1,1)

            cx1 =   ay1 * bz1 - az1 * by1
            cy1 = -(ax1 * bz1 - az1 * bx1)
            cz1 =   ax1 * by1 - ay1 * bx1

            ca1 = sqrt( cx1*cx1 + cy1*cy1 + cz1*cz1 )

            xnorms(i,j) = cx1/ca1
            ynorms(i,j) = cy1/ca1
            znorms(i,j) = cz1/ca1

          enddo
        enddo

!       normp (normal to cell face in p-direction)

        do j = 1,nfp1
          do i = 1,nz

            ax1 = x(i,j,2) - x(i+1,j,1)
            ay1 = y(i,j,2) - y(i+1,j,1)
            az1 = z(i,j,2) - z(i+1,j,1)

            bx1 = x(i+1,j,2) - x(i,j,1)
            by1 = y(i+1,j,2) - y(i,j,1)
            bz1 = z(i+1,j,2) - z(i,j,1)

            cx1 =   ay1 * bz1 - az1 * by1
            cy1 = -(ax1 * bz1 - az1 * bx1)
            cz1 =   ax1 * by1 - ay1 * bx1

            ca1 = sqrt( cx1*cx1 + cy1*cy1 + cz1*cz1 )

            xnormp(i,j) = cx1/ca1
            ynormp(i,j) = cy1/ca1
            znormp(i,j) = cz1/ca1
          enddo
        enddo

!       norml (normal to cell face in phi-direction: longitude)

!        do j = 1,nf
!          do i = 1,nzp1

!            ax1 = x(i+1,j+1,1) - x(i,j,1)
!            ay1 = y(i+1,j+1,1) - y(i,j,1)
!            az1 = z(i+1,j+1,1) - z(i,j,1)

!            bx1 = x(i,j+1,1) - x(i+1,j,1)
!            by1 = y(i,j+1,1) - y(i+1,j,1)
!            bz1 = z(i,j+1,1) - z(i+1,j,1)

!            cx1 =   ay1 * bz1 - az1 * by1
!            cy1 = -(ax1 * bz1 - az1 * bx1)
!            cz1 =   ax1 * by1 - ay1 * bx1

!            ca1 = sqrt( cx1*cx1 + cy1*cy1 + cz1*cz1 )

!            xnorml(i,j) = cx1/ca1
!            ynorml(i,j) = cy1/ca1
!            znorml(i,j) = cz1/ca1

!          enddo
!        enddo

        return
        end

*******************************************
*******************************************

!            vnormal

*******************************************
*******************************************

      subroutine vnormal ( qs,blon0,brs )

!     calculate e x b velocity direction
!     change p but keep q constant

      include 'param-1.00.inc'
      include 'com-1.00.inc'

      real qs(nz,nf),qss(nzp1,nf),pss(nzp1,nf)
      real brs1(nzp1,nf),blats1(nzp1,nf)
      real brs2(nzp1,nf),blats2(nzp1,nf)
      real xs1(nzp1,nf),ys1(nzp1,nf),zs1(nzp1,nf)
      real xs2(nzp1,nf),ys2(nzp1,nf),zs2(nzp1,nf)
      real brs(nz,nf)

!      open (unit=12, file='xs1.dat')
!      open (unit=13, file='ys1.dat')
!      open (unit=15, file='xs2.dat')
!      open (unit=16, file='ys2.dat')

      delps = .001

      do j = 1,nf
        do i = 1,nz-1
          qss(i+1,j) = 0.5 * ( qs(i,j) + qs(i+1,j) )
        enddo
        qss(1,j)    = qs(1,j)  + 0.5 * ( qs(1,j) - qs(2,j) )
        qss(nzp1,j) = qs(nz,j) + 0.5 * ( qs(nz,j) - qs(nz-1,j) )
      enddo

      do j = 1,nf
        do i = 1,nz
          pss(i+1,j) = ps(i,j)
        enddo
      enddo
      do j = 1,nf
        pss(1,j) = pss(2,j)
      enddo

      do j = 1,nf
        do i = 1,nz+1

!         grid 1 (at pss(i,j))

          qtmp   = qss(i,j)
          pvalue = pss(i,j)
          r = re*qp_solve(qtmp,pvalue)

          brs1(i,j)      = r
          br_norm        = r / re
          blats1(i,j)    = asin ( qss(i,j) * br_norm ** 2 ) * rtod

          xs1(i,j)    = brs1(i,j) * cos( blats1(i,j)*po180 ) *
     .                              sin( blon0      * po180 )
          ys1(i,j)    = brs1(i,j) * sin( blats1(i,j)*po180 )
          zs1(i,j)    = brs1(i,j) * cos( blats1(i,j) * po180 ) *
     .                              cos( blon0      * po180 )

!         grid 2 (at pss(i,j)+delps)

          qtmp   = qss(i,j)
          pvalue = pss(i,j) + delps
          r = re*qp_solve(qtmp,pvalue)

          brs2(i,j)      = r
          br_norm        = r / re
          blats2(i,j)    = asin ( qss(i,j) * br_norm ** 2 ) * rtod

!         grid in cartesian coordinates

          xs2(i,j)    = brs2(i,j) * cos( blats2(i,j)*po180 ) *
     .                              sin( blon0      * po180 )
          ys2(i,j)    = brs2(i,j) * sin( blats2(i,j)*po180 )
          zs2(i,j)    = brs2(i,j) * cos( blats2(i,j) * po180 ) *
     .                              cos( blon0       * po180 )

        enddo
      enddo

!     direction of e x b velocity

      do j = 1,nf
        do i = 1,nzp1
          ax1 = xs2(i,j) - xs1(i,j)
          ay1 = ys2(i,j) - ys1(i,j)
          az1 = zs2(i,j) - zs1(i,j)
          a1  = sqrt ( ax1*ax1 + ay1*ay1 + az1*az1 )
          vnx(i,j) = ax1 / a1
          vny(i,j) = ay1 / a1
          vnz(i,j) = az1 / a1
        enddo
      enddo

      do i = 1,nzp1
        vnx(i,nfp1) = vnx(i,nf)
        vny(i,nfp1) = vny(i,nf)
        vnz(i,nfp1) = vnz(i,nf)
      enddo

      return
      end

********************************************************

        function qp_solve(q,p)

        real qp_solve,q,p
        real term1,term2

!       MS: Functionally unnecessary, but kind of neat.  To get r from
!       (q,p) one needs to find the root of a fourth-order polynomial.
!       The code did this with a standard root-finding method,
!       but, of course, there is a closed-form solution.  Since the
!       polynomial has no second or third order terms, the answer
!       isn't completely ugly, and the code below gives the exact
!       result.  Note the special case for q = 0, i.e. points on
!       the magnetic equator.

!        if (p .eq. 0.) then
!          print *,'p = 0 in qp_solve'
!          stop
!        endif

!        if (q .eq. 0.) then
!          qp_solve = p
!        else
!        term1 = ( (sqrt(27.) + sqrt(27.+256.*q**2*p**4)) /
!     .             (16.*abs(q)*p**2) ) ** (1./3.)
!        term2 = 2./(abs(q)*sqrt(3.)) * (term1 - 1./term1)
!        qp_solve = 0.5*sqrt(term2)*(sqrt(2./(p*q**2*sqrt(term2**3)) -
!     .                       1.) - 1.)
!      endif

! MS: The formula is actually my third attempt.  The first had huge
! cancellation problems when q was small, the second had smaller
! problems when term0 was large. This should be well-behaved
! everywhere.  Massive algebra was involved along the way.

        term0 = 256./27.*q**2*p**4
        term1 = ( 1. + sqrt(1. + term0) )**(2./3.)
        term2 = term0**(1./3.)
        term3 = 0.5 * ( (term1**2 + term1*term2 +
     .                term2**2)/term1 )**(3./2.)
        qp_solve = p * (4.*term3) / (1.+term3) /
     .                 ( 1.+sqrt(2.*term3-1.) )

        end

!    IGRF_SUB.F
!    contributed by Paul Bernhardt (NRL)

       subroutine igrf_sub(phi0)

       include 'param-1.00.inc'
       include 'com-1.00.inc'

       phi = phi0 * pie / 180.

       x0  = -353.8674146847993 -
     .        80.0112901196596*Cos(phi) +
     .        165.9453553374371*Cos(2*phi) +
     .        96.90976911496519*Cos(3*phi) -
     .        41.234237919645906*Cos(4*phi) -
     .        5.647625148538237*Cos(5*phi) +
     .        29.81898130515779*Cos(6*phi) +
     .        27.310156089558348*Cos(7*phi) +
     .        7.017930973481185*Cos(8*phi) -
     .        7.511889133487608*Cos(9*phi) -
     .        6.720946815524082*Cos(10*phi) -
     .        0.5967604884638268*Cos(11*phi) +
     .        3.9340250001598633*Cos(12*phi) +
     .        3.232375199053169*Cos(13*phi) +
     .        0.1418342157838855*Cos(14*phi) -
     .        1.7399223106946722*Cos(15*phi) -
     .        1.2396197007828484*Cos(16*phi) +
     .        0.24319862206273882*Cos(17*phi) +
     .        1.0043588925588962*Cos(18*phi) +
     .        0.6050359806604826*Cos(19*phi) -
     .        0.2175679881607454*Cos(20*phi) -
     .        0.5970638247030868*Cos(21*phi) -
     .        0.3095886741276035*Cos(22*phi) +
     .        0.24106351981877516*Cos(23*phi) -
     .        39.488705071918595*Sin(phi) -
     .        35.12750391544765*Sin(2*phi) -
     .        52.61217403191678*Sin(3*phi) +
     .        35.372208218215505*Sin(4*phi) +
     .        62.517949955768565*Sin(5*phi) +
     .        33.70515829546757*Sin(6*phi) -
     .        6.367067330825964*Sin(7*phi) -
     .        16.84168815021257*Sin(8*phi) -
     .        8.77073458124951*Sin(9*phi) +
     .        2.4509417362535326*Sin(10*phi) +
     .        5.964643505586716*Sin(11*phi) +
     .        2.7430884455100206*Sin(12*phi) -
     .        1.469044512466587*Sin(13*phi) -
     .        2.58647152886297*Sin(14*phi) -
     .        0.9057773601920639*Sin(15*phi) +
     .        0.9500051070510113*Sin(16*phi) +
     .        1.2595308671288294*Sin(17*phi) +
     .        0.30670958561561906*Sin(18*phi) -
     .        0.6041955011153372*Sin(19*phi) -
     .        0.669980882605414*Sin(20*phi) -
     .        0.08600912796198656*Sin(21*phi) +
     .        0.4475307568863845*Sin(22*phi) +
     .        0.4538479001719651*Sin(23*phi)

       y0 = 357.71572803711445 +
     .      179.10000315892776*Cos(phi) +
     .      0.48625012094252057*Cos(2*phi) -
     .      18.222451520290683*Cos(3*phi) -
     .      74.19071886155683*Cos(4*phi) -
     .      36.523066435387676*Cos(5*phi) -
     .      1.3391741086911595*Cos(6*phi) +
     .      21.128815605008093*Cos(7*phi) +
     .      12.30309152643703*Cos(8*phi) +
     .      0.7072007210810681*Cos(9*phi) -
     .      4.253716697410791*Cos(10*phi) -
     .      1.9043150223068182*Cos(11*phi) +
     .      1.5369582665526593*Cos(12*phi) +
     .      2.322273352128386*Cos(13*phi) +
     .      0.8660859526717453*Cos(14*phi) -
     .      0.6517850740970705*Cos(15*phi) -
     .      0.8233191201701169*Cos(16*phi) -
     .      0.06587490166894897*Cos(17*phi) +
     .      0.5590511561227545*Cos(18*phi) +
     .      0.4857427492770479*Cos(19*phi) +
     .      0.003220351289183923*Cos(20*phi) -
     .      0.31210096865413983*Cos(21*phi) -
     .      0.19887427994494658*Cos(22*phi) +
     .      0.13672910587356485*Cos(23*phi) +
     .      416.65605532707093*Sin(phi) +
     .      97.64092937463448*Sin(2*phi) +
     .      89.90576641057676*Sin(3*phi) +
     .      19.563789031651595*Sin(4*phi) +
     .      46.61385007838707*Sin(5*phi) +
     .      46.23624804545088*Sin(6*phi) +
     .      19.328984217463248*Sin(7*phi) +
     .      2.4598324280936414*Sin(8*phi) -
     .      3.4043388007500432*Sin(9*phi) +
     .      2.111663037458501*Sin(10*phi) +
     .      4.723429973505645*Sin(11*phi) +
     .      3.7939791153184026*Sin(12*phi) +
     .      0.783837405015646*Sin(13*phi) -
     .      0.8951808216876126*Sin(14*phi) -
     .      0.5853291438883146*Sin(15*phi) +
     .      0.49374732010654354*Sin(16*phi) +
     .      0.9299824387332063*Sin(17*phi) +
     .      0.48650275322307124*Sin(18*phi) -
     .      0.15694935678709498*Sin(19*phi) -
     .      0.36063573138052973*Sin(20*phi) -
     .      0.08389592287106021*Sin(21*phi) +
     .      0.25573941044924925*Sin(22*phi) +
     .      0.2852013162673946*Sin(23*phi)

       z0 = 250.43165707492653 +
     .      340.66122483264905*Cos(phi) +
     .      69.95634836232071*Cos(2*phi) +
     .      44.71168620310048*Cos(3*phi) +
     .      5.572752558766106*Cos(4*phi) -
     .      1.4339461502279616*Cos(5*phi) -
     .      3.217930812956232*Cos(6*phi) -
     .      1.7751137903566045*Cos(7*phi) -
     .      1.5857495900135135*Cos(8*phi) -
     .      0.7461261912547706*Cos(9*phi) +
     .      0.14563859321432496*Cos(10*phi) +
     .      0.3315895855881013*Cos(11*phi) +
     .      0.11379367536280918*Cos(12*phi) -
     .      0.15395572288202422*Cos(13*phi) -
     .      0.19461671510187845*Cos(14*phi) -
     .      0.06361599694133094*Cos(15*phi) +
     .      0.055766203925038796*Cos(16*phi) +
     .      0.06413644582150604*Cos(17*phi) +
     .      0.0016750432363388439*Cos(18*phi) -
     .      0.04182349014199574*Cos(19*phi) -
     .      0.031182905113287293*Cos(20*phi) +
     .      0.001753974707015248*Cos(21*phi) +
     .      0.013973913699992777*Cos(22*phi) +
     .      0.001281301856688131*Cos(23*phi) -
     .      30.176529635127274*Sin(phi) +
     .      3.4340205951414946*Sin(2*phi) -
     .      6.8123984277392795*Sin(3*phi) -
     .      11.668750059562694*Sin(4*phi) -
     .      0.5810580595917775*Sin(5*phi) -
     .      1.5870410791195877*Sin(6*phi) -
     .      1.5244529990716116*Sin(7*phi) -
     .      0.559954677080775*Sin(8*phi) +
     .      0.5194385178951049*Sin(9*phi) +
     .      0.415495066076052*Sin(10*phi) +
     .      0.04521386570333395*Sin(11*phi) -
     .      0.24535417687477926*Sin(12*phi) -
     .      0.18203325048713365*Sin(13*phi) +
     .      0.018878146825439902*Sin(14*phi) +
     .      0.12360222004213386*Sin(15*phi) +
     .      0.07257334462042853*Sin(16*phi) -
     .      0.022737330338338836*Sin(17*phi) -
     .      0.05776054926047047*Sin(18*phi) -
     .      0.023788034168565193*Sin(19*phi) +
     .      0.018325492914075404*Sin(20*phi) +
     .      0.02376235757596904*Sin(21*phi) +
     .      0.002628561346701197*Sin(22*phi) -
     .      0.00943276960175865*Sin(23*phi)

       bb0 = 29229.88363541812 -
     .      1115.0281463875144*Cos(phi) +
     .      211.9965588481871*Cos(2*phi) -
     .      145.87863848381699*Cos(3*phi) -
     .      185.94991131071188*Cos(4*phi) -
     .      94.76388090038132*Cos(5*phi) -
     .      21.781240606511467*Cos(6*phi) +
     .      10.608481026851262*Cos(7*phi) +
     .      26.445003680471622*Cos(8*phi) +
     .      9.033803335738316*Cos(9*phi) -
     .      5.121597505263254*Cos(10*phi) -
     .      8.699077903900074*Cos(11*phi) -
     .      2.3758078601113084*Cos(12*phi) +
     .      3.003659043878612*Cos(13*phi) +
     .      3.3447365817081867*Cos(14*phi) +
     .      0.46321746934859276*Cos(15*phi) -
     .      1.6640994428565439*Cos(16*phi) -
     .      1.3856408382560013*Cos(17*phi) +
     .      0.06865975473204816*Cos(18*phi) +
     .      0.8650543900984624*Cos(19*phi) +
     .      0.4841643209453366*Cos(20*phi) -
     .      0.22968004188392457*Cos(21*phi) -
     .      0.38869426768468657*Cos(22*phi) -
     .      0.00489072283424378*Cos(23*phi) -
     .      343.1161166643772*Sin(phi) -
     .      214.46870730780248*Sin(2*phi) -
     .      218.4373355152967*Sin(3*phi) -
     .      139.32266779124643*Sin(4*phi) -
     .      40.479773535930704*Sin(5*phi) +
     .      21.102185585291537*Sin(6*phi) +
     .      27.763390345914974*Sin(7*phi) +
     .      1.6622501124443485*Sin(8*phi) -
     .      13.49400563169258*Sin(9*phi) -
     .      10.573966878144494*Sin(10*phi) +
     .      0.0010493967978287462*Sin(11*phi) +
     .      5.211690149433707*Sin(12*phi) +
     .      3.1952678133195556*Sin(13*phi) -
     .      1.1454324071564117*Sin(14*phi) -
     .      2.728164501523932*Sin(15*phi) -
     .      1.2850993411366254*Sin(16*phi) +
     .      0.7013933219589532*Sin(17*phi) +
     .      1.1699781168734327*Sin(18*phi) +
     .      0.301352221634545*Sin(19*phi) -
     .      0.5291265885714342*Sin(20*phi) -
     .      0.49114080833339335*Sin(21*phi) +
     .      0.050614256681022596*Sin(22*phi) +
     .      0.28881769341899144*Sin(23*phi)

       thn = 0.19622243171731607 -
     .      0.006649098360752184*Cos(phi) -
     .      0.01956517122663376*Cos(2*phi) -
     .      0.008660216931195966*Cos(3*phi) -
     .      0.006137819025589696*Cos(4*phi) -
     .      0.0075964729910884795*Cos(5*phi) -
     .      0.0016461860924103853*Cos(6*phi) +
     .      0.0007294454409291066*Cos(7*phi) +
     .      0.0023685663829071343*Cos(8*phi) +
     .      0.0003889280477075606*Cos(9*phi) -
     .      0.0006979450055063492*Cos(10*phi) -
     .      0.000885537681495029*Cos(11*phi) -
     .      0.00014878907992479928*Cos(12*phi) +
     .      0.0003679941251949079*Cos(13*phi) +
     .      0.0003130319474168302*Cos(14*phi) -
     .      0.000016633051122379423*Cos(15*phi) -
     .      0.00020931570552896087*Cos(16*phi) -
     .      0.0001303595920774504*Cos(17*phi) +
     .      0.0000388948168387536*Cos(18*phi) +
     .      0.00010454833096758408*Cos(19*phi) +
     .      0.00003792479727719437*Cos(20*phi) -
     .      0.00004613565179395065*Cos(21*phi) -
     .      0.00005090440921535995*Cos(22*phi) +
     .      7.923475578463114e-6*Cos(23*phi) -
     .      0.05678320125406135*Sin(phi) -
     .      0.023022365710781457*Sin(2*phi) -
     .      0.018030195256608213*Sin(3*phi) -
     .      0.013827655583181043*Sin(4*phi) -
     .      0.0061990077703202836*Sin(5*phi) +
     .      0.0017061925607212282*Sin(6*phi) +
     .      0.0023352555360747727*Sin(7*phi) -
     .      0.0002812394123411336*Sin(8*phi) -
     .      0.0014632759741499833*Sin(9*phi) -
     .      0.0009892215014679607*Sin(10*phi) +
     .      0.00012383435618247*Sin(11*phi) +
     .      0.0005411721277811578*Sin(12*phi) +
     .      0.00024600375235788244*Sin(13*phi) -
     .      0.0002043069395354885*Sin(14*phi) -
     .      0.0003033173283262694*Sin(15*phi) -
     .      0.00010053028847115513*Sin(16*phi) +
     .      0.0001067032704224899*Sin(17*phi) +
     .      0.00012248942798267554*Sin(18*phi) +
     .      5.7013723087397106e-6*Sin(19*phi) -
     .      0.00007782884382598787*Sin(20*phi) -
     .      0.0000523747335689761*Sin(21*phi) +
     .      0.000021462140784394788*Sin(22*phi) +
     .      0.00004552716184377109*Sin(23*phi)

       phn = -1.1840115989900575 -
     .      0.2501052377367187*Cos(phi) -
     .      0.0870108497349657*Cos(2*phi) +
     .      0.05524474825558274*Cos(3*phi) +
     .      0.06571634574967282*Cos(4*phi) +
     .      0.054497177647939085*Cos(5*phi) -
     .      0.010914652319742967*Cos(6*phi) -
     .      0.017494764756306118*Cos(7*phi) -
     .      0.001677474136905055*Cos(8*phi) +
     .      0.009094152785717784*Cos(9*phi) +
     .      0.00778609301072576*Cos(10*phi) -
     .      0.000017569943641948104*Cos(11*phi) -
     .      0.00427665106371633*Cos(12*phi) -
     .      0.0031248120532568428*Cos(13*phi) +
     .      0.0004018612433600332*Cos(14*phi) +
     .      0.002025672163296362*Cos(15*phi) +
     .      0.001115968407376618*Cos(16*phi) -
     .      0.0005287555269570307*Cos(17*phi) -
     .      0.0011232678915695546*Cos(18*phi) -
     .      0.0005000380314326315*Cos(19*phi) +
     .      0.00036696679782018067*Cos(20*phi) +
     .      0.000649640645505565*Cos(21*phi) +
     .      0.0002778176080716778*Cos(22*phi) -
     .      0.0002750904248179207*Cos(23*phi) +
     .      0.03756894419604881*Sin(phi) -
     .      0.005007302851653333*Sin(2*phi) +
     .      0.09805567331135825*Sin(3*phi) -
     .      0.015507188585201165*Sin(4*phi) -
     .      0.05239031186985234*Sin(5*phi) -
     .      0.01826942178430312*Sin(6*phi) +
     .      0.0022905823113473396*Sin(7*phi) +
     .      0.018062815520885563*Sin(8*phi) +
     .      0.00816576532607775*Sin(9*phi) -
     .      0.0018070882685123845*Sin(10*phi) -
     .      0.005789120556186378*Sin(11*phi) -
     .      0.0018598667086740701*Sin(12*phi) +
     .      0.0023061252186550446*Sin(13*phi) +
     .      0.002787011388680176*Sin(14*phi) +
     .      0.0006208190253855384*Sin(15*phi) -
     .      0.001298615463680044*Sin(16*phi) -
     .      0.0012840198090914765*Sin(17*phi) -
     .      0.00011250361770319753*Sin(18*phi) +
     .      0.00076532973526117*Sin(19*phi) +
     .      0.0006628511134355898*Sin(20*phi) -
     .      6.212956128694389e-6*Sin(21*phi) -
     .      0.000505326773031062*Sin(22*phi) -
     .      0.0004597641346084901*Sin(23*phi)

         bb0 = 1.e-5 * bb0 ! convert to gauss

!         print *,'phi0,x0,y0,z0,thn,phn,b0',phi0,x0,y0,z0,thn,phn,bb0

       plat = 0.5 * pie - thn
       plon = 2. * pie + phn

       return
       end

*************************************************

!            vector

! MS 9/15/05
! This subroutine finds the coefficients needed to calculate the
! parallel components of vectors.  The math is shamelessly copied from
! the CTIP paper in the red book.  Inputs are the magnetic coordinates
! of a point, outputs are the coefficients at that point.

      subroutine vector(brad,blond,blatd,varg,vathg,vaphig)

      include 'param-1.00.inc'
      include 'com-1.00.inc'

      real brad,blond,blatd,blonr,grad,glat,glon
      real thetaprime,cosphiprime,sinphiprime
      real sinomega,cosomega,coplat,theta
      real arm,athm,aphim,ax,ay,az,arp,athp,aphip,varg,vathg,vaphig

!     Convert blatd and plat to colatitude, blon to radians

      theta  = pie/2. - blatd*po180
      coplat = pie/2. - plat
      blonr  = blond*po180

      arm = 2.* cos (theta) /  sqrt ( ( 2. * cos(theta) ) ** 2 +
     .        sin(theta)   ** 2 )
      athm =     sin (theta) / sqrt ( ( 2. * cos(theta) ) ** 2 +
     .        sin(theta)   ** 2 )
       aphim = 0.

      ax = arm * sin(theta)*cos(blonr) + athm*cos(theta)*cos(blonr)
      ay = arm * sin(theta)*sin(blonr) + athm*cos(theta)*sin(blonr)
      az = arm * cos(theta) - athm*sin(theta)

!     Find geographic coordinates of point

      call btog(brad,blond,blatd,grad,glat,glon)
      glat = pie/2. - glat*po180
      glon = glon * po180

!     See p.243 of red book

      thetaprime = acos(cos(coplat)*cos(glat) +
     .             sin(coplat)*sin(glat)*cos(glon-plon))

!     Account for possibility of no tilt.

      if (coplat .eq. 0.) then
        cosphiprime = 1.
        sinphiprime = 0.
      else
        cosphiprime = ( cos(coplat)*cos(thetaprime) - cos(glat) ) /
     .              (sin(coplat) * sin(thetaprime))
        sinphiprime = sin(glat) * sin(glon-plon) / sin(thetaprime)
      endif

      arp = ax * sin(thetaprime)*cosphiprime +
     .      ay*sin(thetaprime)*sinphiprime +
     .      az*cos(thetaprime)
      athp = ax * cos(thetaprime)*cosphiprime +
     .        ay*cos(thetaprime)*sinphiprime -
     .        az*sin(thetaprime)
      aphip = -ax*sinphiprime + ay*cosphiprime

      sinomega = -sin(coplat)*sin(glon - plon) / sin(thetaprime)

!     MS: This formula is wrong in the red book.  Corrected here.

      cosomega = (cos(coplat) - cos(glat)*cos(thetaprime)) /
     .           (sin(thetaprime) * sin(glat))

      varg  = arp
      vathg = athp*cosomega + aphip*sinomega

!     Because of a typographical error this formula does not appear in
!     the red book.

      vaphig = -athp*sinomega + aphip*cosomega

      return
      end


text:
module constants

use parameters

implicit none

!Universal constants
real(kind=dp), parameter :: Newton_g = 6.67408d-11
real(kind=dp), parameter :: Msolar = 1.989d30
real(kind=dp), parameter :: mu = Newton_g*MBH*Msolar
real(kind=dp), parameter :: light_c = 3.0d8
real(kind=dp), parameter :: convert_m = light_c**2/(Newton_g*MBH*Msolar) !Multiply by this to go TO Natural units
real(kind=dp), parameter :: convert_s = light_c**3/(Newton_g*MBH*Msolar) !Multiply by this to go TO Natural units
real(kind=dp), parameter :: convert_spin= light_c/(Newton_g*(MBH*Msolar)**2.0_dp) !Multiply by this to go TO Natural units
real(kind=dp), parameter :: Rhor = 1.0_dp+sqrt(1.0_dp-a**2) + 1.0d-2 !Horizon + eps
real(kind=dp), PARAMETER :: electron_charge = 4.80320425D-10 !CGS
real(kind=dp), PARAMETER :: electron_mass = 9.10938356D-28 !CGS
real(kind=dp), PARAMETER :: pc = 3.0860e16_dp !meters

!MPD
real(kind=dp), parameter :: KeplerianPeriodSeconds = KeplerianPeriod*365.0_dp*24.0_dp*3600.0_dp
real(kind=dp), parameter :: SM3 =(mu*KeplerianPeriodSeconds**2.0_dp)/(4.0_dp * PI**2.0_dp)
real(kind=dp), parameter :: semi_major = convert_m*SM3**(1.0_dp/3.0_dp)
real(kind=dp), parameter :: semi_latus = semi_major * (1.0_dp - eccentricity**2.0_dp)
real(kind=dp), parameter :: ra = semi_latus/(1.0_dp - eccentricity)
real(kind=dp), parameter :: rp = semi_latus/(1.0_dp + eccentricity)
real(kind=dp), parameter :: theta_min = (90.0_dp - iota) * PI/180.0_dp !Minimum latitude reached in radians
real(kind=dp), parameter :: zMinus = cos(theta_min)**2.0_dp
real(kind=dp), parameter :: m0 = MPSR/MBH
real(kind=dp), parameter :: inertia = 0.40_dp*(MPSR*Msolar)*(RPSR*1.0d3)**2.0_dp !SI units
real(kind=dp), parameter :: s0 = convert_spin*2.0_dp*PI*inertia/p0 !magnitude of spin spatial vector in natural units
integer(kind=dp), parameter :: entries = 13 !Number of differetnai eqns 4x(position,spin,momentum + proper time)
real(kind=dp), parameter :: hs_natural = hs * convert_s
real(kind=dp), parameter :: SpinPeriod = p0 * convert_s

real(kind=dp), parameter :: RObs = 1e4_dp !RO*1e3*pc * convert_m

!MPD I/O
integer(kind=dp), parameter :: nrows = 1d4
character(len=200) :: path, MPDBinaryData, MPDFormatData, RTPath
integer(kind=dp) :: OrbitNrows, OrbitNcols

!Ray tracing

!Dark matter
real(kind=dp) :: r0 = r0_SI*1e3*pc * convert_m
real(kind=dp) :: kDM = r0_SI**3 * rho0  !Check units
integer(kind=dp) :: DM !turn DM halo on /off

!Numerical integration
!Cash-Karp parameters and integration constants
real(kind = dp) :: B21=1.0_dp/5.0_dp
real(kind = dp) :: B31 = 3.0_dp/40.0_dp , B32 = 9.0_dp/40.0_dp
real(kind = dp) :: B41 = 3.0_dp/10.0_dp, B42 = -9.0_dp/10.0_dp, B43 = 6.0_dp/5.0_dp
real(kind = dp) :: B51 = -11.0_dp/54.0_dp, B52 = 5.0_dp/2.0_dp, B53 = -70.0_dp/27.0_dp
real(kind = dp) :: B54 = 35.0_dp/27.0_dp
real(kind = dp) :: B61 = 1631.0_dp/55296.0_dp, B62 = 175.0_dp/512.0_dp, B63 = 575.0_dp/13824.0_dp
real(kind = dp) :: B64 = 44275.0_dp/110592.0_dp, B65 = 253.0_dp/4096.0_dp

real(kind = dp) :: c1 = 37.0_dp/378.0_dp, c3 = 250.0_dp/621.0_dp, c4 = 125.0_dp/594.0_dp
real(kind = dp) :: c6=512.0_dp/1771.0_dp
real(kind = dp) :: cbar1 = 2825.0_dp/27648.0_dp, cbar3 = 18575.0_dp/48384.0_dp
real(kind = dp) :: cbar4=13525.0_dp/55296.0_dp, cbar5 = 277.0_dp/14336.0_dp, cbar6 = 1.0_dp/4.0_dp

real(kind = dp) :: escal = 1.0d19 !SHOULD WE CHANGE THIS? SPEED VS ACCURACY....
real(kind = dp) :: PSHRINK = -0.25_dp,PGROW=-0.20_dp, S=0.90_dp
real(kind=dp), parameter :: dx_eps = 1.0e-22_dp !How small do we push this param?

!Some globally defined MPD variables
real(kind=dp) :: m_sq, s_sq ! mass sqaures and s squared from initial condiitons module
real(kind=dp) :: tau !Proper time. Related to the stepsize

!Christoffel sybols
real(kind=dp) ::    G0_00, G0_01, G0_02, G0_03, G0_11, G0_12, G0_13, G0_22, G0_23, G0_33,       &
                    G1_00, G1_01, G1_02, G1_03, G1_11, G1_12, G1_13, G1_22, G1_23, G1_33,       &
                    G2_00, G2_01, G2_02, G2_03, G2_11, G2_12, G2_13, G2_22, G2_23, G2_33,       &
                    G3_00, G3_01, G3_02, G3_03, G3_11, G3_12, G3_13, G3_22, G3_23, G3_33

!Riemann
REAL(KIND=dp) ::                                                                &
    R_0101, R_0102, R_0103, R_0112, R_0113, R_0123,                             &
    R_0201, R_0202, R_0203, R_0212, R_0213, R_0223,                             &
    R_0301, R_0302, R_0303, R_0312, R_0313, R_0323,                             &
    R_1201, R_1202, R_1203, R_1212, R_1213, R_1223,                             &
    R_1301, R_1302, R_1303, R_1312, R_1313, R_1323,                             &
    R_2301, R_2302, R_2303, R_2312, R_2313, R_2323,                             &
!
    R0_001, R0_002, R0_003, R0_012, R0_013, R0_023,                             &
    R0_101, R0_102, R0_103, R0_112, R0_113, R0_123,                             &
    R0_201, R0_202, R0_203, R0_212, R0_213, R0_223,                             &
    R0_301, R0_302, R0_303, R0_312, R0_313, R0_323,                             &
!
    R1_001, R1_002, R1_003, R1_012, R1_013, R1_023,                             &
    R1_101, R1_102, R1_103, R1_112, R1_113, R1_123,                             &
    R1_201, R1_202, R1_203, R1_212, R1_213, R1_223,                             &
    R1_301, R1_302, R1_303, R1_312, R1_313, R1_323,                             &
!
    R2_001, R2_002, R2_003, R2_012, R2_013, R2_023,                             &
    R2_101, R2_102, R2_103, R2_112, R2_113, R2_123,                             &
    R2_201, R2_202, R2_203, R2_212, R2_213, R2_223,                             &
    R2_301, R2_302, R2_303, R2_312, R2_313, R2_323,                             &
!
    R3_001, R3_002, R3_003, R3_012, R3_013, R3_023,                             &
    R3_101, R3_102, R3_103, R3_112, R3_113, R3_123,                             &
    R3_201, R3_202, R3_203, R3_212, R3_213, R3_223,                             &
    R3_301, R3_302, R3_303, R3_312, R3_313, R3_323

!Spin tensor

real(kind=dp) :: S_01, S_02, S_03, S_12, S_13, S_23

end module constants

ext: f

max_stars_repo_name: tomkimpson/DarkMatterRayTracing

max_issues_repo_name: tomkimpson/DarkMatterRayTracing

max_forks_repo_name: tomkimpson/DarkMatterRayTracing


text:
!##############################################################################
!# ****************************************************************************
!# <name> burgers1dcallback </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains callback functions for the Burgers-1D problem that are
!# used during the matrix/vector assembly for specifying analytical data.
!# There are three callback functions involved, which may be called depending
!# on the situation. All of them correspond to a specific interface for
!# callback functions, defined in 'intf_xxxx.inc' files.
!#
!# 1.) coeff_burgers
!#     -> Returns the coefficients for the System matrix. As the Burgers
!#        equation needs to build a linearised nonlinear matrix,
!#        this routine has to evaluate a FE function and use the FE values
!#        as coefficients in the system matrix.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientMatrixSc.inc'
!#
!# 2.) coeff_RHS
!#     -> Returns analytical values for the right hand side of the Laplace
!#        equation.
!#     -> Corresponds to the interface defined in the file
!#        'intf_coefficientVectorSc.inc'
!#
!# 3.) getBoundaryValues
!#     -> Returns analitical values on the (Dirichlet) boundary of the
!#        problem to solve.
!#     -> Corresponds to the interface defined in the file
!#        'intf_bcassembly.inc'
!#
!# </purpose>
!##############################################################################

module burgers1d_callback

  use fsystem
  use storage
  use linearsolver
  use boundary
  use derivatives
  use spatialdiscretisation
  use linearsystemscalar
  use linearsystemblock
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use feevaluation
  use collection

  implicit none

contains

! ***************************************************************************
  !<subroutine>

  subroutine coeff_burgers (rdiscretisationTrial,rdiscretisationTest,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTrial,IdofsTest,rdomainIntSubset, &
                  Dcoefficients,rcollection)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration

  !<description>
    ! This subroutine is called during the matrix assembly. It has to compute
    ! the coefficients in front of the terms of the bilinear form.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the bilinear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the bilinear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.; trial space.
    type(t_spatialDiscretisation), intent(in)                   :: rdiscretisationTrial

    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.; test space.
    type(t_spatialDiscretisation), intent(in)                   :: rdiscretisationTest

    ! The bilinear form which is currently being evaluated:
    type(t_bilinearForm), intent(in)                            :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(in)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(in)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(in)  :: Dpoints

    ! An array accepting the DOF`s on all elements in the trial space.
    ! DIMENSION(#local DOF`s in trial space,nelements)
    integer, dimension(:,:), intent(in) :: IdofsTrial

    ! An array accepting the DOF`s on all elements trial in the trial space.
    ! DIMENSION(#local DOF`s in test space,nelements)
    integer, dimension(:,:), intent(in) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It is usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(in)              :: rdomainIntSubset

    ! Optional: A collection structure to provide additional
    ! information to the coefficient routine.
    type(t_collection), intent(inout), optional      :: rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the bilinear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the bilinear form.
    real(DP), dimension(:,:,:), intent(out)                      :: Dcoefficients
  !</output>

  !</subroutine>

    ! local variables
    type(t_vectorBlock), pointer :: p_rsolution
    integer(I32) :: ieltyp
    integer :: ieldistr
    real(DP), dimension(npointsPerElement,nelements) :: Dvalues
    real(DP) :: nu

    ! Equation: 1*u_t + 1*u*u_x - nu*u_xx = 0
    Dcoefficients(1,:,:) = 1.0_DP

    nu = 0.01    ! nu = diffusion coefficient; because of the weak formulation,
                 !      the '-' disappears!
    Dcoefficients(3,:,:) = nu

    nu = 0.00    ! No diffusion in y=t direction
    Dcoefficients(4,:,:) = nu

    ! From the collection, get the current solution vector
    p_rsolution => collct_getvalue_vec (rcollection, 'RX')

    ! In the first subvector of p_rsolution, there is our scalar solution u.
    !
    ! Evaluate the solution in all the points where we need it.

    call fevl_evaluate_sim (p_rsolution%RvectorBlock(1), rdomainIntSubset,&
        DER_FUNC,Dvalues)

    ! Transfer the values to the coefficient vector.
    ! DEBUG: Dcoefficients(2,1:npointsPerElement,1:nelements) = 0.0
    Dcoefficients(2,1:npointsPerElement,1:nelements) = Dvalues(1:npointsPerElement,1:nelements)

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine coeff_RHS (rdiscretisation,rform, &
                  nelements,npointsPerElement,Dpoints, &
                  IdofsTest,rdomainIntSubset,&
                  Dcoefficients,rcollection)

    use basicgeometry
    use triangulation
    use collection
    use scalarpde
    use domainintegration

  !<description>
    ! This subroutine is called during the vector assembly. It has to compute
    ! the coefficients in front of the terms of the linear form.
    !
    ! The routine accepts a set of elements and a set of points on these
    ! elements (cubature points) in real coordinates.
    ! According to the terms in the linear form, the routine has to compute
    ! simultaneously for all these points and all the terms in the linear form
    ! the corresponding coefficients in front of the terms.
  !</description>

  !<input>
    ! The discretisation structure that defines the basic shape of the
    ! triangulation with references to the underlying triangulation,
    ! analytic boundary boundary description etc.
    type(t_spatialDiscretisation), intent(in)                   :: rdiscretisation

    ! The linear form which is currently to be evaluated:
    type(t_linearForm), intent(in)                              :: rform

    ! Number of elements, where the coefficients must be computed.
    integer, intent(in)                                         :: nelements

    ! Number of points per element, where the coefficients must be computed
    integer, intent(in)                                         :: npointsPerElement

    ! This is an array of all points on all the elements where coefficients
    ! are needed.
    ! Remark: This usually coincides with rdomainSubset%p_DcubPtsReal.
    ! DIMENSION(dimension,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), intent(in)  :: Dpoints

    ! An array accepting the DOF`s on all elements trial in the trial space.
    ! DIMENSION(#local DOF`s in test space,nelements)
    integer, dimension(:,:), intent(in) :: IdofsTest

    ! This is a t_domainIntSubset structure specifying more detailed information
    ! about the element set that is currently being integrated.
    ! It is usually used in more complex situations (e.g. nonlinear matrices).
    type(t_domainIntSubset), intent(in)              :: rdomainIntSubset

    ! Optional: A collection structure to provide additional
    ! information to the coefficient routine.
    type(t_collection), intent(inout), optional      :: rcollection

  !</input>

  !<output>
    ! A list of all coefficients in front of all terms in the linear form -
    ! for all given points on all given elements.
    !   DIMENSION(itermCount,npointsPerElement,nelements)
    ! with itermCount the number of terms in the linear form.
    real(DP), dimension(:,:,:), intent(out)                      :: Dcoefficients
  !</output>

  !</subroutine>

    Dcoefficients = 0.0_DP

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine getBoundaryValues (Icomponents,rdiscretisation,rboundaryRegion,ielement, &
                                   cinfoNeeded,iwhere,dwhere, Dvalues, rcollection)

  use collection
  use spatialdiscretisation
  use discretebc

!<description>
  ! This subroutine is called during the discretisation of boundary
  ! conditions. It calculates a special quantity on the boundary, which is
  ! then used by the discretisation routines to generate a discrete
  ! 'snapshot' of the (actually analytic) boundary conditions.
!</description>

!<input>
  ! Component specifier.
  ! For Dirichlet boundary:
  !   Icomponents(1) defines the number of the boundary component, the value
  !   should be calculated for (e.g. 1=1st solution component, e.g. X-velocitry,
  !   2=2nd solution component, e.g. Y-velocity,...)
  integer, dimension(:), intent(in)                           :: Icomponents

  ! The discretisation structure that defines the basic shape of the
  ! triangulation with references to the underlying triangulation,
  ! analytic boundary boundary description etc.
  type(t_spatialDiscretisation), intent(in)                   :: rdiscretisation

  ! Boundary region that is currently being processed.
  type(t_boundaryRegion), intent(in)                          :: rboundaryRegion

  ! The element number on the boundary which is currently being processed
  integer, intent(in)                                         :: ielement

  ! The type of information, the routine should calculate. One of the
  ! DISCBC_NEEDxxxx constants. Depending on the constant, the routine has
  ! to return one or multiple information value in the result array.
  integer, intent(in)                                         :: cinfoNeeded

  ! A reference to a geometric object where information should be computed.
  ! cinfoNeeded=DISCBC_NEEDFUNC :
  !   iwhere = number of the point in the triangulation or
  !          = 0, if only the parameter value of the point is known; this
  !               can be found in dwhere,
  ! cinfoNeeded=DISCBC_NEEDDERIV :
  !   iwhere = number of the point in the triangulation or
  !          = 0, if only the parameter value of the point is known; this
  !               can be found in dwhere,
  ! cinfoNeeded=DISCBC_NEEDINTMEAN :
  !   iwhere = number of the edge where the value integral mean value
  !            should be computed
  integer, intent(in)                                          :: iwhere

  ! A reference to a geometric object where information should be computed.
  ! cinfoNeeded=DISCBC_NEEDFUNC :
  !   dwhere = parameter value of the point where the value should be computed,
  ! cinfoNeeded=DISCBC_NEEDDERIV :
  !   dwhere = parameter value of the point where the value should be computed,
  ! cinfoNeeded=DISCBC_NEEDINTMEAN :
  !   dwhere = 0 (not used)
  real(DP), intent(in)                                        :: dwhere

  ! Optional: A collection structure to provide additional
  ! information to the coefficient routine.
  type(t_collection), intent(inout), optional                 :: rcollection

!</input>

!<output>
  ! This array receives the calculated information. If the caller
  ! only needs one value, the computed quantity is put into Dvalues(1).
  ! If multiple values are needed, they are collected here (e.g. for
  ! DISCBC_NEEDDERIV: Dvalues(1)=x-derivative, Dvalues(2)=y-derivative,...)
  real(DP), dimension(:), intent(out)                         :: Dvalues
!</output>

!</subroutine>

    ! To get the X/Y-coordinates of the boundary point, use:
    !
    ! REAL(DP) :: dx,dy
    !
    ! CALL boundary_getCoords(rdiscretisation%p_rboundary, &
    !     rboundaryRegion%iboundCompIdx, dwhere, dx, dy)

    ! In this problem, we are definitely only asked for point values
    ! on the boundary.
    if (cinfoNeeded .ne. DISCBC_NEEDFUNC) then
      print *,'getBoundaryValues: Only point values implemented!'
      stop
    end if

    ! Return zero Dirichlet boundary values for all situations.
    select case (rboundaryRegion%iboundSegIdx)
    case (1)
      ! Initial conditions on segment 1: sin(Pi*x).
      ! dwhere=x here since we are on the unit square on the first segment,
      ! i.e. the parameter value coincides with the x-coordinate. Otherwise,
      ! we would have to ask the analytic boundary for the forrect x-position.
      Dvalues(1) = sin(dwhere * SYS_PI)
    case (2,3,4)
      ! Boundary values on segment 2,4 is =0.
      ! Boundary segment 3 should not occur.
      Dvalues(1) = 0.0_DP
    ! elsewhere: does not occur.
    end select

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      PROGRAM FILTERFIX
C-------------------------
C Fix up HST filter discriptions
C Read in wavelength and throughput
C Interpolate between the bins to give even bin sizes
C--------------------------
C
      IMPLICIT NONE
C
      INTEGER NFILT,IFILT
      PARAMETER (NFILT=8)
c
      CHARACTER*20 FNAME(NFILT),FNAMEI(NFILT)
      DATA FNAME /
     & 'F110W_NIC2.dat',
     & 'F160W_NIC2.dat',
     & 'F606W_ACS.dat',
     & 'F675W_WFPC2.dat',
     & 'F775W_ACS.dat',
     & 'F814W_WFPC2.dat',
     & 'F850LP_ACS.dat',
     & 'F850lp_WFPC2.dat'/

      INTEGER NLAM,ILAM,JLAM,KLAM,NLAMI
      PARAMETER (NLAM=15000)

      REAL LAM(NLAM),THROUG(NLAM)
      REAL LAMX,LAMN,LAMS
C
      REAL LAMI(NLAM),THROUGI(NLAM)
C
      DO IFILT = 1,NFILT
        OPEN (UNIT=30,FILE=FNAME(IFILT),FORM='FORMATTED')
        ILAM = 0
 100    CONTINUE
          ILAM = ILAM + 1
          READ (30,'(11X,F9.3,9X,F12.10)',END=101) LAM(ILAM),
     &                                             THROUG(ILAM)
        GOTO 100
 101    CONTINUE
        CLOSE(30)

C OK fix up

        JLAM = ILAM -1
        LAMN = LAM(1)
        LAMX = LAM(JLAM)
        LAMS = 5.0
        NLAMI = NINT((LAMX-LAMN)/LAMS) + 1
        DO ILAM = 1,NLAMI
          LAMI(ILAM) = LAMN + FLOAT(ILAM-1)*LAMS
          DO KLAM = 1,JLAM-1
            THROUGI(ILAM) = 0.0
            IF (LAMI(ILAM).GE.LAM(KLAM).AND.
     &          LAMI(ILAM).LT.LAM(KLAM+1)) THEN
              THROUGI(ILAM) = THROUG(KLAM) +
     &                        ((LAMI(ILAM)-LAM(KLAM))/
     &                         (LAM(KLAM+1)-LAM(KLAM)))*
     &                        (THROUG(KLAM+1) - THROUG(KLAM))
              GOTO 200
            ENDIF
          ENDDO
 200      CONTINUE
        ENDDO

C Write out the fixed up

        FNAMEI(IFILT) = 'new'//FNAME(IFILT)
        OPEN (UNIT=30,FILE=FNAMEI(IFILT),FORM='FORMATTED')
        DO ILAM = 1,NLAMI
          WRITE(30,'(5X,F8.2,5X,F8.6)') LAMI(ILAM),THROUGI(ILAM)
        ENDDO
        CLOSE(30)

      ENDDO

      END

max_stars_repo_name: temuller/PISCoLA

max_issues_repo_name: temuller/PISCoLA

max_forks_repo_name: temuller/PISCoLA


text:
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+//
!                                                                                //
! mc.f                                                                                //
!                                                                                //
! D. C. Groothuizen Dijkema - January, 2020                                                                                //
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+//

! Approximating mathematical constants with Monte Carlo methods

module monte_carlo
  implicit none
  ! please ignore the incongruity of having both an analytical definition of pi to use in one function, and another function to calculate
  ! the value of pi itself
  real(kind=8), parameter :: pi=4.*atan(1.)
contains

  function simulate_e() result(cnt)
    !
    ! Count the number of uniformly distributed random real numbers on [0,1] which are needed for their sum to become strictly greater than
    ! one.
    ! Euler's number is the expected value of this function.
    !
    ! returns
    ! -------
    ! cnt : integer
    !   - The count of random numbers.
    !

    integer :: cnt
    real(kind=8) :: sum,rand_num

    cnt=0
    sum=0.
    ! produce the count
    do while (sum.LE.1)
      call random_number(rand_num)
      sum=sum+rand_num
      cnt=cnt+1
    end do
  end function simulate_e

  function simulate_pi() result(within)
    !
    ! Randomly generate a point in the unit square and determine if it falls within the unit circle.
    !
    ! returns
    ! -------
    ! within : integer
    !   - 1 if the randomly generated point is within the unit circle, 0 otherwise. This should be a logical, but using an integer allows
    !     simple addition.
    !
    integer :: within

    real(kind=8) :: distance
    real(kind=8), dimension(2) :: point

    ! determine the point and its distance to the origin
    call random_number(point)
    distance=point(1)**2+point(2)**2

    if (distance.LE.1) then
      within=1
    else
      within=0
    end if
  end function simulate_pi

  function simulate_chord_length() result(distance)
    !
    ! Randomly generate a point on the radius of the unit circle, and determine its distance to the point (1,0). This is,
    ! effectively, the distance between any two random points, as those points could be rotated until one of them was (1,0)
    !
    ! returns
    ! -------
    ! distance : real(kind=8)
    !   - The chord length of two random points on the unit circle.
    !
    real(kind=8) :: distance,angle

    call random_number(angle)
    ! the range [0,1] to [0,2*pi]
    angle=angle*2.*pi
    distance=2.*sin(0.5*angle) ! formula derived simply from the definition of the sine function
  end function simulate_chord_length

  subroutine approximate_e(n_itr,means)
    !
    ! Use Monte Carlo simulations to produce successfully more accurate approximations of Euler's number.
    !
    ! parameters
    ! ----------
    ! n_itr : integer, intent(in)
    !   - The number of simulations to run.
    ! means : real(kind=8), intent(inout), dimension(n_itr)
    !   - The array in which to store the moving mean of the simulations.
    !
    integer, intent(in) :: n_itr
    real(kind=8), intent(inout), dimension(n_itr) :: means

    real(kind=8) :: mean
    integer :: itr,approx

    ! initialise the mean
    mean=simulate_e()
    means(1)=mean
    ! produce a sample, update the moving mean, and output
    do itr=2,n_itr
      approx=simulate_e()
      mean=mean+(approx-mean)/(itr+1)
      means(itr)=mean
    end do
  end subroutine approximate_e

  subroutine approximate_pi_integration(n_itr,means)
    !
    ! Use Monte Carlo simulations to produce successfully more accurate approximations of pi.
    !
    ! parameters
    ! ----------
    ! n_itr : integer, intent(in)
    !   - The number of simulations to run.
    ! means : real(kind=8), intent(inout), dimension(n_itr)
    !   - The array in which to store the moving mean of the simulations.
    !
    integer, intent(in) :: n_itr
    real(kind=8), intent(inout), dimension(n_itr) :: means

    integer :: itr,cnt

    ! initialise
    cnt=0
    ! produce a sample, update the moving mean, and output
    do itr=1,n_itr
      cnt=cnt+simulate_pi()
      means(itr)=4*real(cnt)/itr
    end do
  end subroutine approximate_pi_integration

  subroutine approximate_pi_chord_length(n_itr,means)
    !
    ! Use Monte Carlo simulations to produce successfully more accurate approximations of pi.
    !
    ! parameters
    ! ----------
    ! n_itr : integer, intent(in)
    !   - The number of simulations to run.
    ! means : real(kind=8), intent(inout), dimension(n_itr)
    !   - The array in which to store the moving mean of the simulations.
    !
    integer, intent(in) :: n_itr
    real(kind=8), intent(inout), dimension(n_itr) :: means

    real(kind=8) :: mean,approx
    integer :: itr

    ! initialise
    mean=4/simulate_chord_length()
    means(1)=mean
    ! produce a sample, update the moving mean, and output
    do itr=2,n_itr
      approx=simulate_chord_length()
      mean=mean+(approx-mean)/(itr+1)
      means(itr)=4/mean
    end do
  end subroutine approximate_pi_chord_length

end module monte_carlo


text:
c
c     This is a test for MPI_ALLGATHERV for a maximum of 1024 processes.
c
      program testmpi
      integer           mnprocs, lcwk1
      double precision  accurc
      parameter         ( mnprocs = 1024, lcwk1 = mnprocs*3 )
      parameter         ( accurc = 1d-9 )
      integer           comm, rc, myid, nprocs, ierr, i,
     &     recvts(0:mnprocs-1), displs(0:mnprocs-1)
      double precision  wrkbuf(3), cwk1(lcwk1)
      double precision  lobrdr, hibrdr, elem
      logical           error, gerr
      include           'mpif.h'
c
      call MPI_INIT( ierr )
      comm = MPI_COMM_WORLD
      call MPI_COMM_RANK( comm, myid, ierr )
      call MPI_COMM_SIZE( comm, nprocs, ierr )
c
c     check that the maximum number of processes is not exceeded
c
      if (nprocs.gt.mnprocs) then
         print *, 'Can only be run with a maximum of', mnprocs,
     &        ' processes!'
         call MPI_ABORT( comm, 1 )
      end if
c
c     initialize receive buffer with -10
c
      do i = 1, lcwk1
         cwk1(i) = -10
      end do
c
c     initialize send buffer with rank in MPI_COMM_WORLD
c
      do i=1,3
         wrkbuf(i) = myid
      end do
c
c     initialize receive count and displacements
c
      do i = 0, nprocs-1
         recvts(i) = 3
         displs(i) = 3 * i
      end do
c
      call MPI_ALLGATHERV( wrkbuf, recvts(myid), MPI_DOUBLE_PRECISION,
     &     cwk1, recvts, displs, MPI_DOUBLE_PRECISION,
     &     comm, ierr )
c
c     check gathered data
c
      error = .false.
      do i = 1, nprocs
         lobrdr = DBLE(i-1) - accurc
         hibrdr = DBLE(i-1) + accurc
         do j = 1, 3
c     check that received data is equal to expected data within
c     a reasonable level of accuracy
            elem = cwk1((i-1)*3 + j)
            if (elem.lt.lobrdr) then
               error = .true.
            else if (elem.gt.hibrdr) then
               error = .true.
            end if
         end do
      end do
c
c     build global error indicator via MPI_REDUCE
c
      call MPI_REDUCE( error, gerr, 1, MPI_LOGICAL, MPI_LOR, 0, comm,
     &     ierr )
c
c     print message on process 0
c
      if( myid.eq.0 ) then
         if ( gerr ) then
            print *, 'Some errors have occurred during MPI_ALLGATHERV'
         else
            print *, ' No Errors'
         end if
      end if

      call MPI_FINALIZE(rc)
c
      end
c

max_stars_repo_name: RWTH-OS/MP-MPICH

max_issues_repo_name: RWTH-OS/MP-MPICH

max_forks_repo_name: RWTH-OS/MP-MPICH


!
  subroutine make_grid_rect(tetra,verts_rphiz,grid_size,Rmin,Rmax,Zmin,Zmax)
!
  use tetra_grid_mod, only: ntetr,nvert,tetrahedron_grid
  use constants, only: pi
!
  implicit none
!
  type(tetrahedron_grid), dimension(ntetr), intent(out) :: tetra
  double precision,  dimension(3,nvert), intent(out) :: verts_rphiz
  integer, dimension(3), intent(in) :: grid_size
  double precision, intent(in) :: Rmin,Rmax,Zmin,Zmax
  integer :: ir,iphi,iz,iv,itetr,i,j,k,l
  integer :: ind_tetr,ind_tetr1,ind_tetr2,iface1,iface2
  integer :: nvertinner,nper
  integer, dimension(4) :: knots1,knots2
  integer, dimension(3,4,3) :: ivert_prism
  double precision :: hr,hphi,hz,r,phi,z,x,y,s
  double precision :: R_c,Z_c
  double precision, dimension(4)   :: sp,pp,zp
  integer,          dimension(:,:,:),       allocatable :: inodes,itetrbeg
  integer :: nr,nphi,nz
!
  nr = grid_size(1)
  nphi = grid_size(2)
  nz = grid_size(3)
!
! vertices of the first tetrahedron:
! point 1:
  ivert_prism(1,1,1)=0
  ivert_prism(2,1,1)=0
  ivert_prism(3,1,1)=0
! point 2:
  ivert_prism(1,2,1)=1
  ivert_prism(2,2,1)=0
  ivert_prism(3,2,1)=0
! point 3:
  ivert_prism(1,3,1)=0
  ivert_prism(2,3,1)=1
  ivert_prism(3,3,1)=0
! point 4:
  ivert_prism(1,4,1)=0
  ivert_prism(2,4,1)=0
  ivert_prism(3,4,1)=1
!
! vertices of the second tetrahedron:
! point 1:
  ivert_prism(:,1,2)=ivert_prism(:,2,1)
! point 2:
  ivert_prism(:,2,2)=ivert_prism(:,3,1)
! point 3:
  ivert_prism(:,3,2)=ivert_prism(:,4,1)
! point 4:
  ivert_prism(1,4,2)=1
  ivert_prism(2,4,2)=0
  ivert_prism(3,4,2)=1
!
! vertices of the third tetrahedron:
! point 1:
  ivert_prism(:,1,3)=ivert_prism(:,3,1)
! point 2:
  ivert_prism(:,2,3)=ivert_prism(:,4,1)
! point 3:
  ivert_prism(:,3,3)=ivert_prism(:,4,2)
! point 4:
  ivert_prism(1,4,3)=0
  ivert_prism(2,4,3)=1
  ivert_prism(3,4,3)=1
!
  hr=(Rmax-Rmin)/nr
  hphi=2.d0*pi/nphi
  hz=(Zmax-Zmin)/nz
!
  R_c=0.5d0*(Rmax+Rmin)
  Z_c=0.5d0*(Zmax+Zmin)
  nper=0
!
!
  allocate(inodes(0:nr,0:nz,0:nphi))
!
  iv=0
!
  do iphi=0,nphi
    phi=hphi*iphi
    do ir=0,nr
      do iz=0,nz
        r=rmin+hr*ir
        z=zmin+hz*iz
!
        x=(r-R_c)
        y=(z-Z_c)
        r=R_c+x*cos(nper*phi)+y*sin(nper*phi)
        z=Z_c-x*sin(nper*phi)+y*cos(nper*phi)! +hz/3000
        s=r
!
        iv=iv+1
        inodes(ir,iz,iphi)=iv
!
        verts_rphiz(:,iv) = [s,phi,z]
!
      enddo
    enddo
  enddo
!
  allocate(itetrbeg(nr,nz,nphi))
!
  ind_tetr=0
!
  do iphi=1,nphi
    do ir=1,nr
      do iz=1,nz
!
        itetrbeg(ir,iz,iphi)=ind_tetr
!
        do itetr=1,3
          ind_tetr=ind_tetr+1
          do i=1,4
            iv=inodes(ir-1+ivert_prism(1,i,itetr),iz-1+ivert_prism(2,i,itetr),iphi-1+ivert_prism(3,i,itetr))
            tetra(ind_tetr)%ind_knot(i)=iv
            tetra(ind_tetr)%neighbour_tetr(i)=-1
            tetra(ind_tetr)%neighbour_face(i)=-1
            tetra(ind_tetr)%neighbour_perbou_phi(i)=0
          enddo
        enddo
!
        do itetr=1,3
          ind_tetr=ind_tetr+1
          do i=1,4
            iv=inodes(ir-ivert_prism(1,i,itetr),iz-ivert_prism(2,i,itetr),iphi-ivert_prism(3,i,itetr))
            tetra(ind_tetr)%ind_knot(i)=iv
            tetra(ind_tetr)%neighbour_tetr(i)=-1
            tetra(ind_tetr)%neighbour_face(i)=-1
            tetra(ind_tetr)%neighbour_perbou_phi(i)=0
          enddo
        enddo
!
        do i=1,6
          ind_tetr1=itetrbeg(ir,iz,iphi)+i
          do j=1,6
            if(j.eq.i) cycle
            ind_tetr2=itetrbeg(ir,iz,iphi)+j
            call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                 tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
            if(iface1.ne.-1) then
              tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
              tetra(ind_tetr1)%neighbour_face(iface1)=iface2
            endif
          enddo
        enddo
      enddo
    enddo
  enddo
!
!
  do iphi=1,nphi
    do ir=1,nr
      do iz=1,nz
!
        if(ir.gt.1) then
          do i=1,3
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=4,6
              ind_tetr2=itetrbeg(ir-1,iz,iphi)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
        if(ir.lt.nr) then
          do i=4,6
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=1,3
              ind_tetr2=itetrbeg(ir+1,iz,iphi)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
        if(iz.gt.1) then
          do i=1,3
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=4,6
              ind_tetr2=itetrbeg(ir,iz-1,iphi)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
        if(iz.lt.nz) then
          do i=4,6
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=1,3
              ind_tetr2=itetrbeg(ir,iz+1,iphi)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
        if(iphi.gt.1) then
          do i=1,3
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=1,3
              ind_tetr2=itetrbeg(ir,iz,iphi-1)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
          do i=4,6
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=4,6
              ind_tetr2=itetrbeg(ir,iz,iphi-1)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
        if(iphi.lt.nphi) then
          do i=1,3
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=1,3
              ind_tetr2=itetrbeg(ir,iz,iphi+1)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
          do i=4,6
            ind_tetr1=itetrbeg(ir,iz,iphi)+i
            do j=4,6
              ind_tetr2=itetrbeg(ir,iz,iphi+1)+j
              call check_neighbour(tetra(ind_tetr1)%ind_knot(:),              &
                                   tetra(ind_tetr2)%ind_knot(:),iface1,iface2)
              if(iface1.ne.-1) then
                tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
                tetra(ind_tetr1)%neighbour_face(iface1)=iface2
              endif
            enddo
          enddo
        endif
!
      enddo
    enddo
  enddo
!
! Neighbours through the periodic boundary
!
  nvertinner=(nr+1)*(nz+1)*nphi
!
  do ir=1,nr
    do iz=1,nz
!
      do i=1,3
        ind_tetr1=itetrbeg(ir,iz,1)+i
        knots1(:)=tetra(ind_tetr1)%ind_knot(:)
        do j=1,3
          ind_tetr2=itetrbeg(ir,iz,nphi)+j
          knots2(:)=tetra(ind_tetr2)%ind_knot(:)
          knots2(:)=modulo(knots2(:),nvertinner)
          call check_neighbour(knots1,knots2,iface1,iface2)
          if(iface1.ne.-1) then
            tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
            tetra(ind_tetr1)%neighbour_face(iface1)=iface2
            tetra(ind_tetr1)%neighbour_perbou_phi(iface1)=-1
            tetra(ind_tetr2)%neighbour_tetr(iface2)=ind_tetr1
            tetra(ind_tetr2)%neighbour_face(iface2)=iface1
            tetra(ind_tetr2)%neighbour_perbou_phi(iface2)=1
          endif
        enddo
      enddo
!
      do i=4,6
        ind_tetr1=itetrbeg(ir,iz,1)+i
        knots1(:)=tetra(ind_tetr1)%ind_knot(:)
        do j=4,6
          ind_tetr2=itetrbeg(ir,iz,nphi)+j
          knots2(:)=tetra(ind_tetr2)%ind_knot(:)
          knots2(:)=modulo(knots2(:),nvertinner)
          call check_neighbour(knots1,knots2,iface1,iface2)
          if(iface1.ne.-1) then
            tetra(ind_tetr1)%neighbour_tetr(iface1)=ind_tetr2
            tetra(ind_tetr1)%neighbour_face(iface1)=iface2
            tetra(ind_tetr1)%neighbour_perbou_phi(iface1)=-1
            tetra(ind_tetr2)%neighbour_tetr(iface2)=ind_tetr1
            tetra(ind_tetr2)%neighbour_face(iface2)=iface1
            tetra(ind_tetr2)%neighbour_perbou_phi(iface2)=1
          endif
        enddo
      enddo
!
    enddo
  enddo
!
! Check the number of tetrahedrons at the XY boundary
  j=0
!
b: do ind_tetr=1,ntetr
  do i=1,4
     if(tetra(ind_tetr)%neighbour_tetr(i).lt.1) then
       j=j+1
       cycle b
     endif
   enddo
  enddo b
!
  deallocate(inodes,itetrbeg)
!
  end subroutine make_grid_rect
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
  subroutine check_neighbour(knots1,knots2,iface1,iface2)
!
  implicit none
!
  logical :: result
  integer :: iface1,iface2,i,j,matches
  integer, dimension(4) :: knots1,knots2
  integer, dimension(3) :: match1,match2
  logical, dimension(4) :: unused
!
  iface1=-1
  iface2=-1
  matches=0
!
  a: do i=1,4
    do j=1,4
      if(knots1(i).eq.knots2(j)) then
        matches=matches+1
        match1(matches)=i
        match2(matches)=j
        if(matches.eq.3) exit a
      endif
    enddo
    if(i-matches.gt.1) return
  enddo a
!
  unused= .true.
  do i=1,3
    unused(match1(i))= .false.
  enddo
  do i=1,4
    if(unused(i)) then
      iface1=i
      exit
    endif
  enddo
!
  unused= .true.
  do i=1,3
    unused(match2(i))= .false.
  enddo
  do i=1,4
    if(unused(i)) then
      iface2=i
      exit
    endif
  enddo
!
  end subroutine check_neighbour
!
!ccccccccccccccccccccccccccccccccccccccccccc
!