text:
    use dflib
    !Author : globalpolicy
    !Blog   : c0dew0rth.blogspot.com
    !GitHub : globalpolicy
    !Date   : 06:05PM - 12th August, 2019
    !Purpose: To draw flownet (revised)
    !Changes:
    !-Head at dam heel is water depth as should be(was assumed 0 before)
    !-Dam was fixed to be the size of a single grid box. Now, width of dam can be specified.
    !-Dam's heel and toe are Dirichlet boundaries i.e. known heads but any node between are treated as unknowns.
    !Parameters : contourInterval, numDys, numDxs, DxsTillDam, DamWidthInDxs, reservoirHead, downstreamBoundaryHead

    !Built by observing structure of an example coefficient matrix from my Fortran Simulation notes
    !https://photos.app.goo.gl/BS5C2tp4MNMFp9Dh6 -> example grid
    !https://photos.app.goo.gl/GecBS62BS1hTPNri9 -> matrix construction
    !https://photos.app.goo.gl/WgKvzMWf63CRNBwZA -> final matrix
    !Also, during revision of code, the following problems were taken as a reference:
    !https://photos.app.goo.gl/1p7NWiWeYKTRTBEh6
    !https://photos.app.goo.gl/LF51G7jisALXweCL7
    !Coefficient matrices were constructed for them and solved in Excel:
    !https://drive.google.com/open?id=1I513Eq-tPaADIw_MF5oiTIGo4es4k2f8
    !Great lectures I made my notes off of:
    !https://www.youtube.com/watch?v=FgR3KS51mbA - Prof. Wen Shen, Penn State University
    !https://www.youtube.com/watch?v=i0f2DegLBN8 - Sandip Mazumder

    implicit none

    integer i,j,col,row
    character(10) char
    real tmp,distance,smallestDistanceIndex,smallestDistance
    integer numDys,numDxs,numUnknowns,DxsTillDam,DamWidthInDxs
    integer(2) viewPortTopLeftx,viewPortTopLefty,viewPortBottomRightx,viewPortBottomRighty,viewPortWidth !viewPortWidth=viewPortHeight
    real reservoirHead,downstreamBoundaryHead
    real A(500,500),B(500,500),X(500,500),AInverse(500,500)
    logical status
    TYPE(wxycoord) wxy
    integer ret
    real H(500,500)
    real EquipotentialX(500,500),EquipotentialY(500,500) !one column for one head, same column's row element from each matrix give the equipotential point x,y
    real searchHead,currentHead,rightHead,upHead,contourInterval
    real interpolatedX,interpolatedY
    integer numberOfContours, contourId
    integer contourPointsCounters(500) !stores the number of equipotential points obtained for each contour

    contourInterval=10

    numDys=3 !total grid boxes along Y axis, enclosed within boundaries top and bottom
    numDxs=4 !total grid boxes along X axis, enclosed within boundaries left and right
    DxsTillDam=2
    DamWidthInDxs=2

    reservoirHead=100
    downstreamBoundaryHead=50

    viewPortWidth = 400 !width(=height) of the viewport
    viewPortTopLeftx = 200
    viewPortTopLefty = 80

    numUnknowns=numDys*numDxs+(DamWidthInDxs-1)
    numberOfContours=reservoirHead/contourInterval + 1

    !generate coefficient matrix A
    !the rows, columns and nodes in the comments are referencing the full domain grid including the boundaries(https://photos.app.goo.gl/nToWZ9p3znyGxvwy8), not the array/matrix used here
    do i=1,numUnknowns !row index
        do j=1,numUnknowns !column index
            if(i.EQ.j) then

                A(i,j)=-4

                !for immediate right element of diagonal's row (right blob of Laplace computational stencil)
                if(i.LT.numDxs*numDys.AND.MODULO(i,numDxs).EQ.1) then !if current node is on u/s boundary
                    A(i,j+1)=2
                elseif(i.LE.numDxs*numDys.AND.MODULO(i,numDxs).EQ.0) then !if current node is on second-last column (before d/s boundary) and below the surface row
                    A(i,j+1)=0
                    B(i,1)=B(i,1)-downstreamBoundaryHead
                elseif(i.GT.numDxs*numDys.AND.i.EQ.numUnknowns) then !if current node is the final unknown i.e. surface node right beneath dam, right before toe of dam(d/s boundary)
                    A(i,j+1)=0
                    B(i,1)=B(i,1)-0
                else
                    A(i,j+1)=1
                endif

                !for immediate left element of diagonal's row (left blob of Laplace computational stencil)
                if(j.GT.1) then !if current node is not the bottom-left one
                    if(i.LT.numDxs*numDys.AND.MODULO(i,numDxs).EQ.1) then !if current node is on u/s boundary
                        A(i,j-1)=0
                    elseif(DamWidthInDxs.GE.2.AND.i.EQ.numDxs*numDys+1) then !if current node is right beneath the dam, one node to right of heel
                        A(i,j-1)=0
                        B(i,1)=B(i,1)-reservoirHead
                    else
                        A(i,j-1)=1
                    endif
                endif

                !for far right element of diagonal's row (top blob of Laplace computational stencil)
                if(i.LE.numDxs) then !if current node is in the bottom row
                    A(i,j+numDxs)=2
                elseif(j+numDxs.LE.numDxs*numDys) then !if current row is below the second-top row
                    A(i,j+numDxs)=1
                elseif(j.LE.numDxs*numDys) then !if current row is second-top row(row immediately below surface)
                    if(j.LE.numUnknowns-(DamWidthInDxs-1+numDxs-DxsTillDam-1)) then !if this node is u/s of heel of dam, heel inclusive
                        B(i,1)=B(i,1)+(-reservoirHead)
                    elseif(j.GT.numDxs*numDys-(DamWidthInDxs-1).AND.j.LE.numDxs*numDys) then !if node is between(exclusive) heel and toe
                        A(i,j+(DamWidthInDxs-1))=1
                    endif
                else !if current row is topmost row(unknown nodes directly beneath the dam)
                    B(i,1)=B(i,1)+0
                endif

                !for far left element of diagonal's row (bottom blob of Laplace computational stencil)
                if(j.GT.numDxs) then
                    if(j.LE.numDxs*numDys) then
                        A(i,j-numDxs)=1
                    else
                        A(i,j-(DamWIdthInDxs-1))=2
                    endif
                endif

            endif
        end do
    end do

    !generated coefficient matrix A written to file for inspection
    open(1,file='A.txt')
    do i=1,numUnknowns
        write(1,'(F5.1,$)') (A(i,j),j=1,numUnknowns)
        write(1,*)
    end do
    close(1)

    !generated B vector written to file for inspection
    open(1,file='B.txt')
    do i=1,numUnknowns
        write(1,*) B(i,1)
    end do
    close(1)

    call Inverse(A,AInverse,numUnknowns)
    X=matmul(AInverse,B)

    !calculated X vector
    open(1,file='X.txt')
    do i=1,numUnknowns
        write(1,*) X(i,1)
    end do
    close(1)

    !setup viewport
    viewPortBottomRightx = int2(viewPortTopLeftx+real(numDxs)/real(numDys)*viewPortWidth)
    viewPortBottomRighty = int2(viewPortTopLefty+viewPortWidth)
    call SETVIEWPORT(200,80,viewPortBottomRightx,viewPortBottomRighty)
    status = SETCOLOR(4) !red
    status = FLOODFILL(51,81,4) !flood fill viewport
    status = SETWINDOW(.TRUE.,real(0,8),real(0,8),real(numDxs,8),real(numDys,8))
    status = SETCOLOR(1)

    !draw horizontal lines from bottom to top
    do i=0,numDys
        call MOVETO_W(0.0_8,real(i,8),wxy)
        ret = LINETO_W(real(numDxs,8),real(i,8))
    end do

    !draw vertical lines from left to right
    do i=0,numDxs
        call MOVETO_W(real(i,8),0.0_8,wxy)
        ret = LINETO_W(real(i,8),real(numDys,8))
    end do

    !create Head matrix for whole domain from calculated heads at nodes and the boundary condition values
    !populate non-last row(below surface) of Head matrix from the computed X vector
    do i=1,numDys
        do j=1,numDxs+1
            if(j.EQ.numDxs+1) then !if final column
                H(i,j)=downstreamBoundaryHead
            else
                H(i,j)=X((i-1)*numDxs+j,1)
            endif
        end do
    end do
    !populate the last row (corresponding to the surface heads)
    i=1
    do j=1,numDxs+1
        if(j.LE.DxsTillDam+1) then
            H(numDys+1,j)=reservoirHead
        elseif(j.GT.DxsTillDam+1.AND.j.LT.numDxs+1) then
            H(numDys+1,j)=X(numDxs*numDys+i,1)
            i=i+1
        else
            H(numDys+1,j)=0
        endif
    end do

    !output H matrix to file for inspection
    !observe that it is of the form (taking the referenced example at the top)
    !           h1           h2             h3          h4    downstreamBoundaryHead
    !           h5           h6             h7          h8    downstreamBoundaryHead
    !           h9          h10             h11         h12   downstreamBoundaryHead
    !   reservoirHead   reservoirHead  reservoirHead    h13     0
    !where h1,h2,... are computed heads from X vector
    !also observe that H matrix has exactly the same number of elements as nodes in our whole domain grid(i.e. numDys+1 by numDxs+1)
    open(1,file='H.txt')
    do i=1,numDys+1
        write(1,'(F5.1,$,"  ")') (H(i,j),j=1,numDxs+1)
        write(1,*)
    end do
    close(1)

    contourId=0
    do searchHead=0,reservoirHead,contourInterval
        contourId=contourId+1
        !iterate through grid nodes one column at a time, left to right, looking for same potential(searchHead) right and top
        !note that correspondence between our grid coordinate and H matrix's element is
        !Head at coordinate (col,row) of grid  = (row+1,col+1) element of H matrix
        !Makes use of the fact that in a row, head decreases from left to right.
        !Depending on column on grid, in a column, head decreases top to bottom for col<damX and increases for col>=damX
        do col=0,numDxs
10            do row=0,numDys
                currentHead=H(row+1,col+1)
                if(currentHead.NE.searchHead) then
                    if(col<numDxs) then !look toward east node
                        rightHead=H(row+1,col+2)
                        if(searchHead.LE.currentHead.AND.searchHead.GE.rightHead) then
                            interpolatedX=col+1/(rightHead-currentHead)*(searchHead-currentHead)
                            contourPointsCounters(contourId)=contourPointsCounters(contourId)+1
                            EquipotentialX(contourPointsCounters(contourId),contourId)=interpolatedX
                            EquipotentialY(contourPointsCounters(contourId),contourId)=row
                        endif

                    endif
                    if(row<numDys) then !look toward north node
                        upHead=H(row+2,col+1)
                        if(col<DxsTillDam) then !if we're left of dam
                            !head increases towards up
                            if(searchHead.GE.currentHead.AND.searchHead.LE.upHead) then
                                interpolatedY=row+1/(upHead-currentHead)*(searchHead-currentHead)
                                contourPointsCounters(contourId)=contourPointsCounters(contourId)+1
                                EquipotentialY(contourPointsCounters(contourId),contourId)=interpolatedY
                                EquipotentialX(contourPointsCounters(contourId),contourId)=col
                            endif
                        else !if we're at or right of dam's heel
                            !head decreases towards up
                            if(searchHead.LE.currentHead.AND.searchHead.GE.upHead) then
                                interpolatedY=row+1/(upHead-currentHead)*(searchHead-currentHead)
                                contourPointsCounters(contourId)=contourPointsCounters(contourId)+1
                                EquipotentialY(contourPointsCounters(contourId),contourId)=interpolatedY
                                EquipotentialX(contourPointsCounters(contourId),contourId)=col
                            endif
                        endif
                    endif
                else
                    !we're standing at our searchHead!!
                    contourPointsCounters(contourId)=contourPointsCounters(contourId)+1
                    EquipotentialY(contourPointsCounters(contourId),contourId)=row
                    EquipotentialX(contourPointsCounters(contourId),contourId)=col
                endif
            end do
        end do
    end do

    !!----DEPRECATED EQUIPOTENTIAL POINT SORTING METHOD----
    !!sort equipotential points in ascending (descending works too) order of x coordinate for proper equipotential line
    !!sketching (left to right in ascending, right to left in descending).
    !!EDIT : Equipotential curves that are improper functions (two y's for a x) don't behave nicely with
    !!this method though. Only a couple of curves, generally towards the center of the flownet are as such
    !!in my experience.
    !do contourId=1,numberOfContours
    !    do i=1,contourPointsCounters(contourId)
    !        do j=i,contourPointsCounters(contourId)
    !            if(EquipotentialX(j,contourId).LT.EquipotentialX(i,contourId)) then
    !                tmp=EquipotentialY(i,contourId)
    !                EquipotentialY(i,contourId)=EquipotentialY(j,contourId)
    !                EquipotentialY(j,contourId)=tmp
    !
    !                tmp=EquipotentialX(i,contourId)
    !                EquipotentialX(i,contourId)=EquipotentialX(j,contourId)
    !                EquipotentialX(j,contourId)=tmp
    !            endif
    !        end do
    !    end do
    !end do

    !equipotential points sorting for proper equipotential line sketching
    !first sort equipotential points in descending order of Y coordinate(topmost to bottom-most)
    !then sort based on shortest consecutive distance.
    !better than sorting by x coordinate but not totally devoid of some weird artifacts
    do contourId=1,numberOfContours
        !sort descending in Y coordinate
        do i=1,contourPointsCounters(contourId)
            do j=i,contourPointsCounters(contourId)
                if(EquipotentialY(j,contourId).GT.EquipotentialY(i,contourId)) then
                    !swap Y coordinate
                    tmp=EquipotentialY(i,contourId)
                    EquipotentialY(i,contourId)=EquipotentialY(j,contourId)
                    EquipotentialY(j,contourId)=tmp
                    !swap X coordinate
                    tmp=EquipotentialX(i,contourId)
                    EquipotentialX(i,contourId)=EquipotentialX(j,contourId)
                    EquipotentialX(j,contourId)=tmp
                endif
            end do
        end do
        !sort based on shortest consecutive distance (basically ordering the shortest path downward, downward coz the array is descending inY)
        do i=1,contourPointsCounters(contourId)-1
            smallestDistance=10000 !very large number
            do j=i+1,contourPointsCounters(contourId)
                !distance between i and j point
                distance = (EquipotentialX(i,contourId)-EquipotentialX(j,contourId))**2 + (EquipotentialY(i,contourId)-EquipotentialY(j,contourId))**2
                if(distance.LT.smallestDistance) then
                    smallestDistance=distance
                    smallestDistanceIndex=j
                endif
            end do
            !swap i+1 point and point at smallestDistanceIndex
            !X
            tmp = EquipotentialX(i+1,contourId)
            EquipotentialX(i+1,contourId) = EquipotentialX(smallestDistanceIndex,contourId)
            EquipotentialX(smallestDistanceIndex,contourId) = tmp
            !Y
            tmp = EquipotentialY(i+1,contourId)
            EquipotentialY(i+1,contourId) = EquipotentialY(smallestDistanceIndex,contourId)
            EquipotentialY(smallestDistanceIndex,contourId) = tmp
        end do
    end do

    !write equipotential points to file for observation
    open(1,file='Isobars.txt')
        do contourId=1,numberOfContours
            write(1,*) (contourId-1)*contourInterval
            write(1,'("(",F5.2,",",F5.2,")      ",$)') (EquipotentialX(i,contourId),EquipotentialY(i,contourId),i=1,contourPointsCounters(contourId))
            write(1,*) ''
        end do
    close(1)

    !draw equipotential lines
    do contourId=1,numberOfContours
        call MOVETO_W(real(EquipotentialX(1,contourId),8),real(EquipotentialY(1,contourId),8),wxy)

        !draw a circle on the equipotential point(not necessary)
        status = SETCOLOR(6)
        ret = ELLIPSE_W($GFILLINTERIOR,real(EquipotentialX(1,contourId)-0.015,8),real(EquipotentialY(1,contourId)+0.015,8),real(EquipotentialX(1,contourId)+0.015,8),real(EquipotentialY(1,contourId)-0.015,8))
        status = SETCOLOR(2)

        do i=2,contourPointsCounters(contourId)
            !draw a circle on the equipotential point(not necessary)
            status = SETCOLOR(6)
            ret = ELLIPSE_W($GFILLINTERIOR,real(EquipotentialX(i,contourId)-0.015,8),real(EquipotentialY(i,contourId)+0.015,8),real(EquipotentialX(i,contourId)+0.015,8),real(EquipotentialY(i,contourId)-0.015,8))
            status = SETCOLOR(2)

            !avoids jumping lines
            call GETCURRENTPOSITION_W(wxy)
            if(abs(wxy.wx-EquipotentialX(i,contourId)).GT.4.OR.abs(wxy.wy-EquipotentialY(i,contourId)).GT.4) then
                call MOVETO_W(real(EquipotentialX(i,contourId),8),real(EquipotentialY(i,contourId),8),wxy)
            endif

            ret = LINETO_W(real(EquipotentialX(i,contourId),8),real(EquipotentialY(i,contourId),8))

        end do
    end do

    !draw one grid length horizontal line on grid top representing base of dam
    status = SETCOLOR(8)
    call MOVETO_W(real(DxsTillDam,8),real(numDys,8),wxy)
    ret = LINETO_W(real(DxsTillDam+DamWidthInDxs,8),real(numDys,8))

    !write heads on grid nodes
    !note that correspondence between our grid coordinate and H matrix's element is
    !head at coordinate (col,row) of grid  = (row+1,col+1) element of H matrix
    status = SETCOLOR(8)
    status = INITIALIZEFONTS()
    status = SETFONT('t''Arial''h12w8pvb')
    !loop for except the bottom and right boundaries and bottom-right node (else they're clipped)
    do i=0,numDxs-1
        do j=1,numDys
            call MOVETO_W(real(i,8),real(j,8),wxy)
            write(char,'(F6.2)') H(j+1,i+1)
            call OUTGTEXT(char)
        end do
    end do
    !for bottom boundary
    do i=0,numDxs-1
        call MOVETO_W(real(i,8),real(0,8)+0.1,wxy)
        write(char,'(F6.2)') H(0+1,i+1)
        call OUTGTEXT(char)
    end do
    !for right boundary
    do i=1,numDys
        call MOVETO_W(real(numDxs,8)-0.3,real(i,8),wxy)
        write(char,'(F6.2)') H(i+1,numDxs+1)
        call OUTGTEXT(char)
    end do
    !for bottom right node
    call MOVETO_W(real(numDxs,8)-0.3,real(0,8)+0.1,wxy)
    write(char,'(F6.2)') H(0+1,numDxs+1)
    call OUTGTEXT(char)

    end program

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

originalX=X

do i=1,size
    do j=1,size
        if(i.EQ.j) then
            XI(i,j)=1
        else
            XI(i,j)=0 !this is for good measure coz by default, XI may have been changed by the caller itself though it is 0 by default
        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

!restore X to original X since upon inversion completion, X has been reduced to an Identity matrix. good practice
X=originalX

end subroutine

max_stars_repo_name: globalpolicy/Flownet

max_issues_repo_name: globalpolicy/Flownet

max_forks_repo_name: globalpolicy/Flownet


text:
great_ne(z1,ii1).
great_ne(ff1,i1).
great_ne(cc1,f1).
great_ne(kk1,w1).
great_ne(jj1,aa1).
great_ne(u1,f1).
great_ne(p1,hh1).
great_ne(ee1,i1).
great_ne(ee1,ii1).
great_ne(cc1,ff1).
great_ne(jj1,k1).
great_ne(k1,x1).
great_ne(bb1,m1).
great_ne(ff1,c1).
great_ne(bb1,k1).
great_ne(p1,a1).
great_ne(hh1,m1).
great_ne(m1,y1).
great_ne(kk1,y1).
great_ne(bb1,b1).
great_ne(k1,c1).
great_ne(hh1,dd1).
great_ne(p1,h1).
great_ne(kk1,dd1).
great_ne(u1,bb1).
great_ne(u1,c1).
great_ne(aa1,bb1).
great_ne(p1,ff1).
great_ne(aa1,hh1).
great_ne(k1,b1).
great_ne(aa1,k1).
great_ne(q1,m1).
great_ne(t1,ee1).
great_ne(cc1,k1).
great_ne(q1,x1).
great_ne(q1,aa1).
great_ne(v1,b1).
great_ne(jj1,ee1).
great_ne(dd1,v1).
great_ne(jj1,ff1).
great_ne(bb1,i1).
great_ne(kk1,b1).
great_ne(p1,dd1).
great_ne(t1,u1).
great_ne(z1,ee1).
great_ne(z1,hh1).
great_ne(t1,dd1).
great_ne(x1,v1).
great_ne(x1,ii1).
great_ne(aa1,y1).
great_ne(v1,w1).
great_ne(w1,a1).
great_ne(cc1,u1).
great_ne(q1,ii1).
great_ne(p1,b1).
great_ne(ii1,f1).
great_ne(aa1,f1).
great_ne(ii1,c1).
great_ne(t1,hh1).
great_ne(jj1,f1).
great_ne(ee1,hh1).
great_ne(cc1,jj1).
great_ne(u1,x1).
great_ne(q1,i1).
great_ne(bb1,c1).
great_ne(cc1,hh1).
great_ne(z1,h1).
great_ne(hh1,i1).
great_ne(hh1,f1).
great_ne(t1,jj1).
great_ne(m1,i1).
great_ne(jj1,m1).
great_ne(p1,i1).
great_ne(cc1,w1).
great_ne(z1,bb1).
great_ne(aa1,w1).
great_ne(jj1,dd1).
great_ne(u1,jj1).
great_ne(v1,h1).
great_ne(u1,b1).
great_ne(cc1,bb1).
great_ne(dd1,c1).
great_ne(p1,x1).
great_ne(t1,a1).
great_ne(bb1,x1).
great_ne(u1,y1).
great_ne(k1,a1).
great_ne(ee1,k1).
great_ne(q1,bb1).
great_ne(z1,t1).
great_ne(jj1,w1).
great_ne(dd1,m1).
great_ne(cc1,b1).
great_ne(bb1,a1).
great_ne(hh1,bb1).
great_ne(hh1,c1).
great_ne(p1,kk1).
great_ne(z1,c1).
great_ne(b1,ff1).
great_ne(v1,y1).
great_ne(jj1,bb1).
great_ne(q1,ee1).
great_ne(cc1,a1).
great_ne(z1,w1).
great_ne(kk1,f1).
great_ne(t1,ff1).
great_ne(hh1,ff1).
great_ne(b1,i1).
great_ne(w1,ff1).
great_ne(t1,k1).
great_ne(f1,i1).
great_ne(z1,cc1).
great_ne(y1,b1).
great_ne(w1,f1).
great_ne(aa1,ee1).
great_ne(dd1,bb1).
great_ne(f1,a1).
great_ne(kk1,v1).
great_ne(x1,i1).
great_ne(x1,h1).
great_ne(jj1,ii1).
great_ne(u1,kk1).
great_ne(dd1,f1).
great_ne(aa1,h1).
great_ne(k1,f1).
great_ne(ee1,bb1).
great_ne(ii1,y1).
great_ne(m1,c1).
great_ne(m1,ii1).
great_ne(v1,f1).
great_ne(dd1,ff1).
great_ne(cc1,i1).
great_ne(q1,dd1).
great_ne(u1,m1).
great_ne(w1,i1).
great_ne(cc1,h1).
great_ne(z1,f1).
great_ne(ee1,dd1).
great_ne(aa1,m1).
great_ne(ee1,h1).
great_ne(x1,w1).
great_ne(cc1,kk1).
great_ne(aa1,ii1).
great_ne(q1,k1).
great_ne(u1,dd1).
great_ne(t1,kk1).
great_ne(jj1,a1).
great_ne(u1,v1).
great_ne(p1,jj1).
great_ne(cc1,ii1).
great_ne(hh1,ii1).
great_ne(p1,ee1).
great_ne(q1,y1).
great_ne(y1,w1).
great_ne(aa1,i1).
great_ne(v1,ff1).
great_ne(bb1,h1).
great_ne(ii1,w1).
great_ne(dd1,b1).
great_ne(aa1,v1).
great_ne(p1,v1).
great_ne(u1,ff1).
great_ne(bb1,ff1).
great_ne(jj1,i1).
great_ne(t1,i1).
great_ne(jj1,kk1).
great_ne(q1,b1).
great_ne(v1,c1).
great_ne(f1,c1).
great_ne(b1,a1).
great_ne(u1,w1).
great_ne(jj1,hh1).
great_ne(u1,i1).
great_ne(q1,t1).
great_ne(aa1,dd1).
great_ne(ii1,a1).
great_ne(y1,i1).
great_ne(p1,bb1).
great_ne(aa1,kk1).
great_ne(dd1,ii1).
great_ne(t1,v1).
great_ne(z1,jj1).
great_ne(cc1,c1).
great_ne(z1,v1).
great_ne(z1,x1).
great_ne(w1,h1).
great_ne(dd1,y1).
great_ne(jj1,b1).
great_ne(z1,aa1).
great_ne(kk1,ii1).
great_ne(q1,v1).
great_ne(cc1,x1).
great_ne(cc1,y1).
great_ne(jj1,y1).
great_ne(y1,a1).
great_ne(m1,v1).
great_ne(z1,y1).
great_ne(q1,ff1).
great_ne(w1,b1).
great_ne(q1,f1).
great_ne(u1,h1).
great_ne(k1,w1).
great_ne(z1,a1).
great_ne(hh1,b1).
great_ne(q1,p1).
great_ne(f1,h1).
great_ne(v1,a1).
great_ne(dd1,a1).
great_ne(z1,b1).
great_ne(hh1,a1).
great_ne(t1,y1).
great_ne(cc1,t1).
great_ne(b1,h1).
great_ne(m1,w1).
great_ne(p1,c1).
great_ne(ee1,f1).
great_ne(x1,c1).
great_ne(hh1,w1).
great_ne(q1,c1).
great_ne(ee1,c1).
great_ne(kk1,c1).
great_ne(v1,i1).
great_ne(m1,a1).
great_ne(ee1,y1).
great_ne(dd1,k1).
great_ne(p1,u1).
great_ne(q1,kk1).
great_ne(jj1,c1).
great_ne(y1,f1).
great_ne(ff1,h1).
great_ne(t1,m1).
great_ne(dd1,w1).
great_ne(t1,w1).
great_ne(dd1,h1).
great_ne(ff1,b1).
great_ne(b1,bb1).
great_ne(i1,dd1).
great_ne(u1,t1).
great_ne(ii1,u1).
great_ne(dd1,ee1).
great_ne(ee1,q1).
great_ne(dd1,z1).
great_ne(w1,jj1).
great_ne(ii1,hh1).
great_ne(x1,u1).
great_ne(f1,ee1).
great_ne(a1,t1).
great_ne(i1,cc1).
great_ne(v1,ee1).
great_ne(h1,w1).
great_ne(ff1,q1).
great_ne(dd1,kk1).
great_ne(k1,t1).
great_ne(i1,jj1).
great_ne(b1,y1).
great_ne(h1,u1).
great_ne(h1,ee1).
great_ne(x1,dd1).
great_ne(b1,z1).
great_ne(bb1,u1).
great_ne(b1,ee1).
great_ne(ff1,cc1).
great_ne(bb1,hh1).
great_ne(x1,p1).
great_ne(u1,z1).
great_ne(k1,q1).
great_ne(ff1,ee1).
great_ne(h1,q1).
great_ne(a1,f1).
great_ne(jj1,p1).
great_ne(hh1,jj1).
great_ne(b1,jj1).
great_ne(m1,kk1).
great_ne(c1,w1).
great_ne(h1,ii1).
great_ne(y1,cc1).
great_ne(aa1,jj1).
great_ne(y1,v1).
great_ne(hh1,u1).
great_ne(ff1,f1).
great_ne(bb1,aa1).
great_ne(ee1,jj1).
great_ne(h1,ff1).
great_ne(ii1,p1).
great_ne(y1,aa1).
great_ne(ff1,u1).
great_ne(hh1,kk1).
great_ne(m1,u1).
great_ne(b1,hh1).
great_ne(y1,kk1).
great_ne(v1,p1).
great_ne(x1,hh1).
great_ne(i1,bb1).
great_ne(f1,ii1).
great_ne(w1,p1).
great_ne(i1,y1).
great_ne(c1,k1).
great_ne(w1,x1).
great_ne(f1,t1).
great_ne(v1,q1).
great_ne(i1,x1).
great_ne(w1,t1).
great_ne(y1,t1).
great_ne(h1,m1).
great_ne(a1,z1).
great_ne(t1,p1).
great_ne(dd1,q1).
great_ne(h1,bb1).
great_ne(bb1,jj1).
great_ne(a1,cc1).
great_ne(aa1,p1).
great_ne(v1,t1).
great_ne(a1,hh1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
module mod_dft
#include <messenger.h>
    use mod_kinds,      only: ik, rk
    use mod_constants,  only: PI, ZERO, ONE, TWO, DIR_2, &
                              XI_MIN, XI_MAX, ETA_MIN, ETA_MAX, ZETA_MIN, ZETA_MAX
    use type_domain,    only: domain_t
    use type_point,     only: point_t
    use DNAD_D
    implicit none

contains

    !>  Compute the points for a discrete Fourier transform across a boundary.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/17/2016
    !!
    !!
    !!
    !----------------------------------------------------------------------------------
    function compute_dft_points(domain,elems,faces,periodicity) result(points)
        type(domain_t), intent(in)   :: domain
        integer(ik),    intent(in)   :: elems(:)
        integer(ik),    intent(in)   :: faces(:)
        real(rk),       intent(in)   :: periodicity

        type(point_t),  allocatable  :: points(:)
        integer(ik)                  :: nmodes, npoints, ierr, ielem_bc, ielem, var, mode, ipnt, min_y_element, elem_min, iface
        integer(ik)                  :: itime
        integer(ik)                  :: min_y_loc
        real(rk)                     :: xi, eta, zeta, min_y, dy, xloc, zloc

        real(rk),       allocatable  :: mean_y_coordinates(:)
        integer(ik),    allocatable  :: mean_y_elements(:)
        real(rk),       dimension(2) :: ylocs, xlocs, zlocs

        !
        ! TODO: fix hard coded face orientation
        !
        iface = faces(1)

        !
        ! Number of Fourier harmonics to compute
        !
        nmodes = 10

        !
        ! Compute number of points in dft
        !
        npoints = 1 + 2*nmodes

        !
        ! Compute distance between points
        !
        dy = periodicity/real((npoints), rk)

        !
        ! Allocate points
        !
        allocate(points(npoints), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Allocate storage for mean coordinates
        !
        allocate(mean_y_coordinates(size(elems)), &
                 mean_y_elements(   size(elems)), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Coordinate min/max on the boundary in direction of periodicity.
        !
        do ielem_bc = 1,size(elems)

            !
            ! Get element index in block domain
            !
            ielem = elems(ielem_bc)

            !
            ! Get mean y-coordinate from element. y-coordinate variable = 2. Mode = 1
            !
            var  = DIR_2    ! y-coordinate index
            mode = 1        ! Coordinate average
            mean_y_coordinates(ielem_bc) = domain%elems(ielem)%coords%getterm(var,mode,itime)
            mean_y_elements(ielem_bc)    = ielem

        end do ! ielem

        !
        ! Find element with minimum mean y-coordinate
        !
        min_y_element = minloc(mean_y_coordinates, DIM=1)

        elem_min = mean_y_elements(min_y_element)

        !
        ! For the element with the minimum mean y-coordinate, find minimum y-coordinate
        ! on the boundary face. So, for the time being, XI_MIN, ETA_MIN, ZETA=0
        !
!        !xi   = -ONE ! INLET
!        xi   = -ONE  ! OUTLET
!        eta  = -ONE
!        zeta = domain%elems(elem_min)%coords%getterm(3,1) ! get mean z-coordinate
!
!        min_y = domain%elems(elem_min)%y(xi,eta,zeta) + 0.0000001_rk

        if ( (iface == XI_MIN) ) then
            ylocs(1) = domain%elems(elem_min)%y(-ONE, ONE,ZERO)
            ylocs(2) = domain%elems(elem_min)%y(-ONE,-ONE,ZERO)

            xlocs(1) = domain%elems(elem_min)%x(-ONE, ONE,ZERO)
            xlocs(2) = domain%elems(elem_min)%x(-ONE,-ONE,ZERO)

            zlocs(1) = domain%elems(elem_min)%z(-ONE, ONE,ZERO)
            zlocs(2) = domain%elems(elem_min)%z(-ONE,-ONE,ZERO)

            min_y_loc = minloc(ylocs,1)
            min_y     = ylocs(min_y_loc) + 0.000001_rk
            xloc      = xlocs(min_y_loc)
            zloc      = zlocs(min_y_loc)

        else if ( (iface == XI_MAX) ) then
            ylocs(1) = domain%elems(elem_min)%y(ONE, ONE,ZERO)
            ylocs(2) = domain%elems(elem_min)%y(ONE,-ONE,ZERO)

            xlocs(1) = domain%elems(elem_min)%x(ONE, ONE,ZERO)
            xlocs(2) = domain%elems(elem_min)%x(ONE,-ONE,ZERO)

            zlocs(1) = domain%elems(elem_min)%z(ONE, ONE,ZERO)
            zlocs(2) = domain%elems(elem_min)%z(ONE,-ONE,ZERO)

            min_y_loc = minloc(ylocs,1)
            min_y     = ylocs(min_y_loc) + 0.000001_rk
            xloc      = xlocs(min_y_loc)
            zloc      = zlocs(min_y_loc)

        else if ( (iface == ETA_MIN) ) then
            ylocs(1) = domain%elems(elem_min)%y( ONE,-ONE,ZERO)
            ylocs(2) = domain%elems(elem_min)%y(-ONE,-ONE,ZERO)

            xlocs(1) = domain%elems(elem_min)%x( ONE,-ONE,ZERO)
            xlocs(2) = domain%elems(elem_min)%x(-ONE,-ONE,ZERO)

            zlocs(1) = domain%elems(elem_min)%z( ONE,-ONE,ZERO)
            zlocs(2) = domain%elems(elem_min)%z(-ONE,-ONE,ZERO)

            min_y_loc = minloc(ylocs,1)
            min_y     = ylocs(min_y_loc) + 0.000001_rk
            xloc      = xlocs(min_y_loc)
            zloc      = zlocs(min_y_loc)

        else if ( (iface == ETA_MAX) ) then
            ylocs(1) = domain%elems(elem_min)%y( ONE,ONE,ZERO)
            ylocs(2) = domain%elems(elem_min)%y(-ONE,ONE,ZERO)

            xlocs(1) = domain%elems(elem_min)%x( ONE,ONE,ZERO)
            xlocs(2) = domain%elems(elem_min)%x(-ONE,ONE,ZERO)

            zlocs(1) = domain%elems(elem_min)%z( ONE,ONE,ZERO)
            zlocs(2) = domain%elems(elem_min)%z(-ONE,ONE,ZERO)

            min_y_loc = minloc(ylocs,1)
            min_y     = ylocs(min_y_loc) + 0.000001_rk
            xloc      = xlocs(min_y_loc)
            zloc      = zlocs(min_y_loc)

        else if ( (iface == ZETA_MIN) .or. (iface == ZETA_MAX) ) then

            call chidg_signal(FATAL,"compute_dft_ponts: assumes boundary is on a XI or ETA face.")

        end if

        !
        ! Create array of points from (min_y) to (min_y + periodicity)
        !
        !call points(1)%set(xloc, min_y, ZERO)
        call points(1)%set(xloc, min_y, zloc)
        do ipnt = 2,npoints

            !
            ! Copy previous point
            !
            points(ipnt) = points(ipnt - 1)

            !
            ! Add dy
            !
            call points(ipnt)%add_y(dy)

        end do ! ipnt

    end function compute_dft_points
    !**********************************************************************************

    !>  Computes the discrete Fourier transform of an input data set.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes X[i] = SUM[ x(i)e^(-i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine dft(input,outreal,outimag)
        type(AD_D), intent(in)                 :: input(:)
        type(AD_D), intent(inout), allocatable :: outreal(:)
        type(AD_D), intent(inout), allocatable :: outimag(:)

        type(AD_D) :: sreal, simag, angle
        integer      :: n, k, t

        n = size(input)

        !
        ! Initialize derivative arrays
        !
        outreal = input
        outimag = input
        sreal   = input(1)
        simag   = input(1)
        angle   = input(1)

        outreal = ZERO
        outimag = ZERO
        sreal   = ZERO
        simag   = ZERO
        angle   = ZERO

        !
        ! Loop over output modes
        !
        !do k = 1,n+1
        do k = 1,n

            !
            ! Get contribution from input points
            !
            sreal = ZERO
            simag = ZERO
            do t = 1,n

                angle = TWO * PI * real(t-1,rk) * real(k-1,rk) / real(n,rk)

                sreal = sreal + input(t) * cos(-angle)
                simag = simag + input(t) * sin(-angle)

            end do

            outreal(k) = sreal
            outimag(k) = simag

        end do ! k

        outreal = outreal / real(n,rk)
        outimag = outimag / real(n,rk)

    end subroutine dft
    !******************************************************************

    !>  Computes the inverse discrete Fourier transform of an input data set.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes x[i] = (1/N) SUM[ X(i)e^( i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine idft(inreal,inimag,output)
        type(AD_D),   intent(in)                   :: inreal(:)
        type(AD_D),   intent(in)                   :: inimag(:)
        type(AD_D),   intent(inout), allocatable   :: output(:)

        type(AD_D)    :: s
        real(rk)       :: angle
        integer       :: k, t, nmodes, npoints

        nmodes   = size(inreal)
        npoints  = nmodes

        !
        ! Initialize derivative arrays
        !
        output = inreal
        s      = inreal(1)

        output = ZERO
        s      = ZERO

        !
        ! Loop over output points
        !
        do k = 1,npoints

            !
            ! Loop over input modes
            !
            output(k) = ZERO
            s         = ZERO
            do t = 1,nmodes

                angle = TWO * PI * real(t-1,rk) * real(k-1,rk) / real(nmodes,rk)

                s = s + inreal(t)*cos(angle) - inimag(t)*sin(angle)

            end do

            output(k) = s

        end do

        !
        ! Normalization
        !
        !output = output / real(nmodes,8)

    end subroutine idft
    !******************************************************************

    !>  Computes the inverse discrete Fourier transform of an input data at
    !!  user specified angles.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes x[i] = (1/N) SUM[ X(i)e^( i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine idft_angle(inreal, inimag, angles, output)
        type(AD_D),       intent(in)      :: inreal(:)
        type(AD_D),       intent(in)      :: inimag(:)
        type(AD_D),       intent(in)      :: angles(:)
        type(AD_D),       intent(inout), allocatable   :: output(:)

        type(AD_D)            :: sreal, simag, s, freq, mag, phase
        integer               :: iang, imode, nmodes, nangles, quadrant

        nangles  = size(angles)
        nmodes   = size(inreal)

        !
        ! Initialize derivative arrays
        !
        output = angles
        sreal  = inreal(1)
        simag  = inreal(1)
        s      = inreal(1)
        freq   = inreal(1)

        output = ZERO
        sreal  = ZERO
        simag  = ZERO
        s      = ZERO
        freq   = ZERO

        !
        ! Loop over angles
        !
        do iang = 1,nangles

            !
            ! Accumulate contribution from each mode.
            !
            output(iang) = ZERO
            s            = ZERO
            do imode = 1,(nmodes-1)/2

                !
                ! Compute current frequency
                !
                freq  = TWO * PI * real(imode-1,rk)

!                !
!                ! Compute magnitude and phase from the complex mode
!                !
!                mag   = sqrt( inreal(imode)**TWO + inimag(imode)**TWO )
!                phase = asin( inimag(imode) / mag )
!
!                !
!                ! Check quadrant
!                !
!                if      ( (inreal(imode) >= ZERO) .and. (inimag(imode) >= ZERO) ) then
!                    quadrant = 1
!                else if ( (inreal(imode) < ZERO)  .and. (inimag(imode) >= ZERO) ) then
!                    quadrant = 2
!                else if ( (inreal(imode) < ZERO)  .and. (inimag(imode) < ZERO) ) then
!                    quadrant = 3
!                else if ( (inreal(imode) >= ZERO)  .and. (inimag(imode) < ZERO) ) then
!                    quadrant = 4
!                end if
!
!
!
!
!                !
!                ! Account for quadrant
!                !
!                if      ( quadrant == 1 ) then
!                    phase = phase
!                else if ( quadrant == 2) then
!                    phase = PI - phase
!                else if ( quadrant == 3) then
!                    phase = PI - phase
!                else if ( quadrant == 4) then
!                    phase = phase
!                end if
!
!
!
!
!                if ( imode > 1 ) then
!                    s = s +  TWO * mag*cos( freq * angles(iang) + phase )
!                else
!                    s = s +  mag*cos( freq * angles(iang) + phase )
!                end if

                if ( imode > 1 ) then
                    s = s + TWO * (inreal(imode)*cos(freq*angles(iang)) - inimag(imode)*sin(freq*angles(iang)))
                else
                    s = s + inreal(imode)*cos(freq*angles(iang)) - inimag(imode)*sin(freq*angles(iang))
                end if

            end do ! accumulate modes

            !
            ! Store value for point k.
            !
            output(iang) = s

        end do  ! evaluate points

        !
        ! Normalization
        !
        !output = output / real(nmodes,8)

    end subroutine idft_angle
    !******************************************************************

    !>  Computes the inverse discrete Fourier transform of an input data at
    !!  user specified angles.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes x[i] = (1/N) SUM[ X(i)e^( i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine idft_mode_angle(inreal, inimag, imode, angles, output)
        type(AD_D),       intent(in)      :: inreal
        type(AD_D),       intent(in)      :: inimag
        integer(ik),      intent(in)      :: imode
        real(rk),         intent(in)      :: angles(:)
        type(AD_D),       intent(inout)   :: output(:)

        type(AD_D)      :: sreal, simag, s, mag, phase
        real(rk)        :: freq
        integer         :: iang, nangles, quadrant

        nangles  = size(angles)

        !
        ! Initialize derivative arrays
        !
!        output = angles
        sreal  = inreal
        simag  = inreal
        s      = inreal
!        freq   = inreal

        output = ZERO
        sreal  = ZERO
        simag  = ZERO
        s      = ZERO
        freq   = ZERO

        !
        ! Loop over angles
        !
        do iang = 1,nangles

            !
            ! Accumulate contribution from each mode.
            !
            output(iang) = ZERO
            s            = ZERO

            !
            ! Compute current frequency
            !
            freq  = TWO * PI * real(imode-1,rk)

!            !
!            ! Compute magnitude and phase from the complex mode
!            !
!            print*, 'idft_mode_angle - inreal, inimag'
!            print*, inreal%xp_ad_
!            print*, inimag%xp_ad_
!            read*,
!            mag = sqrt( inreal**TWO + inimag**TWO )
!
!
!            if ( mag == ZERO ) then
!                phase = mag*ZERO ! avoid division by zero
!            else
!
!                arg = inimag / mag
!
!                if ( arg == 1 ) then
!
!                phase = asin( inimag / mag )
!            end if
!
!
!
!            !
!            ! Check quadrant
!            !
!            if      ( (inreal >= ZERO)  .and. (inimag >= ZERO) ) then
!                quadrant = 1
!            else if ( (inreal < ZERO)   .and. (inimag >= ZERO) ) then
!                quadrant = 2
!            else if ( (inreal < ZERO)   .and. (inimag < ZERO)  ) then
!                quadrant = 3
!            else if ( (inreal >= ZERO)  .and. (inimag < ZERO)  ) then
!                quadrant = 4
!            end if
!
!
!            !
!            ! Account for quadrant
!            !
!            if      ( quadrant == 1 ) then
!                phase = phase
!            else if ( quadrant == 2) then
!                phase = PI - phase
!            else if ( quadrant == 3) then
!                phase = PI - phase
!            else if ( quadrant == 4) then
!                phase = phase
!            end if
!
!
!
!            print*, 'idft_mode_angle - s before'
!            print*, s%xp_ad_
!            print*, mag%xp_ad_
!            print*, phase%xp_ad_
!
!
!            if ( imode > 1 ) then
!                s = s +  TWO * mag*cos( freq * angles(iang) + phase )
!            else
!                s = s +  mag*cos( freq * angles(iang) + phase ) ! DC offset
!            end if
!
!            print*, 'idft_mode_angle - s after'
!            print*, s%xp_ad_
!            read*,

            if ( imode > 1 ) then
                s = s + TWO * (inreal*cos(freq*angles(iang)) - inimag*sin(freq*angles(iang)))
            else
                s = s + inreal*cos(freq*angles(iang)) - inimag*sin(freq*angles(iang))
            end if

            !
            ! Store value for point k.
            !
            output(iang) = s

        end do  ! evaluate points

        !
        ! Normalization
        !
        !output = output / real(nmodes,8)

    end subroutine idft_mode_angle
    !******************************************************************

    !>  Computes the inverse discrete Fourier transform of an input data at
    !!  user specified points along the period, instead of the original evenly
    !!  spaced periodic points.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes x[i] = (1/N) SUM[ X(i)e^( i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine idft_points(ymin, periodicity, inreal, inimag, points, output)
        real(rk),         intent(in)  :: ymin
        real(rk),         intent(in)  :: periodicity
        type(AD_D),       intent(in)  :: inreal(:)
        type(AD_D),       intent(in)  :: inimag(:)
        type(AD_D),       intent(in)  :: points(:)
        type(AD_D),       intent(inout), allocatable  :: output(:)

        type(AD_D)   :: location(size(points))
        integer      :: npts, k

        npts    = size(points)

        !
        ! Compute location of each point in the period: 0-1
        !
        do k = 1,npts
            location(k) = (points(k)-ymin)/periodicity     ! 0-1.   0 = beginning of period.    1 = end of period.
        end do

        !
        ! Evaluate modes at specified locations
        !
        call idft_angle(inreal,inimag,location,output)

    end subroutine idft_points
    !******************************************************************

    !>  Computes the inverse discrete Fourier transform of an input data at
    !!  user specified points along the period, instead of the original evenly
    !!  spaced periodic points.
    !!
    !!  Adapted from Public Domain code:
    !!  www.nayuki.io/res/how-to-implement-the-discrete-fourier-transform/dft.py
    !!
    !!  Computes x[i] = (1/N) SUM[ X(i)e^( i*2pi*n*k*/N) ]
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/12/2016
    !!
    !!
    !-----------------------------------------------------------------
    subroutine idft_mode_points(ymin, periodicity, inreal, inimag, imode, points, output)
        real(rk),         intent(in)    :: ymin
        real(rk),         intent(in)    :: periodicity
        type(AD_D),       intent(in)    :: inreal
        type(AD_D),       intent(in)    :: inimag
        integer(ik),      intent(in)    :: imode
        real(rk),         intent(in)    :: points(:)
        type(AD_D),       intent(inout) :: output(:)

        real(rk)   :: angles(size(points))
        integer      :: npts, k

        npts = size(points)

        !
        ! Compute location of each point in the period: 0-1
        !
        do k = 1,npts
            angles(k) = (points(k)-ymin)/periodicity     ! 0-1.   0 = beginning of period.    1 = end of period.
        end do

        !
        ! Evaluate modes at specified locations
        !
        call idft_mode_angle(inreal,inimag,imode,angles,output)

    end subroutine idft_mode_points
    !******************************************************************

end module mod_dft

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


program main
  implicit none
  call sub_main('./input/HC_1.dat','ene_HC_1_03-06.dat',3,6)
  call sub_main('./input/HC_2.dat','ene_HC_2_06-07.dat',6,7)
  call sub_main('./input/HC_3.dat','ene_HC_3_07-08.dat',7,8)
  call sub_main('./input/K_1.dat','ene_K_1_07-08.dat',7,8)
  call sub_main('./input/K_2.dat','ene_K_2_10-12.dat',10,12)
  call sub_main('./input/K_3.dat','ene_K_3_05-06.dat',5,6)
  call sub_main('./input/K_4.dat','ene_K_4_02-07.dat',2,7)
  call sub_main('./input/K_4.dat','ene_K_4_06-07.dat',6,7)
  call sub_main('./input/K_4.dat','ene_K_4_02-06.dat',2,6)
  call sub_main('./input/S_1.dat','ene_S_1_09-10.dat',9,10)
  call sub_main('./input/S_1.dat','ene_S_1_15-16.dat',15,16)
  call sub_main('./input/S_2.dat','ene_S_2_10-14.dat',10,14)
  call sub_main('./input/S_3.dat','ene_S_3_11-15.dat',11,15)
  call sub_main('./input/S_4.dat','ene_S_4_06-07.dat',6,7)
  call sub_main('./input/S_4.dat','ene_S_4_10-11.dat',10,11)
  call sub_main('./input/S_5.dat','ene_S_5_07-11.dat',7,11)
  call sub_main('./input/S_6.dat','ene_S_6_06-10.dat',6,10)
contains
  subroutine sub_main(ifile,ofile,p1,p2)
    use omp_lib
    use prod_gate
    use set_gate_list
    implicit none
    character(*), intent(in) :: ifile, ofile
    integer, intent(in) :: p1, p2
    integer :: nol, nog, tq, noth
    character(2) :: cnol
    real(8) :: e(3), st, en, mtime
    real(8), allocatable :: ths(:,:)
    integer, allocatable :: list_p(:,:)
    complex(8), allocatable :: v1(:), v2(:)
    complex(8) :: zdotc
    integer :: lv, i
    call set_mat
    mtime = 1.0d9
    call read_2qgate_list(list_p,tq,nol,nog,ifile)
    write(cnol,'(i2.2)') nol
    call read_so4_theta_list(ths,noth,nol,'./input/thetas_'//cnol//'.dat')
    lv = 2**tq
    allocate(v1(lv), v2(lv))
    open(20,file=ofile)
    write(20,*) "nol, nog, noth, tq, lv", nol, nog, noth, tq, lv
    do i = 1, noth
       st = omp_get_wtime()
       call phi_1_i(v1,lv,list_p,nog,ths(:,i),nol)
       call gate2_vec_o(v2,v1,lv,pxx,p1,p2)
       e(1)=dble(zdotc(lv,v1,1,v2,1))
       call gate2_vec_o(v2,v1,lv,pyy,p1,p2)
       e(2)=dble(zdotc(lv,v1,1,v2,1))
       call gate2_vec_o(v2,v1,lv,pzz,p1,p2)
       e(3)=dble(zdotc(lv,v1,1,v2,1))
       en = omp_get_wtime()
       if(en-st < mtime) mtime = en-st
       write(20,'(3es23.15)') e
    end do
99  write(20,*) "Min. of elapsed time in second is:", mtime
    close(10)
    close(20)
    deallocate(ths,list_p,v1,v2)
  end subroutine sub_main
  !
end program main


text:
c
c   Dint – version 2.0  is licensed under the Apache License, Version 2.0 (the "License");
c   you may not use Dint – version 2.0 except in compliance with the License.
c   You may obtain a copy of the License at
c       http://www.apache.org/licenses/LICENSE-2.0
c   The license is also given in the LICENSE file.
c   Unless required by applicable law or agreed to in writing, software
c   distributed under the License is distributed on an "AS IS" BASIS,
c   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
c   See the License for the specific language governing permissions and limitations under the License.
c
c -------------------------------------------------------------------------------------------
c  Dint : Direct Nonadiabatic Trajectories A code for non-Born–Oppenheimer molecular dynamics
c
c  version 2.0
c
c  A. W. Jasper
c  Argonne National Laboratory
c
c  Rui Ming Zhang
c  Tsinghua University
c
c  and
c
c  D. G. Truhlar
c  University of Minnesota
c
c  copyright  2020
c  Donald G. Truhlar and Regents of the University of Minnesota
c----------------------------------------------------------------------------------------------

      subroutine prediatom(vi,ji,rguess,mm,nsurf,symb,
     &    tau,rin,rout)

c Precompute some information for the atom-diatom initial conditions
c (INITx = 4).

      implicit none
      include 'param.f'

      integer nsurf,arri
      double precision mm(mnat),escat,rguess,
     & einti,eroti,evibi,eelec,rmin,rin,rout,ecoli,rmass,term,
     & mdum,ppreli,bmax,bmin,tau,etmp,
     & ji,vi
      character*2 symb(mnat)

      write(6,*)"Asymptotic analysis"
      write(6,*)"   Emin (eV)      Rmin (A)"
      call diamin2(rmin,etmp,-1.d0,mm,nsurf,symb,rguess)
      write(6,50)etmp*autoev,rmin*autoang
 50   format(2f15.5)

      write(6,*)"Preparing the diatom ",
     &  "with v = ",vi," j = ",ji," surface = ",nsurf
      write(6,*)
      eelec = etmp
      write(6,*)"Energetics (relative to zero of E)"
      write(6,*)"At the classical minimum:"
      write(6,100)"Re    = ",rmin*autoang," A"
      write(6,100)"Eelec = ",eelec*autoev," eV"
      call ewkb2(0.d0,vi,mm,evibi,rin,rout,nsurf,symb,rmin)
      call ewkb2(ji,vi,mm,einti,rin,rout,nsurf,symb,rmin)
      eroti = einti-evibi

      write(6,102)"Eint  = E(v=",vi,",j=",ji,") =",einti*autoev," eV"
      write(6,102)"Evib  = E(v=",vi,",j=",0.,") =",evibi*autoev," eV"
c     & ,2.d0*evibi*autocmi," cm-1"
      write(6,103)"Erot  = Eint  - Evib   = ",eroti*autoev," eV"
      write(6,101)"Turning points = ",rin*autoang," and ",
     & rout*autoang," A"
 100  format(1x,a,20x,f13.5,a)
 101  format(1x,a,f13.5,a,f13.5,a)
 102  format(1x,a,f6.2,a,f6.2,a,4x,f13.5,a,f13.5,a)
 103  format(1x,a,3x,f13.5,a)

      call period2(einti,ji,mm,rin,rout,tau,nsurf,symb)
      write(6,103)"Period for the vibration = ",tau*autofs," fs"
      write(6,*)

      return

      end

max_stars_repo_name: kevinmooreiii/DiNT

max_issues_repo_name: kevinmooreiii/DiNT

max_forks_repo_name: kevinmooreiii/DiNT


text:
c This module/program is part of DaliLite (c) L. Holm 1999
c
	program fssp
c
c   usage: ... fsspfilter.perl | pipeforwarder 1ppt dbsize dalidatpath
c		| fssp.$CPUARC | html.perl > 1ppt.html
c
c
	implicit none
	integer maxres
	parameter(maxres=6000)
c
c	step 1: read DCCP file -- assumed filtered, sorted
c	step 2: fetch DSSP info
c	step 3: statistics + topology
c	step 4: merge split blocks
c	step 5: write summary-line
c	step 6: write alignment-strings
c	step 7: write range-block
c	step 8: write u()t()-block
c
	character*5 cd1,cd2
	character*80 dalidatpath_1,dalidatpath_2,compnd1,compnd
	integer nprot2,nblock,l1(maxres),l2(maxres),r1(maxres),r2(maxres)
	real score,zscore
	integer i,j,k,ierr
c
	character seq(maxres),stru(maxres),aligned(maxres),seq1(maxres)
	character stru1(maxres),c
	integer nres,nchain,nres1,ali(maxres)
	real ca1(3,maxres),ca(3,maxres)
	character*5 resno1(0:maxres),resno(0:maxres)
	character*3 onetothree
	integer from1,from2,to1,to2
	character*6 strid1,strid2
	character*24 flags
	logical lpermut(maxres),lrevers(maxres),lincons(maxres)
	integer nrevers,npermut,lenali,ier,n,dbsize,ierr1
	real x(3,maxres),y(3,maxres),u(3,3),t(3),rmsd,ssq,w(maxres),ide
        real z(3,maxres)
	logical used(maxres)
	character*80 dssptemplate1,dssptemplate2,dssptemplate3

	character rotateaxis(3)
	real angle(3),x_deg,y_deg,z_deg
c
c	input parameters -- reads dccp file from stdin
c
	write(*,*) 'enter code, dbsize, dalidatpath_1, dalidatpath_2'
	read(*,500) cd1
	read(*,*) dbsize
        read(*,510) dalidatpath_1
	read(*,510) dalidatpath_2
c
c	step 0: write out dbsize
c
	write(*,470) dbsize
c
c	step 1: process query -- nchain, nres, cmpnd, seqlength
c
	call getdalidata(dalidatpath_1,cd1,99,nres1,seq1,stru1,compnd1,
     $		ca1,resno1)
	if(cd1(5:5).ne.' ') then
		strid1='code-X'
		strid1(6:6)=cd1(5:5)
	else
		strid1='code  '
	end if
	strid1(1:4)=cd1(1:4)
	write(*,480) strid1
	write(*,485) nres1,(resno1(i),i=1,nres1)
c
c	step 2: process DCCP file
c
	ierr=0
	nprot2=0
	do while (ierr.eq.0)
         call read_nextDCCP(cd1,cd2,score,zscore,nblock,l1,r1,l2,r2,
     $		ierr,5,maxres)
	 write(*,*) 'read_nextDCCP returned ',cd1,cd2,zscore,nblock,ierr
         if(ierr.eq.0) then
	  if(cd2(5:5).ne.' ') then
		strid2='code-X'
		strid2(6:6)=cd2(5:5)
	  else
		strid2='code  '
	  end if
	  strid2(1:4)=cd2(1:4)
          call getdalidata(dalidatpath_2,cd2,99,nres,seq,stru,compnd,
     $			ca,resno)
		write(*,*) '//'
		write(*,490) cd2
		write(*,521) '-compnd ','"',compnd,'"'
		do i=1,nres1
			ali(i)=0
		end do
		do i=1,nblock
			do j=l1(i),r1(i)
				k=l2(i)+j-l1(i)
				ali(j)=abs(k)
			end do
		end do
		! aligned sequence string, to compute pide
		do i=1,nres1
			k=ali(i)
			if(k.eq.0) then
				aligned(i)='.'
			else
				aligned(i)=seq(k)
			end if
		end do
        write(*,530) '-sequen ','"',(seq(i),i=1,nres),'"' ! unaligned sequence
		! sequence identity
		ide=0.0
		do i=1,nres1
		  if(ali(i).ne.0) then
			if(aligned(i).ge.'a'.and.seq1(i).ge.'a') then
				ide=ide+100.0
			else if(aligned(i).ge.'a'.and.seq1(i).eq.'C') then
				ide=ide+100.0
			else if(aligned(i).eq.'C'.and.seq1(i).ge.'a') then
				ide=ide+100.0
			else if(aligned(i).eq.seq1(i)) then
				ide=ide+100.0
			end if
		  end if
		end do
        write(*,530) '-struct ','"',(stru(i),i=1,nres),'"' ! unaligned dssp-sequence
		! merge blocks
		do i=2,nblock
			if(l1(i).eq.r1(i-1)+1.and.l2(i).eq.r2(i-1)+1) then
				write(*,*) 'merge blocks'
				l1(i)=l1(i-1)
				l2(i)=l2(i-1)
				l1(i-1)=0
				l2(i-1)=0
			end if
		end do
		k=nblock
		nblock=0
		do i=1,k
			if(l1(i).ne.0) then
				nblock=nblock+1
				l1(nblock)=l1(i)
				r1(nblock)=r1(i)
				l2(nblock)=l2(i)
				r2(nblock)=r2(i)
			end if
		end do
		call new_topology(l1,r1,l2,r2,lrevers,lpermut,lincons,
     $			nrevers,npermut,nblock,lenali,maxres,used)
		! range lines
		do i=1,nblock
			from1=l1(i)
			to1=r1(i)
			from2=l2(i)
			to2=r2(i)
			flags="                        "
			if(lrevers(i)) flags(1:8)=  'REVERSED'
			if(lpermut(i)) flags(10:17)='PERMUTED'
			if(lincons(i)) flags(19:24)='DOUBLE'
            write(*,550) strid1,strid2,
     $				from1,to1,from2,to2,
     $				onetothree(seq1(from1)),resno1(from1),
     $				onetothree(seq1(to1)),resno1(to1),
     $				onetothree(seq(from2)),resno(from2),
     $				onetothree(seq(to2)),resno(to2),flags
		end do
		! summary line:
		! z,rmsd,lali,lseq2,%ide,revers,permut,nfrag,topo,compnd
		n=0
		do i=1,nres1
		  k=abs(ali(i))
		  if(k.ne.0) then
			n=n+1
			do j=1,3
				x(j,n)=ca(j,k) 	! superimpose second on query
				y(j,n)=ca1(j,i)
			end do
			w(n)=1.0
			!write(*,*) 'pair: ',n,i,k,(x(j,n),j=1,3),(y(j,n),j=1,3)
		  end if
		end do
		do i=1,3
			t(i)=0.0
			do j=1,3
				u(i,j)=0.0
			end do
		end do
		if(n.gt.0) then
			call u3b(w,x,y,n,1,ssq,u,t,ier)
			rmsd=sqrt(ssq/n)
			!write(*,*) 'u3b gives ',n,ssq,rmsd,ier,u,t
		else
			rmsd=0.0
		end if

		do i=1,3
            write(*,570)  strid1,strid2,i,(u(i,j),j=1,3),t(i)
		end do

                ! convert u() to rotate x|y|z
                call u2xyzrotate(u,rotateaxis,x_deg,y_deg,z_deg)
                if(rotateaxis(1).eq.'x') then
                        angle(1)=x_deg
                        angle(2)=y_deg
                        angle(3)=z_deg
                else
                        angle(1)=y_deg
                        angle(2)=z_deg
                        angle(3)=x_deg
                end if
                write(*,580) (rotateaxis(i),angle(i),i=1,3),
     $                  (t(i),i=1,3)
580     format('-transrotate',3(1x,a1,f10.4),1x,3f10.3)

		if(nrevers.eq.0.and.npermut.eq.0) then
			c='S'
		else
			c='N'
		end if
        write(*,560) strid1,strid2,zscore,rmsd,lenali,nres,
     $		  nint(ide/max(1,n)),nrevers,npermut,nblock,c,
     $		  compnd
	 end if
       end do

470	format('-dbsize ',i5)
480	format('-query ',a6)
485	format('-nres ',i5,<maxres>(1x,a5))
486	format('-nchain ',i5)
490	format('-protein ',a5)
500	format(a5)
510	format(a80)
521	format(a8,1x,a1,a70,a1)
520	format(a8,1x,a1,a80,a1)
530	format(a8,1x,<maxres+2>a1)
540	format(a8,1x,a1,a33,a1)
550	format('-ranges  "',a6,1x,a6,2x,i4,' -',i4,' <=> ',i4,' -',i4,
     $		'   (',a3,1x,a5,' - ',a3,1x,a5,' <=> ',a3,1x,a5,' - ',
     $		a3,1x,a5,')',2x,a24,'"')
560	format('-summar  "',a6,1x,a6,f5.1,f5.1,i5,i6,i5,i7,i7,i6,1x,a1,
     $		4x,a70,'"')
570	format('-matrix  "',a6,1x,a6,'  U(',i1,',.) ',3f10.6,f20.6,'"')

	return
	end
c
c-------------------------------------------------------------------------------
c
        subroutine read_nextDCCP(cd1,cd2,score,zscore,
     $          nblock,l1,r1,l2,r2,ierr,iunit,maxres)
        implicit none
c
c       reads dali-scores & Z-scores from domainparser output !
c
        character*5 cd1,cd2,c1,c2
        character*9 line
        real score,zscore
	integer maxres,nblock,l1(maxres),l2(maxres),r1(maxres),r2(maxres)
        integer i,ierr,iunit
        integer ide,lali
        real rmsd
	character*132 longline
c
cDCCP   1   1651.8 0.0 175    24.6         100      1                1hgeB 1hgeBc
100	read(iunit,500,end=999) longline
	read(longline,601,err=100) line,score,rmsd,lali,zscore,
     $          	ide,nblock,c1,c2
        if(line(1:9).eq.'DCCP   1 ') then
                if(c1.eq.cd1) then      ! expected order
10                      read(iunit,605,err=100,end=999) line
                        if(line(1:9).ne.'alignment') goto 10
                        read(iunit,610,err=100) (l1(i),r1(i),i=1,nblock)
                        read(iunit,610,err=100) (l2(i),r2(i),i=1,nblock)
                        cd2=c2
                else if(c2.eq.cd1) then ! reverse pair !
20                      read(iunit,605,err=100,end=999) line
                        if(line(1:9).ne.'alignment') goto 20
                        read(iunit,610,err=100) (l2(i),r2(i),i=1,nblock)
                        read(iunit,610,err=100) (l1(i),r1(i),i=1,nblock)
                        cd2=c1
                else
                        goto 100
                end if
        else
                goto 100
        end if

500	format(a132)
600     format(1x,a9,f8.1,f4.1,i4,16x,i4,4x,i3,6x,2i4,2x,a5,1x,a5)
601     format(1x,a9,f8.1,f4.1,i4,f8.1,8x,i4,4x,i3,16x,a5,1x,a5)
605     format(1x,a9)
610     format(8(i4,2x,i4))

!       normal exit
        ierr=0
        return

!       error exit
999     ierr=-1
        return
	end
c
c-----------------------------------------------------------------------------
c
	subroutine new_topology(l1,r1,l2,r2,lrevers,lpermut,lincons,nrevers,
     $		npermut,nblock,lenali,maxres,used)
	implicit none
	integer maxres,l1(maxres),l2(maxres),r1(maxres),r2(maxres)
	integer nblock,lenali
	logical lpermut(maxres),lrevers(maxres),lincons(maxres),used(maxres)
	integer nrevers,npermut
c
	integer j,k,l,m

		lenali=0
		nrevers=0
		npermut=0
		do k=1,maxres
			used(k)=.false.
		end do
		do j=1,nblock
			lrevers(j)=.false.
			lpermut(j)=.false.
			lincons(j)=.false.
			if(l2(j).le.r2(j)) then
				l=l2(j)
				m=r2(j)
			else
				l=r2(j)
				m=l2(j)
			end if
			do k=l,m
				if(.not.used(k)) then
					used(k)=.true.
				else
					lincons(j)=.true.
				end if
			end do
			lenali=lenali+r1(j)-l1(j)+1
			if(l2(j).gt.r2(j)) then
				nrevers=nrevers+1
				lrevers(j)=.true.
			end if
			if(j.gt.1) then
				if(.not.lrevers(j-1).and.lrevers(j))
     $					lpermut(j)=.true.
				if((l2(j).lt.l2(j-1)).and.
     $	(.not.lrevers(j-1))) then
					npermut=npermut+1
					lpermut(j)=.true.
				end if
				if((l2(j).lt.l2(j-1)).and.
     $	(lrevers(j-1)).and.(.not.lrevers(j))) then
					npermut=npermut+1
					lpermut(j)=.true.
				end if
				if((l2(j).gt.l2(j-1)).and.
     $	(lrevers(j-1))) then
					npermut=npermut+1
					lpermut(j)=.true.
				end if
			end if
		end do

	return
	end

max_stars_repo_name: gkrash/lcspromatch

max_issues_repo_name: gkrash/lcspromatch

max_forks_repo_name: gkrash/lcspromatch


text:
C
C $Id: mapbd.f,v 1.31 2008-09-18 12:31:28 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MAPBD
C
C Calling this do-nothing subroutine forces "ld" to load the following
C block data routine (but only if they are in the same ".f" file).
C
        RETURN
C
      END
CNOSPLIT
      BLOCKDATA MAPBDX
C
        PARAMETER (MNAI=8000)
C
C The common block MAPCM0 contains mathematical constants.
C
        COMMON /MAPCM0/  COS1,DTOR,DTRH,OOPI,PI,PIOF,PIOT,RTDD,RTOD,
     +                   SROT,SIN1,TOPI,TSRT
        DOUBLE PRECISION COS1,DTOR,DTRH,OOPI,PI,PIOF,PIOT,RTDD,RTOD,
     +                   SROT,SIN1,TOPI,TSRT
        SAVE   /MAPCM0/
C
C The common block MAPCM1 contains transformation constants.
C
        COMMON /MAPCM1/  COSO,COSR,SINO,SINR,IPRJ,IROD
        DOUBLE PRECISION COSO,COSR,SINO,SINR
        INTEGER          IPRJ,IROD
        SAVE   /MAPCM1/
C
C The common block MAPCM2 contains area-specification variables.
C
        COMMON /MAPCM2/  BLAM,BLOM,PEPS,SLAM,SLOM,UCEN,UMAX,UMIN,UOFF,
     +                   URNG,VCEN,VMAX,VMIN,VOFF,VRNG,ISSL
        DOUBLE PRECISION BLAM,BLOM,PEPS,SLAM,SLOM,UCEN,UMAX,UMIN,UOFF,
     +                   URNG,VCEN,VMAX,VMIN,VOFF,VRNG
        INTEGER          ISSL
        SAVE   /MAPCM2/
C
C The common block MAPCM3 contains parameters having to do with reading
C the data for outlines.
C
        COMMON /MAPCM3/  ITPN,NOUT,NPTS,IGID,IDLS,IDRS,BLAG,SLAG,BLOG,
     +                   SLOG,PNTS(200),IDOS(4)
        INTEGER          ITPN,NOUT,NPTS,IGID,IDLS,IDRS,IDOS
        REAL             BLAG,SLAG,BLOG,SLOG,PNTS
        SAVE   /MAPCM3/
C
C The common block MAPCM4 contains most of the input parameters.
C
        COMMON /MAPCM4/  GRDR,GRID,GRLA,GRLO,GRPO,OTOL,PDRE,PLA1,PLA2,
     +                   PLA3,PLA4,PLB1,PLB2,PLB3,PLB4,PLNO,PLTO,ROTA,
     +                   SRCH,XLOW,XROW,YBOW,YTOW,IDOT,IDSH,IDTL,ILCW,
     +                   ILTS,JPRJ,ELPF,INTF,LBLF,PRMF
        DOUBLE PRECISION GRDR,GRID,GRLA,GRLO,GRPO,OTOL,PDRE,PLA1,PLA2,
     +                   PLA3,PLA4,PLB1,PLB2,PLB3,PLB4,PLNO,PLTO,ROTA,
     +                   SRCH,XLOW,XROW,YBOW,YTOW
        INTEGER          IDOT,IDSH,IDTL,ILCW,ILTS,JPRJ
        LOGICAL          ELPF,INTF,LBLF,PRMF
        SAVE   /MAPCM4/
C
C The common block MAPCM5 contains various lists ("dictionaries") of
C two-character codes required by EZMAP for parameter-setting.
C
        COMMON /MAPCM5/  DDCT(5),DDCL(5),LDCT(6),LDCL(6),PDCT(19),
     +                   PDCL(19)
        CHARACTER*2      DDCT,DDCL,LDCT,LDCL,PDCT,PDCL
        SAVE   /MAPCM5/
C
C The common block MAPCM6 contains quantities needed by the routines
C MDPTRA and MDPTRI.
C
        COMMON /MAPCM6/  UCNM,UMNM,UMXM,URNM,VCNM,VMNM,VMXM,VRNM,ELPM
        DOUBLE PRECISION UCNM,UMNM,UMXM,URNM,VCNM,VMNM,VMXM,VRNM
        LOGICAL ELPM
        SAVE   /MAPCM6/
C
C The common block MAPCM7 contains parameters describing the portion of
C the plotter frame being used.
C
        COMMON /MAPCM7/  ULOW,UROW,VBOW,VTOW
        DOUBLE PRECISION ULOW,UROW,VBOW,VTOW
        SAVE   /MAPCM7/
C
C The common block MAPCM8 contains parameters set by MDPTRN and used by
C MAPIT in handling "cross-over" problems.
C
        COMMON /MAPCM8/  P,Q,R
        DOUBLE PRECISION P,Q,R
        SAVE   /MAPCM8/
C
C The common block MAPCMA contains values which are used to position
C dots along dotted outlines and to avoid drawing vectors which are
C too short.
C
        COMMON /MAPCMA/  DATL,DBTD,DDTS,DPLT,DPSQ,DSCA,DSSQ
        DOUBLE PRECISION DATL,DBTD,DDTS,DPLT,DPSQ,DSCA,DSSQ
        SAVE   /MAPCMA/
C
C The common block MAPCMC contains variables for MDPBLA.
C
        COMMON /MAPCMC/  IGI1,IGI2,NCRA,NOVS,XCRA(100),YCRA(100)
        INTEGER          IGI1,IGI2,NCRA,NOVS
        REAL             XCRA,YCRA
        SAVE   /MAPCMC/
C
C The common block MAPCMP contains the buffers in which the X and Y
C coordinates of points are collected for an eventual call to POINTS.
C
        COMMON /MAPCMP/  NPTB,XPTB(50),YPTB(50)
        INTEGER          NPTB
        REAL             XPTB,YPTB
        SAVE   /MAPCMP/
C
C The common block MAPCMQ contains quantities specifying the color
C indices to be used for various portions of the plot.
C
        COMMON /MAPCMQ/  ICIN(8)
        INTEGER          ICIN
        SAVE   /MAPCMQ/
C
C The common block MAPCMW contains variables affecting the definitions
C of the equirectangular, cylindrical equal-area, and Winkel tripel
C projections.
C
        COMMON /MAPCMW/  CSLS,CSLT,SLTD,ISLT
        DOUBLE PRECISION CSLS,CSLT,SLTD
        INTEGER ISLT
        SAVE  /MAPCMW/
C
C The common blocks MAPCMX and MAPCMY contain variables used by the new
C routines in EZMAPB (implementing access to the new dataset created in
C early 1998).  MAPCMX contains integer variables and MAPCMY contains
C character variables.
C
        COMMON /MAPCMX/  IATY(MNAI),ISCI(MNAI),IPAR(MNAI)
        INTEGER          IATY,ISCI,IPAR
        SAVE   /MAPCMX/
C
        COMMON /MAPCMY/  NAME(MNAI),FLNS
        CHARACTER*64     NAME
        CHARACTER*512    FLNS
        SAVE   /MAPCMY/
C
C The common block MAPRGD contains two arrays of color indices and
C a couple of point-interpolation values used by the routines that plot
C data from the RANGS/GSHHS dataset.
C
        COMMON /MAPRGD/  ICOL(5),ICSF(5),IDPF,LCRA,NILN,NILT,OLAT,OLON
        INTEGER          ICOL,ICSF,IDPF,LCRA,NILN,NILT
        REAL             OLAT,OLON
        SAVE   /MAPRGD/
C
C The common block MAPSAT contains parameters for the satellite-view
C projection.
C
        COMMON /MAPSAT/  ALFA,BETA,DCSA,DCSB,DSNA,DSNB,SALT,SSMO,SRSS
        DOUBLE PRECISION ALFA,BETA,DCSA,DCSB,DSNA,DSNB,SALT,SSMO,SRSS
        SAVE   /MAPSAT/
C
C Below are descriptions of the variables in each of the common blocks,
C together with data statements giving default values to those variables
C which need default variables.
C
C Variables in MAPCM0:
C
C COS1 is the cosine of 1 degree; DTOR is a conversion constant to go
C from degrees to radians; DTRH is a conversion constant to go from
C degrees to radians, halved; OOPI is one over pi; PI is pi; PIOT is
C pi over two; PIOF is pi over four; RTDD is a conversion constant to go
C from radians to degrees, doubled; RTOD is a conversion constant to go
C from radians to degrees; SIN1 is the sine of 1 degree; SROT is the
C square root of two; TOPI is two over pi; and TSRT is two times the
C square root of two.
C
        DATA COS1 / .999847695156390D0 /
        DATA DTOR / .017453292519943D0 /
        DATA DTRH / .008726646259971D0 /
        DATA OOPI / .318309886183790D0 /
        DATA PI   / 3.14159265358979D0 /
        DATA PIOT / 1.57079632679489D0 /
        DATA PIOF / .785398163397445D0 /
        DATA RTDD / 114.591559026165D0 /
        DATA RTOD / 57.2957795130823D0 /
        DATA SIN1 / .017452406437283D0 /
        DATA SROT / 1.41421356237310D0 /
        DATA TOPI / .636619772367581D0 /
        DATA TSRT / 2.82842712474619D0 /
C
C Variables in MAPCM1:
C
C IPRJ is an integer between 0 and 25, specifying what projection is
C currently in use.  The following table defines the meanings of the
C various possible values of IPRJ.  In the table, the word "arbitrary"
C is used to indicate the the projection may be centered at any point
C on the globe and the word "fast-path" is used to indicate that the
C projection is centered at (0,0) on the globe (which usually speeds
C up the evaluation of the equations).  Projections 16 through 24 are
C fast-path versions of projections 7 through 15.
C
C     IPRJ    Projection Selected                         Type
C     ----    ------------------------------------        -----------
C       0     (UT) USGS transformations                   Various
C       1     (LC) Lambert Conformal                      Conical
C       2     (ST) Stereographic                          Azimuthal
C       3     (OR or SA) Orthographic or Satellite        Azimuthal
C       4     (LE) Lambert Equal-Area                     Azimuthal
C       5     (GN) Gnomonic                               Azimuthal
C       6     (AE) Azimuthal Equidistant                  Azimuthal
C       7     (CE) Cylindrical Equidistant (arbitrary)    Cylindrical
C       8     (ME) Mercator (arbitrary)                   Cylindrical
C       9     (MT) Mollweide-type (arbitrary)             Cylindrical
C      10     (RO) Robinson (arbitrary)                   Cylindrical
C      11     (EA) Cylindrical Equal-Area (arbitrary)     Cylindrical
C      12     (AI) Aitoff (arbitrary)                     Mixed
C      13     (HA) Hammer (arbitrary)                     Mixed
C      14     (TM) True Mollweide (arbitrary)             Mixed
C      15     (WT) Winkel tripel (arbitrary)              Mixed
C      16     (CE) Cylindrical Equidistant (fast-path)    Cylindrical
C      17     (ME) Mercator (fast-path)                   Cylindrical
C      18     (MT) Mollweide-type (fast-path)             Cylindrical
C      19     (RO) Robinson (fast-path)                   Cylindrical
C      20     (EA) Cylindrical Equal-Area (fast-path)     Cylindrical
C      21     (AI) Aitoff (fast-path)                     Mixed
C      22     (HA) Hammer (fast-path)                     Mixed
C      23     (MO) True Mollweide (fast-path)             Mixed
C      24     (WT) Winkel tripel (fast-path)              Mixed
C      25     (RM) Rotated Mercator (fast-path)           Cylindrical
C
C IROD is a flag which, if non-zero, says that we have to use double
C precision at selected points.  SINO, COSO, SINR, and COSR are
C projection variables computed by MDPINT for use by MDPTRN.
C
      DATA IROD / 0 /
C
C Variables in MAPCM2:
C
C UMIN, UMAX, VMIN, and VMAX specify the limits of the rectangle to be
C drawn, in projection space.  UCEN, VCEN, URNG, and VRNG are computed
C by MDPINT for use when the map perimeter is made elliptical (by
C setting the flag ELPF).  BLAM, SLAM, BLOM, and SLOM are respectively
C the biggest latitude, the smallest latitude, the biggest longitude,
C and the smallest longitude on the map.  They are used in MDPGRD and
C in MDPLOT to make the drawing of grids and outlines more efficient.
C ISSL is a flag set by MDPINT to indicate the success (ISSL=1) or
C failure (ISSL=0) of the code setting BLAM, SLAM, BLOM, and SLOM;
C failure indicates that default values were used and that it is not
C safe to take the efficient paths in MDPGRD.  PEPS is set by MDPINT
C for use in MAPIT in testing for cross-over problems.  UMIN and UMAX
C are given default values to prevent code in MDSETI and MDSETR from
C blowing up when PDRE is set prior to the first call to MDPINT.
C (08/07/2001) UOFF and VOFF are being added to this common block to
C support the use of U/V coordinates from which an offset has been
C subtracted to avoid loss of resolution when zooming in on a small
C portion of a projection.
C
      DATA UMIN,UMAX / 0.D0,1.D0 /
C
C Variables in MAPCM3:
C
C ITPN is the unit number of the "tape" from which outline data is to
C be read.  NOUT is the number of the outline to be used; the values 0
C through 5 imply 'NO', 'CO', 'US', 'PS', and 'PO', respectively; thus,
C if NOUT is zero, no outlines are to be used, and, if it is non-zero,
C it is the number of the "file" to be read from unit ITPN.  NPTS, just
C after a read, is the number of elements read into PNTS; it is then
C divided by 2 to become the number of points defined by the group just
C read.  IGID is an identifier for the group, so that, for example, one
C can distinguish a group belonging to a international boundary from
C one belonging to a U.S. state boundary.  IDLS and IDRS are integers
C identifying the areas to the left and right of the edge defined by
C the group of points just read.  BLAG, SLAG, BLOG, and SLOG specify
C the biggest and smallest latitude and the biggest and smallest
C longitude of the points in the group, so that, in some cases at least,
C one can decide quickly not to bother with the group.  PNTS contains
C NPTS coordinate pairs, each consisting of a latitude and a longitude,
C in degrees.  IDOS is an array of area-identifier offsets; for a given
C EZMAP dataset with index NOUT, the area identifiers in the dataset
C are biased by the addition of IDOS(NOUT) in order to generate a
C unique value.
C
      DATA ITPN,NOUT / 1,1 /
C
      DATA IDOS / 0,222,439,1003 /
C
C Variables in MAPCM4:
C
C INTF is a flag whose value at any given time indicates whether the
C package EZMAP is in need of initialization (.TRUE.) or not (.FALSE.).
C JPRJ is an integer between 0 and 13 indicating the type of projection
C currently in use.  PLTO, PLNO, and ROTA are the pole latitude and
C longitude and the rotation angle specified by the last user call to
C MAPROJ.  ILTS is an integer between 1 and 6, specifying how the limits
C of the map are to be chosen.  PLA1-4 and PLB1-4 are the values given
C by the user for PLM1(1), PLM2(1), ..., PLM1(2), PLM2(2), ..., in the
C last call to MAPSET.  PDRE is the plotter resolution - effectively,
C the number of addressable points in the x direction.  GRID, if zero,
C turns off the drawing of the lat/lon grid; if non-zero, it specifies
C the desired default spacing between grid lines, in degrees.  IDSH is
C the desired dash pattern (16-bit binary) for grid lines.  IDOT is a
C flag selecting solid outlines (0) or dotted outlines (1).  LBLF is
C a logical flag indicating whether the international date line, the
C equator, the Greenwich meridian, and the poles are to be labelled or
C not.  PRMF is a logical flag indicating whether or not a perimeter
C is to be drawn.  ELPF is a logical flag indicating whether the map
C perimeter is to be rectangular (.FALSE.) or elliptical (.TRUE.).
C XLOW, XROW, YBOW, and YTOW are fractions between 0. and 1. specifying
C the position of area of the plotter frame in which the map is to be
C put; the map is centered in this area and made as large as possible.
C IDTL is a flag specifying that MAPIT should draw solid outlines (0)
C or dotted outlines (1).  GRDR and SRCH are measured in degrees and
C lie in the range from .001 to 10.  GRDR specifies the resolution with
C which the grid is to be drawn and SRCH the accuracy with which the
C latitude/longitude limits of the map are to be found.  ILCW is the
C character width for characters in the label, as required for use in a
C call to PWRIT.  GRLA and GRLO specify the desired spacings of latitude
C and longitude grid lines, respectively (in degrees); if either is less
C than or equal to zero, the value of GRID is used instead.  GRPO
C specifies the way in which the lat/lon grid is to be modified near
C the poles on projections which map the poles into single points; it
C has a value of the form 1000*GRP1+GRP2, where GRP1 is an integer
C between 0 and 90 and GRP2 is a positive real between 0 and 360.  If
C GRP1 is zero, all the latitude lines of the grid are drawn; if GRP1
C is non-zero, latitude lines at latitudes from GRP1 to 90, inclusive
C (South or North) are omitted.   If GRP2 is zero, no longitude lines
C are drawn near the poles; if GRP2 is non-zero, only longitude lines
C of the grid at multiples of GRP2 are drawn near the poles.  Examples:
C GRPO=0 says "draw all the latitude lines of the grid; omit longitude
C lines of the grid near the poles."  GRPO=1 says "draw entire grid."
C GRPO=75045 says "suppress latitude lines of the grid above 75N and
C below 75S; near the poles, draw longitude lines of the grid only at
C multiples of 45 degrees."  The default is GRPO=90, which says "draw
C all latitude lines of the grid; near the poles, omit longitude lines
C of the grid except those at multiples of 90 degrees".  OTOL is an
C "offset tolerance" (new in 08/08/2001): when the absolute value of the
C difference of the U values at the left and right edges of a map is
C less than OTOL times the absolute value of their sum, or the absolute
C value of the difference of the V values at the bottom and top edges of
C a map is less than OTOL times the absolute value of their sum, MDPINT
C will initialize the projection to use offset values of U and V, thus
C avoiding a loss of precision.
C
      DATA INTF,JPRJ,PLTO,PLNO,ROTA,ILTS,PLA1,PLA2,PLA3,PLA4,PLB1,PLB2 /
     1   .TRUE.,   7,0.D0,0.D0,0.D0,   1,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0 /
C
      DATA PLB3,PLB4,    PDRE, GRID, IDSH,IDOT, LBLF , PRMF ,  ELPF  /
     1     0.D0,0.D0,32768.D0,10.D0,21845,   0,.TRUE.,.TRUE.,.FALSE. /
C
      DATA IDTL,  OTOL /
     1        0,.001D0 /
C
      DATA XLOW,XROW,YBOW,YTOW / .05D0,.95D0,.05D0,.95D0 /
C
      DATA GRDR,SRCH / 1.D0,1.D0 /
C
      DATA ILCW / 1 /
C
      DATA GRLA,GRLO,GRPO / 0.D0,0.D0,90.D0 /
C
C Variables in MAPCM5:
C
C DDCT is the dictionary of available datasets, LDCT the dictionary of
C map limit definition types, and PDCT the dictionary of map projection
C names.  DDCL, LDCL, and PDCL are lower-case equivalents.
C
      DATA DDCT
     +        / 'NO','CO','US','PS','PO' /
      DATA DDCL
     +        / 'no','co','us','ps','po' /
C
      DATA LDCT
     +        / 'MA','CO','PO','AN','LI','GR' /
      DATA LDCL
     +        / 'ma','co','po','an','li','gr' /
C
      DATA PDCT
     +        / 'UT','LC','ST','OR','LE','GN','AE','CE','ME','MT',
     +          'RO','EA','AI','HA','MO','WT','SV','ER','RM' /
      DATA PDCL
     +        / 'ut','lc','st','or','le','gn','ae','ce','me','mt',
     +          'ro','ea','ai','ha','mo','wt','sv','er','rm' /
C
C Variables in MAPCM6:
C
C All of the variables in this common block are basically copies of
C variables in other common blocks (perhaps somewhat modified) for use
C by transformation routines.  ELPM is a copy of the logical variable
C ELPF, which indicates whether the map perimeter is to be rectangular
C or elliptical.  UMNM, UMXM, UMXM, and VMNM specify a rectangular
C perimeter just a little outside the one specified by UMIN, UMAX, VMIN,
C and VMAX.  UCNM, VCNM, URNM, and VRNM are just copies of the values
C UCEN, VCEN, URNG, and VRNG, respectively; they are used when the map
C perimeter is elliptical.  None of these parameters needs a default
C value.
C
C Variables in MAPCM7:
C
C ULOW, UROW, VBOW, and VTOW define the fraction of the plotter frame
C to be occupied by the map - they may be thought of as the first four
C arguments of the SET call or, in the GKS scheme, as the viewport.
C They are computed by MDPINT.  ULOW and UROW are given default values
C to prevent code in MDSETI and MDSETR from blowing up when PDRE is
C set prior to the first call to MDPINT.
C
      DATA ULOW,UROW / 0.D0,1.D0 /
C
C Variables in MAPCM8:
C
C P, Q, and R are set by MDPTRN each time it maps (RLAT,RLON) to (U,V).
C Q is always equal to V, but P is not always equal to U.  Instead, it
C is a value of U from an intermediate step in the projection process.
C For the Lambert conformal conic, P is the distance, in longitude, from
C the central meridian.  For the cylindrical projections, P is a value
C of U prior to multiplication by a function of V shrinking the map
C toward a vertical bisector.  They are all used by MAPIT, while drawing
C lines from point to point, to detect "cross-over" (a jump from one
C side of the map to the other, caused by the projection's having slit
C the globe along some half of a great circle and laid it open with the
C two sides of the slit at opposite ends of the map).
C
C Variables in MAPCMA:
C
C DPLT is the mimimum vector length; MAPIT requires two points to be at
C least DPLT plotter units apart before it will join them with a vector.
C DDTS is the desired distance in plotter units between dots in a dotted
C outline.  These values are relative to the "plotter resolution" PDRE;
C DPLT/PDRE is a fraction of the plotter frame.  DSCA is the ratio of
C the length of a vector, measured in plotter units, to the length of
C the same vector, measured in the u/v plane.  Thus, given a vector of
C length D in the u/v plane, D*DSCA is its length in plotter units.
C DPSQ and DSSQ are the squares of DPLT and DSCA, respectively.  DBTD
C is the distance, in the u/v plane, between two dots DDTS plotter
C units apart.  DPLT and DDTS have the values given below and are not
C reset by the code; DSCA, DPSQ, DSSQ, and DBTD are computed by MDPINT.
C DSCA is given a default value only to keep the routines MDSETI and
C MDSETR from blowing up when DDTS is set prior to any call to MDPINT.
C DATL is used by MAPIT and MDPVP to keep track of where the next point
C along a curve should go.
C
      DATA DPLT,DDTS,DSCA / 1.D0,96.D0,1.D0 /
C
C Variables in MAPCMC:
C
C IGI1 is the group identifier for the principal edge group - the one
C containing all the real map boundaries.  IGI2 is the group identifier
C for the dummy group that splits the map into vertical strips.  NOVS
C is the number of vertical strips to be used.  XCRA and YCRA hold X
C and Y coordinates for delivery to AREDAM.  NCRA is the number of
C coordinates in XCRA and YCRA.
C
      DATA IGI1,IGI2,NOVS,NCRA / 1,2,1,0 /
C
C Variables in MAPCMP:
C
C NPTB is the number of points whose coordinates have been collected in
C the arrays XPTB and YPTB for eventual output by a call to POINTS.
C
      DATA NPTB / 0 /
C
C Variables in MAPCMQ:
C
C The array ICIN specifies color indices to be used for the perimeter,
C for the grid, for labelling, for limbs, for the continental outlines,
C for the U.S. state outlines, and for international political outlines.
C See the routine MDPCHI.  The default values are all negative, to say
C that color should not be set.
C
      DATA ICIN / -1,-1,-1,-1,-1,-1,-1,-1 /
C
C Variables in MAPCMW:
C
C SLTD is a "standard latitude" parameter affecting the behavior of
C three projections: the equirectangular, the cylindrical equal-area,
C and the Winkel tripel.  CSLT is the cosine of SLTD and CSLS is the
C square of CSLT.  When IPRJ is 7, ISLT is 0 if the cylindrical
C equidistant projection is to be used or 1 if the equirectangular
C projection is to be used.
C
      DATA CSLS,CSLT,SLTD,ISLT / 1.D0,1.D0,-1.D0,0 /
C
C Variables in MAPCMX and MAPCMY:
C
C For each I from 1 to MNAI, IATY(I), if non-zero, is the type of the
C area with area identifier I, ISCI(I) is a suggested color index for
C that area, IPAR(I) is the index of the parent area, and NAME(I) is
C the name of the area.  FLNS is the full name of the file from which
C the data in IATY, ISCI, IPAR, and NAME were read; it allows us to
C avoid re-reading the data when that isn't necessary.
C
      DATA (IATY(I),I=1,MNAI) / MNAI*0   /
      DATA (ISCI(I),I=1,MNAI) / MNAI*0   /
      DATA (IPAR(I),I=1,MNAI) / MNAI*0   /
      DATA (NAME(I),I=1,MNAI) / MNAI*' ' /
C
      DATA FLNS / ' ' /
C
C Variables in MAPRGD:
C
C The array ICOL contains five color indices to be used in drawing
C the outlines of polygons of the five different types defined by the
C RANGS/GSHHS data, and the array ICSF contains five color indices to
C be used in filling the polygons.  The routine MDRGSC may be called
C by the user to set these colors and the routine MDRGGC may be called
C to retrieve them.  The indices of each of these arrays are assigned
C to features as follows: 1 => ocean; 2 => land; 3 => lake; 4 =>
C island in lake; and 5 => pond on island.  NILN and NILT determine
C how many points are to be interpolated per degree of longitude or
C degree of latitude, respectively, along those parts of the filled
C polygons lying along the edges of 1-degree cells.  The object of this
C interpolation is to reduce the probability of having little gaps
C along the seams between adjacent cells.  (If it is not done and a
C projection is used that maps the 1-degree cells into patches whose
C edges are neither horizontal nor vertical, the gaps can be quite
C noticeable.)  OLAT and OLON are temporary variables in which some
C values are passed from MPRGSX to MPRGFA; they do not need to be
C initialized.  IDPF is a flag indicating how the RANGS/GSHHS data
C to be processed.  IDPF = 0 means that the data are to be used in
C are unedited form (MDRGSQ is called).  IDPF = 1 means that the data
C are to be edited and that polygon fills are to be done using direct
C calls to GFA, without using an area map (MDRGSX is called).  IDPF =
C 2 means that the data are to be edited and that fills are to be done
C using an area map (MDRGSX is called).
C
      DATA ICOL / 1 , 1 , 1 , 1 , 1 /
      DATA ICSF / 0 , 1 , 0 , 1 , 0 /
C
      DATA IDPF / 1 /
C
      DATA NILN,NILT / 50 , 50 /
C
C Variables in MAPSAT:
C
C The absolute value of SALT, if greater than 1, serves as a flag that
C a satellite-view projection is to be used in place of an orthographic
C projection; its value is the distance of the satellite from the center
C of the earth, in units of earth radii.  In this case, SSMO is the
C square of SALT minus 1 and SRSS is the square root of SSMO.  If ALFA
C is zero, the projection shows the view seen by a satellite looking
C straight at the center of the earth; call this the basic satellite
C view.  If ALFA is non-zero, it and BETA are angles, in degrees,
C determining where the line of sight of the projection is.  If E is
C at the center of the earth, S is at the satellite, and P is a point
C along the line of sight, then ALFA measures the angle ESP.  If O is
C the point at the origin of the basic satellite view and P is the
C projection of the line of sight, then BETA measures the angular
C distance from the positive u axis to the line OP, positive if
C measured counter-clockwise.  DSNA, DCSA, DSNB, and DCSB are sines
C and cosines of ALFA and BETA.  The sign of SALT says whether a normal
C projection (positive) or an extended projection (negative) is to be
C used.  The latter makes it easier to overlay CONREC output on one of
C these projections, by projecting points out of sight around the limb
C to points just outside the limb on the projected view.
C
      DATA SALT,ALFA,BETA / 0.D0,0.D0,0.D0 /
C
      DATA DSNA,DCSA,DSNB,DCSB / 0.D0,1.D0,0.D0,1.D0 /
C
      END
CNOSPLIT
      BLOCKDATA MAPBDQ
C
C The common block MAQCMN contains only those variables that are needed
C for MDQTRA, MDQTRI, and MDQTRN to carry out the transformation that
C was in effect at the time MDQINI was called.  All of them are copies
C of variables in other common blocks and are described above, in the
C commenting for BLOCKDATA MAPBDX.  They are given default values here
C mostly to protect the transformation routines from blowing up if they
C are erroneously called prior to a call to MDQINI.
C
        COMMON /MAQCMN/  ALFA,COSO,COSR,CSLS,CSLT,DCSA,DCSB,DSNA,DSNB,
     +                   DTOR,DTRH,OOPI,PLNO,  PI,PIOT,ROTA,RTDD,RTOD,
     +                   SALT,SINO,SINR,SRSS,SSMO,TOPI,UCNM,UMNM,UMXM,
     +                   UOFF,URNM,VCNM,VMNM,VMXM,VOFF,VRNM,UTPA,IPRF,
     +                   IPRJ,IROD,ELPM
        DOUBLE PRECISION ALFA,COSO,COSR,CSLS,CSLT,DCSA,DCSB,DSNA,DSNB,
     +                   DTOR,DTRH,OOPI,PLNO,  PI,PIOT,ROTA,RTDD,RTOD,
     +                   SALT,SINO,SINR,SRSS,SSMO,TOPI,UCNM,UMNM,UMXM,
     +                   UOFF,URNM,VCNM,VMNM,VMXM,VOFF,VRNM,UTPA(15)
C
        INTEGER IPRF,IPRJ,IROD
C
        LOGICAL ELPM
C
        SAVE   /MAQCMN/
C
        DATA ALFA / 0.D0 /
        DATA COSO / 1.D0 /
        DATA COSR / 1.D0 /
        DATA CSLS / 1.D0 /
        DATA CSLT / 1.D0 /
        DATA DCSA / 1.D0 /
        DATA DCSB / 1.D0 /
        DATA DSNA / 0.D0 /
        DATA DSNB / 0.D0 /
        DATA DTOR / .017453292519943D0 /
        DATA DTRH / .008726646259971D0 /
        DATA ELPM / .FALSE. /
        DATA IPRJ / 16 /
        DATA IROD / 1 /
        DATA OOPI / .318309886183790D0 /
        DATA PLNO / 0.D0 /
        DATA PI   / 3.14159265358979D0 /
        DATA PIOT / 1.57079632679489D0 /
        DATA ROTA / 0.D0 /
        DATA RTDD / 114.591559026165D0 /
        DATA RTOD / 57.2957795130823D0 /
        DATA SALT / 0.D0 /
        DATA SINO / 0.D0 /
        DATA SINR / 0.D0 /
        DATA SRSS / 0.D0 /
        DATA SSMO / 0.D0 /
        DATA TOPI / .636619772367581D0 /
        DATA UCNM / 0.D0 /
        DATA UMNM / -180.000360D0 /
        DATA UMXM / +180.000360D0 /
        DATA UOFF / 0.D0 /
        DATA URNM / 0.D0 /
        DATA VCNM / 0.D0 /
        DATA VMNM / -90.000180D0 /
        DATA VMXM / +90.000180D0 /
        DATA VOFF / 0.D0 /
        DATA VRNM / 0.D0 /
C
        DATA IPRF / 0 /
C
        DATA UTPA / 15*0.D0 /
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


program intel_dgesvd

use,intrinsic:: iso_fortran_env, only: sp=>real32, dp=>real64, compiler_version, stderr=>error_unit

implicit none (type, external)

external :: dgesvd

! Intel(R) Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64,
! Version 19.0.1.144 Build 20181018
! 64 bits: error mag:   0.000000000000E+00  0.000000000000E+00  0.000000000000E+00

! GCC version 7.3.0  (system Lapack)
! 64 bits: error mag:   0.000000000000E+00  2.220446049250E-16  4.996003610813E-16

! GCC version 7.3.0  (MKL 19)
! 64 bits: error mag:   6.952923663893-310  6.941716482845-310  4.940656458412-324

! pgf90 18.10-0 (system Lapack
! 64 bits: error mag:   0.000000000000E+00  2.220446049250E-16  4.996003610813E-16

! pgf90 18.10-0 (MKL 19)
! 64 bits: error mag:   0.000000000000E+00  0.000000000000E+00  0.000000000000E+00

real(dp):: A(3,3) = reshape( &
        [ 1, 0, 1, &
         -1,-2, 0, &
          0, 1,-1 ], &
           shape(A))

integer :: svdinfo
integer, parameter :: M = size(a,1), N = size(a,2), Lratio = 5
integer, parameter :: LWORK = LRATIO*M
real(dp) :: U(M,M),VT(M,M), SWORK(LRATIO*M)

real(dp) :: ss64(M, N), s64(M), e64(M), maxerr
real(dp), parameter :: s64ref(3) = [2.460504870018764_dp, 1.699628148275318_dp, 0.239123278256554_dp]

call dgesvd('A','N',M,M,A,M,SS64,U,M,VT,M, SWORK, LWORK,svdinfo)
if (svdinfo /=0) then
  write(stderr,*) 'svd error: ',svdinfo
  error stop
endif

s64 = ss64(:,1)  ! NOT DIAG!!!!

!e64 = s64 - S64ref

print *,compiler_version()
print '(I3,A,3ES20.12)',storage_size(s64),' bits: error mag: ',e64

maxerr=maxval(abs(e64))

if (maxerr > 1e-7_dp) error stop 'excessive singular value error'

end program