text:
!*******************************************************************************
!>  @file bmw_parallel_context.f
!>  @brief Contains module @ref bmw_parallel_context.
!
!  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 bmw_parallel_context_class. This
!>  contains the state variables needed by BMW for parallel computation.
!*******************************************************************************
      MODULE bmw_parallel_context
      USE mpi_inc
      USE profiler

      IMPLICIT NONE

!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) bmw parallel context base class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Base class representing a bmw parallel context. This contains all memory
!>  needed parameters needed to parallel processing.
!-------------------------------------------------------------------------------
      TYPE :: bmw_parallel_context_class
!>  Problem space offset.
         INTEGER :: offset
!>  Problem space stride.
         INTEGER :: stride

!>  Number of threads.
         INTEGER :: num_threads

#if defined (MPI_OPT)
!>  MPI communicator reference for computation. BMW can be used as an embeded
!>  library. Multiple instances can be running in parallel on separate problems
!>  so do not assume the world comm.
         INTEGER :: comm = MPI_COMM_NULL
!>  Flag to mark if bmw intialized the MPI context.
         LOGICAL :: initialized_mpi
#endif
      CONTAINS
         PROCEDURE, PASS ::                                                    &
     &      set_threads => bmw_parallel_context_set_threads
         PROCEDURE, PASS :: report => bmw_parallel_context_report
         PROCEDURE, PASS :: reduce1 => bmw_parallel_context_reduce1
         PROCEDURE, PASS :: reduce2 => bmw_parallel_context_reduce2
         PROCEDURE, PASS :: reduce3 => bmw_parallel_context_reduce3
         PROCEDURE, PASS :: reduce4 => bmw_parallel_context_reduce4
         GENERIC         :: reduce => reduce1, reduce2, reduce3, reduce4
         PROCEDURE, PASS :: start => bmw_parallel_context_start
         PROCEDURE, PASS :: end => bmw_parallel_context_end
         FINAL           :: bmw_parallel_context_destruct
      END TYPE

!*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Interface for the bmw_parallel_context constructor.
!-------------------------------------------------------------------------------
      INTERFACE bmw_parallel_context_class
         MODULE PROCEDURE bmw_parallel_context_construct
      END INTERFACE

      CONTAINS

!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Construct a @ref bmw_context_class object.
!>
!>  Allocates memory and initializes a @ref bmw_parallel_context_class object.
!>
!>  @param[in] comm        MPI communicator. Only used if MPI support has been
!>                         compiled in.
!>  @returns A pointer to a constructed @ref bmw_parallel_context_class object.
!-------------------------------------------------------------------------------
      FUNCTION bmw_parallel_context_construct(                                 &
#if defined (MPI_OPT)
     &                                        comm                             &
#endif
     &                                       )

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), POINTER ::                           &
     &   bmw_parallel_context_construct
#if defined (MPI_OPT)
      INTEGER, INTENT(in)                        :: comm
#endif

!  local variables
#if defined (MPI_OPT)
      LOGICAL                                    :: isinit
      INTEGER                                    :: status
#endif
      REAL (rprec)                               :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      ALLOCATE(bmw_parallel_context_construct)

#if defined (MPI_OPT)
      CALL MPI_INITIALIZED(isinit, status)

      IF (.not.isinit) THEN
         bmw_parallel_context_construct%initialized_mpi = .true.
         CALL MPI_INIT(status)
      END IF

      bmw_parallel_context_construct%comm = comm
      CALL MPI_COMM_RANK(bmw_parallel_context_construct%comm,                  &
     &                   bmw_parallel_context_construct%offset,                &
     &                   status)
      CALL MPI_COMM_SIZE(bmw_parallel_context_construct%comm,                  &
     &                   bmw_parallel_context_construct%stride,                &
     &                   status)

#else
      bmw_parallel_context_construct%offset = 0
      bmw_parallel_context_construct%stride = 1
#endif

!  Configure the number of threads to use.
      bmw_parallel_context_construct%num_threads = 1

      CALL profiler_set_stop_time('bmw_parallel_context_construct',            &
     &                            start_time)

      END FUNCTION

!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Deconstruct a @ref bmw_parallel_context_class object.
!>
!>  Deallocates memory and uninitializes a @ref bmw_parallel_context_class
!>  object.
!>
!>  @param[inout] this     A @ref bmw_parallel_context_class instance.
!>  @param[in]    finalize Flag to call MPI_FINALIZE. Only used if MPI support
!>                         has been compiled in.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_destruct(this)

      IMPLICIT NONE

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

!  local variables
#if defined (MPI_OPT)
      INTEGER                                          :: status

!  Start of executable code
      IF (this%initialized_mpi) THEN
         CALL MPI_FINALIZE(status)
      END IF
#endif

      END SUBROUTINE

!*******************************************************************************
!  SETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Set the number of threads.
!>
!>  Sets the number of OpenMP threads to use.
!>
!>  @param[in] this        A @ref bmw_parallel_context_class instance.
!>  @param[in] num_threads Number of threads to use.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_set_threads(this, num_threads)
!$    USE omp_lib

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(inout) :: this
      INTEGER, INTENT(in)                               :: num_threads

!  local variables
      REAL (rprec)                                      :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

!  Configure the number of threads to use.
      this%num_threads = num_threads

!$    IF (this%num_threads .gt. 0) THEN
!$       CALL OMP_SET_NUM_THREADS(this%num_threads)
!$    END IF
!$OMP PARALLEL
!$    this%num_threads = OMP_GET_MAX_THREADS()
!$OMP END PARALLEL

      CALL profiler_set_stop_time('bmw_parallel_context_set_threads',          &
     &                            start_time)

      END SUBROUTINE

!*******************************************************************************
!  UTILITY SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Abort the entire program.
!>
!>  Aborts everything.
!>
!>  @param[in] status Error number to abort with.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_abort(status)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER, INTENT(in) :: status

!  local variables
      INTEGER             :: error

!  Start of executable code
#if defined (MPI_OPT)
      CALL MPI_ABORT(MPI_COMM_WORLD, status, error)
#else
      CALL EXIT(status)
#endif

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Report the number of parallel processes and threads.
!>
!>  Provides a summary report if the number of threads and processes currently
!>  configured.
!>
!>  @param[in] this    A @ref bmw_parallel_context_class instance.
!>  @param[in] io_unit Unit number to write messages to.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_report(this, io_unit)

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(inout) :: this
      INTEGER, INTENT(in)                               :: io_unit

!  local variables
      REAL (rprec)                                      :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      IF (this%offset .eq. 0) THEN
#if defined(MPI_OPT)
         WRITE (io_unit,1000) this%stride
#endif
         WRITE (io_unit,1001) this%num_threads
      END IF

      CALL profiler_set_stop_time('bmw_parallel_context_report',               &
     &                            start_time)

1000  FORMAT('Using ',i4,' processes.')
#if defined (MPI_OPT)
1001  FORMAT('Using ',i4,' threads per process.')
#else
1001  FORMAT('Using ',i4,' threads.')
#endif

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Reduce parallel buffers.
!>
!>  Reduce a buffer from all parallel processes. In the single processes case,
!>  this reduces to nothing. This reduces a 1D buffer.
!>
!>  @param[in]    this   A @ref bmw_parallel_context_class instance.
!>  @param[inout] buffer Buffer to reduce.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_reduce1(this, buffer)
      USE stel_kinds, ONLY: rprec

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(in) :: this
      REAL (rprec), DIMENSION(:), INTENT(inout)      :: buffer

!  local variables
      REAL (rprec)                                   :: start_time
      INTEGER                                        :: status

!  Start of executable code
      start_time = profiler_get_start_time()

#if defined (MPI_OPT)
      CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, SIZE(buffer), MPI_REAL8,        &
     &                   MPI_SUM, this%comm, status)
#endif

      CALL profiler_set_stop_time('bmw_parallel_context_reduce1',              &
     &                            start_time)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Reduce parallel buffers.
!>
!>  Reduce a buffer from all parallel processes. In the single processes case,
!>  this reduces to nothing. This reduces a 2D buffer.
!>
!>  @param[in]    this    A @ref bmw_parallel_context_class instance.
!>  @param[inout] buffer Buffer to reduce.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_reduce2(this, buffer)
      USE stel_kinds, ONLY: rprec

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(in) :: this
      REAL (rprec), DIMENSION(:,:), INTENT(inout)    :: buffer

!  local variables
      REAL (rprec)                                   :: start_time
      INTEGER                                        :: status

!  Start of executable code
      start_time = profiler_get_start_time()

#if defined (MPI_OPT)
      CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, SIZE(buffer), MPI_REAL8,        &
     &                   MPI_SUM, this%comm, status)
#endif

      CALL profiler_set_stop_time('bmw_parallel_context_reduce2',              &
     &                            start_time)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Reduce parallel buffers.
!>
!>  Reduce a buffer from all parallel processes. In the single processes case,
!>  this reduces to nothing. This reduces a 3D buffer.
!>
!>  @param[in]    this   A @ref bmw_parallel_context_class instance.
!>  @param[inout] buffer Buffer to reduce.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_reduce3(this, buffer)
      USE stel_kinds, ONLY: rprec

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(in) :: this
      REAL (rprec), DIMENSION(:,:,:), INTENT(inout)  :: buffer

!  local variables
      REAL (rprec)                                   :: start_time
      INTEGER                                        :: status

!  Start of executable code
      start_time = profiler_get_start_time()

#if defined (MPI_OPT)
      CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, SIZE(buffer), MPI_REAL8,        &
     &                   MPI_SUM, this%comm, status)
#endif

      CALL profiler_set_stop_time('bmw_parallel_context_reduce3',              &
     &                            start_time)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Reduce parallel buffers.
!>
!>  Reduce a buffer from all parallel processes. In the single processes case,
!>  this reduces to nothing. This reduces a 4D buffer.
!>
!>  @param[in]    this   A @ref bmw_parallel_context_class instance.
!>  @param[inout] buffer Buffer to reduce.
!-------------------------------------------------------------------------------
      SUBROUTINE bmw_parallel_context_reduce4(this, buffer)
      USE stel_kinds, ONLY: rprec

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bmw_parallel_context_class), INTENT(in)  :: this
      REAL (rprec), DIMENSION(:,:,:,:), INTENT(inout) :: buffer

!  local variables
      REAL (rprec)                                    :: start_time
      INTEGER                                         :: status

!  Start of executable code
      start_time = profiler_get_start_time()

#if defined (MPI_OPT)
      CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, SIZE(buffer), MPI_REAL8,        &
     &                   MPI_SUM, this%comm, status)
#endif

      CALL profiler_set_stop_time('bmw_parallel_context_reduce4',              &
     &                            start_time)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Compute the start index of a flat array.
!>
!>  The total problem size is divided by the number of processes. Any remaining
!>  elements are added to the first few processes.
!>
!>  @param[in] this  A @ref bmw_parallel_context_class instance.
!>  @param[in] total Total number of elements in the loop.
!>  @returns The Starting index for parallel multiprocess computation.
!-------------------------------------------------------------------------------
      PURE FUNCTION bmw_parallel_context_start(this, total)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER :: bmw_parallel_context_start
      CLASS (bmw_parallel_context_class), INTENT(in) :: this
      INTEGER, INTENT(in)                            :: total

!  local variables
      INTEGER                                        :: n_per_process
      INTEGER                                        :: n_left

!  Start of executable code
!  Need to divide the problem space evenly by the number of processes.
      n_per_process = total/this%stride

!  The number of items may not evenly divide by the number of processes.
      n_left = MOD(total, this%stride)

      bmw_parallel_context_start = 1 + this%offset*n_per_process

!  Acound for the remaining elements in the first n_left processes.
      IF (this%offset .lt. n_left) THEN
         bmw_parallel_context_start = bmw_parallel_context_start               &
     &                              + this%offset
      ELSE
         bmw_parallel_context_start = bmw_parallel_context_start               &
     &                              + n_left
      END IF

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Compute the end index of a flat array.
!>
!>  The total problem size is divided by the number of processes. Any remaining
!>  elements are added to the first few processes.
!>
!>  @param[in] this  A @ref bmw_parallel_context_class instance.
!>  @param[in] total Total number of elements in the loop.
!>  @returns The Ending index for parallel multiprocess computation.
!-------------------------------------------------------------------------------
      PURE FUNCTION bmw_parallel_context_end(this, total)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER :: bmw_parallel_context_end
      CLASS (bmw_parallel_context_class), INTENT(in) :: this
      INTEGER, INTENT(in)                            :: total

!  local variables
      INTEGER                                        :: n_per_process
      INTEGER                                        :: n_left

!  Start of executable code
!  Need to divide the problem space evenly by the number of processes.
      n_per_process = total/this%stride

!  The number of items may not evenly divide by the number of processes.
      n_left = MOD(total, this%stride)

      bmw_parallel_context_end = (this%offset + 1)*n_per_process

!  Acound for the remaining elements in the first n_left processes.
      IF (this%offset .lt. n_left) THEN
         bmw_parallel_context_end = bmw_parallel_context_end               &
     &                            + this%offset + 1
      ELSE
         bmw_parallel_context_end = bmw_parallel_context_end               &
     &                            + n_left
      END IF

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Compute the i index of a flat array.
!>
!>  The total problem size represents a three dimensional space. From the flat
!>  index, compute the i index.
!>
!>  @param[in] index Flat index.
!>  @param[in] num_i Size of the ith dimension.
!>  @returns The index of the i dimension.
!-------------------------------------------------------------------------------
      PURE FUNCTION bmw_parallel_context_i(index, num_i)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER :: bmw_parallel_context_i
      INTEGER, INTENT(in) :: index
      INTEGER, INTENT(in) :: num_i

!  Start of executable code
      bmw_parallel_context_i = MOD(index - 1, num_i) + 1

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Compute the j index of a flat array.
!>
!>  The total problem size represents a three dimensional space. From the flat
!>  index, compute the j index.
!>
!>  @param[in] index Flat index.
!>  @param[in] num_i Size of the ith dimension.
!>  @param[in] num_j Size of the jth dimension.
!>  @returns The index of the j dimension.
!-------------------------------------------------------------------------------
      PURE FUNCTION bmw_parallel_context_j(index, num_i, num_j)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER :: bmw_parallel_context_j
      INTEGER, INTENT(in) :: index
      INTEGER, INTENT(in) :: num_i
      INTEGER, INTENT(in) :: num_j

!  Start of executable code
      bmw_parallel_context_j = MOD((index - 1)/num_i, num_j) + 1

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Compute the k index of a flat array.
!>
!>  The total problem size represents a three dimensional space. From the flat
!>  index, compute the k index.
!>
!>  @param[in] index Flat index.
!>  @param[in] num_i Size of the ith dimension.
!>  @param[in] num_j Size of the jth dimension.
!>  @returns The index of the k dimension.
!-------------------------------------------------------------------------------
      PURE FUNCTION bmw_parallel_context_k(index, num_i, num_j)

      IMPLICIT NONE

!  Declare Arguments
      INTEGER :: bmw_parallel_context_k
      INTEGER, INTENT(in) :: index
      INTEGER, INTENT(in) :: num_i
      INTEGER, INTENT(in) :: num_j

!  Start of executable code
      bmw_parallel_context_k = (index - 1)/(num_j*num_i) + 1

      END FUNCTION

      END MODULE

max_stars_repo_name: ORNL-Fusion/BMW

max_issues_repo_name: ORNL-Fusion/BMW

max_forks_repo_name: ORNL-Fusion/BMW


text:
Fantastical mustachioedmanager of the Bike Barn. Knows everything about bikes and their illustrious history. Stop by the Bike Barn to meet the legend, or take a class through the Experimental College on bike maintenance or wheel building.

His interests include solitaire, spidersolitaire, and a more recent majhong. But dont let his mustache fool you as he is a wellversed master of mechanics. Having served in the military prior to his experiences in Davis, he has quite a few stories to tell.

A favorite band of his to listen to is the White Stripes as can often be heard emanating from the office where he does his work.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c-----------------------------------------------------------------------
c This code was modified from a subroutine from the gamerf package at
c http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html .
c The original copyright statement states:
c   Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
c   You may use, copy, modify this code for any purpose and
c   without fee. You may distribute this ORIGINAL package.
c
c Permission to distribute this modified gamma function code
c with the FFTLog package has been granted
c (email from Takuya Ooura to Andrew Hamilton dated 16 March 1999).
c
c Original gamerf2a.doc documentation states:
c
c   Gamma(z)=sqrt(2*pi)*(z+r)^(z-1/2)*exp(-z-r)
c            *(a_0
c            +a_1*(z-1)/z
c            +a_2*(z-1)*(z-2)/z/(z+1)
c            +a_3*(z-1)*(z-2)*(z-3)/z/(z+1)/(z+2)
c            +...)
c        a_n= f_n*(2*n)*(2*n-1)*(2*n-2)*...*(n+1)/1/2/3/.../n
c            -a_0*(2*n)*(2*n-1)*(2*n-2)*...*(n+1)/1/2/3/.../n
c            -a_1*(2*n)*(2*n-1)*(2*n-2)*...*(n+2)/1/2/3/.../(n-1)
c            -...
c            -a_(n-1)*(2*n)/1
c        f_n=1/sqrt(2*pi)*(1*2*3*...*n)*(n+1+r)^(-n-1/2)*exp(n+1+r)

c   C.Lanczos,A Precision Approximation of the Gamma Function,
c   J.SIAM Numer.Anal.Ser.B,Vol.1,1964
c
c Modified 28 Oct 98 by Andrew J S Hamilton
c http://casa.colorado.edu/~ajsh/
c (1) to return ln[Gamma(x)] with the correct phase,
c     as well as Gamma(x), and
c (2) to remain accurate for large absolute values of input x,
c     and for x near 0 or negative integers.
c-----------------------------------------------------------------------
      complex*16 function cdgamma(x,l)
      integer l
      complex*16 x
c *
c * Complex Gamma function in double precision.
c *
c * l = 0: Gamma(x)
c *     1: ln[Gamma(x)]
c *
      real*8 pi,pv,pu,pr,p1,p2,p3,p4,p5,p6,q1,q2,q3,q4,q5,q6
      parameter (
     &    pi = 3.14159265358979324d+00,
     &    pv = 7.31790632447016203d+00,
     &    pu = 3.48064577727581257d+00,
     &    pr = 3.27673720261526849d-02,
     &    p1 = 1.05400280458730808d+01,
     &    p2 = 4.73821439163096063d+01,
     &    p3 = 9.11395751189899762d+01,
     &    p4 = 6.62756400966213521d+01,
     &    p5 = 1.32280130755055088d+01,
     &    p6 = 2.93729529320536228d-01)
      parameter (
     &    q1 = 9.99999999999975753d-01,
     &    q2 = 2.00000000000603851d+00,
     &    q3 = 2.99999999944915534d+00,
     &    q4 = 4.00000003016801681d+00,
     &    q5 = 4.99999857982434025d+00,
     &    q6 = 6.00009857740312429d+00)
      real*8 big
      parameter (big=1.d20)
      real*8 t,ui,ur,vi,vr,wr,wi,xi,xr,yi,yr,zero
      data zero /0.d0/
c
      xr = dble(x)
      xi = dimag(x)
c---x = 0, -1, -2
      if (xr .eq. aint(xr) .and. xr .le. 0.d0 .and. xi .eq. 0.d0) then
c...Gamma
        if (l .eq. 0) then
          wr = xr / 2.d0
c   +Infinity at even negative integers
          if (wr .eq. aint(wr)) then
            yr = 1.d0 / zero
c   -Infinity at odd negative integers
          else
            yr = -1.d0 / zero
          endif
          yi = 0.d0
c...lnGamma
        elseif (l .eq. 1) then
c   real part is +Infinity
          yr = 1.d0 / zero
          yi = pi * aint(xr)
        endif
        goto 200
      endif
c---Re(x) < 1/2 : use reflection formula
      if (xr .lt. .5d0) then
          wr = 1.d0 - xr
          wi = -xi
      else
          wr = xr
          wi = xi
      end if
c---large |x|: keep only leading term of rational function
      t = wr * wr + wi * wi
      if (t .gt. big * big) then
c   Rational function v
        vr = wr / t + pr
        vi = - wi / t
c   ln(overall factor)
c   u = ln(x + pv) - 1
        yr = wr + pv
        ur = yr
        if (ur .lt. 0.d0) ur = - ur
        ui = wi
        if (ui .lt. 0.d0) ui = - ui
        if (ur.ge.ui) then
          t = wi / yr
          ur = log(ur) + log(1.d0 + t * t) / 2.d0 - 1.d0
        else
          t = yr / wi
          ur = log(ui) + log(1.d0 + t * t) / 2.d0 - 1.d0
        endif
        ui = atan2(wi, yr)
c---not large |x|
      else
c   u = u(x) = x + q6 = O(x)
        ur = wr + q6
c   v = v(x,u) = (x + q5) u = (x + q5)(x + q6) = O(x^2)
        vr = ur * (wr + q5) - wi * wi
        vi = wi * (wr + q5) + ur * wi
c   y = y(x,u,v) = p6 + p5 u + p4 v = p4 x^2 + ... = O(x^2)
        yr = p6 + (p5 * ur + p4 * vr)
        yi = p5 * wi + p4 * vi
c   u = u(x,v) = (x + q4) v = (x + q4)(x + q5)(x + q6) = O(x^3)
        ur = vr * (wr + q4) - vi * wi
        ui = vi * (wr + q4) + vr * wi
c   v = v(x,u) = (x + q3) u = (x + q3)(x + q4)...(x + q6) = O(x^4)
        vr = ur * (wr + q3) - ui * wi
        vi = ui * (wr + q3) + ur * wi
c   y = y(y,u,v) = y + p3 u + p2 v = p2 x^4 + ... = O(x^4)
        yr = yr + (p3 * ur + p2 * vr)
        yi = yi + (p3 * ui + p2 * vi)
c   u = u(x,v) = (x + q2) v = (x + q2)(x + q3)...(x + q6) = O(x^5)
        ur = vr * (wr + q2) - vi * wi
        ui = vi * (wr + q2) + vr * wi
c   v = v(x,u) = (x + q1) u = (x + q1)(x + q2)...(x + q6) = O(x^6)
        vr = ur * (wr + q1) - ui * wi
        vi = ui * (wr + q1) + ur * wi
c   Numerator
c   y = y(y,u,v) = y + p1 u + v = x^6 + ... = O(x^6)
c     = (x+q1)...(x+q6) + p1 (x+q2)...(x+q6) + ... + p5 (x+q6) + p6
        yr = yr + (p1 * ur + vr)
        yi = yi + (p1 * ui + vi)
c   Denominator
c   u = x v = x(x + q1)(x + q2)...(x + q6) = O(x^7)
        ur = vr * wr - vi * wi
        ui = vi * wr + vr * wi
c   t = |u|^2
        t = ur * ur + ui * ui
c   Rational function v = y u*/|u|^2 + pr = y/u + pr = pr + 1/x + ...
c     = pr + 1/x ( 1 + 1/(x+q1) ( p1 + 1/(x+q2) ( p2 + ...
        vr = (yr * ur + yi * ui) / t + pr
        vi = (yi * ur - yr * ui) / t
c   Overall factor
c   u = ln(x + pv) - 1
        yr = wr + pv
        ur = log(yr * yr + wi * wi) / 2.d0 - 1.d0
        ui = atan2(wi, yr)
      endif
c---lnGamma
c   y = u(x - .5) - pu
c     = (x - .5) [ln(x + pv) - 1] - pu
      yr = ur * (wr - 0.5d0) - ui * wi - pu
      yi = ui * (wr - 0.5d0) + ur * wi
c   y = y + ln(v)
c     = (x - .5) [ln(x + pv) - 1] - pu + ln(Rational)
c     = lnGamma(x)
      yr = yr + log(vr * vr + vi * vi) / 2.d0
      yi = yi + atan2(vi, vr)
c---Reflection formula Gamma(x) Gamma(1-x) = pi/sin(pi x)
c   sign of Gamma
      t = 1.d0
      if (xr .lt. .5d0) then
        wi = anint(xr)
        wr = xr - wi
        if (wi .gt. xr) wi = wi - 1.d0
c   case of real x
        if (xi .eq. 0.d0) then
c   w = ln[sin(pi x)]
          wr = log(sin(pi * abs(wr)))
          if (l .eq. 0) then
            if (wi .ne. 2.d0 * aint(wi / 2)) t = -1.d0
            wi = 0.d0
          elseif (l .eq. 1) then
            wi = - pi * wi
          endif
c   case where imaginary part of x is < 1 in absolute value
        elseif (abs(xi) .lt. 1.d0) then
          if (l .eq. 0) then
            if (wi .ne. 2.d0 * aint(wi / 2.d0)) t = -1.d0
            ui = 0.d0
          elseif (l .eq. 1) then
            ui = -pi * wi
            if (xi .lt. 0.d0) ui = -ui
          endif
          wr = pi * wr
          wi = pi * xi
          vr = sin(wr) * cosh(wi)
          vi = cos(wr) * sinh(wi)
          if (wr .lt. 0.d0) then
            vr = -vr
            vi = -vi
          endif
c   w = ln[sin(pi x)]
          wr = log(vr * vr + vi * vi) / 2.d0
          wi = ui + atan2(vi, vr)
c   case where imaginary part of x is >= 1 in absolute value
        else
          if (l .eq. 0) then
            if (wi .ne. 2.d0 * aint(wi / 2)) t = -1.d0
            if (wr .ge. 0.d0) then
              ui = pi * (.5d0 - wr)
            else
              ui = pi * (- .5d0 - wr)
            endif
          elseif (l .eq. 1) then
            ui = pi * (.5d0 - xr)
          endif
          wi = exp(- 2.d0 * pi * abs(xi))
          wr = 2.d0 * pi * wr
          vr = (1.d0 - cos(wr) * wi) / 2.d0
          vi = - sin(wr) * wi / 2.d0
          ur = pi * xi
c   w = ln[sin(pi x)]
          if (xi .gt. 0.d0) then
            wr = ur + log(vr * vr + vi * vi) / 2.d0
            wi = ui + atan2(vi, vr)
          elseif (xi .lt. 0.d0) then
            wr = - ur + log(vr * vr + vi * vi) / 2.d0
            wi = - ui - atan2(vi, vr)
          endif
        endif
c   y = ln[Gamma(x)]
        yr = log(pi) - yr - wr
        yi = - yi - wi
      endif
c---Gamma
      if (l .eq. 0) then
        ur = exp(yr)
        if (xi .eq. 0.d0) then
          yr = t * ur
          yi = 0.d0
        else
          yr = t * ur * cos(yi)
          yi = ur * sin(yi)
        endif
      endif
c---finish
  200 cdgamma = dcmplx(yr, yi)
      end
c

max_stars_repo_name: sfarrens/cosmostat


text:
      subroutine postele20(ix, ele, ntip_edge)

	  parameter (nume=40000)
	  parameter (nume2=20000)
	  common/bk06/nprnt,mprint,itmpop,numelt,jprint,idump,locstr
	  common /pcracktip/ connect(4,nume2), node(2,nume2), &
     	       penta(nume2),ndflag(2,nume2), numnpt, numeltu, ndc

	  dimension ix(4,*)
      integer ix

	  integer numelt
	  integer connect, numeltu, penta
	  integer ele, ntip_edge, ndc

	  if(ntip_edge==1) then
	      call cracknode1(1,ele,ix)

		  numeltu=numeltu+1
		  connect(1,numeltu)=ndc
		  connect(2,numeltu)=connect(2,ele)
		  connect(3,numeltu)=connect(3,ele)
		  connect(4,numeltu)=connect(3,ele)
		  penta(numeltu-numelt)=ele

		  connect(2,ele)=ndc

	  else if(ntip_edge==2) then
	      call cracknode1(2,ele,ix)

		  numeltu=numeltu+1
		  connect(1,numeltu)=ndc
		  connect(2,numeltu)=ndc
		  connect(3,numeltu)=connect(3,ele)
		  connect(4,numeltu)=connect(4,ele)
		  penta(numeltu-numelt)=ele

		  connect(3,ele)=ndc

	  else if(ntip_edge==3) then
	      call cracknode1(3,ele,ix)

		  numeltu=numeltu+1
		  connect(1,numeltu)=connect(1,ele)
		  connect(2,numeltu)=ndc
		  connect(3,numeltu)=ndc
		  connect(4,numeltu)=connect(4,ele)
		  penta(numeltu-numelt)=ele

		  connect(4,ele)=ndc

	  else if(ntip_edge==4) then
	      call cracknode1(4,ele,ix)

		  numeltu=numeltu+1
		  connect(1,numeltu)=connect(1,ele)
		  connect(2,numeltu)=connect(2,ele)
		  connect(3,numeltu)=ndc
		  connect(4,numeltu)=ndc
		  penta(numeltu-numelt)=ele

		  connect(1,ele)=ndc

	  end if

      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


! This file was generated by gpufort

!
! Hints:
! Device variables in scope:
!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       character, intent(in)                                 :: uplo

!       complex(8), parameter                                 :: cone = cmplx(1.,0,8), chalf = cmplx(0.5,0,8)

!       real(8), parameter                                    :: one = 1.0_8

!       integer                                               :: i, j

!       integer                                               :: i, j

!       integer                                               :: k, kb, istat

!       integer                                               :: k, kb, istat

!       integer                                               :: k, kb, istat

function test_launch_krnl_185097_0_auto()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use zhegst_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  integer(c_int),intent(IN) :: sharedMem
  type(c_ptr),value,intent(IN) :: stream
  TODO declaration not found :: a
  integer,value :: k
  integer,value :: kb
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_krnl_185097_0_auto(0,c_null_ptr,sharedMem,stream,a,k,kb) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_krnl_185097_0_cpu(0,c_null_ptr,sharedMem,stream,a,k,kb)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

!
! Hints:
! Device variables in scope:
!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       integer, intent(in)                                   :: itype, n, lda, ldb, nb

!       character, intent(in)                                 :: uplo

!       complex(8), parameter                                 :: cone = cmplx(1.,0,8), chalf = cmplx(0.5,0,8)

!       real(8), parameter                                    :: one = 1.0_8

!       integer                                               :: i, j

!       integer                                               :: i, j

!       integer                                               :: k, kb, istat

!       integer                                               :: k, kb, istat

!       integer                                               :: k, kb, istat

function test_launch_krnl_9c3ecb_1_auto()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use zhegst_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  integer(c_int),intent(IN) :: sharedMem
  type(c_ptr),value,intent(IN) :: stream
  TODO declaration not found :: a
  integer,value :: k
  integer,value :: kb
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_krnl_9c3ecb_1_auto(0,c_null_ptr,sharedMem,stream,a,k,kb) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_krnl_9c3ecb_1_cpu(0,c_null_ptr,sharedMem,stream,a,k,kb)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

program test_zhegst_gpu_kernels
  implicit none
  integer :: globalErrorCode = 0, errorCode, fails = 0, tests = 0
  ! declare test functions and return type
  integer :: test_launch_krnl_185097_0_auto
  integer :: test_launch_krnl_9c3ecb_1_auto
  write(*,*) "SUITE test_zhegst_gpu_kernels run ..."
  errorCode = test_launch_krnl_185097_0_auto()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_krnl_185097_0_auto ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_krnl_185097_0_auto ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode
  errorCode = test_launch_krnl_9c3ecb_1_auto()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_krnl_9c3ecb_1_auto ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_krnl_9c3ecb_1_auto ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode

  IF (globalErrorCode > 0) THEN
    write(*,*) "SUITE test_zhegst_gpu_kernels ... FAILURE passed:",(tests-fails)," failed:",fails," total:",tests
  ELSE
    write(*,*) "SUITE test_zhegst_gpu_kernels ... SUCCESS passed:",(tests-fails)," failed:",fails," total:",tests
  END IF
end program test_zhegst_gpu_kernels


text:
*DK kDtreeselect
      subroutine kDtreeselect(ielmtyp, xicelm,yicelm,zicelm, linktree,
     &                        safebox, numfound, itfound, iflag)
C
C######################################################################
C
C     PURPOSE -
C
C        The Purpose of this subroutine is to take a k-D-R+ tree found
C        in linktree and provide an entry point into the routines that
C        find candidates based on the type of element (currently only
C        two)
C
C     NOTES -
C
C
C     INPUT ARGUMENTS -
C        ielmtyp                    The type of intersecting element,
C                                   used to decide what routine to use.
C        xicelm(),yicelm(),zicelm() Coordinates of the intersecting
C                                   element.
C        linktree()                 k-D-R+ tree structure
C        safebox(2,3,*)             bounding boxes corresponding to
C                                   elements in the k-D-R+ tree
C
C     OUTPUT ARGUMENTS -
C        numfound                   number of candidate elements found
C        itfound                    array of candidate elements
C        iflag                      error flag, 0 if successful.
C
C
C     $Log: kdtreeselect.f,v $
C     Revision 2.00  2007/11/05 19:45:59  spchu
C     Import to CVS
C
CPVCS
CPVCS       Rev 1.3   28 Mar 2000 14:09:20   dcg
CPVCS    remove include 'machine.h'
CPVCS
CPVCS       Rev 1.2   22 Mar 2000 08:39:04   dcg
CPVCS    use local_epsilon in place of epsilon
CPVCS
CPVCS       Rev 1.1   08 Feb 2000 08:38:24   dcg
CPVCS    remove comdict
CPVCS
CPVCS       Rev 1.0   Wed Aug 04 10:00:34 1999   bap
CPVCS    Initial revision.
C
C
C######################################################################
C
C     Variable Declarations
C
      implicit none
      include "local_element.h"
C
C     Arguments
      integer ielmtyp
      real*8 xicelm(1000000), yicelm(1000000), zicelm(1000000)
      integer linktree(1000000)
      real*8 safebox(2,3,1000000)
      integer numfound
      integer itfound(1000000)
      integer iflag
C
C     Indices
      integer i
C
C     Bounds
      real*8 xboundhi, yboundhi, zboundhi
      real*8 xboundlo, yboundlo, zboundlo
C
C     Other stuff
      real*8 local_epsilon
      real*8 releps

C
C######################################################################
C
C     Initialize the variables
C
      local_epsilon = 1.0e-10
C
C######################################################################
C
      if (ielmtyp.eq.ifelmlin) then
         releps = sqrt((xicelm(2)-xicelm(1))**2+
     &                 (yicelm(2)-yicelm(1))**2+
     &                 (zicelm(2)-zicelm(1))**2)*local_epsilon
         call lineseg_inter(xicelm(1),yicelm(1),zicelm(1),
     &                      xicelm(2),yicelm(2),zicelm(2),
     &                      linktree, safebox, releps,
     &                      numfound, itfound, iflag)
      else
         xboundhi = xicelm(1)
         xboundlo = xicelm(1)
         yboundhi = yicelm(1)
         yboundlo = yicelm(1)
         zboundhi = zicelm(1)
         zboundlo = zicelm(1)
         do i = 2,nelmnen(ielmtyp)
            xboundhi = max(xicelm(i),xboundhi)
            xboundlo = min(xicelm(i),xboundlo)
            yboundhi = max(yicelm(i),yboundhi)
            yboundlo = min(yicelm(i),yboundlo)
            zboundhi = max(zicelm(i),zboundhi)
            zboundlo = min(zicelm(i),zboundlo)
         enddo
         call bbox_inter(xboundlo, yboundlo, zboundlo,
     &                   xboundhi, yboundhi, zboundhi,
     &                   linktree, safebox,
     &                   numfound, itfound, iflag)
      endif
      return
      end

*DK bbox_inter
      subroutine bbox_inter(xboundlo, yboundlo, zboundlo,
     &                      xboundhi, yboundhi, zboundhi,
     &                      linktree, safebox,
     &                      numfound, itfound, iflag)
C
C######################################################################
C
C     PURPOSE -
C
C        The Purpose of this subroutine is to take a k-D-R+ tree found
C        in linktree and find candidates based on the bounding boxes
C        provided in safebox and defined by <boundlo> and <boundhi>
C
C     NOTES -
C
C        This code is currently NONFUNCTIONAL
C
C     INPUT ARGUMENTS -
C        xboundlo,yboundlo,zboundlo Lower Bounding Box Coordinates of
C                                   the intersecting element.
C        xboundhi,yboundhi,zboundhi Upper Bounding Box Coordinates of
C                                   the intersecting element.
C        linktree()                 k-D-R+ tree structure
C        safebox(2,3,*)             bounding boxes corresponding to
C                                   elements in the k-D-R+ tree
C
C     OUTPUT ARGUMENTS -
C        numfound                   number of candidate elements found
C        itfound                    array of candidate elements
C        iflag                      error flag, 0 if successful.
C
C######################################################################
C
C     Variable Declarations
C
      implicit none
C
C     Arguments
      real*8 xboundlo,yboundlo,zboundlo
      real*8 xboundhi,yboundhi,zboundhi
      integer linktree(1000000)
      real*8 safebox(2,3,1000000)
      integer numfound
      integer itfound(1000000)
      integer iflag
C
C     Indices
      integer i, j

C
C     Stack Variables
      integer stack1(500)
      integer top, node, index
C
C     Useful max and min arrays
      real*8 maxcoord(3)
      real*8 mincoord(3)
C
C     Epsilon factor things
      real*8 releps
C
C######################################################################
C
C     Initialize some important variables...
C
      numfound = 0
      releps = 0.0
C
C     *****************************************************************
C     Put the maxes & mins in to tbe arrays
      maxcoord(1)=xboundhi
      maxcoord(2)=yboundhi
      maxcoord(3)=zboundhi
      mincoord(1)=xboundlo
      mincoord(2)=yboundlo
      mincoord(3)=zboundlo
C
C######################################################################
C
C     Ensure that we have a small, but valid, bounding box. (Find the
C     proper slop factor for the bounding boxes)
C
      do i=1,3
         releps = max(maxcoord(i)-mincoord(i),releps)
      enddo
      if (releps.eq.0) then
         releps = 1.0e-10
      else
         releps = 1.0e-10*releps
      endif
C
C######################################################################
C
C     Now the fun begins...
C
C     *****************************************************************
C     Is the element within the bounding box of the kdtree?
      do i=1,3
         if (mincoord(i).gt.(safebox(2,i,1)+releps)) goto 9999
         if (maxcoord(i).lt.(safebox(1,i,1)-releps)) goto 9999
      enddo
C
C     *****************************************************************
C     Is the root node a leaf?
      if (linktree(1).lt.0) then
         numfound = 1
         itfound(1)= -linktree(1)
         goto 9999
      endif
C
C     *****************************************************************
C     Prep the stack
      top=1
      stack1(top)=1
C
C     *****************************************************************
C     Traverse the k-D tree via the stack
      do while (top.gt.0)
         node=stack1(top)
         top = top-1
         index=linktree(node)
C
C     *****************************************************************
C     See if node's children intersect the bounding box of the element
         do j=0,1
C
C     *****************************************************************
C     If they're leaves, then the best we can do is add them to the
C     candidate list, (There's no way to check the boxes of these
C     nodes without working with the itet arrays)
            if (linktree(index+j).lt.0) then
               numfound=numfound+1
               itfound(numfound) = -linktree(index+j)
C
C     *****************************************************************
C     Otherwise, if they are branch nodes, figure out if they are
C     within the bounding box.
            else
               do i=1,3
                  if (mincoord(i).gt.(safebox(2,i,index+j)+releps))
     &                 goto 10
                  if (maxcoord(i).lt.(safebox(1,i,index+j)-releps))
     &                 goto 10
               enddo
               top = top+1
               stack1(top) = index+j
 10            continue
            endif
         enddo
      enddo
 9999 continue
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
      SUBROUTINE GRATSI(DIX,DIY)
C
C          ------------------------------------------------
C          ROUTINE NO. ( 221)   VERSION (A8.6)    24:NOV:86
C          ------------------------------------------------
C
C          THIS DRAWS A GRATICULE WITH NO ANNOTATION,
C          WITH INTERVALS SET BY <DIX> AND <DIY>.
C
C
C          [X-LINEAR, Y-LINEAR MARKING].
C
C
      REAL    EDGE(2)
      REAL    RDATA(1)
      INTEGER IDATA(1)
C
      COMMON /T0ADIX/ DIVLX,NTIKLX,NTIKHX
      COMMON /T0ADIY/ DIVLY,NTIKLY,NTIKHY
      COMMON /T0NOTA/ NOTATA
      COMMON /T0PPOS/ XPLOT0,YPLOT0
      COMMON /T0TRAC/ IPRINT
      COMMON /T0TRAR/ RTRAC1,RTRAC2,RTRAC3,RTRAC4
      COMMON /T0WNDO/ X1WND0,X2WND0,Y1WND0,Y2WND0
      COMMON /T3LIMS/ IMAXI,RMAXI,RMINI
C
      DATA RDATA /0.0/
C
C
      RTRAC1= DIX
      RTRAC2= DIY
      IF (IPRINT.EQ.1) CALL G0MESG(99,2)
C
C          IF THE WINDOW AREA IS WRONG, NOTHING MORE IS DONE.
C          THE ANNOTATION MARKER IS SAVED AND THEN SET = 0
C          SO THAT ANNOTATION IS NOT DONE WITH A GRATICULE.
C          THE CURRENT PLOTTING POSITION IS THEN SAVED.
C
      IF (ABS(X1WND0-X2WND0).LT.RMINI.OR.
     &    ABS(Y1WND0-Y2WND0).LT.RMINI) RETURN
C
      NTATSV= NOTATA
      NOTATA= 0
      XHERE= XPLOT0
      YHERE= YPLOT0
C
C          THE X-AXIS AND Y-AXIS GRATICULE LINES ARE DRAWN.
C
      IDATA(1)= 1
      CALL G3LINK(0,14,-1,IDATA,RDATA)
      CALL G0DIVS(1,DIX)
      NTICKS= NTIKHX-NTIKLX+1
      IF (NTICKS.LE.0) GO TO 1
C
      EDGE(1)= Y1WND0
      EDGE(2)= Y2WND0
      INDEDG= 1
C
      DO 100 ITICK= 1,NTICKS
        XPOS= DIVLX*(NTIKLX+ITICK-1)
        CALL POSITN(XPOS,EDGE(INDEDG))
        INDEDG= 3-INDEDG
        CALL JOIN(XPOS,EDGE(INDEDG))
  100 CONTINUE
C
    1 CALL G0DIVS(2,DIY)
      NTICKS= NTIKHY-NTIKLY+1
      IF (NTICKS.LE.0) GO TO 2
C
      EDGE(1)= X1WND0
      EDGE(2)= X2WND0
      INDEDG= 1
C
      DO 200 ITICK= 1,NTICKS
        YPOS= DIVLY*(NTIKLY+ITICK-1)
        CALL POSITN(EDGE(INDEDG),YPOS)
        INDEDG= 3-INDEDG
        CALL JOIN(EDGE(INDEDG),YPOS)
  200 CONTINUE
C
    2 IDATA(1)= 0
      CALL G3LINK(0,14,-1,IDATA,RDATA)
C
C          THE ANNOTATION MARKER AND THE CURRENT
C          PLOTTING POSITION ARE THEN RESTORED.
C
      NOTATA= NTATSV
      CALL POSITN(XHERE,YHERE)
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
      SUBROUTINE ZEROF(F,DX,DFDX)
*
*-----------------------------------------------------------------------
*
* FIND DX=-F/DFDX, TO MAKE F ZERO. IF DX=0 AT ENTRY, THEN DFDX IS A
* START APPROXIMATION, AND IT IS THE FIRST CALL. OTHERWISE USE OLD
* INFO.  780926/NORDLUND.
*
* Export version  1988-03-24  ********* Olof Morell *** Uppsala
*
*-----------------------------------------------------------------------
*
      IF(DX.NE.0.) DFDX=(F-FOLD)/DXOLD
      DX=-F/DFDX
      FOLD=F
      DXOLD=DX
*
      RETURN
      END

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
c
c     file tcud24.f
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1999 by UCAR                 .
c  .                                                             .
c  .       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                      MUDPACK version 5.0                    .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... author and specialist
c
c          John C. Adams (National Center for Atmospheric Research)
c          email: johnad@ucar.edu, phone: 303-497-1213

c ... For MUDPACK 5.0 information, visit the website:
c     (http://www.scd.ucar.edu/css/software/mudpack)
c
c
c ... purpose
c
c     test program for the complex MUDPACK solver cud2
c
c ... required MUDPACK files
c
c     cud2.f, cudcom.f
c
c
c *********************************************************
c *********************************************************
c
c     sample program/test driver for cud2
c
c **********************************************************
c **********************************************************
c
c     a sample program/test driver for cud2 is below.  it can be
c     executed as an initial test.  the output is listed for the
c     test case described.
c
c     test the driver below by solving the complex linear elliptic pde
c
c        cmplx(1.+x*x,1.+y*y)*pxx +
c
c        cmplx(exp(-x),exp(-y))*(pyy - py) +
c
c        cmplx(y,x)*p(x,y) = r(x,y)
c
c     on a 49 X 65 grid superimposed on the unit square with
c     specified boundary conditions at xb=1.0, yc=1.0 and
c     mixed boundary conditions
c
c           dp/dx - cmplx(y,y)*p(xa,y) = g(y) at x = xa
c     and
c
c           dp/dy + cmplx(x,x)*p(x,yd) = h(x) at y = yd.
c
c     the exact solution
c
c          p(x,y) = cmplx(x**5,y**5) + 1.0
c
c     is used for testing.  three multigrid cycles (no initial guess)
c     with red/black gauss-seidel point relaxation and the default multigrid
c     options are used to compute the second-order approximation with cud2.
c     Then cud24 is called to produce a fourth order estimate.
c
c ******************************************************
c     output (64 bit floating point arithmetic)
c *******************************************************
c
c     cud2 test
c
c     integer input arguments
c     intl =  0 nxa =  2 nxb =  1 nyc =  1 nyd =  2
c     ixp =  3 jyq =  2 iex =  5 jey =  6
c     nx =  49 ny =  65 iguess =  0 maxcy =  3
c     method =  0 work space estimate =   34192
c
c     multigrid option arguments
c     kcycle =  2
c     iprer =  2
c     ipost =  1
c     intpol =  3
c
c     floating point input parameters
c     xa =  0.000 xb =  1.000 yc =  0.000 yd =  1.000
c     tolerance (error control) =    0.000E+00
c
c     discretization call to cud2 intl =  0
c     ierror =  0 minimum work space =   34192
c
c     approximation call to cud2
c     intl =  1 method =  0 iguess =  0 maxcy =  3
c     ierror =  0
c     maximum error  =   0.640E-03
c
c     cud24 test  ierror =  0
c     maximum error  =   0.346E-05
c
c ************************************************************
c     end of output
c ************************************************************
c
      program tcud24
      implicit none
      integer iixp,jjyq,iiex,jjey,nnx,nny,llw
c
c     set grid sizes with parameter statements
c
      parameter (iixp = 3 , jjyq = 2 , iiex =5, jjey = 6)
      parameter (nnx=iixp*2**(iiex-1)+1, nny=jjyq*2**(jjey-1)+1)
c
c     set exact minimal work space required (see tcud2.f)
c
      parameter (llw=34192)
c
c     dimension solution,right hand side, and work arrays
c
      complex p(nnx,nny),r(nnx,nny),w(llw)
c     put integer and floating point parameter names in contiguous
c     storeage for labelling purposes
c
      integer iprm(16),mgopt(4)
      real fprm(6)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,
     +              iguess,maxcy,method,nwork,lwrkqd,itero
      common/itcud2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,
     +              iguess,maxcy,method,nwork,lwrkqd,itero
      real xa,xb,yc,yd,tolmax,relmax
      common/ftcud2/xa,xb,yc,yd,tolmax,relmax
      equivalence(intl,iprm)
      equivalence(xa,fprm)
      integer i,j,ierror
      complex cxx,cyy,cx,cy,ce,pxx,pyy,px,py,pe
      real dlx,dly,x,y,errmax
c
c     declare coefficient and boundary condition input subroutines external
c
      external cof,bndc
c
c
c     set input integer arguments
c
      intl = 0
c
c     set boundary condition flags
c
      nxa = 2
      nxb = 1
      nyc = 1
      nyd = 2
c
c     set grid sizes from parameter statements
c
      ixp = iixp
      jyq = jjyq
      iex = iiex
      jey = jjey
      nx = nnx
      ny = nny
c
c     set for three multigrid cycle
c
      maxcy = 3
c
c     set work space length approximation from parameter statement
c
      nwork = llw
c
c     set point relaxation
c
      method = 0
c
c     flag no initial guess (this sets full multigrid cycling)
c
      iguess = 0
c
c     set end points of solution rectangle in (x,y) space
c
      xa = 0.0
      xb = 1.0
      yc = 0.0
      yd = 1.0
c
c     set mesh increments
c
      dlx = (xb-xa)/float(nx-1)
      dly = (yd-yc)/float(ny-1)
c
c     set for no error control
c
      tolmax = 0.0
c
c     set right hand side in r
c     initialize p to zero
c
      do i=1,nx
	x = xa+float(i-1)*dlx
	do j=1,ny
	  y = yc+float(j-1)*dly
	  call cof(x,y,cxx,cyy,cx,cy,ce)
	  call exact(x,y,pxx,pyy,px,py,pe)
	  r(i,j) = cxx*pxx+cyy*pyy+cx*px+cy*py+ce*pe
	  p(i,j) = (0.0,0.0)
	end do
      end do
c
c     set specified boundaries in p
c
      x = xb
      do j=1,ny
	y = yc+float(j-1)*dly
	call exact(x,y,pxx,pyy,px,py,pe)
	p(nx,j) = pe
      end do
      y = yc
      do i=1,nx
	x = xa+float(i-1)*dlx
	call exact(x,y,pxx,pyy,px,py,pe)
	p(i,1) = pe
      end do
c
c     set default multigrid opitons
c
      mgopt(1) = 2
      mgopt(2) = 2
      mgopt(3) = 1
      mgopt(4) = 3
c
c     print input parameters (except multigrid options which are default)
c
      write(6,100)
  100 format(//' cud2 test ')
      write (*,101) (iprm(i),i=1,15)
  101 format(/' integer input arguments ',
     +/'intl = ',i2,' nxa = ',i2,' nxb = ',i2,' nyc = ',i2,' nyd = ',i2,
     +/' ixp = ',i2,' jyq = ',i2,' iex = ',i2,' jey = ',i2
     +/' nx = ',i3,' ny = ',i3,' iguess = ',i2,' maxcy = ',i2,
     +/' method = ',i2, ' work space estimate = ',i7)
      write (*,102) (mgopt(i),i=1,4)
  102 format(/' multigrid option arguments ',
     +/' kcycle = ',i2,
     +/' iprer = ',i2,
     +/' ipost = ',i2
     +/' intpol = ',i2)
      write(*,103) xa,xb,yc,yd,tolmax
  103 format(/' floating point input parameters ',
     +/' xa = ',f6.3,' xb = ',f6.3,' yc = ',f6.3,' yd = ',f6.3,
     +/' tolerance (error control) =   ',e10.3)
c
c     intiialization call
c
      write(*,104) intl
  104 format(/' discretization call to cud2', ' intl = ', i2)
      call cud2(iprm,fprm,w,cof,bndc,r,p,mgopt,ierror)
c
c     print error parameter and minimum work space requirement
c
      write (*,105) ierror,iprm(16)
  105 format(' ierror = ',i2, ' minimum work space = ',i7)
      if (ierror.gt.0) call exit(0)
c
c     attempt solution
c
      intl = 1
      write(*,106) intl,method,iguess,maxcy
  106 format(/' approximation call to cud2',
     +/' intl = ',i2, ' method = ',i2,' iguess = ',i2, ' maxcy = ',i2)
      call cud2(iprm,fprm,w,cof,bndc,r,p,mgopt,ierror)
      write (*,107) ierror
  107 format(' ierror = ',i2)
c
c     compute and print exact maximum error
c
      errmax = 0.0
      do j=1,ny
	y = yc+(j-1)*dly
	do i=1,nx
	  x = xa+(i-1)*dlx
	  call exact(x,y,pxx,pyy,px,py,pe)
	  errmax = amax1(errmax,cabs((p(i,j)-pe)))
	end do
      end do
      write(*,108) errmax
  108 format(' maximum error  =  ',e10.3)
c
c      attempt to improve approximation to fourth order
c
      call cud24(w,p,ierror)
      write (*,109) ierror
  109 format(/' cud24 test ', ' ierror = ',i2)
      if (ierror.gt.0) call exit(0)
      errmax = 0.0
      do j=1,ny
	y = yc+(j-1)*dly
	do i=1,nx
	  x = xa+(i-1)*dlx
	  call exact(x,y,pxx,pyy,px,py,pe)
	  errmax = amax1(errmax,abs((p(i,j)-pe)))
	end do
      end do
      write(*,108) errmax
      end

      subroutine cof(x,y,cxx,cyy,cx,cy,ce)
c
c     input pde coefficients at any grid point (x,y) in the solution region
c     (xa.le.x.le.xb,yc.le.y.le.yd) to cud2
c
      implicit none
      real x,y
      complex cxx,cyy,cx,cy,ce
      cxx = cmplx(1.+x*x,1.+y*y)
      cyy = cmplx(exp(-x),exp(-y))
      cx = (0.,0.)
      cy = -cyy
      ce = -cmplx(y,x)
      return
      end

      subroutine bndc(kbdy,xory,alfa,gbdy)
c
c     input mixed derivative b.c. to cud2
c
      implicit none
      integer kbdy
      real xory,x,y
      complex alfa,gbdy
      real xa,xb,yc,yd,tolmax,relmax
      common/ftcud2/xa,xb,yc,yd,tolmax,relmax
      complex pe,px,py,pxx,pyy
      if (kbdy.eq.1) then
c
c     x=xa boundary (nxa must equal 2)
c     b.c. has the form px + alfxa(y)*pe = gbdxa(y)
c     alfa and gbdy corresponding to alfxa(y),gbdxa(y)
c     must be output
c
	y = xory
	x = xa
	call exact(x,y,pxx,pyy,px,py,pe)
	alfa = -cmplx(y,y)
	gbdy = px + alfa*pe
	return
      end if
      if (kbdy.eq.4) then
c
c     y = yd boundary (nyd must equal 2)
c     b.c. has the form py + alfyd(x)*pe = gbdyd(x)
c     alfa and gbdy corresponding to alfyd(x),gbdyd(x)
c     must be output
c
	y = yd
	x = xory
	call exact(x,y,pxx,pyy,px,py,pe)
	alfa = cmplx(x,x)
	gbdy = py + alfa*pe
      return
      end if
      end

      subroutine exact(x,y,pxx,pyy,px,py,pe)
c
c     this subroutine is used to set an exact solution for testing cud2
c
      implicit none
      real x,y
      complex pxx,pyy,px,py,pe
      pe = cmplx(x**5,y**5)+1.0
      px = cmplx(5*x**4,0.)
      py = cmplx(0.,5*y**4)
      pxx = cmplx(20.*x**3,0.)
      pyy = cmplx(0.,20.*y**3)
      return
      end

max_stars_repo_name: utastudents/selalib

max_issues_repo_name: utastudents/selalib

max_forks_repo_name: utastudents/selalib


text:
!=========================================================
subroutine rp1(maxmx, meqn, mwaves, maux, mbc, mx, ql, qr, auxl, auxr, wave, s, amdq, apdq)
!=========================================================
!
! Roe-solver for the 1D ideal MHD equations
!
! waves:     7
! equations: 8
!
! Conserved quantities:
!       1 density
!       2 momentum, x
!       3 momentum, y
!       4 momentum, z
!       5 energy
!       6 magnetic field, x
!       7 magnetic field, y
!       8 magnetic field, z
!
! On input, ql contains the state vector at the left edge of each cell
!           qr contains the state vector at the right edge of each cell
!
! On output, wave contains the waves,
!            s the speeds,
!            amdq the  left-going flux difference  A^- \Delta q
!            apdq the right-going flux difference  A^+ \Delta q
!            (where the flux difference is f(qr(i-1)) - f(ql(i)))
!
! With the Roe solver we have
!   amdq  =  A^- \Delta q    and    apdq  =  A^+ \Delta q
! where A is the Roe matrix.
! An entropy fix can also be incorporated into the flux differences.
!
! Note that the i'th Riemann problem has left state qr(:, i-1)
!                                    and right state ql(:, i)
! From the basic clawpack routine step1, rp is called with ql = qr = q.

    implicit none

    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
    double precision, dimension(meqn, 1-mbc:maxmx+mbc),         intent(in)  :: ql, qr
    double precision, dimension(maux, 1-mbc:maxmx+mbc),         intent(in)  :: auxl, auxr
    double precision, dimension(meqn, mwaves, 1-mbc:maxmx+mbc), intent(out) :: wave
    double precision, dimension(meqn, 1-mbc:maxmx+mbc),         intent(out) :: amdq, apdq
    double precision, dimension(mwaves, 1-mbc:maxmx+mbc),       intent(out) :: s

    ! meqn == 8, mwaves == 7
    double precision, dimension(meqn)        :: fl, fr, df
    double precision, dimension(4)           :: eigl, eigr
    double precision, dimension(mwaves)      :: stmp
    double precision, dimension(meqn,mwaves) :: rr
    double precision, dimension(mwaves,meqn) :: rl
    double precision, dimension(mwaves)      :: alpha, beta

    integer ixy, i, mu1, mu2, mb1, mb2, m1, m2, mm, mflag
    double precision gamma, gamma1
    double precision rho, rhol, rhor, u1, u1l, u1r, u2, u2l, u2r, u3, u3l, u3r, El, Er, p, pl, pr
    double precision B1, B1l, B1r, B2, B2l, B2r, B3, B3l, B3r
    double precision sl, sr, qs

    logical efix

    common /cparam/ gamma, gamma1

    ! use entropy fix for transonic rarefactions
    data efix /.true./

    ixy = 1
    if (ixy == 1) then
        mu1 = 2
        mu2 = 3
        mb1 = 6
        mb2 = 7
    elseif (ixy == 2) then
        mu1 = 3
        mu2 = 2
        mb1 = 7
        mb2 = 6
    endif

    do i = 2-mbc, mx+mbc
        ! Left and right states
        rhol = qr(1, i-1)
        rhor = ql(1, i  )

        ! There is an issue with uninitialized states at the far left and right ghost cells for sharpclaw.
        ! David Ketcheson will open an issue in pyclaw
        if (i > 3-mbc .and. i < mx+mbc-1 .and. (rhol <= 1.d-15 .or. isnan(rhol) .or. rhor <= 1.d-15 .or. isnan(rhor))) then
            print*, ' '
            print*, ' ERROR in RPN:  rhol = ', rhol
            print*, ' ERROR in RPN:  rhor = ', rhor
            print*, '                Cell = ', i
            print*, '             maxmx   = ', maxmx
            print*, '             meqn    = ', meqn
            print*, '             mwaves  = ', mwaves
            print*, '             maux    = ', maux
            print*, '             mbc     = ', mbc
            print*, '             mx      = ', mx
            print*, ' '
            stop
        endif

        u1l = qr(mu1, i-1) / rhol
        u2l = qr(mu2, i-1) / rhol
        u3l = qr(4,   i-1) / rhol
        El  = qr(5,   i-1)
        B1l = qr(mb1, i-1)
        B2l = qr(mb2, i-1)
        B3l = qr(8,   i-1)
        pl  = gamma1 * (El - 0.5d0 * rhol * (u1l**2 + u2l**2 + u3l**2) - 0.5d0 * (B1l**2 + B2l**2 + B3l**2))

        u1r = ql(mu1, i  ) / rhor
        u2r = ql(mu2, i  ) / rhor
        u3r = ql(4,   i  ) / rhor
        Er  = ql(5,   i  )
        B1r = ql(mb1, i  )
        B2r = ql(mb2, i  )
        B3r = ql(8,   i  )
        pr  = gamma1 * (Er - 0.5d0 * rhor * (u1r**2 + u2r**2 + u3r**2) - 0.5d0 * (B1r**2 + B2r**2 + B3r**2))

        if (i > 3-mbc .and. i < mx+mbc-1 .and. (pl <= 1.d-15 .or. isnan(pl) .or. pr <= 1.d-15 .or. isnan(pr))) then
            print*, ' '
            print*, ' ERROR in RPN:  pl = ', pl
            print*, ' ERROR in RPN:  pr = ', pr
            print*, '                Cell = ', i
            print*, '             maxmx   = ', maxmx
            print*, '             meqn    = ', meqn
            print*, '             mwaves  = ', mwaves
            print*, '             maux    = ', maux
            print*, '             mbc     = ', mbc
            print*, '             mx      = ', mx
            print*, ' '
            stop
        endif

        ! Average states
        rho = 0.5d0 * (rhol + rhor)
        u1  = 0.5d0 * (u1l + u1r)
        u2  = 0.5d0 * (u2l + u2r)
        u3  = 0.5d0 * (u3l + u3r)
        p   = 0.5d0 * (pl + pr)
        B1  = 0.5d0 * (B1l + B1r)
        B2  = 0.5d0 * (B2l + B2r)
        B3  = 0.5d0 * (B3l + B3r)

        ! Delta f
        call set_flux_fun(mu1, mu2, mb1, mb2, rhol, u1l, u2l, u3l, pl, B1l, B2l, B3l, El, fl)
        call set_flux_fun(mu1, mu2, mb1, mb2, rhor, u1r, u2r, u3r, pr, B1r, B2r, B3r, Er, fr)
        do m1 = 1, meqn
            df(m1) = fr(m1) - fl(m1)
        enddo

        ! Wave Speeds
        do mm = 1, mwaves
            stmp(mm) = 0.d0
        enddo
        call set_wave_spd(rhol, u1l, u2l, u3l, pl, B1l, B2l, B3l, gamma, stmp)
        eigl(1) = stmp(1)
        eigl(2) = stmp(3)
        eigl(3) = stmp(5)
        eigl(4) = stmp(7)
        call set_wave_spd(rhor, u1r, u2r, u3r, pr, B1r, B2r, B3r, gamma, stmp)
        eigr(1) = stmp(1)
        eigr(2) = stmp(3)
        eigr(3) = stmp(5)
        eigr(4) = stmp(7)
        call set_wave_spd(rho, u1, u2, u3, p, B1, B2, B3, gamma, stmp)
        do m1 = 1, mwaves
          s(m1, i) = stmp(m1)
        enddo

        ! Right Eigenvectors (row,column)
        call set_rght_eig(mu1, mu2, mb1, mb2, rho, u1, u2, u3, p, B1, B2, B3, gamma, rr)

        ! Left Eigenvectors (row,column)
        call set_left_eig(mu1, mu2, mb1, mb2, rho, u1, u2, u3, p, B1, B2, B3, gamma, rl)

        ! Betas and Alphas
        do m1 = 1, mwaves
            beta(m1) = 0.d0
            do m2 = 1, meqn
                beta(m1) = beta(m1) + rl(m1, m2) * df(m2)
            enddo

            if (dabs(s(m1, i)) > 1.d-15) then
                alpha(m1) = beta(m1) / s(m1, i)
            else
                s(m1, i)  = 0.d0
                alpha(m1) = 0.d0
            endif
        enddo

        ! Waves
        do m1 = 1, mwaves
            do m2 = 1, meqn
                wave(m2, m1, i) = alpha(m1) * rr(m2, m1)
            enddo
        enddo

        ! compute left-going and right-going flux differences:
        !------------------------------------------------------
        !
        ! amdq = SUM s*wave   over left-going waves
        ! apdq = SUM s*wave   over right-going waves

        mflag = 0
        if (efix) then
            ! check if entropy fix is needed
            do mm = 1, 4
                if (eigl(mm) * eigr(mm) < 0.d0) then
                    mflag = 1
                endif
            enddo
        endif

        ! Solution update
        if (mflag == 0) then ! # WAVE PROPAGATION METHOD

            do m1 = 1, meqn
                amdq(m1, i) = 0.d0
                apdq(m1, i) = 0.d0

                do m2 = 1, mwaves
                    if (s(m2, i) < -1.d-15) then
                        amdq(m1, i) = amdq(m1, i) + s(m2, i) * wave(m1, m2, i)
                    elseif (s(m2, i) > 1.d-15) then
                        apdq(m1, i) = apdq(m1, i) + s(m2, i) * wave(m1, m2, i)
                    else
                        if (dabs(beta(m2)) > 1.d-15) then
                            print*, ' '
                            print*, ' In RPN: Conservation fix '
                            print*, ' '
                            print*, ' beta(m2) = ', beta(m2)
                            read*
                        endif
                        amdq(m1, i) = amdq(m1, i) + 0.5d0*beta(m2)*rr(m1,m2)
                        apdq(m1, i) = apdq(m1, i) + 0.5d0*beta(m2)*rr(m1,m2)
                    endif
                enddo
            enddo

        else !# HLL ENTROPY FIX

            sl = eigl(1)
            sr = eigr(4)

            do m1 = 1, meqn
                if (dabs(sl-sr) < 1.d-15) then
                    amdq(m1, i) = 0.d0
                    apdq(m1, i) = 0.d0
                else
                    qs = (df(m1) + sl * qr(m1, i-1) - sr * ql(m1, i)) / (sl - sr)

                    amdq(m1, i) = sl * (qs - qr(m1, i-1))
                    apdq(m1, i) = sr * (ql(m1, i) - qs)
                endif

                do m2 = 1, mwaves
                    wave(m1, m2, i) = 0.d0
                enddo
            enddo

        endif

    enddo

    return
end

!=========================================================
subroutine set_flux_fun(mu1, mu2, mb1, mb2, rho, u1, u2, u3, p, B1, B2, B3, E, f)
!=========================================================

    implicit double precision (a-h,o-z)
    double precision, dimension(8), intent(out) :: f

    Bm = 0.5d0 * (B1**2 + B2**2 + B3**2)
    Bu = u1 * B1 + u2 * B2 + u3 * B3

    f(1)   = rho * u1
    f(mu1) = rho * u1**2 + p + Bm - B1**2
    f(mu2) = rho * u1*u2 - B1*B2
    f(4)   = rho * u1*u3 - B1*B3
    f(5)   = u1 * (E + p + Bm) - B1 * Bu
    f(mb1) = 0.d0
    f(mb2) = u1*B2 - u2*B1
    f(8)   = u1*B3 - u3*B1

    return
end

!=========================================================
subroutine set_left_eig(mu1, mu2, mb1, mb2, rho, u1, u2, u3, p, B1, B2, B3, gamma, rl)
!=========================================================

    implicit none
    integer, intent(in) :: mu1, mu2, mb1, mb2
    double precision, intent(in) :: rho, u1, u2, u3, p, B1, B2, B3, gamma
    double precision, dimension(7,8), intent(out) :: rl
    double precision a2, a, d, ca, cf, cs, beta1, beta2, beta3, alphaf, alphas, g1, g2, um2, th1, th2

    a2 = gamma * p / rho
    a = dsqrt(a2)
    d = a2 + (B1**2 + B2**2 + B3**2) / rho

    ca = dsqrt(B1**2 / rho)
    cf = dsqrt(0.5d0 * ( d + dsqrt(d**2 - 4.d0 * a2 * B1 * B1 / rho)))
    cs = dsqrt(0.5d0 * ( d - dsqrt(d**2 - 4.d0 * a2 * B1 * B1 / rho)))

    beta1 = dsign(1.d0, B1)
    if ( (dabs(B2) + dabs(B3)) <= 1.d-15 ) then
        beta2 = 1.d0 / dsqrt(2.d0)
        beta3 = 1.d0 / dsqrt(2.d0)
    else
        beta2 = B2 / dsqrt(B2*B2 + B3*B3)
        beta3 = B3 / dsqrt(B2*B2 + B3*B3)
    endif

    if ( (dabs(B2) + dabs(B3)) <= 1.d-15 .and. dabs(a2 - B1*B1/rho) <= 1.d-15 ) then
        alphaf = 1.d0
        alphas = 1.d0
    else
        alphaf = dsqrt(dabs(cf * cf - ca*ca)) / dsqrt(dabs(cf * cf - cs * cs))
        alphas = dsqrt(dabs(cf * cf - a2   )) / dsqrt(dabs(cf * cf - cs * cs))
    endif

    g1 = gamma - 1.d0
    g2 = gamma - 2.d0

    um2 = u1**2 + u2**2 + u3**2
    th1 = alphaf*alphaf*a2*(cf*cf - (g2/g1)*a2) + alphas*alphas*cf*cf*(cs*cs - (g2/g1)*a2)
    th1 = 1.d0 / th1
    th2 = alphaf*alphaf*cf*a*beta1 + alphas*alphas*cs*ca*beta1
    th2 = 1.d0 / th2

    ! 1 - left eigenvector
    rl(1,1)   = th1*(alphaf/4.d0)*a2*um2 + 0.5d0*th2*(alphaf*a*u1*beta1 - alphas*cs*(beta2*u2 + beta3*u3))
    rl(1,mu1) = -0.5d0*th1*alphaf*a2*u1 - 0.5d0*th2*alphaf*a*beta1
    rl(1,mu2) = -0.5d0*th1*alphaf*a2*u2 + 0.5d0*th2*alphas*cs*beta2
    rl(1,4)   = -0.5d0*th1*alphaf*a2*u3 + 0.5d0*th2*alphas*cs*beta3
    rl(1,5)   = 0.5d0*th1*alphaf*a2
    rl(1,mb1) = 0.d0
    rl(1,mb2) = 0.5d0*th1*alphas*beta2*cf*(cs*cs - (g2/g1)*a2)*dsqrt(rho)
    rl(1,8)   = 0.5d0*th1*alphas*beta3*cf*(cs*cs - (g2/g1)*a2)*dsqrt(rho)

    ! 2 - left eigenvector
    rl(2,1)   = -0.5d0*(beta3*u2-beta2*u3)*beta1
    rl(2,mu1) = 0.d0
    rl(2,mu2) = 0.5d0*beta3*beta1
    rl(2,4)   = -0.5d0*beta2*beta1
    rl(2,5)   = 0.d0
    rl(2,mb1) = 0.d0
    rl(2,mb2) = 0.5d0*beta3*dsqrt(rho)
    rl(2,8)   = -0.5d0*beta2*dsqrt(rho)

    ! 3 - left eigenvector
    rl(3,1)   = th1*(alphas/4.d0)*cf*cf*um2 + 0.5d0*th2*(alphas*ca*u1*beta1 + alphaf*cf*(beta2*u2 + beta3*u3))
    rl(3,mu1) = -0.5d0*th1*alphas*cf*cf*u1 - 0.5d0*th2*alphas*ca*beta1
    rl(3,mu2) = -0.5d0*th1*alphas*cf*cf*u2 - 0.5d0*th2*alphaf*cf*beta2
    rl(3,4)   = -0.5d0*th1*alphas*cf*cf*u3 - 0.5d0*th2*alphaf*cf*beta3
    rl(3,5)   = 0.5d0*th1*alphas*cf*cf
    rl(3,mb1) = 0.d0
    rl(3,mb2) = -0.5d0*th1*alphaf*beta2*cf*(cf*cf - (g2/g1)*a2)*dsqrt(rho)
    rl(3,8)   = -0.5d0*th1*alphaf*beta3*cf*(cf*cf - (g2/g1)*a2)*dsqrt(rho)

    ! 4 - left eigenvector
    rl(4,1)   = 1.d0 - 0.5d0*th1*(alphaf*alphaf*a2 + alphas*alphas*cf*cf)*um2
    rl(4,mu1) = th1*(alphaf*alphaf*a2 + alphas*alphas*cf*cf)*u1
    rl(4,mu2) = th1*(alphaf*alphaf*a2 + alphas*alphas*cf*cf)*u2
    rl(4,4)   = th1*(alphaf*alphaf*a2 + alphas*alphas*cf*cf)*u3
    rl(4,5)   = -th1*(alphaf*alphaf*a2 + alphas*alphas*cf*cf)
    rl(4,mb1) = 0.d0
    rl(4,mb2) = th1*alphaf*alphas*beta2*cf*(cf*cf-cs*cs)*dsqrt(rho)
    rl(4,8)   = th1*alphaf*alphas*beta3*cf*(cf*cf-cs*cs)*dsqrt(rho)

    ! 5 - left eigenvector
    rl(5,1)   = th1*(alphas/4.d0)*cf*cf*um2 - 0.5d0*th2*(alphas*ca*u1*beta1 + alphaf*cf*(beta2*u2 + beta3*u3))
    rl(5,mu1) = -0.5d0*th1*alphas*cf*cf*u1 + 0.5d0*th2*alphas*ca*beta1
    rl(5,mu2) = -0.5d0*th1*alphas*cf*cf*u2 + 0.5d0*th2*alphaf*cf*beta2
    rl(5,4)   = -0.5d0*th1*alphas*cf*cf*u3 + 0.5d0*th2*alphaf*cf*beta3
    rl(5,5)   = 0.5d0*th1*alphas*cf*cf
    rl(5,mb1) = 0.d0
    rl(5,mb2) = -0.5d0*th1*alphaf*beta2*cf*(cf*cf - (g2/g1)*a2)*dsqrt(rho)
    rl(5,8)   = -0.5d0*th1*alphaf*beta3*cf*(cf*cf - (g2/g1)*a2)*dsqrt(rho)

    ! 6 - left eigenvector
    rl(6,1)   = 0.5d0*(beta3*u2-beta2*u3)*beta1
    rl(6,mu1) = 0.d0
    rl(6,mu2) = -0.5d0*beta3*beta1
    rl(6,4)   = 0.5d0*beta2*beta1
    rl(6,5)   = 0.d0
    rl(6,mb1) = 0.d0
    rl(6,mb2) = 0.5d0*beta3*dsqrt(rho)
    rl(6,8)   = -0.5d0*beta2*dsqrt(rho)

    ! 7 - left eigenvector
    rl(7,1)   = th1*(alphaf/4.d0)*a2*um2 - 0.5d0*th2*(alphaf*a*u1*beta1 - alphas*cs*(beta2*u2 + beta3*u3))
    rl(7,mu1) = -0.5d0*th1*alphaf*a2*u1 + 0.5d0*th2*alphaf*a*beta1
    rl(7,mu2) = -0.5d0*th1*alphaf*a2*u2 - 0.5d0*th2*alphas*cs*beta2
    rl(7,4)   = -0.5d0*th1*alphaf*a2*u3 - 0.5d0*th2*alphas*cs*beta3
    rl(7,5)   = 0.5d0*th1*alphaf*a2
    rl(7,mb1) = 0.d0
    rl(7,mb2) = 0.5d0*th1*alphas*beta2*cf*(cs*cs - (g2/g1)*a2)*dsqrt(rho)
    rl(7,8)   = 0.5d0*th1*alphas*beta3*cf*(cs*cs - (g2/g1)*a2)*dsqrt(rho)

    return
end

!=========================================================
subroutine set_rght_eig(mu1, mu2, mb1, mb2, rho, u1, u2, u3, p, B1, B2, B3, gamma, rr)
!=========================================================

    implicit none
    integer, intent(in) :: mu1, mu2, mb1, mb2
    double precision, intent(in) :: rho, u1, u2, u3, p, B1, B2, B3, gamma
    double precision, dimension(8,7), intent(out) :: rr
    double precision a2, a, d, ca, cf, cs, beta1, beta2, beta3, alphaf, alphas, g1, g2, hmf, hpf, hms, hps

    a2 = gamma * p / rho
    a = dsqrt(a2)
    d = a2 + (B1**2 + B2**2 + B3**2) / rho

    ca = dsqrt(B1**2 / rho)
    cf = dsqrt(0.5d0 * ( d + dsqrt(d**2 - 4.d0*a2*B1*B1/rho)))
    cs = dsqrt(0.5d0 * ( d - dsqrt(d**2 - 4.d0*a2*B1*B1/rho)))

    beta1 = dsign(1.d0, B1)
    if ( (dabs(B2) + dabs(B3)) <= 1.d-15 ) then
        beta2 = 1.d0 / dsqrt(2.d0)
        beta3 = 1.d0 / dsqrt(2.d0)
    else
        beta2 = B2 / dsqrt(B2*B2 + B3*B3)
        beta3 = B3 / dsqrt(B2*B2 + B3*B3)
    endif

    if ( (dabs(B2) + dabs(B3)) <= 1.d-15 .and. dabs(a2 - B1*B1/rho) < 1.d-15 ) then
        alphaf = 1.d0
        alphas = 1.d0
    else
        alphaf = dsqrt(dabs(cf*cf - ca*ca)) / dsqrt(dabs(cf*cf-cs*cs))
        alphas = dsqrt(dabs(cf*cf - a2)) / dsqrt(dabs(cf*cf-cs*cs))
    endif

    g1 = gamma - 1.d0
    g2 = gamma - 2.d0
    hmf = (alphaf*cf*cf/g1) - alphaf*cf*u1 + alphas*ca*beta1*(beta2*u2 + beta3*u3) + (g2/g1)*alphaf*(cf*cf-a2)
    hpf = (alphaf*cf*cf/g1) + alphaf*cf*u1 - alphas*ca*beta1*(beta2*u2 + beta3*u3) + (g2/g1)*alphaf*(cf*cf-a2)
    hms = (alphas*cs*cs/g1) - alphas*cs*u1 - alphaf*a*beta1*(beta2*u2 + beta3*u3) + (g2/g1)*alphas*(cs*cs-a2)
    hps = (alphas*cs*cs/g1) + alphas*cs*u1 + alphaf*a*beta1*(beta2*u2 + beta3*u3) + (g2/g1)*alphas*(cs*cs-a2)

    ! 1 - right eigenvector
    rr(1,1)   = alphaf
    rr(mu1,1) = alphaf*(u1-cf)
    rr(mu2,1) = alphaf*u2 + alphas*beta2*ca*beta1
    rr(4,1)   = alphaf*u3 + alphas*beta3*ca*beta1
    rr(5,1)   = 0.5d0*alphaf*(u1**2 + u2**2 + u3**2) + hmf
    rr(mb1,1) = 0.d0
    rr(mb2,1) = alphas*beta2*cf/dsqrt(rho)
    rr(8,1)   = alphas*beta3*cf/dsqrt(rho)

    ! 2 - right eigenvector
    rr(1,2)   = 0.d0
    rr(mu1,2) = 0.d0
    rr(mu2,2) = beta3*beta1
    rr(4,2)   = -beta2*beta1
    rr(5,2)   = (beta3*u2 - beta2*u3)*beta1
    rr(mb1,2) = 0.d0
    rr(mb2,2) = beta3/dsqrt(rho)
    rr(8,2)   = -beta2/dsqrt(rho)

    ! 3 - right eigenvector
    rr(1,3)   = alphas
    rr(mu1,3) = alphas*(u1-cs)
    rr(mu2,3) = alphas*u2 - alphaf*beta2*a*beta1
    rr(4,3)   = alphas*u3 - alphaf*beta3*a*beta1
    rr(5,3)   = 0.5d0*alphas*(u1**2 + u2**2 + u3**2) + hms
    rr(mb1,3) = 0.d0
    rr(mb2,3) = -alphaf*beta2*a2/(cf*dsqrt(rho))
    rr(8,3)   = -alphaf*beta3*a2/(cf*dsqrt(rho))

    ! 4 - right eigenvector
    rr(1,4)   = 1.d0
    rr(mu1,4) = u1
    rr(mu2,4) = u2
    rr(4,4)   = u3
    rr(5,4)   = 0.5d0*(u1**2 + u2**2 + u3**2)
    rr(mb1,4) = 0.d0
    rr(mb2,4) = 0.d0
    rr(8,4)   = 0.d0

    ! 5 - right eigenvector
    rr(1,5)   = alphas
    rr(mu1,5) = alphas*(u1+cs)
    rr(mu2,5) = alphas*u2 + alphaf*beta2*a*beta1
    rr(4,5)   = alphas*u3 + alphaf*beta3*a*beta1
    rr(5,5)   = 0.5d0*alphas*(u1**2 + u2**2 + u3**2) + hps
    rr(mb1,5) = 0.d0
    rr(mb2,5) = -alphaf*beta2*a2/(cf*dsqrt(rho))
    rr(8,5)   = -alphaf*beta3*a2/(cf*dsqrt(rho))

    ! 6 - right eigenvector
    rr(1,6)   = 0.d0
    rr(mu1,6) = 0.d0
    rr(mu2,6) = -beta3*beta1
    rr(4,6)   = beta2*beta1
    rr(5,6)   = -(beta3*u2 - beta2*u3)*beta1
    rr(mb1,6) = 0.d0
    rr(mb2,6) = beta3/dsqrt(rho)
    rr(8,6)   = -beta2/dsqrt(rho)

    ! 7 - right eigenvector
    rr(1,7)   = alphaf
    rr(mu1,7) = alphaf*(u1+cf)
    rr(mu2,7) = alphaf*u2 - alphas*beta2*ca*beta1
    rr(4,7)   = alphaf*u3 - alphas*beta3*ca*beta1
    rr(5,7)   = 0.5d0*alphaf*(u1**2 + u2**2 + u3**2) + hpf
    rr(mb1,7) = 0.d0
    rr(mb2,7) = alphas*beta2*cf/dsqrt(rho)
    rr(8,7)   = alphas*beta3*cf/dsqrt(rho)

    return
end

!=========================================================
subroutine set_wave_spd(rho, u1, u2, u3, p, B1, B2, B3, gamma, s)
!=========================================================

    implicit none
    double precision, intent(in) :: rho, u1, u2, u3, p, B1, B2, B3, gamma
    double precision, dimension(7), intent(out) :: s
    double precision  a2, d, ca, cf, cs

    a2 = gamma * p / rho
    d = a2 + (B1**2 + B2**2 + B3**2) / rho

    ca = dsqrt(B1**2/rho)
    cf = dsqrt(0.5d0 * ( d + dsqrt(d**2 - 4.d0*a2*B1*B1/rho)))
    cs = dsqrt(0.5d0 * ( d - dsqrt(d**2 - 4.d0*a2*B1*B1/rho)))

    s(1) = u1 - cf
    s(2) = u1 - ca
    s(3) = u1 - cs
    s(4) = u1
    s(5) = u1 + cs
    s(6) = u1 + ca
    s(7) = u1 + cf

    return
end

max_stars_repo_name: ranocha/riemann

max_issues_repo_name: ranocha/riemann

max_forks_repo_name: ranocha/riemann


text: #redirect Yolo Community Care Continuum


!--------------------------------------------------------------
! This program is designed for try FORTRAN's Output Format
!--------------------------------------------------------------
! Author: Yang Yang
! Date: 2017/04/05
!--------------------------------------------------------------
!--------------------------------------------------------------

!--------------------------------------------------------------
PROGRAM MAIN
IMPLICIT NONE
! program parameters
INTEGER,PARAMETER:: M = 4, N = 4

! vars in main function
REAL*8 :: A(M,N)
INTEGER i,j

! Define value of array A
DO j = 1,N
  DO i = 1,M
       A(i,j) = (i-1)*N + j + 1.D-2
  END DO
END DO

! Output array to file by calling subroutine
CALL dataOutput(A,M,N)

END PROGRAM MAIN
!--------------------------------------------------------------
!--------------------------------------------------------------

!--------------------------------------------------------------
SUBROUTINE dataOutput(A,M,N)
IMPLICIT NONE

CHARACTER(len=20),PARAMETER:: outputFileName = 'result.txt'   ! Output file name
! Subroutine's input vars
INTEGER M,N
REAL*8 A(M,N)

! Subroutine's locaL vars
INTEGER i,j,ierror
CHARACTER(len=100) M_Temp
WRITE(M_Temp,*) M

! Open OutputFile
OPEN(unit = 10, file = TRIM(outputFileName), IOSTAT = ierror)
IF (ierror.EQ.0) THEN
  PRINT*, 'File result.txt successfully opened!...'
END IF

! Output data to file
DO i = 1,M
  DO j = 1,N
        WRITE(10,100) A(i,j)
  END DO
  WRITE(10,200)
END DO

100 FORMAT(1x,ES24.15,\)          ! Write data to file without enter
200 FORMAT(/)                     ! Wirte enter to file
! Close Output File
CLOSE(10)

END SUBROUTINE dataOutput
!--------------------------------------------------------------
!--------------------------------------------------------------


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 TO SET,CALCULATE OR CALIBRATE NATURAL MORTALITY    *
!  *                       RATE CONSTANT                               *
!  *********************************************************************
!
!  971217 MvdV extra mortality above critical temp for ulvae
!  0895 MvdV output of grazer biomass extended for multiple grazer types

      SUBROUTINE NATMOR(DEATH,ZOOD,TEMP,LCOUPL)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'size.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'cal1.inc'
      INCLUDE 'graas.inc'
      INCLUDE 'ioblck.inc'
      REAL*8 ZOOD(0:MG)
!
!  IF LCAL=1 DEATH IS PUT EQUAL TO THE INPUT VALUE
!
!  IF LCAL=4 DEATH IS CALCULATED AS RMORT1(I) * RMORT2(I) ** TEMP
!
!  DETERMINE HOW DEATH IS TO BE ESTABLISHED
!
      IF (LCAL .EQ. 4) GO TO 70
!
!  LCAL=1
!
      DO 10 I = 1,NUSPEC
   10 RMORT(I) = DEATH
      IF (NREP .GT. 1) GO TO 110
      IF (FLUSH .LT. 1.0D-6) GO TO 30
      IF (IDUMP .EQ. 1) WRITE(IOU(6),20) FLUSH
   20 FORMAT(2X,'Input mortality rates +',2X,F6.2,2X,'are used')
      GO TO 50
   30 CONTINUE
      IF (IDUMP .EQ. 1) WRITE(IOU(6),40)
   40 FORMAT(2X,'Input natural mortality rates are used')
   50 CONTINUE
      GO TO 110
!
!  LCAL=4. STORE MAXIMUM DEATH RATE IN DEATH.
!  Compute mortatlity rate constant or set to BASMOR if option TEMPLIM
!  is on and TEMP is below TEMLIM.
!
   70 CONTINUE
      TEMP2 = TEMP
      IF (LTLIM .EQ. 0) GO TO 90
      IF (TEMP .GE. TEMLIM) GO TO 90
      IF ((NREP .EQ. 1).AND.(LCOUPL.EQ.0)) THEN
         TEMP2 = TEMP + 5.0
         GO TO 90
      END IF
         DO 80 I = 1,NUSPEC
         RMORT(I) = BASMOR
   80    CONTINUE
         DEATH = BASMOR
         GO TO 110
   90 DEATH = 0.0
      DO 100 I = 1,NUSPEC
!       MvdV 971217 addition for ulvae
!       critical temperature in RMORT2
!       RMORT2 (temperature correction) = 1.0
!       temperature muliplication factor in RMORT3
        IF (RMORT2(I).GE.0) THEN
          TMPCOR = RMORT2(I)
        ELSE
          TMPCOR = 1.0
        ENDIF
        RMORT(I) = RMORT1(I) * TMPCOR ** TEMP2
        IF ((RMORT2(I).LT.0.).AND.(TEMP.GT.-1.*RMORT2(I)))
     1    RMORT(I) = MAX(RMORT(I),(TEMP+RMORT2(I)) * RMORT3(I))
        IF (RMORT(I) .GT. DEATH) DEATH = RMORT(I)
  100 CONTINUE
  110 CONTINUE
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


module stpspdmod

!$$$ module documentation block
!           .      .    .                                       .
! module:   stpspdmod    module for stpspd and its tangent linear stpspd_tl
!  prgmmr:
!
! abstract: module for stpspd and its tangent linear stpspd_tl
!
! program history log:
!   2005-05-20  Yanqiu zhu - wrap stpspd and its tangent linear stpspd_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-12-02  Todling - remove stpspd_tl
!   2009-08-12  lueken - update documentation
!   2010-05-13  todling - uniform interface across stp routines
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!
! subroutine included:
!   sub stpspd
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC stpspd

contains

subroutine stpspd(spdhead,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stpspd  calculate penalty and stepsize terms
!                for wind speed, with nonlinear qc.
!   prgmmr: derber           org: np23                date: 1991-02-26
!
! abstract: calculate penalty and stepsize terms for wind speed
!
! program history log:
!   1991-02-26  derber
!   1998-02-03  derber
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-07-29  treadon - add only to module use, add intent in/out
!   2004-10-08  parrish - add nonlinear qc option
!   2005-04-11  treadon - merge stpspd and stpspd_qc into single routine
!   2005-08-02  derber  - modify for variational qc parameters for each ob
!   2005-09-28  derber  - consolidate location and weight arrays
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!                       - unify NL qc
!   2006-09-18  derber  - modify output b1 and b3
!   2007-03-19  tremolet - binning of observations
!   2007-05-10  tremolet - add opt to run as linear procedure
!   2007-06-04  derber  - use quad precision to get reproducability over number of processors
!   2008-12-03  todling - changed handling of ptr%time
!   2009-01-19  todling - re-implement Tremolet's linearization for q1fy10
!   2010-01-04  zhang,b - bug fix: accumulate penalty for multiple obs bins
!   2010-05-13  todling  - update to use gsi_bundle
!   2010-09-25  todling  - fix linearlization
!
!   input argument list:
!     spdhead
!     ru       - search direction for u
!     rv       - search direction for v
!     su       - analysis increment for u
!     sv       - analysis increment for v
!     sges     - step size estimates (nstep)
!     nstep    - number of stepsizes  (==0 means use outer iteration values)
!
!   output argument list
!     out(1:nstep)   - contribution to penalty from wind speed sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_quad
  use qcmod, only: nlnqc_iter,varqc_iter
  use constants, only: zero,half,one,two,tiny_r_kind,cg_term,zero_quad,r3600
  use gsi_4dvar, only: ltlint
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use m_obsNode, only: obsNode
  use m_spdNode, only: spdNode
  use m_spdNode, only: spdNode_typecast
  use m_spdNode, only: spdNode_nextcast
  implicit none

! Declare passed variables
  class(obsNode), pointer             ,intent(in   ) :: spdhead
  integer(i_kind)                     ,intent(in   ) :: nstep
  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges
  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
  type(gsi_bundle)                    ,intent(in   ) :: rval,sval

! Declare local variables
  integer(i_kind) j1,j2,j3,j4,kk,ier,istatus
  real(r_kind) w1,w2,w3,w4,time_spd
  real(r_kind) valu,valv,ucur,vcur,spdnl,spdtl,uu,vv,spd
  real(r_kind),dimension(max(1,nstep)):: pen
  real(r_kind) cg_spd,pencur,wgross,wnotgross
  real(r_kind) pg_spd,pentl
  real(r_kind),pointer,dimension(:) :: su,sv
  real(r_kind),pointer,dimension(:) :: ru,rv
  type(spdNode), pointer :: spdptr

  out=zero_quad

!  If no spd data return
  if(.not. associated(spdhead))return

  time_spd=zero
! Retrieve pointers
! Simply return if any pointer not found
  ier=0
  call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
  if(ier/=0)return

  spdptr => spdNode_typecast(spdhead)
  do while (associated(spdptr))

     if(spdptr%luse)then
        if(nstep > 0)then
           j1 = spdptr%ij(1)
           j2 = spdptr%ij(2)
           j3 = spdptr%ij(3)
           j4 = spdptr%ij(4)
           w1 = spdptr%wij(1)
           w2 = spdptr%wij(2)
           w3 = spdptr%wij(3)
           w4 = spdptr%wij(4)

           valu=w1* ru(j1)+w2* ru(j2)+w3* ru(j3)+w4* ru(j4)
           valv=w1* rv(j1)+w2* rv(j2)+w3* rv(j3)+w4* rv(j4)
           ucur=w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4)+spdptr%uges
           vcur=w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4)+spdptr%vges

           if (ltlint) then
              spd=sqrt(ucur*ucur+vcur*vcur)-spdptr%res
              pencur=spd*spd*spdptr%err2
              do kk=1,nstep
                 spdnl=sqrt(ucur*ucur+vcur*vcur)
                 spdtl=ucur*valu+vcur*valv
                 if (spdnl>tiny_r_kind*100._r_kind) then
                    spdtl=spdtl/spdnl
                 else
                    spdtl=zero
                 endif
                 pentl  =two*spdtl*spd*spdptr%err2
                 pen(kk)=pencur+sges(kk)*pentl
              end do
           else
              do kk=1,nstep
                 uu=ucur+sges(kk)*valu
                 vv=vcur+sges(kk)*valv
                 spd=sqrt(uu*uu+vv*vv)-spdptr%res
                 pen(kk)=spd*spd*spdptr%err2
              end do
           end if
        else
           pen(1)=spdptr%res*spdptr%res*spdptr%err2
        end if

!  Modify penalty term if nonlinear QC
        if (nlnqc_iter .and. spdptr%pg > tiny_r_kind .and. &
                             spdptr%b  > tiny_r_kind) then
           pg_spd=spdptr%pg*varqc_iter
           cg_spd=cg_term/spdptr%b
           wnotgross= one-pg_spd
           wgross = pg_spd*cg_spd/wnotgross
           pencur = -two*log((exp(-half*pencur) + wgross)/(one+wgross))
           do kk=1,max(1,nstep)
              pen(kk) = -two*log((exp(-half*pen(kk)  ) + wgross)/(one+wgross))
           enddo
        endif

        out(1) = out(1)+pen(1)*spdptr%raterr2
        do kk=2,nstep
           out(kk) = out(kk)+(pen(kk)-pen(1))*spdptr%raterr2
        end do

     end if

     spdptr => spdNode_nextcast(spdptr)

  end do
  return
end subroutine stpspd

end module stpspdmod


text:
****************************************************
*
* This is a program to do the triplet mapping for
* droplet inside the linear eddy domain.
*
*******************************************************

	subroutine drop_map(n,i1,m1,ncell
     1  ,flag,index,ibound1,ibound2)

	dimension index(ibound1:ibound2)
	integer j1,j2

c.
c. if this is not the 1st time enter this
c. subroutine, then it is not necessary to
c. set the I.C.s and B.C.s again, go to 1110
c.

	if (flag.eq.1.)goto 1110
c.
c. Do triplet mapping of droplet
c.
c. set I.C.s and B.C.s
c.
	do j=1,ncell
	index(j)=j
	index(j+ncell)=j
Cindex(j-ncell)=j
	end do

1110	continue
c.
c. Be careful not confused here, N is amount of
c. eddy grid points, I1-M1+1 is the 1st point
c. of triplet mapping (because that index(j)=j,
c. so I1 in LEM is equal I1-M1+1 in DROP_MAP.f
c. here. Similarily NN (last point of mapping)
c. is I1-M1+N. ALso think about if there is
c. difference of X(i) in LEM and DROP)
c.
	CALL TRIPLET1(N,I1-M1+1,index)
c write(*,*)N,I1-M1+1
c do j=1,ncell
c write(*,*)j,index(j)
c enddo

	NN=I1-M1+N

c
c check if index(j) is negative
c
	do j=I1-M1+1,I1-M1+N

        if (index(j).le.0)then
	write(*,*)n,I1
        write(*,*)'line 128 m',j,'index(m)',index(j)
        end if

	end do
c.
c. set the 1st point (j1)
c.
	j1=I1-M1+1
	j2=NN
Cwrite(95,*)'jk=',jk,N
Cdo i=j1,j2
Cwrite(95,*)i,index(i)
Cenddo

C
C return to index to the initial value
C

Cdo mp=j1,j2
Cj=mp
Cindex(j-ncell)=index(j)
Cend do

	if(j2.gt.ncell)then
  	do mp=j1,j2
	j=mp
     	index(j-ncell)=index(j)
        enddo
	else
         do mp=j1,j2
	 j=mp
     	 index(j+ncell)=index(j)
 	 enddo
	endif

        return

        end

*********************************************************
*
*       This is a subroutine to perform triplet mapping
*
*********************************************************

        SUBROUTINE TRIPLET1(KK,K1,A)
	parameter(ibound1=0,ibound2=50000)
        Dimension A(ibound1:ibound2),E(ibound1:ibound2)

c
c       kk is the grid number for the eddy, due to the triplet
c       mapping, M is the 1/3 eddy grid pts
c       the logic of triplet mapping can be described as following
c       example
c       initial scalar: 1,2,3,4,5,6,7,8,9
c       after mapping:/  1,4,7,8,5,2,3,6,9
c
        M=kk/3
        MA=(3*M-1)
        MB=3*M-2
        MC=3*M-3

c**************************************************
c
c       STEP1
c
c**************************************************

        DO K=K1,K1+M-1
        E(K)=A(3*M-MA+K1-1)
        MA=MA-3
        END DO

c********************************************************
c
c       STEP2 and STEP2.1 2nd segment and Block Inversion
c
c*********************************************************
        DO K=M+K1, 2*M+K1-1
        E(K)=A(3*M-MB+K1-1)
        MB=MB-3
        END DO
c
c       block inversion of the 2nd segment
c
        MB1=INT((M-1)/2)+(M+K1)
        MW=2*M+K1-1
        DO K=M+K1,MB1
        B1=E(K)
        E(K)=E(MW)

        E(MW)=B1
        MW=MW-1
        END DO

c**********************************************************
c
c       step3 3rd segment
c
c**********************************************************
        DO K=2*M+K1,3*M+K1-1
        E(K)=A(3*M-MC+K1-1)
        MC=MC-3
        END DO

        DO K=K1,KK+K1-1
        A(K)=E(K)
        END DO

        RETURN
        END

max_stars_repo_name: phaustin/empm

max_issues_repo_name: phaustin/empm

max_forks_repo_name: phaustin/empm