text:
!*******************************************************************************
!>  @file unprimed_grid.f
!>  @brief Contains module @ref unprimed_grid.
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  Defines the base class of the type @ref unprimed_grid_class. This contains
!>  the state variables to define the currents and positions of the volumn
!>  integral.
!*******************************************************************************
      MODULE unprimed_grid
      USE stel_kinds, ONLY: rprec
      USE profiler
      USE bmw_parallel_context
      USE m_grid
      USE primed_grid

      IMPLICIT NONE

!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) unprimed grid base class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Base class representing a unprimed grid. This is grid the volume integral
!>  will be summed over.
!-------------------------------------------------------------------------------
      TYPE :: unprimed_grid_class
!>  Vector potential in the R direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: a_r => null()
!>  Vector potential in the Phi direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: a_p => null()
!>  Vector potential in the Z direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: a_z => null()

!>  B field in the R direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: b_r => null()
!>  B field in the Phi direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: b_p => null()
!>  B field in the Z direction.
         REAL (rprec), DIMENSION(:,:,:), POINTER :: b_z => null()
      CONTAINS
         FINAL :: unprimed_grid_destruct
      END TYPE

!-------------------------------------------------------------------------------
!>  Interface to constructors.
!-------------------------------------------------------------------------------
      INTERFACE unprimed_grid_class
         MODULE PROCEDURE unprimed_grid_construct_m,                           &
     &                    unprimed_grid_construct_a
      END INTERFACE

      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Construct a @ref unprimed_grid_class object.
!>
!>  Allocates memory and initializes a @ref unprimed_grid_class object. This
!>  computes the vector potential, magnetic field and positions on the unprimed
!>  grid. Magnetic field components are provided on a truncated grid. Due to the
!>  central differencing, field componetents cannot be provided for first and
!>  last points on the r and z grids.
!>
!>  @param[in] mgrid    A @ref m_grid::m_grid_class object.
!>  @param[in] pgrid    A @ref primed_grid::primed_grid_class object.
!>  @param[in] p_start  Starting phi index to compute fields to.
!>  @param[in] p_end    Ending phi index to compute fields to.
!>  @param[in] parallel @ref bmw_parallel_context::bmw_parallel_context_class
!>                      object instance.
!>  @param[in] io_unit  Unit number to write messages to.
!>  @returns A pointer to a constructed @ref unprimed_grid_class object.
!-------------------------------------------------------------------------------
      FUNCTION unprimed_grid_construct_m(mgrid, pgrid, p_start, p_end,         &
     &                                   parallel, io_unit)
!$    USE omp_lib
      USE bmw_state_flags, ONLY: clear_screen, progress
      USE, INTRINSIC :: iso_fortran_env, Only : output_unit

      IMPLICIT NONE

!  Declare Arguments
      CLASS (unprimed_grid_class), POINTER :: unprimed_grid_construct_m
      CLASS (m_grid_class), INTENT(in)               :: mgrid
      CLASS (primed_grid_class), INTENT(in)          :: pgrid
      INTEGER, INTENT(in)                            :: p_start
      INTEGER, INTENT(in)                            :: p_end
      CLASS (bmw_parallel_context_class), INTENT(in) :: parallel
      INTEGER, INTENT(in)                            :: io_unit

!  local variables
      REAL (rprec)                                   :: start_time
      INTEGER                                        :: i
      INTEGER                                        :: ri
      INTEGER                                        :: zi
      INTEGER                                        :: vi
      INTEGER                                        :: num_r
      INTEGER                                        :: num_p
      INTEGER                                        :: num_z
      REAL (rprec)                                   :: x
      REAL (rprec)                                   :: y
      REAL (rprec)                                   :: ax
      REAL (rprec)                                   :: ay
      INTEGER                                        :: k_p
      INTEGER                                        :: k_m
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: r
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: z
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: r_p
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: r_m
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: cosv
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: sinv
      REAL (rprec)                                   :: total
      REAL (rprec)                                   :: done
      REAL (rprec)                                   :: current
      REAL (rprec), DIMENSION(:,:,:), ALLOCATABLE    :: rp
      REAL (rprec)                                   :: ar_p
      REAL (rprec)                                   :: ar_m
      REAL (rprec)                                   :: ap_p
      REAL (rprec)                                   :: ap_m
      REAL (rprec)                                   :: az_p
      REAL (rprec)                                   :: az_m

!  Start of executable code
      start_time = profiler_get_start_time()

      ALLOCATE(unprimed_grid_construct_m)

      num_r = SIZE(mgrid%a_r, 1)
      num_p = SIZE(mgrid%a_p, 3)
      num_z = SIZE(mgrid%a_z, 2)

      ALLOCATE(unprimed_grid_construct_m%a_r(num_r,num_z,num_p))
      ALLOCATE(unprimed_grid_construct_m%a_p(num_r,num_z,num_p))
      ALLOCATE(unprimed_grid_construct_m%a_z(num_r,num_z,num_p))

      ALLOCATE(unprimed_grid_construct_m%b_r(num_r,num_z,num_p))
      ALLOCATE(unprimed_grid_construct_m%b_p(num_r,num_z,num_p))
      ALLOCATE(unprimed_grid_construct_m%b_z(num_r,num_z,num_p))

      ALLOCATE(cosv(num_p))
      ALLOCATE(sinv(num_p))

      ALLOCATE(r(num_r))
      ALLOCATE(z(num_z))
      ALLOCATE(r_p(num_r))
      ALLOCATE(r_m(num_r))

      total = parallel%end(num_r*num_z*num_p)
      total = CEILING(total/parallel%num_threads)

!$OMP PARALLEL
!$OMP& DEFAULT(SHARED)
!$OMP& PRIVATE(i, ri, zi, vi, x, y, ax, ay, rp, k_p, k_m,                      &
!$OMP&         ar_p, ar_m, ap_p, ap_m, az_p, az_m, current)

!  Multi process will do an all reduce so these arrays need to be initalized.
      IF (parallel%stride .gt. 1) THEN
!$OMP WORKSHARE
         unprimed_grid_construct_m%a_r = 0.0
         unprimed_grid_construct_m%a_p = 0.0
         unprimed_grid_construct_m%a_z = 0.0
         unprimed_grid_construct_m%b_r = 0.0
         unprimed_grid_construct_m%b_p = 0.0
         unprimed_grid_construct_m%b_z = 0.0
         cosv = 0.0
         sinv = 0.0
         r = 0.0
         z = 0.0
         r_p = 0.0
         r_m = 0.0
!$OMP END WORKSHARE
      END IF

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_r), parallel%end(num_r)
         r(i) = (i - 1.0)*mgrid%dr + mgrid%rmin
         r_p(i) = r(i) + mgrid%dr
         r_m(i) = r(i) - mgrid%dr
      END DO
!$OMP END DO

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_z), parallel%end(num_z)
         z(i) = (i - 1.0)*mgrid%dz + mgrid%zmin
      END DO
!$OMP END DO

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_p), parallel%end(num_p)
         x = (i - 1.0)*pgrid%dv
         cosv(i) = COS(x)
         sinv(i) = SIN(x)
      END DO
!$OMP END DO

!$OMP SINGLE
      IF (parallel%stride .gt. 1) THEN
         CALL parallel%reduce(cosv)
         CALL parallel%reduce(sinv)
         CALL parallel%reduce(r)
         CALL parallel%reduce(z)
         CALL parallel%reduce(r_p)
         CALL parallel%reduce(r_m)
      END IF
!$OMP END SINGLE

      current = 0.0

      ALLOCATE(rp(SIZE(pgrid%x, 1),                                            &
     &            SIZE(pgrid%x, 2),                                            &
     &            SIZE(pgrid%x, 3)))

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_r*num_z*num_p),                                &
     &       parallel%end(num_r*num_z*num_p)
         ri = bmw_parallel_context_i(i, num_r)
         zi = bmw_parallel_context_j(i, num_r, num_z)
         vi = bmw_parallel_context_k(i, num_r, num_z)

         IF (parallel%offset .eq. 0) THEN
!$          IF (OMP_GET_THREAD_NUM() .eq. 0) THEN
               current = current + 1.0
               done = 100.0*current/total

               WRITE (io_unit,1000,ADVANCE='NO')                               &
     &            clear_screen, progress(MOD(INT(current),4)), done

               IF (io_unit .ne. output_unit) THEN
                  BACKSPACE (io_unit)
               END IF
!$          END IF
         END IF

         x = r(ri)*cosv(vi)
         y = r(ri)*sinv(vi)

         rp = SQRT((pgrid%x - x)**2.0 + (pgrid%y - y)**2.0 +                   &
     &             (pgrid%z - z(zi))**2.0)

         ax = SUM(pgrid%j_x/rp)*pgrid%dvol
         ay = SUM(pgrid%j_y/rp)*pgrid%dvol

         unprimed_grid_construct_m%a_r(ri,zi,vi) =                             &
     &       x/r(ri)*ax + y/r(ri)*ay + mgrid%a_r(ri,zi,vi)
         unprimed_grid_construct_m%a_p(ri,zi,vi) =                             &
     &      -y/r(ri)*ax + x/r(ri)*ay + mgrid%a_p(ri,zi,vi)
         unprimed_grid_construct_m%a_z(ri,zi,vi) =                             &
     &      SUM(pgrid%j_z/rp)*pgrid%dvol + mgrid%a_z(ri,zi,vi)
      END DO
!$OMP END DO

      DEALLOCATE(rp)

!$OMP SINGLE
!  Multi process did not fill out the entire array. Get the missing pieces from
!  the other processes.
      IF (parallel%stride .gt. 1) THEN
         CALL parallel%reduce(unprimed_grid_construct_m%a_r)
         CALL parallel%reduce(unprimed_grid_construct_m%a_p)
         CALL parallel%reduce(unprimed_grid_construct_m%a_z)
      END IF
!$OMP END SINGLE

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_p*num_r*num_z),                                &
     &       parallel%end(num_p*num_r*num_z)
         ri = bmw_parallel_context_i(i, num_r)
         zi = bmw_parallel_context_j(i, num_r, num_z)
         vi = bmw_parallel_context_k(i, num_r, num_z)

         IF (ri .eq. 1 .or. ri .eq. num_r .or.                                 &
     &       zi .eq. 1 .or. zi .eq. num_z) THEN
            unprimed_grid_construct_m%b_r(ri,zi,vi) = 0.0
            unprimed_grid_construct_m%b_p(ri,zi,vi) = 0.0
            unprimed_grid_construct_m%b_z(ri,zi,vi) = 0.0
         ELSE
            IF (num_p .eq. 1) THEN
               k_p = 1
               k_m = 1
            ELSE
               IF (vi .eq. 1) THEN
                  k_p = 2
                  k_m = num_p
               ELSE IF (vi .eq. num_p) THEN
                  k_p = 1
                  k_m = num_p - 1
               ELSE
                  k_p = vi + 1
                  k_m = vi - 1
               END IF
            END IF

! 1/rdazdp - dapdz
            az_p = unprimed_grid_construct_m%a_z(ri,zi,k_p)
            az_m = unprimed_grid_construct_m%a_z(ri,zi,k_m)
            ap_p = unprimed_grid_construct_m%a_p(ri,zi + 1,vi)
            ap_m = unprimed_grid_construct_m%a_p(ri,zi - 1,vi)
            unprimed_grid_construct_m%b_r(ri,zi,vi) =                          &
     &         (az_p - az_m)/(2.0*pgrid%dv*r(ri)) -                            &
     &         (ap_p - ap_m)/(2.0*mgrid%dz)

! dardz - dazdr
            ar_p = unprimed_grid_construct_m%a_r(ri,zi + 1,vi)
            ar_m = unprimed_grid_construct_m%a_r(ri,zi - 1,vi)
            az_p = unprimed_grid_construct_m%a_z(ri + 1,zi,vi)
            az_m = unprimed_grid_construct_m%a_z(ri - 1,zi,vi)
            unprimed_grid_construct_m%b_p(ri,zi,vi) =                          &
     &         (ar_p - ar_m)/(2.0*mgrid%dz) -                                  &
     &         (az_p - az_m)/(2.0*mgrid%dr)

! 1/r(drapdr - dardp)
            ar_p = unprimed_grid_construct_m%a_r(ri,zi,k_p)
            ar_m = unprimed_grid_construct_m%a_r(ri,zi,k_m)
            ap_p = unprimed_grid_construct_m%a_p(ri + 1,zi,vi)
            ap_m = unprimed_grid_construct_m%a_p(ri - 1,zi,vi)
            unprimed_grid_construct_m%b_z(ri,zi,vi) =                          &
     &         (r_p(ri)*ap_p - r_m(ri)*ap_m)/(2.0*mgrid%dr*r(ri)) -            &
     &         (ar_p - ar_m)/(2.0*pgrid%dv*r(ri))

         END IF
      END DO
!$OMP END DO
!$OMP END PARALLEL

      IF (parallel%stride .gt. 1) THEN
         CALL parallel%reduce(unprimed_grid_construct_m%b_r)
         CALL parallel%reduce(unprimed_grid_construct_m%b_p)
         CALL parallel%reduce(unprimed_grid_construct_m%b_z)
      END IF

      DEALLOCATE(cosv)
      DEALLOCATE(sinv)
      DEALLOCATE(r)
      DEALLOCATE(z)
      DEALLOCATE(r_p)
      DEALLOCATE(r_m)

      IF (parallel%offset .eq. 0) THEN
         WRITE (io_unit,1001) clear_screen
      END IF

      CALL profiler_set_stop_time('unprimed_grid_construct_m',                 &
     &                            start_time)

1000  FORMAT(a,a,f6.2,' % Finished')
1001  FORMAT(a,'Unprimed Grid Finished')

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Construct a @ref unprimed_grid_class object at specfied points.
!>
!>  Allocates memory and initializes a @ref unprimed_grid_class object. This
!>  computes the vector potential, magnetic field and positions on the unprimed
!>  grid. The primed grid is supplied as arbitrary r and z positions with a
!>  uniformly spaced phi grid. Since the primed grid is arbitrary, this function
!>  does not compute the magnetic field components. The middle index represents
!>  phi positions.
!>
!>  @param[in] mgrid    A @ref m_grid::m_grid_class object.
!>  @param[in] pgrid    A @ref primed_grid::primed_grid_class object.
!>  @param[in] r_grid   3D Array of r positions.
!>  @param[in] z_grid   3D Array of z positions.
!>  @param[in] dphi     Phi grid spacing.
!>  @param[in] parallel @ref bmw_parallel_context::bmw_parallel_context_class
!>                      object instance.
!>  @param[in] io_unit  Unit number to write messages to.
!>  @returns A pointer to a constructed @ref unprimed_grid_class object.
!-------------------------------------------------------------------------------
      FUNCTION unprimed_grid_construct_a(mgrid, pgrid, r_grid, z_grid,         &
     &                                   dphi, parallel, io_unit)
!$    USE omp_lib
      USE bmw_state_flags, ONLY: clear_screen, progress
      USE, INTRINSIC :: iso_fortran_env, Only : output_unit

      IMPLICIT NONE

!  Declare Arguments
      CLASS (unprimed_grid_class), POINTER :: unprimed_grid_construct_a
      CLASS (m_grid_class), INTENT(in)               :: mgrid
      CLASS (primed_grid_class), INTENT(in)          :: pgrid
      REAL (rprec), DIMENSION(:,:,:), INTENT(in)     :: r_grid
      REAL (rprec), DIMENSION(:,:,:), INTENT(in)     :: z_grid
      REAL (rprec), INTENT(in)                       :: dphi
      CLASS (bmw_parallel_context_class), INTENT(in) :: parallel
      INTEGER, INTENT(in)                            :: io_unit

!  local variables
      REAL (rprec)                                   :: start_time
      INTEGER                                        :: num_1
      INTEGER                                        :: num_2
      INTEGER                                        :: num_3
      REAL (rprec)                                   :: current
      REAL (rprec)                                   :: done
      REAL (rprec)                                   :: total
      INTEGER                                        :: i
      INTEGER                                        :: ai
      INTEGER                                        :: aj
      INTEGER                                        :: ak
      REAL (rprec), DIMENSION(:,:,:), ALLOCATABLE    :: rp
      REAL (rprec)                                   :: ax
      REAL (rprec)                                   :: ay
      REAL (rprec)                                   :: x
      REAL (rprec)                                   :: y
      REAL (rprec)                                   :: z
      REAL (rprec)                                   :: r
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: phi
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: cosv
      REAL (rprec), DIMENSION(:), ALLOCATABLE        :: sinv

!  Start of executable code
      start_time = profiler_get_start_time()

      IF (MAXVAL(r_grid) .gt. mgrid%rmax .or.                                  &
     &    MINVAL(r_grid) .lt. mgrid%rmin .or.                                  &
     &    MAXVAL(z_grid) .gt. mgrid%zmax .or.                                  &
     &    MINVAL(z_grid) .lt. mgrid%zmin) THEN
         IF (parallel%offset .eq. 0) THEN
            WRITE (io_unit, 2000)
         END IF
         CALL bmw_parallel_context_abort(-1)
      END IF

      ALLOCATE(unprimed_grid_construct_a)

      num_1 = SIZE(r_grid, 1)
      num_2 = SIZE(r_grid, 2)
      num_3 = SIZE(r_grid, 3)

      ALLOCATE(unprimed_grid_construct_a%a_r(num_1,num_2,num_3))
      ALLOCATE(unprimed_grid_construct_a%a_p(num_1,num_2,num_3))
      ALLOCATE(unprimed_grid_construct_a%a_z(num_1,num_2,num_3))

      ALLOCATE(cosv(num_2))
      ALLOCATE(sinv(num_2))
      ALLOCATE(phi(num_2))

      total = parallel%end(num_1*num_2*num_3)
      total = CEILING(total/parallel%num_threads)

!$OMP PARALLEL
!$OMP& DEFAULT(SHARED)
!$OMP& PRIVATE(i, ai, aj, ak, r, x, y, z, ax, ay, current, rp)

!  Multi process will do an all reduce so these arrays need to be initalized.
      IF (parallel%stride .gt. 1) THEN
!$OMP WORKSHARE
         unprimed_grid_construct_a%a_r = 0.0
         unprimed_grid_construct_a%a_p = 0.0
         unprimed_grid_construct_a%a_z = 0.0
         cosv = 0.0
         sinv = 0.0
         phi = 0.0
!$OMP END WORKSHARE
      END IF

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_2), parallel%end(num_2)
         phi(i) = (i - 1.0)*dphi
         cosv(i) = COS(phi(i))
         sinv(i) = SIN(phi(i))
      END DO
!$OMP END DO

!$OMP SINGLE
      IF (parallel%stride .gt. 1) THEN
         CALL parallel%reduce(cosv)
         CALL parallel%reduce(sinv)
         CALL parallel%reduce(phi)
      END IF
!$OMP END SINGLE

      current = 0.0

      ALLOCATE(rp(SIZE(pgrid%x, 1),                                            &
     &            SIZE(pgrid%x, 2),                                            &
     &            SIZE(pgrid%x, 3)))

!$OMP DO
!$OMP& SCHEDULE(STATIC)
      DO i = parallel%start(num_1*num_2*num_3),                                &
     &       parallel%end(num_1*num_2*num_3)
         ai = bmw_parallel_context_i(i, num_1)
         aj = bmw_parallel_context_j(i, num_1, num_2)
         ak = bmw_parallel_context_k(i, num_1, num_2)

         IF (parallel%offset .eq. 0) THEN
!$          IF (OMP_GET_THREAD_NUM() .eq. 0) THEN
               current = current + 1.0
               done = 100.0*current/total

               WRITE (io_unit,1000,ADVANCE='NO')                               &
     &            clear_screen, progress(MOD(INT(current),4)), done

               IF (io_unit .ne. output_unit) THEN
                  BACKSPACE (io_unit)
               END IF
!$          END IF
         END IF

         r = r_grid(ai,aj,ak)
         z = z_grid(ai,aj,ak)
         x = r*cosv(aj)
         y = r*sinv(aj)

         rp = SQRT((pgrid%x - x)**2.0 +                                        &
     &             (pgrid%y - y)**2.0 +                                        &
     &             (pgrid%z - z)**2.0)

         ax = SUM(pgrid%j_x/rp)*pgrid%dvol
         ay = SUM(pgrid%j_y/rp)*pgrid%dvol

!  Interpolate from the mgrid grid to the unprimed R Z grid to add the vacuum
!  vector potential.
         CALL mgrid%interpolate(r, phi(aj), z,                                 &
     &           unprimed_grid_construct_a%a_r(ai,aj,ak),                      &
     &           unprimed_grid_construct_a%a_p(ai,aj,ak),                      &
     &           unprimed_grid_construct_a%a_z(ai,aj,ak))

         unprimed_grid_construct_a%a_r(ai,aj,ak) =                             &
     &      unprimed_grid_construct_a%a_r(ai,aj,ak) + x/r*ax + y/r*ay
         unprimed_grid_construct_a%a_p(ai,aj,ak) =                             &
     &      unprimed_grid_construct_a%a_p(ai,aj,ak) - y/r*ax + x/r*ay
         unprimed_grid_construct_a%a_z(ai,aj,ak) =                             &
     &      unprimed_grid_construct_a%a_z(ai,aj,ak) +                          &
     &      SUM(pgrid%j_z/rp)*pgrid%dvol
      END DO
!$OMP END DO

      DEALLOCATE(rp)

!$OMP END PARALLEL

!  Multi process did not fill out the entire array. Get the missing pieces from
!  the other processes.
      IF (parallel%stride .gt. 1) THEN
         CALL parallel%reduce(unprimed_grid_construct_a%a_r)
         CALL parallel%reduce(unprimed_grid_construct_a%a_p)
         CALL parallel%reduce(unprimed_grid_construct_a%a_z)
      END IF

      IF (parallel%offset .eq. 0) THEN
         WRITE (io_unit,1001) clear_screen
      END IF

      CALL profiler_set_stop_time('unprimed_grid_construct_a',                 &
     &                            start_time)

1000  FORMAT(a,a,f6.2,' % Finished')
1001  FORMAT(a,'Unprimed Grid Finished')

2000  FORMAT('Error: Unprimed grid extends beyond the mgrid grid.')

      END FUNCTION

!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Deconstruct a @ref unprimed_grid_class object.
!>
!>  Deallocates memory and uninitializes a @ref unprimed_grid_class object.
!>
!>  @param[inout] this A @ref unprimed_grid_class instance.
!-------------------------------------------------------------------------------
      SUBROUTINE unprimed_grid_destruct(this)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (unprimed_grid_class), INTENT(inout) :: this

!  Start of executable code
      IF (ASSOCIATED(this%a_r)) THEN
         DEALLOCATE(this%a_r)
         this%a_r => null()
      END IF

      IF (ASSOCIATED(this%a_p)) THEN
         DEALLOCATE(this%a_p)
         this%a_p => null()
      END IF

      IF (ASSOCIATED(this%a_z)) THEN
         DEALLOCATE(this%a_z)
         this%a_z => null()
      END IF

      IF (ASSOCIATED(this%b_r)) THEN
         DEALLOCATE(this%b_r)
         this%b_r => null()
      END IF

      IF (ASSOCIATED(this%b_p)) THEN
         DEALLOCATE(this%b_p)
         this%b_p => null()
      END IF

      IF (ASSOCIATED(this%b_z)) THEN
         DEALLOCATE(this%b_z)
         this%b_z => null()
      END IF

      END SUBROUTINE

      END MODULE

max_stars_repo_name: ORNL-Fusion/BMW

max_issues_repo_name: ORNL-Fusion/BMW

max_forks_repo_name: ORNL-Fusion/BMW


text:
!==============================================================================!
  module Flow_Mod
!------------------------------------------------------------------------------!
!   Module for basic flow field plus temperature.                              !
!   It is a bit of a mumbo-jumbo at this moment, it will furhter have to       !
!   differentiate into numerical and physica parts.                            !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Var_Mod
  use Bulk_Mod
  use Matrix_Mod
!------------------------------------------------------------------------------!
  implicit none
!==============================================================================!

  ! Velocity components
  type(Var_Type) :: u
  type(Var_Type) :: v
  type(Var_Type) :: w

  ! Temperature
  type(Var_Type) :: t

  ! Pressure
  type(Var_Type) :: p
  type(Var_Type) :: pp

  ! Mass fluxes throught cell faces
  real, allocatable :: flux(:)

  ! Variables determining if we are dealing with heat transfer and buoyancy
  logical :: heat_transfer
  logical :: buoyancy

  ! Geometrical staff
  real,allocatable :: f_coef(:)  ! face coefficient
  real,allocatable :: fw(:)      ! weight factors for the fluid phase

  ! Right hand side for velocity and pressure equations
  type(Matrix_Type) :: a  ! system matrix for all variables
  real, allocatable :: b(:)

  ! Mass fluxes, bulk velocities and pressure drops
  type(Bulk_Type) :: bulk

  ! Physical properties
  real :: viscosity, density, conductivity, diffusivity, capacity

  ! Angular velocity
  real :: omega_x, omega_y, omega_z, omega

  ! Reference temperature
  real :: t_ref, t_inf

  ! Heat and heat flux to the domain
  real :: heat, heat_flux

  ! Gravity
  real :: grav_x, grav_y, grav_z

  end module

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


module combinations_module
   use iso_fortran_env
implicit none
contains
!
!----------------------------------
!combinations - function:
!----------------------------------
!
integer function combination_x_of_p(x,p)
implicit none
!
integer(kind=int64), intent(in) :: x, p
integer(kind=int64) :: diff, prod, x_fact, p_x_fact, i
!
diff = p - x
write(*,*) 'x', x, 'p', p
write(*,*) 'diff', diff
!
!--------------------
! p-x > x :
!--------------------
!
if (diff .ge. x) then
   write(*,*) 'p - x > x'
   prod = 1
   do i = 1, x
      prod = prod * (diff + i)
      write(*,*) 'prod', prod
   end do
   !
   if (x .ne. 0) then
      x_fact = 1
      do i = 1, x
         x_fact = x_fact * i
      end do
   end if
   if (x .eq. 0) then
      x_fact = 1
   end if
   !
   write(*,*) 'prod', prod
   write(*,*) 'x_fact', x_fact
   combination_x_of_p = prod/x_fact
   write(*,*) 'combination x of p', combination_x_of_p
!
!--------------------------
! p-x < x :
!--------------------------
!
else if (diff .lt. x) then
   write(*,*) 'x > p - x'
   prod = 1
   do i = 1, diff
      prod = prod * (x+i)
   end do
   !
   p_x_fact = 1
   do i = 1, p - x
      p_x_fact = p_x_fact * i
   end do
   !
   combination_x_of_p = prod/p_x_fact
end if
!
!-----------------------------------------------------------------
! example -> x .lt. p-x
!-----------------------------------------------------------------
!
!   p = 8
!   x = 2
!   p - x = 6
!
!   x < p-x
!
!   6! 7 8            p-x+1 p-x+2
! -------------- = --------------------
!   6! 2!               x!
!
!-----------------------------------------------------------------
! example -> x .gt. p=x
!-----------------------------------------------------------------
!
!   p = 8
!   x = 5
!   p - x = 3
!
!   5! 6 7 8        x+1 x+2 x+3
! ------------ = ------------------
!   5! 3!           p - x!
!
end function combination_x_of_p

end module combinations_module


text:
module tensor_arithmetic_primitives_dev
    use :: tensor_arithmetic_primitives_api

    use :: intrinsic_tap_module, only : intrinsic_tap

    implicit none
    public
end module tensor_arithmetic_primitives_dev

max_stars_repo_name: cheshyre/ntcl-tensor

max_issues_repo_name: cheshyre/ntcl-tensor

max_forks_repo_name: cheshyre/ntcl-tensor


! { dg-do run }
! { dg-options "-fcoarray=single" }
!
! Tests if team_number intrinsic fucntion works
!
  use iso_fortran_env, only : team_type
  implicit none
  type(team_type) :: team
  integer, parameter :: standard_initial_value=-1

  associate(new_team => mod(this_image(),2)+1)
    form team (new_team,team)
      change team (team)
    end team
  end associate

  if (team_number()/=standard_initial_value) STOP 1
end


text:
subroutine do_initialization() bind(C)

  use microphysics_module

  implicit none

  character (len=32) :: probin_file
  integer :: probin_pass(32)
  integer :: i

  probin_file = "probin"
  do i = 1, len(trim(probin_file))
     probin_pass(i) = ichar(probin_file(i:i))
  enddo

  call runtime_init(probin_pass(1:len(trim(probin_file))), len(trim(probin_file)))

  call microphysics_init()

end subroutine do_initialization

max_stars_repo_name: ChrisDeGrendele/Microphysics

max_issues_repo_name: Youhichka/Microphysics

max_forks_repo_name: Youhichka/Microphysics


text:
      function dotprd (v1,v2)
c     implicit double precision (a-h,o-z)                                    dp
      common/bk03/numdc,imassn,idampn,irller,penstf
      common/bk17/dn1,dn2,nwebuf,ntime,numnp,neq,ibar,mthsol
      dimension v1(*),v2(*)
      common/main_block/ a(1)
c
      dotprd=fdot(v1,v2,neq)
      if (numdc.eq.0) return
      k73=igtpnt(73)
      k57=igtpnt(57)
      k72=igtpnt(72)
      change=0.0
      call moddot(a(k57),a(k72),a(k73),numdc,v1,v2,change)
      dotprd=dotprd-change
      return
c
      end

ext: f

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
*----------------------------------------------------------------------*
      subroutine contr_xarc_sort(contr)
*----------------------------------------------------------------------*
*     if possible, reorder vertices such that they contribute to xarcs
*     in ascending order
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'ifc_operators.h'

      integer, parameter ::
     &     ntest = 00

      type(contraction), intent(inout), target ::
     &     contr

      integer ::
     &     nvtx, narc, nxarc, nj, ivtx, iarc, hlp
      integer ::
     &     svertex_hlp(contr%nvtx)
      type(cntr_arc), pointer ::
     &     xarc(:), arc(:)
      type(cntr_vtx), pointer ::
     &     vtx_hlp(:)
      logical ::
     &     changed, swap
      integer, pointer ::
     &     svmap(:,:), xl4vtx(:), ireo(:), ireo2(:)

      if (ntest.ge.100)
     &     call write_title(lulog,wst_dbg_subr,'contr_xarc_sort')

      nxarc = contr%nxarc

      ! if there are no xarcs, we can quit
      if (nxarc.eq.0) return

      nvtx = contr%nvtx

      ! only one vertex: nothing to do
      if (nvtx.le.1) return

      nj = 0
      xarc => contr%xarc
      do iarc = 1, nxarc
        nj = max(nj,xarc(iarc)%link(2))
      end do

      ! also this is an uninteresting case
      if (nj.le.1) return

      narc = contr%narc
      arc => contr%arc

      allocate(xl4vtx(nvtx),svmap(nvtx,nj),ireo(nvtx),ireo2(nvtx))

      xl4vtx(1:nvtx) = 0
      svmap(1:nvtx,1:nj) = 0
      do iarc = 1, nxarc
        ivtx = xarc(iarc)%link(1)
        xl4vtx(ivtx) = xl4vtx(ivtx)+1
        svmap(ivtx,xl4vtx(ivtx)) = xarc(iarc)%link(2)
      end do

      do ivtx = 1, nvtx
        if (xl4vtx(ivtx).gt.2)
     &   call quit(1,'contr_xarc_sort','xl4vtx(ivtx).gt.2: Test this!')
        if (xl4vtx(ivtx).eq.2) then
          if (svmap(ivtx,2).le.svmap(ivtx,1))
     &   call quit(1,'contr_xarc_sort','not ascending svmap!')
        end if
      end do

      do ivtx = 1, nvtx
        ireo(ivtx) = ivtx
      end do

      changed = .true.
      do while (changed)
!     dbg
!        write(lulog,*) 'current ireo: ',ireo
!     dbg

        changed = .false.
        do ivtx = 2, nvtx
          if (svmap(ireo(ivtx-1),1).eq.svmap(ireo(ivtx),1)
     &         .and.svmap(ireo(ivtx),1)==0) cycle
          swap = svmap(ireo(ivtx-1),1).gt.svmap(ireo(ivtx),1)
          if(.not.swap.and.svmap(ireo(ivtx-1),1).eq.svmap(ireo(ivtx),1))
     &         then
            swap = svmap(ireo(ivtx-1),2).gt.svmap(ireo(ivtx),2)
          end if
!     dbg
!          if (swap) write(lulog,*) 'request swap of ',
!     &         ireo(ivtx-1),ireo(ivtx)
!     dbg
          if (swap.and.may_commute(ireo(ivtx-1),ireo(ivtx))) then
            changed = .true.
!     dbg
!            write(lulog,*) 'request accepted'
!     dbg

            hlp = ireo(ivtx)
            ireo(ivtx)=ireo(ivtx-1)
            ireo(ivtx-1) = hlp
          end if
        end do
      end do

      ! note: ireo(idx) = oldidx
      !   set ireo2(oldidx) = idx
      do ivtx = 1, nvtx
        ireo2(ireo(ivtx)) = ivtx
      end do
!     dbg
!      write(lulog,*) 'final reo:  ',ireo
!      write(lulog,*) 'final reo2: ',ireo2
!     dbg

      allocate(vtx_hlp(nvtx))
      vtx_hlp = contr%vertex

      do ivtx = 1, nvtx
        contr%vertex(ivtx) = vtx_hlp(ireo(ivtx))
      end do

      deallocate(vtx_hlp)

      svertex_hlp = contr%svertex
      do ivtx = 1, nvtx
        contr%svertex(ivtx) = svertex_hlp(ireo(ivtx))
      end do

      do iarc = 1, narc
        contr%arc(iarc)%link(1) = ireo2(contr%arc(iarc)%link(1))
        contr%arc(iarc)%link(2) = ireo2(contr%arc(iarc)%link(2))
      end do

      do iarc = 1, nxarc
        contr%xarc(iarc)%link(1) = ireo2(contr%xarc(iarc)%link(1))
      end do

      call update_svtx4contr(contr)

      deallocate(svmap,xl4vtx,ireo,ireo2)

      return

      contains

      logical function may_commute(ivtx1,ivtx2)

      implicit none
      integer ivtx1, ivtx2

      integer iarc

      may_commute = .true.
      do iarc = 1, narc
!     dbg
!        write(lulog,*) 'checking: ',iarc
!     dbg

        if (arc(iarc)%link(1)==ivtx1 .and.
     &       arc(iarc)%link(2)==ivtx2) then
          may_commute = .false.
          exit
        end if
      end do

      return

      end function

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!##############################################################################
!# ****************************************************************************
!# <name> triatest3d </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module reads in a 3D mesh, optionally refines it a few times and
!# creates a standard mesh from it. Afterwards (nearly) all arrays from the
!# triangulation are printed onto the screen and to the log file.
!# </purpose>
!##############################################################################
module triatest3d

  use fsystem
  use genoutput
  use boundary
  use triangulation
  use storage

  implicit none

contains

  subroutine tria_test3d(NLMAX)
  integer, optional, intent(IN) :: NLMAX

  type(t_triangulation) :: rtria
  integer, dimension(:), pointer :: p_Idata,p_Iidx
  integer, dimension(:,:), pointer :: p_Idata2
  real(DP), dimension(:), pointer :: p_Ddata
  real(DP), dimension(:,:), pointer :: p_Ddata2
  integer :: i,j

    ! Print a header
    call output_separator(OU_SEP_STAR)
    call output_line('T R I A _ T E S T 3 D')
    call output_separator(OU_SEP_STAR)

    ! Now read in the basic triangulation.
    call tria_readTriFile3D (rtria, './pre/CUBE.tri')

    ! Optionally refine the mesh
    if(present(NLMAX)) then
      call tria_quickRefine2LevelOrdering (NLMAX-1,rtria)
    end if

    ! Standardise the mesh
    call tria_initStandardMeshFromRaw (rtria)

    ! First of all, print out the basic information
    call output_line('Triangulation Scalars')
    call output_line('----------------------------------')
    call output_line('NVT....................: ' // trim(sys_si(rtria%NVT,8)))
    call output_line('NMT....................: ' // trim(sys_si(rtria%NMT,8)))
    call output_line('NAT....................: ' // trim(sys_si(rtria%NAT,8)))
    call output_line('NEL....................: ' // trim(sys_si(rtria%NEL,8)))
    call output_line('NBCT...................: ' // trim(sys_si(rtria%NBCT,8)))
    call output_line('NblindBCT..............: ' // trim(sys_si(rtria%NblindBCT,8)))
    call output_line('NVBD...................: ' // trim(sys_si(rtria%NVBD,8)))
    call output_line('NMBD...................: ' // trim(sys_si(rtria%NMBD,8)))
    call output_line('NABD...................: ' // trim(sys_si(rtria%NABD,8)))
    call output_line('NNVE...................: ' // trim(sys_si(rtria%NNVE,8)))
    call output_line('NNEE...................: ' // trim(sys_si(rtria%NNEE,8)))
    call output_line('NNAE...................: ' // trim(sys_si(rtria%NNAE,8)))
    call output_line('NNelAtVertex...........: ' // trim(sys_si(rtria%NNelAtVertex,8)))
    call output_line('NNelAtEdge.............: ' // trim(sys_si(rtria%NNelAtEdge,8)))
    call output_line('nverticesPerEdge.......: ' // trim(sys_si(rtria%nverticesPerEdge,8)))
    call output_line('nverticesOnAllEdges....: ' // trim(sys_si(rtria%nVerticesOnAllEdges,8)))
    call output_line('nverticesInEachElement.: ' // trim(sys_si(rtria%nverticesInEachElement,8)))
    call output_line('nverticesInAllElements.: ' // trim(sys_si(rtria%nverticesInAllElements,8)))
    call output_line('nadditionalVertices....: ' // trim(sys_si(rtria%NMT,8)))

    ! Print the vertice coordinates
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_double2d(rtria%h_DvertexCoords, p_Ddata2)
    call output_line('p_DvertexCoords')
    do i = lbound(p_Ddata2,2), ubound(p_Ddata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ' // trim(sys_sdEP(p_Ddata2(1,i),20,12)) //&
              ' , ' // trim(sys_sdEP(p_Ddata2(2,i),20,12)) // &
              ' , ' // trim(sys_sdEP(p_Ddata2(3,i),20,12)) // ' /)')
    end do

    ! Print vertices-at-Element
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IverticesAtElement,p_Idata2)
    call output_line('p_IverticesAtElement')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print edges-at-Element
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IedgesAtElement,p_Idata2)
    call output_line('p_IedgesAtElement')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print faces-at-Element
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IfacesAtElement,p_Idata2)
    call output_line('h_IfacesAtElement')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print neighbours-at-Element
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IneighboursAtElement,p_Idata2)
    call output_line('h_IneighboursAtElement')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print vertices-at-edge
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IverticesAtEdge,p_Idata2)
    call output_line('h_IverticesAtEdge')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print faces at edge
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IfacesAtEdgeIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IfacesAtEdgeIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IfacesAtEdge,p_Idata)
      call output_line('h_IfacesAtEdgeIdx / h_IfacesAtEdge')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IfacesAtEdgeIdx = ST_NOHANDLE')
    end if

    ! Print elements at edge
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IelementsAtEdgeIdx3d .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IelementsAtEdgeIdx3d,p_Iidx)
      call storage_getbase_int(rtria%h_IelementsAtEdge3d,p_Idata)
      call output_line('h_IelementsAtEdgeIdx3d / h_IelementsAtEdge3d')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IelementsAtEdgeIdx3d = ST_NOHANDLE')
    end if

    ! Print vertices at face
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IverticesAtFace,p_Idata2)
    call output_line('h_IverticesAtFace')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print edges at face
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IedgesAtFace,p_Idata2)
    call output_line('h_IedgesAtFace')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print elements at face
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int2d(rtria%h_IelementsAtFace,p_Idata2)
    call output_line('h_IelementsAtFace')
    do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
      call output_line(trim(sys_si(i,8)) // ': (/ ', bnolinebreak = .true.)
      do j = lbound(p_Idata2,1), ubound(p_Idata2,1)-1
        call output_line(trim(sys_si(p_Idata2(j,i),8)) // ' , ', bnolinebreak = .true.)
      end do
      call output_line(trim(sys_si(p_Idata2(ubound(p_Idata2,1),i),8)) // ' /)')
    end do

    ! Print edges at vertices
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IedgesAtVertexIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IedgesAtVertexIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IedgesAtVertex,p_Idata)
      call output_line('h_IedgesAtVertexIdx / h_IedgesAtVertex')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IedgesAtVertexIdx = ST_NOHANDLE')
    end if

    ! Print elements at vertex
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int(rtria%h_IelementsAtVertexIdx,p_Iidx)
    call storage_getbase_int(rtria%h_IelementsAtVertex,p_Idata)
    call output_line('h_IelementsAtVertexIdx / h_IelementsAtVertex')
    do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
      call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                       ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
      do j = p_Iidx(i)+1, p_Iidx(i+1)-1
        call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
      end do
    end do

    ! Print nodal property
    call output_separator(OU_SEP_MINUS)
    call storage_getbase_int(rtria%h_InodalProperty,p_Idata)
    call output_line('h_InodalProperty')
    do i = lbound(p_Idata,1), ubound(p_Idata,1)
      call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_si(p_Idata(i),8)))
    end do

    ! Print maco nodal property
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_ImacroNodalProperty .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_ImacroNodalProperty,p_Idata)
      call output_line('h_ImacroNodalProperty')
      do i = lbound(p_Idata,1), ubound(p_Idata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_si(p_Idata(i),8)))
      end do
    else
      call output_line('h_ImacroNodalProperty = ST_NOHANDLE')
    end if

    ! Print element volume
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_DelementVolume .ne. ST_NOHANDLE) then
      call storage_getbase_double(rtria%h_DelementVolume,p_Ddata)
      call output_line('h_DelementVolume')
      do i = lbound(p_Ddata,1), ubound(p_Ddata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_sdEP(p_Ddata(i),20,12)))
      end do
    else
      call output_line('h_DelementVolume = ST_NOHANDLE')
    end if

    ! Print refinement patch
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IrefinementPatchIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IrefinementPatchIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IrefinementPatch,p_Idata)
      call output_line('h_IrefinementPatchIdx / h_IrefinementPatch')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IrefinementPatchIdx = ST_NOHANDLE')
    end if

    ! Print coarse grid element
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IcoarseGridElement .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IcoarseGridElement,p_Idata)
      call output_line('h_IcoarseGridElement')
      do i = lbound(p_Idata,1), ubound(p_Idata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_si(p_Idata(i),8)))
      end do
    else
      call output_line('h_IcoarseGridElement = ST_NOHANDLE')
    end if

    ! Print boundary vertices
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IboundaryCpIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IboundaryCpIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IverticesAtBoundary,p_Idata)
      call output_line('h_IboundaryCpIdx / h_IverticesAtBoundary')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IboundaryCpIdx = ST_NOHANDLE')
    end if

    ! Print boundary edges
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IboundaryCpEdgesIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IboundaryCpEdgesIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IedgesAtBoundary,p_Idata)
      call output_line('h_IboundaryCpEdgesIdx / h_IedgesAtBoundary')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IboundaryCpEdgesIdx = ST_NOHANDLE')
    end if

    ! Print boundary faces
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IboundaryCpFacesIdx .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IboundaryCpFacesIdx,p_Iidx)
      call storage_getbase_int(rtria%h_IfacesAtBoundary,p_Idata)
      call output_line('h_IboundaryCpFacesIdx / h_IfacesAtBoundary')
      do i = lbound(p_Iidx,1), ubound(p_Iidx,1)-1
        call output_line(trim(sys_si(i,8)) // ' , '// trim(sys_si(p_Iidx(i),8)) //&
                         ' : ' // trim(sys_si(p_Idata(p_Iidx(i)),8)))
        do j = p_Iidx(i)+1, p_Iidx(i+1)-1
          call output_line('           ' // trim(sys_si(j,8)) // ' : ' // trim(sys_si(p_Idata(j),8)))
        end do
      end do
    else
      call output_line('h_IboundaryCpFacesIdx = ST_NOHANDLE')
    end if

    ! Print boundary elements
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IelementsAtBoundary .ne. ST_NOHANDLE) then
      call storage_getbase_int(rtria%h_IelementsAtBoundary,p_Idata)
      call output_line('h_IelementsAtBoundary')
      do i = lbound(p_Idata,1), ubound(p_Idata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_si(p_Idata(i),8)))
      end do
    else
      call output_line('h_IelementsAtBoundary = ST_NOHANDLE')
    end if

    ! Print vertex parameter values
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_DvertexParameterValue .ne. ST_NOHANDLE) then
      call storage_getbase_double(rtria%h_DvertexParameterValue,p_Ddata)
      call output_line('h_DvertexParameterValue')
      do i = lbound(p_Ddata,1), ubound(p_Ddata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_sdEP(p_Ddata(i),20,12)))
      end do
    else
      call output_line('h_DvertexParameterValue = ST_NOHANDLE')
    end if

    ! Print edge parameter values
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_DedgeParameterValue .ne. ST_NOHANDLE) then
      call storage_getbase_double(rtria%h_DedgeParameterValue,p_Ddata)
      call output_line('h_DedgeParameterValue')
      do i = lbound(p_Ddata,1), ubound(p_Ddata,1)
        call output_line(trim(sys_si(i,8)) // ': ' // trim(sys_sdEP(p_Ddata(i),20,12)))
      end do
    else
      call output_line('h_DedgeParameterValue = ST_NOHANDLE')
    end if

    ! Print boundary vertex positions
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IboundaryVertexPos .ne. ST_NOHANDLE) then
      call storage_getbase_int2d(rtria%h_IboundaryVertexPos, p_Idata2)
      call output_line('h_IboundaryVertexPos')
      do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
        call output_line(trim(sys_si(i,8)) // ': (/ ' // trim(sys_si(p_Idata2(1,i),8)) //&
                ' , ' // trim(sys_si(p_Idata2(2,i),8)) // ' /)')
      end do
    else
      call output_line('h_IboundaryVertexPos = ST_NOHANDLE')
    end if

    ! Print boundary edge positions
    call output_separator(OU_SEP_MINUS)
    if (rtria%h_IboundaryEdgePos .ne. ST_NOHANDLE) then
      call storage_getbase_int2d(rtria%h_IboundaryEdgePos, p_Idata2)
      call output_line('h_IboundaryEdgePos')
      do i = lbound(p_Idata2,2), ubound(p_Idata2,2)
        call output_line(trim(sys_si(i,8)) // ': (/ ' // trim(sys_si(p_Idata2(1,i),8)) //&
                ' , ' // trim(sys_si(p_Idata2(2,i),8)) // ' /)')
      end do
    else
      call output_line('h_IboundaryEdgePos = ST_NOHANDLE')
    end if

    ! Release the mesh
    call tria_done(rtria)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
C { dg-do run }

      INCLUDE "omp_lib.h"

      DOUBLE PRECISION :: D, E
      LOGICAL :: L
      INTEGER (KIND = OMP_LOCK_KIND) :: LCK
      INTEGER (KIND = OMP_NEST_LOCK_KIND) :: NLCK

      D = OMP_GET_WTIME ()

      CALL OMP_INIT_LOCK (LCK)
      CALL OMP_SET_LOCK (LCK)
      IF (OMP_TEST_LOCK (LCK)) CALL ABORT
      CALL OMP_UNSET_LOCK (LCK)
      IF (.NOT. OMP_TEST_LOCK (LCK)) CALL ABORT
      IF (OMP_TEST_LOCK (LCK)) CALL ABORT
      CALL OMP_UNSET_LOCK (LCK)
      CALL OMP_DESTROY_LOCK (LCK)

      CALL OMP_INIT_NEST_LOCK (NLCK)
      IF (OMP_TEST_NEST_LOCK (NLCK) .NE. 1) CALL ABORT
      CALL OMP_SET_NEST_LOCK (NLCK)
      IF (OMP_TEST_NEST_LOCK (NLCK) .NE. 3) CALL ABORT
      CALL OMP_UNSET_NEST_LOCK (NLCK)
      CALL OMP_UNSET_NEST_LOCK (NLCK)
      IF (OMP_TEST_NEST_LOCK (NLCK) .NE. 2) CALL ABORT
      CALL OMP_UNSET_NEST_LOCK (NLCK)
      CALL OMP_UNSET_NEST_LOCK (NLCK)
      CALL OMP_DESTROY_NEST_LOCK (NLCK)

      CALL OMP_SET_DYNAMIC (.TRUE.)
      IF (.NOT. OMP_GET_DYNAMIC ()) CALL ABORT
      CALL OMP_SET_DYNAMIC (.FALSE.)
      IF (OMP_GET_DYNAMIC ()) CALL ABORT

      CALL OMP_SET_NESTED (.TRUE.)
      IF (.NOT. OMP_GET_NESTED ()) CALL ABORT
      CALL OMP_SET_NESTED (.FALSE.)
      IF (OMP_GET_NESTED ()) CALL ABORT

      CALL OMP_SET_NUM_THREADS (5)
      IF (OMP_GET_NUM_THREADS () .NE. 1) CALL ABORT
      IF (OMP_GET_MAX_THREADS () .NE. 5) CALL ABORT
      IF (OMP_GET_THREAD_NUM () .NE. 0) CALL ABORT
      CALL OMP_SET_NUM_THREADS (3)
      IF (OMP_GET_NUM_THREADS () .NE. 1) CALL ABORT
      IF (OMP_GET_MAX_THREADS () .NE. 3) CALL ABORT
      IF (OMP_GET_THREAD_NUM () .NE. 0) CALL ABORT
      L = .FALSE.
C$OMP PARALLEL REDUCTION (.OR.:L)
      L = OMP_GET_NUM_THREADS () .NE. 3
      L = L .OR. (OMP_GET_THREAD_NUM () .LT. 0)
      L = L .OR. (OMP_GET_THREAD_NUM () .GE. 3)
C$OMP MASTER
      L = L .OR. (OMP_GET_THREAD_NUM () .NE. 0)
C$OMP END MASTER
C$OMP END PARALLEL
      IF (L) CALL ABORT

      IF (OMP_GET_NUM_PROCS () .LE. 0) CALL ABORT
      IF (OMP_IN_PARALLEL ()) CALL ABORT
C$OMP PARALLEL REDUCTION (.OR.:L)
      L = .NOT. OMP_IN_PARALLEL ()
C$OMP END PARALLEL
C$OMP PARALLEL REDUCTION (.OR.:L) IF (.TRUE.)
      L = .NOT. OMP_IN_PARALLEL ()
C$OMP END PARALLEL

      E = OMP_GET_WTIME ()
      IF (D .GT. E) CALL ABORT
      D = OMP_GET_WTICK ()
C Negative precision is definitely wrong,
C bigger than 1s clock resolution is also strange
      IF (D .LE. 0 .OR. D .GT. 1.) CALL ABORT
      END


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis

       function zenith(lat,lon,slat,slon)

c      computes zenith angle give earth lat,lon and satellite subpoint
c      slat, slon.

c      returns zenith angle in degrees for use in SSEC routines

c      input coordinates are all in radians.  E longitude =+

c      authored by Dan Birkenheuer   12/1/95

        real lat,lon,slat,slon
        real zenith
        real cos_of_zen

        cos_of_zen = sin (slat)*sin(lat)+cos(slat)*cos(lat)*
     1      cos(abs(slon-lon))

        zenith = acos(cos_of_zen)

        zenith = zenith + atan2 ( sin(zenith),6.6166-cos_of_zen )

c       convert to degrees

        zenith = zenith * 180./acos(-1.0)

        return

      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
MODULE hsort_module

  implicit none

  PRIVATE
  PUBLIC :: indexx

CONTAINS

  SUBROUTINE indexx( n, arr, indx )
    implicit none
    integer,intent(IN)  :: n
    real(8),intent(IN)  :: arr(n)
    integer,intent(OUT) :: indx(n)
    call indexx_0( n, arr, indx )
    !call indexx_1( n, arr, indx )
    !call ascending_sort_indexx_2( n, arr, indx )
  END SUBROUTINE indexx

  SUBROUTINE indexx_0(n,arr,indx)
    implicit none
    integer,intent(IN)  :: n
    real(8),intent(IN)  :: arr(n)
    integer,intent(OUT) :: indx(n)
    integer,parameter :: m=7,nstck=50
    integer :: i,j,k,ll,jstck,jndx,lndx,ir,istck(nstck),loop0,loop1
    real(8) :: arrj,arrl

    do i=1,n
       indx(i)=i
    end do
    jstck=0
    ll=1
    ir=n
    do loop0=1,n
    if ( ir-ll < m ) then
       do j=ll+1,ir
          jndx=indx(j)
          arrj=arr(jndx)
          do i=j-1,1,-1
             if ( arr(indx(i)) <= arrj ) exit
             indx(i+1)=indx(i)
          end do
          if ( i < 1 ) i=0
          indx(i+1)=jndx
       end do
       if ( jstck == 0 ) return
       ir=istck(jstck)
       ll=istck(jstck-1)
       jstck=jstck-2
    else
       k=(ll+ir)/2
       call iswap( indx(k), indx(ll+1) )
       if ( arr(indx(ll+1)) > arr(indx(ir)) ) then
          call iswap( indx(ll+1), indx(ir) )
       end if
       if ( arr(indx(ll)) > arr(indx(ir)) ) then
          call iswap( indx(ll), indx(ir) )
       end if
       if ( arr(indx(ll+1)) > arr(indx(ll)) ) then
          call iswap( indx(ll+1), indx(ll) )
       end if
       i=ll+1
       j=ir
       lndx=indx(ll)
       arrl=arr(lndx)
       do loop1=1,n
          i=i+1
          if ( arr(indx(i)) < arrl ) cycle
          j=j-1
          do while ( arr(indx(j)) > arrl )
             j=j-1
          end do
          if ( j < i ) exit
          call iswap( indx(i), indx(j) )
       end do ! loop1
       indx(ll)=indx(j)
       indx(j )=lndx
       jstck=jstck+2
       if ( jstck > nstck ) call stop_program( "stop@indexx_0" )
       if ( ir-i+1 >= j-ll ) then
          istck(jstck)  =ir
          istck(jstck-1)=i
          ir=j-1
       else
          istck(jstck)  =j-1
          istck(jstck-1)=ll
          ll=i
       end if
    end if
    end do ! loop0
  END SUBROUTINE indexx_0

  SUBROUTINE indexx_1( ndat, arry, indx )
    implicit none
    integer,intent(IN)  :: ndat
    real(8),intent(IN)  :: arry(ndat)
    integer,intent(OUT) :: indx(ndat)
    integer :: ilayer,iboss,jlayer,jboss,kboss
    integer :: i,i1,i2,max_layer,loop,ndat_max
    integer,allocatable :: irslt(:)
    logical :: flag_boss_change
    real(8) :: arry_max

    arry_max = maxval( arry )

    max_layer = log( dble(ndat+1) )/log( 2.0d0 )

    ndat_max = 2**max_layer - 1
    if ( ndat_max < ndat ) max_layer = max_layer + 1
    ndat_max = 2**max_layer - 1
    if ( ndat_max < ndat ) stop "stop(1)"

    do i=1,ndat
       indx(i)=i
    end do

! ---

    do ilayer=2,max_layer

       loop_boss : do iboss=2**(ilayer-2),2**(ilayer-1)-1

          if ( iboss > ndat ) exit loop_boss

          flag_boss_change = .false.

          do i=2*iboss,2*iboss+1

             if ( i > ndat ) exit loop_boss

             if ( arry(indx(i)) < arry(indx(iboss)) ) then
                call iswap( indx(iboss), indx(i) )
                flag_boss_change = .true.
             end if

             if ( flag_boss_change ) then
                kboss = iboss
                do jlayer=ilayer-1,2,-1
                   do jboss=2**(jlayer-2),2**(jlayer-1)-1
                      if ( arry(indx(kboss)) < arry(indx(jboss)) ) then
                         call iswap( indx(jboss), indx(kboss) )
                         kboss = jboss
                         exit
                      end if
                   end do ! jboss
                end do ! jlayer
                if ( arry(indx(kboss)) < arry(indx(1)) ) then
                   call iswap( indx(1), indx(kboss) )
                end if
             end if

          end do ! i

       end do loop_boss

    end do ! ilayer

! ---

    allocate( irslt(ndat) ) ; irslt=0

    do loop=1,ndat

       irslt(loop)=indx(1)
       indx(1)=-1

       jboss=1
       do ilayer=2,max_layer

          i1=2*jboss
          i2=2*jboss+1

          if ( i1 > ndat ) then
             exit
          else if ( i1 <= ndat .and. i2 > ndat ) then
             call iswap( indx(jboss), indx(i1) )
             exit
          end if

          if ( indx(i1) < 0 .and. indx(i2) < 0 ) then
             exit
          else if ( indx(i1) > 0 .and. indx(i2) < 0 ) then
             call iswap( indx(jboss), indx(i1) )
             jboss=i1
          else if ( indx(i1) < 0 .and. indx(i2) > 0 ) then
             call iswap( indx(jboss), indx(i2) )
             jboss=i2
          else if ( arry(indx(i1)) < arry(indx(i2)) ) then
             call iswap( indx(jboss), indx(i1) )
             jboss=i1
          else
             call iswap( indx(jboss), indx(i2) )
             jboss=i2
          end if

       end do ! ilayer

    end do ! loop

    indx(:)=irslt(:)

    deallocate( irslt )

  END SUBROUTINE indexx_1

  SUBROUTINE iswap( i, j )
    implicit none
    integer,intent(INOUT) :: i,j
    integer :: k
    k=i
    i=j
    j=k
  END SUBROUTINE iswap

  SUBROUTINE ascending_sort_indexx_2( ndat, arry, indx )
    implicit none
    integer,intent(IN)  :: ndat
    real(8),intent(IN)  :: arry(ndat)
    integer,intent(OUT) :: indx(ndat)
    integer :: ipos, jpos,i,j,ndati,ndatj,loop,ipivo
    real(8) :: apivo

    do i=1,ndat
       indx(i) = i
    end do

    if ( ndat <= 1 ) then
       return
    else if ( ndat == 2 ) then
       if ( arry(2) < arry(1) ) then
          indx(1) = 2
          indx(2) = 1
       end if
       return
    end if

    apivo = arry(1)
    if ( apivo < arry(2) ) apivo = arry(2)

    ipos = 1
    jpos = ndat

! ---

    do loop=1,ndat

    if ( ipos == jpos ) then
       if ( arry(indx(ipos)) >= apivo ) then
          ipos = ipos - 1
       else if ( arry(indx(jpos)) < apivo ) then
          jpos = jpos + 1
       end if
       exit
    end if

    do i=ipos,jpos-1
       if ( arry(indx(i)) >= apivo ) then
          ipos = i
          exit
       end if
    end do

    do j=jpos,ipos+1,-1
       if ( arry(indx(j)) < apivo ) then
          jpos = j
          exit
       end if
    end do

    if ( i /= ipos .and. j == jpos ) then
       ipos = jpos
       jpos = jpos + 1
       exit
    else if ( i == ipos .and. j /= jpos ) then
       ipos = ipos - 1
       jpos = ipos + 1
       exit
    else if ( i /= ipos .and. j /= jpos ) then
       jpos = ipos + 1
       exit
    end if

    if ( arry(indx(ipos)) > arry(indx(jpos)) ) then
       i=indx(ipos)
       indx(ipos)=indx(jpos)
       indx(jpos)=i
    end if

    if ( ipos+1 == jpos ) exit

    ipos = ipos + 1
    jpos = jpos - 1

    end do

    if ( ipos+1 == jpos ) then
       ndati = ipos
       ndatj = ndat-jpos+1
       ipivo = 0
       call indexx_2_sub( ipivo, ndat, ndati, arry, indx )
       call indexx_2_sub( ipivo, ndat, ndatj, arry, indx(jpos) )
    else
       stop "error!:stop@ascending_sort_indexx_2"
    end if

  END SUBROUTINE ascending_sort_indexx_2

  RECURSIVE SUBROUTINE indexx_2_sub( ipivo_in, ndat, ndat1, arry, indx )
    implicit none
    integer,intent(IN) :: ipivo_in, ndat, ndat1
    real(8),intent(IN) :: arry(ndat)
    integer,intent(INOUT) :: indx(ndat1)
    integer :: i,j,ipos,jpos,ndati,ndatj,loop,ipivo
    real(8) :: apivo,apivo0,apivo1,err

    if ( ndat1 <= 1 ) then
       return
    else if ( ndat1 == 2 ) then
       if ( arry(indx(2)) < arry(indx(1)) ) then
          i=indx(1)
          indx(1) = indx(2)
          indx(2) = i
       end if
       return
    end if

    if ( ipivo_in == 0 ) then
       apivo = arry(indx(1))
       if ( apivo < arry(indx(2)) ) apivo = arry(indx(2))
    else if ( ipivo_in == 1 ) then
       apivo0 = 1.d100
       apivo1 =-1.d100
       do i=1,ndat1
          apivo = arry(indx(i))
          apivo0 = min( apivo, apivo0 )
          apivo1 = max( apivo, apivo1 )
       end do
       apivo = 0.5d0*( apivo0 + apivo1 )
    else
       stop "error!: stop@indexx_sub_2(1)"
    end if

    ipos = 1
    jpos = ndat1

! ---

    do loop=1,ndat1

    if ( ipos == jpos ) then
       if ( arry(indx(ipos)) >= apivo ) then
          ipos = ipos - 1
       else if ( arry(indx(jpos)) < apivo ) then
          jpos = jpos + 1
       end if
       exit
    end if

    do i=ipos,jpos-1
       if ( arry(indx(i)) >= apivo ) then
          ipos = i
          exit
       end if
    end do

    do j=jpos,ipos+1,-1
       if ( arry(indx(j)) < apivo ) then
          jpos = j
          exit
       end if
    end do

    if ( i /= ipos .and. j == jpos ) then
       ipos = jpos
       jpos = jpos + 1
       exit
    else if ( i == ipos .and. j /= jpos ) then
       ipos = ipos - 1
       jpos = ipos + 1
       exit
    else if ( i /= ipos .and. j /= jpos ) then
       jpos = ipos + 1
       exit
    end if

    if ( arry(indx(ipos)) > arry(indx(jpos)) ) then
       i=indx(ipos)
       indx(ipos)=indx(jpos)
       indx(jpos)=i
    end if

    if ( ipos+1 == jpos ) exit

    ipos = ipos + 1
    jpos = jpos - 1

    end do

    if ( ipos+1 == jpos ) then

       ndati = ipos
       ndatj = ndat1-jpos+1
       ipivo = 0

       if ( ipos == 0 .or. jpos == 0 ) then
          do i=1,ndat1
             err=abs(arry(indx(i))-apivo)
             if ( err > 1.d-12 ) exit
          end do
          if ( i > ndat1 ) return
          ipivo=1
       end if

       if ( ipos > 1 ) then
          call indexx_2_sub( ipivo,ndat,ndati,arry,indx )
       end if

       if ( jpos >= 1 .and. jpos < ndat1 ) then
          call indexx_2_sub( ipivo,ndat,ndatj,arry,indx(jpos) )
       end if

    else

       stop "error!: stop@indexx_2_sub"

    end if

  END SUBROUTINE indexx_2_sub

END MODULE hsort_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
 MODULE GDSWZD_EQUID_CYLIND_MOD
!$$$  MODULE DOCUMENTATION BLOCK
!
! $Revision: 74917 $
!
! MODULE:  GDSWZD_EQUID_CYLIND_MOD  GDS WIZARD MODULE FOR EQUIDISTANT
!                                   CYLINDRICAL GRIDS
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: - CONVERT FROM EARTH TO GRID COORDINATES OR VICE VERSA.
!           - COMPUTE VECTOR ROTATION SINES AND COSINES.
!           - COMPUTE MAP JACOBIANS.
!           - COMPUTE GRID BOX AREA.
!
! PROGRAM HISTORY LOG:
!   2015-01-21  GAYNO   INITIAL VERSION FROM A MERGER OF
!                       ROUTINES GDSWIZ00 AND GDSWZD00.
!   2015-09-17  GAYNO   RENAME MODULE TO "GDSWZD_EQUID_CYLIND_MOD"
!
! USAGE:  "USE GDSWZD_EQUID_CYLIND_MOD"  THEN CALL THE PUBLIC DRIVER
!         ROUTINE "GDSWZD_EQUID_CYLIND"
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
!
 IMPLICIT NONE

 PRIVATE

 PUBLIC                  :: GDSWZD_EQUID_CYLIND

 REAL                    :: DLAT ! GRID RESOLUTION IN DEGREES N/S DIRECTION
 REAL                    :: DLON ! GRID RESOLUTION IN DEGREES E/W DIRECTION
 REAL                    :: RERTH

 CONTAINS

 SUBROUTINE GDSWZD_EQUID_CYLIND(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
                                XPTS,YPTS,RLON,RLAT,NRET,  &
                                CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  GDSWZD_EQUID_CYLIND   GDS WIZARD FOR EQUIDISTANT CYLINDRICAL
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB 2 GRID DEFINITION
!           TEMPLATE (PASSED IN INTEGER FORM AS DECODED BY THE
!           NCEP G2 LIBRARY) AND RETURNS ONE OF THE FOLLOWING:
!             (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
!             (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
!           WORKS FOR EQUIDISTANT CYLINDRICAL PROJECTIONS.
!           IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
!           BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
!           OUTPUT ELEMENTS ARE SET TO FILL VALUES.  THE ACTUAL
!           NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.
!           OPTIONALLY, THE VECTOR ROTATIONS, MAP JACOBIANS AND
!           GRID BOX AREAS MAY BE RETURNED.  TO COMPUTE THE VECTOR
!           ROTATIONS, THE OPTIONAL ARGUMENTS 'SROT' AND 'CROT'
!           MUST BE PRESENT.  TO COMPUTE THE MAP JACOBIANS, THE
!           OPTIONAL ARGUMENTS 'XLON', 'XLAT', 'YLON', 'YLAT' MUST
!           BE PRESENT. TO COMPUTE THE GRID BOX AREAS, THE OPTIONAL
!           ARGUMENT 'AREA' MUST BE PRESENT.
!
! PROGRAM HISTORY LOG:
!   96-04-10  IREDELL
!   97-10-20  IREDELL  INCLUDE MAP OPTIONS
! 2015-01-21  GAYNO    MERGER OF GDSWIZ00 AND GDSWZD00.  MAKE
!                      CROT,SORT,XLON,XLAT,YLON,YLAT AND AREA
!                      OPTIONAL ARGUMENTS.  MAKE PART OF A MODULE.
!                      MOVE VECTOR ROTATION, MAP JACOBIAN AND GRID
!                      BOX AREA COMPUTATIONS TO SEPARATE SUBROUTINES.
! 2015-07-13  GAYNO    CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAY
!                      WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAY.
!                      RENAME ROUTINE AS "GDSWZD_EQUID_CYLIND".
!
! USAGE:    CALL GDSWZD_EQUID_CYLIND(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,
!    &                  FILL,XPTS,YPTS,RLON,RLAT,NRET,
!    &                  CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
!
!   INPUT ARGUMENT LIST:
!     IGDTNUM  - INTEGER GRID DEFINITION TEMPLATE NUMBER.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                MUST BE "0" FOR EQUIDISTANT CYLINDRICAL
!                PROJECTIONS.
!     IGDTMPL  - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY.
!                CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE FOR SECTION
!                THREE:
!                 (1):  SHAPE OF EARTH, OCTET 15
!                 (2):  SCALE FACTOR OF SPHERICAL EARTH RADIUS,
!                       OCTET 16
!                 (3):  SCALED VALUE OF RADIUS OF SPHERICAL EARTH,
!                       OCTETS 17-20
!                 (4):  SCALE FACTOR OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 21
!                 (5):  SCALED VALUE OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 22-25
!                 (6):  SCALE FACTOR OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 26
!                 (7):  SCALED VALUE OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 27-30
!                 (8):  NUMBER OF POINTS ALONG A PARALLEL, OCTS 31-34
!                 (9):  NUMBER OF POINTS ALONG A MERIDIAN, OCTS 35-38
!                 (10): BASIC ANGLE OF INITIAL PRODUCTION DOMAIN,
!                       OCTETS 39-42.
!                 (11): SUBDIVISIONS OF BASIC ANGLE, OCTETS 43-46
!                 (12): LATITUDE OF FIRST GRID POINT, OCTETS 47-50
!                 (13): LONGITUDE OF FIRST GRID POINT, OCTETS 51-54
!                 (14): RESOLUTION AND COMPONENT FLAGS, OCTET 55
!                 (15): LATITUDE OF LAST GRID POINT, OCTETS 56-59
!                 (16): LONGITUDE OF LAST GRID POINT, OCTETS 60-63
!                 (17): I-DIRECTION INCREMENT, OCTETS 64-67
!                 (18): J-DIRECTION INCREMENT, OCTETS 68-71
!                 (19): SCANNING MODE, OCTET 72
!     IGDTLEN  - INTEGER NUMBER OF ELEMENTS (19) OF THE GRID DEFINITION
!                TEMPLATE ARRAY.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     IOPT     - INTEGER OPTION FLAG
!                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
!                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
!     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
!     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
!                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
!                (ACCEPTABLE RANGE: -360. TO 360.)
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
!                (ACCEPTABLE RANGE: -90. TO 90.)
!
!   OUTPUT ARGUMENT LIST:
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
!     NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
!     CROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES
!     SROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION SINES
!                (UGRID=CROT*UEARTH-SROT*VEARTH;
!                 VGRID=SROT*UEARTH+CROT*VEARTH)
!     XLON     - REAL, OPTIONAL (NPTS) DX/DLON IN 1/DEGREES
!     XLAT     - REAL, OPTIONAL (NPTS) DX/DLAT IN 1/DEGREES
!     YLON     - REAL, OPTIONAL (NPTS) DY/DLON IN 1/DEGREES
!     YLAT     - REAL, OPTIONAL (NPTS) DY/DLAT IN 1/DEGREES
!     AREA     - REAL, OPTIONAL (NPTS) AREA WEIGHTS IN M**2
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER,             INTENT(IN   ) :: IGDTNUM, IGDTLEN
 INTEGER,             INTENT(IN   ) :: IGDTMPL(IGDTLEN)
 INTEGER,             INTENT(IN   ) :: IOPT, NPTS
 INTEGER,             INTENT(  OUT) :: NRET
!
 REAL,                INTENT(IN   ) :: FILL
 REAL,                INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
 REAL,                INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
 REAL, OPTIONAL,      INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
 REAL, OPTIONAL,      INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
 REAL, OPTIONAL,      INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
!
 INTEGER                            :: IM, JM, ISCAN, N, ISCALE
!
 LOGICAL                            :: LROT, LMAP, LAREA
!
 REAL                               :: ECCEN_SQUARED
 REAL                               :: HI, RLAT1, RLON1, RLAT2, RLON2
 REAL                               :: XMAX, XMIN, YMAX, YMIN
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IF(PRESENT(CROT)) CROT=FILL
 IF(PRESENT(SROT)) SROT=FILL
 IF(PRESENT(XLON)) XLON=FILL
 IF(PRESENT(XLAT)) XLAT=FILL
 IF(PRESENT(YLON)) YLON=FILL
 IF(PRESENT(YLAT)) YLAT=FILL
 IF(PRESENT(AREA)) AREA=FILL
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! IS THIS AN EQUIDISTANT CYCLINDRICAL GRID?
 IF(IGDTNUM/=0)THEN
   CALL EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
   RETURN
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 CALL EARTH_RADIUS(IGDTMPL,IGDTLEN,RERTH,ECCEN_SQUARED)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  ENSURE EARTH RADIUS IS DEFINED.
 IF(RERTH<0.) THEN
   CALL EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
   RETURN
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IM=IGDTMPL(8)
 JM=IGDTMPL(9)
 ISCALE=IGDTMPL(10)*IGDTMPL(11)
 IF(ISCALE==0) ISCALE=10**6
 RLAT1=FLOAT(IGDTMPL(12))/FLOAT(ISCALE)
 RLON1=FLOAT(IGDTMPL(13))/FLOAT(ISCALE)
 RLAT2=FLOAT(IGDTMPL(15))/FLOAT(ISCALE)
 RLON2=FLOAT(IGDTMPL(16))/FLOAT(ISCALE)
 ISCAN=MOD(IGDTMPL(19)/128,2)
 HI=(-1.)**ISCAN
 DLON=HI*(MOD(HI*(RLON2-RLON1)-1+3600,360.)+1)/(IM-1)
 DLAT=(RLAT2-RLAT1)/(JM-1)
 XMIN=0
 XMAX=IM+1
 IF(IM.EQ.NINT(360/ABS(DLON))) XMAX=IM+2
 YMIN=0
 YMAX=JM+1
 NRET=0
 IF(PRESENT(CROT).AND.PRESENT(SROT))THEN
   LROT=.TRUE.
 ELSE
   LROT=.FALSE.
 ENDIF
 IF(PRESENT(XLON).AND.PRESENT(XLAT).AND.PRESENT(YLON).AND.PRESENT(YLAT))THEN
   LMAP=.TRUE.
 ELSE
   LMAP=.FALSE.
 ENDIF
 IF(PRESENT(AREA))THEN
   LAREA=.TRUE.
 ELSE
   LAREA=.FALSE.
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
 IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN
!$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
   DO N=1,NPTS
     IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND. &
        YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
       RLON(N)=MOD(RLON1+DLON*(XPTS(N)-1)+3600,360.)
       RLAT(N)=MIN(MAX(RLAT1+DLAT*(YPTS(N)-1),-90.),90.)
       NRET=NRET+1
       IF(LROT)  CALL EQUID_CYLIND_VECT_ROT(CROT(N),SROT(N))
       IF(LMAP)  CALL EQUID_CYLIND_MAP_JACOB(XLON(N),XLAT(N),YLON(N),YLAT(N))
       IF(LAREA) CALL EQUID_CYLIND_GRID_AREA(RLAT(N),AREA(N))
     ELSE
       RLON(N)=FILL
       RLAT(N)=FILL
     ENDIF
   ENDDO
!$OMP END PARALLEL DO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
 ELSEIF(IOPT.EQ.-1) THEN
!$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
   DO N=1,NPTS
     IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN
       XPTS(N)=1+HI*MOD(HI*(RLON(N)-RLON1)+3600,360.)/DLON
       YPTS(N)=1+(RLAT(N)-RLAT1)/DLAT
       IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND.  &
          YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
         NRET=NRET+1
         IF(LROT)  CALL EQUID_CYLIND_VECT_ROT(CROT(N),SROT(N))
         IF(LMAP)  CALL EQUID_CYLIND_MAP_JACOB(XLON(N),XLAT(N),YLON(N),YLAT(N))
         IF(LAREA) CALL EQUID_CYLIND_GRID_AREA(RLAT(N),AREA(N))
       ELSE
         XPTS(N)=FILL
         YPTS(N)=FILL
       ENDIF
     ELSE
       XPTS(N)=FILL
       YPTS(N)=FILL
     ENDIF
   ENDDO
!$OMP END PARALLEL DO
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE GDSWZD_EQUID_CYLIND
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 SUBROUTINE EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  EQUID_CYLIND_ERROR   ERROR HANDLER
!   PRGMMR: GAYNO       ORG: W/NMC23       DATE: 2015-07-13
!
! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
!           A "FILL" VALUE TO THE OUTPUT FIELDS.
!
! PROGRAM HISTORY LOG:
! 2015-07-13  GAYNO     INITIAL VERSION
!
! USAGE:    CALL EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
!
!   INPUT ARGUMENT LIST:
!     IOPT     - INTEGER OPTION FLAG
!                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
!                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
!     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
!     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
!                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
!   OUTPUT ARGUMENT LIST:
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER, INTENT(IN   ) :: IOPT, NPTS
!
 REAL,    INTENT(IN   ) :: FILL
 REAL,    INTENT(  OUT) :: RLAT(NPTS),RLON(NPTS)
 REAL,    INTENT(  OUT) :: XPTS(NPTS),YPTS(NPTS)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IF(IOPT>=0) THEN
   RLON=FILL
   RLAT=FILL
 ENDIF
 IF(IOPT<=0) THEN
   XPTS=FILL
   YPTS=FILL
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE EQUID_CYLIND_ERROR
!
 SUBROUTINE EQUID_CYLIND_VECT_ROT(CROT,SROT)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  EQUID_CYLIND_VECT_ROT   VECTOR ROTATION FIELDS FOR
!                                      EQUIDISTANT CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE VECTOR ROTATION SINES AND
!           COSINES FOR AN EQUIDISTANT CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "EQUID_CYLIND_VECT_ROT"
!
! USAGE:    CALL EQUID_CYLIND_VECT_ROT(CROT,SROT)
!
!   INPUT ARGUMENT LIST:
!     NONE
!
!   OUTPUT ARGUMENT LIST:
!     CROT     - CLOCKWISE VECTOR ROTATION COSINES (REAL)
!     SROT     - CLOCKWISE VECTOR ROTATION SINES (REAL)
!                (UGRID=CROT*UEARTH-SROT*VEARTH;
!                 VGRID=SROT*UEARTH+CROT*VEARTH)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE

 REAL,                INTENT(  OUT) :: CROT, SROT

 CROT=1.0
 SROT=0.0

 END SUBROUTINE EQUID_CYLIND_VECT_ROT
!
 SUBROUTINE EQUID_CYLIND_MAP_JACOB(XLON,XLAT,YLON,YLAT)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  EQUID_CYLIND_MAP_JACOB  MAP JACOBIANS FOR
!                                      EQUIDISTANT CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE MAP JACOBIANS FOR
!           AN EQUIDISTANT CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "EQUID_CYLIND_MAP_JACOB"
!
! USAGE:  CALL EQUID_CYLIND_MAP_JACOB(XLON,XLAT,YLON,YLAT)
!
!   INPUT ARGUMENT LIST:
!     NONE
!
!   OUTPUT ARGUMENT LIST:
!     XLON     - DX/DLON IN 1/DEGREES (REAL)
!     XLAT     - DX/DLAT IN 1/DEGREES (REAL)
!     YLON     - DY/DLON IN 1/DEGREES (REAL)
!     YLAT     - DY/DLAT IN 1/DEGREES (REAL)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$

 REAL,                INTENT(  OUT) :: XLON,XLAT,YLON,YLAT

 XLON=1.0/DLON
 XLAT=0.
 YLON=0.
 YLAT=1.0/DLAT

 END SUBROUTINE EQUID_CYLIND_MAP_JACOB
!
 SUBROUTINE EQUID_CYLIND_GRID_AREA(RLAT,AREA)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  EQUID_CYLIND_GRID_AREA  GRID BOX AREA FOR
!                                      EQUIDISTANT CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE GRID BOX AREA FOR
!           AN EQUIDISTANT CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "EQUID_CYLIND_GRID_AREA"
!
! USAGE:  CALL EQUID_CYLIND_GRID_AREA(RLAT,AREA)
!
!   INPUT ARGUMENT LIST:
!     RLAT     - LATITUDE OF GRID POINT IN DEGREES (REAL)
!
!   OUTPUT ARGUMENT LIST:
!     AREA     - AREA WEIGHTS IN M**2 (REAL)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE

 REAL,                INTENT(IN   ) :: RLAT
 REAL,                INTENT(  OUT) :: AREA

 REAL,                PARAMETER     :: PI=3.14159265358979
 REAL,                PARAMETER     :: DPR=180./PI

 REAL                               :: DSLAT, RLATU, RLATD

 RLATU=MIN(MAX(RLAT+DLAT/2,-90.),90.)
 RLATD=MIN(MAX(RLAT-DLAT/2,-90.),90.)
 DSLAT=SIN(RLATU/DPR)-SIN(RLATD/DPR)
 AREA=RERTH**2*ABS(DSLAT*DLON)/DPR

 END SUBROUTINE EQUID_CYLIND_GRID_AREA

 END MODULE GDSWZD_EQUID_CYLIND_MOD

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM


module fcm

use mpi
use settings
use partitioner, only: mesh,m
implicit none

! Delta function = (2*pi*sigma**2)^(-3/2) * exp( (x-Yn)**2/(2*sigma**2) )
! to be calculated once at each timestep for all particles
double precision, dimension(:,:,:,:), allocatable :: delta_fn

contains

subroutine calc_delta_fn(Yp)
    implicit none
    double precision, dimension(Np,3) :: Yp
    double precision :: temp2
    double precision, dimension(nz,ny,m) :: dist_sqr, temp1
    double precision :: sigma
    integer :: i

    if (.not. allocated(delta_fn)) allocate(delta_fn(nz,ny,1:m,Np))

    sigma = a/sqrt(pi)

    temp2 = (2*pi*sigma**2)**(-1.5)

    do i=1,Np
        dist_sqr = (mesh%x-Yp(i,1))**2 + (mesh%y-Yp(i,2))**2 + (mesh%z-Yp(i,3))**2
	    temp1 = -dist_sqr/(2*sigma**2)
        delta_fn(:,:,:,i) = temp2*exp(temp1)
    end do
end subroutine calc_delta_fn

subroutine updateParticles(Yp,Vp,Vp0,u,v,w)

	implicit none

	double precision, dimension(nz,ny,0:m+1) :: u,v,w
	double precision, dimension(Np,3) :: Yp, Vp, Vp0

    double precision :: dv
    double precision, dimension(Np,3) :: V_loc

	integer :: i

    integer :: id, Nproc, ierr

    !call mpi_comm_rank(MPI_COMM_WORLD, id, ierr)
    !call mpi_comm_size(MPI_COMM_WORLD, Nproc, ierr)

    if (.not. allocated(delta_fn)) then
        !print*, 'Error (forceCoupling): call subroutine calc_delta_fn() before calling forceCoupling()'
        call calc_delta_fn(Yp)
    end if

	Vp0 = Vp

    dv = dx*dy*dz

    do i=1,Np
	!do i=id*(Np/Nproc)+1,min((id+1)*(Np/Nproc),Np)

		!dist_sqr = (x-Yp(i,1))**2 + (y-Yp(i,2))**2 + (z-Yp(i,3))**2

		!temp1 = -dist_sqr/(2*sigma**2)

		V_loc(i,1) = sum(u(:,:,1:m)*delta_fn(:,:,:,i))*dv
		V_loc(i,2) = sum(v(:,:,1:m)*delta_fn(:,:,:,i))*dv
		V_loc(i,3) = sum(w(:,:,1:m)*delta_fn(:,:,:,i))*dv

	end do

    !call mpi_bcast(Vp(id*(Np/Nproc)+1:min((id+1)*(Np/Nproc),Np),:),3*(min((id+1)*(Np/Nproc)-id*(Np/Nproc)+1,Np)),MPI_DOUBLE_PRECISION,id,MPI_COMM_WORLD,ierr)

    call mpi_allreduce(V_loc(:,1),Vp(:,1),Np,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
    call mpi_allreduce(V_loc(:,2),Vp(:,2),Np,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
    call mpi_allreduce(V_loc(:,3),Vp(:,3),Np,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)

    !print*, V_loc(1,1), V_loc(1,2), V_loc(1,3)

	Yp = Yp + dt*Vp

end subroutine updateParticles

subroutine forceCoupling(Gx,Gy,Gz,Yp)

	double precision, dimension(nz,ny,1:m) :: Gx,Gy,Gz
	double precision, dimension(Np,3) :: Yp

	double precision, dimension(Np,3) :: F, F_stokes, F_part_coll, F_wall_coll, F_vwf, F_pivkin

	integer :: i

	!integer :: error, Nproc, id

    !print*, 'entered forceCoupling sub'

    call calc_delta_fn(Yp)

	!call mpi_comm_size(MPI_COMM_WORLD, Nproc, error)
    !call mpi_comm_rank(MPI_COMM_WORLD, id, error)

    !if (id==0) then
    !    print*, 'nprocs=', Nproc
    !end if
    !print*, id

	! sigma = a/sqrt(pi)

	! assume particles are neutrally buoyant. therefore no inertia term

	!F =
	!F = 1e-15;

	!temp2 = (2*pi*sigma**2)**(-1.5)

	! Stokes force to ensure terminal velocity of U0
	!F_stokes(:,1) = 6*pi*mu*a*U0
	F_stokes(:,1) = 0
	F_stokes(:,2:3) = 0

	!call particleCollisionForce(F_part_coll,Yp)
	!call wallCollisionForce(F_wall_coll,Yp)
	!call vwfForce(F_vwf,Yp)

	! Exclude other forces for now
	F_part_coll = 0
	F_wall_coll = 0
	F_vwf = 0

	!call pivkinForceField(F_pivkin,Yp)
	F_pivkin = 0

	F = F_stokes + F_part_coll + F_wall_coll + F_vwf + F_pivkin
	!print*, F, -6*pi*mu*a*U0
	!F(:,1) = -6*pi*mu*a*U0; F(:,2:3) = 0
	!F = 0

	Gx = 0; Gy  =0; Gz = 0;

    !print*, id*(Np/Nproc)+1,min((id+1)*(Np/Nproc),Np)

    ! local sum
	!do i=id*(Np/Nproc)+1,min((id+1)*(Np/Nproc),Np)
	do i=1,Np

		!dist_sqr = (x-Yp(i,1))**2 + (y-Yp(i,2))**2 + (z-Yp(i,3))**2

		!temp1 = -dist_sqr/(2*sigma**2)

		Gx = Gx + F(i,1)*delta_fn(:,:,:,i)/rho
		Gy = Gy + F(i,2)*delta_fn(:,:,:,i)/rho
		Gz = Gz + F(i,3)*delta_fn(:,:,:,i)/rho

	end do

    !print *, 'before allreduce'

    !call mpi_allreduce(Gx_loc,Gx,nx*ny*nz,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,error)
    !call mpi_allreduce(Gy_loc,Gy,nx*ny*nz,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,error)
    !call mpi_allreduce(Gz_loc,Gz,nx*ny*nz,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,error)

    !print*, 'after allreduce'

end subroutine forceCoupling

subroutine pivkinForceField(F,Yp)

	double precision, dimension(Np,3) :: Yp
	double precision, dimension(Np,3) :: F

	double precision, dimension(3) :: r, r_wall

	double precision :: r_mag, r_wall_mag

	integer :: i,j

	! parameters
	double precision, parameter :: R_piv = 1.0
	double precision, parameter :: L_piv = 1.5
	double precision, parameter :: B_piv = 3.5
	double precision, parameter :: M_piv = (L_piv+B_piv)/2
	double precision, parameter :: alpha_piv = 4e-12

	do i=1,Np
		F(i,:) = 0

		! Wall
		r_wall_mag = a - sqrt(Yp(i,2)**2+Yp(i,3)**2)
		r_wall(1)=0; r_wall(2)=Yp(i,2);	r_wall(3)=Yp(i,3);
		r_wall = -r_wall_mag*r_wall/sqrt(sum(r_wall**2))	! directed away from wall
		if (r_wall_mag/a <= R_piv) then
			F(i,:) = F(i,:) + 2*alpha_piv*(1-(r_wall_mag/a)/R_piv)*r_wall/r_wall_mag
		end if

		! other particles
		do j=1,Np
			if (i/=j) then
				r = Yp(i,:)-Yp(j,:)	! directed j -> i
				r_mag = sqrt(sum(r**2))
				if (r_mag/a <= R_piv) then
					F(i,:) = F(i,:) + 2*alpha_piv*(1-(r_mag/a)/R_piv)*r/r_mag ! repulsive
				else if ( (r_mag/a >= R_piv) .and. (r_mag/a <= L_piv) ) then
					F(i,:) = F(i,:) + 0
				else if ( (r_mag/a >= L_piv) .and. (r_mag/a <= M_piv) ) then
					F(i,:) = F(i,:) - alpha_piv*((r_mag/a)/L_piv-1)*r/r_mag ! attractive
				else if ( (r_mag/a >= M_piv) .and. (r_mag/a <= B_piv) ) then
					F(i,:) = F(i,:) - alpha_piv*(M_piv/L_piv-1)*((B_piv-(r_mag/a))/(B_piv-M_piv))*r/r_mag ! attractive
				end if
			end if
		end do

	end do

end subroutine pivkinForceField

subroutine particleCollisionForce(F,Yp)

	double precision, dimension(Np,3) :: Yp

	double precision :: R_ref
	double precision :: F_ref
	double precision, dimension(3) :: r			! r vector from j to i
	double precision, dimension(Np,3) :: F		! force vector for each particle
	double precision :: r_mag

	integer :: i, j

	! Parameters
	R_ref = 2.001*a
	F_ref = 1e-20

	do i=1,Np
		F(i,:) = 0

		! other particles
		do j=1,Np
			if (i/=j) then
				r = Yp(i,:)-Yp(j,:)
				r_mag = sqrt(sum(r**2))
				if (r_mag < R_ref) then
					F(i,:) = F(i,:) + F_ref*( (R_ref**2 - r_mag**2)/(R_ref**2 - 4*a**2) )**2 * r / r_mag
				end if
			end if
		end do

	end do

end subroutine particleCollisionForce

subroutine wallCollisionForce(F,Yp)

	double precision, dimension(Np,3) :: F		! force vector for each particle
	double precision, dimension(Np,3) :: Yp

	double precision, dimension(3) :: r					! axis -> Yp
	double precision :: r_mag

	double precision :: R_ref		! close to channel RADIUS
	double precision :: F_ref		! force

	integer :: i

	! Parameters
	R_ref = 0.999*RADIUS
	F_ref = 1e-10

	F = 0

	do i=1,Np
		r(1) = 0; r(2) = Yp(i,2)-Ly/2; r(3) = Yp(i,3)-Lz/2
		r_mag = sqrt(sum(r**2))
		if (r_mag > R_ref) then
			F(i,:) = F_ref*( (r_mag**2 - R_ref**2)/(RADIUS**2 - R_ref**2) )**2 * r / r_mag
		end if
	end do

end subroutine wallCollisionForce

subroutine vwfForce(F,Yp)
	double precision, dimension(Np,3) :: F		! force vector for each particle
	double precision, dimension(Np,3) :: Yp

	double precision :: x0, y0, z0 	! location of injury

	double precision, dimension(3) :: r			! r vector from j to i
	double precision :: r_mag
	double precision :: vwfRange

	double precision :: F_ref

	integer :: i

	x0 = 30e-6; y0=0.0; z0 = Lz/2;
	vwfRange = 5e-6
	F_ref = 1e-10

	F = 0

	do i=1,Np

		r(1) = Yp(i,1) - x0; r(2) = Yp(i,2) - y0; r(3) = Yp(i,3) - z0;
		r_mag = sqrt(sum(r**2))

		if (r_mag<vwfRange) then
			F(i,:) = -F_ref* r
		end if
	end do

end subroutine vwfForce

subroutine fbgForce(F,Yp)

	! depends on local fbg concentration
	! To be added when biochem is implemented

	double precision, dimension(Np,3) :: F		! force vector for each particle
	double precision, dimension(Np,3) :: Yp

	double precision, dimension(3) :: r			! r vector from j to i
	double precision :: r_mag

	integer :: i, j

	do i=1,Np
		F(i,:) = 0

	end do

end subroutine fbgForce

end module fcm


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE DLWQPP ( NOTOT  , NOSYS  , NOSEG  , NOPA   , NOSFUN ,
     *                    ITIME  , IMFLAG , IDFLAG , IHFLAG , MONAME ,
     *                    SYNAME , DUNAME , WANAME , IDUMP  , NODUMP ,
     *                    IWASTE , NOWAST , CONC   , CONS   , PARAM  ,
     *                    FUNC   , SEGFUN , VOLUME , WASTE  , BOUND  ,
     *                    NOBND  , ITSTRT , ITSTOP , NX     , NY     ,
     *                    LGRID  , NODISP , NOVELO , NOQ    , NOQ1   ,
     *                    NOQ2   , NOQ3   , DISPER , VELO   , ASMASS ,
     *                    IBFLAG , NOCONS , NOFUN  , CONAME , PANAME ,
     *                    FUNAME , SFNAME , BONAME )
!
!
!     Deltares      SECTOR WATERRESOURCES AND ENVIRONMENT
!
!     CREATED             : june 1988 by L. Postma
!
!     FUNCTION            : Parameter list and header for user supplied
!                           subroutine for POSTPROCESSING.
!
!     LOGICAL UNITS       : none explicitly, the user may use any unit
!                           number if NOT between 10 and 35 !!!!!!!!!!
!
!     SUBROUTINES CALLED  : none explicitly
!
!     PARAMETERS          :
!
!     NAME    KIND     LENGTH     FUNCT.  DESCRIPTION
!     ----    -----    ------     ------- -----------
!     NOTOT   INTEGER       1     INPUT   Total number of substances
!     NOSYS   INTEGER       1     INPUT   Number of active substances
!     NOSEG   INTEGER       1     INPUT   Nr. of computational elements
!     NOPA    INTEGER       1     INPUT   Number of parameters
!     NOSFUN  INTEGER       1     INPUT   Number of segment functions
!     ITIME   INTEGER       1     INPUT   Time in system clock units
!     IMFLAG  LOGICAL       1     INPUT   .TRUE. if DELWAQ sets a
!                                                     monitoring step
!     IDFLAG  LOGICAL       1     INPUT   as IMFLAG for dump actions
!     IHFLAG  LOGICAL       1     INPUT   as IMFLAG for history actions
!     MONAME  CHAR*40       4     INPUT   Model and run names
!     SYNAME  CHAR*20    NOTOT    INPUT   names of systems
!     DUNAME  CHAR*20    NODUMP   INPUT   names of dump locations
!     WANAME  CHAR*20    NOWAST   INPUT   names of waste locations
!     IDUMP   INTEGER    NODUMP   INPUT   dump segment numbers
!     NODUMP  INTEGER       1     INPUT   number of dump locations
!     IWASTE  INTEGER    NOWAST   INPUT   waste segment numbers
!     NOWAST  INTEGER       1     INPUT   number of waste locations
!     CONC    REAL   NOTOT,NOSEG  INPUT   Model concentrations
!     CONS    REAL          *     IN/OUT  Model constants
!     PARAM   REAL    NOPA,NOSEG  IN/OUT  Model parameters
!     FUNC    REAL          *     IN/OUT  Model functions at ITIME
!     SEGFUN  REAL   NOSEG,NOSFUN IN/OUT  Segment functions at ITIME
!     VOLUME  REAL      NOSEG     INPUT   Segment volumes
!     WASTE   REAL NOTOT+1,NOWAST INPUT   waste loads
!     BOUND   REAL   NOSYS,NOBND  INPUT   boundary concentrations
!     NOBND   INTEGER     1       INPUT   number of boundary conditions
!     ITSTRT  INTEGER     1       INPUT   model start time in units
!     ITSTOP  INTEGER     1       INPUT   model stop  time in units
!     NX      INTEGER     1       INPUT   width of grid
!     NY      INTEGER     1       INPUT   depth of grid
!     LGRID   INTEGER     NX*NY   INPUT   grid-layout
!     NOQ     INTEGER       1     INPUT   Total number of exchanges
!     NOQ1    INTEGER       1     INPUT   number of exchanges 1st direction
!     NOQ2    INTEGER       1     INPUT   number of exchanges 2nd direction
!     NOQ3    INTEGER       1     INPUT   number of exchanges 3rd direction
!     NOVELO  INTEGER       1     INPUT   Number of user-flows
!     NODISP  INTEGER       1     INPUT   Number of user-dispersions
!     DISPER  REAL   NODISP*NOQ   OUTPUT  User defined dispersion
!     VELO    REAL   NOVELO*NOQ   OUTPUT  User defined flows
!     ASMASS  REAL  NOTOT*NOSEG*? IN/OUT  mass balance per comp. elems.
!                                         if IBFLAG = 1
!     IBFLAG  INTEGER    1        INPUT   if 1 then mass balance p.c.e.
!     NOCONS  INTEGER       1     INPUT   Number of constants used
!     NOFUN   INTEGER       1     INPUT   Number of functions ( user )
!     CONAME  CHAR*20   NOCONS    INPUT   Constant names
!     PANAME  CHAR*20   NOPA      INPUT   Parameter names
!     FUNAME  CHAR*20   NOFUN     INPUT   Function names
!     SFNAME  CHAR*20   NOSFUN    INPUT   Segment function names
!     BONAME  CHAR*20   NOBND     INPUT   Boundary names
!
!     ==================================================================
!
!
!
      CHARACTER*40 MONAME(4)
      CHARACTER*20 SYNAME (NOTOT), DUNAME (NODUMP), WANAME (NOWAST),
     &             CONAME (*)    , PANAME (*)     , FUNAME (*)     ,
     &             SFNAME (*)    , BONAME (*)
      DIMENSION    IDUMP (NODUMP), IWASTE (NOWAST)
      DIMENSION    VOLUME (NOSEG), CONS        (*),FUNC        (*),
     &             CONC   (NOTOT,NOSEG)           ,PARAM   (NOPA,NOSEG),
     &             WASTE (0:NOTOT,NOWAST)         ,SEGFUN(NOSEG,NOSFUN),
     &             BOUND ( NOSYS, NOBND )           ,LGRID (NX    ,NY ),
     &             DISPER( NODISP, NOQ  )           ,VELO  (NOVELO,NOQ),
     &             ASMASS( NOTOT,NOSEG,*)
      LOGICAL IMFLAG, IDFLAG, IHFLAG
!
!
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd