text:
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       this is the end of the debugging code, and the beginning
c       of the python version of quaplot routines - and then some...
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
        subroutine zpyplot(iw,zs,n,itype,title)
        implicit real *8 (a-h,o-z)
        real *8 zs(2,1)
        real *8, allocatable :: xs(:),ys(:)
        character *1 title(1)
c
        allocate(xs(n), ys(n))
c
        do 1600 i=1,n
        xs(i)=zs(1,i)
        ys(i)=zs(2,i)
 1600 continue
c
        call pyplot(iw,xs,ys,n,itype,title)
c
        return
        end
c
c
c
c
c
        subroutine zpyplot2(iw, zs, n, itype, zs2, n2,
     1      itype2, title)
        implicit real *8 (a-h,o-z)
        real *8 zs(2,1),zs2(2,1)
        real *8, allocatable :: xs(:),ys(:),xs2(:),ys2(:)
        character *1 title(1)
c
        allocate(xs(n), ys(n), xs2(n2), ys2(n2))
c
        do 1600 i=1,n
        xs(i)=zs(1,i)
        ys(i)=zs(2,i)
 1600 continue
c
        do 1800 i=1,n2
        xs2(i)=zs2(1,i)
        ys2(i)=zs2(2,i)
 1800 continue
c
        call pyplot2(iw,xs,ys,n,itype,xs2,ys2,n2,itype2,title)
c
        return
        end
c
c
c
c
c
        subroutine zpyplot3(iw, zs, n, itype, zs2, n2, itype2,
     1      zs3, n3, itype3, title)
        implicit real *8 (a-h,o-z)
        real *8 zs(2,1), zs2(2,1), zs3(2,1)
        real *8, allocatable :: xs(:),ys(:),xs2(:),ys2(:)
        real *8, allocatable :: xs3(:), ys3(:)
        character *1 title(1)
c
        allocate(xs(n), ys(n), xs2(n2), ys2(n2))
        allocate(xs3(n), ys3(n))
c
        do i=1,n
          xs(i)=zs(1,i)
          ys(i)=zs(2,i)
        enddo
c
        do i=1,n2
          xs2(i)=zs2(1,i)
          ys2(i)=zs2(2,i)
        enddo
c
        do i=1,n3
          xs3(i) = zs3(1,i)
          ys3(i) = zs3(2,i)
        enddo
c
        call pyplot3(iw, xs, ys, n, itype, xs2, ys2, n2, itype2,
     1      xs3, ys3, n3, itype3, title)
c
        return
        end
c
c
c
c
c
        subroutine pyplot(iw,x,y,n,itype,title)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1)
        character *1 a1,a10,file1(10),file11(9),title(1),file1p(10),
     1      temp(32)
        character *2 ls1
        character *9 file4
        character *10 file8
        character *32 title2
        character *10 file8p
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1        (file1p,file8p)

c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat'
        file1(5)=a10
        file1(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1

        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
        call quamesslen3(title,nchar)
c
        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue
c
        ls1='k.'
        if (itype .eq. 2) ls1='kx'
        if (itype .eq. 3) ls1='k-'
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/2,2)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
c
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
cccc        write(iun87, '(a)') 'pt.axes().set_aspect("equal")'
ccc        write(iun87,'(a)') 'pt.show()'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'

c
c       now print out the data file, plotiw.dat
c
        iun88=88
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i)
 1600 continue

c
        return
        end
c
c
c
c
c
        subroutine pyplot2(iw,x,y,n,itype,x2,y2,n2,itype2,title)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1),x2(1),y2(1)
        character *1 a1,a10,file1(11),file2(11),file11(9),
     1      title(1),temp(32)
        character *2 ls1,ls2
        character *9 file4
        character *11 file8,file9
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1      (file2,file9)
c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)
c
        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
        call quamesslen3(title,nchar)
c

        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue
c

        ls1='k.'
        if (itype .eq. 2) ls1='kx'
        if (itype .eq. 3) ls1='k-'
c
        ls2='g.'
        if (itype2 .eq. 2) ls2='gx'
        if (itype2 .eq. 3) ls2='g-'
c

        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/2,2)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a)') 'y = y.reshape(np.size(y)/2,2)'
        write(iun87,'(a)') 'y1=y[:,0]'
        write(iun87,'(a)') 'y2=y[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.plot(y1,y2,"',ls2,'")'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87, '(a)') 'pt.axes().set_aspect("equal")'
        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i)
 1600 continue
c
        iun89=899
        open(unit=iun89,file=file9)
        do 2600 i=1,n2
        write(iun89,*) x2(i), y2(i)
 2600 continue

c
        return
        end
c
c
c
c
c
        subroutine pyplot3(iw, x, y, n, itype, x2, y2, n2, itype2,
     1      x3, y3, n3, itype3, title)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1),x2(1),y2(1),x3(1),y3(1)
        character *1 a1,a10,file1(11),file2(11),file3(11),
     1      file11(9),title(1),temp(32),file1p(10)
        character *2 ls1,ls2,ls3
        character *9 file4
        character *11 file8,file9,file10
        character *32 title2
        character *10 file8p
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1      (file2,file9), (file3,file10), (file1p,file8p)
c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)
c
        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1
c
        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1

        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
        call quamesslen3(title,nchar)
c

        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue
c

        ls1='k.'
        if (itype .eq. 2) ls1='kx'
        if (itype .eq. 3) ls1='k-'
c
        ls2='g.'
        if (itype2 .eq. 2) ls2='gx'
        if (itype2 .eq. 3) ls2='g-'
c
        ls3='r.'
        if (itype3 .eq. 2) ls2='rx'
        if (itype3 .eq. 3) ls2='r-'
c

        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/2,2)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a)') 'y = y.reshape(np.size(y)/2,2)'
        write(iun87,'(a)') 'y1=y[:,0]'
        write(iun87,'(a)') 'y2=y[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'z = np.loadtxt("',file10,'")'
        write(iun87,'(a)') 'z = z.reshape(np.size(z)/2,2)'
        write(iun87,'(a)') 'z1=z[:,0]'
        write(iun87,'(a)') 'z2=z[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.plot(y1,y2,"',ls2,'")'
        write(iun87,'(a,a,a)') 'pt.plot(z1,z2,"',ls3,'")'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87, '(a)') 'pt.axes().set_aspect("equal")'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i)
 1600 continue
c
        iun89=899
        open(unit=iun89,file=file9)
        do 2600 i=1,n2
        write(iun89,*) x2(i), y2(i)
 2600 continue
c
        iun90=900
        open(unit=iun90,file=file10)
        do i=1,n3
          write(iun90,*) x3(i), y3(i)
        enddo

c
        return
        end
c
c
c
        subroutine pyplot3b(iw, x, y, n, itype, x2, y2, n2, itype2,
     1      x3, y3, n3, itype3,xmin,xmax,ymin,ymax)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1),x2(1),y2(1),x3(1),y3(1)
        character *1 a1,a10,file1(11),file2(11),file3(11),
     1      file11(9),title(1),temp(32),file1p(10)
        character *2 ls1,ls2,ls3
        character *9 file4
        character *11 file8,file9,file10
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1      (file2,file9), (file3,file10),(file1p,file8p)
c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)
c
        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.eps'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1
c
        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
c

        ls1='k.'
        if (itype .eq. 2) ls1='kx'
        if (itype .eq. 3) ls1='k-'
c
        ls2='g.'
        if (itype2 .eq. 2) ls2='gx'
        if (itype2 .eq. 3) ls2='g-'
c
        ls3='r.'
        if (itype3 .eq. 2) ls2='rx'
        if (itype3 .eq. 3) ls2='r-'
c

        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/2,2)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a)') 'y = y.reshape(np.size(y)/2,2)'
        write(iun87,'(a)') 'y1=y[:,0]'
        write(iun87,'(a)') 'y2=y[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'z = np.loadtxt("',file10,'")'
        write(iun87,'(a)') 'z = z.reshape(np.size(z)/2,2)'
        write(iun87,'(a)') 'z1=z[:,0]'
        write(iun87,'(a)') 'z2=z[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.plot(y1,y2,"',ls2,'")'
        write(iun87,'(a,a,a)') 'pt.plot(z1,z2,"',ls3,'")'
        write(iun87,'(a,e22.16,a,e22.16,a)')
     1    'pt.xlim((',xmin,',',xmax,'))'
        write(iun87,'(a,e22.16,a,e22.16,a)')
     1    'pt.ylim((',ymin,',',ymax,'))'
        write(iun87, '(a)') 'pt.axis("off")'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'",dpi=100)'
c        write(iun87,'(a)') 'pt.show()'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i)
 1600 continue
c
        iun89=899
        open(unit=iun89,file=file9)
        do 2600 i=1,n2
        write(iun89,*) x2(i), y2(i)
 2600 continue
c
        iun90=900
        open(unit=iun90,file=file10)
        do i=1,n3
          write(iun90,*) x3(i), y3(i)
        enddo

c
        return
        end
c
c
c
        subroutine pyplot4b(iw, x, y, n, itype, x2, y2, n2, itype2,
     1      x3, y3, n3, itype3, x4,y4,n4,itype4,xmin,xmax,ymin,ymax)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1),x2(1),y2(1),x3(1),y3(1),x4(1),y4(1)
        character *1 a1,a10,file1(11),file2(11),file3(11),file5(11),
     1      file11(9),title(1),temp(32),file1p(10)
        character *2 ls1,ls2,ls3,ls4
        character *9 file4
        character *11 file8,file9,file10,file12
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1      (file2,file9), (file3,file10),(file5,file12),
     2      (file1p,file8p)
c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)
c
        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.eps'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1
c
        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file12='plotiw.dat4'
        file5(5)=a10
        file5(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
c

        ls1='k.'
        if (itype .eq. 2) ls1='kx'
        if (itype .eq. 3) ls1='k-'
c
        ls2='g.'
        if (itype2 .eq. 2) ls2='gx'
        if (itype2 .eq. 3) ls2='g-'
c
        ls3='r.'
        if (itype3 .eq. 2) ls3='rx'
        if (itype3 .eq. 3) ls3='r-'
c
        ls4='c.'
        if (itype4 .eq. 2) ls4='cx'
        if (itype4 .eq. 3) ls4='c-'
c

        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/2,2)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a)') 'y = y.reshape(np.size(y)/2,2)'
        write(iun87,'(a)') 'y1=y[:,0]'
        write(iun87,'(a)') 'y2=y[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'z = np.loadtxt("',file10,'")'
        write(iun87,'(a)') 'z = z.reshape(np.size(z)/2,2)'
        write(iun87,'(a)') 'z1=z[:,0]'
        write(iun87,'(a)') 'z2=z[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'w = np.loadtxt("',file12,'")'
        write(iun87,'(a)') 'w = w.reshape(np.size(w)/2,2)'
        write(iun87,'(a)') 'w1=w[:,0]'
        write(iun87,'(a)') 'w2=w[:,1]'
c
        write(iun87,'(a)') ''
        write(iun87,'(a)') ''
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.plot(y1,y2,"',ls2,'")'
        write(iun87,'(a,a,a)') 'pt.plot(z1,z2,"',ls3,'")'
        write(iun87,'(a,a,a)') 'pt.plot(w1,w2,"',ls4,'")'
        write(iun87,'(a,e22.16,a,e22.16,a)')
     1    'pt.xlim((',xmin,',',xmax,'))'
        write(iun87,'(a,e22.16,a,e22.16,a)')
     1    'pt.ylim((',ymin,',',ymax,'))'
        write(iun87, '(a)') 'pt.axis("off")'
c        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'",dpi=100)'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i)
 1600 continue
c
        iun89=899
        open(unit=iun89,file=file9)
        do 2600 i=1,n2
        write(iun89,*) x2(i), y2(i)
 2600 continue
c
        iun90=900
        open(unit=iun90,file=file10)
        do i=1,n3
          write(iun90,*) x3(i), y3(i)
        enddo

        iun91=901
        open(unit=iun91,file=file12)
        do i=1,n4
          write(iun91,*) x4(i), y4(i)
        enddo

c
        return
        end
c
c
c
c
        subroutine pyimage(iw,m,n,vals,title)
        implicit real *8 (a-h,o-z)
        real *8 vals(m,n),target(2,m,n)
        character *1 a1,a10,file1(10),file11(9),title(1),
     1      temp(32), file1p(10)
        character *9 file4
        character *10 file8, file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1      (file1p,file8p)
c
c       output the data stored in vals to a file which python can
c       read and plot as a heat map using matplotlib
c
c       this routine requires a very special formatting of the inputs
c       target and vals - vals(1,1) is the value in the upper left
c       hand corner of the plot, vals(m,1) is the value in the lower
c       left hand corner, vals(1,n) is the value in the upper right
c       hand corner, and vals(m,n) is the value in the lower right
c       hand corner.
c
c       input:
c
c         iw - file number, same as in quaplot
c         m - number of y values on the grid
c         n - number of x values on the grid
c         vals - 'z' value at the x,y grid points
c         title - the title of the plot, should be of the form 'title*',
c             and can only be a maximum of 32 characters long
c
c       note that on exit, two files will be created - plotiw.py
c       and plotiw.dat
c
c
c       first construct the file names using iw
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, gniw
c
        iun87=877
        open(unit=iun87,file=file4)
c
        call quamesslen3(title,nchar)
c
        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
        write(iun87,'(a)') 'pt.rc("font", size=16)'
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a,i5,a,i5,a)')
     1      'x = x.reshape(', m, ',', n, ', order="F")'
        write(iun87,'(a)') 'pt.imshow(x, aspect="auto")'
        write(iun87,'(a)') 'cb=pt.colorbar(shrink=.9)'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 j=1,n
        do 1400 i=1,m
        write(iun88,*) vals(i,j)
 1400 continue
 1600 continue
c
        return
        end
c
c
c
        subroutine pyimage2(iw,m,n,vals,x,y,npts,x1,x2,
     1      y1,y2,title)
        implicit real *8 (a-h,o-z)
        real *8 vals(m,n),target(2,m,n),x(1),y(1)
        character *1 a1,a10,file1(11),file11(9),title(1),
     1      temp(32),file2(11),file1p(10)
        character *9 file4
        character *11 file8,file9
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2)
        equivalence (file2,file9), (file1p,file8p)
c
c       output the data stored in vals to a file which python can
c       read and plot as a heat map using matplotlib - this routine
c       also plots the scatterer, contains in x,y, as a filled
c       in polygon
c
c       this routine requires a very special formatting of the inputs
c       target and vals - vals(1,1) is the value in the upper left
c       hand corner of the plot, vals(m,1) is the value in the lower
c       left hand corner, vals(1,n) is the value in the upper right
c       hand corner, and vals(m,n) is the value in the lower right
c       hand corner.
c
c       input:
c
c         iw - file number, same as in quaplot
c         m - number of y values on the grid
c         n - number of x values on the grid
c         vals - 'z' value at the x,y grid points
c         x,y - points on the scatterer
c         npts - length of x and y
c         title - the title of the plot, should be of the form 'title*',
c             and can only be a maximum of 32 characters long
c
c       note that on exit, two files will be created - plotiw.py
c       and plotiw.dat
c
c
c       first construct the file names using iw
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, gniw
c
        iun87=877
        open(unit=iun87,file=file4)
c
c     Option 2 - local + global qbx
        call quamesslen3(title,nchar)
c
        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue

c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
        write(iun87,'(a)') 'pt.rc("font", size=16)'
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a,i5,a,i5,a)')
     1      'x = x.reshape(', m, ',', n, ', order="F")'
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,'])'
        write(iun87,'(a)') 'cb=pt.colorbar(shrink=.9)'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a)') 'pt.fill(y[:,0],y[:,1],facecolor="w")'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat1 and plotiw.dat2
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 j=1,n
        do 1400 i=1,m
        write(iun88,*) vals(i,j)
 1400 continue
 1600 continue
c
        call prinf('inside pyimage2, npts = *', npts, 1)
c
        iun89=889
        open(unit=iun89,file=file9)
        do 2600 i=1,npts
        write(iun89,*) x(i), y(i)
 2600 continue
c
        return
        end
c
c
c
c
        subroutine pyimage3(iw,m,n,vals,nd,x,y,npts,x1,x2,
     1      y1,y2)
        implicit real *8 (a-h,o-z)
        integer npts(1)
        real *8 vals(m,n),target(2,m,n),x(1),y(1)
        character *1 a1,a10,file1(11),file11(9),title(1),
     1      temp(32),file2(11),file1p(10),file3(11)
        character *9 file4
        character *11 file8,file9,file10
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2)
        equivalence (file2,file9), (file1p,file8p),(file3,file10)
c
c       output the data stored in vals to a file which python can
c       read and plot as a heat map using matplotlib - this routine
c       also plots the scatterer, contains in x,y, as a filled
c       in polygon
c
c       this routine requires a very special formatting of the inputs
c       target and vals - vals(1,1) is the value in the upper left
c       hand corner of the plot, vals(m,1) is the value in the lower
c       left hand corner, vals(1,n) is the value in the upper right
c       hand corner, and vals(m,n) is the value in the lower right
c       hand corner.
c
c       input:
c
c         iw - file number, same as in quaplot
c         m - number of y values on the grid
c         n - number of x values on the grid
c         vals - 'z' value at the x,y grid points
c         nd - number of scatterers
c         x,y - points on the scatterer
c         npts(nd) - length of x and y
c         x1,x2,y1,y2 - x and y limits for the plot
c
c       note that on exit, two files will be created - plotiw.py
c       and plotiw.dat
c
c
c       first construct the file names using iw
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1

        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, gniw
c
        iun87=877
        open(unit=iun87,file=file4)
c
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
        write(iun87,'(a)') 'pt.rc("font", size=16)'
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a,i5,a,i5,a)')
     1      'x = x.reshape(', m, ',', n, ', order="F")'
ccc        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a,a,a)') 'f=open("',file10,'","r")'
        write(iun87,'(a)') 'iss = 0'
        write(iun87,'(a,i5)') 'nd = ',nd
        write(iun87,'(a)') 'for i in range(nd):'
        write(iun87,'(4x,a)') 'nn = int(f.readline())'
        write(iun87,'(4x,a)')
     1   'pt.fill(y[iss:iss+nn,0],y[iss:iss+nn,1],facecolor="w")'
        write(iun87,'(4x,a)') 'iss=iss+nn'
        write(iun87,*)
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,'])'
        write(iun87,'(a)') 'cb=pt.colorbar(shrink=.9)'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat1 and plotiw.dat2
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 j=1,n
        do 1400 i=1,m
        write(iun88,*) vals(i,j)
 1400 continue
 1600 continue
c
        call prinf('inside pyimage2, npts = *', npts, 1)
c
c
c
        iun89=889
        open(unit=iun89,file=file9)
        iun90=890
        open(unit=iun90,file=file10)
        ntot = 0

        do 2550 i=1,nd

        ntot = ntot + npts(i)
        write(iun90,*) npts(i)

 2550   continue

ccc        call prinf('ntot=*',ntot,1)
ccc        call prin2('x=*',x,ntot)
ccc        call prin2('y=*',y,ntot)
ccc        call prinf('npts=*',npts,nd)

        do 2600 i=1,ntot
        write(iun89,*) x(i), y(i)
 2600 continue
c
        return
        end
c
c
c
c
c
        subroutine pyimage4(iw,m,n,vals,nd,x,y,npts,x1,x2,
     1      y1,y2,icmp)
        implicit real *8 (a-h,o-z)
        integer npts(1)
        real *8 vals(m,n),target(2,m,n),x(1),y(1)
        character *1 a1,a10,file1(11),file11(9),title(1),
     1      temp(32),file2(11),file1p(10),file3(11)
        character *9 file4
        character *11 file8,file9,file10
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2)
        equivalence (file2,file9), (file1p,file8p),(file3,file10)
c
c       output the data stored in vals to a file which python can
c       read and plot as a heat map using matplotlib - this routine
c       also plots the scatterer, contains in x,y, as a filled
c       in polygon
c
c       this routine requires a very special formatting of the inputs
c       target and vals - vals(1,1) is the value in the upper left
c       hand corner of the plot, vals(m,1) is the value in the lower
c       left hand corner, vals(1,n) is the value in the upper right
c       hand corner, and vals(m,n) is the value in the lower right
c       hand corner.
c
c       input:
c
c         iw - file number, same as in quaplot
c         m - number of y values on the grid
c         n - number of x values on the grid
c         vals - 'z' value at the x,y grid points
c         nd - number of scatterers
c         x,y - points on the scatterer
c         npts(nd) - length of x and y
c         x1,x2,y1,y2 - x and y limits for the plot
c
c       note that on exit, two files will be created - plotiw.py
c       and plotiw.dat
c
c
c       first construct the file names using iw
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1

        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, gniw
c
        iun87=877
        open(unit=iun87,file=file4)
c
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
	write(iun87,'(a)') 'import colormaps as cmaps'
        write(iun87,'(a)') ''
        write(iun87,'(a)') 'pt.rc("font", size=16)'
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a,i5,a,i5,a)')
     1      'x = x.reshape(', m, ',', n, ', order="F")'
ccc        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a,a,a)') 'f=open("',file10,'","r")'
        write(iun87,'(a)') 'iss = 0'
        write(iun87,'(a,i5)') 'nd = ',nd
        write(iun87,'(a)') 'for i in range(nd):'
        write(iun87,'(4x,a)') 'nn = int(f.readline())'
        write(iun87,'(4x,a)')
     1   'pt.fill(y[iss:iss+nn,0],y[iss:iss+nn,1],facecolor="w")'
        write(iun87,'(4x,a)') 'iss=iss+nn'
        write(iun87,*)

	if (icmp .eq. 1) then
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = "hot")'

	elseif(icmp .eq. 2) then
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.inferno)'

	elseif(icmp .eq. 3) then
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.plasma)'

	elseif(icmp .eq. 4) then
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.viridis)'

	elseif(icmp .eq. 5) then
        write(iun87,'(a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.parula)'
	endif

        write(iun87,'(a)') 'cb=pt.colorbar(shrink=.9)'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat1 and plotiw.dat2
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 j=1,n
        do 1400 i=1,m
        write(iun88,*) vals(i,j)
 1400 continue
 1600 continue
c
        call prinf('inside pyimage2, npts = *', npts, 1)
c
c
c
        iun89=889
        open(unit=iun89,file=file9)
        iun90=890
        open(unit=iun90,file=file10)
        ntot = 0

        do 2550 i=1,nd

        ntot = ntot + npts(i)
        write(iun90,*) npts(i)

 2550   continue

ccc        call prinf('ntot=*',ntot,1)
ccc        call prin2('x=*',x,ntot)
ccc        call prin2('y=*',y,ntot)
ccc        call prinf('npts=*',npts,nd)

        do 2600 i=1,ntot
        write(iun89,*) x(i), y(i)
 2600 continue
c
        return
        end
c--------------------------------------------------------
c
        subroutine pyimage5(iw,m,n,vals,nd,x,y,npts,x1,x2,
     1      y1,y2,icmp)
        implicit real *8 (a-h,o-z)
        integer npts(1)
        real *8 vals(m,n),target(2,m,n),x(1),y(1)
        character *1 a1,a10,file1(11),file11(9),title(1),
     1      temp(32),file2(11),file1p(10),file3(11)
        character *9 file4
        character *11 file8,file9,file10
        character *10 file8p
        character *32 title2
c
        equivalence (file1,file8), (file11,file4), (temp,title2)
        equivalence (file2,file9), (file1p,file8p),(file3,file10)
c
c       output the data stored in vals to a file which python can
c       read and plot as a heat map using matplotlib - this routine
c       also plots the scatterer, contains in x,y, as a filled
c       in polygon
c
c       this routine requires a very special formatting of the inputs
c       target and vals - vals(1,1) is the value in the upper left
c       hand corner of the plot, vals(m,1) is the value in the lower
c       left hand corner, vals(1,n) is the value in the upper right
c       hand corner, and vals(m,n) is the value in the lower right
c       hand corner.
c
c       input:
c
c         iw - file number, same as in quaplot
c         m - number of y values on the grid
c         n - number of x values on the grid
c         vals - 'z' value at the x,y grid points
c         nd - number of scatterers
c         x,y - points on the scatterer
c         npts(nd) - length of x and y
c         x1,x2,y1,y2 - x and y limits for the plot
c
c       note that on exit, two files will be created - plotiw.py
c       and plotiw.dat
c
c
c       first construct the file names using iw
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat1'
        file1(5)=a10
        file1(6)=a1
c
        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
        file9='plotiw.dat2'
        file2(5)=a10
        file2(6)=a1

        file10='plotiw.dat3'
        file3(5)=a10
        file3(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1
c
c       print out the contents of the scripting file, gniw
c
        iun87=877
        open(unit=iun87,file=file4)
c
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import numpy as np'
	write(iun87,'(a)') 'import colormaps as cmaps'
        write(iun87,'(a)') ''
        write(iun87,'(a)') 'pt.rc("font", size=16)'
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a,i5,a,i5,a)')
     1      'x = x.reshape(', m, ',', n, ', order="F")'
ccc        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
        write(iun87,'(a,a,a)') 'y = np.loadtxt("',file9,'")'
        write(iun87,'(a,a,a)') 'f=open("',file10,'","r")'
        write(iun87,'(a)') 'iss = 0'
        write(iun87,'(a,i5)') 'nd = ',nd
        write(iun87,'(a)') 'for i in range(nd):'
        write(iun87,'(4x,a)') 'nn = int(f.readline())'
        write(iun87,'(4x,a)')
     1   'pt.plot(y[iss:iss+nn,0],y[iss:iss+nn,1],"kx")'
        write(iun87,'(4x,a)') 'iss=iss+nn'
        write(iun87,*)

	if (icmp .eq. 1) then
        write(iun87,'(a,e11.5,a,e11.5,a,e11.5,a,e11.5,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = "hot")'

	elseif(icmp .eq. 2) then
        write(iun87,'(a,e11.5,a,e11.5,a,e11.5,a,e11.5,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.inferno)'

	elseif(icmp .eq. 3) then
        write(iun87,'(a,e11.5,a,e11.5,a,e11.5,a,e11.5,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.plasma)'

	elseif(icmp .eq. 4) then
        write(iun87,'(a,e11.5,a,e11.5,a,e11.5,a,e11.5,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.viridis)'

	elseif(icmp .eq. 5) then
        write(iun87,'(a,e11.5,a,e11.5,a,e11.5,a,e11.5,a)')
     1      'pt.imshow(x, extent=[',x1,',',x2,',',y1,',',y2,
     2      '],cmap = cmaps.parula)'
	endif

        write(iun87,'(a)') 'cb=pt.colorbar(shrink=.9)'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'
c        write(iun87,'(a)') 'pt.show()'

c
c       now print out the data file, plotiw.dat1 and plotiw.dat2
c
        iun88=888
        open(unit=iun88,file=file8)
        do 1600 j=1,n
        do 1400 i=1,m
        write(iun88,*) vals(i,j)
 1400 continue
 1600 continue
c
        call prinf('inside pyimage2, npts = *', npts, 1)
c
c
c
        iun89=889
        open(unit=iun89,file=file9)
        iun90=890
        open(unit=iun90,file=file10)
        ntot = 0

        do 2550 i=1,nd

        ntot = ntot + npts(i)
        write(iun90,*) npts(i)

 2550   continue

ccc        call prinf('ntot=*',ntot,1)
ccc        call prin2('x=*',x,ntot)
ccc        call prin2('y=*',y,ntot)
ccc        call prinf('npts=*',npts,nd)

        do 2600 i=1,ntot
        write(iun89,*) x(i), y(i)
 2600 continue
c
        return
        end
c
c
        SUBROUTINE quamesslen3(MES,nchar)
        CHARACTER *1 MES(1),AST
        DATA AST/'*'/
C
C         DETERMINE THE LENGTH OF THE MESSAGE
C
        I=0
        DO 1400 I=1,10000
        IF(MES(I).EQ.AST) GOTO 1600
        I1=I
 1400 CONTINUE
 1600 CONTINUE
c
        nchar=i1
c
        RETURN
        END
c
c
c
c
c
        subroutine int2char2(n,a)
        implicit real *8 (a-h,o-z)
        character *1 a
c
        if (n .eq. 0) a='0'
        if (n .eq. 1) a='1'
        if (n .eq. 2) a='2'
        if (n .eq. 3) a='3'
        if (n .eq. 4) a='4'
        if (n .eq. 5) a='5'
        if (n .eq. 6) a='6'
        if (n .eq. 7) a='7'
        if (n .eq. 8) a='8'
        if (n .eq. 9) a='9'
c
        return
        end

c
c
c
c
        subroutine pyplot3d(iw,x,y,z,n,itype,title)
        implicit real *8 (a-h,o-z)
        real *8 x(1),y(1),z(1)
        character *1 a1,a10,file1(10),file11(9),title(1),file1p(10),
     1      temp(32)
        character *2 ls1
        character *9 file4
        character *10 file8
        character *32 title2
        character *10 file8p
c
        equivalence (file1,file8), (file11,file4), (temp,title2),
     1        (file1p,file8p)

c
c       plots the points determined by x,y using itype style
c
        i1=mod(iw,10)
        i10=(iw-i1)/10
        call int2char2(i1,a1)
        call int2char2(i10,a10)

        file8='plotiw.dat'
        file1(5)=a10
        file1(6)=a1
c
        file4='plotiw.py'
        file11(5)=a10
        file11(6)=a1

        file8p='plotiw.pdf'
        file1p(5)=a10
        file1p(6)=a1
c
c
c       print out the contents of the scripting file, plotiw.py
c
        iun87=877
        open(unit=iun87,file=file4)
c
        call quamesslen3(title,nchar)
c
        title2=''
        do 1200 i=1,nchar
        temp(i)=title(i)
 1200 continue
c
        ls1='marker=x'
        if (itype .eq. 2) ls1='marker=o'
        if (itype .eq. 3) ls1='marker=^'
c
        write(iun87,'(a)') '#!/usr/bin/python'
        write(iun87,'(a)') 'import matplotlib.pyplot as pt'
        write(iun87,'(a)') 'import mpl_toolkits.mplot3d as Axes3D'
        write(iun87,'(a)') 'import numpy as np'
        write(iun87,'(a)') ''
c
        write(iun87,'(a,a,a)') 'x = np.loadtxt("',file8,'")'
        write(iun87,'(a)') 'x = x.reshape(np.size(x)/3,3)'
        write(iun87,'(a)') 'x1=x[:,0]'
        write(iun87,'(a)') 'x2=x[:,1]'
        write(iun87,'(a)') 'x3=x[:,2]'
c
        write(iun87,'(a)') 'fig=pt.figure()'
        write(iun87,'(a)') 'ax=fig.add_subplot(111,projection="3d")'
        write(iun87,'(a,a,a)') 'pt.plot(x1,x2,x3,"',ls1,'")'
        write(iun87,'(a,a,a)') 'pt.title("', title2, '")'
cccc        write(iun87, '(a)') 'pt.axes().set_aspect("equal")'
ccc        write(iun87,'(a)') 'pt.show()'
        write(iun87,'(a,a,a)') 'pt.savefig("',file8p,'")'

c
c       now print out the data file, plotiw.dat
c
        iun88=88
        open(unit=iun88,file=file8)
        do 1600 i=1,n
        write(iun88,*) x(i), y(i), z(i)
 1600 continue

c
        return
        end
c
c

max_stars_repo_name: askhamwhat/biharm-evals

max_issues_repo_name: askhamwhat/biharm-evals

max_forks_repo_name: askhamwhat/biharm-evals


text:
! --------------------------------------------------------------------
!
! Copyright (C) 2015 Rocco Meli
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! ---------------------------------------------------------------------

PROGRAM classicalMD

    USE DYNAMICS, only: classical_step
    USE OUTPUT

    IMPLICIT NONE

    INTEGER, PARAMETER :: steps = 1000
    REAL*8, PARAMETER:: dt = 1e-3

    INTEGER :: i ! Loop index

    INTEGER, PARAMETER :: Nn = 3
    CHARACTER(len=2), dimension(Nn) :: atoms
    REAL*8, dimension(Nn,3) :: pos
    REAL*8, dimension(Nn,3) :: vel
    REAL*8, dimension(Nn,3) :: force
    REAL*8, dimension(Nn) :: mass

    atoms = (/"H","H"/)

    pos(1,:) = (/0.0D0, 0.0D0, 0.0D0/)
    pos(2,:) = (/1.2*1.0D0, 0.0D0, 0.0D0/)
    pos(3,:) = (/0.5*1.0D0, 1.3*1.0D0, 0.0D0/)

    vel(1,:) = (/0.0D0, 0.0D0, 0.0D0/)
    vel(2,:) = (/0.0D0, 0.0D0, 0.0D0/)
    vel(3,:) = (/0.0D0, 0.0D0, 0.0D0/)

    force = (-1) * dU(Nn,pos) ! Compute the initial force

    mass = (/1.0D0, 1.0D0, 10.0D0/)

    OPEN(unit=123,file="classicalMD.xyz",form="formatted",status="new",action="write") ! Open file 123

    DO i = 1, steps
        CALL append_xyz(Nn,atoms,pos,123)
        CALL classical_step(Nn,mass,pos,vel,force,dU,dt)
    END DO

    CLOSE(unit=123) ! Close file 123

    CONTAINS
        FUNCTION dU(Nn,pos)

            ! INPUT
            INTEGER, intent(in) :: Nn                   ! Number of particles
            REAL*8, dimension(Nn,3), intent(in) :: pos  ! Particle's position

            ! OUTPUT
            REAL*8, dimension(Nn,3) :: dU               ! Gradient of the PES

            ! PARAMETERS
            REAL*8, PARAMETER :: k = 500 ! Spring constant (TODO: different force constant for different interactions)
            REAL*8, PARAMETER :: d = 1 ! Equilibrium distance (TODO: different equilibrium distance for different interactions)

            ! INTERMEDIATE VARIABLES
            REAL*8 :: r, dx, dy, dz
            INTEGER :: i, j ! Loop indices

            dU(:,:) = 0.0D0

            DO i = 1, Nn-1
                DO j = i+1, Nn
                    dx = pos(i,1) - pos(j,1)
                    dy = pos(i,2) - pos(j,2)
                    dz = pos(i,3) - pos(j,3)

                    r = DSQRT(dx**2 + dy**2 + dz**2)

                    dU(i,1) = dU(i,1) - k * (d - r) / r * dx
                    dU(i,2) = dU(i,2) - k * (d - r) / r * dy
                    dU(i,3) = dU(i,3) - k * (d - r) / r * dz

                    dU(j,1) = dU(j,1) + k * (d - r) / r * dx
                    dU(j,2) = dU(j,2) + k * (d - r) / r * dy
                    dU(j,3) = dU(j,3) + k * (d - r) / r * dz
                END DO
            END DO

        END FUNCTION

END PROGRAM classicalMD

max_stars_repo_name: RMeli/Hartree-Fock

max_issues_repo_name: RMeli/Hartree-Fock

max_forks_repo_name: RMeli/Hartree-Fock


!-------------------------------------------------------------------------------
!> module ATMOSPHERE / Physics Cloud Microphysics
!!
!! @par Description
!!          Cloud Microphysics wrapper
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2013-12-06 (S.Nishizawa) [new]
!!
!<
!-------------------------------------------------------------------------------
module mod_atmos_phy_mp_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ATMOS_PHY_MP_driver_setup
  public :: ATMOS_PHY_MP_driver

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
!  real(RP), private, allocatable :: DENS_t(:,:,:)
!  real(RP), private, allocatable :: MOMZ_t(:,:,:)
!  real(RP), private, allocatable :: MOMX_t(:,:,:)
!  real(RP), private, allocatable :: MOMY_t(:,:,:)
!  real(RP), private, allocatable :: RHOT_t(:,:,:)
!  real(RP), private, allocatable :: QTRC_t(:,:,:,:)
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup Cloud Microphysics
  !-----------------------------------------------------------------------------
  subroutine ATMOS_PHY_MP_driver_setup( MP_TYPE )
    use scale_process, only: &
       PRC_MPIstop
    use scale_atmos_phy_mp, only: &
       ATMOS_PHY_MP_setup
    use mod_atmos_vars, only: &
       DENS, &
       RHOT, &
       QTRC
    implicit none

    character(len=H_SHORT), intent(in) :: MP_TYPE
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '+++ Module[Cloud Microphisics]/Categ[ATMOS]'

!    allocate( DENS_t(KA,IA,JA) )
!    allocate( MOMZ_t(KA,IA,JA) )
!    allocate( MOMX_t(KA,IA,JA) )
!    allocate( MOMY_t(KA,IA,JA) )
!    allocate( RHOT_t(KA,IA,JA) )
!    allocate( QTRC_t(KA,IA,JA,QA) )

    call ATMOS_PHY_MP_setup( MP_TYPE )

!    call ATMOS_PHY_MP_driver( .true., .false. )

    return
  end subroutine ATMOS_PHY_MP_driver_setup

  !-----------------------------------------------------------------------------
  !> Cloud Microphysics
  !-----------------------------------------------------------------------------
  subroutine ATMOS_PHY_MP_driver( update_flag, history_flag )
    use scale_atmos_phy_mp, only: &
       ATMOS_PHY_MP
    use mod_atmos_vars, only: &
       ATMOS_vars_fillhalo, &
       ATMOS_vars_total, &
       DENS, &
       MOMZ, &
       MOMX, &
       MOMY, &
       RHOT, &
       QTRC
    implicit none
    logical, intent(in) :: update_flag
    logical, intent(in) :: history_flag

    if ( update_flag ) then
       call ATMOS_PHY_MP( &
            DENS, &
            MOMZ, &
            MOMX, &
            MOMY, &
            RHOT, &
            QTRC )

       call ATMOS_vars_fillhalo

       call ATMOS_vars_total
    end if

    return
  end subroutine ATMOS_PHY_MP_driver

end module mod_atmos_phy_mp_driver


text:
      SUBROUTINE HIRECT(IM)
*
*
*       Rectification of hierarchical binary.
*       -------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
*
*
*       Include diagnostic error check.
      RB = 0.0
      UPR2 = 0.0
      TD2 = 0.0
      DO 5 K = 1,4
          RB = RB + UM(K,IM)**2
          UPR2 = UPR2 + UMDOT(K,IM)**2
          TD2 = TD2 + 2.0*UM(K,IM)*UMDOT(K,IM)
    5 CONTINUE
*
      ZMB = CM(1,IM) + CM(2,IM)
      HI = (2.0*UPR2 - ZMB)/RB
      ERR = (HI - HM(IM))/HI
      ZMU = CM(1,IM)*CM(2,IM)/ZMB
      DB = ZMU*(HI - HM(IM))
      IF (ABS(DB).GT.1.0D-08) THEN
      SEMI = -0.5*ZMB/HM(IM)
      ECC2 = (1.0 - RB/SEMI)**2 + TD2**2/(ZMB*SEMI)
      ECC = SQRT(ECC2)
      RA = RB/SEMI
      WRITE (16,3) TTOT, NAMEG(IM), KSTARM(IM), ECC, RA, HM(IM), DB, ERR
    3 FORMAT (' HIRECT:   T NM K* E R/A H DB DH/H ',
     &                    F9.3,I6,I4,F8.4,F8.4,F9.2,1P,2E10.1)
      CALL FLUSH(16)
      END IF
*       Initialize iteration counter for difficult case (SJA 10/97).
      ITER = 0
*
*       Form square regularized velocity for the explicit binding energy.
   10 RB = 0.0D0
      UPR2 = 0.0
      DO 15 K = 1,4
          RB = RB + UM(K,IM)**2
          UPR2 = UPR2 + UMDOT(K,IM)**2
   15 CONTINUE
*
*       Form KS scaling factors from energy and angular momentum relation.
      A1 = 0.25D0*ZMB/UPR2
*       Solve for C1 from H = (2*U'*U'*C1**2 - M)/(U*U*C2**2) with C2 = 1/C1.
      A2 = A1**2 + 0.5D0*HM(IM)*RB/UPR2
*
*       Avoid negative round-off value on second try (NB! no change in CK).
      IF (ITER.EQ.2.AND.A2.LT.0.0) A2 = 0.0D0
*
*       Check for undefined case (circular orbit or eccentric anomaly = 90).
      IF (A2.GE.0.0D0) THEN
          IF (A1.LT.1.0) THEN
*       Choose square root sign from eccentric anomaly (e*cos(E) = 1 - R/a).
              C1 = SQRT(A1 + SQRT(A2))
          ELSE
              C1 = SQRT(A1 - SQRT(A2))
          END IF
          CK = 1.0
      ELSE
*       Adopt C1*C2 = CK for difficult case (Seppo's suggestion of 1991).
          C1 = 1.0
          CK = ZMB/SQRT(-8.0D0*HM(IM)*RB*UPR2)
          WRITE (6,20)  IM, KSTARM(IM), RB, HM(IM), UPR2, A2, CK-1.0
   20     FORMAT (' WARNING!    HIRECT    IM K* R H UPR2 A2 CK-1 ',
     &                                    2I4,1P,5E10.2)
          ITER = ITER + 1
      END IF
*
*       Specify KS coordinate scaling from angular momentum conservation.
      C2 = CK/C1
*
*       Transform KS variables to yield the prescribed elements.
      DO 25 K = 1,4
          UM(K,IM) = C2*UM(K,IM)
          UMDOT(K,IM) = C1*UMDOT(K,IM)
   25 CONTINUE
*
*       Improve solution by second iteration in case of CK procedure.
      ITER = ITER + 1
      IF (ITER.EQ.2) GO TO 10
*
      RETURN
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
We are (or were) in Davis because our significant other (boyfriend, girlfriend, husband, wife, life partner, etc.) is a student or faculty at UC Davis.

   Users/BillKendrick
   Users/HollywoOd
   Users/PenguiN42
   Users/NicoleBush
   Users/SaraGoldband
   Users/DanielHoang
   Users/DorrieSiobhan
   Users/condemned2bfree
   Users/bear
   Users/JabberWokky was living in Sacramento when he met the woman who is Users/SarahEdwards now his wife while doing Rocky Horror in Davis.  They moved in together in Davis.
   Users/mheathershaw

Feel free to add yourself if you are (or were) trapped in this (admittedly nice) little town thanks to someone else.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C (C) Copyright 1996-2016 ECMWF.
C
C This software is licensed under the terms of the Apache Licence Version 2.0
C which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation nor
C does it submit to any jurisdiction.
C

      program  magics

C Defines the colours table for the shading
      character*25 colour
      dimension colour(24)

C Defines the lines for the title
      character*100 title
      dimension title(4)

      data colour /
     +    'rgb(0,0,0.5)',
     +    'rgb(0,0,0.5)',
     +    'rgb(0,0,0.85)',
     +    'rgb(0,0,0.5)',
     +    'rgb(0,0,0.85)',
     +    'rgb(0.25,0,1)',
     +    'blue_purple',
     +    'greenish_blue',
     +    'blue_green',
     +    'bluish_green',
     +    'yellow_green',
     +    'greenish_yellow',
     +    'yellow',
     +    'orangish_yellow',
     +    'orange_yellow',
     +    'yellowish_orange',
     +    'orange',
     +    'reddish_orange',
     +    'red_orange',
     +    'orangish_red',
     +    'red',
     +    'magenta',
     +    'magenta',
     +    'magenta' /

      data title /
     +    "<grib_info id='z500' key='base-date' ",
     +    "<grib_info id='z500' key='valid-date' ",
     +    "parameter:<grib_info id='z500' key='name'/>",
     +    "parameter:<grib_info id='t850' key='name'/>"/

C open magics
      call popen

C Setting of the output file name
      call ???('???', '???')

C Setting the coordinates of the geographical area
C  For example take Europe
      call psetr
      call psetr
	  .....
      call psetc('subpage_map_projection', '???)

C Import the t850 data
      call psetc('grib_id', 't850')
      call psetc('grib_input_file_name', '../t850.grb')
      call ???

C Define the shading for t850
C find the parameters to set
	  call psetc...
	  call psetr...
	  call pseti...
      call pcont

C Import the z500 data
      call psetc('???', 'z500')
      call psetc('????', '../z500.grb')
      call pgrib

C Define the contouring for z500
C find the parameters to set
	  call psetc...
	  call psetr...
	  call pseti...
      call pcont

C Add the coastlines
      call psetc('map_grid_colour', 'tan')
      call psetc('map_grid', 'on')
      call psetc('map_coastline_colour', 'tan')
      call pcoast

C add a title
      call pset1c('text_lines', title, 4)
      call psetc('text_justification', 'left')
      call psetc('text_html', 'true')
      call psetc('text_colour', 'black')
      call ptext

C In a second step: try to position text whereever you want
C     call psetc('text_mode', 'positional')
C     call psetr('text_box_y_length', ???)
C     call psetc('text_border', '???')
C     call psetr('text_box_x_length', ???)
C     call psetr('text_box_x_position', ???)
C     call psetc('text_box_blanking', 'on')
C     call psetc('text_border_colour', 'black')
C     call psetr('text_box_y_position', ???)

C Add a legend
      call psetc('legend', 'on')
      call psetc('legend_display_type', 'continuous')

C In a second step:  the Try to move it around!
C     call psetr('legend_box_x_position', ???)
C     call psetr('legend_box_y_position', ??)
C     call psetr('legend_box_x_length', ???)
C     call psetc('legend_box_mode', 'positional')
C     call psetr('legend_box_y_length', ???)

      call pclose

      end


! Copyright 2016 The Go Authors. All rights reserved.
! Use of this source code is governed by a BSD-style
! license that can be found in the LICENSE file.

function the_answer() result(j) bind(C)
  use iso_c_binding, only: c_int
  integer(c_int) :: j ! output
  j = 42
end function the_answer


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: vector_field_grad.f90 8157 2014-10-09 13:02:44Z sfilippo $
!
! Description:
!    Computes the gradient of a VectorField object.
!    The result is a VECTOR object with rank 1.
!
!    REMARK: remember to allocate the LHS before assigning the function result
!
SUBMODULE(op_grad) vector_field_grad_implementation

    IMPLICIT NONE

    CONTAINS

        MODULE PROCEDURE vector_field_grad
            USE class_psblas
            USE class_connectivity
            USE class_face
            USE class_least_squares
            USE class_mesh
            USE class_vector_field
            USE class_vector

            IMPLICIT NONE
            !
            INTEGER :: i, ib, ibf, ic, icoo, IF, im, info, ib_offset
            INTEGER :: n, nb, nbc, ncd, ncells
            INTEGER, POINTER :: ic2c(:) => NULL(), if2b(:) => NULL()
            REAL(psb_dpk_) :: xc, yc, zc, xf, yf, zf, xnb, ynb, znb
            REAL(psb_dpk_), ALLOCATABLE :: phi_x(:,:)
            REAL(psb_dpk_), ALLOCATABLE :: phi_b(:,:)
            REAL(psb_dpk_), ALLOCATABLE :: rhs(:,:,:)
            TYPE(mesh), POINTER :: msh => NULL()

            ! Gets internal and boundary values
            CALL phi%get_x(phi_x)
            CALL phi%get_bx(phi_b)

            ! Face-centered fields are not supported
            IF(phi%on_faces_()) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            ! Ncells = local + halos
            ncells = SIZE(phi_x,1)

            IF(.NOT.ALLOCATED(grad)) THEN
                ALLOCATE(grad(3,ncells),stat=info)
                IF(info /= 0) THEN
                    WRITE(*,200)
                    CALL abort_psblas
                END IF
            END IF
            ! REMARK: compiler-dependent behaviour in SCALAR_FIELD_GRAD
            ! Intel:   allocated(grad) = F
            ! Gfortran: allocated(grad) = T

            ! Gets PHI mesh
            CALL phi%get_mesh(msh)
            ncd = msh%ncd

            ! Allocates RHS for least squares regression
            ALLOCATE(rhs(ncd+1,3,ncells),stat=info)
            IF(info /= 0) THEN
                WRITE(*,200)
                CALL abort_psblas
            END IF

            ! Initialization
            grad(:,:) = vector_(0.d0,0.d0,0.d0)

            ! Next loops sweep only strictly local cells => NCELLS redefinition.
            ncells = psb_cd_get_local_rows(msh%desc_c)
            nbc = msh%nbc

            rhs = 0.d0

            ! ----- Builds RHS -----

            DO ic = 1, ncells

                ! --- IC ---
                xc = msh%cell_cntr(ic)%x_()
                yc = msh%cell_cntr(ic)%y_()

                DO icoo = 1, 3
                    rhs(1,icoo,ic) = phi_x(ic,icoo)
                    rhs(2,icoo,ic) = phi_x(ic,icoo) * xc
                    rhs(3,icoo,ic) = phi_x(ic,icoo) * yc
                END DO

                ! --- IC's neighbors ---
                CALL msh%c2c%get_ith_conn(ic2c,ic)
                n = SIZE(ic2c)

                DO i = 1, n
                    nb = ic2c(i)
                    xnb = msh%cell_cntr(nb)%x_()
                    ynb = msh%cell_cntr(nb)%y_()

                    DO icoo = 1, 3
                        rhs(1,icoo,ic) = rhs(1,icoo,ic) + phi_x(nb,icoo)
                        rhs(2,icoo,ic) = rhs(2,icoo,ic) + phi_x(nb,icoo) * xnb
                        rhs(3,icoo,ic) = rhs(3,icoo,ic) + phi_x(nb,icoo) * ynb
                    END DO
                END DO

            END DO

            ! Center(s) of possible boundary faces
            DO ib = 1, nbc
                CALL msh%f2b%get_ith_conn(if2b,ib)
                n = SIZE(if2b)
                ib_offset = COUNT(msh%faces%flag_() > 0 .AND. msh%faces%flag_() < ib)

                DO i = 1, n
                    IF = if2b(i)
                    ibf = ib_offset + i
                    im = msh%faces(IF)%master_()

                    xf = msh%face_cntr(IF)%x_()
                    yf = msh%face_cntr(IF)%y_()

                    DO icoo = 1, 3
                        rhs(1,icoo,im) = rhs(1,icoo,im) + phi_b(ibf,icoo)
                        rhs(2,icoo,im) = rhs(2,icoo,im) + phi_b(ibf,icoo) * xf
                        rhs(3,icoo,im) = rhs(3,icoo,im) + phi_b(ibf,icoo) * yf
                    END DO
                END DO
            END DO

            ! Extra elements for 3D case only
            IF(ncd == 3) THEN

                DO ic = 1, ncells

                    ! --- IC ---
                    zc = msh%cell_cntr(ic)%z_()

                    DO icoo = 1, 3
                        rhs(4,icoo,ic) = phi_x(ic,icoo) * zc
                    END DO

                    ! --- IC's neighbors ---
                    CALL msh%c2c%get_ith_conn(ic2c,ic)
                    n = SIZE(ic2c)

                    DO i = 1, n
                        nb = ic2c(i)
                        znb = msh%cell_cntr(nb)%z_()

                        DO icoo = 1, 3
                            rhs(4,icoo,ic) = rhs(4,icoo,ic) + phi_x(nb,icoo) * znb
                        END DO
                    END DO
                END DO

                ! Center(s) of possible boundary faces
                DO ib = 1, nbc
                    CALL msh%f2b%get_ith_conn(if2b,ib)
                    n = SIZE(if2b)
                    ib_offset = COUNT(msh%faces%flag_() > 0 .AND. msh%faces%flag_() < ib)

                    DO i = 1, n
                        IF = if2b(i)
                        ibf = ib_offset + i
                        im = msh%faces(IF)%master_()

                        zf = msh%face_cntr(IF)%z_()

                        DO icoo = 1, 3
                            rhs(4,icoo,im) = rhs(4,icoo,im) + phi_b(ibf,icoo) * zf
                        END DO
                    END DO
                END DO

            END IF

            ! Solves least squares problem
            DO ic = 1, ncells
                DO icoo = 1, 3
                    CALL msh%lsr(ic)%solve_least_squares(rhs(:,icoo,ic))
                END DO
            END DO

            ! Sets grad values
            SELECT CASE(ncd)
            CASE(2)
                DO ic = 1, ncells
                    DO icoo = 1, 3
                        grad(icoo,ic) = vector_(rhs(2,icoo,ic),rhs(3,icoo,ic),0.d0)
                    END DO
                END DO
            CASE(3)
                DO ic = 1, ncells
                    DO icoo = 1, 3
                        grad(icoo,ic) = vector_(rhs(2,icoo,ic),rhs(3,icoo,ic),rhs(4,icoo,ic))
                    END DO
                END DO
            END SELECT

            ! Updates halo elements
            DO icoo = 1, 3
                CALL update_vector_halo(grad(icoo,:),msh%desc_c)
            END DO

            NULLIFY(if2b)
            NULLIFY(ic2c)
            DEALLOCATE(rhs)
            NULLIFY(msh)
            DEALLOCATE(phi_b,phi_x)

        100 FORMAT(' ERROR! Face-centered fields not supported by VECTOR_FIELD_GRAD')
        200 FORMAT(' ERROR! Memory allocation failure in VECTOR_FIELD_GRAD')

        END PROCEDURE vector_field_grad

END SUBMODULE vector_field_grad_implementation

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


! { dg-do compile }

program test_failed_images_1
  implicit none

  integer, allocatable :: fi(:)
  real :: r
  integer :: i

  fi = failed_images()         ! OK
  fi = failed_images(TEAM=1)   ! { dg-error "'team' argument of 'failed_images' intrinsic at \\(1\\) not yet supported" }
  fi = failed_images(KIND=1)   ! OK
  fi = failed_images(KIND=4)   ! OK
  fi = failed_images(KIND=0)   ! { dg-error "'kind' argument of 'failed_images' intrinsic at \\\(1\\\) must be positive" }
  fi = failed_images(KIND=r)   ! { dg-error "'kind' argument of 'failed_images' intrinsic at \\\(1\\\) must be INTEGER" }
  fi = failed_images(KIND=i)   ! { dg-error "Constant expression required at \\\(1\\\)" }
  fi = failed_images(KIND=42)  ! { dg-error "'kind' argument of 'failed_images' intrinsic at \\\(1\\\) shall specify a valid integer kind" }

end program test_failed_images_1


text:
!-----------------------------------------------------------------------
!
! Header with QFunctions
!
      include 't502-operator-f.h'
!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceed/fortran.h'

      integer ceed,err,i,j
      integer stridesu_small(3),stridesu_large(3)
      integer erestrictx,erestrictu
      integer erestrictui_small,erestrictui_large
      integer bx_small,bu_small,bx_large,bu_large
      integer qf_setup,qf_mass
      integer op_setup_small,op_mass_small,op_setup_large,op_mass_large
      integer qdata_small,qdata_large,x,u,v
      integer nelem,p,q,scale
      parameter(nelem=15)
      parameter(p=5)
      parameter(q=8)
      parameter(scale=3)
      integer nx,nu
      parameter(nx=nelem+1)
      parameter(nu=nelem*(p-1)+1)
      integer indx(nelem*2)
      integer indu(nelem*p)
      real*8 arrx(nx)
      integer*8 voffset,xoffset

      real*8 hu(nu*2),hv(nu*2)
      real*8 total1,total2

      character arg*32

      external setup,mass

      call getarg(1,arg)
      call ceedinit(trim(arg)//char(0),ceed,err)

      do i=0,nx-1
        arrx(i+1)=i/(nx-1.d0)
      enddo
      do i=0,nelem-1
        indx(2*i+1)=i
        indx(2*i+2)=i+1
      enddo

      call ceedelemrestrictioncreate(ceed,nelem,2,1,1,nx,ceed_mem_host,&
     & ceed_use_pointer,indx,erestrictx,err)

      do i=0,nelem-1
        do j=0,p-1
          indu(p*i+j+1)=2*(i*(p-1)+j)
        enddo
      enddo

      call ceedelemrestrictioncreate(ceed,nelem,p,2,1,2*nu,ceed_mem_host,&
     & ceed_use_pointer,indu,erestrictu,err)
      stridesu_small=[1,q,q]
      call ceedelemrestrictioncreatestrided(ceed,nelem,q,1,q*nelem,&
     & stridesu_small,erestrictui_small,err)
      stridesu_large=[1,q*scale,q*scale]
      call ceedelemrestrictioncreatestrided(ceed,nelem,q*scale,1,&
     & q*nelem*scale,stridesu_large,erestrictui_large,err)

      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx_small,err)
      call ceedbasiscreatetensorh1lagrange(ceed,1,2,p,q,ceed_gauss,bu_small,err)
      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q*scale,&
     & ceed_gauss,bx_large,err)
      call ceedbasiscreatetensorh1lagrange(ceed,1,2,p,q*scale,&
     & ceed_gauss,bu_large,err)

! Common QFunctions

      call ceedqfunctioncreateinterior(ceed,1,setup,&
     &SOURCE_DIR&
     &//'t502-operator.h:setup'//char(0),qf_setup,err)
      call ceedqfunctionaddinput(qf_setup,'weight',1,ceed_eval_weight,err)
      call ceedqfunctionaddinput(qf_setup,'x',1,ceed_eval_grad,err)
      call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err)

      call ceedqfunctioncreateinterior(ceed,1,mass,&
     &SOURCE_DIR&
     &//'t502-operator.h:mass'//char(0),qf_mass,err)
      call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err)
      call ceedqfunctionaddinput(qf_mass,'u',2,ceed_eval_interp,err)
      call ceedqfunctionaddoutput(qf_mass,'v',2,ceed_eval_interp,err)

      call ceedvectorcreate(ceed,nx,x,err)
      xoffset=0
      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)

! Small operator

      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
     & ceed_qfunction_none,op_setup_small,err)
      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
     & ceed_qfunction_none,op_mass_small,err)

      call ceedvectorcreate(ceed,nelem*q,qdata_small,err)

      call ceedoperatorsetfield(op_setup_small,'weight',&
     & ceed_elemrestriction_none,bx_small,ceed_vector_none,err)
      call ceedoperatorsetfield(op_setup_small,'x',erestrictx,&
     & bx_small,ceed_vector_active,err)
      call ceedoperatorsetfield(op_setup_small,'rho',erestrictui_small,&
     & ceed_basis_collocated,ceed_vector_active,err)
      call ceedoperatorsetfield(op_mass_small,'rho',erestrictui_small,&
     & ceed_basis_collocated,qdata_small,err)
      call ceedoperatorsetfield(op_mass_small,'u',erestrictu,&
     & bu_small,ceed_vector_active,err)
      call ceedoperatorsetfield(op_mass_small,'v',erestrictu,&
     & bu_small,ceed_vector_active,err)

! Large operator

      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
     & ceed_qfunction_none,op_setup_large,err)
      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
     & ceed_qfunction_none,op_mass_large,err)

      call ceedvectorcreate(ceed,nelem*q*scale,qdata_large,err)

      call ceedoperatorsetfield(op_setup_large,'weight',&
     & ceed_elemrestriction_none,bx_large,ceed_vector_none,err)
      call ceedoperatorsetfield(op_setup_large,'x',erestrictx,&
     & bx_large,ceed_vector_active,err)
      call ceedoperatorsetfield(op_setup_large,'rho',erestrictui_large,&
     & ceed_basis_collocated,ceed_vector_active,err)
      call ceedoperatorsetfield(op_mass_large,'rho',erestrictui_large,&
     & ceed_basis_collocated,qdata_large,err)
      call ceedoperatorsetfield(op_mass_large,'u',erestrictu,&
     & bu_large,ceed_vector_active,err)
      call ceedoperatorsetfield(op_mass_large,'v',erestrictu,&
     & bu_large,ceed_vector_active,err)

! Setup U, V

      call ceedvectorcreate(ceed,2*nu,u,err)
      call ceedvectorgetarraywrite(u,ceed_mem_host,hu,voffset,err)
      do i=1,nu
        hu(voffset+2*i-1)=1.
        hu(voffset+2*i)=2.
      enddo
      call ceedvectorrestorearray(u,hu,voffset,err)
      call ceedvectorcreate(ceed,2*nu,v,err)

! Small apply

      call ceedoperatorapply(op_setup_small,x,qdata_small,&
     & ceed_request_immediate,err)
      call ceedoperatorapply(op_mass_small,u,v,ceed_request_immediate,err)

      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
      total1=0.
      total2=0.
      do i=1,nu
        total1=total1+hv(voffset+2*i-1)
        total2=total2+hv(voffset+2*i)
      enddo
      if (abs(total1-1.)>1.0d-10) then
! LCOV_EXCL_START
        write(*,*) 'Computed Area: ',total1,' != True Area: 1.0'
! LCOV_EXCL_STOP
      endif
      if (abs(total2-2.)>1.0d-10) then
! LCOV_EXCL_START
        write(*,*) 'Computed Area: ',total2,' != True Area: 2.0'
! LCOV_EXCL_STOP
      endif
      call ceedvectorrestorearrayread(v,hv,voffset,err)

! Large apply

      call ceedoperatorapply(op_setup_large,x,qdata_large,&
     & ceed_request_immediate,err)
      call ceedoperatorapply(op_mass_large,u,v,ceed_request_immediate,err)

      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
      total1=0.
      total2=0.
      do i=1,nu
        total1=total1+hv(voffset+2*i-1)
        total2=total2+hv(voffset+2*i)
      enddo
      if (abs(total1-1.)>1.0d-10) then
! LCOV_EXCL_START
        write(*,*) 'Computed Area: ',total1,' != True Area: 1.0'
! LCOV_EXCL_STOP
      endif
      if (abs(total2-2.)>1.0d-10) then
! LCOV_EXCL_START
        write(*,*) 'Computed Area: ',total2,' != True Area: 2.0'
! LCOV_EXCL_STOP
      endif
      call ceedvectorrestorearrayread(v,hv,voffset,err)

      call ceedvectordestroy(qdata_small,err)
      call ceedvectordestroy(qdata_large,err)
      call ceedvectordestroy(x,err)
      call ceedvectordestroy(u,err)
      call ceedvectordestroy(v,err)
      call ceedoperatordestroy(op_mass_small,err)
      call ceedoperatordestroy(op_setup_small,err)
      call ceedoperatordestroy(op_mass_large,err)
      call ceedoperatordestroy(op_setup_large,err)
      call ceedqfunctiondestroy(qf_mass,err)
      call ceedqfunctiondestroy(qf_setup,err)
      call ceedbasisdestroy(bu_small,err)
      call ceedbasisdestroy(bx_small,err)
      call ceedbasisdestroy(bu_large,err)
      call ceedbasisdestroy(bx_large,err)
      call ceedelemrestrictiondestroy(erestrictu,err)
      call ceedelemrestrictiondestroy(erestrictx,err)
      call ceedelemrestrictiondestroy(erestrictui_small,err)
      call ceedelemrestrictiondestroy(erestrictui_large,err)
      call ceeddestroy(ceed,err)
      end
!-----------------------------------------------------------------------

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
!> \brief Computes a Diagonal Lumped Boundary Mass matrix for the Current Mesh. <BR>
!> The Matrix is stored inverted (1/mass), so that no divisions are required when applying it
!! @param bvmass Boundary Inverted Mass matrix
subroutine nsc_CalculateBoundaryVmass(a)
   use typre
   use def_parame
   use Mod_Memor
   use Mod_Mesh
   use Mod_Element
   use Mod_NSCompressible

   implicit none
   class(NSCompressibleProblem) :: a

   integer(ip)             :: iboun,nboun,npoin
   integer(ip)             :: ipoin,ndime,igaub
   integer(ip)             :: calculate_statistics

   class(FiniteElement), pointer :: e => NULL()
   real(rp)    :: auxshape(a%Mesh%mnode)
   real(rp)    :: dsurf

   !Initializations
   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNboun(nboun)

   !FOR BVMASS
   call a%Memor%alloc(npoin,a%bvmass,'BVMASS','nsc_BoundaryVmass')
   a%bvmass=0.0_rp

   !We force a closed integration rule for computing vmass
   !Loop over elements and compute vmass
   call a%Mesh%ElementAlloc(e,a%Memor,'ForceClosedRule','nsc_BoundaryVmass')

   ! Loop over boundaries
   boundaries: do iboun=1,nboun
      !Load Element
      call a%Mesh%BoundaryLoad(iboun,e)

      call e%elmdel

      !Gauss Point Loop
      gauss_points : do igaub=1,e%pgaub
        e%igaub = igaub
        !Calculate exterior Normal
        call e%bounor

        auxshape = 0.0_rp
        dsurf=e%weigb(e%igaub)*e%eucta
        auxshape(e%lboel(1:e%pnodb)) = dsurf*e%shapb(1:e%pnodb,e%igaub)
        call a%Mesh%AssemblyToArray(e,1_ip,auxshape,a%bvmass)
      enddo   gauss_points

   end do boundaries

   call a%Mesh%ElementDealloc(e,a%Memor,'ForceClosedRule','nsc_BoundaryVmass')

   !Communicate between subdomains
   call a%Mesh%ArrayCommunicator%GhostCommunicate(1_ip,a%bvmass)

   !Store the inverse of vmass
   do ipoin=1,npoin
      if (a%bvmass(ipoin) /= 0.0_rp) then
         a%bvmass(ipoin)=1.0_rp/a%bvmass(ipoin)
         if (a%bvmass(ipoin) < 0.0_rp) call runend('Negative bvmass value')
      endif
   end do

end subroutine nsc_CalculateBoundaryVmass

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
c******************************************************************************
c
c      Rebasing
c
c******************************************************************************
c
       SUBROUTINE SJERB(NAMAX,NSMAX,NIW,NRW,ATOM,CHEM,IW,RW)
c
c      Rebases the system following a change in the active phases.
c      Checks for permissibility of the demise of phase(s).
c------------------------------------------------------------------------------
c      At call:
c
c          KPCX contains the proposed active phase list
c
c      On return:
c
c          If successful, KERR = 0 and the base system is reset.
c
c          KERR = 1 if the phase change is unfeasible.
c------------------------------------------------------------------------------
c      Sets up a provisional system excluding the species in the absent
c      phase(s). Uses SJISMP to find the base set for this system.  If the
c      populations are impossible in this subset of species, or if the
c      number of bases are changed, the rebasing is not allowed.
c------------------------------------------------------------------------------
c      Nomenclature:
c
c      Variables in the argument list:
c @        NAMAX       maximum number of atoms
c @        NSMAX       maximum number of species
c @        NIW         dimension of work array IW
c @        NRW         dimension of work array IR
c @        ATOM(I)     CHARACTER*2 name of Ith atom
c @        CHEM(J)     CHARACTER*8 name of Jth species
c          IW(I)       integer work array
c          RW(I)       real work array
c
c      Variables in the integer work array IW:
c @#       IB(K) = I   if the Ith atom is the Kth independent system atom
c -        IBO(K) = I  saved IB
c @#       JB(K) = J   if the Jth system species is the Kth base species
c -        JBO(K)      saved JB
c @#       KB(J)       basis control
c -        KBO(J)      saved KB
c #        KERR        error flag
c @        KPCX(K)     trial phase control; 0 if phase empty, 1 if populated
c @        N(I,J)      number of Ith atoms in Jth species
c @        NA          number of atom types
c @        NB          number of base species
c @        NP          number of phases
c @        NS          number of species
c @        NSP(M)      number of species in Mth phase
c
c      Variables used only internally:
c -        NBO         NB for current system
c------------------------------------------------------------------------------
       IMPLICIT    REAL*8  (A-H,O-Z)
       CHARACTER*2 ATOM
       CHARACTER*8 CHEM
c------------------------------------------------------------------------------
       DIMENSION   ATOM(NAMAX),CHEM(NSMAX),IW(NIW),RW(NRW)
c------------------------------------------------------------------------------

c    pointers
       COMMON /SJEPTR/
     ;   IoKERR,IoKMON,IoKTRE,IoKUMO,IoNA,IoNB,IoNP,IoNS,IoIB,IoIBO,
     ;   IoJB,IoJBAL,IoJBA,IoJBB,IoJBO,IoJBX,IoJS2,IoKB,IoKB2,IoKBA,
     ;   IoKBB,IoKBO,IoKPC,IoKPCX,IoLB2,IoMPA,IoMPJ,IoN,LoN,IoNSP,
     ;   IoFRND,IoHUGE,IoR1,IoR2,IoR3,IoA,LoA,IoB,LoB,IoBBAL,
     ;   IoCM,LoCM,IoD,LoD,IoDC,LoDC,IoDPML,IoDLAM,IoDLAY,IoE,
     ;   LoE,IoEEQN,IoELAM,IoELMA,IoELMB,IoF,IoG,IoHA,IoHC,IoPA,
     ;   IoPC,IoPMOL,IoQ,LoQ,IoQC,LoQC,IoRC,LoRC,IoRL,IoRP,
     ;   IoSMOA,IoSMOB,IoSMOO,IoSMOL,IoSMUL,IoW,IoX,IoXO,IoY,IoZ
c------------------------------------------------------------------------------
c    get parameters
       NP = IW(IoNP)
       NS = IW(IoNS)
c
c    save the present base parameters and load trials
       NBO = IW(IoNB)
       DO 9 K=1,NBO
           IW(IoIBO+K) = IW(IoIB+K)
           IW(IoJBO+K) = IW(IoJB+K)
9          CONTINUE
       DO 19 J=1,NS
           IW(IoKBO+J) = IW(IoKB+J)
19         CONTINUE
c
c    set basis control for rebasing
       J2 = 0
       DO 39 M=1,NP
           J1 =J2 + 1
           J2 = J2 + IW(IoNSP+M)
           DO 29 J=J1,J2
               IF (IW(IoKPCX+M).NE.0)  THEN
c                   phase will remain allow all non-excluded species
                       IF (IW(IoKB+J).NE.4)  IW(IoKB+J) = 0
                   ELSE
c                    phase will be absent
                       IW(IoKB+J) = 4
                   ENDIF
29             CONTINUE
39         CONTINUE
c
c    obtain the provisional base (include phase considerations)
       K = 1
       CALL SJISMP(NAMAX,NSMAX,NIW,NRW,ATOM,CHEM,IW,RW,K)
       IF (IW(IoKERR).NE.0)  RETURN
       NB = IW(IoNB)
c
c    check for illegal rebasing
       IF ((IW(IoKERR).NE.0).OR.(IW(IoNB).NE.NBO))  THEN
c        rebasing failed; restore old bases
           IW(IoNB) = NBO
           DO 43 K=1,NB
               IW(IoIB+K) = IW(IoIBO+K)
               IW(IoJB+K) = IW(IoJBO+K)
43             CONTINUE
           DO 49 J=1,NS
               IW(IoKB+J) = IW(IoKBO+J)
49             CONTINUE
           IW(IoKERR) = 1
           RETURN
           ENDIF
c
c    set the base controls
50     DO 59 J=1,NS
c        remove false exclusion flags
           IF ((IW(IoKB+J).EQ.4).AND.(IW(IoKBO+J).NE.4)) IW(IoKB+J) = 0
59         CONTINUE
c
c    calculate conditioning matrix and conditioned populations
       L = 0
       CALL SJICPC(NIW,NRW,IW,RW,L)
       RETURN
       END
c******************************************************************************

max_stars_repo_name: kjfergu/esp

max_issues_repo_name: fergu/esp

max_forks_repo_name: fergu/esp


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE USRBC1 ( VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                    XBD, YBD, NNODES, NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*), V2BD(
     1*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRBC2 ( VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                    XBD, YBD, NNODES, NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*),
     *          V2BD(*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRBC3 ( VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                    XBD, YBD, NNODES, NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*),
     *          V2BD(*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRBC4 (VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                   XBD, YBD, NNODES, NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*),
     *          V2BD(*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRBC5 ( VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                    XBD, YBD, NNODES, 1NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*),
     *          V2BD(*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRBC6 ( VALUE, UBD, VBD, PBD, TBD, V1BD, V2BD,
     *                    XBD, YBD, NNODES, NELEM, TIME, KCYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE A TIME-DEPENDENT BOUNDARY CONDITION
C
C     ******************************************************************
C
      DIMENSION VALUE(*), UBD(*), VBD(*), PBD(*), TBD(*), V1BD(*),
     *          V2BD(*), XBD(*), YBD(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRCON ( AK, T, S, VAR1, VAR2, NNODES, MAT)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE THERMAL CONDUCTIVITY
C
C     ******************************************************************
C
      DIMENSION AK(*), T(*), S(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE  USRCP (RHOCP, T, VAR1, VAR2, NNODES, MAT, RHO)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE SPECIFIC HEAT
C
C     ******************************************************************
C
      DIMENSION RHOCP(*), T(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRDIF ( AD, T, VAR1, VAR2, NNODES, MAT, IEQN)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE DIFFUSION COEFFICIENT (AUXILIARY EQN)
C
C     ******************************************************************
C
      DIMENSION AD(*), T(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USREX ( ZETA, T, VAR1, VAR2, NNODES, MAT, IEQN)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE EXPANSION COEFFICIENT (AUXILIARY EQN)
C
C     ******************************************************************
C
      DIMENSION ZETA(*), T(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USREXT ( BETA, T, VAR1, VAR2, NNODES, MAT)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE THERMAL EXPANSION COEFFICIENT
C
C     ******************************************************************
C
      DIMENSION BETA(*), T(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRGRV ( GX, GY, X, Y, NNODES, NELEM)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE VARIABLY ORIENTED GRAVITY VECTOR
C
C     ******************************************************************
C
      DIMENSION X(*), Y(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRHTC ( HT, TSURF, TREF, XSURF, YSURF, TIME, NELEM,
     *                    IVALUE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE VARIABLE HEAT TRANSFER COEFFICIENT
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRMSH ( X, Y, MAXI, MAXJ)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO GENERATE MESH POINT LOCATIONS
C
C     ******************************************************************
C
      DIMENSION X(MAXI,*), Y(MAXI,*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRVHS ( QVALUE, T, VAR1, VAR2, X, Y, NNODES, MAT,
     *                    NELEM, TIME, ICYCLE)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE VOLUMETRIC HEAT SOURCE
C
C     ******************************************************************
C
      DIMENSION QVALUE(*), T(*), X(*), Y(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END
      SUBROUTINE USRVIS ( AMU, T, S, VAR1, VAR2, NNODES, MAT, ICYCLE)
      DIMENSION AMU(*), T(*), S(*), VAR1(*), VAR2(*)
C
C     ***************************************************************
C
C     USER SUBROUTINE TO EVALUATE DYNAMIC VISCOSITY
C
C     ***************************************************************
C
      RETURN
      END
      SUBROUTINE USRVS ( QVALUE, T, VAR1, VAR2, X, Y, NNODES, MAT,
     *                   NELEM, TIME, ICYCLE, IEQN)
C
C     ******************************************************************
C
C     USER SUBROUTINE TO EVALUATE VOLUMETRIC SOURCE (AUXILIARY EQN)
C
C     ******************************************************************
C
      DIMENSION QVALUE(*), T(*), X(*), Y(*), VAR1(*), VAR2(*)
C
C     ******************************************************************
C
      RETURN
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine extrapolate_ad_h_comp(nface,ielfa,xrlfa,adv,advfa,
     &                  hel,hfa)
!
!     inter/extrapolation of adv at the center of the elements
!     to the center of the faces
!
!     inter/extrapolation of h at the center of the elements
!     to the center of the faces; division through advfa to obtain
!     the face velocity
!
      implicit none
!
      integer nface,ielfa(4,*),ipo1,iel2,ipo3,i,j
!
      real*8 xrlfa(3,*),xl1,advfa(*),adv(*),hel(3,*),hfa(3,*)
!
      do i=1,nface
         ipo1=ielfa(1,i)
         xl1=xrlfa(1,i)
         iel2=ielfa(2,i)
         if(iel2.gt.0) then
!
!           internal face
!
            advfa(i)=1.d0/(xl1/adv(ipo1)+xrlfa(2,i)/adv(iel2))
            do j=1,3
               hfa(j,i)=(xl1*hel(j,ipo1)/adv(ipo1)
     &                   +xrlfa(2,i)*hel(j,iel2)/adv(iel2))
            enddo
         elseif(ielfa(3,i).gt.0) then
!
!           external face; linear extrapolation
!
            ipo3=ielfa(3,i)
            advfa(i)=1.d0/(xl1/adv(ipo1)+xrlfa(3,i)/adv(ipo3))
            do j=1,3
               hfa(j,i)=(xl1*hel(j,ipo1)/adv(ipo1)+
     &                   xrlfa(3,i)*hel(j,ipo3)/adv(ipo3))
            enddo
         else
!
!           external face: constant extrapolation (only one adjacent
!           element layer)
!
            advfa(i)=adv(ipo1)
            do j=1,3
               hfa(j,i)=hel(j,ipo1)/adv(ipo1)
            enddo
         endif
c         write(*,*) 'extrapolate_ad ',ielfa(1,i),ielfa(4,i),advfa(i)
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
subroutine NPZD_N2  ! Add N2 fixation, use P as the limiting element
use bio_MOD
implicit none
integer :: k
real    :: NO3,PO4,PHY, ZOO, DETn,DETp,DIA
real    :: par_, muNet,muDIA
real    :: Kp, gmax, mz,KPO4
real    :: QN, Qp, bI0, KFe_, N2P_phy, theta_DIA
real    :: pp_NDn,pp_PDp,Res_DIA,Mort_DIA
real, parameter :: N2P_dia = 45. ! N:P ratio of diazotrophs
real, parameter :: N2C_dia = 0.2  ! N:C ratio of diazotrophs
real, parameter :: P2C_dia = N2C_dia/N2P_dia ! P:C ratio of diazotrophs
real, parameter :: Q0p     = 3D-3 ! Subsistence P quota of ordinary phyto.
real, parameter :: Lnifdet = 0.05 ! Diazotroph mortality   rate
real, parameter :: RDN     = 0.1  ! Regeneration rate of detritus to DIN
real            :: Lnifdin = 0.17*16

Lnifdin=params(iLnifp)  ! Diazotroph respiration rate
Kp  = params(iKPHY)
gmax= params(igmax)
mz  = params(imz)
!bI0 = 10**params(ibI0B)
bI0 = 0.
KPO4= params(iKN)/16d0
if (do_IRON) KFe = params(iKFe)

DO k = 1, nlev
   tf_p= TEMPBOL(Ep,Temp(k))
   NO3 = Vars(iNO3,   k)
   PO4 = Vars(iPO4,   k)
   PHY = Vars(iPHY(1),k)
   ZOO = Vars(iZOO,   k)
   DETn= Vars(iDET,   k)  ! Detritus N
   DETp= Vars(iDETp,  k)  ! Detritus P

   ! Add diazotrophs:
   DIA = Vars(iDIA,   k)

   ! Check whether in the MLD or not
   if (k .lt. N_MLD) then
      par_ = PAR(k)
   else
      par_ = PARavg
   endif
   Varout(oPAR_,k)=par_
   ! Calculate nondiazotroph growth rate if not considering mixing
   CALL MONOD(Temp(k), PAR(k), NO3,PO4,params(imu0),params(iQ0N),   &
               Q0p, params(iaI0_C),                                 &
               bI0, params(iKN), KPO4, DFe(k), KFe_,         &
               muNet, QN, QP, Varout(oTheta(1),k),                  &
               Varout(oSI(1),k), Varout(oLno3(1),k))

   ! Calculate diazotroph growth rate if not considering mixing
   call nif_growth(PO4,PAR(k),Temp(k),wstr0(1),muDIA,theta_DIA)

   ! Unit of PHY: mM P
   Varout(oPPt,k)=(PHY*muNet/QP + DIA*muDIA/P2C_dia)*12d0

   ! Calculate nondiazo growth rate, theta, and QN based on environmental conditions
   CALL MONOD(Temp(k), par_, NO3, PO4, params(imu0),params(iQ0N),  &
        Q0p,  params(iaI0_C), &
        bI0,  params(iKN), KPO4, DFe(k), KFe_,              &
        muNet,Varout(oQN(1),k),Varout(oQP(1),k),Varout(oTheta(1),k), &
        Varout(oSI(1),k), Varout(oLno3(1),k))

   N2P_phy = Varout(oQN(1),k)/Varout(oQP(1),k)  !N:P of normal phytoplankton

   ! Calculate diazotroph growth rate if not considering mixing
   call nif_growth(PO4,par_,Temp(k),wstr0(1),muDIA,theta_DIA)

   ! Correct the unit of growth rate
   muNet = muNet*dtdays

   Varout(oDIAu,k)=muDIA
   muDIA = muDIA*dtdays

   ! The total amount of phytoplankton grazed by zooplankton (molN;gmax is the maximal specific ingestion rate!)
   tf_z  = TEMPBOL(Ez,Temp(k))  ! Temperature effect of zooplankton
   gbar  = grazing(grazing_formulation,Kp,PHY) !Functional response of zooplankton
   INGES = gmax*tf_z*dtdays * gbar   !Specific ingestion rate of zooplankton
   Zmort = ZOO*ZOO*dtdays*mz*tf_z  !Mortality term for ZOO

   !Zooplankton excretion rate (-> DOM)
   RES = INGES*(1d0-GGE-unass)

   !ZOOPLANKTON EGESTION (-> POM)
   EGES = INGES*unass

  ! For production/destruction matrix:
  pp_NDn = dtdays*RDN*DETn*tf_z
  pp_PDp = dtdays*params(iRDN_P)*DETp*tf_z
  pp_NZ  = ZOO*RES
  pp_DZ  = ZOO*EGES+Zmort
  pp_ZP  = ZOO*INGES

  ! Respiration of diazotrophs (DIA -> DIP):
  Res_DIA  = Lnifdin*tf_p*dtdays*DIA**2

  ! Mortality of diazotrophs (DIA -> POP):
  Mort_DIA = Lnifdet*tf_p*dtdays*DIA

  ! N as the unit:
  Varout(oDET,k)  = (DETn + pp_DZ*N2P_phy) - pp_NDn + Mort_DIA*N2P_dia

  ! P as the unit
  Varout(oDETp,k) = (DETp + pp_DZ)- pp_PDp + Mort_DIA

  Varout(oNO3,k)  = (NO3+pp_NDn+pp_NZ*N2P_phy+Res_DIA*N2P_dia) &
                  - PHY*muNet*N2P_phy   ! N as the unit

  Varout(oPO4,k)  = (PO4+pp_PDp+pp_NZ        +Res_DIA) &
                  - PHY*muNet-DIA*muDIA ! P as the unit

  Varout(oPHY(1),k)    = PHY*(1d0 + muNet)-pp_ZP
  Varout(oZOO,k)       = (ZOO+pp_ZP)-pp_DZ-pp_NZ
  Varout(oDIA,k)       = DIA*(1d0 + muDIA)-Mort_DIA-Res_DIA

  Varout(omuNet(1), k) = muNet/dtdays
  Varout(oGraz(1), k)  = pp_ZP/PHY/dtdays
  Varout(oZ2N,k)   = pp_NZ/dtdays
  Varout(oD2N,k)   = pp_ND/dtdays
  Varout(oCHLt,k)  = PHY/Varout(oQP(1),k)*Varout(otheta(1),k) &
                   + DIA/P2C_dia*theta_DIA
  Varout(oCHL(1),k)= Varout(oCHLt,k)

  ! Total particulate organic nitrogen (PON)
  Varout(oPON,k)=(Varout(oZOO,k)+Varout(oPHY(1),k))*N2P_phy &
            +     Varout(oDET,k)+Varout(oDIA,k)*N2P_dia

  ! Total POP:
  Varout(oPOP,k)=Varout(oZOO,k)+Varout(oPHY(1),k) + Varout(oDETp,k) + Varout(oDIA,k)
Enddo
return
End subroutine NPZD_N2

max_stars_repo_name: BingzhangChen/NPZDFeCONT

max_issues_repo_name: BingzhangChen/NPZDFeCONT

max_forks_repo_name: BingzhangChen/NPZDFeCONT


text:
      SUBROUTINE SELBUF(NBUFFR)
C
C          ------------------------------------------------
C          ROUTINE NO. (  12)   VERSION (A7.4)    11:FEB:85
C          ------------------------------------------------
C
C          THIS SELECTS BUFFER NO. 'NBUFFR'.
C
C
      REAL    RDATA(1)
      INTEGER IDATA(1)
      LOGICAL ERRON
C
      COMMON /T0BUFN/ KBUFR0
      COMMON /T0TRAC/ IPRINT
      COMMON /T0TRAI/ ITRAC1,ITRAC2,ITRAC3,ITRAC4
      COMMON /T3ERRS/ ERRON,NUMERR
C
      DATA RDATA /0.0/
C
C
      CALL G3INIT(2)
C
      ITRAC1= NBUFFR
      IF (IPRINT.EQ.1) CALL G0MESG(114,5)
      IF (NBUFFR.LE.0.OR.NBUFFR.GT.16) GO TO 901
C
      KBUFR0= NBUFFR
      IDATA(1)= NBUFFR
      CALL G3LINK(4,1,-1,IDATA,RDATA)
      RETURN
C
  901 NUMERR= 6
      IF (ERRON) CALL G0ERMS
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
program main
  use gobject_aux
  use iso_c_binding

  implicit none

  block
  end block
end program main

max_stars_repo_name: dram/fortran-gobject

max_issues_repo_name: dram/fortran-gobject

max_forks_repo_name: dram/fortran-gobject