text:
c     Subroutine to identify complete months in record
c     AJ_Kettle, Nov09/2017

      SUBROUTINE identify_complete_months(l_prod,s_day_date)

      IMPLICIT NONE
c************************************************************************
      INTEGER             :: l_prod
      CHARACTER(LEN=8)    :: s_day_date(100000)

      INTEGER             :: i,j,k,ii,jj,kk

      INTEGER             :: n_monthdays_expect(12)
      INTEGER             :: n_monthdays_expect2(12)
      INTEGER             :: i_vec_year(300)
      REAL                :: f_vec_year(300)
      REAL                :: f_result
      INTEGER             :: i_refyear
      INTEGER             :: i_matrix_monthdays(300,12)
      INTEGER             :: i_matrix_truedaycnt(300,12)
      INTEGER             :: i_loc_year
      INTEGER             :: i_loc_month
      INTEGER             :: ii1,ii2,ii3

      INTEGER             :: l_cnt_month
      INTEGER             :: l_cnt_month_bad
      INTEGER             :: l_cnt_month_good

      CHARACTER(LEN=4)    :: s_day_year(100000)
      CHARACTER(LEN=2)    :: s_day_month(100000)
      CHARACTER(LEN=2)    :: s_day_day(100000)
      REAL                :: f_day_year(100000)
      REAL                :: f_day_month(100000)
      REAL                :: f_day_day(100000)
      INTEGER             :: i_day_year(100000)
      INTEGER             :: i_day_month(100000)
      INTEGER             :: i_day_day(100000)

      CHARACTER(LEN=8)    :: s_date_single
      CHARACTER(LEN=2)    :: s_mon_ref

      INTEGER             :: i_monend(10000)
      INTEGER             :: i_monstart(10000)
      INTEGER             :: i_monlength(10000)
      INTEGER             :: i_yearid(10000)
      INTEGER             :: i_monthid(10000)
      INTEGER             :: n_month

      CHARACTER(LEN=2)    :: s_test2
      CHARACTER(LEN=4)    :: s_test4
      REAL                :: f_test
      INTEGER             :: i_test

c************************************************************************
c     Declare expected number of days in each month

      n_monthdays_expect(1)=31
      n_monthdays_expect(2)=28
      n_monthdays_expect(3)=31
      n_monthdays_expect(4)=30
      n_monthdays_expect(5)=31
      n_monthdays_expect(6)=30
      n_monthdays_expect(7)=31
      n_monthdays_expect(8)=31
      n_monthdays_expect(9)=30
      n_monthdays_expect(10)=31
      n_monthdays_expect(11)=30
      n_monthdays_expect(12)=31

      n_monthdays_expect2(1)=31
      n_monthdays_expect2(2)=29
      n_monthdays_expect2(3)=31
      n_monthdays_expect2(4)=30
      n_monthdays_expect2(5)=31
      n_monthdays_expect2(6)=30
      n_monthdays_expect2(7)=31
      n_monthdays_expect2(8)=31
      n_monthdays_expect2(9)=30
      n_monthdays_expect2(10)=31
      n_monthdays_expect2(11)=30
      n_monthdays_expect2(12)=31

      ii1=0
      ii2=0
      ii3=0

      i_refyear=300

      DO i=1,i_refyear
       i_vec_year(i)=1750+(i-1)
       f_vec_year(i)=1750.0+(i-1.0)

       f_result=MODULO(f_vec_year(i),4.0)
c       print*,'f_result=',f_result

c      Ordinary years
       IF (f_result.NE.0.0) THEN
        DO j=1,12
         i_matrix_monthdays(i,j)=n_monthdays_expect(j)
        ENDDO
        ii1=ii1+1
       ENDIF

c      Leap year condition
       IF (f_result.EQ.0.0) THEN
        DO j=1,12
         i_matrix_monthdays(i,j)=n_monthdays_expect2(j)
        ENDDO
        ii2=ii2+1
       ENDIF

c      Hard wire century years
       IF (i_vec_year(i).EQ.1800.OR.i_vec_year(i).EQ.1900) THEN
        DO j=1,12
         i_matrix_monthdays(i,j)=n_monthdays_expect(j)
        ENDDO
        ii3=ii3+1
       ENDIF

      ENDDO

c      print*,'i_vec_year=',(i_vec_year(i),i=1,5)
c      print*,'f_vec_year=',(f_vec_year(i),i=1,5)
c      print*,'ii1,ii2,ii3=',ii1,ii2,ii3
c************************************************************************
c     Extract component fields of date

      DO i=1,l_prod      !cycle through length of file
       s_date_single=s_day_date(i)

       s_day_year(i) =s_date_single(1:4)
       s_day_month(i)=s_date_single(5:6)
       s_day_day(i)  =s_date_single(7:8)
c******
c      Convert string to float
       s_test4=s_day_year(i)
       READ(s_test4,*) f_test
       f_day_year(i)=f_test

       s_test2=s_day_month(i)
       READ(s_test2,*) f_test
       f_day_month(i)=f_test

       s_test2=s_day_year(i)
       READ(s_test2,*) f_test
       f_day_day(i)=f_test
c******
c      Convert string to integer
       s_test4=s_day_year(i)
       READ(s_test4,*) i_test
       i_day_year(i)=i_test

       s_test2=s_day_month(i)
       READ(s_test2,*) i_test
       i_day_month(i)=i_test

       s_test2=s_day_year(i)
       READ(s_test2,*) i_test
       i_day_day(i)=i_test
c******
      ENDDO

c      print*,'f_day_year=',(f_day_year(i),i=1,10)
c************************************************************************
c     Count days on standard matrix

c     Initial i_matrix_truedaycnt
      DO j=1,300
       DO k=1,12
        i_matrix_truedaycnt(j,k)=0
       ENDDO
      ENDDO

c     Increment counter
      DO i=1,l_prod      !cycle through length of file
       i_loc_year =i_day_year(i)-1750+1
       i_loc_month=i_day_month(i)

       i_matrix_truedaycnt(i_loc_year,i_loc_month)=
     +    i_matrix_truedaycnt(i_loc_year,i_loc_month)+1

       IF (i_matrix_truedaycnt(i_loc_year,i_loc_month).GT.31) THEN
        print*,'emergency stop at point, month days exceeded',
     +   i_day_year(i),i_day_month(i),
     +   i_matrix_truedaycnt(i_loc_year,i_loc_month)
        CALL SLEEP(50)
       ENDIF
      ENDDO

c     Test to see in number of days in true counts exceeds 31
      DO j=1,300
       DO k=1,12
        IF (i_matrix_truedaycnt(j,k).GT.31) THEN
         print*,'emergency stop, month days exceeded',j,k,
     +     i_matrix_truedaycnt(j,k)
c         CALL SLEEP(5)
        ENDIF
       ENDDO
      ENDDO

c     Count number of months with deficient days
      l_cnt_month     =0
      l_cnt_month_bad =0
      l_cnt_month_good=0
      DO j=1,300
       DO k=1,12
        IF (i_matrix_truedaycnt(j,k).GT.0) THEN
         l_cnt_month=l_cnt_month+1

         IF (i_matrix_truedaycnt(j,k).EQ.i_matrix_monthdays(j,k)) THEN
          l_cnt_month_good=l_cnt_month_good+1
         ENDIF
         IF (i_matrix_truedaycnt(j,k).LT.i_matrix_monthdays(j,k)) THEN
          l_cnt_month_bad =l_cnt_month_bad+1
         ENDIF
        ENDIF
       ENDDO
      ENDDO

c      print*,'matrix comparison l_cnt_month...',
c     +  l_cnt_month,l_cnt_month_good,l_cnt_month_bad
c      IF (l_cnt_month_bad.GT.0) THEN
c       CALL SLEEP(5)
c      ENDIF
c************************************************************************
c     Test if first year before 1750
      IF (i_day_year(1).LT.1750) THEN
       print*,'emergency stop; 1st year before 1750',i_day_year(1)
       CALL SLEEP(100)
      ENDIF
c************************************************************************
c     Count number of month changes in record
      ii=1
      i_monstart(ii)=ii         !first month start indicator
      s_mon_ref=s_day_month(1)
      i_yearid(ii) =i_day_year(1)
      i_monthid(ii)=i_day_month(1)

      DO i=2,l_prod

c      Condition for month switch
       IF (s_day_month(i).NE.s_mon_ref) THEN
        s_mon_ref=s_day_month(i)
        i_monend(ii)=i-1
        ii=ii+1
        i_monstart(ii)=i
        i_yearid(ii) =i_day_year(i)
        i_monthid(ii)=i_day_month(i)

        IF (ii.GT.10000) THEN
         print*,'Month storage array exceeded'
         CALL SLEEP(10)
        ENDIF
       ENDIF
      ENDDO
      i_monend(ii)=l_prod      !last month end indicator

      n_month=ii

c     Find lengths of months
      DO i=1,n_month
       i_monlength(i)=i_monend(i)-i_monstart(i)+1

c       IF (i_monlength(i).GT.31) THEN
c        print*,'i_monlength exceeded',i_monlength(i)
c        print*,'i_monstart(i)=',i_monstart(i)
c        print*,'i_monend(i)=',  i_monend(i)
c       ENDIF
      ENDDO

c      print*,'ii nmonth=',n_month
c      print*,'i_monstart start=',(i_monstart(i),i=1,5)
c      print*,'i_monend   start=',(i_monend(i),i=1,5)
c      print*,'i_monstart end=',  (i_monstart(i),i=n_month-5,n_month)
c      print*,'i_monend   end=',  (i_monend(i),i=n_month-5,n_month)
c      print*,'i_monlength=',     (i_monlength(i),i=1,5)

c      print*,'i_yearid=',        (i_yearid(i),i=1,14)
c      print*,'i_monthid=',       (i_monthid(i),i=1,14)

c      print*,'s_day_month=',(s_day_month(i),i=1,5)

c************************************************************************
c     Test to see if day counts complete

      DO i=1,n_month         !cycle through sequential month list

c      Cycle through leap year list
       DO j=1,i_refyear
c       test if sequential year matches vec_year
        IF (i_yearid(i).EQ.i_vec_year(j)) THEN
c        Cycle through matrix months
         DO k=1,12
c         Test if match between months
          IF (i_monthid(i).EQ.k) THEN
c          Test if number of days good
           IF (i_monlength(i).EQ.i_matrix_monthdays(j,k)) THEN
c            print*,'good match',j,k,i_matrix_monthdays(j,k)
c            CALL SLEEP(1)
           ENDIF

c           IF (i_monlength(i).NE.i_matrix_monthdays(j,k)) THEN
c            print*,'months incomplete'
c            print*,'i_matrix_monthdays ref=',i_matrix_monthdays(j,k)
c            print*,'year/month index=',j,k
c            print*,'i_monlength test=',      i_monlength(i)
cc            CALL SLEEP(5)
c           ENDIF

          ENDIF
         ENDDO
        ENDIF
       ENDDO                 !close j
      ENDDO

c      CALL SLEEP(5)
c************************************************************************
      RETURN
      END

max_stars_repo_name: ajkettle/glamod-nuim

max_issues_repo_name: ajkettle/glamod-nuim

max_forks_repo_name: ajkettle/glamod-nuim


text:
Module Draw

   use Image

   implicit none

   type point
      integer :: x, y
   end type point

   type bucket
      integer :: yMax, yMin, x, sign, dx, dy, summ
   end type bucket

   interface draw_line
      module procedure draw_lineRGB
      module procedure draw_lineRGBA
   end interface draw_line

   interface draw_rectangle
      module procedure draw_rectangleRGB
      module procedure draw_rectangleRGB_A
      module procedure draw_rectangleRGBA
   end interface

   interface draw_circle
      module procedure draw_circleRGB
      ! module procedure draw_circleRGB_A
      module procedure draw_circleRGBA
   end interface

   interface flood_fill
      module procedure flood_fillRGB
      module procedure flood_fillRGBA
   end interface

Contains

   subroutine draw_lineRGB(img, p1, p2, colour)

      implicit none

      type(RGBimage),      intent(INOUT) :: img
      type(RGB), optional, intent(IN)    :: colour
      type(point),         intent(IN)    :: p1, p2
      type(RGB)                          :: c
      type(point)                        :: p1t, p2t
      integer                            :: x, y, dx, dy, error, ystep
      logical                            :: steep

      if(.not.present(colour))then
         c = RGB(0, 0, 0)
      else
         c = colour
      end if

      p1t = p1
      p2t = p2

      steep = (abs(p2t%y - p1t%y) > abs(p2t%x -p1t%x))

      if(steep)then
         call swap(p1t%x, p1t%y)
         call swap(p2t%x, p2t%y)
      end if

      if(p1t%x > p2t%x)then
         call swap(p1t%x, p2t%x)
         call swap(p1t%y, p2t%y)
      end if

      dx = p2t%x - p1t%x
      dy = abs(p2t%y - p1t%y)
      error = dx/2
      y = p1t%y

      if(p1t%y < p2t%y)then
         ystep = 1
      else
         ystep = -1
      end if

    do x = p1t%x, p2t%x
       if (steep) then
          call set_pixel(img, y, x, c)
       else
          call set_pixel(img, x, y, c)
       end if
       error = error - dy
       if ( error < 0 ) then
          y = y + ystep
          error = error + dx
       end if
    end do
   end subroutine draw_lineRGB

   subroutine draw_lineRGBA(img, p1, p2, colour)

      implicit none

      type(RGBAimage),      intent(INOUT) :: img
      type(RGBA), optional, intent(IN)    :: colour
      type(point),         intent(IN)    :: p1, p2
      type(RGBA)                          :: c
      type(point)                        :: p1t, p2t
      integer                            :: x, y, dx, dy, error, ystep
      logical                            :: steep

      if(.not.present(colour))then
         c = RGBA(0, 0, 0, 0)
      else
         c = colour
      end if

      p1t = p1
      p2t = p2

      steep = (abs(p2t%y - p1t%y) > abs(p2t%x -p1t%x))

      if(steep)then
         call swap(p1t%x, p1t%y)
         call swap(p2t%x, p2t%y)
      end if

      if(p1t%x > p2t%x)then
         call swap(p1t%x, p2t%x)
         call swap(p1t%y, p2t%y)
      end if

      dx = p2t%x - p1t%x
      dy = abs(p2t%y - p1t%y)
      error = dx/2
      y = p1t%y

      if(p1t%y < p2t%y)then
         ystep = 1
      else
         ystep = -1
      end if

    do x = p1t%x, p2t%x
       if (steep) then
          call set_pixel(img, y, x, c)
       else
          call set_pixel(img, x, y, c)
       end if
       error = error - dy
       if ( error < 0 ) then
          y = y + ystep
          error = error + dx
       end if
    end do
   end subroutine draw_lineRGBA

   subroutine draw_rectangleRGB(img, p1, p2, colour, fill)

      implicit none

      type(RGBimage),    intent(INOUT) :: img
      type(RGB),         intent(IN)    :: colour
      type(point),       intent(INOUT) :: p1, p2
      logical, optional, intent(IN)    :: fill
      integer                          :: i, j
      logical                          :: flag, xf, yf

      if(present(fill))then
         if(fill)then
            flag = .TRUE.
         else
            flag = .FALSE.
         end if
      else
         flag = .FALSE.
      end if

      if(p1%x > p2%x)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
         xf = .TRUE.
      end if

      if(p1%y > p2%y.and..not.xf)then
         call swap(p1%y, p2%y)
         yf = .TRUE.
      end if

      if(flag)then
         do j = p1%y , p2%y
            do i = p1%x , p2%x
               call set_pixel(img, i, j, colour)
            end do
         end do
      else
         do i = p1%x , p2%x
            call set_pixel(img, i, p1%y, colour)
            call set_pixel(img, i, p2%y, colour)
         end do

         do i = p1%y , p2%y
            call set_pixel(img, p1%x, i, colour)
            call set_pixel(img, p2%x, i, colour)
         end do
      end if

      if(xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      elseif(yf.and..not.xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      end if
   end subroutine draw_rectangleRGB

   subroutine draw_rectangleRGB_A(img, p1, p2, colour, fill)

      implicit none

      type(RGBAimage),    intent(INOUT) :: img
      type(RGB),         intent(IN)     :: colour
      type(point),       intent(INOUT)  :: p1, p2
      logical, optional, intent(IN)     :: fill
      type(RGBA)                        :: c_a
      integer                           :: i, j
      logical                           :: flag, xf, yf

      c_a = RGBA(colour%red, colour%green, colour%green, 255)

      if(present(fill))then
         if(fill)then
            flag = .TRUE.
         else
            flag = .FALSE.
         end if
      else
         flag = .FALSE.
      end if

      if(p1%x > p2%x)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
         xf = .TRUE.
      end if

      if(p1%y > p2%y.and..not.xf)then
         call swap(p1%y, p2%y)
         yf = .TRUE.
      end if

      if(flag)then
         do j = p1%y , p2%y
            do i = p1%x , p2%x
               call set_pixel(img, i, j, c_a)
            end do
         end do
      else
         do i = p1%x , p2%x
            call set_pixel(img, i, p1%y, c_a)
            call set_pixel(img, i, p2%y, c_a)
         end do

         do i = p1%y , p2%y
            call set_pixel(img, p1%x, i, c_a)
            call set_pixel(img, p2%x, i, c_a)
         end do
      end if

      if(xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      elseif(yf.and..not.xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      end if
   end subroutine draw_rectangleRGB_A

   subroutine draw_rectangleRGBA(img, p1, p2, colour, fill, blend)

      implicit none

      type(RGBAimage),   intent(INOUT) :: img
      type(RGBA),        intent(IN)    :: colour
      type(point),       intent(INOUT) :: p1, p2
      logical, optional, intent(IN)    :: fill, blend
      type(RGBA)                       :: c1, c2
      integer                          :: i, j
      logical                          :: flag_f, flag_b, xf, yf

      if(present(blend))then
         if(blend)then
            flag_b = .TRUE.
         else
            flag_b = .FALSE.
         end if
      else
         flag_b = .FALSE.
      end if

      if(present(fill))then
         if(fill)then
            flag_f = .TRUE.
         else
            flag_f = .FALSE.
         end if
      else
         flag_f = .FALSE.
      end if

      if(p1%x > p2%x)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
         xf = .TRUE.
      end if

      if(p1%y > p2%y.and..not.xf)then
         call swap(p1%y, p2%y)
         yf = .TRUE.
      end if

      if(flag_f.and.flag_b)then
         do j = p1%y , p2%y
            do i = p1%x , p2%x
               call get_pixel(img, i, j, c1)
               c2 = alpha_comp(colour, c1)
               call set_pixel(img, i, j, c2)
            end do
         end do
      elseif(flag_f.and..not.flag_b)then
         do j = p1%y , p2%y
            do i = p1%x , p2%x
               call set_pixel(img, i, j, colour)
            end do
         end do
      else
         do i = p1%x , p2%x
            call set_pixel(img, i, p1%y, colour)
            call set_pixel(img, i, p2%y, colour)
         end do

         do i = p1%y , p2%y
            call set_pixel(img, p1%x, i, colour)
            call set_pixel(img, p2%x, i, colour)
         end do
      end if

      if(xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      elseif(yf.and..not.xf)then
         call swap(p1%x, p2%x)
         call swap(p1%y, p2%y)
      end if
   end subroutine draw_rectangleRGBA

   subroutine draw_circleRGB(img, p, radius, colour, fill)

      implicit none

      type(RGBimage),    intent(INOUT) :: img
      type(RGB),         intent(IN)    :: colour
      type(point),       intent(IN)    :: p
      integer,           intent(IN)    :: radius
      logical, optional, intent(IN)    :: fill
      integer                          :: x, y, error

      x = radius
      y = 0
      error = 0

      do while( x >= y)

         call set_pixel(img, p%x + x, p%y + y , colour)
         call set_pixel(img, p%x + y, p%y + x , colour)
         call set_pixel(img, p%x - y, p%y + x , colour)
         call set_pixel(img, p%x - x, p%y + y , colour)
         call set_pixel(img, p%x - x, p%y - y , colour)
         call set_pixel(img, p%x - y, p%y - x , colour)
         call set_pixel(img, p%x + y, p%y - x , colour)
         call set_pixel(img, p%x + x, p%y - y , colour)

         y = y + 1
         error = error + 1 + 2*y
         if(2*(error-x) + 1 > 0)then
            x = x - 1
            error = error + (1 - 2*x)
         end if
      end do

      if(present(fill))then
         if(fill)then
            call flood_fill(img, p%x, p%y, colour, RGB(img%Red(x, y), img%Green(x, y), img%Blue(x, y)))
         end if
      end if

   end subroutine draw_circleRGB

   subroutine draw_circleRGBA(img, p, radius, colour, fill, blend)

      implicit none

      type(RGBAimage),   intent(INOUT) :: img
      type(RGBA),        intent(IN)    :: colour
      type(point),       intent(IN)    :: p
      integer,           intent(IN)    :: radius
      logical, optional, intent(IN)    :: fill, blend
      type(RGBA)                       :: c1, c2
      integer                          :: x, y, error,i,j

      x = radius
      y = 0
      error = 0
      c1 = colour

      if(present(fill).and.present(blend))then
         if(fill.and.blend)then

            do i = -radius, radius, 1
               do j = -radius, radius, 1
                  if(i*i + j*j <= radius*radius)then
                     call get_pixelRGBA(img, p%x+i, p%y+j, c2)
                     c1 = alpha_comp(colour, c2)
                     call set_pixel(img, p%x+i, p%y+j, c1)
                  end if
               end do
            end do
         elseif(fill)then
            do i = -radius, radius, 1
               do j = -radius, radius, 1
                  if(i*i + j*j <= radius*radius)then
                     call set_pixel(img, p%x+i, p%y+j, c1)
                  end if
               end do
            end do
         elseif(blend)then
            do while(x >= y)

               call get_pixelRGBA(img, p%x+x, p%y+y, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x + x, p%y + y , c1)

               call get_pixelRGBA(img, p%x+y, p%y+x, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x + y, p%y + x , c1)

               call get_pixelRGBA(img, p%x-y, p%y+x, c2)
               c1 = alpha_comp(colour,c2)
               call set_pixel(img, p%x - y, p%y + x , c1)

               call get_pixelRGBA(img, p%x-x, p%y+y, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x - x, p%y + y , c1)

               call get_pixelRGBA(img, p%x-x, p%y-y, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x - x, p%y - y , c1)

               call get_pixelRGBA(img, p%x-y, p%y-x, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x - y, p%y - x , c1)

               call get_pixelRGBA(img, p%x+y, p%y-x, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x + y, p%y - x , c1)

               call get_pixelRGBA(img, p%x+x, p%y-y, c2)
               c1 = alpha_comp(colour, c2)
               call set_pixel(img, p%x + x, p%y - y , c1)

               y = y + 1
               error = error + 1 + 2*y
               if(2*(error-x) + 1 > 0)then
                  x = x - 1
                  error = error + (1 - 2*x)
               end if
            end do
         else
            do while(x >= y)

               call set_pixel(img, p%x + x, p%y + y , c1)
               call set_pixel(img, p%x + y, p%y + x , c1)
               call set_pixel(img, p%x - y, p%y + x , c1)
               call set_pixel(img, p%x - x, p%y + y , c1)
               call set_pixel(img, p%x - x, p%y - y , c1)
               call set_pixel(img, p%x - y, p%y - x , c1)
               call set_pixel(img, p%x + y, p%y - x , c1)
               call set_pixel(img, p%x + x, p%y - y , c1)

               y = y + 1
               error = error + 1 + 2*y
               if(2*(error-x) + 1 > 0)then
                  x = x - 1
                  error = error + (1 - 2*x)
               end if
            end do
         end if
      elseif(present(fill))then
         if(fill)then
            do i = -radius, radius, 1
               do j = -radius, radius, 1
                  if(i*i + j*j <= radius*radius)then
                     call get_pixelRGBA(img, p%x+i, p%y+j, c2)
                     c1 = alpha_comp(colour, c2)
                     call set_pixel(img, p%x+i, p%y+j, c1)
                  end if
               end do
            end do
         else
            do while( x >= y)

               call set_pixel(img, p%x + x, p%y + y , colour)
               call set_pixel(img, p%x + y, p%y + x , colour)
               call set_pixel(img, p%x - y, p%y + x , colour)
               call set_pixel(img, p%x - x, p%y + y , colour)
               call set_pixel(img, p%x - x, p%y - y , colour)
               call set_pixel(img, p%x - y, p%y - x , colour)
               call set_pixel(img, p%x + y, p%y - x , colour)
               call set_pixel(img, p%x + x, p%y - y , colour)

               y = y + 1
               error = error + 1 + 2*y
               if(2*(error-x) + 1 > 0)then
                  x = x - 1
                  error = error + (1 - 2*x)
               end if
            end do
         end if
      else
         do while( x >= y)

               call set_pixel(img, p%x + x, p%y + y , colour)
               call set_pixel(img, p%x + y, p%y + x , colour)
               call set_pixel(img, p%x - y, p%y + x , colour)
               call set_pixel(img, p%x - x, p%y + y , colour)
               call set_pixel(img, p%x - x, p%y - y , colour)
               call set_pixel(img, p%x - y, p%y - x , colour)
               call set_pixel(img, p%x + y, p%y - x , colour)
               call set_pixel(img, p%x + x, p%y - y , colour)

               y = y + 1
               error = error + 1 + 2*y
               if(2*(error-x) + 1 > 0)then
                  x = x - 1
                  error = error + (1 - 2*x)
               end if
         end do
      end if

   end subroutine draw_circleRGBA

   subroutine draw_polygon(img, colour, p1, p2, p3, fill)

      implicit none

      type(RGBaimage),    intent(INOUT) :: img
      type(RGBa),         intent(IN)    :: colour
      type(point),       intent(INOUT) :: p1, p2, p3
      logical, optional, intent(IN)    :: fill
      integer                          :: x, y
      type(bucket)                     :: buck(3)

      call draw_line(img, p1, p2, colour)
      call draw_line(img, p2, p3, colour)
      call draw_line(img, p3, p1, colour)

      do x = 1 , 3
         buck(x) = bucket(max(p1%y,p2%y), min(p1%y,p2%y), x, 1, abs(p2%x-p1%x), abs(p2%y-p1%y), 0)
      end do

      if(present(fill))then
         if(fill)then
            x = (p1%x + p2%x + p3%x)/3
            y = (p1%y + p2%y + p3%y)/3
            print*,x,y
!            call set_pixel(img, x,y,RGB(255,255,255))
            call flood_fill(img, x, y, colour, RGBa(img%Red(x, y), img%Green(x, y), img%Blue(x, y),img%alpha(x,y)))
         end if
      end if

   end subroutine draw_polygon

   subroutine swap(a, b)

      implicit none

      integer, intent(INOUT) :: a, b
      integer                :: tmp

      tmp = a
      a = b
      b = tmp

   end subroutine swap

   ! subroutine fill_poly(img, colour)

   !    implicit none

   !    type(RGBAimage), intent(INOUT) :: img
   !    type(RGBA),      intent(IN)    :: colour

   ! end subroutine fill_poly

   recursive subroutine flood_fillRGB(img, x, y, colour, old)

      implicit none

      type(RGBimage), intent(INOUT) :: img
      type(RGB),      intent(IN)    :: colour, old
      integer,        intent(IN)    :: x, y
      integer                       :: x1

      if(old == colour)return
      if( (RGB(img%Red(x, y), img%Green(x, y), img%Blue(x, y)) /= old))return

      x1 = x
      do while(x1 < img%width .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                             img%Blue(x1, y)) == old)
         call set_pixel(img, x1, y, colour)
         x1 = x1 + 1
      end do

      x1 = x - 1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                                img%Blue(x1, y)) == old)
            call set_pixel(img, x1, y, colour)

            x1 = x1 - 1
            if(x1==0)exit
         end do
      end if

!      move up
      x1 = x
      do while(x1 < img%width .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                             img%Blue(x1, y)) == colour)
         if(y==1)exit
         if(y > 1 .and. RGB(img%Red(x1, y-1), img%Green(x1, y-1), &
                                             img%Blue(x1, y-1))==old)then
            call flood_fill(img, x, y-1, colour, old)
         end if
         x1 = x1 +1

      end do

      x1 = x - 1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                                img%Blue(x1, y)) == colour)
            if(x1==0 .or. y==1)exit
            if(y > 1 .and. RGB(img%Red(x1, y-1), img%Green(x1, y-1), &
                                                img%Blue(x1, y-1))==old)then
               call flood_fill(img, x, y-1, colour, old)
            end if
            x1 = x1 - 1
            if(x1==0 .or. y==0)exit
         end do
      end if
      !move down

      x1 = x
      do while(x1 < img%width .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                             img%Blue(x1, y)) == colour)
         if(y==img%height)exit
         if(y < img%height-1 .and. RGB(img%Red(x1, y+1), img%Green(x1, y+1), &
                                             img%Blue(x1, y+1))==old)then
            call flood_fill(img, x, y+1, colour, old)
         end if
         x1 = x1 +1
      end do

      x1 = x-1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGB(img%Red(x1, y), img%Green(x1, y), &
                                                img%Blue(x1, y)) == colour)
            if(y==img%height)exit
            if(y < img%height .and. RGB(img%Red(x1, y+1), img%Green(x1, y+1), &
                                                img%Blue(x1, y+1))==old)then
               call flood_fill(img, x, y+1, colour, old)
            end if
            x1 = x1 - 1
            if(x1==0.or.y==0)exit
         end do
      end if
   end subroutine flood_fillRGB

   recursive subroutine flood_fillRGBA(img, x, y, colour, old)

      implicit none

      type(RGBAimage), intent(INOUT) :: img
      type(RGBA),      intent(IN)    :: colour, old
      integer,         intent(IN)    :: x, y
      type(RGBA)                     :: c
      integer                        :: x1

      if(old == colour)return
      call get_pixelRGBA(img, x, y, c)
      if(c /= old)return
      x1 = x
      do while(x1 < img%width .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                         img%Blue(x1, y), img%alpha(x1, y)) == old)
         call set_pixel(img, x1, y, colour)
         x1 = x1 + 1
      end do

      x1 = x - 1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                     img%Blue(x1, y), img%alpha(x1, y)) == old)
            call set_pixel(img, x1, y, colour)

            x1 = x1 - 1
            if(x1==0)exit
         end do
      end if

!      move up
      x1 = x
      do while(x1 < img%width .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                         img%Blue(x1, y), img%alpha(x1, y)) == colour)
         if(y==1)exit
         if(y > 1 .and. RGBA(img%Red(x1, y-1), img%Green(x1, y-1), &
                             img%Blue(x1, y-1), img%alpha(x1, y-1)) == old)then
            call flood_fill(img, x, y-1, colour, old)
         end if
         x1 = x1 +1

      end do

      x1 = x - 1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                     img%Blue(x1, y), img%alpha(x1, y)) == colour)
            if(x1==0 .or. y==1)exit
            if(y > 1 .and. RGBA(img%Red(x1, y-1), img%Green(x1, y-1), &
                                img%Blue(x1, y-1), img%alpha(x1, y-1)) == old)then
               call flood_fill(img, x, y-1, colour, old)
            end if
            x1 = x1 - 1
            if(x1==0 .or. y==0)exit
         end do
      end if
      !move down

      x1 = x
      do while(x1 < img%width .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                         img%Blue(x1, y), img%alpha(x1, y)) == colour)
         if(y==img%height)exit
         if(y < img%height-1 .and. RGBA(img%Red(x1, y+1), img%Green(x1, y+1), &
                                        img%Blue(x1, y+1), img%alpha(x1, y+1)) == old)then
            call flood_fill(img, x, y+1, colour, old)
         end if
         x1 = x1 +1
      end do

      x1 = x-1
      if(x1/=0)then
         do while(x1 >= 1 .and. RGBA(img%Red(x1, y), img%Green(x1, y), &
                                     img%Blue(x1, y), img%alpha(x1, y)) == colour)
            if(y==img%height)exit
            if(y < img%height .and. RGBA(img%Red(x1, y+1), img%Green(x1, y+1), &
                                         img%Blue(x1, y+1), img%alpha(x1, y+1)) == old)then
               call flood_fill(img, x, y+1, colour, old)
            end if
            x1 = x1 - 1
            if(x1==0.or.y==0)exit
         end do
      end if

   end subroutine flood_fillRGBA

   function alpha_comp(ca, cb) result(co)

        implicit none

        type(RGBA), intent(IN) :: cb, ca
        type(RGBA)             :: co
        real                   :: a_tmp, b_tmp

        a_tmp = ca%alpha/255.
        b_tmp = cb%alpha/255.

         co%red =  clampInt(int(ca%red * a_tmp) + int(cb%red * (1. - a_tmp)), 0, 255)
         co%green =  clampInt(int(ca%green * a_tmp) + int(cb%green * (1. - a_tmp)), 0, 255)
         co%blue =  clampInt(int(ca%blue * a_tmp) + int(cb%blue * (1. - a_tmp)), 0, 255)
         co%alpha = clampInt(int((b_tmp + (1. - b_tmp) * a_tmp)*255.), 0, 255)

   end function alpha_comp

end module Draw

max_stars_repo_name: lewisfish/FortranForm

max_issues_repo_name: lewisfish/FortranForm

max_forks_repo_name: lewisfish/FortranForm


text:
      module ptrd_mod
        use stel_kinds
        USE v3_utilities, ONLY: assert
        implicit none

#if defined(MPI_OPT)
        contains

        subroutine pdtrdf(m,nblock,Amat,Bmat,Cmat,ipivmat, desc)
        implicit none
!
!       perform LU factorization of block tridiagonal system
!
!       [A1  C1       ]   [ L1          ] [ I U1        ]
!       [B1  A2 C2    ]   [ B1 L2       ] [   I  U2     ]
!       [    B2 A3 C3 ] = [    B2 L3    ] [      I   U3 ]
!       [       B4 A4 ]   [       B3 L4 ] [          I  ]
!
!
! L1 = A1
! L1*U1 = C1 => U1 = L1\C1
! B1*U1 + L2 = A2 =>  L2 = A2 - B1*U1
! L2*U2 = C2 =>  U2 = L2\C2
!
! Lk = Ak - B_{k-1}*U_{k-1}
! Uk = Lk\Ck
!
!
        REAL(rprec), parameter :: one = 1.0d0
        REAL(rprec), parameter :: zero = 0.0d0

        integer, intent(in) :: m, nblock
        REAL(rprec), target, dimension(:,:), intent(inout) :: Amat
        REAL(rprec), target, dimension(:,:), intent(inout) :: Bmat
        REAL(rprec), target, dimension(:,:), intent(inout) :: Cmat
        integer, dimension(:,:), target, intent(inout) :: ipivmat
        integer, dimension(:), intent(in) :: desc

        integer :: k, ia,ja, ib,jb, iu,ju, ic,jc
        integer :: info, nrhs, mm,nn
        integer, dimension(:), pointer :: ipiv
        REAL(rprec) :: alpha, beta
        REAL(rprec), dimension(:), pointer :: Ak, Bkm1, Ukm1
        REAL(rprec), dimension(:), pointer :: Lk, Ck

        nullify( Ak )
        nullify( Bkm1 )
        nullify( Ukm1 )
        nullify( Lk )
        nullify( Ck )

        do k=1,nblock
          if (k .ge. 2) then
!            -------------------------
!            Ak = Ak - B_{k-1}*U_{k-1}
!            -------------------------
             Ak => Amat(:,k)
             Bkm1 => Bmat(:,k-1)
             Ukm1 => Cmat(:,k-1)
             alpha = -one
             beta = one
             ib = 1
             jb = 1
             iu = 1
             ju = 1
             ia = 1
             ja = 1

             call pdgemm( 'N', 'N', m,m,m,                                    &
     &                   alpha, Bkm1,ib,jb,desc,                              &
     &                          Ukm1,iu,ju,desc,                              &
     &                   beta,  Ak, ia,ja, desc )
          endif

!         ------------
!         Uk = Lk \ Ck
!         ------------
          Lk => Amat(:,k)
          ipiv => ipivmat(:,k)
          ia = 1
          ja = 1
          info = 0
          mm = m
          nn = m
          call pdgetrf(mm,nn,Lk,ia,ja,desc,ipiv,info)
          call assert(info.eq.0,'pdtrdf: pdgetrf return info != 0')

          if (k .le. (nblock-1)) then
             nrhs = m
             ic = 1
             jc = 1
             Ck => Cmat(:,k)
             call pdgetrs('No',m,nrhs,Lk,ia,ja,desc,ipiv,               &
     &                    Ck,ic,jc,desc, info)
             call assert(info.eq.0,'pdtrdf: pdgetrs return info != 0')
          endif

        enddo
        return
        end subroutine pdtrdf

        subroutine pdtrds(m,nblock,Amat,Bmat,Cmat,ipivmat,desc,         &
     &                    nrhs, Rrhs_in,ir,jr,descRrhs_in)
        use descriptor_mod, ir1=>ir, jr1=>jr, mm1=>mm, info1=>info
        implicit none

!
!       use LU factorization of block tridiagonal system
!
!       [A1  C1       ]   [ L1          ] [ I U1        ]
!       [B1  A2 C2    ]   [ B1 L2       ] [   I  U2     ]
!       [    B2 A3 C3 ] = [    B2 L3    ] [      I   U3 ]
!       [       B4 A4 ]   [       B3 L4 ] [          I  ]
!
!
!
!
        integer, parameter :: idebug = 0
        REAL(rprec), parameter :: one = 1.0d0
        REAL(rprec), parameter :: zero = 0.0d0

        integer, intent(in) :: m, nblock
        REAL(rprec), dimension(:,:), target, intent(inout) :: Amat
        REAL(rprec), dimension(:,:), target, intent(inout) :: Bmat
        REAL(rprec), dimension(:,:), target, intent(inout) :: Cmat
        integer, dimension(:,:), target, intent(inout) :: ipivmat
        integer, dimension(:), intent(in) :: desc

        integer, intent(in) :: nrhs
        integer, intent(in) :: ir,jr
        REAL(rprec), dimension(:) :: Rrhs_in
        integer, dimension(:), intent(in) :: descRrhs_in

        integer :: k, info, mm,nn,kk
        integer :: ia,ja, irk, jrk, iy,jy, ix,jx, iu,ju, ib,jb
        REAL(rprec) :: alpha, beta

        REAL(rprec), pointer, dimension(:) :: Lk,Bkm1,Uk
        REAL(rprec), pointer, dimension(:) :: Rrhsk,Rrhskm1,Rrhskp1
        integer, pointer, dimension(:) :: ipiv

	integer :: inc1,inc2,iblock,irhs

	logical, parameter :: use_pdcopy = .false.
	REAL(rprec), dimension(:,:), target, allocatable :: Rrhs
	integer, dimension(DLEN_) :: descRrhs
    INTEGER           :: NUMROC
    EXTERNAL          :: NUMROC

!	-----------------------------------------------------
!	change storage to nblock copies of m by nrhs matrices
!	also to avoid alignment constraints in scalapack
!	-----------------------------------------------------
	icontxt = descRrhs_in(CTXT_)
	mb = descRrhs_in(MB_)
	nb = 1
	rsrc = 0
	csrc = 0

	call blacs_gridinfo( icontxt, nprow,npcol,myrow,mycol)
	Locp = numroc( m, mb, myrow,rsrc,nprow)
	Locq = numroc( nrhs,nb,mycol,csrc,npcol)
	lld = max(1,Locp)
	call descinit(descRrhs,m,nrhs,mb,nb,rsrc,csrc,icontxt,lld,info)
	call assert(info.eq.0,'pdtrds: descinit return info != 0')

	ineed = max(1,Locp*Locq)
	allocate( Rrhs(ineed,nblock),stat=ierr)
	call assert(ierr.eq.0,'pdtrds: alloc Rrhs,ierr != 0')
        Rrhs = 0.0d0

	nullify( Rrhsk )
	nullify( Rrhskm1 )
	nullify( Rrhskp1 )

!	---------
!	copy data
!	---------

	if (use_pdcopy) then

	do iblock=1,nblock
	do irhs=1,nrhs
	  ia = (ir-1) + 1 + (iblock-1)*m
	  ja = (jr-1) + irhs
	  inc1 = 1

	  ib = 1
	  jb = irhs
	  inc2 = 1

	  Rrhsk => Rrhs(:,iblock)
	  call pdcopy(m, Rrhs_in,ia,ja,descRrhs_in,inc1,                    &
     &                   Rrhsk,ib,jb,descRrhs,inc2 )
	enddo
	enddo

	else

	do iblock=1,nblock
	  irhs = 1
	  ia = (ir-1) + 1 + (iblock-1)*m
	  ja = (jr-1) + irhs

	  ib = 1
	  jb = irhs
	  Rrhsk => Rrhs(:,iblock)
	  alpha = 1.0d0
	  beta = 0.0d0
	  call pdgeadd( 'N',m,nrhs,alpha,Rrhs_in,ia,ja,descRrhs_in,           &
     &            beta, Rrhsk,ib,jb,descRrhs )
	enddo

	endif

!   -----------------
!   L*U*x = r
!   (1) solve L*y = r
!   (2) solve U*x = y
!   -----------------
        nullify(Lk)
        nullify(Bkm1)
        nullify(Uk)
        nullify(ipiv)

        icontxt = desc(CTXT_)
	call blacs_gridinfo(icontxt, nprow,npcol,myrow,mycol)
        isroot = (myrow.eq.0).and.(mycol.eq.0)

!   --------------------------------
!  (1) solve L*y = r
!
!   [L1             ] [ y1 ]   [ r1 ]
!   [B1  L2         ] [ y2 ]   [ r2 ]
!   [    B2  L3     ] [ y3 ] = [ r3 ]
!   [        B3  L4 ] [ y4 ]   [ r4 ]
!
!
!   y1 = L1\r1
!   y2 = L2\( r2 - B1*y1 )
!   y3 = L3\( r3 - B2*y2 )
!   y4 = L4\( r4 - B3*y3 )
!
!   yk = Lk\( rk - B_{k-1}*y_{k-1} )
!   --------------------------------
          if (isroot .and. (idebug.ge.1)) then
             write(*,*) 'pdtrds 77: m,nblock,nrhs ',m,nblock,nrhs
             write(*,*) 'descRrhs_in(M_) ',descRrhs_in(M_)
             write(*,*) 'descRrhs_in(N_) ',descRrhs_in(N_)
             write(*,*) 'descRrhs_in(MB_) ',descRrhs_in(MB_)
             write(*,*) 'descRrhs_in(NB_) ',descRrhs_in(NB_)
          endif

        do k=1,nblock
          if (k .ge. 2) then
!           --------------------------
!           rk <- rk - B_{k-1}*y_{k-1}
!           --------------------------
            Bkm1 => Bmat(:,k-1)

            alpha = -one
            beta = one
            mm = m
            nn = nrhs
            kk = m
            ib = 1
            jb = 1

            iy = 1
            jy = 1
            irk = 1
            jrk = 1
	    Rrhsk => Rrhs(:,k)
	    Rrhskm1 => Rrhs(:,k-1)

            call pdgemm( 'N', 'N', mm,nn,kk,                                 &
     &          alpha,  Bkm1, ib,jb, desc,                                   &
     &                  Rrhskm1, iy,jy, descRrhs,                            &
     &          beta,   Rrhsk, irk,jrk,descRrhs )
          endif

!         ------------
!         yk = Lk \ rk
!         ------------

          Lk => Amat(:,k)
          ia = 1
          ja = 1
          ipiv => ipivmat(:,k)

	  Rrhsk => Rrhs(:,k)
          jrk = 1
          irk = 1

          info = 0
          if (isroot .and. (idebug.ge.1)) then
             write(*,*) 'pdtrds 106: k,irk,jrk ',k,irk,jrk
          endif
          call pdgetrs( 'N', m,nrhs, Lk,ia,ja,desc, ipiv,               &
     &                  Rrhsk,irk,jrk,descRrhs,info)
          call assert(info.eq.0,'pdtrds: pdgetrs return info != 0')

        enddo

!  (2)  solve U*x = y
!
! [I  U1         ]   [ x1 ]   [ y1 ]
! [   I   U2     ]   [ x2 ]   [ y2 ]
! [       I   U3 ]   [ x3 ] = [ y3 ]
! [           I  ]   [ x4 ]   [ y4 ]
!
!
! x4 = y4
! x3 = y3 - U3*y4
! x2 = y2 - U2*y3
! x1 = y1 - U1*y2
!
! xk = yk - Uk*x_{k+1}
!

        do k=nblock-1,1,-1
!          --------------------
!          xk = yk - Uk*x_{k+1}
!          --------------------
           alpha = -one
           beta = one

           mm = m
           nn = nrhs
           kk = m
           Uk => Cmat(:,k)
           iu = 1
           ju = 1

           ix = 1
           jx = 1
           irk = 1
           jrk = 1

	   Rrhskp1 => Rrhs(:,k+1)
	   Rrhsk => Rrhs(:,k)

           call pdgemm( 'N','N', mm,nn,kk,                                     &
     &                 alpha,  Uk,iu,ju,desc,                                  &
     &                         Rrhskp1,ix,jx,descRrhs,                         &
     &                 beta,   Rrhsk,irk,jrk,descRrhs )
         enddo

!	----------------
!	copy results out
!	----------------

	if (use_pdcopy) then

	do iblock=1,nblock
	do irhs=1,nrhs
	  ia = (ir-1) + 1 + (iblock-1)*m
	  ja = (jr-1) + irhs
	  inc1 = 1

	  ib = 1
	  jb = irhs
	  inc2 = 1

	  Rrhsk => Rrhs(:,iblock)
          call pdcopy(m,Rrhsk,ib,jb,descRrhs,inc2,                          &
     &                  Rrhs_in,ia,ja,descRrhs_in,inc1)
	enddo
	enddo

	else

	do iblock=1,nblock
	  irhs = 1
	  ia = (ir-1) + 1 + (iblock-1)*m
	  ja = (jr-1) + irhs

	  ib = 1
	  jb = irhs
	  Rrhsk => Rrhs(:,iblock)
	  alpha = 1.0d0
	  beta = 0.0d0

	  call pdgeadd('N',m,nrhs,alpha,Rrhsk,ib,jb,descRrhs,               &
     &                            beta,Rrhs_in,ia,ja,descRrhs_in)
	enddo

	endif
	deallocate( Rrhs, stat=ierr)
	call assert(ierr.eq.0,'pdtrds: dealloc Rhs ')

        return
        end subroutine pdtrds

#if defined(NEED_TOOLS)
      SUBROUTINE PDGETRS( TRANS, N, NRHS, A, IA, JA, DESCA, IPIV, B,    &
     &                    IB, JB, DESCB, INFO )
!
!  -- ScaLAPACK routine (version 1.7) --
!     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
!     and University of California, Berkeley.
!     May 1, 1997
!
!     .. Scalar Arguments ..
      CHARACTER          TRANS
      INTEGER            IA, IB, INFO, JA, JB, N, NRHS
!     ..
!     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCB( * ), IPIV( * )
      REAL(dp)   A( * ), B( * )
!     ..
!
!  Purpose
!  =======
!
!  PDGETRS solves a system of distributed linear equations
!
!                   op( sub( A ) ) * X = sub( B )
!
!  with a general N-by-N distributed matrix sub( A ) using the LU
!  factorization computed by PDGETRF.
!  sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1), op( A ) = A or A**T and
!  sub( B ) denotes B(IB:IB+N-1,JB:JB+NRHS-1).
!
!  Notes
!  =====
!
!  Each global data object is described by an associated description
!  vector.  This vector stores the information required to establish
!  the mapping between an object element and its corresponding process
!  and memory location.
!
!  Let A be a generic term for any 2D block cyclicly distributed array.
!  Such a global array has an associated description vector DESCA.
!  In the following comments, the character _ should be read as
!  "of the global array".
!
!  NOTATION        STORED IN      EXPLANATION
!  --------------- -------------- --------------------------------------
!  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
!                                 DTYPE_A = 1.
!  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
!                                 the BLACS process grid A is distribu-
!                                 ted over. The context itself is glo-
!                                 bal, but the handle (the integer
!                                 value) may vary.
!  M_A    (global) DESCA( M_ )    The number of rows in the global
!                                 array A.
!  N_A    (global) DESCA( N_ )    The number of columns in the global
!                                 array A.
!  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
!                                 the rows of the array.
!  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
!                                 the columns of the array.
!  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
!                                 row of the array A is distributed.
!  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
!                                 first column of the array A is
!                                 distributed.
!  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
!                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
!
!  Let K be the number of rows or columns of a distributed matrix,
!  and assume that its process grid has dimension p x q.
!  LOCr( K ) denotes the number of elements of K that a process
!  would receive if K were distributed over the p processes of its
!  process column.
!  Similarly, LOCc( K ) denotes the number of elements of K that a
!  process would receive if K were distributed over the q processes of
!  its process row.
!  The values of LOCr() and LOCc() may be determined via a call to the
!  ScaLAPACK tool function, NUMROC:
!          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
!          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
!  An upper bound for these quantities may be computed by:
!          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
!          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
!
!  This routine requires square block data decomposition ( MB_A=NB_A ).
!
!  Arguments
!  =========
!
!  TRANS   (global input) CHARACTER
!          Specifies the form of the system of equations:
!          = 'N':  sub( A )    * X = sub( B )  (No transpose)
!          = 'T':  sub( A )**T * X = sub( B )  (Transpose)
!          = 'C':  sub( A )**T * X = sub( B )  (Transpose)
!
!  N       (global input) INTEGER
!          The number of rows and columns to be operated on, i.e. the
!          order of the distributed submatrix sub( A ). N >= 0.
!
!  NRHS    (global input) INTEGER
!          The number of right hand sides, i.e., the number of columns
!          of the distributed submatrix sub( B ). NRHS >= 0.
!
!  A       (local input) REAL(dp) pointer into the local
!          memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
!          On entry, this array contains the local pieces of the factors
!          L and U from the factorization sub( A ) = P!L!U; the unit
!          diagonal elements of L are not stored.
!
!  IA      (global input) INTEGER
!          The row index in the global array A indicating the first
!          row of sub( A ).
!
!  JA      (global input) INTEGER
!          The column index in the global array A indicating the
!          first column of sub( A ).
!
!  DESCA   (global and local input) INTEGER array of dimension DLEN_.
!          The array descriptor for the distributed matrix A.
!
!  IPIV    (local input) INTEGER array, dimension ( LOCr(M_A)+MB_A )
!          This array contains the pivoting information.
!          IPIV(i) -> The global row local row i was swapped with.
!          This array is tied to the distributed matrix A.
!
!  B       (local input/local output) REAL(dp) pointer into the
!          local memory to an array of dimension
!          (LLD_B,LOCc(JB+NRHS-1)).  On entry, the right hand sides
!          sub( B ). On exit, sub( B ) is overwritten by the solution
!          distributed matrix X.
!
!  IB      (global input) INTEGER
!          The row index in the global array B indicating the first
!          row of sub( B ).
!
!  JB      (global input) INTEGER
!          The column index in the global array B indicating the
!          first column of sub( B ).
!
!  DESCB   (global and local input) INTEGER array of dimension DLEN_.
!          The array descriptor for the distributed matrix B.
!
!  INFO    (global output) INTEGER
!          = 0:  successful exit
!          < 0:  If the i-th argument is an array and the j-entry had
!                an illegal value, then INFO = -(i!100+j), if the i-th
!                argument is a scalar and had an illegal value, then
!                INFO = -i.
!
!  =====================================================================
!
!     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,   &
     &                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,   &
     &                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,  &
     &                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL(dp)   ONE
      PARAMETER          ( ONE = 1.0D+0 )
!     ..
!     .. Local Scalars ..
      LOGICAL            NOTRAN
      INTEGER            IAROW, IBROW, ICOFFA, ICTXT, IROFFA, IROFFB,    &
     &                   MYCOL, MYROW, NPCOL, NPROW
!     ..
!     .. Local Arrays ..
      INTEGER            DESCIP( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
!     ..
!     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCSET, PCHK2MAT,     &
     &                   PDLAPIV, PDTRSM, PXERBLA
!     ..
!     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2P, NUMROC
      EXTERNAL           INDXG2P, LSAME, NUMROC
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC          ICHAR, MOD
!     ..
!     .. Executable Statements ..
!
!     Get grid parameters
!
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
!
!     Test the input parameters
!
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(700+CTXT_)
      ELSE
         NOTRAN = LSAME( TRANS, 'N' )
         CALL CHK1MAT( N, 2, N, 2, IA, JA, DESCA, 7, INFO )
         CALL CHK1MAT( N, 2, NRHS, 3, IB, JB, DESCB, 12, INFO )
         IF( INFO.EQ.0 ) THEN
            IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),   &
     &                       NPROW )
            IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ),   &
     &                       NPROW )
            IROFFA = MOD( IA-1, DESCA( MB_ ) )
            ICOFFA = MOD( JA-1, DESCA( NB_ ) )
            IROFFB = MOD( IB-1, DESCB( MB_ ) )
            IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.  &
     &         LSAME( TRANS, 'C' ) ) THEN
               INFO = -1
            ELSE IF( IROFFA.NE.0 ) THEN
               INFO = -5
            ELSE IF( ICOFFA.NE.0 ) THEN
               INFO = -6
            ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(700+NB_)
!            ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
!               INFO = -10
            ELSE IF( DESCB( MB_ ).NE.DESCA( NB_ ) ) THEN
               INFO = -(1200+NB_)
            ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
               INFO = -(1200+CTXT_)
            END IF
         END IF
         IF( NOTRAN ) THEN
            IDUM1( 1 ) = ICHAR( 'N' )
         ELSE IF( LSAME( TRANS, 'T' ) ) THEN
            IDUM1( 1 ) = ICHAR( 'T' )
         ELSE
            IDUM1( 1 ) = ICHAR( 'C' )
         END IF
         IDUM2( 1 ) = 1
         CALL PCHK2MAT( N, 2, N, 2, IA, JA, DESCA, 7, N, 2, NRHS, 3,    &
     &                  IB, JB, DESCB, 12, 1, IDUM1, IDUM2, INFO )
      END IF

      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PDGETRS', -INFO )
         RETURN
      END IF
!
!     Quick return if possible
!
      IF( N.EQ.0 .OR. NRHS.EQ.0 ) RETURN

      CALL DESCSET( DESCIP, DESCA( M_ ) + DESCA( MB_ )*NPROW, 1,        &
     &              DESCA( MB_ ), 1, DESCA( RSRC_ ), MYCOL, ICTXT,      &
     &              DESCA( MB_ ) + NUMROC( DESCA( M_ ), DESCA( MB_ ),   &
     &              MYROW, DESCA( RSRC_ ), NPROW ) )

      IF( NOTRAN ) THEN
!
!        Solve sub( A ) * X = sub( B ).
!
!        Apply row interchanges to the right hand sides.
!
         CALL PDLAPIV( 'Forward', 'Row', 'Col', N, NRHS, B, IB, JB,     &
     &                 DESCB, IPIV, IA, 1, DESCIP, IDUM1 )
!
!        Solve L*X = sub( B ), overwriting sub( B ) with X.
!
         CALL PDTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, &
     &                ONE, A, IA, JA, DESCA, B, IB, JB, DESCB )
!
!        Solve U*X = sub( B ), overwriting sub( B ) with X.
!
         CALL PDTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,   &
     &                NRHS, ONE, A, IA, JA, DESCA, B, IB, JB, DESCB )
      ELSE
!
!        Solve sub( A )' * X = sub( B ).
!
!        Solve U'*X = sub( B ), overwriting sub( B ) with X.
!
         CALL PDTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, &
     &                ONE, A, IA, JA, DESCA, B, IB, JB, DESCB )
!
!        Solve L'*X = sub( B ), overwriting sub( B ) with X.
!
         CALL PDTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS,    &
     &                ONE, A, IA, JA, DESCA, B, IB, JB, DESCB )
!
!        Apply row interchanges to the solution vectors.
!
         CALL PDLAPIV( 'Backward', 'Row', 'Col', N, NRHS, B, IB, JB,    &
     &                 DESCB, IPIV, IA, 1, DESCIP, IDUM1 )

      END IF

      RETURN

!     End of PDGETRS

      END SUBROUTINE PDGETRS

      SUBROUTINE PDELGET( SCOPE, TOP, ALPHA, A, IA, JA, DESCA )
      use descriptor_mod, desca_x=>desca
      implicit none
!
!  -- ScaLAPACK tools routine (version 1.7) --
!     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
!     and University of California, Berkeley.
!     May 1, 1997
!
!     .. Scalar Arguments ..
      CHARACTER*1        SCOPE, TOP
      INTEGER            IA, JA
      REAL(dp)   ALPHA
!     ..
!     .. Array arguments ..
      INTEGER            DESCA( * )
      REAL(dp)   A( * )
!     ..
!
!  Purpose
!  =======
!
!  PDELGET sets alpha to the distributed matrix entry A( IA, JA ).
!  The value of alpha is set according to the scope.
!
!  Notes
!  =====
!
!  Each global data object is described by an associated description
!  vector.  This vector stores the information required to establish
!  the mapping between an object element and its corresponding process
!  and memory location.
!
!  Let A be a generic term for any 2D block cyclicly distributed array.
!  Such a global array has an associated description vector DESCA.
!  In the following comments, the character _ should be read as
!  "of the global array".
!
!  NOTATION        STORED IN      EXPLANATION
!  --------------- -------------- --------------------------------------
!  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
!                                 DTYPE_A = 1.
!  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
!                                 the BLACS process grid A is distribu-
!                                 ted over. The context itself is glo-
!                                 bal, but the handle (the integer
!                                 value) may vary.
!  M_A    (global) DESCA( M_ )    The number of rows in the global
!                                 array A.
!  N_A    (global) DESCA( N_ )    The number of columns in the global
!                                 array A.
!  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
!                                 the rows of the array.
!  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
!                                 the columns of the array.
!  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
!                                 row of the array A is distributed.
!  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
!                                 first column of the array A is
!                                 distributed.
!  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
!                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
!
!  Let K be the number of rows or columns of a distributed matrix,
!  and assume that its process grid has dimension p x q.
!  LOCr( K ) denotes the number of elements of K that a process
!  would receive if K were distributed over the p processes of its
!  process column.
!  Similarly, LOCc( K ) denotes the number of elements of K that a
!  process would receive if K were distributed over the q processes of
!  its process row.
!  The values of LOCr() and LOCc() may be determined via a call to the
!  ScaLAPACK tool function, NUMROC:
!          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
!          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
!  An upper bound for these quantities may be computed by:
!          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )!MB_A
!          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )!NB_A
!
!  Arguments
!  =========
!
!  SCOPE   (global input) CHARACTER*1
!          The BLACS scope in which alpha is updated.
!          If SCOPE = 'R', alpha is updated only in the process row
!                          containing A( IA, JA ),
!          If SCOPE = 'C', alpha is updated only in the process column
!                          containing A( IA, JA ),
!          If SCOPE = 'A', alpha is updated in all the processes of the
!                          grid,
!          otherwise alpha is updated only in the process containing
!           A( IA, JA ).
!
!  TOP     (global input) CHARACTER!1
!          The topology to be used if broadcast is needed.
!
!  ALPHA   (global output) DOUBLE PRECISION, the scalar alpha.
!
!  A       (local input) REAL(dp) pointer into the local memory
!          to an array of dimension (LLD_A,!) containing the local
!          pieces of the distributed matrix A.
!
!  IA      (global input) INTEGER
!          The row index in the global array A indicating the first
!          row of sub( A ).
!
!  JA      (global input) INTEGER
!          The column index in the global array A indicating the
!          first column of sub( A ).
!
!  DESCA   (global and local input) INTEGER array of dimension DLEN_.
!          The array descriptor for the distributed matrix A.
!
!  =====================================================================
!
!     .. Parameters ..
!      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
!     $                   LLD_, MB_, M_, NB_, N_, RSRC_
!      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
!     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
!     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL(dp)   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
!     ..
!     .. Local Scalars ..
      INTEGER            IACOL, IAROW, ICTXT, IIA, IOFFA, JJA  !, MYCOL, MYROW, NPCOL, NPROW
!     ..
!     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DGEBR2D, DGEBS2D, INFOG2L
!     ..
!     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
!     ..
!     .. Executable Statements ..
!
!     Get grid parameters.
!
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
!
      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,&
                    IAROW, IACOL )
!
      ALPHA = ZERO
!
      IF( LSAME( SCOPE, 'R' ) ) THEN
         IF( MYROW.EQ.IAROW ) THEN
            IF( MYCOL.EQ.IACOL ) THEN
               IOFFA = IIA+(JJA-1)*DESCA( LLD_ )
               CALL DGEBS2D( ICTXT, SCOPE, TOP, 1, 1, A( IOFFA ), 1 )
               ALPHA = A( IOFFA )
            ELSE
               CALL DGEBR2D( ICTXT, SCOPE, TOP, 1, 1, ALPHA, 1,         &
                             IAROW, IACOL )
            END IF
         END IF
      ELSE IF( LSAME( SCOPE, 'C' ) ) THEN
         IF( MYCOL.EQ.IACOL ) THEN
            IF( MYROW.EQ.IAROW ) THEN
               IOFFA = IIA+(JJA-1)*DESCA( LLD_ )
               CALL DGEBS2D( ICTXT, SCOPE, TOP, 1, 1, A( IOFFA ), 1 )
               ALPHA = A( IOFFA )
            ELSE
               CALL DGEBR2D( ICTXT, SCOPE, TOP, 1, 1, ALPHA, 1,         &
                             IAROW, IACOL )
            END IF
         END IF
      ELSE IF( LSAME( SCOPE, 'A' ) ) THEN
         IF( ( MYROW.EQ.IAROW ).AND.( MYCOL.EQ.IACOL ) ) THEN
            IOFFA = IIA+(JJA-1)*DESCA( LLD_ )
            CALL DGEBS2D( ICTXT, SCOPE, TOP, 1, 1, A( IOFFA ), 1 )
            ALPHA = A( IOFFA )
         ELSE
            CALL DGEBR2D( ICTXT, SCOPE, TOP, 1, 1, ALPHA, 1,            &
                          IAROW, IACOL )
         END IF
      ELSE
         IF( MYROW.EQ.IAROW .AND. MYCOL.EQ.IACOL )                      &
            ALPHA = A( IIA+(JJA-1)*DESCA( LLD_ ) )
      END IF
!
      RETURN
!
!     End of PDELGET
!
      END

      SUBROUTINE PDELSET( A, IA, JA, DESCA, ALPHA )
      use descriptor_mod, desca_x=>desca
      implicit none
!
!  -- ScaLAPACK tools routine (version 1.7) --
!     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
!     and University of California, Berkeley.
!     May 1, 1997
!
!     .. Scalar Arguments ..
      INTEGER            IA, JA
      REAL(dp)   ALPHA
!     ..
!     .. Array arguments ..
      INTEGER            DESCA( * )
      REAL(dp)   A( * )
!     ..
!
!  Purpose
!  =======
!
!  PDELSET sets the distributed matrix entry A( IA, JA ) to ALPHA.
!
!  Notes
!  =====
!
!  Each global data object is described by an associated description
!  vector.  This vector stores the information required to establish
!  the mapping between an object element and its corresponding process
!  and memory location.
!
!  Let A be a generic term for any 2D block cyclicly distributed array.
!  Such a global array has an associated description vector DESCA.
!  In the following comments, the character _ should be read as
!  "of the global array".
!
!  NOTATION        STORED IN      EXPLANATION
!  --------------- -------------- --------------------------------------
!  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
!                                 DTYPE_A = 1.
!  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
!                                 the BLACS process grid A is distribu-
!                                 ted over. The context itself is glo-
!                                 bal, but the handle (the integer
!                                 value) may vary.
!  M_A    (global) DESCA( M_ )    The number of rows in the global
!                                 array A.
!  N_A    (global) DESCA( N_ )    The number of columns in the global
!                                 array A.
!  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
!                                 the rows of the array.
!  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
!                                 the columns of the array.
!  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
!                                 row of the array A is distributed.
!  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
!                                 first column of the array A is
!                                 distributed.
!  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
!                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
!
!  Let K be the number of rows or columns of a distributed matrix,
!  and assume that its process grid has dimension p x q.
!  LOCr( K ) denotes the number of elements of K that a process
!  would receive if K were distributed over the p processes of its
!  process column.
!  Similarly, LOCc( K ) denotes the number of elements of K that a
!  process would receive if K were distributed over the q processes of
!  its process row.
!  The values of LOCr() and LOCc() may be determined via a call to the
!  ScaLAPACK tool function, NUMROC:
!          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
!          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
!  An upper bound for these quantities may be computed by:
!          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
!          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
!
!  Arguments
!  =========
!
!  A       (local output) REAL(dp) pointer into the local memory
!          to an array of dimension (LLD_A,*) containing the local
!          pieces of the distributed matrix A.
!
!  IA      (global input) INTEGER
!          The row index in the global array A indicating the first
!          row of sub( A ).
!
!  JA      (global input) INTEGER
!          The column index in the global array A indicating the
!          first column of sub( A ).
!
!  DESCA   (global and local input) INTEGER array of dimension DLEN_.
!          The array descriptor for the distributed matrix A.
!
!  ALPHA   (local input) DOUBLE PRECISION
!          The scalar alpha.
!
!  =====================================================================
!
!     .. Parameters ..
!      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
!     $                   LLD_, MB_, M_, NB_, N_, RSRC_
!      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
!     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
!     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
!     ..
!     .. Local Scalars ..
      INTEGER            IACOL, IAROW, IIA, JJA   !, MYCOL, MYROW, NPCOL,  NPROW
!     ..
!     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, INFOG2L
!     ..
!     .. Executable Statements ..
!
!     Get grid parameters.
!
      CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
!
      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,&
                    IAROW, IACOL )
!
      IF( MYROW.EQ.IAROW .AND. MYCOL.EQ.IACOL )                         &
         A( IIA+(JJA-1)*DESCA( LLD_ ) ) = ALPHA
!
      RETURN
!
!     End of PDELSET
!
      END
#endif

#endif

	  end module ptrd_mod

max_stars_repo_name: ORNL-Fusion/SIESTA

max_issues_repo_name: ORNL-Fusion/SIESTA

max_forks_repo_name: ORNL-Fusion/SIESTA


text:
!>
!! \brief This module sets elemental abundances
!!
!! Module for Capreole / C2-Ray (f90)
!!
!! \b Author: Garrelt Mellema
!!
!! \b Date: 2008-05-27
!!
!! This module is also accepted by the F compiler (Dec 9, 2003)
!<

module abundances

  use precision, only: dp

  implicit none

  private

  !> Helium abundance (by number)
  real(kind=dp),public,parameter :: abu_he=0.1_dp
  !> Carbon abundance (by number)
  real(kind=dp),public,parameter :: abu_c=3.3e-4_dp
  !> Oxygen abundance (by number)
  real(kind=dp),parameter,public :: abu_o=3.0e-4_dp

  !> Mean molecular weight (assuming metallicities of order 10^-4).
  real(kind=dp),public,parameter :: mu=(1.0-abu_he)+4.0*abu_he

  !-----------------------------------------------------------------------
  !     set heavy element abundances
  !     from Osterbrock (1989), Table 5.13
  !-----------------------------------------------------------------------
  !parameter(abu_he=0.1)
  !abu_he=0.074              !cosmic abundance of He (24.2% by mass)
  !parameter(abu_c=7.1e-7)
  !abu_n=1.3e-7
  !abu_o=4.2e-7
  !abu_ne=1.0e-7
  !abu_s=1.0e-8

  !abu_c=3.3e-4
  !   abu_n=0.0
  !   abu_o=6.6e-4
  !   abu_ne=8.3e-5
  !   abu_s=1.6e-5

end module abundances

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!!*****************************************************
!
!         Turbo AMG (ramg) interface functions
!         1. V_cycle
!         2. plain
!         3. driver
!
!*****************************************************

      !*******************************************
      !    ramg_V_cycle:
      !     One V_cycle call.
      !     Now use SAMG, later will use own code
      !*******************************************
      recursive subroutine ramg_V_cycle(
     &               ramg_sol,ramg_rhs,clevel,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper
     &           )!ramg_V_cycle
      use ramg_data
      include "common.h"
      include "mpif.h"

      !Variable Declaration
      ! arguments
      integer, intent(in)         :: clevel
      real(kind=8),intent(inout),dimension(amg_nshg(clevel))
     &               :: ramg_sol,ramg_rhs

      integer,intent(in),dimension(nshg+1)       :: colm
      integer,intent(in),dimension(nnz_tot)      :: rowp
      real(kind=8),intent(in),dimension(9,nnz_tot)  :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)  :: lhsP
      integer,intent(in),dimension(nlwork)          :: ilwork
      integer,intent(in),dimension(nshg)            :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC

      ! local
      real(kind=8),dimension(amg_nshg(clevel))  :: myvF,myvE
      real(kind=8),dimension(:),allocatable :: myvC,myvCS
      real(kind=8)          :: cpusec(10)

      ! In scale

      if (clevel.ne.1) then
      ramg_rhs = ramg_rhs*amg_ppeDiag(clevel)%p
      endif

      if (ramg_levelx-clevel .eq. 0) then
      ! finest level solver
          call ramg_direct_solve(ramg_rhs,ramg_sol,
     &    colm,rowp,lhsK,lhsP,
     &    ilwork,BC,iBC,iper)
      !**********************
      ! higher level smoother
      !**********************
      else if (ramg_levelx-clevel .gt. 0) then
          allocate(myvC(amg_nshg(clevel+1)))
          allocate(myvCS(amg_nshg(clevel+1)))

          ! smoothing (pre)
          myvF = 0 ! initial guess is zero
          call ramg_smoother(clevel,ramg_sol,myvF,ramg_rhs,
     &              colm,rowp,lhsK,lhsP,
     &              ilwork,BC,iBC,iper,2) !dx1

          ! restriction
          call ramg_calcAv_g(clevel,myvF,ramg_sol,colm,rowp,lhsK,lhsP,
     &               ilwork,BC,iBC,iper,1)
          myvF =  myvF - ramg_rhs! r2
          call ramg_calcIvFC(clevel,clevel+1,myvF,myvC)
          ! up one level, solve
          myvCS = 0
          call ramg_V_cycle(myvCS,myvC,clevel+1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
          ! prolongation
          call ramg_calcIvCF(clevel,clevel+1,myvCS,myvF,
     &         ilwork,BC,iBC,iper) ! v2hat
          ramg_sol = ramg_sol - myvF

          ! smoothing (post)
          call ramg_smoother(clevel,myvF,ramg_sol,ramg_rhs,
     &              colm,rowp,lhsK,lhsP,
     &              ilwork,BC,iBC,iper,1)
          ramg_sol =  myvF

          deallocate(myVC)
          deallocate(myVCS)
      end if

      ! Out scale
      if (clevel.ne.1) then
      ramg_sol = ramg_sol * amg_ppeDiag(clevel)%p
      endif

      return

      end subroutine ramg_V_cycle

!*****************************************************
!      ramg Cg solver
!      using ramg_V_cycle as preconditioner
!*****************************************************
      subroutine ramg_CG(
     &              sol,rhs,level,
     &              colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper,iterNum)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer,intent(in) :: level
      real(kind=8),intent(in),dimension(amg_nshg(level)) :: rhs
      real(kind=8),intent(inout),dimension(amg_nshg(level)) :: sol

      integer,intent(in),dimension(nshg+1)       :: colm
      integer,intent(in),dimension(nnz_tot)      :: rowp

      real(kind=8),intent(in),dimension(9,nnz_tot)    :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)    :: lhsP

      integer,intent(in),dimension(nlwork)   :: ilwork
      integer,intent(in),dimension(nshg)     :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC
      integer,intent(inout)                    :: iterNum

      ! local
      real(kind=8),dimension(amg_nshg(level))   :: cgP,cgQ,cgZ,cgR
      real(kind=8)                       :: rz,rz_0,pq,alpha,beta
      real(kind=8)               :: tmp,norm_0,norm_1,norm_e,norm_c
      real(kind=8)               :: mres_0,mres_n

      cgR = rhs
      call ramg_L2_norm_p(level,cgR,1,norm_0)
      norm_c = norm_0
      norm_e = norm_0*1e-7

      cgZ = cgR

      call ramg_dot_p(level,cgZ,cgR,1,rz)
      rz_0 = rz
      cgP = cgZ
      sol = 0

      iterNum = 1

      do while ( (norm_c .gt. norm_e).and.(iterNum.lt.500) )
         call ramg_calcAv_g(level,cgQ,cgP,colm,rowp,lhsK,lhsP,
     &                      ilwork,BC,iBC,iper,1)
         call ramg_dot_p(level,cgP,cgQ,1,pq)
         alpha = rz/pq
         sol = sol + alpha*cgP
         cgR = cgR - alpha*cgQ
         call ramg_L2_norm_p(level,cgR,1,tmp)
         norm_c = tmp
         cgZ = cgR
         tmp = rz
         call ramg_dot_p(level,cgZ,cgR,1,rz)
         beta = rz/tmp
         cgP = beta*cgP+cgZ
         iterNum = iterNum + 1
      enddo

      end subroutine ! ramg_CG

      subroutine ramg_PCG(
     &              sol,rhs,run_mode,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper
     &           )!ramg_CG
      use ramg_data
      include "common.h"
      include "mpif.h"

      real(kind=8),intent(inout),dimension(nshg) :: sol
      real(kind=8),intent(in),dimension(nshg)    :: rhs
      integer,intent(in)                                :: run_mode

      integer,intent(in),dimension(nshg+1)       :: colm
      integer,intent(in),dimension(nnz_tot)      :: rowp
      real(kind=8),intent(in),dimension(9,nnz_tot)  :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)  :: lhsP
      integer,intent(in),dimension(nlwork)          :: ilwork
      integer,intent(in),dimension(nshg)            :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC

      ! local
      real(kind=8),dimension(amg_nshg(1))   :: cgP,cgQ,cgZ,cgR
      real(kind=8)                       :: rz,rz_0,pq,alpha,beta
      real(kind=8)               :: tmp,norm_0,norm_p,norm_e,norm_c
      real(kind=8)               :: mres_0,mres_n
      real(kind=8),dimension(nshg,4)               :: diag
      integer                    :: iterNum

      ! controls

      real(kind=8),dimension(maxiters) :: rconv
      real(kind=8) :: avgconv
      logical :: vcycle,gcycle,restart
      character cflagv,cflagg

      cgR = rhs
      call ramg_L2_norm_p(1,cgR,1,norm_0)
      norm_c = norm_0
      norm_e = norm_0 * epstol(2)
      write(*,*)'norm_0 ', norm_0

      vcycle = .true.
      gcycle = .true.
      restart = .false.

      ! amg preconditioning control

      if (vcycle) then
         call ramg_V_cycle(cgZ,cgR,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         if (gcycle) then
         call ramg_G_cycle(cgZ,cgR,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         endif
      else
          cgZ = cgR
      endif

      call ramg_dot_p(1,cgZ,cgR,1,rz)
      rz_0 = rz

      cgP = cgZ
      sol = 0

      iterNum = 1
      rconv(1) = 1.0
      norm_p=norm_c
      avgconv = 0

      do while ( (norm_c.gt.norm_e).and.(iterNum.lt.maxiters))
      !do while ((iterNum.lt.maxiters))
         call ramg_calcAv_g(1,cgQ,cgP,colm,rowp,lhsK,lhsP,
     &                      ilwork,BC,iBC,iper,1)
         call ramg_dot_p(1,cgP,cgQ,1,pq)
         alpha = rz/pq
         sol = sol + alpha*cgP
         cgR = cgR - alpha*cgQ
         call ramg_L2_norm_p(1,cgR,1,tmp)
         norm_c = tmp
         if (iterNum.eq.1) rconv(1) = norm_p/norm_c

      ! amg preconditioning control
      if (vcycle) then
         call ramg_V_cycle(cgZ,cgR,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         if (gcycle) then
         call ramg_G_cycle(cgZ,cgR,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         endif
      else
          cgZ = cgR
      endif

         tmp = rz
         call ramg_dot_p(1,cgZ,cgR,1,rz)
         beta = rz/tmp
         cgP = beta*cgP + cgZ
         iterNum = iterNum + 1
         rconv(iterNum) = norm_p/norm_c
         norm_p = norm_c
         cflagv='n'
         cflagg='n'
         if (vcycle) cflagv='v'
         if (gcycle) cflagg='g'
         if (myrank.eq.master) then
         write(*,"((A7)(I4)(E14.4)(T30)(F8.4)(T40)(A1)(A1))")
     &    'AMGCG: ',
     &    iterNum,norm_c/norm_0,rconv(iterNum),cflagv,cflagg
         end if
      enddo

      end subroutine !ramg_PCG

!*******************************************
!      ramg Interface
!      Interface for libLES
!*******************************************
      subroutine ramg_interface(
     &           colm,rowp,lhsK,lhsP,flowDiag,
     &           mcgR,mcgZ,
     &           ilwork,BC,iBC,iper
     &           )
      use ramg_data
      include "common.h"
      include "mpif.h"
      integer,intent(in),dimension(nshg+1)             :: colm
      integer,intent(in),dimension(nnz_tot)            :: rowp
      real(kind=8),intent(in),dimension(9,nnz_tot)     :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)     :: lhsP
      real(kind=8),intent(in),dimension(nshg)          :: mcgR
      real(kind=8),intent(inout),dimension(nshg)       :: mcgZ
      integer, intent(in), dimension(nlwork)           :: ilwork
      integer, intent(in),dimension(nshg)              :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC)   :: BC
      real(kind=8),dimension(nshg)                     :: myR
      real(kind=8),intent(in),dimension(nshg,4)        :: flowDiag

      integer::i,j,k,amgmode
      logical :: precflag
      !character*10                                     :: fname

      mcgZ = 0
      amgmode = 11

      ramg_winct = mod(ramg_winct+1,2)

      if (irun_amg_prec.eq.0) then
          precflag = .false.
      else if ((irun_amg_prec.ge.2).and.(ramg_flag.eq.1)) then
      ! first solve attempt with CG if using smart solver
          precflag=.false.
      else
          precflag = .true.
      endif

      if (precflag) then !.and.(ramg_redo.ne.0)) then
          myR = mcgR*amg_ppeDiag(1)%p
          call ramg_driver(colm,rowp,lhsK,lhsP,
     &           nshg,nnz_tot,nflow,
     &           myR,
     &           nlwork,ilwork,ndofBC,BC,iBC,iper,mcgZ,1,
     &           ramg_eps,amgmode)
          mcgZ = mcgZ*amg_ppeDiag(1)%p
      else
          mcgZ = mcgR
      endif

      end subroutine ramg_interface

!*******************************************
!        ramg Driver
!        1. extract PPE / system
!        2. call coarsening
!        3. solve
!********************************************

      !*************************************
      !    ramg_driver
      !     Input:  global matrix,# of systems
      !             control params
      !     Output: solution
      !*************************************
      subroutine ramg_driver(
     &           colm,rowp,lhsK,lhsP,
     &           nshg,nnz_tot,nflow,
     &           pperhs,
     &           nlwork,ilwork,ndofBC,BC,iBC,iper,
     &           ramg_sol,n_sol,
     &           amg_eps,amg_mode
     &           )

      use ramg_data
      implicit none

      !***********parameters**************
      !the matrix
      integer,intent(in)                               :: nshg
      integer,intent(in)                        :: nnz_tot
      integer,intent(in)                               :: nflow
      !the matrix
      integer,intent(in),dimension(nshg+1)             :: colm
      integer,intent(in),dimension(nnz_tot)            :: rowp
      real(kind=8),intent(in),dimension(9,nnz_tot)     :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)     :: lhsP
      !the forcing term, rhs
      real(kind=8),intent(in),dimension(nshg)       :: pperhs
      ! the boundary info
      integer, intent(in)                              :: nlwork
      integer, intent(in), dimension(nlwork)           :: ilwork
      integer, intent(in)                              :: ndofBC
      integer, intent(in),dimension(nshg)              :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC)   :: BC
      !the output solution
      integer, intent(in)                              :: n_sol
      real(kind=8),intent(inout),dimension(nshg,n_sol) :: ramg_sol

      !the tolerance
      real(kind=8),intent(in)                          :: amg_eps
      !control run mode
      integer,intent(inout)                            :: amg_mode
      ! my AMG parameter
      !*********parameters end**************

      !*****local variables*****************
      real(kind=8)                      :: cpusec(10)

      call cpu_time(cpusec(1))

      if (amg_mode .eq. 11 ) then     ! solve PPE les Precondition ramg
         call cpu_time(cpusec(2))
         call ramg_V_cycle(ramg_sol,pperhs,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         call ramg_G_cycle(ramg_sol,pperhs,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
         call cpu_time(cpusec(3))
         ramg_time(4)=ramg_time(4)+cpusec(3)-cpusec(2)
      else if (amg_mode .eq. 3) then ! solve PPE CG
          call cpu_time(cpusec(2))
          call ramg_PCG(ramg_sol,pperhs,1,
     &      colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
          call cpu_time(cpusec(3))
          ramg_time(4)=ramg_time(4)+cpusec(3)-cpusec(2)
      end if

      return

      end subroutine ramg_driver
!*******************************************
!      <EOF> ramg Driver
!*******************************************

ext: f


text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE INTSTR (NNUM, LNGSTR, INUM, ISTR, LSTR)
C=======================================================================
C$Id: intstr.f,v 1.2 2009/03/25 12:46:02 gdsjaar Exp $
C$Log: intstr.f,v $
CRevision 1.2  2009/03/25 12:46:02  gdsjaar
CAdd copyright and license notice to all files.
C
CRevision 1.1.1.1  1990/08/14 16:15:14  gdsjaar
CTesting
C
c Revision 1.1  90/08/14  16:15:12  gdsjaar
c Initial revision
c
c Revision 1.1  90/08/09  13:39:32  gdsjaar
c Initial revision
c

C   --*** INTSTR *** (STRLIB) Convert integer numbers to strings
C   --   Written by Amy Gilkey - revised 03/14/88
C   --
C   --INTSTR converts a set of integer numbers into a consistent set of
C   --strings of the same length (right justified).
C   --
C   --Parameters:
C   --   NNUM - IN - the number of integers numbers in the set
C   --   LNGSTR - IN - the number of digits in the string; <= 0 if length
C   --      of longest integer
C   --   INUM - IN - the array of integers to be converted
C   --   ISTR - OUT - the set of integer strings
C   --   LSTR - OUT - the maximum length of the number strings

      INTEGER NNUM
      INTEGER LNGSTR
      INTEGER INUM(*)
      CHARACTER*(*) ISTR(*)
      INTEGER LSTR

      CHARACTER*20 SCISTR
      CHARACTER*5 FFMT

      IF ((NNUM .EQ. 1) .AND. (LNGSTR .LE. 0)) THEN

C      --Handle special case of single number, smallest length

         WRITE (FFMT, 10000, IOSTAT=IDUM) LEN (ISTR(1))

         WRITE (ISTR(1), FFMT, IOSTAT=IDUM) INUM(1)
         CALL SQZSTR (ISTR(1), LSTR)

      ELSE

C      --Determine smallest length of largest number string

         IF (LNGSTR .LE. 0) THEN
            MINE = 0
            MAXE = 0
            DO 100 I = 1, NNUM
               MINE = MIN (MINE, INUM(I))
               MAXE = MAX (MAXE, INUM(I))
  100       CONTINUE

            MM = MAX (IABS (MAXE), IABS (MINE))
            IF (MM .EQ. IABS (MINE)) MM = -MM
            WRITE (SCISTR, '(I12)', IOSTAT=IDUM) MM
            CALL SQZSTR (SCISTR, LSTR)

         ELSE
            LSTR = LNGSTR
         END IF

         LSTR = MIN (LSTR, LEN (ISTR(1)))

C      --Convert all integers to string of given length

         WRITE (FFMT, 10000, IOSTAT=IDUM) LSTR
10000     FORMAT ('(I', I2.2, ')')

         DO 110 I = 1, NNUM
            WRITE (ISTR(I), FFMT, IOSTAT=IDUM) INUM(I)
  110    CONTINUE
      END IF

      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text: #redirect Buddhist Meditation Group at the Unitarian Church

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
James Cothern served as Trustee of the Davis JUSD School Board from 1977 to 1981.  At the time Cothern worked as an economist for the UC Davis agricultural extension service.

Source:

Sias, Spencer; 9 March 1977, Cothern, Damesyn, Darling, School critics triumph,  Davis Enterprise, pg. 1.


text:
module test2008_15_m2
  integer i
  contains
    subroutine ffff(i,pppp)
      integer i
      integer pppp(i)
       integer x
       x = i+2
    end subroutine ffff
end module test2008_15_m2

subroutine u1
 use test2008_15_m2
 integer,parameter::cc=12
 integer,allocatable::a(:)
 allocate(a(cc))
 call ffff(cc,a)
end subroutine u1

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


SUBROUTINE CmpXqtl(xqtl2, DMFTrans, nbands, nind, maxdim2)
  ! Creates xqtl from transformation DMFTrans
  IMPLICIT NONE
  COMPLEX*16, intent(inout):: xqtl2(nbands,maxdim2,maxdim2)
  COMPLEX*16, intent(in)   :: DMFTrans(nbands,nbands,maxdim2,maxdim2)
  INTEGER, intent(in)      :: nbands, nind, maxdim2
  ! local variables
  INTEGER    :: ind1, ind2, i
  do ind1=1,nind
     do ind2=1,nind
        do i=1,nbands ! over bands
           xqtl2(i,ind1,ind2) = xqtl2(i,ind1,ind2) + DMFTrans(i,i,ind1,ind2)
        enddo
     enddo
  enddo
END SUBROUTINE CmpXqtl


text: The Yolo Plankroad was a raised road of wooden planks that spanned a distance of four and a half miles across the Yolo Basin.  Built in 1855, it extended from the Sacramento River to a portion of the Davis Wetlands located about five miles to the southeast of Woodland.  The western endpoint of the road was marked by the Tule House, a hotel and stage stop built on stilts.  The main advantage of the unique road was that it could be used by travelers in the winter and spring, seasons which normally made the area impassable due to flooding.  However, the road required a high degree of upkeep, and, despite its purpose, was susceptible to flooding during some winter months.  The Yolo Plankroad was eventually replaced by the Tule Jake Road, a regular road which was only travelable in the dry season.


text:
      SUBROUTINE SUBA
C
C     UNSTEADY FLOW ANAYSIS OF A SUPERSONIC CASCADE
C
C     LIFT AND MOMENT COEFICIENT
C
      DIMENSION       PRES1(21),PRES2(21),PRES3(21),PRES4(21),QRES4(21),
     1                SBKDE1(201),SBKDE2(201),SUMSV1(201),SUMSV2(201),
     2                SVKL1(201),SVKL2(201),XLSV1(21),XLSV2(21),
     3                XLSV3(21),XLSV4(21)
      COMPLEX         SBKDE1,SBKDE2,F4,F4S,AM4,F5S,F6S,AM4TST,SUM3,SUM4,
     1                AM5TT,AM6,SUMSV1,SUMSV2,SVKL1,SVKL2,F5,F5T,AM5,
     2                AM5T,AI,A,B,BSYCON,ALP,F1,AM1,ALN,BLKAPM,BKDEL3,
     3                F1S,C1,C2P,C2N,C2,AMTEST,FT2,BLAM1,FT3,AM2,SUM1,
     4                SUM2,F2,BLAM2,FT2T,C1T,FT3T,F2P,AM2P,SUM1T,SUM2T,
     5                C1P,C1N,BKDEL1,BKDEL2,BLKAP1,ARG,ARG2,FT3TST,BC,
     6                BC2,BC3,BC4,BC5,CA1,CA2,CA3,CA4,CLIFT,CMOMT,
     7                PRES1,PRES2,PRES3,PRES4,QRES4,FQA,FQB,T1,T2,T3,T4,
     8                GUSAMP,FQ7,CEXP3,CEXP4,CEXP5,CONST,C1A,C2A
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /SYSTEM/ SYSBUF,IBBOUT
      COMMON /BLK1  / SCRK,SPS,SNS,DSTR,AI,PI,DEL,SIGMA,BETA,RES
      COMMON /BLK2  / BSYCON
      COMMON /BLK3  / SBKDE1,SBKDE2,F4,F4S,AM4,F5S,F6S,AM4TST,SUM3,SUM4,
     1                AM5TT,AM6,SUMSV1,SUMSV2,SVKL1,SVKL2,F5,F5T,AM5,
     2                AM5T,A,B,ALP,F1,AM1,ALN,BLKAPM,BKDEL3,F1S,C1,C2P,
     3                C2N,C2,AMTEST,FT2,BLAM1,FT3,AM2,SUM1,SUM2,F2,
     4                BLAM2,FT2T,C1T,FT3T,F2P,AM2P,SUM1T,SUM2T,C1P,C1N,
     5                BKDEL1,BKDEL2,BLKAP1,ARG,ARG2,FT3TST,BC,BC2,BC3,
     6                BC4,BC5,CA1,CA2,CA3,CA4,CLIFT,CMOMT,PRES1,PRES2,
     7                PRES3,PRES4,QRES4,FQA,FQB,FQ7
      COMMON /BLK4  / I,R,Y,A1,B1,C4,C5,GL,I6,I7,JL,NL,RI,RT,R5,SN,SP,
     1                XL,Y1,AMU,GAM,IDX,INX,NL2,RL1,RL2,RQ1,RQ2,XL1,
     2                ALP1,ALP2,GAMN,GAMP,INER,IOUT,REDF,STAG,STEP,
     3                AMACH,BETNN,BETNP,BKAP1,XLSV1,XLSV2,XLSV3,XLSV4,
     4                ALPAMP,AMOAXS,GUSAMP,DISAMP,PITAXS,PITCOR
C
      S1   = SPS - SNS
      S2   = SPS*DEL - SIGMA
      S3   = SPS/(DSTR**2)
      S4   = SNS/DSTR
      S0   = 2.0 - SPS + SNS
      T1   = CEXP(-AI*SIGMA)
      T2   = CEXP(AI*SIGMA)
      A1   = 2.0*PI/S1
      B1   = S2/S1
      GAM  = S2
      C1P  = GAM/DSTR - SCRK
      C1N  = GAM/DSTR + SCRK
      ALP  = GAM*S3 + S4*CSQRT(C1P)*CSQRT(C1N)
      BC   = -B1/ALP*BSYCON/SIN(PI*B1/A1)
      T3   = ALP - DEL
      F1   = (ALP-AMU)/T3*AI*SNS/(BETA*(GAM-ALP*SPS))
      ARG2 = DEL
      CALL AKAPM (ARG2,BKDEL1)
      ARG  = DEL - GL
      CALL AKAPM  (ARG,BKDEL2)
      CALL DLKAPM (ARG2,BLKAP1)
      INX  = 0
      CALL DRKAPM (ALP,INX,BLKAPM)
      F1   = F1*BKDEL1/BLKAPM*(-T3/(T3+GL)*A*AI*BKDEL2/BKDEL1 +
     1       B*BLKAP1+B/T3)
      F1S  = F1
      NL   = 10
      RL1  = NL - 1
      CEXP3 = CEXP(-AI*T3/RL1*S1)
      PRES1(1) = F1S
      NNL1 = NL - 1
      DO 453 JL = 1,NNL1
      PRES1(JL+1) = PRES1(JL)*CEXP3
  453 CONTINUE
      F1   = F1*AI/T3*(CEXP(-AI*T3*S1)-1.0)
      AM1  = F1/(AI*T3)-F1S/(AI*T3)*S1*CEXP(-AI*T3*S1)
      AMTEST = 0.0
      FQB  = BKDEL1/(BETA*BC)*CEXP(AI*S2/2.0)*
     1      (-A*AI*BKDEL2/BKDEL1+B*BLKAP1)
      DO 20 I = 1,200
      R    = I
      GAMP = 2.0*PI*R + S2
      GAMN =-2.0*PI*R + S2
      C1P  = (GAMP/DSTR) - SCRK
      C2P  = (GAMP/DSTR) + SCRK
      ALP  = GAMP*S3 + S4*CSQRT(C1P)*CSQRT(C2P)
      T3   = ALP - DEL
      IDX  = I
      CALL DRKAPM (ALP,IDX,BLKAPM)
      C1   = (ALP-AMU)/T3*AI*SNS/(BETA*(GAMP-ALP*SPS))*BKDEL1/
     1       (BLKAPM)*(-T3/(T3+GL)*A*AI*BKDEL2/BKDEL1+B*BLKAP1+B/T3)
      C1N  = (GAMN/DSTR) - SCRK
      C2N  = (GAMN/DSTR) + SCRK
      ALN  = GAMN*S3 + S4*CSQRT(C1N)*CSQRT(C2N)
      T4   = ALN - DEL
      IDX  =-I
      CALL DRKAPM (ALN,IDX,BLKAPM)
      C2   = (ALN-AMU)/T4*AI*SNS/(BETA*(GAMN-ALN*SPS))*BKDEL1/
     1       (BLKAPM)*(-T4/(T4+GL)*A*AI*BKDEL2/BKDEL1+B*BLKAP1+B/T4)
      F1   = F1+C1*AI/T3*(CEXP(-AI*T3*S1)-1.0)+C2*AI/
     1       T4*(CEXP(-AI*T4*S1)-1.0)
      AM1  = AM1+C1/(AI*T3)*(-S1*CEXP(-AI*T3*S1)+AI/
     1       T3*(CEXP(-AI*T3*S1)-1.0))+C2/(AI*T4)*
     2       (-S1*CEXP(-AI*T4*S1)+AI/T4*(CEXP(-AI*T4*S1)-1.0))
      C2A  = C2
      C1A  = C1
      AA   = S1/RL1
      CEXP3 = CEXP(-AI*T3*AA)
      CEXP4 = CEXP(-AI*T4*AA)
      TEMP  = 2.0*PI*R
      CEXP5 = CEXP(AI*(SIGMA-SNS*DEL)/S1*AA)
      CONST = 4.0*FQB/TEMP
      PRES1(1) = PRES1(1) + C1 + C2
      DO 454 JL = 1,NNL1
      CONST = CONST*CEXP5
      C1A   = C1A*CEXP3
      C2A   = C2A*CEXP4
      PRES1(JL+1) = PRES1(JL+1) + C1A + C2A
      PRES1(JL+1) = PRES1(JL+1) + CONST*SIN(TEMP*JL/RL1)
  454 CONTINUE
      IF (CABS((AM1-AMTEST)/AM1) .LT. 0.0005) GO TO 45
      AMTEST = AM1
   20 CONTINUE
      GO TO 9992
 9992 WRITE  (IBBOUT,3005) UFM
 3005 FORMAT (A23,' FROM AMG MODULE. AM1 LOOP IN SUBROUTINE SUBA DID ',
     1        'NOT CONVERGE.')
      CALL MESAGE (-61,0,0)
   45 CONTINUE
      AA    = S1/RL1
      CEXP3 = CEXP(AI*(SIGMA-SNS*DEL)/RL1)
      CONST = FQB
      TEMP  = 2.0*AA/(SPS-SNS)
      PRES1(1) = PRES1(1) - FQB
      DO 4541 JL = 1,NNL1
      CONST = CONST*CEXP3
      PRES1(JL+1) = PRES1(JL+1) - CONST*(1.0-JL*TEMP)
 4541 CONTINUE
      Y    = 0.0
      Y1   = SNS
      ARG  = DEL - GL
      CALL ALAMDA (ARG,Y,BLAM1)
      CALL ALAMDA (ARG,Y1,BLAM2)
      CALL AKAPPA (ARG,BKAP1)
      FT2  = A*AI*(DEL-GL-AMU)*BLAM1/BKAP1
      FT2T = A*AI*(DEL-GL-AMU)*BLAM2/BKAP1
      ARG  = DEL
      CALL ALAMDA (ARG,Y,BLAM1)
      CALL ALAMDA (ARG,Y1,BLAM2)
      CALL AKAPPA (ARG,BKAP1)
      GAM  = SQRT(DEL**2-SCRK**2)
      S5   = SIN(SNS*GAM)
      S6   = COS(SNS*GAM)
      C1   =-1.0/(BETA*GAM*S5)
      C1T  = C1*(AI*SPS*T2*S6-SNS*DEL/GAM*T2*S5)-BLAM2/BKAP1*DEL/GAM*(S5
     1       +GAM*SNS*S6)/(GAM*S5)
      C1   = C1*(ARG/GAM*SNS*S5+AI*SPS*T2)-BLAM1/BKAP1*DEL/(GAM*S5)*(S5/
     1       GAM+SNS*S6)
      FT3  =-B*(BLAM1/BKAP1+(DEL-AMU)*C1)
      FT3T =-B*(BLAM2/BKAP1+(DEL-AMU)*C1T)
      IF (GL .EQ. 0.0) GO TO 50
      F2   = FT2*(CEXP(2.0*AI*GL)-CEXP(AI*GL*S1))/(AI*GL)+
     1       FT3*S0+B*AI*(DEL-AMU)*BLAM1/BKAP1*(4.0-S1**2)/2.0
      AM2  = FT2*(2.0*CEXP(2.0*AI*GL)/(AI*GL)-S1/(AI*GL)*CEXP(GL*AI*S1)+
     1       (CEXP(2.0*AI*GL)-CEXP(AI*S1*GL))/GL**2)+FT3*(4.0-S1**2)/2.0
     2       +B*AI*(DEL-AMU)*BLAM1/BKAP1*(8.0-S1**3)/3.0
      F2P  = FT2T*T1*CEXP(AI*GL*SNS)/(AI*GL)*(CEXP(2.0*AI*GL)-
     1       CEXP(AI*GL*S1))+FT3T*T1*S0+B*AI*(DEL-AMU)*T1*BLAM2/
     2       BKAP1*(S0**2/2.0+SPS*S0)
      AM2P = FT2T*T1*(CEXP(AI*GL*SPS)/(AI*GL)*S0*CEXP(AI*GL*S0)+
     1       CEXP(AI*GL*SPS)/(GL**2)*(CEXP(AI*GL*S0)-1.0))+
     2       FT3T*T1*S0**2/2.0+B*AI*(DEL-AMU)*T1*BLAM2/BKAP1*(S0**3/3.0+
     3       SPS*S0**2/2.0)
      GO TO 55
   50 CONTINUE
      F2   = FT2*S0+FT3*S0+B*AI*(DEL-AMU)*BLAM1/BKAP1*(4.-S1**2)/2.
      AM2  = FT2*(4.0-S1**2)/2.0+FT3*(4.0-S1**2)/2.0+B*AI*(DEL-AMU)*
     1       BLAM1/BKAP1*(8.0-S1**3)/3.0
      F2P  = FT2T*T1*S0+FT3T*T1*S0+B*AI*(DEL-AMU)*T1*BLAM2/BKAP1*(S0**2
     1       /2.0+SPS*S0)
      AM2P = FT2T*T1*S0**2/2.0+FT3T*T1*S0**2/2.0+B*AI*(DEL-AMU)*T1*BLAM2
     1       /BKAP1*(S0**3/3.0+SPS*S0**2/2.0)
   55 CONTINUE
      NL2  = 20
      RL2  = NL2 - 1
      AA   = SPS - SNS
      CONST = B*AI*(DEL-AMU)*BLAM1/BKAP1
      TEMP = S0/RL2
      C1A  = AI*GL
      CEXP3 = CEXP(C1A*AA)
      CEXP4 = CEXP(C1A*TEMP)
      DO 455 JL = 1,NL2
      XL = AA + TEMP*(JL-1)
      PRES2(JL) = FT2*CEXP3 + FT3+CONST*XL
      CEXP3 = CEXP3*CEXP4
  455 CONTINUE
      CALL SUBBB
      RETURN
      END