text:
c
c***************************************************************************
c
      subroutine sigmad_cql3d(i, j, n, m, rho, rho_a,
     .   gradprlb, bmod, bmod0,
     .   xm, q, xn, xnuomg,
     .   xkt, omgc, omgp2,
     .   lmin, lmax, nzfun, ibessel,
     .   xkxsav, xkysav, nphi, capr,
     .   bx,    by,    bz,
     .   uxx, uxy, uxz,
     .   uyx, uyy, uyz,
     .   uzx, uzy, uzz,
     .   sigxx, sigxy, sigxz,
     .   sigyx, sigyy, sigyz,
     .   sigzx, sigzy, sigzz,
     .   delta0, ndist, nupar, nuper, n_psi,
     .   n_psi_dim,
     .   UminPara, UmaxPara, UPERP, UPARA,
     .   vc_mks, df_cql_uprp, df_cql_uprl, nbessj, ncount)

*     ---------------------------------------------------------
*     This routine uses the modified Z functions Z0, Z1, Z2
*     with the appropriate sign changes for k_parallel < 0.0
*     No rotation is made.  Result is in the Stix frame.
*     ---------------------------------------------------------

      implicit none

      integer lmin, lmax, nzfun, lmaxdim, l, labs, ibessel,
     .    i, j, n, m, nphi, ndist, myid, iflag, nproc, ni, mi

      real xkperp, xkprl, xm, q, xn, xkt, omgc, omgp2, xme
      real xkprl_eff, fgam, y0, sgn_kprl, reson, dzetal
      real dakbdkb, xnuomg, gradprlb, bmod, bmod0, nu_coll, gamma_coll
      real akprl, gammab, rho, alpha, eps0, omgrf, v0i, emax
      real a, b, xnurf, pi, delta0, rhol
      real bx, by, bz, bratio, denom
      real time, t1, dummy, tmin, second1

      real xkxsav, xkysav, capr, argd
      real xkphi
      real xkalp, xkbet, xk0, rgamma

      complex zi, zfunct, fzeta, omgrfc

      complex zfunct0, zeta0, sig3cold, z0, z1, z2, dz0, dz1, dz2

      complex sig0, sig1, sig2, sig3, sig4, sig5
      complex sig0_a, sig1_a, sig2_a, sig3_a, sig4_a, sig5_a
      complex sig0_h, sig1_h, sig2_h, sig3_h, sig4_h, sig5_h

      complex sig0l, sig1l, sig2l, sig3l, sig4l, sig5l

      complex sigxx, sigxy, sigxz,
     1        sigyx, sigyy, sigyz,
     1        sigzx, sigzy, sigzz

      real uxx, uxy, uxz,
     1     uyx, uyy, uyz,
     2     uzx, uzy, uzz

      parameter (lmaxdim = 99)

      complex xil(0: lmaxdim), xilp(0: lmaxdim)
      complex exil(0: lmaxdim), exilp(0: lmaxdim),
     .                    exilovergam(0: lmaxdim)

      complex zetal, zieps0, arg,
     .   al, bl, cl,
     .   gamma, zeta_eff

      integer  :: n_psi_dim
      integer :: nuper, nupar, n_psi

c      integer, parameter :: NBESSJ = 7
      integer NBESSJ, ncount, nc

      real :: UPERP(NUPER),UPARA(NUPAR)
      real :: UPERP0, UPARA0, dfdupar0, dfduper0

      real :: DFDUPER(NUPER,NUPAR),DFDUPAR(NUPER,NUPAR)
      real :: W,K1(3),K2(3),KPER1,KPER2,ENORM,ZSPEC,ASPEC,BMAG,DensSPEC
      integer :: NSBESSJ,IFAIL
      COMPLEX WSPEC(3,3)
      complex :: factor
	real :: UminPara,UmaxPara
      real :: XI1(NUPER),JNXI1(NUPER, NBESSJ)
      real :: XI2(NUPER),JNXI2(NUPER, NBESSJ)

      real :: df_cql_uprp(NUPER, NUPAR, n_psi_dim)
      real :: df_cql_uprl(NUPER, NUPAR, n_psi_dim)
      real :: rho_a(n_psi_dim), vc_mks

      integer :: i_uperp, i_upara, i_psi

      common/sigcom/zi, eps0, v0i, omgrf, xk0

      nu_coll =  .01 * omgrf

      xme = 9.11e-31
      zieps0 = zi * eps0
      alpha = sqrt(2. * xkt / xm)
      rhol = alpha / omgc
      xkphi = nphi / capr
      omgrfc = omgrf * (1. + zi * xnuomg)

      xkalp = uxx * xkxsav + uxy * xkysav + uxz * xkphi
      xkbet = uyx * xkxsav + uyy * xkysav + uyz * xkphi
      xkprl = uzx * xkxsav + uzy * xkysav + uzz * xkphi

      xkperp = sqrt(xkalp**2 + xkbet**2)

      akprl = abs(xkprl)

      sgn_kprl = sign(1.0, xkprl)

      if (akprl .lt. 0.01) then

         xkprl = 0.01 * sgn_kprl
         akprl = abs(xkprl)

      end if

*     -----------------------
*     Maxwellian distribution
*     -----------------------
      if(ndist .eq. 0 .or. ndist .eq. 2)then

         gamma = 0.5 * xkperp**2 * rhol**2
         rgamma = real(gamma)

         if(rgamma .ge. 1.0e-08)
     .      call besiexp(gamma, lmax, exil, exilp, lmaxdim, exilovergam)

         if(rgamma .lt. 1.0e-08)
     .      call bes_expand(gamma, lmax, exil, exilp, lmaxdim,
     .                                                      exilovergam)

         sig0 = 0.0
         sig1 = 0.0
         sig2 = 0.0
         sig3 = 0.0
         sig4 = 0.0
         sig5 = 0.0

         do l = lmin, lmax
            labs = abs(l)

            reson = (omgrf - l * real(omgc)) / omgrf
            if (abs(reson) .lt. 0.02)then
               zetal = (omgrfc - l * omgc) / (xkprl * alpha)
               dzetal = omgrf * xnuomg / (xkprl * alpha)
            else
               zetal = (omgrf  - l * omgc) / (xkprl * alpha)
               dzetal = 0.0
            end if

            gammab = l * omgc / (2.0 * alpha * xkprl**2)
     .                                            * gradprlb / bmod
            gamma_coll = nu_coll / (akprl * alpha)

            if(xm .eq. xme)gammab = 0.0
            if(abs(gammab) .gt. 1000.0) gammab = 1000.0
            if(abs(gammab) .lt. .01)gammab = .01

            if (nzfun .eq. 0) call z_approx(sgn_kprl, zetal, 0.0,
     .                                                     z0, z1, z2)
            if (nzfun .eq. 1) call z_approx(sgn_kprl, zetal, gammab,
     .                                                     z0, z1, z2)
            if (nzfun .eq. 2) call z_smithe(sgn_kprl, zetal, gammab,
     .                                                     z0, z1, z2)
            if (nzfun .eq. 3) call z_table (sgn_kprl, zetal, gammab,
     .                                         gamma_coll, z0, z1, z2)

            al = 1.0 / (xkprl * alpha) * z0
            bl = 1.0 / (xkprl * alpha) * z1
            cl = 1.0 / (xkprl * alpha) * z2

            sig0l = - zieps0 * omgp2 * rhol**2
     .                                 * (exil(labs) - exilp(labs)) * al
            sig1l = - zieps0 * omgp2 * l**2 * exilovergam(labs) * al
            sig2l = - eps0 * omgp2 * l * (exil(labs) - exilp(labs)) * al
            sig3l = - zieps0 * omgp2 * 2.0 * exil(labs) * cl
            sig4l = - zieps0 * omgp2 * rhol * l * exilovergam(labs) * bl
            sig5l = - eps0 * omgp2 * rhol
     .                                 * (exil(labs) - exilp(labs)) * bl

            sig0 = sig0 + sig0l
            sig1 = sig1 + sig1l
            sig2 = sig2 + sig2l
            sig3 = sig3 + sig3l
            sig4 = sig4 + sig4l
            sig5 = sig5 + sig5l

         end do

         sig1 = sig1 + delta0 * eps0 * omgrf * xkperp**2 / xk0**2
         sig3 = sig3 + delta0 * eps0 * omgrf * xkperp**2 / xk0**2

         sig0_a = zi * aimag(sig0)
         sig1_a = zi * aimag(sig1)
         sig2_a =       real(sig2)
         sig3_a = zi * aimag(sig3)
         sig4_a = zi * aimag(sig4)
         sig5_a =       real(sig5)

      end if

*     -----------------------------------
*     Funky METS calls for non Maxwellian:
*     -----------------------------------
      if (ndist .eq. 1 .or. ndist .eq. 2) then

         Emax = 0.5 * xm * vc_mks**2
         Enorm = Emax / 1.6e-19

         W = omgrf
         K1(1) = xkperp
         K1(2) = 0.0
         K1(3) = xkprl
         K2(1) = K1(1)
         K2(2) = K1(2)
         K2(3) = K1(3)
         KPER1 = xkperp
         KPER2 = xkperp

         ZSPEC = q / 1.6e-19
         ASPEC = xm / 1.67e-27
         BMAG = omgc * (xm / q)
         NSBESSJ = 2 * NBESSJ + 8
c         NSBESSJ = NBESSJ
         DensSPEC = xn
         bratio = bmod0 / bmod

!         call dummy_dist(NUPAR, NUPER,
!     &                 UminPara, UmaxPara,
!     &                 UPERP, UPARA, DFDUPER, DFDUPAR)

!        ------------------------------------------------
!        get CQL3D distribution function on the midplane
!        ------------------------------------------------
         call cql3d_dist(nupar, nuper, n_psi,
     .                 n_psi_dim, rho_a, rho,
     .                 UminPara,UmaxPara,
     .                 df_cql_uprp, df_cql_uprl,
     .                 UPERP, UPARA, DFDUPER, DFDUPAR)

!        go to 1000

!        ------------------------------------------------
!        map CQL3D distribution function off the midplane
!        ------------------------------------------------
         if(bratio .gt. 0.0)then
            do ni = 1, nuper
               do mi = 1, nupar

                  dfdupar0 = dfdupar(ni, mi)
                  dfduper0 = dfduper(ni, mi)

                  argd = uperp(ni)**2 * (1. - bratio) + upara(mi)**2
                  if (argd .le. 0.0) argd = 1.0e-06

                  upara0  = sign(1.0, upara(mi)) * sqrt(argd)

                  dfduper(ni, mi) = dfduper0 * sqrt(bratio)
     .               + dfdupar0 * uperp(ni) / upara0 * (1.0 - bratio)

                  dfdupar(ni, mi) = dfdupar0 * upara(mi) / upara0

               end do
            end do
         end if

 1000    continue

         t1 = second1(dummy)

         do nc = 1, ncount

            call WMATPRECALC_AORSA(ZSPEC,ASPEC,ENORM,BMAG,KPER1,UPERP,
     &                   NUPER,NBESSJ,NSBESSJ,XI1,JNXI1,IFAIL)

         end do

         time = second1(dummy) - t1

  	   tmin = time / 60.
         if (myid.eq.0) then
	      write(6 ,833) time
            write(15,833) time
         endif

  833 format('time taken by WMATPRECALC =',f9.3,4h sec)

         t1 = second1(dummy)

         do nc = 1, ncount

            call GETNONMAXSIGMA_AORSA_NEW(W,
     &                          ZSPEC,ASPEC,DensSPEC,BMAG,
     &                          K1,XI1,JNXI1,
     &                          K1,XI1,JNXI1,NBESSJ,
     &                          Enorm,UminPara,UmaxPara,
     &                          NUPAR,NUPER,UPERP,UPARA,
     &			        DFDUPER,DFDUPAR,
     &                          WSPEC,IFAIL)

         end do

         time = second1(dummy) - t1

  	   tmin = time / 60.
         if (myid.eq.0) then
	      write(6 ,834) time
            write(15,834) time
         endif

  834 format('time taken by GETNONMAXSIGMA =',f9.3,4h sec)

         factor = cmplx(0.,-omgrf * eps0)

         sig1 = WSPEC(1,1) * factor
         sig2 = WSPEC(1,2) * factor
         sig3 = WSPEC(3,3) * factor
         sig4 = 0.0
         sig5 = 0.0
         sig0 = 0.0
         if (xkperp .gt. 0.01) then
            sig4 = WSPEC(3,1) / xkperp * factor
            sig5 = WSPEC(3,2) / xkperp * factor
            sig0 = (WSPEC(2,2) * factor - sig1) / xkperp**2
         endif

         sig0_h =       real(sig0)
         sig1_h =       real(sig1)
         sig2_h = zi * aimag(sig2)
         sig3_h =       real(sig3)
         sig4_h =       real(sig4)
         sig5_h = zi * aimag(sig5)

c         write(6, *) "rho = ", rho
c         write(6, *) "xn = ", xn
c         write(6, *) "xkt = ", xkt
c         write(6, *) "xkperp = ", xkperp
c         write(6, *) "xkprl = ", xkprl
c         write(6, *) "bmod = ", bmod
c         write(6, *) "bmod0 = ", bmod0
c         write(6, *) "bratio =", bratio
c         write(6, *) "xm = ", xm
c         write(6, *) "q = ", q
c         write(6, *) "lmax = ", lmax
c         write(6, *) "omgrf = ", omgrf
c         write(6, *) "omgc = ", omgc
         reson = omgrf / omgc
c         write(6, *) "reson = ", reson
c         write(6, *) "vc_mks = ", vc_mks
c         write(6, *) "DFDUPER(10, 30) = ", DFDUPER(10, 30)
c         write(6, *) "DFDUPAR(10, 30) = ", DFDUPAR(10, 30)
c         write(6, *) "UPERP(10) = ", UPERP(10)
c         write(6, *) "UPARA(30) = ", UPARA(30)
c         write(6, *) "Enorm = ", Enorm
c         write(6, *) "UminPara = ", UminPara
c         write(6, *) "UmaxPara = ", UmaxPara
         write(6, *) "NUPER = ", NUPER
         write(6, *) "NUPAR = ", NUPAR

         write(15, *) "NUPER = ", NUPER
         write(15, *) "NUPAR = ", NUPAR

      end if

!      if(ndist .eq. 2) then

!         sig0 = sig0_a + sig0_h
!         sig1 = sig1_a + sig1_h
!         sig2 = sig2_a + sig2_h
!         sig3 = sig3_a + sig3_h
!         sig4 = sig4_a + sig4_h
!         sig5 = sig5_a + sig5_h

!      end if

         write(6, 1312)sig0
         write(6, 1312)sig1
         write(6, 1312)sig2
         write(6, 1312)sig3
         write(6, 1312)sig4
         write(6, 1312)sig5

         write(15, 1312)sig0
         write(15, 1312)sig1
         write(15, 1312)sig2
         write(15, 1312)sig3
         write(15, 1312)sig4
         write(15, 1312)sig5

*     -----------------------------
*     Swanson's rotation (original):
*     -----------------------------
      sigxx = sig1 + sig0 * xkbet**2
      sigxy = sig2 - sig0 * xkbet * xkalp
      sigxz = sig4 * xkalp + sig5 * xkbet

      sigyx = - sig2 - sig0 * xkbet * xkalp
      sigyy =   sig1 + sig0 * xkalp**2
      sigyz =   sig4 * xkbet - sig5 * xkalp

      sigzx = sig4 * xkalp - sig5 * xkbet
      sigzy = sig4 * xkbet + sig5 * xkalp
      sigzz = sig3

      return

  101 format(i10, 1p8e12.4)
 1314 format(4i10, 1p9e12.4)
 1312 format(1p9e12.4)
  100 format('ier = ', i5, 'besic failed')
      end

c
c***************************************************************************
c

      subroutine cql3d_dist(nupar, nuper, n_psi,
     .                       n_psi_dim, rho_a, rho,
     .                       UminPara,UmaxPara,
     .                       df_cql_uprp, df_cql_uprl,
     .                       UPERP, UPARA, DFDUPER, DFDUPAR)

      implicit none
      integer, intent(in) :: nupar, nuper, n_psi
      integer, intent(in) :: n_psi_dim

      real, intent(out) :: UminPara,UmaxPara
      real, intent(out) :: UPERP(NUPER)
      real, intent(out) :: UPARA(NUPAR)
      real, intent(out) :: DFDUPER(NUPER,NUPAR),
     .                      DFDUPAR(NUPER,NUPAR)

      real, intent(in) :: df_cql_uprp(NUPER, NUPAR, n_psi_dim)
      real, intent(in) :: df_cql_uprl(NUPER, NUPAR, n_psi_dim)
      real :: f, u2, rho, rho_a(n_psi_dim)

      integer n, m, i_psi, i_psin

      i_psi = 1

      do i_psin = 1, n_psi - 1
         if(rho .ge. rho_a(i_psin)     .and.
     .      rho .lt. rho_a(i_psin + 1))
     .      i_psi = i_psin
      end do

      if(rho .ge. rho_a(n_psi)) i_psi = n_psi

c      write(6,  *)
c      write(6,  *) "i_psi = ", i_psi

c      write(15, *)
c      write(15, *) "i_psi = ", i_psi

      if(i_psi .lt. n_psi)then
         do n = 1, nuper
            do m = 1, nupar
               DFDUPER(n, m) = df_cql_uprp(n, m, i_psi)
     .            + (df_cql_uprp(n, m, i_psi + 1)
     .                                     - df_cql_uprp(n, m, i_psi))
     .                    * (rho              - rho_a(i_psi))
     .                    / (rho_a(i_psi + 1) - rho_a(i_psi))
               DFDUPAR(n, m) = df_cql_uprl(n, m, i_psi)
     .            + (df_cql_uprl(n, m, i_psi + 1)
     .                                     - df_cql_uprl(n, m, i_psi))
     .                    * (rho              - rho_a(i_psi))
     .                    / (rho_a(i_psi + 1) - rho_a(i_psi))
            end do
         end do
      end if

      if(i_psi .eq. n_psi)then
         do n = 1, nuper
            do m = 1, nupar
               DFDUPER(n, m) = df_cql_uprp(n, m, i_psi)
               DFDUPAR(n, m) = df_cql_uprl(n, m, i_psi)
            end do
         end do
      end if

  310 format(1p6e12.4)
  311 format(i10, 1p6e12.4)

      return
      end subroutine cql3d_dist

c
c***************************************************************************
c
       subroutine dummy_dist(nupar, nuper,
     &                       UminPara, UmaxPara,
     &                       UPERP, UPARA, DFDUPER, DFDUPAR)
       implicit none

       integer, intent(in) :: nupar, nuper
       real, intent(out) :: UminPara,UmaxPara
       real, intent(out) :: UPERP(NUPER)
       real, intent(out) :: UPARA(NUPAR)
       real, intent(out) :: DFDUPER(NUPER, NUPAR),
     &                      DFDUPAR(NUPER, NUPAR)
       real :: pi,fnorm,f,u2
       integer n,m

       pi = 4.0*atan(1.0)
       fnorm = pi**1.5

       do n = 1, NUPER
          UPERP(n) = 3.5*(real(n-1)/real(NUPER-1))
       end do

       UminPara = -3.5
       UmaxPara =  3.5

       do m = 1, NUPAR
          UPARA(m) = 3.5*(-1.0+2.*(real(m-1)/real(NUPAR-1)))
       end do

       do n = 1, NUPER
          do m = 1, NUPAR
             u2 = UPERP(n)**2 + UPARA(m)**2
             f = exp(-u2) / fnorm
             DFDUPER(n, m) = f *(-2.0 * UPERP(n))
             DFDUPAR(n, m) = f *(-2.0 * UPARA(m))
          end do
       end do

       return
       end subroutine dummy_dist

c
c***************************************************************************
c
      subroutine sigmah_stix(i, j, n, m,
     .   gradprlb, bmod,
     .   xm, q, xn, xnuomg,
     .   xkt, omgc, omgp2,
     .   lmin, lmax, nzfun, ibessel,
     .   xkxsav, xkysav, nphi, capr,
     .   bx,    by,    bz,
     .   uxx, uxy, uxz,
     .   uyx, uyy, uyz,
     .   uzx, uzy, uzz,
     .   sigxx, sigxy, sigxz,
     .   sigyx, sigyy, sigyz,
     .   sigzx, sigzy, sigzz,
     .   delta0, iflag_gammab)

*     ---------------------------------------------------------
*     This routine uses the modified Z functions Z0, Z1, Z2
*     with the appropriate sign changes for k_parallel < 0.0
*     No rotation is made.  Result is in the Stix frame.
*     ---------------------------------------------------------

      implicit none

      integer lmin, lmax, nzfun, lmaxdim, l, labs, ibessel,
     .    i, j, n, m, nphi, iflag_gammab

      real xkperp, xkprl, xm, q, xn, xkt, omgc, omgp2, xme, delta0
      real xkprl_eff, fgam, y0, sgn_kprl, reson, dzetal
      real dakbdkb, xnuomg, gradprlb, bmod, nu_coll, gamma_coll
      real akprl, gammab, rho, alpha, eps0, omgrf, v0i
      real a, b
      real bx, by, bz

      real xkxsav, xkysav, capr
      real xkphi
      real xkalp, xkbet, xk0, rgamma

      complex zi, zfunct, fzeta, omgrfc

      complex zfunct0, zeta0, sig3cold, z0, z1, z2, dz0, dz1, dz2

      complex sig0, sig1, sig2, sig3, sig4, sig5
      complex sig0l, sig1l, sig2l, sig3l, sig4l, sig5l

      complex sigxx, sigxy, sigxz,
     1        sigyx, sigyy, sigyz,
     1        sigzx, sigzy, sigzz

      real uxx, uxy, uxz,
     1     uyx, uyy, uyz,
     2     uzx, uzy, uzz

      parameter (lmaxdim = 99)

      complex xil(0: lmaxdim), xilp(0: lmaxdim)
      complex exil(0: lmaxdim), exilp(0: lmaxdim),
     .                    exilovergam(0: lmaxdim)

      complex zetal, zieps0, arg,
     .   al, bl, cl,
     .   gamma, zeta_eff

      common/sigcom/zi, eps0, v0i, omgrf, xk0

      nu_coll =  .01 * omgrf

      xme = 9.11e-31
      zieps0 = zi * eps0
      alpha = sqrt(2. * xkt / xm)
      rho = alpha / omgc
      xkphi = nphi / capr
      omgrfc = omgrf * (1. + zi * xnuomg)

      xkalp = uxx * xkxsav + uxy * xkysav + uxz * xkphi
      xkbet = uyx * xkxsav + uyy * xkysav + uyz * xkphi
      xkprl = uzx * xkxsav + uzy * xkysav + uzz * xkphi

      xkperp = sqrt(xkalp**2 + xkbet**2)

      akprl = abs(xkprl)

      sgn_kprl = sign(1.0, xkprl)

      if (akprl .lt. 0.01) then

         xkprl = 0.01 * sgn_kprl
         akprl = abs(xkprl)

      end if

      gamma = 0.5 * xkperp**2 * rho**2
      rgamma = real(gamma)

      if(rgamma .ge. 1.0e-08)
     .   call besiexp(gamma, lmax, exil, exilp, lmaxdim, exilovergam)

      if(rgamma .lt. 1.0e-08)
     .   call bes_expand(gamma, lmax, exil, exilp, lmaxdim, exilovergam)

      sig0 = 0.0
      sig1 = 0.0
      sig2 = 0.0
      sig3 = 0.0
      sig4 = 0.0
      sig5 = 0.0

      do l = lmin, lmax
         labs = abs(l)

         reson = (omgrf - l * real(omgc)) / omgrf
         if (abs(reson) .lt. 0.02)then
            zetal = (omgrfc - l * omgc) / (xkprl * alpha)
            dzetal = omgrf * xnuomg / (xkprl * alpha)
         else
            zetal = (omgrf  - l * omgc) / (xkprl * alpha)
            dzetal = 0.0
         end if

         gammab = l * omgc / (2.0 * alpha * xkprl**2)
     .                                            * gradprlb / bmod

         gamma_coll = nu_coll / (akprl * alpha)

         if(xm .eq. xme)gammab = 0.0

         if(abs(gammab) .gt. 1000.0) gammab = 1000.0
         if(abs(gammab) .lt. .01)gammab = .01

         if (nzfun .eq. 0) call z_approx(sgn_kprl, zetal, 0.0,
     .                                                     z0, z1, z2)
         if (nzfun .eq. 1) call z_approx(sgn_kprl, zetal, gammab,
     .                                                     z0, z1, z2)
         if (nzfun .eq. 2) call z_smithe(sgn_kprl, zetal, gammab,
     .                                                     z0, z1, z2)
         if (nzfun .eq. 3) call z_table (sgn_kprl, zetal, gammab,
     .                                         gamma_coll, z0, z1, z2)

         al = 1.0 / (xkprl * alpha) * z0
         bl = 1.0 / (xkprl * alpha) * z1
         cl = 1.0 / (xkprl * alpha) * z2

         sig0l = - zieps0 * omgp2 * rho**2 * (exil(labs) - exilp(labs))
     .           * al
         sig1l = - zieps0 * omgp2 * l**2 * exilovergam(labs) * al
         sig2l = - eps0 * omgp2 * l * (exil(labs) - exilp(labs)) * al
         sig3l = - zieps0 * omgp2 * 2.0 * exil(labs) * cl
         sig4l = - zieps0 * omgp2 * rho * l * exilovergam(labs)  * bl
         sig5l = - eps0 * omgp2 * rho * (exil(labs) - exilp(labs)) * bl

         sig0 = sig0 + sig0l
         sig1 = sig1 + sig1l
         sig2 = sig2 + sig2l
         sig3 = sig3 + sig3l
         sig4 = sig4 + sig4l
         sig5 = sig5 + sig5l

      end do

      sig1 = sig1 + delta0 * eps0 * omgrf * xkperp**2 / xk0**2
      sig3 = sig3 + delta0 * eps0 * omgrf * xkperp**2 / xk0**2

*     -----------------------------
*     Swanson's rotation (original):
*     -----------------------------
      sigxx = sig1 + sig0 * xkbet**2
      sigxy = sig2 - sig0 * xkbet * xkalp
      sigxz = sig4 * xkalp + sig5 * xkbet

      sigyx = - sig2 - sig0 * xkbet * xkalp
      sigyy =   sig1 + sig0 * xkalp**2
      sigyz =   sig4 * xkbet - sig5 * xkalp

      sigzx = sig4 * xkalp - sig5 * xkbet
      sigzy = sig4 * xkbet + sig5 * xkalp
      sigzz = sig3

      return

  101 format(i10, 1p8e12.4)
 1314 format(4i10, 1p9e12.4)
 1312 format(1p9e12.4)
  100 format('ier = ', i5, 'besic failed')
      end

c
c***************************************************************************
c

      subroutine sigmac_stix(i, j, n, m,
     .   xm, q, xn, xnuomg,
     .   xkt, omgc, omgp2,
     .   lmin, lmax, nzfun, ibessel,
     .   xkxsav, xkysav, nphi, capr,
     .   bx,    by,    bz,
     .   uxx, uxy, uxz,
     .   uyx, uyy, uyz,
     .   uzx, uzy, uzz,
     .   sigxx, sigxy, sigxz,
     .   sigyx, sigyy, sigyz,
     .   sigzx, sigzy, sigzz,
     .   delta0)

*     ----------------------------------------------------
*     This routine calculates sigma_cold in the Stix frame
*     ----------------------------------------------------

      implicit none

      integer lmin, lmax, nzfun, lmaxdim, l, labs, ibessel,
     .    i, j, n, m, nphi

      real xkperp, xkprl, xm, q, xn, xkt, omgc, omgp2, xkb, akb
      real dakbdkb
      real akprl, gammab, alpha, eps0, omgrf, v0i
      real cosalp, sinalp, a, b
      real bx, by, bz
      real xkxsav, xkysav, capr
      real xkphi
      real xkalp, xkbet, xnuomg, xk0, delta0

      complex zi, omgrfc

      complex sig0, sig1, sig2, sig3, sig4, sig5

      complex sigxx, sigxy, sigxz,
     1        sigyx, sigyy, sigyz,
     1        sigzx, sigzy, sigzz

      real uxx, uxy, uxz,
     1     uyx, uyy, uyz,
     2     uzx, uzy, uzz

      parameter (lmaxdim = 99)

      complex tpl, tml, tplp, tplpp, zetalp, zetalm, zieps0,
     .   bl, gamma

      common/sigcom/zi, eps0, v0i, omgrf, xk0

      zieps0 = zi * eps0
      xkphi = nphi / capr
      omgrfc = omgrf * (1. + zi * xnuomg)

      xkalp = uxx * xkxsav + uxy * xkysav + uxz * xkphi
      xkbet = uyx * xkxsav + uyy * xkysav + uyz * xkphi
      xkprl = uzx * xkxsav + uzy * xkysav + uzz * xkphi

      sig0 = 0.0
      sig1 = zieps0 * omgrfc * omgp2 / (omgrfc**2 - omgc**2)
      sig2 = - eps0 * omgc   * omgp2 / (omgrfc**2 - omgc**2)
      sig3 = zieps0 * omgp2 / omgrfc
      sig4 = 0.0
      sig5 = 0.0

*     -------------------
*     Swanson's rotation:
*     -------------------
      sigxx = sig1 + sig0 * xkbet**2
      sigxy = sig2 - sig0 * xkbet * xkalp
      sigxz = sig4 * xkalp + sig5 * xkbet

      sigyx = - sig2 - sig0 * xkbet * xkalp
      sigyy =   sig1 + sig0 * xkalp**2
      sigyz =   sig4 * xkbet - sig5 * xkalp

      sigzx = sig4 * xkalp - sig5 * xkbet
      sigzy = sig4 * xkbet + sig5 * xkalp
      sigzz = sig3

      return

  101 format(i10, 1p8e12.4)
 1314 format(4i10, 1p9e12.4)
 1312 format(1p9e12.4)
  100 format('ier = ', i5, 'besic failed')
      end

c
c***************************************************************************
c

      subroutine besiexp(gamma, lmax, expbes, expbesp, lmaxdim,
     .   expbesovergam)

*     ----------------------------------------------------------
*     Calculates exp(-gamma) times the modified bessel functions
*     (and derivatives) of order up to lmax
*     ----------------------------------------------------------

      implicit none

      integer lmax, nmax, ier, l, lmaxdim
      real gammod
      complex gamma, expbes(0: lmaxdim), expbesp(0: lmaxdim),
     .   expbesovergam(0: lmaxdim),
     .   xil(0: lmaxdim), xilp(0: lmaxdim), exgam

      complex b(100)

      exgam = exp(-gamma)
      gammod = cabs(gamma)

      if(gammod .le. 700.)then
         nmax = lmax + 1
         call besic(gamma, nmax, b, ier)
         if(ier .ne. 0)write(6,100) ier

         do l = 0, lmax
            xil(l) = b(l+1)
         end do

         do l = 0, lmax
           if(l .eq. 0) xilp(0) = xil(1)
           if(l .ne. 0) xilp(l) = xil(l-1) - l / gamma * xil(l)
           expbes(l) = exgam * xil(l)
           expbesp(l) = exgam * xilp(l)
         end do
      end if

      if(gammod .gt. 700.)then
         do l = 0, lmax
            call bes_asym(gamma, l, expbes(l), expbesp(l))
         end do
      end if

      do l = 0, lmax
         expbesovergam(l) = expbes(l) / gamma
      end do

  100 format('ier = ', i5, 'besic failed')
      return
      end

c
c***************************************************************************
c

      subroutine bes_expand(gamma, lmax, expbes, expbesp, lmaxdim,
     .   expbesovergam)

*-------------------------------------------------------------------
*     Calculates exp(-gamma) times the modified bessel functions
*     (and derivatives) of order up to lmax using second order
*     expansion for small argument
*-------------------------------------------------------------------

      implicit none

      integer lmax, nmax, ier, l, lmaxdim
      real factrl, factl
      complex gamma, expbes(0:lmaxdim), expbesp(0:lmaxdim),
     .   expbesovergam(0:lmaxdim), xilovergam,
     .   xil(0:lmaxdim), xilp(0:lmaxdim), exgam

      complex b(100)

      exgam = 1.0 - gamma + gamma**2 / 2.0

      do l = 0, lmax
         factl = factrl(l)
         xil(l) = gamma**l / (2**l * factl) *
     .                                 ( 1. + gamma**2 / (4. * (l+1)))
         xilp(l) = gamma**(l-1) / (2**l * factl) *
     .                     (l + (l+2) * gamma**2 / (4. * (l+1)))

         xilovergam = gamma**(l-1) / (2**l * factl) *
     .                                 ( 1. + gamma**2 / (4. * (l+1)))

         expbes(l) = exgam * xil(l)
         expbesp(l) = exgam * xilp(l)
         expbesovergam(l) = exgam * xilovergam

c         write(6, 100)l, expbes(l), expbesp(l)
      end do

  100 format(i10, 1p8e12.4)

      return
      end

c
c***************************************************************************
c

      subroutine bes_asym(z, n, exil, exilp)

      implicit none

      integer mu, n
      real pi
      complex z, exil, exilp
      data pi/3.141592654/

      mu = 4 * n**2
      exil =  1.0 / csqrt(2.0 * pi * z)
     1   * (1.0
     1   - (mu - 1)/(8.0 * z)
     1   + (mu - 1) * (mu - 9) / (2.0 * (8.0 * z)**2)
     1   - (mu - 1) * (mu - 9) * (mu - 25) / (6.0 * (8.0 * z)**3)  )
      exilp = 1.0 / csqrt(2.0 * pi * z)
     1   * (1.0
     1   - (mu + 3)/(8.0 * z)
     1   + (mu - 1) * (mu + 15) / (2.0 * (8.0 * z)**2)
     1   - (mu - 1) * (mu - 9) * (mu + 35) / (6.0 * (8.0 * z)**3)  )

      return
      end

c
c***************************************************************************
c

      function fzeta (arg)
      complex a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, arg, aux0
     1   , aux1, term, z, zz, fzeta
      data d1r/0.0/
      data d1i/1.77245385090551/
      data d2r/0.0/
      data d2i/3.54490770181103/
      data d3r/0.0/
      data d3i/7.08981540362206/
      data d4/0.33333333333333/
      data eps/1.0E-07/

c     data d4/0.33333333333333/
c     common/zetcom/eps
cray  code analysis
cray  optimize
c
      i = 0
      z = arg
      zz = z*z
      x = real(z)
      y = aimag(z)
      d1 = cmplx(d1r,d1i)
      d2 = cmplx(d2r,d2i)
      d3 = cmplx(d3r,d3i)
      ymag = abs(y)
      if (ymag - 1.0 .ge. 0.) then
c
c     continued fraction method: abs(y).ge.1.0
c
         y0 = y
         y = ymag
         aux1 = 1.5 - z*z
         aux2 = 0.0
         del = 1.5
         a1 = 0.0
         a2 = -1.0
         b1 = 1.0
         b2 = aux1
         c1 = a2/b2
c
  100    continue
         aux1 = aux1 + 2.0
         aux2 = aux2 - del
         del = del + 2.0
         a3 = aux1*a2 + aux2*a1
         b3 = aux1*b2 + aux2*b1
         c2 = a3/b3
         c3 = c2 - c1
         c3r = real(c3)
         c3i = aimag(c3)
         if (abs(c3r) + abs(c3i) .lt. eps) go to 110
         a1 = a2
         a2 = a3
         b1 = b2
         b2 = b3
         c1 = c2
         go to 100
  110    continue
         if (y0 .lt. 0.) then
            y = y0
            c2 = conjg(c2) - d3*z*exp(-zz)
         endif
         aux0 = -(0.5*c2 + 1.0)/z
      else
c
c     asymptotic series method: abs(x).ge.4.0 and abs(y).lt.1.0
c
         xmag = abs(x)
         if (xmag - 4.0 .lt. 0.) go to 130
         term = 1.0/z
         aux0 = -term
         aux1 = 0.5*term**2
         p = 1.0
         if (y .le. 0.) then
            if (y .ne. 0.) then
               aux0 = aux0 + d2*exp(-zz)
            else
               aux0 = aux0 + d1*exp(-zz)
            endif
         endif
  120    continue
         term = aux1*term*p
         aux0 = aux0 - term
         p = p + 2.0
         termr = real(term)
         termi = aimag(term)
c     if(abs(termr)+abs(termi).lt.eps)30,18
         if (abs(termr) + abs(termi) .lt. eps) go to 160
         go to 120
c
c     power series method: abs(x).lt.4.0 and abs(y).lt.1.0
c
  130    continue
         aux0 = 1.0
         aux1 = -(zz + zz)
         aux2 = eps/(eps + xmag + ymag)
         term = d4*aux1
         p = 3.0
  140    continue
         aux0 = aux0 + term
         termr = real(term)
         termi = aimag(term)
c     if(abs(termr)+abs(termi).lt.aux2)26,24
         if (abs(termr) + abs(termi) .lt. aux2) go to 150
         p = p + 2.0
         term = aux1*term/p
         go to 140
  150    continue
         aux0 = d1*exp(-zz) - 2.0*z*aux0
      endif
  160 continue
      fzeta = aux0
      if (i .le. 0) return
      fzeta = -2.0*(1.0 + arg*aux0)
      return
      end

c
c***************************************************************************
c
      real function second1_old(dummy)

      implicit none

      integer :: v(8)
c      integer mtime, mclock
      real dummy

c*****Fortran 90 standard for wall clock time
c      call date_and_time(values=v)
c      second1_old=(v(5)*3600)+(v(6)*60)+v(7)+0.001d0*v(8)

c*****FALCON:
      double precision mpi_wtime
      external mpi_wtime
      second1_old = MPI_WTIME()

c*****EAGLE:
c      mtime = mclock()
c      second1_old  = 0.01 * mtime

      return
      end

c
c***************************************************************************
c

      real function second1(dummy)
!     returns a monotonically increasing wall-clock time-stamp in seconds
!     with an attempt to cross multiple days between calls
!     (assumes one day if end-of-the-month is spanned between calls)

      implicit none

      real dummy

      integer :: v(8)

      integer,save :: ld
      logical,save :: elapsed=.false.

!*****Fortran 90 standard for wall clock time
      call date_and_time(values=v)

      if(elapsed)then  ! second has been called before
        if(v(3).ge.ld)then ! days have advanced monotonically
         second1=(v(3)-ld)*86400+(v(5)*3600)+(v(6)*60)+v(7)+0.001d0*v(8)
        else  ! day counter reset past end of the month
          print *,'WARNING: SECOND1 TRIPPED OVER THE END OF THE MONTH'
          print *,'         SECOND1 IS ASSUMING ONLY ONE ELAPSED DAY'
!         proper fix requires full blown calendar program
          ld=0
         second1=(v(3)-ld)*86400+(v(5)*3600)+(v(6)*60)+v(7)+0.001d0*v(8)
        endif
      else  ! first call
        second1=(v(5)*3600)+(v(6)*60)+v(7)+0.001d0*v(8)
        ld=v(3)
        elapsed=.true.
      endif

      return
      end

max_stars_repo_name: ORNL-Fusion/aorsa

max_issues_repo_name: ORNL-Fusion/aorsa

max_forks_repo_name: ORNL-Fusion/aorsa


! { dg-do run }

!$omp declare reduction (foo : character(kind=1, len=*) &
!$omp & : omp_out = fn (omp_out, omp_in)) initializer (omp_priv = '')
!$omp declare reduction (bar : character(kind=1, len=:) &
!$omp & : omp_out = fn (omp_in, omp_out)) initializer (omp_priv = '')
!$omp declare reduction (baz : character(kind=1, len=1) &
!$omp & : omp_out = char (ichar (omp_out) + ichar (omp_in) &
!$omp & - ichar ('0'))) initializer (omp_priv = '0')
!$omp declare reduction (baz : character(kind=1, len=2) &
!$omp & : omp_out = char (ichar (omp_out(1:1)) + ichar (omp_in(1:1)) &
!$omp & - ichar ('0')) // char (ichar (omp_out(2:2)) + &
!$omp & ichar (omp_in(2:2)) - ichar ('0'))) initializer (omp_priv = '00')
  interface
    elemental function fn (x, y)
      character (len=64), intent (in) :: x, y
      character (len=64) :: fn
    end function
  end interface
  character(kind=1, len=64) :: c(-3:-2,1:1,7:8), d(2:3,-7:-5)
  character(kind = 1, len=1) :: e(2:4)
  character(kind = 1, len=1+1) :: f(8:10,9:10)
  integer :: i, j, k
  c = ''
  d = ''
  e = '0'
  f = '00'
!$omp parallel do reduction (foo : c) reduction (bar : d) &
!$omp & reduction (baz : e, f) private (j, k)
  do i = 1, 64
    forall (j = -3:-2, k = 7:8) &
      c(j,1,k) = trim(c(j,1,k)) // char (ichar ('0') + i)
    d = char (ichar ('0') + i) // d
    e = char (ichar (e) + mod (i, 3))
    f = char (ichar (f(:,:)(1:1)) + mod (i, 2)) &
&	// char (ichar (f(:,:)(2:2)) + mod (i, 3))
  end do
  do i = 1, 64
    if (any (index (c, char (ichar ('0') + i)) .eq. 0)) call abort
    if (any (index (d, char (ichar ('0') + i)) .eq. 0)) call abort
  end do
  if (any (e.ne.char (ichar ('0') + 64))) call abort
  if (any (f(:,:)(1:1).ne.char (ichar ('0') + 32))) call abort
  if (any (f(:,:)(2:2).ne.char (ichar ('0') + 64))) call abort
end
elemental function fn (x, y)
  character (len=64), intent (in) :: x, y
  character (len=64) :: fn
  fn = trim(x) // y
end function


! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::;
!> \callgraph
!! \callergraph
!!  Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by
!!  **ncol**,  and set the values of any piece of
!!  of the patch which extends outside the physical domain
!!  using the boundary conditions.
!!
!!  ### Standard boundary condition choices for amr2ez in clawpack
!!
!!  At each boundary  k = 1 (left),  2 (right),  3 (bottom), 4 (top):
!!
!!  mthbc(k) =
!!  * 0  for user-supplied BC's (must be inserted!)
!!  * 1  for zero-order extrapolation
!!  * 2  for periodic boundary conditions
!!  * 3  for solid walls, assuming this can be implemented
!!                   by reflecting the data about the boundary and then
!!                   negating the 2'nd (for k=1,2) or 3'rd (for k=3,4)
!!                   component of q.
!!  * 4  for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half
!!
!!  The corners of the grid patch are at
!!     (xlo_patch,ylo_patch)  --  lower left corner
!!     (xhi_patch,yhi_patch) --  upper right corner
!!
!!  The physical domain itself is a rectangle bounded by
!!     (xlower,ylower)  -- lower left corner
!!     (xupper,yupper)  -- upper right corner
!!
!   This figure below does not work with doxygen
!   the picture is the following:
!  ____________________________________________________
!
!                _____________________ (xupper,yupper)
!               |                     |
!           ____|____ (xhi_patch,yhi_patch)
!           |   |    |                |
!           |   |    |                |
!           |   |    |                |
!           |___|____|                |
!  (xlo_patch,ylo_patch) |            |
!               |                     |
!               |_____________________|
!    (xlower,ylower)
!  ____________________________________________________
!!
!!
!>  Any cells that lie outside the physical domain are ghost cells whose
!!  values should be set in this routine.  This is tested for by comparing
!!  xlo_patch with xlower to see if values need to be set at the left
!   as in the figure above,
!
!>  and similarly at the other boundaries.
!!  Patches are guaranteed to have at least 1 row of cells filled
!!  with interior values so it is possible to extrapolate.
!!  Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set.
!!
!!  Make sure the order the boundaries are specified is correct
!!  so that diagonal corner cells are also properly taken care of.
!!
!!  Periodic boundaries are set before calling this routine, so if the
!!  domain is periodic in one direction only you
!!  can safely extrapolate in the other direction.
!!
!!  Don't overwrite ghost cells in periodic directions!
!!
!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**)
!! \param aux data array for auxiliary variables
!! \param nrow number of cells in *i* direction on this grid
!! \param ncol number of cells in *j* direction on this grid
!! \param meqn number of equations for the system
!! \param naux number of auxiliary variables
!! \param hx spacing (mesh size) in *i* direction
!! \param hy spacing (mesh size) in *j* direction
!! \param level AMR level of this grid
!! \param time setting ghost cell values at time **time**
!! \param xlo_patch left bound of the input grid
!! \param xhi_patch right bound of the input grid
!! \param ylo_patch lower bound of the input grid
!! \param yhi_patch upper bound of the input grid
! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;

subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time,   &
     xlo_patch, xhi_patch, ylo_patch, yhi_patch, irr, lstgrd)

  use amr_module, only: mthbc, xlower, ylower, xupper, yupper
  use amr_module, only: xperdom,yperdom,spheredom
  use amr_module, only: xcirr,ycirr

  implicit none

  ! Input/Output
  integer, intent(in) :: nrow, ncol, meqn, naux, level
  integer, intent(in) :: irr(nrow,ncol), lstgrd
  real(kind=8), intent(in) :: hx, hy, time
  real(kind=8), intent(in) :: xlo_patch, xhi_patch
  real(kind=8), intent(in) :: ylo_patch, yhi_patch
  real(kind=8), intent(in out) :: val(meqn, nrow, ncol)
  real(kind=8), intent(in out) :: aux(naux, nrow, ncol)

  ! Local storage
  integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt,kuse
  real(kind=8) :: hxmarg, hymarg
  real(kind=8) :: pr,u,v,rho,xcen,ycen
  real(kind=8) :: state(meqn)

  hxmarg = hx * .01d0
  hymarg = hy * .01d0

  ! Use periodic boundary condition specialized code only, if only one
  ! boundary is periodic we still proceed below
  if (xperdom .and. (yperdom .or. spheredom)) then
     return
  end if

  ! Each check has an initial check to ensure that the boundary is a real
  ! boundary condition and otherwise skips the code.  Otherwise
  !-------------------------------------------------------
  ! Left boundary:
  !-------------------------------------------------------
  if (xlo_patch < xlower-hxmarg) then
     ! number of grid cells from this patch lying outside physical domain:
     nxl = int((xlower + hxmarg - xlo_patch) / hx)

     select case(mthbc(1))
     case(0) ! User defined boundary condition
        do j = 1, ncol
           do i=1, nxl
              kuse = irr(i,j)
              call getCellAvgState(meqn,xlo_patch,ylo_patch,hx,hy,   &
                                   state,lstgrd,kuse,i,j,time)
              !if (kuse .ne. -1 .and. kuse .ne. lstgrd) then
              !   xcen = xcirr(kuse)
              !   ycen = ycirr(kuse)
              !else
              !   ycen = ylo_patch + (dfloat(j)-.5d0)* hy
              !   xcen = xlo_patch + (dfloat(i)-.5d0)* hx
              !endif
              !call channelInit(xcen,ycen,state)
              val(1, i, j) = state(1)
              val(2, i, j) = state(2)
              val(3, i, j) = state(3)
              val(4, i, j) = state(4)
           end do
        end do

     case(1) ! Zero-order extrapolation
        do j = 1, ncol
           do i=1, nxl
              val(:, i, j) = val(:, nxl + 1, j)
           end do
        end do

     case(2) ! Periodic boundary condition
        continue

     case(3) ! Wall boundary conditions
        do j = 1, ncol
           do i=1, nxl
              val(:, i, j) = val(:, 2 * nxl + 1 - i, j)
           end do
        end do
        ! negate the normal velocity:
        do j = 1, ncol
           do i=1, nxl
              val(2, i, j) = -val(2, i, j)
           end do
        end do

     case(4) ! Spherical domain
        continue

     case default
        print *, "Invalid boundary condition requested."
        stop
     end select
  end if

  !-------------------------------------------------------
  ! Right boundary:
  !-------------------------------------------------------
  if (xhi_patch > xupper+hxmarg) then

     ! number of grid cells lying outside physical domain:
     nxr = int((xhi_patch - xupper + hxmarg) / hx)
     ibeg = max(nrow - nxr + 1, 1)

     select case(mthbc(2))
     case(0) ! User defined boundary condition
        ! set exact instead of outflow, since this code used for accuracy test
        do j = 1, ncol
           do i=ibeg, nrow
              kuse = irr(i,j)
              call getCellAvgState(meqn,xlo_patch,ylo_patch,hx,hy,   &
                                   state,lstgrd,kuse,i,j,time)

              !if (kuse .ne. -1 .and. kuse .ne. lstgrd) then
              !   xcen = xcirr(kuse)
              !   ycen = ycirr(kuse)
              !else
              !   ycen = ylo_patch + (dfloat(j)-.5d0)* hy
              !   xcen = xlo_patch + (dfloat(i)-.5d0)* hx
              !endif
              !call channelInit(xcen,ycen,state)
              val(1, i, j) = state(1)
              val(2, i, j) = state(2)
              val(3, i, j) = state(3)
              val(4, i, j) = state(4)
           end do
        end do

     case(1) ! Zero-order extrapolation
        do i = ibeg, nrow
           do j = 1, ncol
              val(:, i, j) = val(:, ibeg - 1, j)
           end do
        end do

     case(2) ! Periodic boundary condition
        continue

     case(3) ! Wall boundary conditions
        do i=ibeg, nrow
           do j = 1, ncol
              val(:, i, j) = val(:, 2 * ibeg - 1 - i, j)
           end do
        end do
        ! negate the normal velocity:
        do i = ibeg, nrow
           do j = 1, ncol
              val(2, i, j) = -val(2, i, j)
           end do
        end do

     case(4) ! Spherical domain
        continue

     case default
        print *, "Invalid boundary condition requested."
        stop

     end select
  end if

  !-------------------------------------------------------
  ! Bottom boundary:
  !-------------------------------------------------------
  if (ylo_patch < ylower - hymarg) then

     ! number of grid cells lying outside physical domain:
     nyb = int((ylower + hymarg - ylo_patch) / hy)

     select case(mthbc(3))
     case(0) ! User defined boundary condition
        ! Replace this code with a user defined boundary condition
        stop "A user defined boundary condition was not provided."

     case(1) ! Zero-order extrapolation
        do j = 1, nyb
           do i = 1, nrow
              val(:,i,j) = val(:, i, nyb + 1)
           end do
        end do

     case(2) ! Periodic boundary condition
        continue

     case(3) ! Wall boundary conditions
        do j = 1, nyb
           do i = 1, nrow
              val(:,i,j) = val(:, i, 2 * nyb + 1 - j)
           end do
        end do
        ! negate the normal velocity:
        do j = 1, nyb
           do i = 1, nrow
              val(3,i,j) = -val(3, i, j)
           end do
        end do

     case(4) ! Spherical domain
        continue

     case default
        print *, "Invalid boundary condition requested."
        stop

     end select
  end if

  !-------------------------------------------------------
  ! Top boundary:
  !-------------------------------------------------------
  if (yhi_patch > yupper + hymarg) then

     ! number of grid cells lying outside physical domain:
     nyt = int((yhi_patch - yupper + hymarg) / hy)
     jbeg = max(ncol - nyt + 1, 1)

     select case(mthbc(4))
     case(0) ! User defined boundary condition
        ! Replace this code with a user defined boundary condition
        stop "A user defined boundary condition was not provided."

     case(1) ! Zero-order extrapolation
        do j = jbeg, ncol
           do i = 1, nrow
              val(:, i, j) = val(:, i, jbeg - 1)
           end do
        end do

     case(2) ! Periodic boundary condition
        continue

     case(3) ! Wall boundary conditions
        do j = jbeg, ncol
           do i = 1, nrow
              val(:, i, j) = val(:, i, 2 * jbeg - 1 - j)
           end do
        end do
        ! negate the normal velocity:
        do j = jbeg, ncol
           do i = 1, nrow
              val(3, i, j) = -val(3, i, j)
           end do
        end do

     case(4) ! Spherical domain
        continue

     case default
        print *, "Invalid boundary condition requested."
        stop

     end select
  end if

end subroutine bc2amr


text:
C **********************************************************************SPLD   1
C ** PROGRAM SPLOAD - LOADS NAD 1927 AND NAD 1983 DIRECT ACCESS FILES **SPLD   2
C ** U. S. GEOLOGICAL SURVEY - CHIPPEAUX AND LINCK           09/19/89 **SPLD   3
C **********************************************************************SPLD   4
      PROGRAM SPLOAD                                                    SPLD   5
C                                                                       SPLD   6
C     PROGRAM TO COPY NAD 1927 AND NAD 1983 SEQUENTIAL FILES TO DIRECT  SPLD   7
C     ACCESS FILES FOR USE BY GCTP                                      SPLD   8
C                                                                       SPLD   9
      REAL*8 PARM(9)                                                    SPLD  10
      INTEGER*4 I, ID, IZONE, LENGTH                                    SPLD  11
      CHARACTER*32 ZNAME, GETF27, GETF83, PUTF27, PUTF83                SPLD  12
      CHARACTER*3 STAT27, STAT83                                        SPLD  13
C                                                                       SPLD  14
C     IPR    = LOGICAL UNIT NUMBER OF PRINTER                           SPLD  15
C     GETF27 = NAME OF NAD 1927 INPUT FILE                              SPLD  16
C     PUTF27 = NAME OF NAD 1927 OUTPUT FILE                             SPLD  17
C     GETF83 = NAME OF NAD 1983 INPUT FILE                              SPLD  18
C     PUTF83 = NAME OF NAD 1983 OUTPUT FILE                             SPLD  19
C     STAT27 = STATUS OF NAD 1927 OUTPUT FILE                           SPLD  20
C     STAT83 = STATUS OF NAD 1983 OUTPUT FILE                           SPLD  21
C     LENGTH = LENGTH IN BYTES OR WORDS OF OUTPUT RECORDS               SPLD  22
C                                                                       SPLD  23
      PARAMETER (IPR = 6)                                               SPLD  24
      PARAMETER (GETF27 = 'nad1927', GETF83 = 'nad1983')                SPLD  25
      PARAMETER (PUTF27 = 'nad27sp', PUTF83 = 'nad83sp')                SPLD  26
      PARAMETER (STAT27 = 'NEW', STAT83 = 'NEW')                        SPLD  27
      PARAMETER (LENGTH = 108)                                          SPLD  28
C                                                                       SPLD  29
      OPEN(UNIT=2,FILE=GETF27,STATUS='OLD',                             SPLD  30
     .     ACCESS='SEQUENTIAL')                                         SPLD  31
      OPEN(UNIT=3,FILE=PUTF27,STATUS=STAT27,RECL=LENGTH,                SPLD  32
     .     ACCESS='DIRECT')                                             SPLD  33
C                                                                       SPLD  34
      OPEN(UNIT=12,FILE=GETF83,STATUS='OLD',                            SPLD  35
     .     ACCESS='SEQUENTIAL')                                         SPLD  36
      OPEN(UNIT=13,FILE=PUTF83,STATUS=STAT83,RECL=LENGTH,               SPLD  37
     .     ACCESS='DIRECT')                                             SPLD  38
C                                                                       SPLD  39
    1 FORMAT(A32,7X,I1,31X,I4)                                          SPLD  40
    2 FORMAT(3D25.16)                                                   SPLD  41
    3 FORMAT(1X,I3,5X,A32,5X,I4,5X,I1)                                  SPLD  42
C                                                                       SPLD  43
      I = 0                                                             SPLD  44
   10 READ(2,1,END=20) ZNAME, ID, IZONE                                 SPLD  45
      I = I + 1                                                         SPLD  46
      READ(2,2) PARM(1), PARM(2), PARM(3)                               SPLD  47
      READ(2,2) PARM(4), PARM(5), PARM(6)                               SPLD  48
      READ(2,2) PARM(7), PARM(8), PARM(9)                               SPLD  49
      WRITE(3,REC=I) ZNAME, ID, PARM                                    SPLD  50
      WRITE(IPR,3) I, ZNAME, IZONE, ID                                  SPLD  51
      GO TO 10                                                          SPLD  52
C                                                                       SPLD  53
   20 I = 0                                                             SPLD  54
   30 READ(12,1,END=40) ZNAME, ID, IZONE                                SPLD  55
      I = I + 1                                                         SPLD  56
      READ(12,2) PARM(1), PARM(2), PARM(3)                              SPLD  57
      READ(12,2) PARM(4), PARM(5), PARM(6)                              SPLD  58
      READ(12,2) PARM(7), PARM(8), PARM(9)                              SPLD  59
      WRITE(13,REC=I)  ZNAME, ID, PARM                                  SPLD  60
      WRITE(IPR,3) I, ZNAME, IZONE, ID                                  SPLD  61
      GO TO 30                                                          SPLD  62
C                                                                       SPLD  63
   40 CLOSE(2)                                                          SPLD  64
      CLOSE(3)                                                          SPLD  65
      CLOSE(12)                                                         SPLD  66
      CLOSE(13)                                                         SPLD  67
C                                                                       SPLD  68
      STOP                                                              SPLD  69
      END                                                               SPLD  70

ext: f

max_stars_repo_name: SCECcode/cvmhlabn

max_issues_repo_name: SCECcode/cvmhlabn

max_forks_repo_name: SCECcode/cvmhlabn


          program demo_sort2020
          use M_strings, only : split2020
          implicit none
          character(len=*),parameter :: gen='(*("[",g0,"]":,","))'

           ! Execution of TOKEN form
           block
             character (len=:), allocatable :: string
             character (len=:), allocatable :: tokens(:)
             character (len=*),parameter :: set = " ,"
             string = 'first,second,third'
             call split2020(string, set, tokens )
             write(*,gen)tokens

           ! assigns the value ['first ','second','third ' ]
           ! to TOKENS.
           endblock

           ! Execution of BOUNDS form

           block
             character (len=:), allocatable :: string
             character (len=*),parameter :: set = " ,"
             integer, allocatable        :: first(:), last(:)
             string =    'first,second,,forth'
             call split2020 (string, set, first, last)
             write(*,gen)first
             write(*,gen)last

           ! will assign the value [ 1, 7, 14, 15 ] to FIRST,
           ! and the value [ 5, 12, 13, 19 ] to LAST.
           endblock

           ! Execution of STEP form
           block
             character (len=:), allocatable :: string
             character (len=*),parameter :: set = " ,"
             integer :: p, istart, iend
             string = " one,   last  example  "
             do while (p < len(string))
               istart = p + 1
               call split2020 (string, set, p)
               iend=p-1
               if(iend.gt.istart)then
                  print '(t3,a,1x,i0,1x,i0)', string (istart:iend),istart,iend
               endif
             enddo
           endblock
       end program demo_sort2020


module ArrayOperationClass
    use MathClass
    implicit none

    interface MergeArray
        module procedure MergeArrayInt, MergeArrayReal
    end interface MergeArray

    interface CopyArray
        module procedure CopyArrayInt, CopyArrayReal,CopyArrayIntVec, CopyArrayRealVec
    end interface CopyArray

    interface TrimArray
        module procedure TrimArrayInt, TrimArrayReal
    end interface TrimArray

    interface ImportArray
        module procedure ImportArrayInt, ImportArrayReal
    end interface ImportArray

    interface ExportArray
        module procedure ExportArrayInt, ExportArrayReal
    end interface ExportArray

    interface ExportArraySize
        module procedure ExportArraySizeInt, ExportArraySizeReal
    end interface ExportArraySize

    interface InOrOut
            module procedure InOrOutReal,InOrOutInt
    end interface

    interface ShowArray
        module procedure ShowArrayInt, ShowArrayReal
    end interface ShowArray

    interface ShowArraySize
        module procedure    ShowArraySizeInt, ShowArraySizeReal
        module procedure    ShowArraySizeIntvec, ShowArraySizeRealvec
        module procedure    ShowArraySizeIntThree, ShowArraySizeRealThree
    end interface ShowArraySize

    interface ExtendArray
        module procedure  :: ExtendArrayReal,ExtendArrayInt
    end interface ExtendArray

    interface insertArray
        module procedure :: insertArrayInt, insertArrayReal
    end interface insertArray

    interface removeArray
        module procedure :: removeArrayReal,removeArrayInt
    end interface removeArray
contains

! ############## Elementary Entities ##############

!=====================================
subroutine MergeArrayInt(a,b,c)
    integer,intent(in)::a(:,:)
    integer,intent(in)::b(:,:)
    integer,allocatable,intent(out)::c(:,:)
    integer i,j,an,am,bn,bm

    if(allocated(c)) deallocate(c)
    an=size(a,1)
    am=size(a,2)

    bn=size(b,1)
    bm=size(b,2)

    if(am/=bm)then
        print *, "ERROR :: MergeArray, size(a,2)/= size(b,2)"
        return
    endif
    allocate(c(an+bn,am) )
    do i=1,an
        c(i,:)=a(i,:)
    enddo
    do i=1,bn
        c(i+an,:)=b(i,:)
    enddo

end subroutine
!=====================================

!=====================================
subroutine MergeArrayReal(a,b,c)
    real(8),intent(in)::a(:,:)
    real(8),intent(in)::b(:,:)
    real(8),allocatable,intent(out)::c(:,:)
    integer i,j,an,am,bn,bm
    if(allocated(c)) deallocate(c)
    an=size(a,1)
    am=size(a,2)

    bn=size(b,1)
    bm=size(b,2)

    if(am/=bm)then
        print *, "ERROR :: MergeArray, size(a,2)/= size(b,2)"
        return
    endif
    allocate(c(an+bn,am) )
    do i=1,an
        c(i,:)=a(i,:)
    enddo
    do i=1,bn
        c(i+an,:)=b(i,:)
    enddo

end subroutine
!=====================================

!=====================================
subroutine CopyArrayInt(a,ac)
    integer,allocatable,intent(inout)::a(:,:)
    integer,allocatable,intent(inout)::ac(:,:)
    integer i,j,n,m

    if(.not.allocated(a) )then
        if(allocated(ac) )then
            deallocate(ac)
            return
        endif
    endif
    n=size(a,1)
    m=size(a,2)
    if(allocated(ac) ) deallocate(ac)

    allocate(ac(n,m) )
    ac(:,:)=a(:,:)

end subroutine
!=====================================

!=====================================
subroutine CopyArrayReal(a,ac)
    real(8),allocatable,intent(inout)::a(:,:)
    real(8),allocatable,intent(inout)::ac(:,:)
    integer i,j,n,m

    if(.not.allocated(a) )then
        if(allocated(ac) )then
            deallocate(ac)
            return
        endif
    endif
    n=size(a,1)
    m=size(a,2)

    if(allocated(ac) ) deallocate(ac)
    allocate(ac(n,m) )
    ac(:,:)=a(:,:)

end subroutine
!=====================================

!=====================================
subroutine CopyArrayIntVec(a,ac)
    integer,allocatable,intent(inout)::a(:)
    integer,allocatable,intent(inout)::ac(:)
    integer i,j,n,m

    if(.not.allocated(a) )then
        if(allocated(ac) )then
            deallocate(ac)

        endif
        return
    endif
    n=size(a,1)
    if(allocated(ac) ) deallocate(ac)

    allocate(ac(n) )
    ac(:)=a(:)

end subroutine
!=====================================

!=====================================
subroutine CopyArrayRealVec(a,ac)
    real(8),allocatable,intent(inout)::a(:)
    real(8),allocatable,intent(inout)::ac(:)
    integer i,j,n,m

    if(.not.allocated(a) )then
        if(allocated(ac) )then
            deallocate(ac)

        endif
        return
    endif
    n=size(a,1)
    if(allocated(ac) ) deallocate(ac)

    allocate(ac(n) )
    ac(:)=a(:)

end subroutine
!=====================================

!=====================================
subroutine TrimArrayInt(a,k)
    integer,allocatable,intent(inout)::a(:,:)
    integer,intent(in)::k
    integer,allocatable::ac(:,:)
    integer :: i,j,n,m

    n=size(a,1)
    m=size(a,2)
    allocate(ac(k,m ))

    do i=1,k
        ac(i,:) = a(i,:)
    enddo
    deallocate(a)
    allocate(a(k,m) )
    a(:,:)=ac(:,:)
    return

end subroutine
!=====================================

!=====================================
subroutine TrimArrayReal(a,k)
    real(8),allocatable,intent(inout)::a(:,:)
    integer,intent(in)::k
    real(8),allocatable::ac(:,:)
    integer :: i,j,n,m

    n=size(a,1)
    m=size(a,2)
    allocate(ac(k,m ))

    do i=1,k
        ac(i,:) = a(i,:)
    enddo

    deallocate(a)
    allocate(a(k,m) )
    a(:,:)=ac(:,:)
    return

end subroutine
!=====================================

!##################################################
subroutine ImportArrayInt(Mat,OptionalFileHandle,OptionalSizeX,OptionalSizeY)
    integer,allocatable,intent(inout)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle,OptionalSizeX,OptionalSizeY

    integer i,j,n,m,fh

    if(present(OptionalSizeX) )then
        n=OptionalSizeX
    elseif(allocated(Mat) )then
        n=size(Mat,1)
    else
        n=1
        print *, "Caution :: ArrayOperationClass/ImportArray >> No size_X is set"
    endif

    if(present(OptionalSizeY) )then
        m=OptionalSizeY
    elseif(allocated(Mat) )then
        m=size(Mat,2)
    else
        m=1
        print *, "Caution :: ArrayOperationClass/ImportArray >> No size_Y is set"
    endif

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(.not.allocated(Mat))then
        allocate(Mat(n,m) )
    endif

    if(size(Mat,1)/=n .or. size(Mat,2)/=m  )then
        deallocate(Mat)
        allocate(Mat(n,m) )
    endif

    do i=1,size(Mat,1)
        read(fh,*) Mat(i,:)
    enddo

    !include "./ImportArray.f90"

end subroutine ImportArrayInt
!##################################################

!##################################################
subroutine ImportArrayReal(Mat,OptionalFileHandle,OptionalSizeX,OptionalSizeY)
    real(8),allocatable,intent(inout)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle,OptionalSizeX,OptionalSizeY

    !include "./ImportArray.f90"

    integer i,j,n,m,fh

    if(present(OptionalSizeX) )then
        n=OptionalSizeX
    elseif(allocated(Mat) )then
        n=size(Mat,1)
    else
        n=1
        print *, "Caution :: ArrayOperationClass/ImportArray >> No size_X is set"
    endif

    if(present(OptionalSizeY) )then
        m=OptionalSizeY
    elseif(allocated(Mat) )then
        m=size(Mat,2)
    else
        m=1
        print *, "Caution :: ArrayOperationClass/ImportArray >> No size_Y is set"
    endif

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(.not.allocated(Mat))then
        allocate(Mat(n,m) )
    endif

    if(size(Mat,1)/=n .or. size(Mat,2)/=m  )then
        deallocate(Mat)
        allocate(Mat(n,m) )
    endif

    do i=1,size(Mat,1)
        read(fh,*) Mat(i,:)
    enddo

end subroutine ImportArrayReal
!##################################################

!##################################################
subroutine ExportArraySizeInt(Mat,RankNum,OptionalFileHandle)
    integer,intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle
    integer,intent(in)::RankNum

    !#include "./ExportArraySize.f90"
    integer :: fh
    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    endif

    write(fh,*) size(Mat,RankNum)

end subroutine ExportArraySizeInt
!##################################################

!##################################################
subroutine ExportArraySizeReal(Mat,RankNum,OptionalFileHandle)
    real(8),intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle
    integer,intent(in)::RankNum

    !#include "./ExportArraySize.f90"
    integer :: fh
    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    endif

    write(fh,*) size(Mat,RankNum)

end subroutine ExportArraySizeReal
!##################################################

!##################################################
subroutine ExportArrayInt(Mat,OptionalFileHandle)
    integer,intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    do i=1,size(Mat,1)
        write(fh,*) Mat(i,:)
    enddo

end subroutine ExportArrayInt
!##################################################

!##################################################
subroutine ExportArrayReal(Mat,OptionalFileHandle)
    real(8),intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    do i=1,size(Mat,1)
        write(fh,*) Mat(i,:)
    enddo

end subroutine ExportArrayReal
!##################################################

!##################################################
subroutine ShowArrayInt(Mat,IndexArray,FileHandle,Name)
    integer,intent(in)::Mat(:,:)
    integer,optional,intent(in) :: IndexArray(:,:)
    integer,optional,intent(in)::FileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i,j,k

    if(present(FileHandle) )then
        fh=FileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(present(IndexArray))then

        do i=1,size(IndexArray,1)
            do j=1,size(IndexArray,2)
                k = IndexArray(i,j)
                if(k <= 0)then
                    cycle
                endif

                if(present(FileHandle) .or. present(Name) )then
                    write(fh,*) Mat(k,:)
                else
                    print *, Mat(k,:)
                endif
            enddo
        enddo
    else

        do j=1,size(Mat,1)

            if(present(FileHandle) .or. present(Name) )then
                write(fh,*) Mat(j,:)
            else
                print *, Mat(j,:)
            endif

        enddo

    endif

    if(present(FileHandle) .or. present(Name) )then
        flush(fh)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArrayReal(Mat,IndexArray,FileHandle,Name)
    real(8),intent(in)::Mat(:,:)
    integer,optional,intent(in) :: IndexArray(:,:)
    integer,optional,intent(in)::FileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i,j,k

    if(present(FileHandle) )then
        fh=FileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(present(IndexArray))then

        do i=1,size(IndexArray,1)
            do j=1,size(IndexArray,2)
                k = IndexArray(i,j)
                if(k <= 0)then
                    cycle
                endif

                if(present(FileHandle) .or. present(Name) )then
                    write(fh,*) Mat(k,:)
                else
                    print *, Mat(k,:)
                endif
            enddo
        enddo
    else

        do j=1,size(Mat,1)

            if(present(FileHandle) .or. present(Name) )then
                write(fh,*) Mat(j,:)
            else
                print *, Mat(j,:)
            endif

        enddo

    endif

    if(present(FileHandle) .or. present(Name) )then
        flush(fh)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeInt(Mat,OptionalFileHandle,Name)
    integer,allocatable,intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle) )then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1), size(Mat,2)
    if(present(OptionalFileHandle) .or. present(Name) )then
        write(fh,*) size(Mat,1), size(Mat,2)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeReal(Mat,OptionalFileHandle,Name)
    real(8),allocatable,intent(in)::Mat(:,:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle) )then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1), size(Mat,2)
    if(present(OptionalFileHandle).or. present(Name) )then
        write(fh,*) size(Mat,1), size(Mat,2)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeIntvec(Mat,OptionalFileHandle,Name)
    integer,allocatable,intent(in)::Mat(:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle).or. present(Name) )then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1)
    if(present(OptionalFileHandle).or. present(Name) )then
        write(fh,*) size(Mat,1)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeRealvec(Mat,OptionalFileHandle,Name)
    real(8),allocatable,intent(in)::Mat(:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle).or. present(Name) )then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1)
    if(present(OptionalFileHandle) .or. present(Name))then
        write(fh,*) size(Mat,1)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeIntThree(Mat,OptionalFileHandle,Name)
    integer,allocatable,intent(in)::Mat(:,:,:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle) .or. present(Name))then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1), size(Mat,2) , size(Mat,3)
    if(present(OptionalFileHandle) .or. present(Name))then
        write(fh,*) size(Mat,1), size(Mat,2), size(Mat,3)
    endif

    if(present(Name) )then
        close(fh)
    endif
end subroutine
!##################################################

!##################################################
subroutine ShowArraySizeRealThree(Mat,OptionalFileHandle,Name)
    real(8),allocatable,intent(in)::Mat(:,:,:,:)
    integer,optional,intent(in)::OptionalFileHandle
    character(*),optional,intent(in)::Name

    !#include "./ExportArray.f90"
    integer :: fh,i

    if(present(OptionalFileHandle) )then
        fh=OptionalFileHandle
    else
        fh=10
    endif

    if(present(Name) )then
        open(fh,file=Name)
    endif

    if(.not. allocated(Mat) )then
        print *, "Not allocated!"
        if(present(OptionalFileHandle) .or. present(Name))then
            write(fh,*) "Not allocated!"
        endif
    endif
    print *, size(Mat,1), size(Mat,2) , size(Mat,3)
    if(present(OptionalFileHandle).or. present(Name) )then
        write(fh,*) size(Mat,1), size(Mat,2), size(Mat,3)
    endif

    if(present(Name) )then
        close(fh)
    endif

end subroutine
!##################################################

!##################################################
function InOrOutReal(x,xmax,xmin,DimNum) result(Inside)
    real(8),intent(in)::x(:)
    real(8),intent(in)::xmax(:),xmin(:)
    integer,optional,intent(in)::DimNum
    integer :: dim_num
    logical ::Inside
    integer :: i,j,n,cout

    cout=0
    if(present(DimNum) )then
        dim_num=DimNum
    else
        dim_num=size(x)
    endif

    do i=1,dim_num
        if(xmin(i) <= x(i) .and. x(i)<=xmax(i) )then
            cout=cout+1
        else
            cycle
        endif
    enddo

    if(cout==dim_num)then
        Inside = .true.
    else
        Inside = .false.
    endif

end function

!##################################################

!##################################################
function InOrOutInt(x,xmax,xmin,DimNum) result(Inside)
    integer,intent(in)::x(:)
    integer,intent(in)::xmax(:),xmin(:)
    integer,optional,intent(in)::DimNum
    integer :: dim_num
    logical ::Inside
    integer :: i,j,n,cout

    cout=0
    if(present(DimNum) )then
        dim_num=DimNum
    else
        dim_num=size(x)
    endif

    do i=1,dim_num
        if(xmin(i) <= x(i) .and. x(i)<=xmax(i) )then
            cout=cout+1
        else
            cycle
        endif
    enddo

    if(cout==dim_num)then
        Inside = .true.
    else
        Inside = .false.
    endif

end function

!##################################################

!##################################################
subroutine ExtendArrayReal(mat,extend1stColumn,extend2ndColumn,DefaultValue)
    real(8),allocatable,intent(inout)::mat(:,:)
    real(8),allocatable :: buffer(:,:)
    logical,optional,intent(in) :: extend1stColumn,extend2ndColumn
    real(8),optional,intent(in) :: DefaultValue
    real(8) :: val
    integer :: n,m

    if( present(DefaultValue) )then
        val = DefaultValue
    else
        val = 0.0d0
    endif

    n=size(mat,1)
    m=size(mat,2)
    if(present(extend1stColumn) )then
        if(extend1stColumn .eqv. .true.)then
            allocate(buffer(n+1, m ) )
            buffer(:,:)=val
            buffer(1:n,1:m)=mat(1:n,1:m)
            deallocate(mat)
            call copyArray(buffer,mat)
            deallocate(buffer)
        endif
    endif

    n=size(mat,1)
    m=size(mat,2)
    if(present(extend2ndColumn) )then
        if(extend2ndColumn .eqv. .true.)then
            allocate(buffer(n, m+1 ) )
            buffer(:,:)=val
            buffer(1:n,1:m)=mat(1:n,1:m)
            deallocate(mat)
            call copyArray(buffer,mat)
            deallocate(buffer)
        endif
    endif

end subroutine
!##################################################

!##################################################
subroutine ExtendArrayInt(mat,extend1stColumn,extend2ndColumn,DefaultValue)
    integer,allocatable,intent(inout)::mat(:,:)
    integer,allocatable :: buffer(:,:)
    logical,optional,intent(in) :: extend1stColumn,extend2ndColumn
    integer,optional,intent(in) :: DefaultValue
    integer :: val
    integer :: i,j,k,n,m

    if( present(DefaultValue) )then
        val = DefaultValue
    else
        val = 0
    endif

    n=size(mat,1)
    m=size(mat,2)
    if(present(extend1stColumn) )then
        if(extend1stColumn .eqv. .true.)then
            allocate(buffer(n+1, m ) )
            buffer(:,:)=val
            buffer(1:n,1:m)=mat(1:n,1:m)
            deallocate(mat)
            call copyArray(buffer,mat)
            deallocate(buffer)
        endif
    endif

    n=size(mat,1)
    m=size(mat,2)
    if(present(extend2ndColumn) )then
        if(extend2ndColumn .eqv. .true.)then
            allocate(buffer(n, m+1 ) )
            buffer(:,:)=val
            buffer(1:n,1:m)=mat(1:n,1:m)
            deallocate(mat)
            call copyArray(buffer,mat)
            deallocate(buffer)
        endif
    endif

end subroutine
!##################################################

!##################################################
subroutine insertArrayInt(mat,insert1stColumn,insert2ndColumn,DefaultValue,NextOf)
    integer,allocatable,intent(inout)::mat(:,:)
    integer,allocatable :: buffer(:,:)
    logical,optional,intent(in) :: insert1stColumn,insert2ndColumn
    integer,optional,intent(in) :: DefaultValue,NextOf
    integer :: val
    integer :: i,nof

    call extendArray(mat,insert1stColumn,insert2ndColumn,DefaultValue)

    if(present(DefaultValue))then
        val =DefaultValue
    else
        val =0
    endif

    if(present(NextOf))then
        nof=NextOf
    else
        if( present(insert1stColumn ))then
            if( insert1stColumn .eqv. .true. )then
                nof=size(mat,1)-1
            endif
        endif

        if( present(insert2ndColumn ))then
            if( insert2ndColumn .eqv. .true. )then
                nof=size(mat,2)-1
            endif
        endif

    endif

    if( present(insert1stColumn ))then
        if( insert1stColumn .eqv. .true. )then
            do i=size(mat,1)-1,nof,-1
                mat(i+1,:)=mat(i,:)
            enddo
            mat(nof+1,:)=val
        endif
    endif
    if( present(insert2ndColumn ))then
        if( insert2ndColumn .eqv. .true. )then
            do i=size(mat,1)-1,nof,-1
                mat(:,i+1)=mat(:,i)
            enddo
            mat(:,nof+1)=val
        endif
    endif
end subroutine
!##################################################

!##################################################
subroutine insertArrayReal(mat,insert1stColumn,insert2ndColumn,DefaultValue,NextOf)
    real(8),allocatable,intent(inout)::mat(:,:)
    real(8),allocatable :: buffer(:,:)
    logical,optional,intent(in) :: insert1stColumn,insert2ndColumn

    integer,optional,intent(in) :: NextOf
    real(8),optional,intent(in) :: DefaultValue
    real(8) :: val
    integer :: i,nof

    call extendArray(mat,insert1stColumn,insert2ndColumn,DefaultValue)

    if(present(DefaultValue))then
        val =DefaultValue
    else
        val =0
    endif

    if(present(NextOf))then
        nof=NextOf
    else
        if( present(insert1stColumn ))then
            if( insert1stColumn .eqv. .true. )then
                nof=size(mat,1)-1
            endif
        endif

        if( present(insert2ndColumn ))then
            if( insert2ndColumn .eqv. .true. )then
                nof=size(mat,2)-1
            endif
        endif

    endif

    if( present(insert1stColumn ))then
        if( insert1stColumn .eqv. .true. )then
            do i=size(mat,1)-1,nof,-1
                mat(i+1,:)=mat(i,:)
            enddo
            mat(nof+1,:)=val
        endif
    endif
    if( present(insert2ndColumn ))then
        if( insert2ndColumn .eqv. .true. )then
            do i=size(mat,1)-1,nof,-1
                mat(:,i+1)=mat(:,i)
            enddo
            mat(:,nof+1)=val
        endif
    endif
end subroutine
!##################################################

!##################################################
subroutine removeArrayInt(mat,remove1stColumn,remove2ndColumn,NextOf)
    integer,allocatable,intent(inout)::mat(:,:)
    integer,allocatable :: buffer(:,:)
    logical,optional,intent(in) :: remove1stColumn,remove2ndColumn

    integer,optional,intent(in) :: NextOf

    integer :: i,nof

    if( present(remove1stColumn ))then
        if( remove1stColumn .eqv. .true. )then
            nof=size(mat,1)
        endif
    endif

    if( present(remove2ndColumn ))then
        if( remove2ndColumn .eqv. .true. )then
            nof=size(mat,2)
        endif
    endif

    if(present(NextOf) )then
        if(NextOf >= nof)then
            return
        endif
    endif
    call copyArray(mat,buffer)

    if(present(NextOf))then
        nof=NextOf
    else
        if( present(remove1stColumn ))then
            if( remove1stColumn .eqv. .true. )then
                nof=size(mat,1)-1
            endif
        endif

        if( present(remove2ndColumn ))then
            if( remove2ndColumn .eqv. .true. )then
                nof=size(mat,2)-1
            endif
        endif
    endif

    deallocate(mat)

    if( present(remove1stColumn ))then
        if( remove1stColumn .eqv. .true. )then
            if(size(buffer,1)-1 == 0)then
                print *, "Array is deleted"
                return
            endif
            allocate(mat(size(buffer,1)-1,size(buffer,2)) )
            mat(:,:)=0.0d0
            do i=1,nof
                mat(i,:)=buffer(i,:)
            enddo

            do i=nof+1,size(buffer,1)-1
                mat(i,:)=buffer(i+1,:)
            enddo

        endif
    endif
    if( present(remove2ndColumn ))then
        if( remove2ndColumn .eqv. .true. )then

            if(size(buffer,2)-1 == 0)then
                print *, "Array is deleted"
                return
            endif
            allocate(mat(size(buffer,1),size(buffer,2)-1) )
            mat(:,:)=0.0d0
            do i=1,nof
                mat(:,i)=buffer(:,i)
            enddo

            do i=nof+1,size(buffer,2)-1
                mat(:,i)=buffer(:,i+1)
            enddo
        endif
    endif

end subroutine
!##################################################

!##################################################
subroutine removeArrayReal(mat,remove1stColumn,remove2ndColumn,NextOf)
    real(8),allocatable,intent(inout)::mat(:,:)
    real(8),allocatable :: buffer(:,:)
    logical,optional,intent(in) :: remove1stColumn,remove2ndColumn

    integer,optional,intent(in) :: NextOf

    integer :: i,nof

    if( present(remove1stColumn ))then
        if( remove1stColumn .eqv. .true. )then
            nof=size(mat,1)
        endif
    endif

    if( present(remove2ndColumn ))then
        if( remove2ndColumn .eqv. .true. )then
            nof=size(mat,2)
        endif
    endif

    if(present(NextOf) )then
        if(NextOf >= nof)then
            return
        endif
    endif

    call copyArray(mat,buffer)

    if(present(NextOf))then
        nof=NextOf
    else
        if( present(remove1stColumn ))then
            if( remove1stColumn .eqv. .true. )then
                nof=size(mat,1)-1
            endif
        endif

        if( present(remove2ndColumn ))then
            if( remove2ndColumn .eqv. .true. )then
                nof=size(mat,2)-1
            endif
        endif
    endif

    deallocate(mat)

    if( present(remove1stColumn ))then
        if( remove1stColumn .eqv. .true. )then
            if(size(buffer,1)-1 == 0)then
                print *, "Array is deleted"
                return
            endif
            allocate(mat(size(buffer,1)-1,size(buffer,2)) )
            mat(:,:)=0.0d0
            do i=1,nof
                mat(i,:)=buffer(i,:)
            enddo

            do i=nof+1,size(buffer,1)-1
                mat(i,:)=buffer(i+1,:)
            enddo

        endif
    endif
    if( present(remove2ndColumn ))then
        if( remove2ndColumn .eqv. .true. )then

            if(size(buffer,2)-1 == 0)then
                print *, "Array is deleted"
                return
            endif
            allocate(mat(size(buffer,1),size(buffer,2)-1) )
            mat(:,:)=0.0d0
            do i=1,nof
                mat(:,i)=buffer(:,i)
            enddo

            do i=nof+1,size(buffer,2)-1
                mat(:,i)=buffer(:,i+1)
            enddo
        endif
    endif

end subroutine
!##################################################

end module ArrayOperationClass