      subroutine cli_pgen(iwgn)

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine generates precipitation data when the user chooses to
!!    simulate or when data is missing for particular days in the weather file

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    j           |none          |HRU number
!!    pcp_stat(:,1,:)|mm/day     |average amount of precipitation falling in
!!                               |one day for the month
!!    pcp_stat(:,2,:)|mm/day     |standard deviation for the average daily
!!                               |precipitation
!!    pcp_stat(:,3,:)|none       |skew coefficient for the average daily
!!                               |precipitation
!!    pr_w(1,:,:) |none          |probability of wet day after dry day in month
!!    pr_w(2,:,:) |none          |probability of wet day after wet day in month
!!    rnd3(:)     |none          |random number between 0.0 and 1.0
!!    rndseed(:,:)|none          |random number seeds
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    rnd3(:)     |none          |random number between 0.0 and 1.0
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Log
!!    SWAT: Aunif, Dstn1

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use basin_module
      use climate_module
      use hydrograph_module
      use time_module

      implicit none

      real :: vv               !none          |random number between 0.0 and 1.0
      real :: pcpgen           !mm H2O        |generated precipitation value for the day
      real :: v8               !none          |random number between 0.0 and 1.0
      real :: r6               !none          |variable to hold intermediate calculation
      real :: xlv              !none          |variable to hold intermediate calculation
      real :: aunif            !              |
      real :: xx               !              |
      real :: cli_dstn1        !              |
      integer :: iwgn          !              |

      pcpgen = 0.
      vv = Aunif(rndseed(idg(1),iwgn))
      if (wst(iwst)%weat%precip_prior_day == "dry")  then
        xx = wgn(iwgn)%pr_wd(time%mo)
      else
        xx = wgn(iwgn)%pr_ww(time%mo)
      endif
      if (vv > xx) then
        pcpgen = 0.
      else
        v8 = Aunif(rndseed(idg(3),iwgn))
        !!skewed rainfall distribution
        r6 = wgn(iwgn)%pcpskw(time%mo) / 6.
        xlv = (cli_Dstn1(rnd3(iwgn),v8) - r6) * r6 + 1.
        xlv = (xlv**3 - 1.) * 2. / wgn(iwgn)%pcpskw(time%mo)
        rnd3(iwgn) = v8
        pcpgen = xlv * wgn(iwgn)%pcpstd(time%mo) + wgn_pms(iwgn)%pcpmean(time%mo)
        pcpgen = pcpgen * wgn_pms(iwgn)%pcf(time%mo)
        if (pcpgen < .1) pcpgen = .1
      end if

      !! precip for the next day
      wst(iwst)%weat%precip_next = pcpgen

      return
      end subroutine cli_pgen


text:
      SUBROUTINE XERBLA( SRNAME, INFO )
      CHARACTER*(*)      SRNAME
      INTEGER            INFO

      END

ext: f

max_stars_repo_name: KavrakiLab/amino

max_issues_repo_name: KavrakiLab/amino

max_forks_repo_name: KavrakiLab/amino


text:
Flip Cup is a skill based game where 2 teams flip their cups after drinking its contents in a relay race fashion.

What is needed to play:
 2 teams of equal size
 A long, clean table
 1 cup per person
 Beer (or other alcohol such as wine or liquor NOTE: using other alcohol can end very dangerously, or very funnily)

Playing
 Each person should fill up their cup to about quarter full
 Match Up: Raise your cup the the person across the table from you to ensure that each team has the same number of players and to make sure each opposing player has the same amount of beer.
 To start, the 2 first people must cheers, touch their cups to the table, and drink.
 When it is their turn, each player must drink the contents of their cup, place it upside down on the edge of the table, and flip it 360 degrees until it lands upside down on the table again and comes to a rest. Each player must do this until all players have done it successfully. Note: many people play Tip Cup where you start with the cup right side up and flip it to upside down. This requires much less skill.
 The next person to go may not touch their cup until the person before has successfully flipped
 The game ends when the last person on either team finishes

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE XTRNLD(I1,I2,KCASE)
*
*
*       External force & derivatives.
*       -----------------------------
*
      INCLUDE 'common6.h'
      REAL*8  XI(3),XIDOT(3),FIRR(3),FREG(3),FD(3),FDR(3)
*
*
*       See whether to include the standard or point-mass tidal force.
      IF (KZ(14).LE.2.AND.KCASE.EQ.1) THEN
*       Include tidal force & first derivative (I1 = I2 for single body).
          DO 10 I = I1,I2
              FI(1,I) = FI(1,I) + TIDAL(4)*XDOT(2,I)
              FI(2,I) = FI(2,I) - TIDAL(4)*XDOT(1,I)
              FR(1,I) = FR(1,I) + TIDAL(1)*X(1,I)
              FR(3,I) = FR(3,I) + TIDAL(3)*X(3,I)
              D1(1,I) = D1(1,I) + TIDAL(4)*(FI(2,I) + FR(2,I))
              D1(2,I) = D1(2,I) - TIDAL(4)*(FI(1,I) + FR(1,I))
              D1R(1,I) = D1R(1,I) + TIDAL(1)*XDOT(1,I)
              D1R(3,I) = D1R(3,I) + TIDAL(3)*XDOT(3,I)
   10     CONTINUE
      END IF
*
      IF (KZ(14).LE.2.AND.KCASE.EQ.2) THEN
*       Add the second and third derivatives due to the tidal field.
          DO 20 I = I1,I2
              D2(1,I) = D2(1,I) + TIDAL(4)*FDOT(2,I)
              D2(2,I) = D2(2,I) - TIDAL(4)*FDOT(1,I)
              D2R(1,I) = D2R(1,I) + TIDAL(1)*F(1,I)
              D2R(3,I) = D2R(3,I) + TIDAL(3)*F(3,I)
              D3(1,I) = D3(1,I) + TIDAL(4)*(D2(2,I) + D2R(2,I))
              D3(2,I) = D3(2,I) - TIDAL(4)*(D2(1,I) + D2R(1,I))
              D3R(1,I) = D3R(1,I) + TIDAL(1)*FDOT(1,I)
              D3R(3,I) = D3R(3,I) + TIDAL(3)*FDOT(3,I)
   20     CONTINUE
      END IF
*
*       Check options for galactic potential (called from FPOLY1).
      IF (KZ(14).GE.3.AND.KCASE.EQ.1) THEN
*       Obtain regular force components for 3D galaxy or Plummer model.
          DO 30 I = I1,I2
              DO 25 K = 1,3
                  XI(K) = X(K,I)
                  XIDOT(K) = XDOT(K,I)
                  FIRR(K) = 0.0
                  FD(K) = 0.0
                  FREG(K) = 0.0
                  FDR(K) = 0.0
   25         CONTINUE
*      Jpetts, added argument I to end of XTRNLF
              CALL XTRNLF(XI,XIDOT,FIRR,FREG,FD,FDR,1,I)
              DO 28 K = 1,3
                  FR(K,I) = FR(K,I) + FREG(K)
                  D1R(K,I) = D1R(K,I) + FDR(K)
   28         CONTINUE
   30     CONTINUE
*
      END IF
*
      RETURN
*
      END


text:
Module m_indmed
Integer, Parameter :: kdp = selected_real_kind(15)
public :: indmed
private :: kdp
private :: R_indmed, I_indmed, D_indmed
private :: r_med, i_med, d_med
Integer, Allocatable, Dimension(:), Private, Save :: IDONT
interface indmed
  module procedure d_indmed, r_indmed, i_indmed
end interface indmed
contains

Subroutine D_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call d_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine D_indmed
   Recursive Subroutine d_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out):: ires_med
! __________________________________________________________
!
      Real (kind=kdp), Parameter :: XHUGE = HUGE (XDATT)
      Real (kind=kdp) :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call d_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine d_med
!
Subroutine R_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Real, Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call r_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine R_indmed
   Recursive Subroutine r_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Real, Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out) :: ires_med
! __________________________________________________________
!
      Real, Parameter :: XHUGE = HUGE (XDATT)
      Real :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call r_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine r_med
Subroutine I_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Integer, Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call i_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine I_indmed
   Recursive Subroutine i_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Integer, Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out) :: ires_med
! __________________________________________________________
!
      Integer, Parameter :: XHUGE = HUGE (XDATT)
      Integer :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call i_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine i_med
end module m_indmed

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      subroutine mnmx(vx,vy,vz,x,y,a,llx,lly,lxini,lxe,lyin,lye,s)
      dimension x(llx),y(lly),a(llx,lly)
      common/gfvp/ v1,v2,v3,xr,yr,infin
      v1=vx
      v2=vy
      v3=vz
      lxin=lxini
      if(lxin.ge.0)   go to 55
      lxin=-lxin
      k=-1
      go to   56
  55  k=1
56    continue
      xmin=x(lxin)
      ymin=y(lyin)
      xmax=x(lxe)
      ymax=y(lye)
C
         if (v1.le.xmin)   goto 11
      if(v1.ge.xmax)   go to 12
      if(v2.le.ymin)   go to 26
      if (v2.ge.ymax)   go to 27
C
C     this is point v inside surface
      do 1 mv=lxin,lxe
      if(v1.le.x(mv))  go to 3
  1   continue
   3  v1=x(mv)
      if(mv.ge.lxe)  go to 28
      do 2 mg=lyin,lye
      if(v2.le.y(mg))   go to 4
  2   continue
  4   v2 = y (mg)
      if(mg.le.lye)   goto  27
C
C     this is south
26    do 7 m=lxin,lxe
      if(v1.le.x(m))   go to 8
  7   continue
  8   v1=x(m)
      go to   30
C
C     this is   north
27    do 31 m=lxin,lxe
      if(v1.le.x(m))   go to 32
 31    continue
 32   v1=x(m)
      go to 30
11    if(v2.le.ymin)   goto 30
      if(v2.ge.ymax)   goto 30
C     this is west
       do 5 m=lyin,lye
      if(v2.le.y(m))   go to 6
   5  continue
6      v2=y(m)
      go to 30
12     if(v2.le.ymin)  go to 30
      if(v2.ge.ymax)   go to 30
C     this is east
 28    do 33 m=lyin,lye
      if (v2.le.y(m))   go to 34
  33  continue
  34  v2=y(m)
       go to   30
C   --------------------------------------------------------------------
30    jbg=0
       do 99 j=lyin,lye
       do 99 i=lxin,lxe
      call map(x(i),y(j),a(i,j))
             if(infin.gt.0)   go to 99
      jbg=jbg+1
       if(jbg.gt.1)   go to 71
      xmin=xr
      xmax = xr
      ymin = yr
      ymax = yr
      go to  99
71    continue
      if (xmin.gt.xr)  xmin=xr
      if (xmax.lt.xr)   xmax=xr
      if (ymin.gt.yr)  ymin=yr
      if (ymax.lt.yr)  ymax=yr
   99 continue
       if(k.ge.0)   goto 52
             if(xmn.gt.xmin)   xmn=xmin
      if(xmx.lt.xmax)   xmx=xmax
             if(ymn.gt.ymin)   ymn=ymin
             if(ymx.lt.ymax)   ymx=ymax
             go to   53
C  ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
52    continue
       xmn=xmin
       xmx=xmax
       ymn=ymin
       ymx=ymax
53    call limits(xmn,xmx,ymn,ymx)
       s=(xmx-xmn)/(ymx-ymn)
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
      SUBROUTINE BFSMAT (ND,NE,NB,NP,NTP,LENGTH,NTOTAL,SCR1,JF,JL,NAS,
     1                   FMACH,YB,ZB,YS,ZS,X,DELX,EE,XIC,SG,CG,AR,RIA,
     2                   NBEA1,NBEA2,NASB,NSARAY,NCARAY,BFS,AVR,CBAR,
     3                   A0,XIS1,XIS2,KR,NSBEA,NT0)
C
C     NOTE:
C     A JUMP (VIA AN ASSIGN STATEMENT) TO 200 AND A JUMP TO 1100 (ALSO
C     VIA AN ASSIGN STATEMENT), INTO THE MIDDLES OF SOME DO LOOPS, ARE
C     ACCEPTABLE ANSI 77 FORTRAN. HOWEVER, IBM COMPILER MAY COMPLAIN.
C     THIS PROBLEM IS NOW ELIMINATED (BY G.C. 9/89)
C
C
C        ND         SYMMETRY FLAG
C        NE         GROUND EFFECTS FLAG
C        NB         NUMBER OF BODIES
C        NP         NUMBER OF PANELS
C        NTP        NUMBER OF LIFTING SURFACE BOXES
C        NTOTAL     NTP + TOTAL NO. OF Y AND Z ORIENTED BODY ELEMENTS
C        LENGTH     NTOTAL + THE TOTAL NUMBER OF Z- AND Y-ORIENTED
C                   SLENDER BODY ELEMENTS
C        SCR1       FILE FOR OUTPUT
C        JF         ROW FOR FIRST ZY BODY
C        JL         ROW FOR LAST  ZY BODY
C        NAS        ARRAY CONTAINING THE NUMBER OF ASSOCIATED BODIES
C                   FOR EACH PANEL
C        FMACH      MACH NUMBER
C        YB         ARRAY OF -Y- COORDINATES OF THE BODIES
C        ZB         ARRAY OF -Z- COORDINATES OF THE BODIES
C        YS         ARRAY OF -Y- COORDINATES OF STRIPS AND BODIES
C        ZS         ARRAY OF -Z- COORDINATES OF STRIPS AND BODIES
C        X          ARRAY OF 3/4 CHORD LOCATIONS OF BOXES AND
C                            1/2 CHORD FOR BODY ELEMENTS
C        DELX       ARRAY OF LENGTHS OF BOXES AND BODY ELEMENTS
C        EE         ARRAY OF THE SEMI-WITH OF STRIPS
C        XIC        ARRAY OF 1/4 CHORD COORDINATES OF BOXES
C        SG         ARRAY OF SINE   OF STRIP DIHEDRAL ANGLE
C        CG         ARRAY OF COSINE OF STRIP DIHEDRAL ANGLE
C        AR         ARRAY OF RATIO OF MAJOR AXES OF BODIES
C        RIA        ARRAY OF RADII OF BODY ELEMENTS
C        NBEA1      ARRAY OF NUMBER OF BODY ELEMENTS PER BODY
C        NBEA2      ARRAY OF THE BODY ORIENTATION FLAGS PER BODY
C        NASB       ARRAY OF THE BODIES ASSOCIATED WITH PANELS
C        NSARAY     ARRAY OF THE NUMBER OF STRIPS PER PANEL
C        NCARAY     ARRAY OF THE NUMBER OF CHORDWISE DIV. PER PANEL
C        BFS        WORK ARRAY FOR TEMPORARY STORAGE OF THE BFS COLS.
C        AVR        ARRAY OF RADII OF BODIES
C        CBAR       REFERENCE CHORD
C        A0         ARRAY OF SLENDER BODY ELEMENT RADII
C        XIS1       ARRAY OF SLENDER BODY ELEMENT LEADING  EDGE COORD.S
C        XIS2       ARRAY OF SLENDER BODY ELEMENT TRAILING EDGE COORD.S
C        KR         REDUCED FREQUENCY
C        NSBEA      ARRAY OF THE NUMBER OF ELEMENTS PER SLENDER BODY
C
      LOGICAL    LAST
      INTEGER    SCR1
      REAL       KR
      COMPLEX    BFS(LENGTH,2) , FWZ, FWY, EIKJ1 , EIKJ2
      DIMENSION  YB(1), ZB(1), YS(1), ZS(1), X(1), DELX(1), EE(1),
     1           XIC(1), SG(1), CG(1), AR(1), RIA(1), AVR(1), XIS1(1),
     2           XIS2(1), A0(1), NAS(1), NASB(1), NBEA1(1), NBEA2(1) ,
     3           NSBEA(1), NCARAY(1), NSARAY(1)
C
C
      BETA2 = 1.0 - FMACH**2
      ICOL  = 0
      KSP   = 1
      GO TO 3000
C
   50 CONTINUE
C
C     -Y- ORIENTED BODIES AS SENDING ELEMENTS
C
      SGS   =-1.0
      CGS   = 0.0
      NASD  = 0
      IZYFLG= 3
      ASSIGN 4010 TO IBODY
      GO TO 1000
  100 CONTINUE
C
C     - LIFTING SURF. BOXES AS SENDING ELEMENTS
C
      IF (NTP .LE. 0) GO TO 800
      J     = 1
      JP1   = J
      IBOX  = 0
      ISTRIP= 0
      ISN   = 0
      KSP   = 1
C
C     LOOP FOR -PANEL-
C
      DO 700 ISP = 1,NP
      NS  = NSARAY(ISP)
      NC  = NCARAY(ISP)
      NS  = (NS-ISN) / NC
      ISN = NSARAY(ISP)
      NASD= NAS(ISP)
C
C     LOOP FOR -STRIP-
C
      DO 600 IS = 1,NS
      ISTRIP= ISTRIP + 1
      DYS   = YS(ISTRIP)
      DZS   = ZS(ISTRIP)
      SGS   = SG(ISTRIP)
      CGS   = CG(ISTRIP)
      WIDTH = 2.0 * EE(ISTRIP)
C
C     LOOP FOR -BOX-
C
      DO 500 IB = 1,NC
      IBOX = IBOX + 1
      DXS  = XIC(IBOX)
C
      ICOL = ICOL + 1
C
      CALL FWMW (ND,NE,SGS,CGS,IRB,DRIA,AR,DXLE,DXTE,YB,ZB,DXS,DYS,
     1           DZS,NASD,NASB(KSP),KR,BETA2,CBAR,AVR,FWZ,FWY)
      BFS(ICOL,1) = FWZ * (DXTE - DXLE)
      BFS(ICOL,2) = FWY * (DXTE - DXLE)
      BFS(ICOL,1) = BFS(ICOL,1) * SCALE
      BFS(ICOL,2) = BFS(ICOL,2) * SCALE
C
      AREA = WIDTH * DELX(IBOX)
      BFS(ICOL,1) = BFS(ICOL,1) * AREA
      BFS(ICOL,2) = BFS(ICOL,2) * AREA
  500 CONTINUE
  600 CONTINUE
      KSP = KSP + NASD
  700 CONTINUE
C
C     -Z-  ORIENTED BODIES AS SENDING ELEMENTS
C
  800 CONTINUE
      SGS   = 0.0
      CGS   = 1.0
      NASD  = 0
      IZYFLG= 1
      ASSIGN   50 TO IBODY
      GO TO  1000
C
C
C     *** LOOP FOR EACH INTERFERENCE BODY SENDING ELEMENT
C
 1000 CONTINUE
      INDEX = NTP
C
C     --ISB-- IS THE SENDING BODY
C
      DO 1900 ISB  = 1,NB
      IF (NBEA2(ISB) .EQ.     2 ) GO TO 1070
      IF (NBEA2(ISB) .NE. IZYFLG) GO TO 1850
 1070 DYS   = YB(ISB)
      NSBE  = NBEA1(ISB)
      JP1   = 1
      LAST  = .FALSE.
      DZS   = ZB(ISB)
      EARG2 = 1.0
C
C     --ISBE-- IS THE ELEMENT OF THE SEND BODY
C
      DO 1800 ISBE = 1,NSBE
      EARG1 = EARG2
      INDEX = INDEX + 1
      DXS   = X (INDEX) - DELX(INDEX) /4.0
      EARG2 = KR * DELX(INDEX) / CBAR
C
C     CALCULATE THIS COLUMN
C
      ICOL  = ICOL + 1
      EIKJ1 = CMPLX(COS(EARG1),-SIN(EARG1))
      EIKJ2 = CMPLX(COS(EARG2), SIN(EARG2))
C
      CALL FWMW (ND,NE,SGS,CGS,IRB,DRIA,AR,DXLE,DXTE,YB,ZB,DXS,DYS,
     1           DZS,NASD,NASB(KSP),KR,BETA2,CBAR,AVR,FWZ,FWY)
      BFS(ICOL,1) = FWZ * (DXTE - DXLE)
      BFS(ICOL,2) = FWY * (DXTE - DXLE)
      BFS(ICOL,1) = BFS(ICOL,1) * SCALE
      BFS(ICOL,2) = BFS(ICOL,2) * SCALE
C
C
C     IS THIS THE FIRST COLUMN, YES  BRANCH
C
      IF (ISBE .EQ. 1) GO TO 1800
      BFS(ICOL-1,1) = BFS(ICOL-1,1)*EIKJ1 - BFS(ICOL,1)*EIKJ2
      BFS(ICOL-1,2) = BFS(ICOL-1,2)*EIKJ1 - BFS(ICOL,2)*EIKJ2
 1800 CONTINUE
      GO TO 1900
 1850 INDEX = INDEX + NBEA1(ISB)
 1900 CONTINUE
C
C     RETURN TO CALLING POINT - EITHER Y OR Z SENDING BODY ELEM
C
C     *** GO  EITHER TO THE  Y-ORIENTED INTERFERENCE BODY ELEMENT LOOP
C         OR  TO THE LOOP FOR SLENDER BODY SENDING ELEMENTS
C
      GO TO IBODY, (50,4010)
C
C
C     CALCULATE EACH ROW OF THE SENDING COLUMN
C
 3000 CONTINUE
      IY   = 0
      JF   = 0
      NW   = LENGTH*2
      IROW = 0
C
C     --IRB-- IS THE RECEIVING BODY
C
      IRB  = 0
 3050 IRB  = IRB + 1
      IF (IRB .GT. NB) GO TO 3900
      NRBE = NSBEA(IRB)
      ITSB = NBEA2(IRB)
C
      XYB  = YB(IRB)
      XZB  = ZB(IRB)
      SCALE= 1.0
      IF (ND.NE.0 .AND. XYB.EQ.0.0) SCALE = .5
      IF (NE.NE.0 .AND. XZB.EQ.0.0) SCALE = SCALE*.5
C
C     --IRBE-- IS THE ELEM. OF THE REC. BODY
C
      IRBE = 0
 3060 IRBE = IRBE + 1
      IF (IRBE .GT. NRBE) GO TO 3800
      IY   = IY   + 1
      IROW = IROW + 1
      DRIA = A0  (IY)
      DXLE = XIS1(IY)
      DXTE = XIS2(IY)
      XX1  = DXLE
      XX2  = DXTE
      XAA  = DRIA
      ICOL = 0
      GO  TO  100
C
 3100 CONTINUE
      GO TO (3110,3120,3130), ITSB
 3110 CALL WRITE (SCR1,BFS(1,1),NW,0)
      GO TO 3140
 3120 CALL WRITE (SCR1,BFS(1,2),NW,0)
      CALL WRITE (SCR1,BFS(1,1),NW,0)
      IF (JF .EQ. 0) JF = IROW
      IROW = IROW + 1
      GO TO 3140
 3130 CALL WRITE (SCR1,BFS(1,2),NW,0)
      IROW = IROW - 1
 3140 CONTINUE
      GO TO 3060
 3800 CONTINUE
      GO TO 3050
 3900 CONTINUE
      JL = IROW
      RETURN
C
C
 4010 CONTINUE
C
C
C     *** LOOP FOR EACH SLENDER BODY SENDING ELEMENT
C
      IZYFLG= 1
      SGS   = 0.0
      CGS   = 1.0
 4050 CONTINUE
      LSBE  =  0
      DO 5000 LSB = 1,NB
C
C     --LSB-- IS THE INDEX OF THE SLENDER SENDING BODY
C
      IF (NSBEA(LSB) .EQ.      0) GO TO 5000
      IF (NBEA2(LSB) .EQ.      2) GO TO 4070
      IF (NBEA2(LSB) .NE. IZYFLG) GO TO 4097
 4070 CONTINUE
      XETA  =  YB(LSB)
      XZETA =  ZB(LSB)
      SCALE2=  SCALE
      MSBE  =  NSBEA(LSB)
      DO 4080 LSBS = 1,MSBE
      LSBE  =  LSBE + 1
      ICOL  =  ICOL + 1
      XXIJ  = .50 * XIS1(LSBE) + .50 * XIS2(LSBE)
      CALL FWMW (ND,NE,SGS,CGS,IRB,DRIA,AR,DXLE,DXTE,YB,ZB,XXIJ,XETA,
     1           XZETA,NASD,NASB,KR,BETA2,CBAR,AVR,FWZ,FWY)
      BFS(ICOL,1) =  FWZ * (DXTE - DXLE)
      BFS(ICOL,2) =  FWY * (DXTE - DXLE)
      BFS(ICOL,1) =  BFS(ICOL,1) * SCALE2
      BFS(ICOL,2) =  BFS(ICOL,2) * SCALE2
 4080 CONTINUE
      GO TO 5000
 4097 LSBE = LSBE + NSBEA(LSB)
 5000 CONTINUE
      IF (IZYFLG .EQ. 3) GO TO  5010
      IZYFLG = 3
      SGS =-1.0
      CGS = 0.0
      GO TO 4050
 5010 CONTINUE
C
      GO TO 3100
C
      END


text:
!----- LGPL --------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This library is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License as published by the Free Software Foundation version 2.1.
!
!  This library 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
!  Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with this library; 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, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.

!  $Id: ec_instance.f90 5725 2016-01-13 11:51:32Z leander $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_lgpl/ec_module/packages/ec_module/src/ec_instance.f90 $

!> This module contains the constructor and destructor for the datatype tEcInstance.
!! @author edwin.bos@deltares.nl
module m_ec_instance
   use m_ec_typedefs
   use m_ec_parameters
   use m_ec_message
   use m_ec_alloc
   use m_ec_converter
   use m_ec_connection
   use m_ec_item
   use m_ec_quantity
   use m_ec_elementSet
   use m_ec_field
   use m_ec_filereader
   use m_ec_bcreader
   use m_ec_netcdf_timeseries

   implicit none

   private

   public :: ecInstanceCreate
   public :: ecInstanceFree
   public :: ecInstanceCreateItem
   public :: ecInstanceCreateQuantity
   public :: ecInstanceCreateElementSet
   public :: ecInstanceCreateField
   public :: ecInstanceCreateFileReader
   public :: ecInstanceCreateBCBlock
   public :: ecInstanceCreateNetCDF
   public :: ecInstanceCreateConverter
   public :: ecInstanceCreateConnection
   public :: ecInstancePrintState
   public :: ecInstanceListSourceItems
   public :: ecInstanceListFileReaders

   contains

      ! =======================================================================

      !> Dynamically allocates memory for a new tEcInstance, making the tEcInstance an unnamed data object.
      !! Subsequently, the new tEcInstance is initialized.
      function ecInstanceCreate(ptr) result (success)
         logical                    :: success !< function status
         type(tEcInstance), pointer :: ptr     !< intent(out)
         integer                    :: istat   !< allocate() status
         !
         success = .false.
         !
         if (associated(ptr)) then
            call setECMessage("ERROR: ec_instance::ecInstanceCreate: Dummy argument ptr is already associated.")
         else
            allocate(ptr, STAT = istat)
            if (istat == 0) then
               success = .true.
               ! Allocate arrays with a default size of 10.
               allocate(ptr%ecConnectionsPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecConnectionsPtr array.")
                  success = .false.
               end if
               ptr%nConnections = 0
               allocate(ptr%ecConvertersPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecConvertersPtr array.")
                  success = .false.
               end if
               ptr%nConverters = 0
               allocate(ptr%ecElementSetsPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecElementSetsPtr array.")
                  success = .false.
               end if
               ptr%nElementSets = 0
               allocate(ptr%ecFieldsPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecFieldsPtr array.")
                  success = .false.
               end if
               ptr%nFields = 0
               allocate(ptr%ecFileReadersPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecFileReadersPtr array.")
                  success = .false.
               end if
               ptr%nFileReaders = 0
               allocate(ptr%ecBCBlocksPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecBCBlocksPtr array.")
                  success = .false.
               end if
               ptr%nNetCDFs = 0
               allocate(ptr%ecNetCDFsPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecNetCDFsPtr array.")
                  success = .false.
               end if
               ptr%nBCBlocks = 0
               allocate(ptr%ecItemsPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecItemsPtr array.")
                  success = .false.
               end if
               ptr%nItems = 0
               allocate(ptr%ecQuantitiesPtr(10), STAT = istat)
               if (istat /= 0) then
                  call setECMessage("ERROR: ec_instance::ecInstanceCreate: Unable to allocate memory for ecQuantitiesPtr array.")
                  success = .false.
               end if
               ptr%nQuantities = 0
               ptr%idCounter = 0
            end if
         end if
      end function ecInstanceCreate

      ! =======================================================================

      !> Frees and deallocates dynamically allocated memory and thereby nullifies the tEcInstance's pointers.
      function ecInstanceFree(ptr) result (success)
         logical                    :: success !< function status
         type(tEcInstance), pointer :: ptr     !< intent(inout)
         !
         integer :: istat !< deallocate() status
         !
         success = .false.
         !
         if (.not. associated(ptr)) then
            call setECMessage("INFO: ec_instance::ecInstanceFree: Dummy argument ptr is already disassociated.")
         else
            ! Delegate Free-and-deallocate call to all constituent data types.
            if (.not.ecConnectionFree1dArray(ptr%ecConnectionsPtr, ptr%nConnections)) then
               return
            endif
            if (.not.ecConverterFree1dArray(ptr%ecConvertersPtr, ptr%nConverters)) then
               return
            endif
            if (.not.ecElementSetFree1dArray(ptr%ecElementSetsPtr, ptr%nElementSets)) then
               return
            endif
            if (.not.ecFieldFree1dArray(ptr%ecFieldsPtr, ptr%nFields)) then
               return
            endif
            if (.not.ecFileReaderFree1dArray(ptr%ecFileReadersPtr, ptr%nFileReaders)) then
               return
            endif
            if (.not.ecBCBlockFree1dArray(ptr%ecBCBlocksPtr, ptr%nBCBlocks)) then
               return
            endif
            if (.not.ecNetCDFFree1dArray(ptr%ecNetCDFsPtr, ptr%nNetCDFs)) then
               return
            endif
            if (.not.ecItemFree1dArray(ptr%ecItemsPtr, ptr%nItems)) then
               return
            endif
            if (.not.ecQuantityFree1dArray(ptr%ecQuantitiesPtr, ptr%nQuantities)) then
               return
            endif
               ! Finally deallocate the tEcInstance pointer.
            deallocate(ptr, stat = istat)
            if (istat == 0) success = .true.
         end if
      end function ecInstanceFree

      ! =======================================================================
      ! Create methods
      ! =======================================================================

      !> Create a new Item and register it with the Instance.
      !! Failure is indicated by: itemId == ec_undef_int
      function ecInstanceCreateItem(instancePtr) result(itemId)
         integer                    :: itemId      !< id of the new Item
         type(tEcInstance), pointer :: instancePtr !< intent(in)
         !
         type(tEcItem), pointer :: itemPtr !< the new Item
         !
         itemPtr => null()
         itemId = ec_undef_int
         itemPtr => ecItemCreate(instancePtr%idCounter + 1)
         if (associated(itemPtr)) then
            ! ensure capacity
            if (instancePtr%nItems == size(instancePtr%ecItemsPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecItemsPtr, instancePtr%nItems)) then
                  return
               end if
            end if
            ! register the Item
            instancePtr%nItems = instancePtr%nItems + 1
            instancePtr%ecItemsPtr(instancePtr%nItems)%ptr => itemPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            itemId = instancePtr%idCounter
         end if
      end function ecInstanceCreateItem

      ! =======================================================================

      !> Create a new Quantity and register it with the Instance.
      !! Failure is indicated by: quantityId == ec_undef_int
      function ecInstanceCreateQuantity(instancePtr) result(quantityId)
         integer                    :: quantityId  !< id of the new Quantity
         type(tEcInstance), pointer :: instancePtr !< intent(in)
         !
         type(tEcQuantity), pointer :: quantityPtr !< the new Quantity
         !
         quantityPtr => null()
         quantityId = ec_undef_int
         quantityPtr => ecQuantityCreate(instancePtr%idCounter + 1)
         if (associated(quantityPtr)) then
            ! ensure capacity
            if (instancePtr%nQuantities == size(instancePtr%ecQuantitiesPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecQuantitiesPtr, instancePtr%nQuantities)) then
                  return
               end if
            end if
            ! register the Quantity
            instancePtr%nQuantities = instancePtr%nQuantities + 1
            instancePtr%ecQuantitiesPtr(instancePtr%nQuantities)%ptr => quantityPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            quantityId = instancePtr%idCounter
         end if
      end function ecInstanceCreateQuantity

      ! =======================================================================

      !> Create a new ElementSet and register it with the Instance.
      !! Failure is indicated by: elementSetId == ec_undef_int
      function ecInstanceCreateElementSet(instancePtr) result(elementSetId)
         integer                    :: elementSetId !< id of the new ElementSet
         type(tEcInstance), pointer :: instancePtr  !< intent(in)
         !
         type(tEcElementSet), pointer :: elementSetPtr !< the new ElementSet
         !
         elementSetPtr => null()
         elementSetId = ec_undef_int
         elementSetPtr => ecElementSetCreate(instancePtr%idCounter + 1)
         if (associated(elementSetPtr)) then
            ! ensure capacity
            if (instancePtr%nElementSets == size(instancePtr%ecElementSetsPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecElementSetsPtr, instancePtr%nElementSets)) then
                  return
               end if
            end if
            ! register the ElementSet
            instancePtr%nElementSets = instancePtr%nElementSets + 1
            instancePtr%ecElementSetsPtr(instancePtr%nElementSets)%ptr => elementSetPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            elementSetId = instancePtr%idCounter
         end if
      end function ecInstanceCreateElementSet

      ! =======================================================================

      !> Create a new Field and register it with the Instance.
      !! Failure is indicated by: fieldId == ec_undef_int
      function ecInstanceCreateField(instancePtr) result(fieldId)
         integer                    :: fieldId     !< id of the new Field
         type(tEcInstance), pointer :: instancePtr !< intent(in)
         !
         type(tEcField), pointer :: fieldPtr !< the new Field
         !
         fieldPtr => null()
         fieldId = ec_undef_int
         fieldPtr => ecFieldCreate(instancePtr%idCounter + 1)
         if (associated(fieldPtr)) then
            ! ensure capacity
            if (instancePtr%nFields == size(instancePtr%ecFieldsPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecFieldsPtr, instancePtr%nFields)) then
                  return
               end if
            end if
            ! register the Field
            instancePtr%nFields = instancePtr%nFields + 1
            instancePtr%ecFieldsPtr(instancePtr%nFields)%ptr => fieldPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            fieldId = instancePtr%idCounter
         end if
      end function ecInstanceCreateField

      ! =======================================================================

      !> Create a new Connection and register it with the Instance.
      !! Failure is indicated by: connectionId == ec_undef_int
      function ecInstanceCreateConnection(instancePtr) result(connectionId)
         integer                    :: connectionId !< id of the new Connection
         type(tEcInstance), pointer :: instancePtr  !< intent(in)
         !
         type(tEcConnection), pointer :: connectionPtr !< the new Connection
         !
         connectionPtr => null()
         connectionId = ec_undef_int
         connectionPtr => ecConnectionCreate(instancePtr%idCounter + 1)
         if (associated(connectionPtr)) then
            ! ensure capacity
            if (instancePtr%nConnections == size(instancePtr%ecConnectionsPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecConnectionsPtr, instancePtr%nConnections)) then
                  return
               end if
            end if
            ! register the Connection
            instancePtr%nConnections = instancePtr%nConnections + 1
            instancePtr%ecConnectionsPtr(instancePtr%nConnections)%ptr => connectionPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            connectionId = instancePtr%idCounter
         end if
      end function ecInstanceCreateConnection

      ! =======================================================================

      !> Create a new Converter and register it with the Instance.
      !! Failure is indicated by: converterId == ec_undef_int
      function ecInstanceCreateConverter(instancePtr) result(converterId)
         integer                    :: converterId !< id of the new Converter
         type(tEcInstance), pointer :: instancePtr !< intent(in)
         !
         type(tEcConverter), pointer :: converterPtr !< the new Converter
         !
         converterPtr => null()
         converterId = ec_undef_int
         converterPtr => ecConverterCreate(instancePtr%idCounter + 1)
         if (associated(converterPtr)) then
            ! ensure capacity
            if (instancePtr%nConverters == size(instancePtr%ecConvertersPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecConvertersPtr, instancePtr%nConverters)) then
                  return
               end if
            end if
            ! register the Converter
            instancePtr%nConverters = instancePtr%nConverters + 1
            instancePtr%ecConvertersPtr(instancePtr%nConverters)%ptr => converterPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            converterId = instancePtr%idCounter
         end if
      end function ecInstanceCreateConverter

      ! =======================================================================

      !> Create a new  FileReader and register it with the Instance.
      !! Failure is indicated by: converterId == ec_undef_int
      function ecInstanceCreateFileReader(instancePtr) result(fileReaderId)
         integer                               :: fileReaderId !< id of the new FileReader
         type(tEcInstance), pointer            :: instancePtr  !< intent(in)
         !
         type(tEcFileReader), pointer :: fileReaderPtr !< the new FileReader
         !
         fileReaderId = ec_undef_int
         fileReaderPtr => ecFileReaderCreate(instancePtr%idCounter + 1)
         if (associated(fileReaderPtr)) then
            ! ensure capacity
            if (instancePtr%nFileReaders == size(instancePtr%ecFileReadersPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecFileReadersPtr, instancePtr%nFileReaders)) then
                  return
               end if
            end if
            ! register the FileReader
            instancePtr%nFileReaders = instancePtr%nFileReaders + 1
            instancePtr%ecFileReadersPtr(instancePtr%nFileReaders)%ptr => fileReaderPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            fileReaderId = instancePtr%idCounter
         end if
      end function ecInstanceCreateFileReader

      !> Create a new BCBlock and register it with the Instance (analog of filereader creation).
      function ecInstanceCreateBCBlock(instancePtr) result(bcblockId)
         integer                               :: bcblockId    !< id of the new BCBlock
         type(tEcInstance), pointer            :: instancePtr  !< intent(in)
         !
         type(tEcBCBlock), pointer             :: bcblockPtr   !< intent(out), the new BCBlock
         bcblockId = ec_undef_int
         bcblockPtr => ecBCBlockCreate(instancePtr%idCounter + 1)
         if (associated(bcblockPtr)) then
            ! ensure capacity
            if (instancePtr%nBCBlocks == size(instancePtr%ecBCBlocksPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecBCBlocksPtr, instancePtr%nBCBlocks)) then
                  return
               end if
            end if
            ! register the BCBlock
            instancePtr%nBCBlocks = instancePtr%nBCBlocks + 1
            instancePtr%ecBCBlocksPtr(instancePtr%nBCBlocks)%ptr => bcblockPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            bcblockId = instancePtr%idCounter
         end if
      end function ecInstanceCreateBCBlock

      ! =======================================================================
      !> Create a new NetCDF instance and register it with the Instance (analog of filereader creation).
      function ecInstanceCreateNetCDF(instancePtr) result(netCDFId)
      use m_ec_netcdf_timeseries
      implicit none
         integer                               :: netCDFId    !< id of the new BCBlock
         type(tEcInstance), pointer            :: instancePtr !< intent(in)
         !
         type(tEcNetCDF), pointer              :: netCDFPtr   !< intent(out), the new BCBlock
         netCDFId = ec_undef_int
         netCDFPtr => ecNetCDFCreate(instancePtr%idCounter + 1)
         if (associated(netCDFPtr)) then
            ! ensure capacity
            if (instancePtr%nNetCDFs == size(instancePtr%ecNetCDFsPtr)) then
               if (.not. ecArrayIncrease(instancePtr%ecNetCDFsPtr, instancePtr%nNetCDFs)) then
                  return
               end if
            end if
            ! register the BCBlock
            instancePtr%nNetCDFs = instancePtr%nNetCDFs + 1
            instancePtr%ecNetCDFsPtr(instancePtr%nNetCDFs)%ptr => netCDFPtr
            instancePtr%idCounter = instancePtr%idCounter + 1
            netCDFId = instancePtr%idCounter
         end if
      end function ecInstanceCreateNetCDF

      ! =======================================================================

      !>
      subroutine ecInstanceListSourceItems(instancePtr,dev)
         implicit none
         ! List source items by quantity and location
         type(tEcInstance), pointer :: instancePtr           !< EC instance
         integer, intent(in)        :: dev
         type(tEcItem), pointer     :: sourceItemPtr
         integer  :: ii
         do ii=1, instancePtr%nItems
            sourceItemPtr => instancePtr%ecItemsPtr(ii)%ptr
            if (sourceItemPtr%role == itemType_source) then
                     write(dev,'(a,i4.4,a,i1,a)') 'Source Item ',sourceItemPtr%id
                     write(dev,'(a,i4.4,a,i1,a)') '  Quantity = '//trim(sourceItemPtr%quantityPtr%name)
                     write(dev,'(a,i4.4,a,i1,a)') '  Location = '//trim(sourceItemPtr%elementsetPtr%name)
                     write(dev,*) ''
            endif
         enddo
      end subroutine ecInstanceListSourceItems

      ! =======================================================================
      subroutine ecInstanceListFileReaders(instancePtr,dev)
         implicit none
         ! List source items by quantity and location
         type(tEcInstance), pointer :: instancePtr           !< EC instance
         integer, intent(in)        :: dev
         type(tEcFileReader), pointer     :: fileReaderPtr
         type(tEcItem)      , pointer     :: sourceItemPtr
         character(len=:)   , allocatable :: filename
         integer  :: ii, jj

         do ii=1, instancePtr%nFileReaders
            fileReaderPtr => instancePtr%ecFileReadersPtr(ii)%ptr
            select case (fileReaderPtr%ofType)
            case (provFile_uniform)
               filename = trim(fileReaderPtr%fileName)
            case (provFile_bc)
               filename = trim(fileReaderPtr%bc%fname)
            case default
               filename = ''
            end select
            write(dev,'(a,i4.4,a,i1,a)') 'Filereader ',fileReaderPtr%id,' ('''//filename//''') provides items: '
            do jj=1, fileReaderPtr%nItems
               sourceItemPtr => fileReaderPtr%items(jj)%ptr
               write(dev,'(a,i4.4,a,i1,a)') '   Item ',sourceItemPtr%id
               write(dev,'(a,i4.4,a,i1,a)') '      Quantity = '//trim(sourceItemPtr%quantityPtr%name)
               write(dev,'(a,i4.4,a,i1,a)') '      Location = '//trim(sourceItemPtr%elementsetPtr%name)
            enddo
         write(dev,*) ''
         enddo
      end subroutine ecInstanceListFileReaders
      ! =======================================================================

      subroutine ecInstancePrintState(instancePtr,messenger,lvl)
      use m_ec_message
      implicit none
         type(tEcInstance), pointer :: instancePtr           !< intent(in)
         interface
            subroutine messenger(lvl,msg)
            integer, intent(in)             :: lvl
            character(len=*), intent(in)    :: msg
            end subroutine
         end interface
         integer, intent(in)        :: lvl

         integer                      :: ii, ic, js, i, j
         type(tEcItem),       pointer :: targetItemPtr
         type(tEcConnection), pointer :: connectionPtr
         type(tEcItem),       pointer :: sourceItemPtr
         type(tEcFileReader), pointer :: FileReaderPtr
         type(tEcBCBlock),    pointer :: BCBlockPtr
         character(len=maxMessageLen) :: line

         call messenger(lvl, '.')
         call messenger(lvl, 'EC MODULE ITEMS:')
         do ii=1, instancePtr%nItems
            ! TODO: This lookup loop of items may be expensive for large models, use a lookup table with ids.
            targetItemPtr => instancePtr%ecItemsPtr(ii)%ptr
            if (targetItemPtr%role == itemType_target) then
               write(line,'(a,i4.4,a,i1,a)') 'Target Item ', targetItemPtr%id, ' (name='//trim(targetItemPtr%quantityPtr%name)//', vectormax=',targetItemPtr%quantityPtr%vectormax,')'
               call messenger(lvl, line)
               do ic=1, targetItemPtr%nConnections
                  connectionPtr => targetItemPtr%connectionsPtr(ic)%ptr
                  write(line,'(a,i4.4)') '   Connection ',connectionPtr%id
                  call messenger(lvl, line)
                  do js=1, connectionPtr%nSourceItems
                     sourceItemPtr => connectionPtr%sourceItemsPtr(js)%ptr
                     write(line,'(a,i4.4,a,i1,a)') '      Source Item ',sourceItemPtr%id, ' (name='//trim(sourceItemPtr%quantityPtr%name)//', vectormax=',sourceItemPtr%quantityPtr%vectormax,')'
                     call messenger(lvl, line)
                     ! Find the FileReader which can update this source Item.
                     frs: do i=1, instancePtr%nFileReaders
                        do j=1, instancePtr%ecFileReadersPtr(i)%ptr%nItems
                           if (instancePtr%ecFileReadersPtr(i)%ptr%items(j)%ptr%id == sourceItemPtr%id) then
                              fileReaderPtr => instancePtr%ecFileReadersPtr(i)%ptr
                              if (associated(fileReaderPtr%bc)) then
                                 BCBlockPtr => fileReaderPtr%bc
                                 write(line,'(a,i4.4,a)') '         File Reader ',fileReaderPtr%id, '(filename='//trim(fileReaderPtr%bc%fname)//')'
                                 call messenger(lvl, line)
                                 write(line,'(a,i4.4)') '            BCBlock ',BCBlockPtr%id
                                 call messenger(lvl, line)
                              else
                                 write(line,'(a,i4.4,a)') '         File Reader ',fileReaderPtr%id, '(filename='//trim(fileReaderPtr%filename)//')'
                                 call messenger(lvl, line)
                              end if
                              if (associated(sourceItemPtr%QuantityPtr)) then
                              !  if (allocated(sourceItemPtr%QuantityPtr%name))  write(dev,'(a)') '            Quantity = '//trim(sourceItemPtr%QuantityPtr%name)
                                 if (len_trim(sourceItemPtr%QuantityPtr%name)>0) then
                                    write(line,'(a)') '            Quantity = '//trim(sourceItemPtr%QuantityPtr%name)
                                    call messenger(lvl, line)
                                 end if
                              end if
                              if (associated(sourceItemPtr%ElementSetPtr)) then
                                 if (len_trim(sourceItemPtr%ElementSetPtr%name)>0) then
                                    write(line,'(a)') '            Location = '//trim(sourceItemPtr%ElementSetPtr%name)
                                    call messenger(lvl, line)
                                 end if
                              end if
                              exit frs ! exits outer named do loop
                           end if
                        end do
                     end do frs
                  enddo ! SOURCE ITEMS
               enddo	! CONNECTIONS
            endif
         enddo ! TARGET ITEMS OF INSTANCE
         call messenger(lvl,'.')
      end subroutine ecInstancePrintState

end module m_ec_instance

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!> sorting and reordering
module modSort
  private

  !> generic sort
  interface sort
    module procedure::qsortI
  end interface

  public::sort

contains

  !> quicksort integer array (driver)
  subroutine qsortI(a,perm)
    integer,intent(inout)::a(:) !< the array
    integer,intent(inout),optional::perm(*) !< the permutation array

    if(present(perm))then
      forall(i=1:size(a))
        perm(i)=i
      end forall
      call qsortI_actual(a,perm=perm,first=1,last=size(a))
    else
      call qsortI_actual(a,first=1,last=size(a))
    end if
  end subroutine

  !> quicksort integer array (actual)
  recursive subroutine qsortI_actual(a,perm,first,last)
    integer,intent(inout)::a(*) !< the array
    integer,intent(inout),optional::perm(*) !< the permutation array
    integer,intent(in)::first,last !< sort section from first to last
    integer::tmp,pivot

    pivot=a((first+last)/2)
    i=first
    j=last
    do while(.true.)
      do while(a(i)<pivot)
        i=i+1
      end do
      do while(pivot<a(j))
        j=j-1
      end do
      if(i>=j) exit
      tmp=a(i)
      a(i)=a(j)
      a(j)=tmp
      if(present(perm))then
        tmp=perm(i)
        perm(i)=perm(j)
        perm(j)=tmp
      end if
      i=i+1
      j=j-1
    end do
    if(present(perm))then
      if(first<i-1) call qsortI_actual(a,perm=perm,first=first,last=i-1)
      if(j+1<last) call qsortI_actual(a,perm=perm,first=j+1,last=last)
    else
      if(first<i-1) call qsortI_actual(a,first=first,last=i-1)
      if(j+1<last) call qsortI_actual(a,first=j+1,last=last)
    end if
  end subroutine

end module

max_stars_repo_name: mianzhi/fosolvers

max_issues_repo_name: mianzhi/fosolvers

max_forks_repo_name: mianzhi/fosolvers