Artificial Viscosity Smoothening

dimension Q(100,100),R(100,100),RTranspose(100,100),RTransposeR(100,100),RTransposeRInverse(100,100),RTransposeQ(100,100),U(100,100)
    integer qrow,qcol,rrow,rcol
    real epsilon_i_plus_1by2,epsilon_i_minus_1by2,nu_i,nu_i_plus_1,K,tmp !variables for artificial smoothening
    K=3.5 !trial artificial viscosity for smoothening. change until the unit hydrograph is smooth enough

    qrow=16 !16 non-zero DRH ordinates
    qcol=1
    rrow=16 !16 rows of Rainfall matrix
    rcol=14 !14 rows of Rainfall matrix (qrow-(no. of ERH ordinates-1)) = 16-(3-1) = 16-2 = 14
    open(5,file='Q.txt')
    do i=1,qrow
        read (5,*)(Q(i,j),j=1,qcol)
    end do
    open(5,file='R.txt')
    do i=1,rrow
        read (5,*)(R(i,j),j=1,rcol)
    end do

    call Transpose(R,RTranspose,rrow,rcol)
    call MatMul(RTranspose,R,RTransposeR,rcol,rrow,rrow,rcol)
    call Inverse(RTransposeR,RTransposeRInverse,rcol)
    call MatMul(RTranspose,Q,RTransposeQ,rcol,rrow,qrow,qcol)
    call MatMul(RTransposeRInverse,RTransposeQ,U,rcol,rcol,rcol,1)

    open(5,file='U.txt')
    do i=1,rcol !U vector has a size of rcol

        !this mess inside the if block below is for artificial smoothening
        if(i.GT.10) then !filter the UH ordinates that you want to smoothen(here, UH ordinates above the 10th)
            if(i.NE.rcol) then !if not the last value
                nu_i_plus_1=abs(U(i+1,1)-2*U(i,1)+U(i-1,1))/(abs(U(i+1,1))+2*abs(U(i,1))+abs(U(i-1,1)))
                nu_i=abs(U(i,1)-2*U(i-1,1)+U(i-2,1))/(abs(U(i,1))+2*abs(U(i-1,1))+abs(U(i-2,1)))
                nu_i_minus_1=abs(U(i-1,1)-2*U(i-2,1)+U(i-3,1))/(abs(U(i-1,1))+2*abs(U(i-2,1))+abs(U(i-3,1)))
            else !if last value
                nu_i_plus_1=abs(U(i,1)+U(i-1,1))/(abs(U(i,1))+abs(U(i-1,1)))
                nu_i=abs(U(i-1,1)+U(i-2,1))/(abs(U(i-1,1))+abs(U(i-2,1)))
                nu_i_minus_1=abs(U(i-2,1)+U(i-3,1))/(abs(U(i-2,1))+abs(U(i-3,1)))
            endif

            if(nu_i_plus_1.GT.nu_i) then
                tmp=nu_i_plus_1
            else
                tmp=nu_i
            endif
            epsilon_i_plus_1by2=K*tmp

            if(nu_i_minus_1.GT.nu_i) then
                tmp=nu_i_minus_1
            else
                tmp=nu_i
            endif
            epsilon_i_minus_1by2=K*tmp

            U(i,1)=U(i,1)+epsilon_i_plus_1by2*(U(i+1,1)-U(i,1))-epsilon_i_minus_1by2*(U(i,1)-U(i-1,1))
        endif

        write (5,*) U(i,1)
    end do

    end program


subroutine Transpose(A,R,rows,cols) !rows,cols = of A
integer rows,cols
dimension A(100,100),R(100,100) !CAUTION: This declaration is crucial
    do i=1,cols !i goes from 1 to cols
        do j=1,rows !and j from 1 to rows because transpose should have the dimension: cols x rows
            R(i,j)=A(j,i) !fills from column of A to row of R
        end do
    end do
    return
end subroutine Transpose


subroutine MatMul(A,B,C,arow,acol,brow,bcol)
integer arow,acol,brow,bcol
!To see how this works, try writing the elements of the product matrix in sigma notation
dimension A(100,100),B(100,100),C(100,100) !CAUTION: This declaration is crucial
do k=1,arow !k runs from 1 to arow
    do l=1,bcol !and l from 1 to bcol because that's the output matrix dimension and we _have_ to populate each element of it
        sum=0
        do i=1,acol !this loop is the essence of the aforementioned sigma notation for the (k,l) element of the product matrix
            sum=sum+A(k,i)*B(i,l)
        end do
        C(k,l)=sum
    end do
end do
return
end subroutine MatMul

subroutine Inverse(X,XI,size)
dimension X(100,100),XI(100,100)
integer rowindex,colindex,size,tmpi,tmpj,i,j
real pivotelement

do i=1,size
    do j=1,size
        if(i.EQ.j) then
            XI(i,j)=1
        endif
    end do
end do


do rowindex=1,size !loop over rows of X

    pivotelement=X(rowindex,rowindex) !this is crucial!

    !divide the row by pivot element
    do colindex=1,size
        X(rowindex,colindex)=X(rowindex,colindex)/pivotelement
    end do
    do colindex=1,size
        XI(rowindex,colindex)=XI(rowindex,colindex)/pivotelement
    end do


    !perform row operations to make non-pivot elements of this column zero
    do i=1,size !i gives the row number being manipulated
        pivotelement=X(i,rowindex) !pivot element for the currently manipulated row. this i crucial

        if(i.NE.rowindex) then !row operations are for non-current rows only

            do j=1,size !column index for ith row of X matrix
                X(i,j)=X(i,j)-pivotelement*X(rowindex,j)
            end do
            do j=1,size !column index for ith row of XI matrix
                XI(i,j)=XI(i,j)-pivotelement*XI(rowindex,j)
            end do


        endif
    end do



end do


end subroutine