!Copyright (c) 2015 M. Zapukhlyak
!THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
!FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
!COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
!IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
!CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

!-----------------------------------------------------------------------
!Module factory_method_module implements factory method design pattern
!-----------------------------------------------------------------------
module factory_method_module
   !use
   implicit none
   private ! all by default
   public :: abstract_product, concrete_product_a, concrete_product_b
   public :: abstract_creator, concrete_creator_a, concrete_creator_b

   ! product definitions
   type, abstract :: abstract_product
   contains
      procedure (abstract_product_procedure), deferred, pass :: get_name
   end type abstract_product

   abstract interface
      function abstract_product_procedure(this) result(name)
         import  :: abstract_product
         class(abstract_product), intent(inout) :: this
         character(len=:), allocatable :: name
      end function abstract_product_procedure
   end interface

   type, extends(abstract_product) :: concrete_product_a
   contains
      procedure, pass :: get_name => get_name_concrete_product_a
   end type concrete_product_a

   type, extends(abstract_product) :: concrete_product_b
   contains
      procedure, pass :: get_name => get_name_concrete_product_b
   end type concrete_product_b

   ! creator definitions
   type, abstract :: abstract_creator
   contains
      procedure (abstract_creator_procedure), deferred, pass :: factory_method
   end type abstract_creator

   abstract interface
      subroutine abstract_creator_procedure(this, prod)
         import  :: abstract_creator
         import  :: abstract_product
         class(abstract_creator), intent(in) :: this
         class(abstract_product), allocatable, intent(out) :: prod
      end subroutine abstract_creator_procedure
   end interface

   type, extends(abstract_creator) :: concrete_creator_a
   contains
      procedure, pass :: factory_method => factory_method_concrete_creator_a
   end type concrete_creator_a

   type, extends(abstract_creator) :: concrete_creator_b
   contains
      procedure, pass :: factory_method => factory_method_concrete_creator_b
   end type concrete_creator_b

contains
   !-----------------------------------------------------------------------
   !Function
   !-----------------------------------------------------------------------
   function  get_name_concrete_product_a(this) result(name)
      implicit none
      class(concrete_product_a), intent(inout) :: this
      character(len=:),allocatable :: name
      name = "concrete_product_a"
   end function get_name_concrete_product_a
   !-----------------------------------------------------------------------
   !Function
   !-----------------------------------------------------------------------
   function  get_name_concrete_product_b(this) result(name)
      implicit none
      class(concrete_product_b), intent(inout) :: this
      character(len=:),allocatable :: name
      name = "concrete_product_b"
   end function get_name_concrete_product_b
   !-----------------------------------------------------------------------
   !Function
   !-----------------------------------------------------------------------
   subroutine  factory_method_concrete_creator_a(this, prod)
      implicit none
      class(concrete_creator_a), intent(in) :: this
      class(abstract_product), allocatable, intent(out) :: prod
      ! prod = concrete_product_a() ! Assignment to an allocatable polymorphic variable at (1) is not yet supported
      allocate(prod, source=concrete_product_a())
   end subroutine factory_method_concrete_creator_a
   !-----------------------------------------------------------------------
   !Function
   !-----------------------------------------------------------------------
   subroutine  factory_method_concrete_creator_b(this, prod)
      implicit none
      class(concrete_creator_b), intent(in) :: this
      class(abstract_product), allocatable, intent(out) :: prod
      ! prod = concrete_product_b() ! Assignment to an allocatable polymorphic variable at (1) is not yet supported
      allocate(prod, source=concrete_product_b())
   end subroutine factory_method_concrete_creator_b

end module factory_method_module

!-----------------------------------------------------------------------
!Main program test_factory_method_module
!-----------------------------------------------------------------------
program    test_factory_method_module
   use factory_method_module
   implicit none
   integer, parameter :: n = 2
   type(concrete_creator_a), target :: creator_a
   type(concrete_creator_b), target :: creator_b

   type creator_container
      class(abstract_creator), pointer :: c
   end type creator_container

   class(creator_container), allocatable, dimension(:) :: creators
   class(abstract_product), allocatable :: prod
   integer :: i

   call creator_a % factory_method(prod)
   write(*,*) prod % get_name()

   select type (prod)
    class is (abstract_product)
      write(*,*) "class is abstract_product???"
    type is (concrete_product_a)
      write(*,*) "type is concrete_product_a"
    type is (concrete_product_b)
      write(*,*) "type is concrete_product_b"
    class default
      write(*,*) "type is unknown?"
   end select

   deallocate(prod)

   call creator_b % factory_method(prod)

   write(*,*) prod % get_name()
   select type (prod)
    class is (abstract_product)
      write(*,*) "class is abstract_product???"
    type is (concrete_product_a)
      write(*,*) "type is concrete_product_a"
    type is (concrete_product_b)
      write(*,*) "type is concrete_product_b"
    class default
      write(*,*) "type is unknown?"
   end select

   deallocate(prod)

   allocate(creators(n))
   creators(1) % c => creator_a
   creators(2) % c => creator_b

   do i = 1, n
      call creators(i) % c % factory_method(prod)
      write(*,*) prod % get_name()
      deallocate(prod)
   enddo

end program test_factory_method_module


text:
PROGRAM my_prog
i = 1
END PROGRAM my_prog

max_stars_repo_name: ouankou/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


! Reads the inputs to create the mesh
subroutine mesh
use Declarations
real(dp) :: delta_x, &    ! Horizontal space step [m]
        &   delta_y       ! Vertical space ste [m]
integer :: choice_nb_layers ! Variable to decide if vertical space step is equal to layer thicknesses (=1) or not (=0)
real(dp), allocatable :: y(:), &        ! Variable used to calculate yn
                     &   thickness(:)   ! Variable used to calculate yn

choice_nb_layers = 0

open(7,file='inputs.txt')
read(7,*) L1
read(7,*) L2
read(7,*) N
read(7,*) M
read(7,*)
read(7,*)
read(7,*) beta
read(7,*) nb_layers
allocate (thickness(nb_layers))
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
read(7,*)
do j = 1,nb_layers
   read(7,*) thickness(j)
enddo

close (7)

if (M .eq. 999) then   ! use layer thicknesses for vertical mesh resolution
   M = nb_layers + 1
   choice_nb_layers = 1
endif

!================== Allocate mesh variables =================================
allocate (y(M),xn(N,M),yn(N,M),xv(N-1,M-1),yv(N-1,M-1), Vol(N-1,M-1))
allocate (Sx_west(N-1,M-1),Sy_west(N-1,M-1),Sx_east(N-1,M-1),Sy_east(N-1,M-1))
allocate (Sx_north(N-1,M-1),Sy_north(N-1,M-1),Sx_south(N-1,M-1),Sy_south(N-1,M-1))

y(1) = 0._dp

if (choice_nb_layers .eq. 0) then
   delta_y = L2 / real(M-1,dp)
   do j = 1,M-1
      y(j+1) = y(j) + delta_y
   enddo
else
   do j=1,M-1
      y(j+1) = y(j)+thickness(j)
   enddo
endif

do j = 1,M
   xn(1,j) = dsin(beta)*y(j)
   yn(1,j) = dcos(beta)*y(j)
enddo

delta_x = L1 / real(N-1,dp)

do j = 1,M
   do i = 2,N
      xn(i,j) = xn(1,j) + real((i-1),dp) * delta_x * dcos(beta)
      yn(i,j) = yn(1,j) - dsin(beta)*real(i-1,dp)*delta_x
   enddo
enddo
do j = 1,M-1
   do i = 1,N-1
      xv(i,j) = (xn(i,j)+xn(i+1,j)+xn(i,j+1)+xn(i+1,j+1))*0.25_dp
      yv(i,j) = (yn(i,j)+yn(i+1,j)+yn(i,j+1)+yn(i+1,j+1))*0.25_dp
      Vol(i,j) = (yn(i,j+1)-yn(i,j))*delta_x/dcos(beta)
      Sx_west(i,j) = dabs(yn(i,j+1)-yn(i,j))
      Sy_west(i,j) = dabs(xn(i,j+1)-xn(i,j))
      Sx_east(i,j) = dabs(yn(i+1,j+1)-yn(i+1,j))
      Sy_east(i,j) = dabs(xn(i+1,j+1)-xn(i+1,j))
      Sx_north(i,j) = dabs(yn(i+1,j+1)-yn(i,j+1))
      Sy_north(i,j) = dabs(xn(i+1,j+1)-xn(i,j+1))
      Sx_south(i,j) = dabs(yn(i+1,j)-yn(i,j))
      Sy_south(i,j) = dabs(xn(i+1,j)-xn(i,j))
    enddo
enddo

! Check if sum of all small volumes is equal to rectangle volume
write(*,*) 'Sum of all small volumes [m2] = ', sum(Vol)
write(*,*) 'Total volume [m2] = ', L1 * L2

deallocate(y, thickness)

return
end


text:
! =============================================================================
!                       Test FFT module
!
!       This unit test checks the FFT module with reference solutions.
! =============================================================================
program test_fft_2
    use unit_test
    use constants, only : pi, twopi, zero, one
    use stafft, only : initfft, forfft
    implicit none

    integer, parameter          :: n = 16
    double precision            :: trig(2 * n)
    integer                     :: factors(5)
    double precision            :: fun(n)
    double precision            :: ref(2 * n)            ! reference solution in Fourier space
    double precision            :: a ! 0 < a < 1
    double precision            :: da = 0.1d0
    double precision            :: err = zero
    double precision, parameter :: atol = dble(1.0e-14)

    ! set up FFT
    call initfft(n, factors, trig)

    ! test function 1
    a = da
    do while (a <= 0.9d0)

        ! initialize function
        fun = func_1(a, n)

        ! compute reference solution
        call dft(n, fun, ref)

        ! forward FFT (Hermitian form) (overwrites fun)
        call forfft(1, n, fun, trig, factors)

        err = max(err, get_max_absolute_error(n, fun, ref))

        a = a + da
    enddo

    ! test function 2
    a = da
    do while (a <= 0.9d0)

        ! initialize function
        fun = func_2(a, n)

        ! compute reference solution
        call dft(n, fun, ref)

        ! forward FFT (Hermitian form) (overwrites fun)
        call forfft(1, n, fun, trig, factors)

        err = max(err, get_max_absolute_error(n, fun, ref))

        a = a + da
    enddo

    ! final check

    call print_result_dp('Test FFT transform', err, atol)

    contains

        ! initialize function in real space from -pi to pi
        ! f(x) = 1 / (1 - a * cos(x + 1))
        ! 0 < a < 1
        function func_1(a, n) result(f1)
            double precision, intent(in) :: a
            integer,          intent(in) :: n
            double precision             :: f1(n)
            double precision             :: dx, xx, phase
            integer                      :: j

            phase = one ! phase offset in radians

            dx = twopi / dble(n)
            do j = 0, n-1
                xx = -pi + dx * dble(j)
                f1(j+1) = one / (one - a * dcos(xx + phase))
            enddo
        end function func_1

        ! initialize function in real space from -pi to pi
        ! f(x) = 1 / (1 - a * sin(x + 1))
        ! 0 < a < 1
        function func_2(a, n) result(f2)
            double precision, intent(in) :: a
            integer,          intent(in) :: n
            double precision             :: f2(n)
            double precision             :: dx, xx, phase
            integer                      :: j

            phase = one ! phase offset in radians

            dx = twopi / dble(n)
            do j = 0, n-1
                xx = -pi + dx * dble(j)
                f2(j+1) = one / (one - a * dsin(xx + phase))
            enddo
        end function func_2

        ! discrete Fourier transform
        subroutine dft(m, f, fhat)
            integer,          intent(in)  :: m
            double precision, intent(in)  :: f(m)       ! function in real space
            double precision, intent(out) :: fhat(2*m)  ! function in Fourier space
            double precision              :: fac, arg
            integer                       :: k, l
            double precision              :: sfac       ! normalization factor

            fhat = 0.d0
            fac = twopi / dble(m)

            do k = 0, m - 1
                do l = 0, m - 1
                    arg = fac * dble(k * l)
                    fhat(k+1)   = fhat(k+1)   + f(l+1) * dcos(arg)
                    fhat(k+1+m) = fhat(k+1+m) - f(l+1) * dsin(arg)
                enddo
            enddo

            sfac = one / dsqrt(dble(m))

            fhat = fhat * sfac

        end subroutine dft

        function get_max_absolute_error(n, fun, ref) result(err)
            integer,          intent(in) :: n
            double precision, intent(in) :: fun(n), ref(2 * n)
            double precision             :: err, err_cos, err_sin
            integer                      :: i, j, k, nw

            nw = n / 2

            ! Cosine coefficients
            err_cos = maxval(abs(ref(1:nw+1) - fun(1:nw+1)))

            ! Sine coefficients
            i = 1
            j = nw - 1
            k = n + 1

            err_sin = maxval(abs(ref(k+i:k+j) - fun(k-i:k-j:-1)))

            err = max(err_cos, err_sin)

        end function get_max_absolute_error

end program test_fft_2

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


subroutine readXsection(k,rmanning,timesDepth)

    use constants_module
    use arrays_module
    use arrays_section_module
    use xsec_attribute_module

    implicit none
    save

    integer, intent(in) :: k
    real, intent(in) :: rmanning, timesDepth

    real xcs(maxTableLength), ycs(maxTableLength), el1(nel),a1(nel),peri1(nel),redi1(nel)
    real conv1(nel), tpW1(nel), diffArea(nel), newI1(nel), deffPere(nel), newdPdA(nel)
    integer i_start(nel), i_end(nel), i_area, i_find, i, j, jj, num
    real el_min, el_max, el_range, el_incr, el_now, x1, y1, x2, y2, x_start, x_end, waterElev
    real f2m, cal_area, cal_peri, cal_topW, cal_dist, cal_tri_area, cal_multi_area, cal_perimeter, diffAreaCenter
    integer i1, i2
    !character*4 file_num

    ! Assign some parameters
    ! f2m is conversion if the data is in feet

    f2m=1.0

!     Open data file

    write(file_num,'(i4.4)')k
        print*, file_num
        open(11,file=trim(xSection_path)//file_num)
!     Skip headings

      read(11,*)

!     Get CS data
      do i=2,maxTableLength
          read(11,*,end=300)x1,y1
          xcs(i)=x1*f2m
          ycs(i)=y1*f2m
      enddo
300   close(11)
      num=i

!     get el_max, el_min, el_range,el_incr nel=100
      el_min=99999.
      el_max=-99999.
      do i=2,num-1
          if(ycs(i).lt.el_min)el_min=ycs(i)
          if(ycs(i).gt.el_max)el_max=ycs(i)
      enddo
      el_range=(el_max-el_min)*timesDepth
      el_incr=el_range/real(nel-1.0)

      xcs(1)=xcs(2)
      ycs(1)=el_min+el_range+1.
      xcs(num)=xcs(num-1)
      ycs(num)=el_min+el_range+1.

!     output cs data
      open(11,file=trim(xSection_path)//file_num//'.dat')
      do i=1,num
          write(11,*)xcs(i),ycs(i) !changing all into m unit
      enddo
      close(11)

      open(11,file=trim(xSection_path)//file_num//'_lines')
      open(22,file=trim(xSection_path)//file_num//'_tab')
      write(22,'(120a)')' Elev(m)    Area(m2)     Peri(m)      Radi(m)   Conv(m3/s)    topWidth(m)    newI1(m3)    dPdA(1/m)'  ! Hu changed

      do j=1,nel
          el_now=el_min+real(j-1)*el_incr

          if(abs(el_now - el_min) < TOLERANCE) then
            el_now=el_now+0.00001
          end if

          i_start(1)=-999
          i_end(1)=-999
          i_area=0
          i_find=0
        do i=1,num-1
          y1=ycs(i)
          y2=ycs(i+1)
          if(el_now.le.y1 .and. el_now.gt.y2 .and. i_find.eq.0)then
              i_find=1
              i_area=i_area+1
              i_start(i_area)=i
          endif
          if(el_now.gt.y1 .and. el_now.le.y2 .and. i_find.eq.1)then
              i_find=0
              i_end(i_area)=i
          endif
        enddo

        cal_area=0.
        cal_peri=0.
        cal_topW=0.
                newI1=0.0 !Hu changed

        do i=1,i_area
            x1=xcs(i_start(i))
            x2=xcs(i_start(i)+1)
            y1=ycs(i_start(i))
            y2=ycs(i_start(i)+1)
            if(y1.eq.y2)then
                x_start=x1
            else
                x_start=x1+(el_now-y1)/(y2-y1)*(x2-x1)
            endif
            write(11,*)x_start,el_now

            x1=xcs(i_end(i))
            x2=xcs(i_end(i)+1)
            y1=ycs(i_end(i))
            y2=ycs(i_end(i)+1)

            if(y1.eq.y2)then
              x_end=x1
            else
              x_end=x1+(el_now-y1)/(y2-y1)*(x2-x1)
            endif

            cal_topW=x_end-x_start+cal_topW

            write(11,*)x_end,el_now
            write(11,*)'NaN NaN'

            i1=i_start(i)
            i2=i_end(i)

            cal_area = cal_area    &
                     +cal_tri_area(el_now,x_start,xcs(i1+1),ycs(i1+1))    &
                     +cal_multi_area(el_now,xcs,ycs,maxTableLength,i1+1,i2)    &
                     +cal_tri_area(el_now,x_end,xcs(i2),ycs(i2))
            cal_peri = cal_peri    &
                    +cal_dist(x_start,el_now,xcs(i1+1),ycs(i1+1))    &
                    +cal_perimeter(xcs,ycs,maxTableLength,i1+1,i2)    &
                    +cal_dist(x_end,el_now,xcs(i2),ycs(i2))
            if(i1.eq.1)cal_peri=cal_peri    &
                     -cal_dist(x_start,el_now,xcs(i1+1),ycs(i1+1))
            if(i2.eq.(num-1))cal_peri=cal_peri    &
                     -cal_dist(x_end,el_now,xcs(i2),ycs(i2))

        end do

        el1(j)=el_now
        a1(j)=cal_area
        peri1(j)=cal_peri
        redi1(j)=a1(j)/peri1(j)
        conv1(j)=1./rmanning*a1(j)*(redi1(j))**(2./3.)
        tpW1(j)=cal_topW

        if(j.eq.1) then
          diffArea(j)=a1(j)
          deffPere(j)=peri1(j)
        else
          diffArea(j)=a1(j)-a1(j-1)
          deffPere(j)=peri1(j)-peri1(j-1)
        endif

        newdPdA(j)=deffPere(j)/diffArea(j)

        waterElev=el1(j)

        do jj=2,j
          diffAreaCenter=el1(jj)-el_incr*0.5
          newI1(j)=newI1(j)+diffArea(jj)*(waterElev-diffAreaCenter)

        enddo

        write(22,10)el1(j),a1(j),peri1(j),redi1(j),conv1(j),    &
                   tpW1(j),newI1(j),newdPdA(j)
10      format(f9.2,3f12.2,2f20.3,2f16.3)

        ! new change 20200107 to make the model faster
        xsec_tab(1,j,k) = el1(j)
        xsec_tab(2,j,k) = a1(j)
        xsec_tab(3,j,k) = peri1(j)
        xsec_tab(4,j,k) = redi1(j)
        xsec_tab(5,j,k) = conv1(j)
        xsec_tab(6,j,k) = tpW1(j)
        xsec_tab(7,j,k) = newI1(j)
        xsec_tab(8,j,k) = newdPdA(j)
      end do

      close(11)
      close(22)

      z(k)= el_min
end subroutine

function cal_tri_area(el,x0,x1,y1)
      cal_tri_area=abs(0.5*(x1-x0)*(el-y1))
      return
end function

function cal_trap_area(el,x1,y1,x2,y2)
      cal_trap_area=abs(0.5*(x2-x1)*(el-y1+el-y2))
      return
end function

function cal_multi_area(el,xx,yy,n,i1,i2)
      real xx(n),yy(n)

      area=0
      do i=i1,i2-1
          x1=xx(i)
          y1=yy(i)
          x2=xx(i+1)
          y2=yy(i+1)
          area=area+cal_trap_area(el,x1,y1,x2,y2)
      enddo
      cal_multi_area=area
      return
end function

function cal_dist(x1,y1,x2,y2)
          cal_dist=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+1.e-32)
      return
end function

function cal_perimeter(xx,yy,n,i1,i2)
        real xx(n),yy(n)

      p=0.
      do i=i1,i2-1
        x1=xx(i)
        y1=yy(i)
        x2=xx(i+1)
        y2=yy(i+1)
        p=p+cal_dist(x1,y1,x2,y2)
      enddo
      cal_perimeter=p
      return
end function


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 06.01.01/04                           *
C  *    TEST TITLE : Behavior of matrix composition        *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      REAL      ACT4X4(4,4), EXP4X4(4,4), AIN4X4(4,4), BIN4X4(4,4)
      REAL      ACT3X3(3,3), EXP3X3(3,3), AIN3X3(3,3), BIN3X3(3,3)
      REAL      INC
      INTEGER   I,J
      LOGICAL   TRNSEQ

      CALL INITGL ('06.01.01/04')

C open PHIGS
      CALL XPOPPH  (ERRFIL, MEMUN)
      CALL SETMSG ('1 2 11', 'The matrix returned by <compose ' //
     1             'matrix 3> should perform the same 3D '      //
     2             'transformation as that resulting from the ' //
     3             'multiplication of the two specified matrices.')

      AIN4X4(1,1) = -2.2
      AIN4X4(1,2) = 3.3
      AIN4X4(1,3) = 4.4
      AIN4X4(1,4) = 5.5
      AIN4X4(2,1) = 0.0023
      AIN4X4(2,2) = -727.12
      AIN4X4(2,3) = 2.09
      AIN4X4(2,4) = 33.33
      AIN4X4(3,1) = 0
      AIN4X4(3,2) = 34
      AIN4X4(3,3) = 54.1
      AIN4X4(3,4) = -0.045
      AIN4X4(4,1) = 7
      AIN4X4(4,2) = 6
      AIN4X4(4,3) = 5
      AIN4X4(4,4) = 4

      INC = 0.0
      DO 401 I = 1,4
      DO 402 J = 1,4
        BIN4X4(I,J) = TAN(-1.56 + INC)
        INC = INC + 3.0/16.0
402   CONTINUE
401   CONTINUE

      CALL PCOM3  (AIN4X4, BIN4X4, ERRIND, ACT4X4)
      CALL CHKINQ ('pcom3', ERRIND)

C  Compute exp4x4 = expected array
      CALL ECOM3  (AIN4X4, BIN4X4, EXP4X4)

C  pass if matrices represent same transformation
      CALL IFPF   (TRNSEQ(4, ACT4X4, EXP4X4))

      CALL SETMSG ('1 2 12', 'The matrix returned by <compose ' //
     1             'matrix> should perform the same 2D '        //
     2             'transformation as that resulting from the ' //
     3             'multiplication of the two specified matrices.')

      AIN3X3(1,1) = -2.2
      AIN3X3(1,2) = 3.3
      AIN3X3(1,3) = 5.5
      AIN3X3(2,1) = 0.0023
      AIN3X3(2,2) = -727.12
      AIN3X3(2,3) = 33.33
      AIN3X3(3,1) = 0
      AIN3X3(3,2) = 34
      AIN3X3(3,3) = -0.045

      INC = 0.0
      DO 301 I = 1,3
      DO 302 J = 1,3
        BIN3X3(I,J) = TAN(-1.56 + INC)
        INC = INC + 3.0/9.0
302   CONTINUE
301   CONTINUE

      CALL PCOM  (AIN3X3, BIN3X3, ERRIND, ACT3X3)
      CALL CHKINQ ('pcom', ERRIND)

C  Compute exp3x3 = expected array
      CALL ECOM  (AIN3X3, BIN3X3, EXP3X3)

C  pass if matrices represent same transformation
      CALL IFPF   (TRNSEQ(3, ACT3X3, EXP3X3))

      CALL ENDIT
      END

ext: f


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
subroutine setlg(lg_sta,lg_end,lg_num,ista_Mx_ori,iend_Mx_ori,inum_Mx_ori,    &
                 Hgs,Nd,rLsize1,imesh_oddeven,iperiodic)
  implicit none
  integer :: iperiodic
  integer :: lg_sta(3),lg_end(3),lg_num(3)
  integer :: ista_Mx_ori(3),iend_Mx_ori(3),inum_Mx_ori(3)
  integer :: Nd
  integer :: imesh_oddeven(3)
  integer :: jj
  real(8) :: Hgs(3)
  real(8) :: rLsize1(3)
  real(8),parameter :: epsilon=1.d-10

  select case(iperiodic)
    case(0)
      iend_Mx_ori(:)=int((rLsize1(:)+epsilon)/2.d0/Hgs(:))+Nd
      lg_end(:)=int((rLsize1(:)+epsilon)/2.d0/Hgs(:))

      do jj=1,3
        select case(imesh_oddeven(jj))
          case(1)
            ista_Mx_ori(jj)=-(int((rLsize1(jj)+epsilon)/2.d0/Hgs(jj))+Nd)
            lg_sta(jj)=-(int((rLsize1(jj)+epsilon)/2.d0/Hgs(jj)))
          case(2)
            ista_Mx_ori(jj)=-(int((rLsize1(jj)+epsilon)/2.d0/Hgs(jj))+Nd)+1
            lg_sta(jj)=-(int((rLsize1(jj)+epsilon)/2.d0/Hgs(jj)))+1
        end select
      end do

    case(3)
      ista_Mx_ori(:)=1-Nd
      lg_sta(:)=1
      iend_Mx_ori(:)=int((rLsize1(:)+epsilon)/Hgs(:))+Nd
      lg_end(:)=int((rLsize1(:)+epsilon)/Hgs(:))
  end select

  inum_Mx_ori(:)=iend_Mx_ori(:)-ista_Mx_ori(:)+1
  lg_num(:)=lg_end(:)-lg_sta(:)+1

end subroutine setlg

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
program foo
  integer, pointer :: iptr1, iptr2
  integer, target :: i1, i2

  iptr1 => i1
  iptr2 => i2

  nullify(iptr1, iptr2)
end program foo

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
module MyMod1
  implicit none

  integer, parameter :: myModVal1 = 1

end module MyMod1

max_stars_repo_path: test cases/fortran/8 module names/mod1.f90


text:
	SUBROUTINE FL_TINQ  ( table, type, exist, fname, iret )
C************************************************************************
C* FL_TINQ								*
C* 									*
C* This subroutine returns the path to table file, if the table exists.	*
C* 									*
C* The TABLE file name is split into the path and filename and the file	*
C* is located by searching in the following order:			*
C* 									*
C*	1. filename (local)						*
C* 	2. path/filename (TABLE as given)				*
C* 	3. $NCDESK/type/filename					*
C* 	4. $NCSITE/type/filename					*
C*	5. $GEMTBL/type/filename					*
C* 									*
C* FL_TINQ  ( TABLE, TYPE, EXIST, FNAME, IRET )				*
C* 									*
C* Input parameters:							*
C*	TABLE		CHAR*		Table file name 		*
C*	TYPE		CHAR*		File name type			*
C*									*
C* Output parameters:							*
C*	EXIST		LOGICAL		Table exists flag		*
C*	FNAME		CHAR*		Full file name for the table	*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 5/01	Copied from FL_TBOP			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	table, type, fname
	LOGICAL		exist
C
	CHARACTER	newfil*132, name*132
	CHARACTER       filnam*132, name3*132, name4*132, name5*132
	CHARACTER	patfil*132, path*132
	LOGICAL		found
C------------------------------------------------------------------------
	iret  = 0
	found = .false.
	fname = ' '
	name  = ' '
C
	IF  ( table .eq. ' ' )  THEN
	    iret  = -1
	    RETURN
	END IF
C
C*	Strip GEMTBL from beginning of the table file name if present.
C
	CALL ST_RMST ( table, '$GEMTBL/', ipos, patfil, ier )
	IF ( ipos .eq. 0 ) THEN
	    CALL ST_RMST ( table, 'GEMTBL:', ipos, patfil, ier )
	END IF
	IF ( ipos .eq. 0 ) patfil = table
C
C*	Extract the filename from the end of the path.
C
	CALL FL_PATH ( patfil, path, filnam, ier )
	CALL ST_LSTR ( filnam, lens1, ier )
	CALL ST_LSTR ( type, lenst, ier )
C
C*	Check to see if the file is local.
C
	CALL FL_INQR ( filnam, found, newfil, iret )
	IF ( found ) THEN
	    name  = filnam (1:lens1)
	END IF
C
C*	Check to see if file is down a path.
C
	IF ( ( .not. found ) .and. ( path .ne. ' ' ) ) THEN
	    CALL FL_INQR ( patfil, found, newfil, iret )
	    IF ( found ) THEN
		name  = patfil
	    END IF
	END IF
C
C*	Check to see if file exists in terms of type in NCDESK.
C
	IF ( ( .not. found ) .and. type .ne. ' ' ) THEN
	    name3 = '$NCDESK/' // type (1:lenst) // '/' //
     +		     filnam (1:lens1)
	    CALL FL_INQR  ( name3, found, newfil, iret )
  	    IF  ( found )  THEN
		name  = name3
	    END IF
	END IF
C
C*	Check to see if file exists in terms of type in NCSITE.
C
	IF ( ( .not. found ) .and. type .ne. ' ' ) THEN
	    name4 = '$NCSITE/' // type (1:lenst) // '/' //
     +		     filnam (1:lens1)
	    CALL FL_INQR  ( name4, found, newfil, iret )
  	    IF  ( found )  THEN
		name  = name4
	    END IF
	END IF
C
C*	Check to see if file exists in terms of type in GEMTBL.
C
	IF ( .not. found .and. type .ne. ' ' ) THEN
	    name5 = '$GEMTBL/' // type (1:lenst) // '/' //
     +		     filnam (1:lens1)
	    CALL FL_INQR ( name5, found, newfil, iret )
	    IF ( found ) THEN
		name  = name5
	    END IF
	END IF
C
	fname = name
	exist = found
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


text:
C++********************************************************************
C Set of routines:
C JLP_FOURIER
C
C JLP
C Version 11-07-91
C--*******************************************************************
C JLP_FOURIER
C Computes the Fourier Transform of a real image (full of zeroes...)
C
C Zero frequency is in RE(1,1), and IM(1,1)
C
C Input:
C IMAGE: Input image
C
C Output:
C RE,IM: Real and imaginary part
C
C RE(KX,KY) = Sum on IX,IY of IMAGE(IX,IY)*COS(2*PI*(IX*(KX-1)/NX+IY*(KY-1)/NY))
C IM(KX,KY) = Sum on IX,IY of IMAGE(IX,IY)*SIN(2*PI*(IX*(KX-1)/NX+IY*(KY-1)/NY))
C
C Tests: Big improvement in accuracy
C          when X0,X1,X0,Y1,PI,ANGLE,MYCOS, and MYSIN are real*8...
C          and when DCOS and DSIN are used.
C*******************************************************************
	SUBROUTINE JLP_FOURIER(IMAGE,RE,IM,NX,NY,IDIM,ISIGN)
	PARAMETER (IDIM1=512, SMALLER_VALUE=1.E-12)
	REAL IMAGE(IDIM,*),RE(IDIM,*),IM(IDIM,*)
	REAL*8 X0,Y0,X1,Y1,SIGN
	REAL*8 PI,ANGLE
	INTEGER NX,NY,ISIGN
	REAL*8 MYCOS(IDIM1),MYSIN(IDIM1)

	PI=3.14159265358979323846
	X0=2.*PI/FLOAT(NX)
	Y0=2.*PI/FLOAT(NY)
	SIGN=FLOAT(ISIGN)

C Resetting arrays to zero:
	DO KY=1,NY
	  DO KX=1,NX
	    RE(KX,KY)=0.
	    IM(KX,KY)=0.
	  END DO
	END DO

        IF(NX.EQ.NY)THEN
C When NX=NY, we have:
C RE(KX,KY) = Sum on IX,IY of IMAGE(IX,IY)*COS((2*PI/NX)*(IX*(KX-1)+IY*(KY-1)))
C IM(KX,KY) = Sum on IX,IY of IMAGE(IX,IY)*SIN((2*PI/NX)*(IX*(KX-1)+IY*(KY-1)))
C
C Cosinus and sinus:
C MYCOS(I) = COS(2*PI*(I-1)/NX)
         DO I=1,NX
	   ANGLE=X0*FLOAT(I-1)
	   MYCOS(I)=DCOS(ANGLE)
	   MYSIN(I)=DSIN(ANGLE)
         END DO

C Fourier Transform with pre-computed cosinus and sinus :
	DO IY=1,NY
	  DO IX=1,NX
	    IF(ABS(IMAGE(IX,IY)).GT.SMALLER_VALUE)THEN
	      DO KY=1,NY
	        DO KX=1,NX
                  I=(IX-1)*(KX-1)+(IY-1)*(KY-1)
                  I=MOD(I,NX)+1
                  RE(KX,KY)=RE(KX,KY)+IMAGE(IX,IY)*MYCOS(I)
                  IM(KX,KY)=IM(KX,KY)+IMAGE(IX,IY)*MYSIN(I)
	        END DO
	      END DO
	    ENDIF
	  END DO
	END DO

        ELSE

C Fourier Transform without pre-computed cosinus and sinus:
	DO IY=1,NY
	  Y1=Y0*FLOAT(IY-1)
	  DO IX=1,NX
	    X1=X0*FLOAT(IX-1)
	    IF(ABS(IMAGE(IX,IY)).GT.SMALLER_VALUE)THEN
	      DO KY=1,NY
	        DO KX=1,NX
	          ANGLE=X1*FLOAT(KX-1)+Y1*FLOAT(KY-1)
                  RE(KX,KY)=RE(KX,KY)+IMAGE(IX,IY)*DCOS(ANGLE)
                  IM(KX,KY)=IM(KX,KY)+IMAGE(IX,IY)*DSIN(ANGLE)
	        END DO
	      END DO
	    ENDIF
	  END DO
	END DO

        ENDIF

	RETURN
	END

ext: for

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib