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 precfd(ne,ipkon,kon,lakon,ipnei,neifa,neiel,ipoface,
     &  nodface,ielfa,nflnei,nface,ifaext,nfaext,
     &  isolidsurf,nsolidsurf,set,nset,istartset,iendset,ialset,
     &  vel,vold,mi,neij,nef,nactdoh,ipkonf,lakonf,ielmatf,ielmat,
     &  ielorienf,ielorien,norien,cs,mcs,tieset,x,y,z,xo,yo,zo,nx,
     &  ny,nz,co,ifatie)
!
!     gathering topological information for computational fluid
!     dynamics calculations
!
      implicit none
!
      character*8 lakon(*),lakonf(*)
      character*81 set(*),noset,tieset(3,*),slavset,mastset
!
      integer ne,ipkon(*),ipnei(*),ipoface(*),nodface(5,*),neifa(*),
     &  ielfa(4,*),nflnei,nface,i,j,k,index,indexe,neiel(*),ithree,
     &  nfaext,ifaext(*),isolidsurf(*),nsolidsurf,indexf,nneigh,
     &  nset,istartset(*),iendset(*),ialset(*),iaux,kflag,ifour,iel2,
     &  ifaceq(8,6),ifacet(7,4),ifacew(8,5),kon(*),nodes(4),iel1,j2,
     &  indexold,ifree,ifreenew,ifreenei,mi(*),neij(*),ifreenei2,nef,
     &  nactdoh(*),ipkonf(*),ielmatf(mi(3),*),ielmat(mi(3),*),nf(5),
     &  nope,ielorien(mi(3),*),ielorienf(mi(3),*),norien,jopposite8(6),
     &  jopposite6(5),itie,nx(*),ny(*),nz(*),noden(1),nelemm,nelems,
     &  n,mcs,l,jfacem,jfaces,islav,imast,ifaces,ifacem,ifatie(*),
     &  nodeinface,nodeoutface,nopes
!
      real*8 vel(nef,0:7),vold(0:mi(2),*),coords(3),cs(17,*),x(*),
     &  y(*),z(*),xo(*),yo(*),zo(*),co(3,*),a(3),b(3),xn(3),p(3),
     &  q(3),c(3,3),dot,dc,ds,dd,theta,pi
!
!     nodes belonging to the element faces
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,11,
     &             1,2,4,5,9,8,12,
     &             2,3,4,6,10,9,13,
     &             1,4,3,8,10,7,14/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
      data jopposite6 /2,1,0,0,0/
      data jopposite8 /2,1,5,6,3,4/
      data nf /3,3,4,4,4/
!
!     transfer element information from ipkon, lakon, ielmat
!     and ielorien to ipkonf, lakonf, ielmatf and ielorienf.
!     The latter fields were created
!     for CFD applications in which the elements are numerated in
!     ascending order (without gaps)
!
      do i=1,ne
c         write(*,*) 'precfd ',i,nactdoh(i)
         if(nactdoh(i).ne.0) then
            ipkonf(nactdoh(i))=ipkon(i)
            lakonf(nactdoh(i))=lakon(i)
            do j=1,mi(3)
               ielmatf(j,nactdoh(i))=ielmat(j,i)
               if(norien.gt.0) ielorienf(j,nactdoh(i))=ielorien(j,i)
            enddo
         endif
      enddo
!
      kflag=1
      ithree=3
      ifour=4
!
!     determining the external element faces of the fluid mesh
!     the faces are catalogued by the three lowes nodes numbers
!     in ascending order. ipoface(i) points to a face for which
!     node i is the lowest node and nodface(1,ipoface(i)) and
!     nodface(2,ipoface(i)) are the next lower ones.
!     nodface(3,ipoface(i)) contains the element number,
!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
!     is a pointer to the next surface for which node i is the
!     lowest node; if there are no more such surfaces the pointer
!     has the value zero
!     An external element face is one which belongs to one element
!     only
!
      ifree=1
      ifreenei=0
!
      do i=1,6*nef
         nodface(5,i)=i+1
      enddo
      do i=1,nef
         indexe=ipkonf(i)
         if(lakonf(i)(4:4).eq.'8') then
!
!           hex element
!
            ipnei(i)=ifreenei
            do j=1,6
               do k=1,4
                  nodes(k)=kon(indexe+ifaceq(k,j))
               enddo
               call isortii(nodes,iaux,ifour,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         else if(lakonf(i)(4:4).eq.'6') then
!
!           wedge element
!
            ipnei(i)=ifreenei
            do j=1,5
               do k=1,nf(j)
                  nodes(k)=kon(indexe+ifacew(k,j))
               enddo
               call isortii(nodes,iaux,nf(j),kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         else
!
!           tet element
!
            ipnei(i)=ifreenei
            do j=1,4
               do k=1,3
                  nodes(k)=kon(indexe+ifacet(k,j))
               enddo
               call isortii(nodes,iaux,ithree,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         endif
      enddo
!
      nflnei=ifreenei
      ipnei(nef+1)=ifreenei
      nface=ifree-1
!
!     check for cyclic symmetry (translational or rotational)
!
      do i=1,mcs
         itie=int(cs(17,i))
         if((tieset(1,itie)(81:81).ne.'P').and.
     &      (tieset(1,itie)(81:81).ne.'Z')) cycle
!
         slavset=tieset(2,itie)
         mastset=tieset(3,itie)
!
         do j=1,nset
            if(set(j).eq.slavset) exit
         enddo
         islav=j
!
         do j=1,nset
            if(set(j).eq.mastset) exit
         enddo
         imast=j
!
!        catalogueing the center of gravity of the master faces
!
         n=0
         do j=istartset(imast),iendset(imast)
            ifacem=ialset(j)
!
            nelemm=int(ifacem/10)
            jfacem=ifacem-nelemm*10
!
!           taking the renumbering of the elements into account
!
            nelemm=nactdoh(nelemm)
!
            indexe=ipkonf(nelemm)
!
!           determination of the face center of the master face
!
            if(lakonf(nelemm)(4:4).eq.'8') then
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,4
                  nodes(k)=kon(indexe+ifaceq(k,jfacem))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/4.d0
               enddo
            elseif(lakonf(nelemm)(4:4).eq.'6') then
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,nf(jfacem)
                  nodes(k)=kon(indexe+ifacew(k,jfacem))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/nf(jfacem)
               enddo
            else
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,3
                  nodes(k)=kon(indexe+ifaceq(k,jfacem))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/3.d0
               enddo
            endif
!
            if(j.eq.istartset(imast)) then
               if(tieset(1,itie)(81:81).eq.'Z') then
!
!                 check the correct oprientation of the axis vector
!
                  if(lakonf(nelemm)(4:4).eq.'8') then
                     nope=8
                     nopes=4
                  elseif(lakonf(nelemm)(4:4).eq.'6') then
                     nope=6
                     nopes=nf(jfacem)
                  else
                     nope=4
                     nopes=3
                  endif
!
!                 two nodes on the axis
!
                  do k=1,3
                     a(k)=cs(5+k,i)
                     b(k)=cs(8+k,i)
                  enddo
!
!                 looking for a master node not on the axis
!
                  do k=1,nopes
                     if((co(1,nodes(k))-a(1))**2+
     &                  (co(2,nodes(k))-a(2))**2+
     &                  (co(3,nodes(k))-a(3))**2.gt.1.d-20) exit
                  enddo
                  nodeinface=nodes(k)
!
!                 looking for a node belonging to the same element
!                 but not in the master face
!
                  loop:do k=1,nope
                     nodeoutface=kon(indexe+k)
                     do l=1,nopes
                        if(nodes(l).eq.nodeoutface) cycle loop
                     enddo
                     exit
                  enddo loop
!
!                 normalized vector along the axis
!
                  dd=dsqrt((b(1)-a(1))**2+
     &                     (b(2)-a(2))**2+
     &                     (b(3)-a(3))**2)
                  do k=1,3
                     xn(k)=(b(k)-a(k))/dd
                  enddo
!
!                 vector connecting a on the axis with
!                 nodeinface and nodeoutface
!
                  do k=1,3
                     p(k)=co(k,nodeinface)-a(k)
                     q(k)=co(k,nodeoutface)-a(k)
                  enddo
!
!                 n.(pxq) must be negative for the axis vector
!                 to have the right orientation (= you turn from the
!                 slave surface to the master surface through the
!                 body in clockwise direction while looking in
!                 the direction of xn)
!
                  dot=xn(1)*(p(2)*q(3)-q(2)*p(3))
     &               +xn(2)*(p(3)*q(1)-q(3)*p(1))
     &               +xn(3)*(p(1)*q(2)-q(1)*p(2))
!
!                 revert the direction of the axis if not correct
!
                  if(dot.gt.0.d0) then
                     do k=1,3
                        b(k)=2.d0*a(k)-b(k)
                        cs(8+k,i)=b(k)
                        xn(k)=-xn(k)
                     enddo
                  endif
!
!                 angle from the master to the slave surface
!
                  pi=4.d0*datan(1.d0)
                  theta=-2.d0*pi/cs(1,i)
!
!                 rotation matrix rotating a vector in the master
!                 surface into a vector in the slave surface
!
                  dc=dcos(theta)
                  ds=dsin(theta)
!
!                 C-matrix from Guido Dhondt, The Finite Element
!                 Method for Three-Dimensional Thermomechanical
!                 Applications p 158
!
                  c(1,1)=dc+(1.d0-dc)*xn(1)*xn(1)
                  c(1,2)=   (1.d0-dc)*xn(1)*xn(2)-ds*xn(3)
                  c(1,3)=   (1.d0-dc)*xn(1)*xn(3)+ds*xn(2)
                  c(2,1)=   (1.d0-dc)*xn(2)*xn(1)+ds*xn(3)
                  c(2,2)=dc+(1.d0-dc)*xn(2)*xn(2)
                  c(2,3)=   (1.d0-dc)*xn(2)*xn(3)-ds*xn(1)
                  c(3,1)=   (1.d0-dc)*xn(3)*xn(1)-ds*xn(2)
                  c(3,2)=   (1.d0-dc)*xn(3)*xn(2)+ds*xn(1)
                  c(3,3)=dc+(1.d0-dc)*xn(3)*xn(3)
               endif
            endif
!
!           storing the face centers (translated from the master
!           face into the slave face) in the vectors x,y and z
!
            n=n+1
            if(tieset(1,itie)(81:81).eq.'P') then
               x(n)=coords(1)-cs(6,i)
               y(n)=coords(2)-cs(7,i)
               z(n)=coords(3)-cs(8,i)
            else
!
!              vector orthogonal to axis (in the master surface)
!
               do k=1,3
                  p(k)=coords(k)-a(k)
               enddo
               dd=p(1)*xn(1)+p(2)*xn(2)+p(3)*xn(3)
               do k=1,3
                  p(k)=p(k)-dd*xn(k)
               enddo
!
!              vector rotated into the slave surface
!
               do k=1,3
                  q(k)=c(k,1)*p(1)+c(k,2)*p(2)+c(k,3)*p(3)
               enddo
!
!              rotated node in the slave surface
!
               x(n)=coords(1)+q(1)-p(1)
               y(n)=coords(2)+q(2)-p(2)
               z(n)=coords(3)+q(3)-p(3)
            endif
         enddo
!
!        preparing the fields for near3d
!
         do j=1,n
            nx(j)=j
            ny(j)=j
            nz(j)=j
            xo(j)=x(j)
            yo(j)=y(j)
            zo(j)=z(j)
         enddo
         kflag=2
         call dsort(x,nx,n,kflag)
         call dsort(y,ny,n,kflag)
         call dsort(z,nz,n,kflag)
!
!        loop over all slave faces
!
         do j=istartset(islav),iendset(islav)
            ifaces=ialset(j)
!
            nelems=int(ifaces/10)
            jfaces=ifaces-nelems*10
!
!           taking the renumbering of the elements into account
!
            nelems=nactdoh(nelems)
!
            indexe=ipkonf(nelems)
!
!           determination of the face center of the slave face
!
            if(lakonf(nelems)(4:4).eq.'8') then
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,4
                  nodes(k)=kon(indexe+ifaceq(k,jfaces))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/4.d0
               enddo
            elseif(lakonf(nelems)(4:4).eq.'6') then
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,nf(jfaces)
                  nodes(k)=kon(indexe+ifacew(k,jfaces))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/nf(jfaces)
               enddo
            else
               do l=1,3
                  coords(l)=0.d0
               enddo
               do k=1,3
                  nodes(k)=kon(indexe+ifaceq(k,jfaces))
                  do l=1,3
                     coords(l)=coords(l)+co(l,nodes(k))
                  enddo
               enddo
               do l=1,3
                  coords(l)=coords(l)/3.d0
               enddo
            endif
!
            nneigh=1
            call near3d(xo,yo,zo,x,y,z,nx,ny,nz,coords(1),
     &           coords(2),coords(3),n,noden,nneigh)
!
            ifacem=ialset(istartset(imast)+noden(1)-1)
!
            nelemm=int(ifacem/10)
            jfacem=ifacem-nelemm*10
!
!           taking the renumbering of the elements into account
!
            nelemm=nactdoh(nelemm)
!
            ielfa(2,neifa(ipnei(nelems)+jfaces))=nelemm
            ielfa(2,neifa(ipnei(nelemm)+jfacem))=nelems
!
            neiel(ipnei(nelems)+jfaces)=nelemm
            neiel(ipnei(nelemm)+jfacem)=nelems
!
            neij(ipnei(nelems)+jfaces)=jfacem
            neij(ipnei(nelemm)+jfacem)=jfaces
!
!           field indicating whether a face (slave: +, master: -)
!           is part of a cyclic symmetry condition i
!
            ifatie(neifa(ipnei(nelems)+jfaces))=i
            ifatie(neifa(ipnei(nelemm)+jfacem))=-i
         enddo
      enddo
!
!     catalogueing external faces
!
      nfaext=0
!
      do i=1,nface
         if(ielfa(2,i).ne.0) cycle
         nfaext=nfaext+1
         ifaext(nfaext)=i
         iel1=ielfa(1,i)
         if(lakonf(iel1)(4:4).eq.'8') then
            j=ielfa(4,i)
            j=jopposite8(j)
         elseif(lakonf(iel1)(4:4).eq.'6') then
            j=ielfa(4,i)
            j=jopposite6(j)
            if(j.eq.0) cycle
         else
            cycle
         endif
         indexf=ipnei(iel1)
         ielfa(3,i)=neiel(indexf+j)
      enddo
!
c      do i=1,nef
c         write(*,*)'precfd neiel',i,ipnei(i),(neiel(ipnei(i)+j),j=1,6)
c         write(*,*)'precfd neij',i,ipnei(i),(neij(ipnei(i)+j),j=1,6)
c         write(*,*)'precfd neifa',i,ipnei(i),(neifa(ipnei(i)+j),j=1,6)
c      enddo
c      do k=1,nface
c         write(*,*) 'precfd ielfa',k,(ielfa(j,k),j=1,4),ifatie(k)
c      enddo
c      do k=1,nfaext
c         write(*,*) 'precfd ifaext',k,ifaext(k)
c      enddo
c      write(*,*)
!
!     faces belonging to solid surfaces
!
      nsolidsurf=0
      noset(1:13)='SOLIDSURFACET'
      do i=1,nset
         if(set(i)(1:13).eq.noset(1:13)) exit
      enddo
      if(i.gt.nset) then
         write(*,*) '*WARNING in precfd: facial surface SOLID SURFACE '
         write(*,*) '         has not been defined.'
         write(*,*)
      else
c         nsolidsurf=0
         do j=istartset(i),iendset(i)
            nsolidsurf=nsolidsurf+1
            isolidsurf(nsolidsurf)=ialset(j)
         enddo
         call isortii(isolidsurf,iaux,nsolidsurf,kflag)
      endif
c      do i=1,nsolidsurf
c         write(*,*) 'precfd solidsurf ',i,isolidsurf(i)
c      enddo
!
!     initial conditions: element values
!     vel was initialized at allocation
!
      do i=1,nef
         indexe=ipkonf(i)
         if(lakonf(i)(4:4).eq.'8') then
            nope=8
         elseif(lakonf(i)(4:4).eq.'6') then
            nope=6
         else
            nope=4
         endif
         do j=0,4
            do k=1,nope
               vel(i,j)=vel(i,j)+vold(j,kon(indexe+k))
            enddo
            vel(i,j)=vel(i,j)/nope
         enddo
      enddo
c      do i=1,nef
c         write(*,*) 'precfd vel ',i,(vel(j,i),j=1,3)
c      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: #redirect Donald G. Brewer

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module alm_blm_clm
contains
  !SUBROUTINE cmp_alm_blm_clm(alm,blm,clm,aalm,bblm,cclm,Aweight,latom,lfirst,jatom,nemin,nemax,n0,nnlo,is,DM_nemin,DM_nemax,DM_nemaxx,nbands,nbandsx,Qforce_j,lmax2lmax2_,nbands_dft_,dim_,lomaxlomax_,nloat_,time_bl,time_bl_w,time_lo,time_lo_w)
  SUBROUTINE cmp_alm_blm_clm(alm,blm,clm,aalm,bblm,cclm,Aweight,latom,lfirst,jatom,nemin,nemax,n0,nnlo,is,DM_nemin,DM_nemax,DM_nemaxx,nbands,nbandsx,Qforce_j,time_bl,time_bl_w,time_lo,time_lo_w)
    USE defs,  ONLY: pi, imag
    USE param, ONLY: lmax2, lomax, nume, nloat, iblock, nmat
    USE xa,    ONLY: bkrloc, fj, dfj
    USE xa3,   ONLY: bk3, k3, As
    USE structure, ONLY: rmt, vol, rotij, br1, pos, tauij, rotloc
    USE atspdt,ONLY: p, dp, pe, dpe
    USE lo,    ONLY: nlo, lapw, ilo
    USE dmfts, ONLY: iso, shft
    USE mod_lo, ONLY: lomain
    IMPLICIT NONE
    COMPLEX*16, intent(out) :: alm((lmax2+1)*(lmax2+1),nume),blm((lmax2+1)*(lmax2+1),nume),clm((lomax+1)*(lomax+1),nume,nloat)
    COMPLEX*16, intent(in)  :: Aweight(nbands,nbands)
    INTEGER, intent(in)     :: latom,lfirst,jatom,nemin,nemax,n0,nnlo,is,DM_nemin,DM_nemax,DM_nemaxx,nbands,nbandsx
    LOGICAL, intent(in)     :: Qforce_j
    REAL*8, intent(inout)   :: time_bl, time_bl_w, time_lo, time_lo_w
    ! These data can not be passed as arguments to the function, because they might not be allocated, hence is should be packed inside the module.
    COMPLEX*16, ALLOCATABLE :: aalm(:,:,:), bblm(:,:,:), cclm(:,:,:,:)
    !INTEGER, intent(in)     :: lmax2lmax2_,nbands_dft_,dim_,lomaxlomax_,nloat_
    !COMPLEX*16, intent(out) :: aalm(lmax2lmax2_,nbands_dft_,dim_), bblm(lmax2lmax2_,nbands_dft_,dim_), cclm(lomaxlomax_,nbands_dft_,nloat_,dim_)
    ! locals
    REAL*8  :: fac, twopi, arg1, argt2, arg2, rmt2
    COMPLEX*16 :: phshel, cfac
    REAL*8  :: bkrot2(3)
    REAL*8  :: time1, time1_w, time2, time2_w
    INTEGER :: ii,ibb,i,i3, index, lmx, l, m, maxx, lda, ldb, ldc, jlo, num, lmax2lmax2, lomaxlomax, i_h_k, ik, ilo_max, ind, nbands_dft, ip
    REAL*8 :: h_al(iblock),h_bl(iblock)
    COMPLEX*16:: h_yl
    COMPLEX*16 :: yl((LMAX2+1)*(LMAX2+1)), h_alyl((LMAX2+1)*(LMAX2+1),iblock), h_blyl((LMAX2+1)*(LMAX2+1),iblock)
    COMPLEX*16, allocatable :: xlm_tmp(:,:)
    REAL*8, allocatable     :: h_k(:,:)
    complex*16, allocatable :: h_ablyl_hk(:,:)
    complex*16, allocatable :: ylm_tmp(:,:)
    LOGICAL :: nonzero_shft
    REAL*8  :: rotloc_x_BR1(3,3), rotloc_x_BR1_x_rotij(3,3)
    INTEGER :: j
    complex*16 :: cone, czero
    cone  = cmplx(1.d0, 0.d0, 8)
    czero = cmplx(0.d0, 0.d0, 8)

    CALL cputim(time1)
    CALL walltim(time1_w)

    lmax2lmax2 = (lmax2+1)*(lmax2+1)
    lomaxlomax = (lomax+1)*(lomax+1)
    nbands_dft = nemax-nemin+1
    twopi = pi*2.d0
    fac=4.0d0*pi*rmt(jatom)**2/SQRT(vol)
    lmx=lmax2
    ilo_max=0
    DO l=0,lmax2
       if (ilo(l).gt.ilo_max) ilo_max=ilo(l)
    ENDDO

    alm(:,:)=0.d0
    blm(:,:)=0.d0
    clm(:,:,:)=0.d0

    if (Qforce_j) then
       allocate( h_ablyl_hk(lmax2lmax2,iblock) )
       allocate( h_k(iblock,3) )

       if (size(aalm,3).NE.3 .or. size(cclm,4).NE.3) print*, 'ERROR alm_blm_clm: Third dimension of aalm should be 3 but is', size(aalm,3)
       if (size(aalm,1).NE.lmax2lmax2 .or. size(cclm,1).NE.lomaxlomax) print*, 'ERROR alm_blm_clm: first dimension of aalm should be', lmax2lmax2, 'but is', size(aalm,1)
       if (size(cclm,3).NE.nloat) print*, 'ERROR alm_blm_clm: Second dimension of cclm should be nloat but is', size(cclm,3)
       if (size(aalm,2).NE.size(cclm,2)) print*, 'ERROR alm_blm_clm: nbands_dft should be dimension in aalm and cclm', size(aalm,3)
       if (size(aalm,2).NE.nbands_dft) print*, 'ERROR alm_blm_clm: size(aalm,2)!=nbands_dft', size(aalm,2), nbands_dft

       aalm(1:lmax2lmax2,1:nbands_dft,1:3)=0.d0
       bblm(1:lmax2lmax2,1:nbands_dft,1:3)=0.d0
       cclm(1:lomaxlomax,1:nbands_dft,1:nloat,1:3)=0.d0
    endif

    rotloc_x_BR1(:,:) = matmul(rotloc(:,:,jatom), BR1)
    rotloc_x_BR1_x_rotij = matmul( rotloc_x_BR1, rotij(:,:,latom) )
    nonzero_shft = sum(abs(shft(latom,:))) .GT. 1e-10

    DO ii=1,n0-nnlo,iblock  !! ii=1,iblock,2*iblock,3*ibloc,... n0-nnlo+d
       ibb=MIN(iblock,n0-nnlo-ii+1)  ! size of this block
       if(ibb.gt.0) then
          if (Qforce_j) h_k(:,:) = 0.d0
          !$OMP PARALLEL DO PRIVATE(bkrot2,bkrloc,yl,arg1,arg2,argt2,phshel)&
          !$OMP& SHARED(h_k,h_alyl)&
          !$OMP& SCHEDULE(STATIC)
          DO i=ii,min(ii+iblock-1,n0-nnlo)
             !---------  rotates ylm(k+K) to ylm(k'+K) where k' is in star of irreducible k. ------------!
             ! bk3(:,i) !-----  reciprocal vector and irreducible vector: G=K+k ----!
             ! BKROT2 = R_n.R_a.(k+K), transformation from the first atom to an equivalent atom
             bkrot2 = matmul(rotij(:,:,latom), bk3(:,i))
             ! ! !---- BR1 transforms integer reciprocal lattice vectors, as given in the VECTORLIST of LAPW1, into cartesian system ----!
             ! ! ! BKROT3 = R_n.R_a.(k+K), but in cartesian coordinate system
             ! ! bkrot3= matmul(br1,bkrot2)
             ! ! !---- BKRLOC = Rotloc.R_g.(k+K),  is rotates Rotloc.R.(k+K) . Rotation Rotloc entered by user
             ! ! ! Here we use rotloc and not crotloc, hence we rather use the transformation
             ! ! ! from the structure file and not the user input. This is because the wave functions are
             ! ! ! compatible with the rotation in the structure file!
             ! ! bkrloc = matmul(rotloc(:,:,jatom),bkrot3)
             bkrloc = matmul( rotloc_x_BR1,bkrot2)
             !---- YLM = Y_{L}(Rotloc.R_g.(k+K))
             CALL YLM (bkrloc,lmax2,yl,lmax2lmax2)
             ! (R_n.R_a.(k+K)) *  R(first) * 2pi
             arg1 = dot_product(bkrot2, pos(:,lfirst))*twopi
             ! ARGT2 = (R_a.(k+K)).tau_n * 2pi
             argt2= dot_product(bk3(:,i), tauij(:,latom))*twopi
             ! ARG2 = (R_a.(k+K)) *  shft * 2pi
             arg2=0.d0
             if (nonzero_shft) arg2 = dot_product(bk3(:,i), shft(latom,:))*twopi
             ! PHSEHL = e^{I*2pi*( (R_a.(k+K))*tau_n + (K+k)*tau(isym) + (R_n.R_a.(k+K)*R(first)))}
             phshel=EXP( imag*(arg1+arg2+argt2) )
             h_alyl(:(lmx+1)**2,i-ii+1) = dconjg(yl(:(lmx+1)**2))*phshel
             IF(Qforce_j) THEN
                h_k(i-ii+1,:) = matmul( rotloc_x_BR1_x_rotij, k3(:,i) )
             ENDIF
          ENDDO
          !$OMP END PARALLEL DO

          rmt2=1.0D0/(rmt(jatom)**2)
          DO l=0,lmx
             ! !$OMP PARALLEL DO PRIVATE(i3) SHARED(h_al,h_bl) SCHEDULE(STATIC)
             DO i=ii,min(ii+iblock-1,n0-nnlo)
                i3=i-ii+1
                if(lapw(l)) then
                   h_al(i3)=dfj(l,i,jatom)*pe(l)-fj(l,i,jatom)*dpe(l)
                   h_bl(i3)=fj(l,i,jatom)*dp(l)-dfj(l,i,jatom)*p(l)
                ELSE
                   h_al(i3) = rmt2*fj(l,i,jatom)/p(l)
                   h_bl(i3) = 0.d0
                ENDIF
             ENDDO
             ! !$OMP END PARALLEL DO
             DO i3=1,ibb ! size of this block
                h_blyl(l**2+1:(l+1)**2,i3)=h_bl(i3)*h_alyl(l**2+1:(l+1)**2,i3)
                h_alyl(l**2+1:(l+1)**2,i3)=h_al(i3)*h_alyl(l**2+1:(l+1)**2,i3)
             ENDDO
          ENDDO
          !
          lda=lmax2lmax2
          ldc=lda
          ldb=nmat
          !---- h_alyl(2*lmax+1,iblock)  contains rotated Apw's, such that chi(r) = (Apw*u(r) + Bpw*udot(r))*Ylm(r)
          !---- h_blyl(2*lmax+1,iblock)  contains rotated Bpw's
          !---- A(:N,iband,is)              contains eigenvectors, also named C(k+G,inad,is)
          !---- alm[:,iband][irf=1,is] += sum_{iK\in block} h_alyl[:,iK][is]*A[iK,iband][is]
          !---- alm[:,iband][irf=2,is] += sum_{iK\in block} h_blyl[:,iK][is]*A[iK,iband][is]
          !----
          !---- The results is:
          !---- alm[lm,iband][1,is] = sum_G Apw(lm,is,K+G) * C(k+G,iband,is)
          !---- alm[lm,iband][2,is] = sum_G Bpw(lm,is,K+G) * C(k+G,iband,is)
          !---- Where C(k+G,iband,is) are eigenvectors, and Apw and Bpw are expansion coefficients defined in Shick et.al., PRB 60, 10763 (1999).
          CALL zgemm('N','N',lmax2lmax2,nemax-nemin+1,ibb,cone,h_alyl,lda,As(ii,nemin,is),ldb,cone,alm(1,nemin),ldc)
          CALL zgemm('N','N',lmax2lmax2,nemax-nemin+1,ibb,cone,h_blyl,lda,As(ii,nemin,is),ldb,cone,blm(1,nemin),ldc)
          IF (Qforce_j) THEN
             do i_h_k=1,3
                !$OMP PARALLEL DO SHARED(h_ablyl_hk,h_alyl,h_k) SCHEDULE(STATIC)
                do i3=1,ibb
                   h_ablyl_hk(:,i3) = h_alyl(:,i3)*h_k(i3,i_h_k)   ! h_ablyl <-  alm(lm,K)*K
                enddo
                !$OMP END PARALLEL DO
                CALL zgemm('N','N',lmax2lmax2,nemax-nemin+1,ibb,cone,h_ablyl_hk,lda,As(ii,nemin,is),ldb,cone,aalm(1,1,i_h_k),ldc)
             enddo
             do i_h_k=1,3
                !$OMP PARALLEL DO SHARED(h_ablyl_hk,h_blyl,h_k) SCHEDULE(STATIC)
                do i3=1,ibb
                   h_ablyl_hk(:,i3) = h_blyl(:,i3)*h_k(i3,i_h_k)   ! h_ablyl <-  blm(lm,K)*K
                enddo
                !$OMP END PARALLEL DO
                CALL zgemm('N','N',lmax2lmax2,nemax-nemin+1,ibb,cone,h_ablyl_hk,lda,As(ii,nemin,is),ldb,cone,bblm(1,1,i_h_k),ldc)
             enddo
          ENDIF
       endif
    ENDDO

    if (Qforce_j) then
       deallocate( h_k )
       deallocate( h_ablyl_hk )
    endif

    CALL cputim(time2)
    CALL walltim(time2_w)

    time_bl=time_bl+time2-time1
    time_bl_w=time_bl_w+time2_w-time1_w

    CALL cputim(time1)
    CALL walltim(time1_w)

    !! output: alm, blm, clm, aalm, bblm, cclm
    if (nlo.ne.0) call lomain(rotloc(:,:,jatom),is,nemin,nemax,lfirst,latom,n0,jatom,alm,blm,clm,Qforce_j,aalm,bblm,cclm,lmax2)

    !$OMP PARALLEL DO PRIVATE(cfac) SHARED(alm,blm) SCHEDULE(STATIC)
    DO l=0,lmax2
       cfac=fac*imag**l
       alm(l**2+1:(l+1)**2,nemin:nemax)=alm(l**2+1:(l+1)**2,nemin:nemax)*cfac
       blm(l**2+1:(l+1)**2,nemin:nemax)=blm(l**2+1:(l+1)**2,nemin:nemax)*cfac
    ENDDO
    !$OMP END PARALLEL DO
    !$OMP PARALLEL DO PRIVATE(cfac) SHARED(clm) SCHEDULE(STATIC)
    DO l=0,lomax
       cfac=fac*imag**l
       clm(l**2+1:(l+1)**2,nemin:nemax,1:ilo(l))=clm(l**2+1:(l+1)**2,nemin:nemax,1:ilo(l))*cfac
    ENDDO
    !$OMP END PARALLEL DO

    IF (Qforce_j) THEN
       !$OMP PARALLEL DO PRIVATE(cfac) SHARED(aalm,bblm) SCHEDULE(STATIC)
       DO l=0,lmax2
          cfac = fac*imag**l
          aalm(l**2+1:(l+1)**2,1:nbands_dft,1:3) = aalm(l**2+1:(l+1)**2,1:nbands_dft,1:3)*cfac
          bblm(l**2+1:(l+1)**2,1:nbands_dft,1:3) = bblm(l**2+1:(l+1)**2,1:nbands_dft,1:3)*cfac
       ENDDO
       !$OMP END PARALLEL DO
       !$OMP PARALLEL DO PRIVATE(cfac) SHARED(cclm) SCHEDULE(STATIC)
       DO l=0,lomax
          cfac = fac*imag**l
          cclm(l**2+1:(l+1)**2,1:nbands_dft,1:ilo(l),1:3) = cclm(l**2+1:(l+1)**2,1:nbands_dft,1:ilo(l),1:3)*cfac
       ENDDO
       !$OMP END PARALLEL DO
    END IF

    ! Chenging alm,blm,clm to make them DMFT-like
    !   alm_{lm,i}     <--  alm_{lm,j} * A^{DMFT}_{ji}
    !   blm_{lm,i}     <--  blm_{lm,j} * A^{DMFT}_{ji}
    !   clm_{lm,i,jlo} <--  clm_{lm,j,jlo} * A^{DMFT}_{ji}
    !
    ALLOCATE( xlm_tmp(lmax2lmax2,nbandsx) )
    CALL zgemm('N','N',lmax2lmax2,nbandsx,nbands,cone,alm(1,DM_nemin),lmax2lmax2,Aweight,nbands,czero,xlm_tmp,lmax2lmax2)
    alm(:,DM_nemin:DM_nemaxx)=xlm_tmp(:,:)
    CALL zgemm('N','N',lmax2lmax2,nbandsx,nbands,cone,blm(1,DM_nemin),lmax2lmax2,Aweight,nbands,czero,xlm_tmp,lmax2lmax2)
    blm(:,DM_nemin:DM_nemaxx)=xlm_tmp(:,:)
    IF (Qforce_j) THEN
       ! Chenging aalm,bblm,cclm to make them DMFT-like
       do ik=1,3
          CALL zgemm('N','N',lmax2lmax2,nbandsx,nbands,cone,aalm(1,1+DM_nemin-nemin,ik),lmax2lmax2,Aweight,nbands,czero,xlm_tmp,lmax2lmax2)
          aalm(:,DM_nemin-nemin+1:DM_nemaxx-nemin+1,ik)=xlm_tmp(:,:)
          CALL zgemm('N','N',lmax2lmax2,nbandsx,nbands,cone,bblm(1,1+DM_nemin-nemin,ik),lmax2lmax2,Aweight,nbands,czero,xlm_tmp,lmax2lmax2)
          bblm(:,DM_nemin-nemin+1:DM_nemaxx-nemin+1,ik)=xlm_tmp(:,:)
       enddo
    ENDIF
    DEALLOCATE( xlm_tmp )
    ALLOCATE( xlm_tmp(lomaxlomax,nbandsx) )
    DO jlo=1,ilo_max
       CALL zgemm('N','N',lomaxlomax,nbandsx,nbands,cone,clm(1,DM_nemin,jlo),lomaxlomax,Aweight,nbands,czero,xlm_tmp,lomaxlomax)
       clm(:,DM_nemin:DM_nemaxx,jlo)=xlm_tmp(:,:)
    ENDDO
    IF (Qforce_j) THEN
       do ik=1,3
          DO jlo=1,ilo_max
             CALL zgemm('N','N',lomaxlomax,nbandsx,nbands,cone,cclm(1,1+DM_nemin-nemin,jlo,ik),lomaxlomax,Aweight,nbands,czero,xlm_tmp,lomaxlomax)
             cclm(:,DM_nemin-nemin+1:DM_nemaxx-nemin+1,jlo,ik)=xlm_tmp(:,:)
          ENDDO
       enddo
    ENDIF
    DEALLOCATE( xlm_tmp )

    CALL cputim(time2)
    CALL walltim(time2_w)

    time_lo=time_lo+time2-time1
    time_lo_w=time_lo_w+time2_w-time1_w

  END SUBROUTINE cmp_alm_blm_clm
END module alm_blm_clm


 MODULE WANG
   implicit none
   integer, parameter:: in=10,im=2,igen=5           !! 'in' is the number of species in each patch, 'im' is the number of patches in the metacommunity
   integer, parameter:: imn=im*in, npar=imn*(imn+3)+im**2,ifact=1  !! npar is the number of parameters, ifact gives number of community pairs
   INTEGER, PARAMETER :: mlyap=0,NEQ=(mlyap+1)*imn,niter=30000,ntrans=20000,ninic=2000,net_Tscl=1000000  !! ntrans is transient time, niter is iteration beyond transients
   real(kind=4), parameter:: d_max=14000.d0, a_max=1000.d0, xL=1.d0
   real(kind=4):: par(npar)

    CONTAINS

      SUBROUTINE DERIVS(NEQ,T,Y,YDOT)
        IMPLICIT NONE
        INTEGER, intent(in):: NEQ
        real(kind=4), intent(in):: T, Y(NEQ)
        real(kind=4), intent(out):: YDOT(NEQ)
        integer(kind=4)::i,j,iti,itj,kc,icptch,iflag
        real(kind=4):: cptsum,cdsum
        real(kind=4):: ar(imn),akc(imn),acpt(imn,imn),ad(imn),xadj(im,im) !! parameters

        !*************************
        kc=0; do i=1,imn                                       !!
                kc=kc+1; ar(i)=par(kc)                         !!
              enddo                                            !! maximal growth rates
!             write(*,24),ar; pause
              do i=1,imn                                       !!
                kc=kc+1; akc(i)=par(kc)                        !!
              enddo                                            !! carrying capacities
!            write(*,24) akc; pause
              acpt=0.0                                         !! competition initialization
              do i=1,imn; do j=1,imn                           !!
                kc=kc+1; acpt(i,j)=par(kc)                     !! competition terms
              enddo; enddo                                     !!
!              write(*,24) acpt; pause
              do i=1,imn                                       !!
                kc=kc+1; ad(i)=par(kc)                         !!
              enddo                                            !! dispersal terms
           !write(*,24),ad; pause
              do i=1,im
               do j=1,im
               kc=kc+1; xadj(i,j)=par(kc)
               enddo
              enddo
            !  print*,xadj
        !*************************
        do iti=1,imn

           cptsum=0.0                                                   !!
           do itj=1,imn                                                 !! competition summation
            cptsum=cptsum+acpt(iti,itj)*Y(itj)                          !!
           enddo                                                        !!

           cdsum=0.0                                                                            !!
           icptch=int(float(iti-1)/float(in))+1                                                 !!
           iflag=mod(iti,in); if(iflag==0) iflag=in                                             !!
           kc=0
           do itj=1,im                               !! dispersal
            if(itj/=icptch) then
              cdsum=cdsum+xadj(itj,icptch)*Y(iflag+(itj-1)*in)                 !!
            endif
           !if(itj/=icptch) write(*,24) dfloat(iti),dfloat(icptch),dfloat(iflag+(itj-1)*in),ad(iti); pause
           enddo

           if (sum(xadj(:,icptch))>0.00001d0) then
             cdsum=cdsum/sum(xadj(:,icptch))
            else
             cdsum=0.d0
            endif

            YDOT(iti)=ar(iti)*Y(iti)*(1.0-(Y(iti)+cptsum)/akc(iti))+ad(iti)*(cdsum-Y(iti))

       enddo
       24  Format(30F15.7)
      RETURN
      END SUBROUTINE DERIVS
    END MODULE WANG

!******************************************************************

    PROGRAM DEMOWANG

      USE WANG
      USE random
!     Type declarations:
      IMPLICIT NONE
      real(kind=4):: h,T,TOUT, xrho1, xrho2, xvar1, xvar2, kppa
      real(kind=4):: ald,aldi,gmd,gmdi,btd, reg_avg
      real(kind=4):: alcv,gmcv,btcv
      real(kind=4):: sum_patch_l,svsum,spsum,tot_mean
      real(kind=4):: p_waxmn, p_mc
      integer(kind=4):: inic,IOUT,i,ii,j,jj,l,kcout,kij
      real(kind=4),allocatable:: Y(:),YDER(:),ctmp(:,:), noi_spec(:), noi_ptch(:), ypopl(:,:)
      real(kind=4),allocatable:: yavg(:,:),yprel(:,:),ysave(:,:,:),lo_simp(:), omg(:), sysame(:)
      real(kind=4),allocatable:: ymean(:), ysd(:), ycov(:,:)
      real(kind=4),allocatable:: adjm(:,:), beta(:), dist_mat(:,:)
      real(kind=4),allocatable:: simp_conc(:),xp_i(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! variables for the random number generator
    INTEGER, ALLOCATABLE  :: ix(:), seed(:)

    REAL, ALLOCATABLE     :: xep(:), xptch(:), xmeanep(:), xmeanptch(:), covep(:), covptch(:), chol_fep(:), &
                            chol_fptch(:), xseed(:)

    INTEGER               :: nep, nptch, ndf, n_binom, k, kep, kptch, ier, pos1, pos2, which, iindx

    REAL                  :: average, stdvn, shape, a, b, mu, p, pop_mean, pop_sd,   &
                            one = 1.0, zero = 0.0, sk, two = 2.0, pcntile(9), &
                            middle, xmax, xmin, xtseed
    LOGICAL               :: first

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   allocate(Y(NEQ),YDER(NEQ),ctmp(imn,imn),noi_spec(imn),noi_ptch(imn)); par=0.0
   allocate(yavg(im,in), yprel(im,in),lo_simp(im),sysame(in),ysave(im,in,niter))
   allocate(ymean(im),ysd(im),ycov(im,im))
   allocate(simp_conc(im),omg(im),xp_i(in))
   allocate(adjm(im,im),beta(im),dist_mat(im,im))    !!!! NETWORK PARAMETERS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! RANDOM NUMBER GENERATOR PREPARATION
!     Set the random number seed.
 CALL RANDOM_SEED(size=k)
 ALLOCATE (seed(k),xseed(k))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FIXING THE SEED
!do i=1,k; call cpu_time(xtseed); xseed(i)=xtseed; enddo
 call random_number(xseed)
 seed=int(50000*xseed)
 CALL RANDOM_SEED(put=seed)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!! STOCHASTIC PARAMETERS FROM PAPER
xrho1=0.0  !!print*, "Between species correlation";
!xrho2=0.0
read*,xrho2; print*, xrho2, "started" !! between patch correlation
xvar1=0.04  !! species specific variance
xvar2=0.04  !! patch specific variance
kppa=8.0    !! kappa

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! allocation for the random number for species specific environmental response
    kep=in                    !! species/patch
    nep = kep*(kep+1)/2
    ALLOCATE( xep(kep), xmeanep(kep), covep(nep), chol_fep(nep))

!!!!!! MEAN VECTOR
    xmeanep=0.0
! !
!!!!!! COVARIANCE MATRIX FOR SPECIES SPECIFIC ENVIRONMENTAL RESPONSE
        pos2 = 0
        DO i = 1, kep
         pos1 = pos2 + 1
         pos2 = pos2 + i
         covep(pos1:pos2)=xrho1*xvar1; covep(pos2)=xvar1
        END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! allocation for the random number for patch specific environmental response
    kptch=im
    nptch = kptch*(kptch+1)/2
    ALLOCATE ( xptch(kptch), xmeanptch(kptch), covptch(nptch), chol_fptch(nptch) )
!!!!!! MEAN VECTOR FOR PATCHES
    xmeanptch=0.0
! !
!!!!!! COVARIANCE MATRIX FOR BETWEEN PATCH ENVIRONMENTAL RESPONSE
    pos2 = 0
    DO i = 1, kptch
        pos1 = pos2 + 1
        pos2 = pos2 + i
        covptch(pos1:pos2)=xrho2*xvar2; covptch(pos2)=xvar2
    END DO

    !!!!!!!!!!!!!!!!! NETWORK STUFF

       do i=1,im
         beta(i)=0.5d0!0.3d0+0.4d0*dfloat(i)/dfloat(NEQ)  ! AREA PROXIES for Waxman
       enddo

     dist_mat=0.d0
      do i=1,im
       do j=i+1,im
        call random_number(dist_mat(j,i)); dist_mat(i,j)=dist_mat(j,i)  !! dist_mat is the distance matrix for the KHs
        enddo
      enddo
      adjm=0.d0; kij=0; write(1,*) 'digraph {'
        do i=1,im                                      !
          do j=1,im
            kij=kij+1                                          !
            p_waxmn=(beta(i)*beta(j))*exp(-dist_mat(i,j)/xL)!(beta(i)*beta(j))*dexp(-dist_mat(i,j)/xL)  !
            call random_number(p_mc)                         ! Adjacency matrix
            !if(p_mc<p_waxmn) adjm(i,j)=1.d0
            adjm(i,j)=1.d0
            if(i==j) adjm(i,j)=0.d0                           !
            if(adjm(i,j)==1) write(1,*) i,'->',j!
          enddo                                              !
         enddo
         write(1,*) '}'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! RANDOM NUMBER GEN. PREPARATION DONE
      do inic=1,ninic!; print*,inic
      par=0.0
!     Set  the problem parameters:
      call random_number(par(1:imn));
      par(1:imn)=0.5+par(1:imn)*1.0              !! M*N MAXIMAL GROWTH RATES

      call random_number(par(imn+1:2*imn))
      par(imn+1:2*imn)=0.5+par(imn+1:2*imn)*1.0  !! M*N CARRYING CAPACITIES
      ctmp=0.0
          do i=1,im
           call random_number(ctmp((i-1)*in+1:i*in,(i-1)*in+1:i*in))
           ctmp((i-1)*in+1:i*in,(i-1)*in+1:i*in)=0.25+ctmp((i-1)*in+1:i*in,(i-1)*in+1:i*in)*0.5
          if (i>1) ctmp((i-1)*in+1:i*in,(i-1)*in+1:i*in)=ctmp(1:in,1:in) !! weird addition following the Shaopeng supplied code - ask?
          enddo
          kcout=0
          do i=1,imn
           do j=1,imn
            kcout=kcout+1
            if(i==j) ctmp(i,j)=0.0
            par(2*imn+kcout)=ctmp(i,j)            !! (M*N)^2 COMPETITION COEFFICIENTS FOR M PATCHES
           enddo
          enddo
          !write(*,24)ctmp; pause
      do i=1,im
        call random_number(par(2*imn+imn**2+1+i*in-igen:2*imn+imn**2+i*in))                                         !! DISPERSAL RATES
       ! print*,par(2*imn+imn**2+1+i*in-igen:2*imn+imn**2+i*in)
        par(2*imn+imn**2+1+i*in-igen:2*imn+imn**2+i*in)=0.1+par(2*imn+imn**2+1+i*in-igen:2*imn+imn**2+i*in)*0.05         !!
        !print*,par(2*imn+imn**2+1+i*in-igen:2*imn+imn**2+i*in)
      enddo
      kcout=2*imn+imn**2+imn
      do i=1,im
        do j=1,im
         kcout=kcout+1; par(kcout)=adjm(i,j)
        enddo
      enddo
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !     Set the initial conditions:
       call random_number(Y); YDER=0.0
! !     Output points:
       T = 0.0
       h=0.01
! !     Perform the integration:
       DO IOUT = 1, ntrans
         call DERIVS(NEQ,T,Y,YDER)
         Y=Y+h*YDER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(mod(IOUT,net_Tscl)==0) then  !! NETWORK RESET
         adjm=0.d0; kij=0 !! INITIALIZE ADJACENCY MATRIX
         do i=1,im                                                             !
          do j=1,im
            kij=kij+1                                                             !
            p_waxmn=(beta(i)*beta(j))*exp(-dist_mat(i,j)/xL)                   !
            call random_number(p_mc)                         ! Adjacency matrix
            !if(p_mc<p_waxmn) adjm(i,j)=1.d0                   !
            adjm(i,j)=1.d0
            if(i==j) adjm(i,j)=0.d0                           !
          enddo                                              !
         enddo
             kcout=2*imn+imn**2+imn
              do i=1,im
               do j=1,im
                kcout=kcout+1; par(kcout)=adjm(i,j)
               enddo
              enddo
       endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       T=T+h
       END DO ! IOUT
       Y=Y+0.1
  !!! TO TAKE CARE OF ZERO POPULATIONS
       first = .true.

       DO IOUT = 1, niter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! EVOLUTION BEGINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! between species random number
        CALL random_mvnorm(kep, xmeanep, covep, chol_fep, first, xep, ier)
        IF (ier .NE. 0) THEN
            WRITE(*, *) '** Covariance matrix1 is not +ve definite **'
            PAUSE; !EXIT replaced since the do loop was removed
        END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! between patch random number
        CALL random_mvnorm(kptch, xmeanptch, covptch, chol_fptch, first, xptch, ier)
        IF (ier .NE. 0) THEN
        WRITE(*, *) '** Covariance matrix2 is not +ve definite **'
        PAUSE; !EXIT replaced since the do loop was removed
        END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        do i=1,im
           noi_spec((i-1)*in+1:i*in)=xep
        enddo         !!!!!!!!!!!!!noise arrangement for integration

        kcout=0
        do i=1,im
          do j=1,in
           kcout=kcout+1
           noi_ptch(kcout)=xptch(i)
          enddo
        enddo
         call DERIVS(NEQ,T,Y,YDER)
         do i=1,imn
           Y(i)=Y(i)+h*YDER(i)+sqrt(h)*Y(i)*(noi_spec(i)+noi_ptch(i) &
           & +kppa*noi_ptch(i)*noi_spec(i)) !!  noise added; T=T+h
         enddo
        kcout=0
        do l=1,im
         do i=1,in
           kcout=kcout+1
           ysave(l,i,IOUT)=Y(kcout)
          enddo
        enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(mod(IOUT,net_Tscl)==0) then  !! NETWORK RESET
         adjm=0.d0; kij=0 !! INITIALIZE ADJACENCY MATRIX
         do i=1,im                                                             !
          do j=1,im
            kij=kij+1                                                             !
            p_waxmn=(beta(i)*beta(j))*exp(-dist_mat(i,j)/xL)                   !
            call random_number(p_mc)                         ! Adjacency matrix
            !if(p_mc<p_waxmn) adjm(i,j)=1.d0                   !
            adjm(i,j)=1.d0
            if(i==j) adjm(i,j)=0.d0                           !
          enddo                                              !
         enddo
             kcout=2*imn+imn**2+imn
              do i=1,im
               do j=1,im
                kcout=kcout+1; par(kcout)=adjm(i,j)
               enddo
              enddo
       endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        T=T+h
       END DO ! IOUT
!!!!!!!!!!!!!!!!!!!!!!!! ENS based
   do l=1,im
     sum_patch_l=sum(ysave(l,:,niter))
    do i=1,in
     yprel(l,i)=ysave(l,i,niter)/sum_patch_l
    enddo
   enddo

  aldi=0.0
    do l=1,im
    simp_conc(l)=sum(yprel(l,:)**2)
    omg(l)=sum(ysave(l,:,niter))/sum(ysave(:,:,niter))
      aldi=aldi+omg(l)*simp_conc(l)
    enddo

    ald=1.0/aldi                              !! alpha_d

    do i=1,in
     xp_i(i)=sum(ysave(:,i,niter))/sum(ysave(:,:,niter))
    enddo

    gmdi=sum(xp_i**2) !! gamma_d inverse
    gmd=1.0/gmdi                                !! gamma_d
    btd=gmd/ald                                  !! beta_d

! !!!!! VARIABILITY BASED MEASURES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      svsum=0.0
        do i=1,im        !!!! temporal mean for each patch
        spsum=0.0
         do IOUT=1,niter
          spsum=spsum+sum(ysave(i,:,IOUT))
         enddo
         ymean(i)=spsum/float(niter) !! temporal mean of community biomass in patch i
         svsum=svsum+ymean(i)
        enddo          !                                       QC done
        tot_mean=svsum!/float(niter)   ! temporal mean of metacommunity biomass
        reg_avg=tot_mean/float(im)    !regional average of temporal means of local community biomass

         ycov=0.d0
            do IOUT=1,niter
             do i=1,im
              do j=1,im
               ycov(i,j)=ycov(i,j)+(sum(ysave(i,:,IOUT))-ymean(i))*(sum(ysave(j,:,IOUT))-ymean(j))
              enddo
             enddo
            enddo
              ycov=ycov/float(niter)  !! covariance matrix between community biomasses
            do i=1,im
             ysd(i)=sqrt(ycov(i,i))
            enddo
      !      print*,ycov; pause
        alcv=(sum(ysd)/tot_mean)**2                  !! alpha_cv
        gmcv=sum(ycov)/(tot_mean)**2                 !! gamma_cv
        btcv=alcv/gmcv  !! is it correct?            !! beta_cv
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        write(11,24) ald, btd, gmd, alcv, btcv, gmcv
        write(*,24) ald, btd, gmd, alcv, btcv, gmcv
   enddo !!!! INIC
!! OUTPUT FILE FORMAT
      24  Format(20F20.7)

END PROGRAM DEMOWANG


text:
subroutine sldsup_turnonModel_NH(a)
   use typre
   use Mod_SUPSolids_NH
   implicit none
   class(SUPSolidsProblem_NH) :: a

       a%sld_type = 'NONLI'
       a%sld_model ='NEOHO'

end subroutine sldsup_turnonModel_NH

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    CompDrivingForce.f90
    !> \brief   Compute the driving force of all pure condensed phases.
    !> \author  M.H.A. Piro
    !> \date    Apr. 26, 2012
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer          Description of change
    !   ----            ----------          ---------------------
    !   04/26/2012      M.H.A. Piro         Original code
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to compute the driving force of all pure condensed phases in
    !! the database.  The driving force is defined as the difference between the standard molar Gibbs energy of
    !! a pure condensed phase and the corresponding value computed from the element potentials.  This value is
    !! used to determine whether a pure condensed phase should be added to the system.  For a more thorough
    !! explanation of the chemical significance of the driving force, refer to the following literature:
    !!
    !! H.L. Lukas, S.G. Fries, B. Sundman, Computational Thermodynamics - The Calphad Method, Cambridge
    !! University Press, New York, 2007.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    !> \param[out]  iMaxDrivingForce    An integer scalar representing the index of the pure condensed phase
    !!                                   with the maximum driving force.  A value of zero is returned by default.
    !> \param[out]  dMaxDrivingForce    A double real scalar representing the maximum driving force of all pure
    !!                                   condensed phases.  A value of zero is returned by default.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine CompDrivingForce(iMaxDrivingForce,dMaxDrivingForce)

    USE ModuleThermo
    USE ModuleGEMSolver

    implicit none

    integer :: i, j, iMaxDrivingForce
    real(8) :: dTemp, dMaxDrivingForce

    ! Initialize variables:
    iMaxDrivingForce = 0
    dMaxDrivingForce = 0D0

    ! Loop through all pure condensed phases:
    do i = nSpeciesPhase(nSolnPhasesSys) + 1, nSpecies - nDummySpecies

        ! Compute the chemical potential of this phase as defined by the element potentials:
        dTemp = 0D0
        do j = 1, nElements
            dTemp = dTemp + dElementPotential(j) * dStoichSpecies(i,j)
        end do

        ! Compute the driving force:
        dTemp = dStdGibbsEnergy(i) - dTemp

        ! Normalize per gram-atom:
        dTemp = dTemp / dSpeciesTotalAtoms(i)

        ! Check for most negative value:
        if (dTemp < dMaxDrivingForce) then
            iMaxDrivingForce = i
            dMaxDrivingForce = dTemp
        end if

    end do

    return

end subroutine CompDrivingForce

max_stars_repo_name: nrc-fuels/thermochimica

max_issues_repo_name: nrc-fuels/thermochimica

max_forks_repo_name: nrc-fuels/thermochimica


text:
       subroutine srfflx( lin, &
                          lout, &
                          new_nq, &
                          new_solsym, &
                          srf_flx_map, &
                          srf_flx_cnt, &
			  tag )

      use var_mod, only : var_lim

      implicit none

!-----------------------------------------------------------------------
!       ... Dummy args
!-----------------------------------------------------------------------
      integer, intent(in)  ::  lin, &       ! input unit number
                               lout, &      ! output unit number
                               new_nq, &    ! species count
                               tag          ! emission or deposition tag ( 1,2 )
      integer, intent(inout) ::  srf_flx_cnt  ! count of species with srf flux
      integer, intent(out) ::  srf_flx_map(*)  ! srf flux "map"
      character(len=8), intent(in) ::  new_solsym(*) ! species names

!-----------------------------------------------------------------------
!       ... Local variables
!-----------------------------------------------------------------------
      integer  ::  nchar
      integer  ::  toklen(20)
      integer  ::  j, k
      integer  ::  no_tokens

      character(len=320) ::  buff
      character(len=320) ::  buffh
      character(len=8)  ::  tokens(20)

      logical  ::  found

      integer, parameter ::  symlen = 8

!-----------------------------------------------------------------------
!       ... Read the surface flux species
!-----------------------------------------------------------------------
      do
         call cardin( lin, buff, nchar )
         buffh = buff
         call upcase( buffh )
         if( tag == 1 .and. buffh == 'ENDSURFACEFLUX' ) then
	    exit
         else if( tag == 2 .and. buffh == 'ENDSURFACEDEPOSITION' ) then
	    exit
	 end if
         call gettokens(  buff,     nchar,    ',',      symlen, &
                          tokens,   toklen,   20,       no_tokens )
         if( no_tokens == 0 ) then
            call errmes( ' SRFFLX: Species input line in error@', lout, buff, 1, buff )
         end if
         do j = 1,no_tokens
            srf_flx_cnt = srf_flx_cnt + 1
            if( srf_flx_cnt > var_lim ) then
               call errmes( ' SRFFLX: Species count exceeds limit@', lout, buff, 1, buff )
            end if
	    found = .false.
	    do k = 1,new_nq
	       if( tokens(j) == new_solsym(k) ) then
		  if( srf_flx_cnt > 1 ) then
		     if( any( srf_flx_map(:srf_flx_cnt) == k ) ) then
			if( tag == 1 ) then
                           call errmes( '# is already in srf emis list@', lout, tokens(j), toklen(j), buff )
			else if( tag == 2 ) then
                           call errmes( '# is already in dry dep list@', lout, tokens(j), toklen(j), buff )
			end if
		     end if
		  end if
	          srf_flx_map(srf_flx_cnt) = k
	          found = .true.
		  exit
	       end if
	    end do
	    if( .not. found ) then
               call errmes( '# is not in solution species list@', lout, tokens(j), toklen(j), buff )
            end if
         end do
      end do

      end subroutine srfflx

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_licenses: NCSA Unlicense


text:
module multid_slope_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  private

  public multid_slope

contains

! :::
! ::: ------------------------------------------------------------------
! :::

  subroutine multid_slope(q, q_lo, q_hi, ncomp, nc, &
                          flatn, &
                          dqx, dqy, qpd_lo, qpd_hi, &
                          dx, dy, &
                          ilo, jlo, ihi, jhi)

    use bl_constants_module

    use amrex_fort_module, only : rt => amrex_real
    implicit none

    integer, intent(in) :: q_lo(3), q_hi(3)
    integer, intent(in) :: qpd_lo(3), qpd_hi(3)
    integer, intent(in) :: ilo, jlo, ihi, jhi
    integer, intent(in) :: ncomp, nc

    real(rt)        , intent(in) :: dx, dy

    real(rt)        , intent(in) ::     q( q_lo(1): q_hi(1), q_lo(2): q_hi(2), ncomp)
    real(rt)        , intent(in) :: flatn( q_lo(1): q_hi(1), q_lo(2): q_hi(2))

    real(rt)        , intent(inout) :: dqx(qpd_lo(1):qpd_hi(1),qpd_lo(2):qpd_hi(2), ncomp)
    real(rt)        , intent(inout) :: dqy(qpd_lo(1):qpd_hi(1),qpd_lo(2):qpd_hi(2), ncomp)

    integer :: i, j, m, niter

    real(rt)        , allocatable :: q_nd(:,:)

    integer, parameter :: ill = 1, ilh = 2, irl = 3, irh = 4
    real(rt)         :: ss(4), ss_temp(4), min_ss(4), max_ss(4), diff(4)
    real(rt)         :: A_x, A_y, A_xy
    real(rt)         :: sumdiff, sgndiff, redfac, redmax, div
    real(rt)        , parameter :: eps = 1.e-10_rt
    integer :: kdp
    integer, parameter :: MAX_ITER = 3

    ! here we do a multidimensional reconstruction of the data, forming
    ! a bilinear polynomial:
    !     q(x,y) = A_xy (x-x_i)(y-y_i) + A_x (x-x_i) + A_y (y-y_i) + A_avg
    !
    ! First find the values at the nodes.  Here q_nd(i,j) will refer to
    ! a_{i-1/2,j-1/2}

    allocate(q_nd(q_lo(1):q_hi(1),q_lo(2):q_hi(2)))

    do j = jlo-2, jhi+3
       do i = ilo-2, ihi+3

          q_nd(i,j) = &
             (         q(i-2,j-2,nc) -  7.0e0_rt*(q(i-1,j-2,nc) + q(i,j-2,nc)) +       q(i+1,j-2,nc) + &
              (-7.0e0_rt)*q(i-2,j-1,nc) + 49.0e0_rt*(q(i-1,j-1,nc) + q(i,j-1,nc)) - 7.0e0_rt*q(i+1,j-1,nc) + &
              (-7.0e0_rt)*q(i-2,j  ,nc) + 49.0e0_rt*(q(i-1,j  ,nc) + q(i,j  ,nc)) - 7.0e0_rt*q(i+1,j  ,nc) + &
                       q(i-2,j+1,nc) -  7.0e0_rt*(q(i-1,j+1,nc) + q(i,j+1,nc)) +       q(i+1,j+1,nc))/ &
             144.0e0_rt

       enddo
    enddo

    ! now get the coefficients of the bilinear polynomial.  Here we loop
    ! ! over zone centers.
    do j = jlo-1, jhi+1
       do i = ilo-1, ihi+1

          ss(ill) = q_nd(i,j)     ! a_{i-1/2,j-1/2}
          ss(ilh) = q_nd(i,j+1)   ! a_{i-1/2,j+1/2}
          ss(irl) = q_nd(i+1,j)   ! a_{i+1/2,j-1/2}
          ss(irh) = q_nd(i+1,j+1) ! a_{i+1/2,j+1/2}

          A_x  = ((ss(irh) + ss(irl)) - (ss(ilh) + ss(ill)))/(2.0e0_rt*dx)
          A_y  = ((ss(ilh) + ss(irh)) - (ss(ill) + ss(irl)))/(2.0e0_rt*dy)
          A_xy = ((ss(irh) - ss(irl)) - (ss(ilh) - ss(ill)))/(dx*dy)

          ! check if we are within the cc-values
          ss_temp(ill) = q(i,j,nc) - HALF*dx*A_x - HALF*dy*A_y + 0.25e0_rt*dx*dy*A_xy
          ss_temp(ilh) = q(i,j,nc) - HALF*dx*A_x + HALF*dy*A_y - 0.25e0_rt*dx*dy*A_xy
          ss_temp(irl) = q(i,j,nc) + HALF*dx*A_x - HALF*dy*A_y - 0.25e0_rt*dx*dy*A_xy
          ss_temp(irh) = q(i,j,nc) + HALF*dx*A_x + HALF*dy*A_y + 0.25e0_rt*dx*dy*A_xy

          min_ss(ill) = min(q(i-1,j-1,nc), q(i,j-1,nc), q(i-1,j,nc), q(i,j,nc))
          max_ss(ill) = max(q(i-1,j-1,nc), q(i,j-1,nc), q(i-1,j,nc), q(i,j,nc))

          min_ss(ilh) = min(q(i-1,j+1,nc), q(i,j+1,nc), q(i-1,j,nc), q(i,j,nc))
          max_ss(ilh) = max(q(i-1,j+1,nc), q(i,j+1,nc), q(i-1,j,nc), q(i,j,nc))

          min_ss(irl) = min(q(i+1,j-1,nc), q(i,j-1,nc), q(i+1,j,nc), q(i,j,nc))
          max_ss(irl) = max(q(i+1,j-1,nc), q(i,j-1,nc), q(i+1,j,nc), q(i,j,nc))

          min_ss(irh) = min(q(i+1,j+1,nc), q(i,j+1,nc), q(i+1,j,nc), q(i,j,nc))
          max_ss(irh) = max(q(i+1,j+1,nc), q(i,j+1,nc), q(i+1,j,nc), q(i,j,nc))

          ! limit
          do m = 1, 4
             ss_temp(m) = max(min(ss_temp(m), max_ss(m)), min_ss(m))
          enddo

          do niter = 1, MAX_ITER
             sumdiff = (ss_temp(ill) + ss_temp(ilh) + &
                        ss_temp(irl) + ss_temp(irh)) - 4.0e0_rt*q(i,j,nc)
             sgndiff = sign(1.e0_rt,sumdiff)

             do m = 1, 4
                diff(m) = (ss_temp(m) - q(i,j,nc))*sgndiff
             enddo

             kdp = 0

             do m = 1, 4
                if (diff(m) > eps) kdp = kdp + 1
             enddo

             do m = 1, 4
                if (kdp < 1) then
                   div = 1.e0_rt
                else
                   div = dble(kdp)
                endif

                if (diff(m) > eps) then
                   redfac = sumdiff*sgndiff/div
                   kdp = kdp-1
                else
                   redfac = 0.0e0_rt
                endif

                if (sgndiff > 0.0e0_rt) then
                   redmax = ss_temp(m) - min_ss(m)
                else
                   redmax = max_ss(m) - ss_temp(m)
                endif

                redfac = min(redfac, redmax)
                sumdiff = sumdiff - redfac*sgndiff
                ss_temp(m) = ss_temp(m) - redfac*sgndiff
             enddo
          enddo

          ! construct the final slopes
          A_x  = ((ss_temp(irh) + ss_temp(irl)) - &
                  (ss_temp(ilh) + ss_temp(ill)))/(2.0e0_rt*dx)

          A_y  = ((ss_temp(ilh) + ss_temp(irh)) - &
                  (ss_temp(ill) + ss_temp(irl)))/(2.0e0_rt*dy)

          A_xy = ((ss_temp(irh) - ss_temp(irl)) - &
                  (ss_temp(ilh) - ss_temp(ill)))/(dx*dy)

          ! now construct the limited differences in x and y
          dqx(i,j,nc) = flatn(i,j)*A_x*dx
          dqy(i,j,nc) = flatn(i,j)*A_y*dy

       enddo
    enddo

    deallocate(q_nd)

  end subroutine multid_slope

end module multid_slope_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
module param
    parameter nu = 10
    real*8 U0, dUdx
    real*8 U(nu)
end module param

module mtr
    parameter (nx = 100)
    parameter nL = 102, nt = 8
    real*8 Nd(-nx - 1:nx + 1), x(-nx - 1:nx + 1)
    real*8 xL(nl), Yl(nl, nt), Y1(nl), ts(nt), xmax
    integer JF(nt), it, its
    real*8 Xx(50, nt), CSD(50, nt)
end module mtr

program main
end program main

max_stars_repo_name: leginee/netbeans

max_issues_repo_name: leginee/netbeans

max_forks_repo_name: leginee/netbeans


!***********************************************************************************************************************************
! UNIT:
!
!  (MODULE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  MODULE MR_MOD_FUNC_DS

    USE MR_KINDS

    USE MR_DEF_CONSTS_N_REF_PARS

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: MR_FUNC_DS

    INTERFACE MR_FUNC_DS
      MODULE PROCEDURE MR_FUNC_DS_PARD_KIND
      MODULE PROCEDURE MR_FUNC_DS_FDRD_KIND
    END INTERFACE

!***********************************************************************************************************************************

  CONTAINS

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  FUNCTION MR_FUNC_DS_PARD_KIND( D0 ) RESULT( DS )

    IMPLICIT NONE

    REAL(PARD_KIND) , INTENT(IN ) :: D0

    REAL(PARD_KIND) :: DS

    DS = D0 / ( DSPAR * ZR )

  END FUNCTION MR_FUNC_DS_PARD_KIND

!***********************************************************************************************************************************
! UNIT:
!
!  (SUBROUTINE)
!
! PURPOSE:
!
!   TO
!
! DEFINITION OF VARIABLES:
!
!
!
! RECORD OF REVISIONS:
!
!      DATE       |    PROGRAMMER    |    DESCRIPTION OF CHANGE
!      ====       |    ==========    |    =====================
!   20XX-XX-XX    |     DR. HYDE     |    ORIGINAL CODE.
!
!***********************************************************************************************************************************
  FUNCTION MR_FUNC_DS_FDRD_KIND( D0 ) RESULT( DS )

    IMPLICIT NONE

    REAL(FDRD_KIND) , INTENT(IN ) :: D0

    REAL(FDRD_KIND) :: DS

    DS = D0 / ( DSPAR * ZR )

  END FUNCTION MR_FUNC_DS_FDRD_KIND

  END MODULE MR_MOD_FUNC_DS


text:
Spiders are ubiquitous in Davis. Cellar spiders, aka Daddy Long Legs, are probably the most common and least harmful, though the Black Widow spiders have also been noted in the area. Brown recluse spiders are not found in Davis; people commonly misidentify brown colored spiders as Brown recluse spiders.

Some people think of spiders as pests, but the truth is that spiders feed on flies, mosquitoes, cockroaches, etc., thus keeping populations of exterminators many others pests in check.

Spiders are arthropods of the class Arachnida and the order Araneae. They are related to mites, ticks and scorpions, although they differ from them in that they have two body segments, all have poison glands in their chelicera, or jaws, and have a spinneret, or silk producing gland on their abdomen. Some types spin webs, some use their webs as bolas, some use their webs to create trap doors. Some spiders rarely use their webs at all.  The largest group of spiders is the Salticidae, or the jumping spiders.

There are many, many different spiders to be found outside, from flower spiders, to jumping spiders, and even fairly large garden spiders, which can be over an inch long from tip of cephalothorax to end of abdomen.

Indoors, the most common spider is probably the Daddy Long Legs (not to be mistaken for the Harvestman arachnid) or Cellar spider, spiders of the family Pholcidae. They can also be identified with the cobwebby, hammocklike webs they commonly spin in houses. You probably have this spider in your apartment. Contrary to urban myth, these spiders are completely harmless to humans, whether alive or ground up into a powder. A possible domestic spider also includes the Black Widow Spider, a member of the genus Latrodectus. Unlike the aforementioned Cellar spider, the females pack a potent venom capable of harming an adult human, though it is not usually fatal. They are identified by their black or brown color, and large abdomen. Females also have the red hourglass on their underside which identifies them as widows. They commonly can inhabit sheds and other generally disused structures, and lay hanging sacs of eggs. Some places I have found widows  include under the cement bike lock structures that are all over the campus, and the old shed near the Bee Biology building. Davis is full of black widows. I have several times found them under the lids of the recycling containers here. Theyre pretty shy, dont like light, and move relatively slowly. Davis is also home to wolf spiders and jumping spiders.

Spiders can be found almost anywhere where insects, their primary prey, can be found. In other words, almost anywhere. Spiderlings can use their silk to balloon long distances, and distribute themselves over a wide range.

 Species

  Black Widow Black Widows
  wiki:WikiPedia:Daddy_longlegs_spider Daddy Long Legs
  wiki:WikiPedia:Phidippus_audax Daring Jumping Spider

 Images

These are all local species, but not all photos are of local spiders.  If you see a good example of one, please snap a picture so we can see our local eight legged beasties.

  Im not sure this is true. Could someone with better knowledge remove nonlocal spiders below? It would be very helpful. Users/NumiaCairaguas
    Davis is in the range for all the spiders listed below.  Ask a professor?

Maybe someone wants to classify these (taken in Users/WesHardaker my yard):

Found this guy creeping on the carpet when I first moved into my apartment last fall  Users/JournaL

I found this generic brown spider in my bathtub.  I feel like these are pretty common here, but I dont know what kind of spider it is.

To learn more about creepy crawlies inhabiting Davis our town, please take a look at our Town Wildlife.

20070327 12:53:45 nbsp My scout master back in the day made us totally fear the brown recluse.  Those things will mess you up.  If you dont die your flesh like turns to liquid.  Nasty things! Users/BradBenedict

20070327 15:18:17 nbsp Professor Kimsey (he teaches the zoology class, BIS 1B) told me that Brown Recluse Spiders dont live in California. I want to see references that say otherwise, or else that spider should be removed from this page. Users/NumiaCairaguas
   I know multiple people who have been bitten by Brown Recluse in Ca Daubert
      Kimsey was right.  See http://en.wikipedia.org/wiki/Brown_recluse_spider#Distribution  , or AlexMandels link in the comment thrad below.   Its a very, very common mistake, so I doubt they actually were the brown recluse spiders. users/edwins ES

20070327 15:18:53 nbsp ...I forgot to say please. :) Users/NumiaCairaguas

20070327 15:21:32 nbsp I looked up info on BRs a while back because my wife thought she saw one...  everything thing I found said they didnt exist in central CA (though some had been supposedly seen in southern CA, those references were really weak) Users/WesHardaker
  I looked this up a while ago too, and had found that the natural range of Brown Recluse did not include much of California if any part (South Eastern near Arizona was the most likely). I did however here recently that they have been spotted in CA as invasive species brought in by accident in some places. Its unclear wether there is a self sustaining population anywhere in the state though. This http://www.ipm.ucdavis.edu/PMG/PESTNOTES/pn7468.html UCR ANR page appears to have the most comprehensive and reliable information. Users/AlexMandel

20070717 09:15:20 nbsp I snapped a shot of that black and white jumping spider. It thought it was really pretty, but it seems there is a better picture up here. Does anyone know if it has a common name? Users/CarlosOverstreet

 20070718 07:53:13 nbsp The wikipedia link under the picture, assuming youre talking about the same one Im thinking of, says its a jumping spider.  And if you dont like them the only place you can live is greenland... Users/WesHardaker

20130209 11:59:12 nbsp Want to participate in a local Citizen Science Project? As of 2012 Explorit Science Center has been doing a Citizen Science Project concerning spider distribution in Yolo County. Anyone who takes photos anywhere in Yolo County is invited to participate and help build up the database. Participants are asked to register at http://www.explorit.org/csp/registration to get an ID number as an official participant and then sign up with Explorits Where Is My Spider Project at iNaturalist.org at http://www.inaturalist.org/projects/whereismyspider. A pilot was run in 2012 and as of 2013 the project is starting over using iNaturalist. There are smart phone apps for recording observations on iNaturalist! Users/AnneHance


text:
!-*- mode: compilation; default-directory: "/tmp/" -*-
!Compilation started at Sun May 19 23:14:14
!
!a=./F && make $a && $a < unixdict.txt
!f95 -Wall -ffree-form F.F -o F
!101
!110010
!10001100101000
!
!Compilation finished at Sun May 19 23:14:14
!
!
!   tobin=: -.&' '@":@#:
!   tobin 5
!101
!   tobin 50
!110010
!   tobin 9000
!10001100101000

program bits
  implicit none
  integer, dimension(3) :: a
  integer :: i
  data a/5,50,9000/
  do i = 1, 3
    call s(a(i))
  enddo

contains

  subroutine s(a)
    integer, intent(in) :: a
    integer :: i
    if (a .eq. 0) then
      write(6,'(a)')'0'
      return
    endif
    do i = 31, 0, -1
      if (btest(a, i)) exit
    enddo
    do while (0 .lt. i)
      if (btest(a, i)) then
        write(6,'(a)',advance='no')'1'
      else
        write(6,'(a)',advance='no')'0'
      endif
      i = i-1
    enddo
    if (btest(a, i)) then
      write(6,'(a)')'1'
    else
      write(6,'(a)')'0'
    endif
  end subroutine s

end program bits

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


submodule (pathlib) pathlib_windows

implicit none (type, external)

contains

module procedure is_absolute

character :: f

character(:), allocatable  :: buf

buf = expanduser(path)

is_absolute = .false.
if(len_trim(buf) < 2) return

f = buf(1:1)

is_absolute = (((f >= "a" .and. f <= "z") .or. (f >= "A" .and. f <= "Z")) .and. &
  buf(2:2) == ":")
!! NEED all these parentheses

end procedure is_absolute

module procedure copyfile

integer :: i,j
!! https://docs.microsoft.com/en-us/windows-server/administration/windows-commands/copy
character(*), parameter :: CMD='copy /y '

character(:), allocatable  :: s,d

s = filesep_windows(expanduser(source))
d = filesep_windows(expanduser(dest))

call execute_command_line(CMD // s // ' ' // d, exitstat=i, cmdstat=j)
if (i /= 0 .or. j /= 0) error stop "could not copy " // source // " => " // dest

end procedure copyfile

module procedure mkdir
!! create a directory, with parents if needed
integer :: i,j
!! https://docs.microsoft.com/en-us/windows-server/administration/windows-commands/md
character(*), parameter :: CMD='mkdir '
character(:), allocatable  :: buf

buf = expanduser(path)

if(directory_exists(buf)) return

call execute_command_line(CMD // filesep_windows(buf), exitstat=i, cmdstat=j)
if (i /= 0 .or. j /= 0) error stop "could not create directory " // path

end procedure mkdir

end submodule pathlib_windows


text:
www.ajspainting.net
http://ajspainting.wikispot.org/Front_Page
                                                                                Making Your  Dreams a Reality
                                                                                Residential                                 Industrial                                    Commercial
Welcome to AJS PAINTING!
If you are in need of professional painting services for your residence or business, please contact AJS PAINTING today.  Our skilled painters know from years of experience how to tackle all types of commercial, industrial or residential painting projects, including interior and exterior house painting. With every project, were devoted to providing the highest quality work and professional painting service in the area. Whatever your needs, big or small, AJS PAINTING will provide a professional quality paint job and an experience youll be thrilled with.