SUBROUTINE MESGER90(MAXRLO,IDIMA,A,IMIN,IMUT,XMIN,XMAX,DRMIN,FAC,NRLOC,NSERIE,IRLOC)

implicit none
INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision
real(DP) :: A(IDIMA),xmax,drminx,drmin,drmixx,iflag,indr
real(DP) :: radx,im,iflagx,b1,b2,fac,xmin
integer :: i,j,nrloc,nserie,idown,idima,maxrlo,imin,imut,IRLOC(MAXRLO)
save

      A(1)=0.
      A(2)=0.
!-Initialise IRLOC(i).

      DO 10 I=1,MAXRLO
       IRLOC(I)=0
  10  CONTINUE
     ! IRLOC(1:MAXRLO)=0

!-Distribute as equally as possible the NSERIE boreholes into the NRLOC
!-radial subregions (zones). There will be NRLOC subregions with each
!-IRLoc[j] borehole(s) connected in serie. At the best, NRLOC = NSERIE and
!-IRLoc[j]=1, j=1 to NRLOC. The finest resolution (smaller IRLoc[j]) is
!-set at the border of the store volume (xmax=radius).
      DO 20 I=1,NSERIE
       J= MOD(I-1,NRLOC) + 1
       IRLOC(J)=IRLOC(J) + 1
  20  CONTINUE
!-Calculate the number of boreholes in serie from the center of the store
!-to the upper limit of each radial subregion (zone).

      DO 30 I=2,NRLOC
       IRLOC(I)=IRLOC(I-1)+IRLOC(I)
  30  CONTINUE
    !   IRLOC(2:NRLOC)=IRLOC(1:NRLOC-1)+IRLOC(2:NRLOC)

!-Calculate the lower radial mesh indice of each radial subregion. The
!-results are stored in IRLoc(i).
      DO 40 I=1,NRLOC
       IDOWN=NRLOC+1-I
       IRLOC(IDOWN+1)=IRLOC(IDOWN)+IMIN
  40  CONTINUE
      IRLOC(1)=IMIN
!-Generate the mesh according to the number of boreholes connected in
!-series; calculate A(i).
      A(IMIN)=XMIN
      IMUT=IMIN+NSERIE
      A(IMUT)=XMAX
      DO 50 I=1,NSERIE
       A(IMIN+I)=XMAX*SQRT(dble(I)/dble(NSERIE))
  50  CONTINUE
!-Finer radial mesh within each radial division if necessary.
      DRMINX=DRMIN
      DRMIXX=0.
!-----IFLAG=1 to indicate that the smallest mesh has to be repeated twice
      IFLAG=1
      DO 130 I=1,NSERIE
       INDR=NSERIE-I+IMIN
       RADX=A(INDR+1)-A(INDR)
       IF (RADX.LT.(1.0*DRMINX)) THEN
!-------No additional mesh.
        IFLAG=0
        DRMIXX=RADX
       ELSE
!-------The smallest mesh DRMIN is not repeated if RADX is <=3*DRMIN.
        IF (RADX.LE.(3.*DRMIN)) IFLAG=0
!-------Calculate IM, the number of additional meshes.
        IM=0
        IFLAGX=IFLAG
        B1=DRMINX
        B2=RADX
  60    B2=B2-B1
         IF (B2.LT.0) GO TO 70
         IM=IM+1
         IF (IFLAGX.EQ.1) THEN
           IFLAGX=0
          ELSE
           B1=B1*FAC
          END IF
         GO TO 60
  70    CONTINUE
!-------Update IRLOC, the indice array defining the radial subregions.
        DO 80 J=1,NRLOC
         IF (IRLOC(J+1).GT.INDR) IRLOC(J+1)=IRLOC(J+1)+IM
  80    CONTINUE
!-------Shift in the radial array A(j), for j>INDR.
        DO 90 J=1,(IMUT-INDR)
         IDOWN=IMUT+1-J
         A(IDOWN+IM)=A(IDOWN)
  90    CONTINUE
!-------Calculate and add the new radial coordinate A(i).
        IF (IM.GT.0) THEN
          B1=DRMINX
          DO 100 J=1,IM
           IDOWN=INDR+IM+1-J
           A(IDOWN)=A(IDOWN+1)-B1
           IF (IFLAG.EQ.1) THEN
             IFLAG=0
            ELSE
             B1=B1*FAC
            END IF
  100     CONTINUE
          IF ( ( A(INDR+1)-A(INDR) ).LT.( A(INDR+2)-A(INDR+1) ) ) THEN
            A(INDR+1)= .5*(A(INDR)+A(INDR+2))
           END IF
         END IF
!-------Update the total number of meshes.
        IMUT=IMUT+IM
        DRMIXX=A(INDR+1)-A(INDR)
       END IF
!------Calculate DRMINX for the next radial division.
       B1=DRMIN
       B2=A(IMUT)-A(INDR)
  110  B2=B2-B1
        IF (B2.LT.0) GO TO 120
        B1=B1*FAC
        GO TO 110
  120  CONTINUE
       DRMINX=B1
       DRMIXX=2.*DRMIXX
       IF (DRMINX.LT.DRMIXX) DRMINX=DRMIXX
  130 CONTINUE

RETURN
END subroutine mesger90