text:
      subroutine table_element(cmohex,
     &                         ipxintrp,ipyintrp,ipzintrp,nump,
     &                         ipielement,
     &                         ierr2)
C
      implicit none
C
      character*132 logmess
C
C
C#######################################################################
C
C      PURPOSE -
C
C         THIS ROUTINE TAKES A SET OF POINTS IN X,Y,Z AND, FOR EACH
C         POINT, FINDS AN INTERPOLATED FIELD VALUE FOR SOME SCALAR
C         FIELD ON A  GRID IN X,Y,Z.  THE METHOD IS TO
C         FIND WHICH ELEMENT IN LIST THE POINT LIES INSIDE AND THEN
C         LINEARLY INTERPOLATE IN THAT ELEMENT TO GET THE RESULTING
C         VALUE.  THE INTERPOLATION IS APPLIED TO EITHER THE FIELD
C         OR TO SOME FUNCTION OF THE FIELD, AS DETERMINED BY ITYPE.
C         IF A FUNCTION OF FIELD IS USED, THE RESULT IS FED INTO THE
C         INVERSE FUNCTION.
C
C         FORMAT:
C
C      INPUT ARGUMENTS -
C
C        cmotable - THE NAME OF THE CMO CONTAINING THE TABLE
C
C        ipxintrp - POINTER TO X-COORDINATES FOR INTERPOLATION POINTS
C
C        ipyintrp - POINTER TO Y-COORDINATES FOR INTERPOLATION POINTS
C
C        ipzintrp - POINTER TO Z-COORDINATES FOR INTERPOLATION POINTS
C
C        nump     - NUMBER OF POINTS TO BE INTERPOLATED
C
C
C     OUTPUT ARGUMENTS -
C
C        ipielement - POINTER TO THE LIST OF OUTPUT ELEMENTS
C
C        ierr2    - ERROR INDICATOR, RETURNS PLACE IN LIST WHERE THERE
C                   WAS NO BOUNDING TET IN THE LIST FOR THIS POINT.
C                   RETURNS 0 IF NO TPROBLEMS, <0 IF SOME OTHER PROBLEM.
C
C     CHANGE HISTORY -
C
C        $Log: table_element.f,v $
C        Revision 2.00  2007/11/09 20:04:04  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.10   30 Jan 2002 10:28:00   dcg
CPVCS    remove duplicate declarations - put declarations in order
CPVCS
CPVCS       Rev 1.9   17 Dec 2001 12:50:42   dcg
CPVCS    implicit none plus clean up
CPVCS
CPVCS       Rev 1.8   Mon Apr 14 17:02:28 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.7   Mon Apr 07 14:54:16 1997   het
CPVCS    Add the element to element interpoation algorithm.
CPVCS
CPVCS       Rev 1.6   Fri Jan 24 13:54:12 1997   het
CPVCS    Add the AMR arrays to the table lookup scheme to only look
CPVCS       at active elements.
CPVCS
CPVCS       Rev 1.5   Mon Dec 09 09:01:44 1996   het
CPVCS    Get the active set of elements.
CPVCS
CPVCS       Rev 1.4   Mon Dec 02 08:56:20 1996   het
CPVCS
CPVCS       Rev 1.3   Tue Nov 26 13:49:52 1996   het
CPVCS    Fix an error with calculating volumes of hybrid pri and pyr elements
CPVCS
CPVCS       Rev 1.2   Thu Nov 21 19:06:54 1996   het
CPVCS    Start rays at the element nearest the centroid of the points set.
CPVCS
CPVCS       Rev 1.1   Mon Nov 11 20:57:04 1996   het
CPVCS    Add pri and pyr element types.
CPVCS
CPVCS       Rev 1.0   Mon Jul 29 15:42:42 1996   dcg
CPVCS    Initial revision.
C
C#######################################################################
C
C
      include 'local_element.h'
      include 'consts.h'
C
      character*32 cmo, isubname
C
      pointer (ipitp1, itp1)
      integer itp1(*)
C
      pointer (ipxic_table, xic_table)
      pointer (ipyic_table, yic_table)
      pointer (ipzic_table, zic_table)
      real*8 xic_table(*), yic_table(*), zic_table(*)
C
      pointer (ipxic_value, xic_value)
      pointer (ipyic_value, yic_value)
      pointer (ipzic_value, zic_value)
      real*8 xic_value(*), yic_value(*), zic_value(*)
C
      pointer (ipitetclr, itetclr)
      pointer (ipitettyp, itettyp)
      pointer (ipitetoff, itetoff)
      pointer (ipjtetoff, jtetoff)
      integer itetclr(*), itettyp(*),
     *        itetoff(*), jtetoff(*)
      pointer (ipitet, itet1)
      pointer (ipjtet, jtet1)
      integer itet1(*), jtet1(*)
C
      pointer (ipitflag, itflag)
      integer itflag(*)
C
      pointer (ipitpath, itpath)
      integer itpath(*)
C
      pointer (ipireal1, ireal1)
      integer ireal1(*)
C
      pointer (ipitetkid, itetkid)
      integer itetkid(*)
C
      pointer (ipitactive, itactive)
      integer itactive(*)
C
      real*8 volface(100)
C
      integer iclrpath
      data iclrpath / 0 /
C
      real*8 xsmall
      data xsmall / 1.0d-06 /
C
      character*32 coption
      data coption / '-active-' /
C
C
C######################################################################
C
      integer ierr2,nnodes,ilen,ityp,ierror,length,ierrdum,it,iclose,
     *  ittest,ismax,index,if4,if1,if3,if2,j,nef,nen,
     *  i1,i2,i3,i4,i5,i6,i7,i8,in_element,ipath,iflag,irestart,
     *  iconvex,itnext,ilost,ntactive,i,idum,icharlnf,mbndry,
     *  nefcmo,nelements,nump,icscode,ierr,itstart,ip,iprintcount
      real*8 dzt,dyt,dxt,dz,dy,dx,volelm,x1,x2,x3,x4,x5,x6,x7,x8,
     *  y1,y2,y3,y4,y5,y6,y7,y8,z1,z2,z3,z4,z5,z6,z7,z8,za,ya,xa,
     *  xavg,yavg,zavg,distmin,xavg_value,yavg_value,zavg_value,
     *  dist,tepsilon
      pointer (ipxintrp, xintrp)
      pointer (ipyintrp, yintrp)
      pointer (ipzintrp, zintrp)
      real*8 xintrp(nump),yintrp(nump),zintrp(nump)
      pointer (ipielement, ielement)
      integer ielement(nump)
      character*(*) cmohex
C
C
C######################################################################
C
C
      isubname='intrphex'
      cmo=cmohex
C
      ierr2 = 0
      tepsilon = 1.0d-10
C
      call cmo_get_info('nnodes',cmo,nnodes,ilen,ityp,ierror)
      call cmo_get_info('nelements',cmo,nelements,ilen,ityp,ierror)
      call cmo_get_info('faces_per_element',cmo,nefcmo,ilen,ityp,ierror)
      call cmo_get_info('mbndry',cmo,mbndry,ilen,ityp,ierror)
C
      call cmo_get_info('itp1',cmo,ipitp1,ilen,ityp,ierror)
      call cmo_get_info('xic',cmo,ipxic_table,ilen,ityp,ierror)
      call cmo_get_info('yic',cmo,ipyic_table,ilen,ityp,ierror)
      call cmo_get_info('zic',cmo,ipzic_table,ilen,ityp,ierror)
C
      length=nnodes
      call mmgetblk("ireal1",isubname,ipireal1,length,1,icscode)
C
      call unpacktp("allreal","set",nnodes,ipitp1,ipireal1,ierrdum)
      if(ierrdum.ne.0) call x3d_error('unpacktp', isubname)
C
      ipxic_value=ipxintrp
      ipyic_value=ipyintrp
      ipzic_value=ipzintrp
C
      call cmo_get_info('itetclr',cmo,
     *                        ipitetclr,ilen,ityp,ierr)
      call cmo_get_info('itettyp',cmo,
     *                        ipitettyp,ilen,ityp,ierr)
      call cmo_get_info('itetoff',cmo,
     *                        ipitetoff,ilen,ityp,ierr)
      call cmo_get_info('jtetoff',cmo,
     *                        ipjtetoff,ilen,ityp,ierr)
      call cmo_get_info('itet',cmo,ipitet,ilen,ityp,ierr)
      call cmo_get_info('jtet',cmo,ipjtet,ilen,ityp,ierr)
C
      call mmfindbk('itetkid',cmo,ipitetkid,length,icscode)
      if(icscode.ne.0) then
         length=nelements
         call mmgetblk('itetkid',isubname,ipitetkid,length,1,icscode)
         do it=1,nelements
            itetkid(it)=0
         enddo
      endif
C
      length=nelements
      call mmgetblk('itactive',isubname,ipitactive,length,1,icscode)
      length=icharlnf(coption)
      if(coption(1:length).eq.'-all-') then
         ntactive=nelements
         do it=1,nelements
            itactive(it)=it
         enddo
      elseif(coption(1:length).eq.'-active-') then
         ntactive=0
         call get_active_elements(cmo,ntactive,ipitactive)
      endif
C
      length=nelements
      call mmgetblk('itflag',isubname,ipitflag,length,2,ierr)
      call mmgetblk('itpath',isubname,ipitpath,length,2,ierr)
C
      xavg_value=zero
      yavg_value=zero
      zavg_value=zero
      do i1=1,nump
         xavg_value=xavg_value+xic_value(i1)
         yavg_value=yavg_value+yic_value(i1)
         zavg_value=zavg_value+zic_value(i1)
      enddo
      xavg_value=xavg_value/nump
      yavg_value=yavg_value/nump
      zavg_value=zavg_value/nump
      itstart=1
      distmin=1.0d+100
      do idum=1,ntactive
         it=itactive(idum)
         xavg=zero
         yavg=zero
         zavg=zero
         do i=1,nelmnen(itettyp(it))
            i1=itet1(itetoff(it)+i)
            xavg=xavg+xic_table(i1)
            yavg=yavg+yic_table(i1)
            zavg=zavg+zic_table(i1)
         enddo
         xavg=xavg/real(nelmnen(itettyp(it)))
         yavg=yavg/real(nelmnen(itettyp(it)))
         zavg=zavg/real(nelmnen(itettyp(it)))
         dist=(xavg-xavg_value)**2+
     *        (yavg-yavg_value)**2+
     *        (zavg-zavg_value)**2
         if(dist.lt.distmin) then
            itstart=it
            distmin=dist
         endif
      enddo
C
      iprintcount=0
C
      do ip=1,nump
         xa=xic_value(ip)
         ya=yic_value(ip)
         za=zic_value(ip)
C
      itnext=itstart
      iconvex=0
 200  continue
      do it=1,nelements
         itflag(it)=0
      enddo
      ilost=0
      irestart=0
      iflag=0
      ipath=0
      dowhile(iflag.eq.0)
         if(ipath.ge.nelements) then
            iprintcount=iprintcount+1
            if(iprintcount.le.30) then
               write(logmess,'(a,a,i10)') 'To many paths: ',
     *                                    cmo(1:icharlnf(cmo)),
     *                                    ip
               call writloga('default',0,logmess,0,ierr)
            endif
            goto 110
         endif
         ipath=ipath+1
         itpath(ipath)=itnext
         itflag(itnext)=1
         in_element=-1
         if(itettyp(itnext).eq.ifelmtri) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            call inside_tri2d(x1,y1,z1,x2,y2,z2,x3,y3,z3,
     *                        xa,ya,za,
     *                        in_element)
         elseif(itettyp(itnext).eq.ifelmqud) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            i4=itet1(itetoff(itnext)+4)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            x4=xic_table(i4)
            y4=yic_table(i4)
            z4=zic_table(i4)
            call inside_quad2d(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                         xa,ya,za,
     *                         in_element)
         elseif(itettyp(itnext).eq.ifelmtet) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            i4=itet1(itetoff(itnext)+4)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            x4=xic_table(i4)
            y4=yic_table(i4)
            z4=zic_table(i4)
            call inside_tet(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                      xa,ya,za,
     *                      in_element)
         elseif(itettyp(itnext).eq.ifelmpyr) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            i4=itet1(itetoff(itnext)+4)
            i5=itet1(itetoff(itnext)+5)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            x4=xic_table(i4)
            y4=yic_table(i4)
            z4=zic_table(i4)
            x5=xic_table(i5)
            y5=yic_table(i5)
            z5=zic_table(i5)
            call inside_pyr(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                      x5,y5,z5,
     *                      xa,ya,za,
     *                      in_element)
         elseif(itettyp(itnext).eq.ifelmpri) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            i4=itet1(itetoff(itnext)+4)
            i5=itet1(itetoff(itnext)+5)
            i6=itet1(itetoff(itnext)+6)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            x4=xic_table(i4)
            y4=yic_table(i4)
            z4=zic_table(i4)
            x5=xic_table(i5)
            y5=yic_table(i5)
            z5=zic_table(i5)
            x6=xic_table(i6)
            y6=yic_table(i6)
            z6=zic_table(i6)
            call inside_pri(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                      x5,y5,z5,x6,y6,z6,
     *                      xa,ya,za,
     *                      in_element)
         elseif(itettyp(itnext).eq.ifelmhex) then
            i1=itet1(itetoff(itnext)+1)
            i2=itet1(itetoff(itnext)+2)
            i3=itet1(itetoff(itnext)+3)
            i4=itet1(itetoff(itnext)+4)
            i5=itet1(itetoff(itnext)+5)
            i6=itet1(itetoff(itnext)+6)
            i7=itet1(itetoff(itnext)+7)
            i8=itet1(itetoff(itnext)+8)
            x1=xic_table(i1)
            y1=yic_table(i1)
            z1=zic_table(i1)
            x2=xic_table(i2)
            y2=yic_table(i2)
            z2=zic_table(i2)
            x3=xic_table(i3)
            y3=yic_table(i3)
            z3=zic_table(i3)
            x4=xic_table(i4)
            y4=yic_table(i4)
            z4=zic_table(i4)
            x5=xic_table(i5)
            y5=yic_table(i5)
            z5=zic_table(i5)
            x6=xic_table(i6)
            y6=yic_table(i6)
            z6=zic_table(i6)
            x7=xic_table(i7)
            y7=yic_table(i7)
            z7=zic_table(i7)
            x8=xic_table(i8)
            y8=yic_table(i8)
            z8=zic_table(i8)
            call inside_hex(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                      x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,
     *                      xa,ya,za,
     *                      in_element)
         endif
         nen=nelmnen(itettyp(itnext))
         nef=nelmnef(itettyp(itnext))
         if(in_element.ge.0) then
            iflag=1
         else
            if(itettyp(itnext).eq.ifelmtri) then
               call volume_tri(x1,y1,z1,x2,y2,z2,x3,y3,z3,
     *                         volelm)
            elseif(itettyp(itnext).eq.ifelmqud) then
               call volume_qud(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                         volelm)
            elseif(itettyp(itnext).eq.ifelmtet) then
               call volume_tet(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                         volelm)
            elseif(itettyp(itnext).eq.ifelmpyr) then
               call volume_hex(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                         x5,y5,z5,x5,y5,z5,x5,y5,z5,x5,y5,z5,
     *                         volelm)
            elseif(itettyp(itnext).eq.ifelmpri) then
               call volume_hex(x1,y1,z1,x2,y2,z2,x2,y2,z2,x3,y3,z3,
     *                         x4,y4,z4,x5,y5,z5,x5,y5,z5,x6,y6,z6,
     *                         volelm)
            else
               call volume_hex(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
     *                         x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,
     *                         volelm)
            endif
            do i=1,nef
               if(itettyp(itnext).eq.ifelmtri .or.
     *            itettyp(itnext).eq.ifelmqud) then
                  if2=itet1(itetoff(itnext) +
     *                ielmface1(1,i,itettyp(itnext)))
                  if3=itet1(itetoff(itnext) +
     *                ielmface1(2,i,itettyp(itnext)))
                  x1=0.0
                  y1=0.0
                  z1=0.0
                  do j=1,nelmnen(itettyp(itnext))
                     if1=itet1(itetoff(itnext)+j)
                     x1=x1+xic_table(if1)
                     y1=y1+yic_table(if1)
                     z1=z1+zic_table(if1)
                  enddo
                  x1=x1/nelmnen(itettyp(itnext))
                  y1=y1/nelmnen(itettyp(itnext))
                  z1=z1/nelmnen(itettyp(itnext))
                  x2=xic_table(if2)
                  y2=yic_table(if2)
                  z2=zic_table(if2)
                  x3=xic_table(if3)
                  y3=yic_table(if3)
                  z3=zic_table(if3)
                  dx=  ((y3-y2)*(za-z2)-(ya-y2)*(z3-z2))
                  dy= -((x3-x2)*(za-z2)-(xa-x2)*(z3-z2))
                  dz=  ((x3-x2)*(ya-y2)-(xa-x2)*(y3-y2))
                  dxt= ((y3-y2)*(z1-z2)-(y1-y2)*(z3-z2))
                  dyt=-((x3-x2)*(z1-z2)-(x1-x2)*(z3-z2))
                  dzt= ((x3-x2)*(y1-y2)-(x1-x2)*(y3-y2))
                  volface(i)=-(dxt*dx+dyt*dy+dzt*dz)
               elseif(itettyp(itnext).eq.ifelmtet) then
                  if2=itet1(itetoff(itnext) +
     *                ielmface1(1,i,itettyp(itnext)))
                  if3=itet1(itetoff(itnext) +
     *                ielmface1(2,i,itettyp(itnext)))
                  if4=itet1(itetoff(itnext) +
     *                ielmface1(3,i,itettyp(itnext)))
                  x2=xic_table(if2)
                  y2=yic_table(if2)
                  z2=zic_table(if2)
                  x3=xic_table(if3)
                  y3=yic_table(if3)
                  z3=zic_table(if3)
                  x4=xic_table(if4)
                  y4=yic_table(if4)
                  z4=zic_table(if4)
                  dx= ((y2-y3)*(z4-z3)-(y4-y3)*(z2-z3))
                  dy=-((x2-x3)*(z4-z3)-(x4-x3)*(z2-z3))
                  dz= ((x2-x3)*(y4-y3)-(x4-x3)*(y2-y3))
                  volface(i)=+((x3-xa)*dx+(y3-ya)*dy+(z3-za)*dz)
               elseif(itettyp(itnext).eq.ifelmpyr) then
                  if2=itet1(itetoff(itnext) +
     *                ielmface1(1,i,itettyp(itnext)))
                  if3=itet1(itetoff(itnext) +
     *                ielmface1(2,i,itettyp(itnext)))
                  if4=itet1(itetoff(itnext) +
     *                ielmface1(3,i,itettyp(itnext)))
                  x2=xic_table(if2)
                  y2=yic_table(if2)
                  z2=zic_table(if2)
                  x3=xic_table(if3)
                  y3=yic_table(if3)
                  z3=zic_table(if3)
                  x4=xic_table(if4)
                  y4=yic_table(if4)
                  z4=zic_table(if4)
                  dx= ((y2-y3)*(z4-z3)-(y4-y3)*(z2-z3))
                  dy=-((x2-x3)*(z4-z3)-(x4-x3)*(z2-z3))
                  dz= ((x2-x3)*(y4-y3)-(x4-x3)*(y2-y3))
                  volface(i)=+((x3-xa)*dx+(y3-ya)*dy+(z3-za)*dz)
               elseif(itettyp(itnext).eq.ifelmpri) then
                  if2=itet1(itetoff(itnext) +
     *                ielmface1(1,i,itettyp(itnext)))
                  if3=itet1(itetoff(itnext) +
     *                ielmface1(2,i,itettyp(itnext)))
                  if4=itet1(itetoff(itnext) +
     *                ielmface1(3,i,itettyp(itnext)))
                  x2=xic_table(if2)
                  y2=yic_table(if2)
                  z2=zic_table(if2)
                  x3=xic_table(if3)
                  y3=yic_table(if3)
                  z3=zic_table(if3)
                  x4=xic_table(if4)
                  y4=yic_table(if4)
                  z4=zic_table(if4)
                  dx= ((y2-y3)*(z4-z3)-(y4-y3)*(z2-z3))
                  dy=-((x2-x3)*(z4-z3)-(x4-x3)*(z2-z3))
                  dz= ((x2-x3)*(y4-y3)-(x4-x3)*(y2-y3))
                  volface(i)=+((x3-xa)*dx+(y3-ya)*dy+(z3-za)*dz)
               elseif(itettyp(itnext).eq.ifelmhex) then
                  if2=itet1(itetoff(itnext) +
     *                ielmface1(1,i,itettyp(itnext)))
                  if3=itet1(itetoff(itnext) +
     *                ielmface1(2,i,itettyp(itnext)))
                  if4=itet1(itetoff(itnext) +
     *                ielmface1(3,i,itettyp(itnext)))
                  x2=xic_table(if2)
                  y2=yic_table(if2)
                  z2=zic_table(if2)
                  x3=xic_table(if3)
                  y3=yic_table(if3)
                  z3=zic_table(if3)
                  x4=xic_table(if4)
                  y4=yic_table(if4)
                  z4=zic_table(if4)
                  dx= ((y2-y3)*(z4-z3)-(y4-y3)*(z2-z3))
                  dy=-((x2-x3)*(z4-z3)-(x4-x3)*(z2-z3))
                  dz= ((x2-x3)*(y4-y3)-(x4-x3)*(y2-y3))
                  volface(i)=+((x3-xa)*dx+(y3-ya)*dy+(z3-za)*dz)
               endif
            enddo
            do j=1,nef
               index=ismax(nef,volface,1)
               if(volface(index).ge.-xsmall*max(xsmall,volelm)) then
                  if(jtet1(jtetoff(itnext)+index).eq.mbndry) then
                     iflag=2
                  elseif(jtet1(jtetoff(itnext)+index).gt.mbndry) then
                     ittest=1+
     *                    (jtet1(jtetoff(itnext)+index)-mbndry-1)/nefcmo
                     if(itflag(ittest).eq.0) then
                        dowhile(itetkid(ittest).ne.0)
                           itflag(ittest)=1
                           ittest=itetkid(ittest)
                        enddo
                        itnext=ittest
                        goto 100
                     else
                        volface(index)=-1.0d+30
                     endif
                  else
                     ittest=1+(jtet1(jtetoff(itnext)+index)-1)/nefcmo
                     if(itflag(ittest).eq.0) then
                        dowhile(itetkid(ittest).ne.0)
                           itflag(ittest)=1
                           ittest=itetkid(ittest)
                        enddo
                        itnext=ittest
                        goto 100
                     else
                        volface(index)=-1.0d+30
                     endif
                  endif
               endif
            enddo
         endif
         if(iflag.eq.0.and.itflag(itnext).ne.0) then
            do i=1,nef
               if(jtet1(jtetoff(itnext)+i).eq.mbndry) then
               elseif(jtet1(jtetoff(itnext)+i).gt.mbndry) then
                  ittest=1+
     *                   (jtet1(jtetoff(itnext)+i)-mbndry-1)/nefcmo
                  if(itflag(ittest).eq.0) then
                     dowhile(itetkid(ittest).ne.0)
                        itflag(ittest)=1
                        ittest=itetkid(ittest)
                     enddo
                     itnext=ittest
                     goto 100
                  endif
               else
                  ittest=1+(jtet1(jtetoff(itnext)+i)-1)/nefcmo
                  if(itflag(ittest).eq.0) then
                     dowhile(itetkid(ittest).ne.0)
                        itflag(ittest)=1
                        ittest=itetkid(ittest)
                     enddo
                     itnext=ittest
                     goto 100
                  endif
               endif
            enddo
C*****      iflag=3
            ilost=ilost+1
            irestart=itnext
            ipath=ipath+1
            itpath(ipath)=-irestart
            do it=1,nelements
               itflag(it)=0
            enddo
         endif
 100     continue
      enddo
C
 110  continue
C
      if(iclrpath.gt.0) then
         do i=1,ipath
            itetclr(abs(itpath(i)))=iclrpath
         enddo
      endif
C
      if(iflag.eq.0) then
         ielement(ip)=0
      elseif(iflag.eq.2) then
         if(iconvex.eq.0) then
            distmin=1.0d+100
            do idum=1,ntactive
               it=itactive(idum)
               xavg=zero
               yavg=zero
               zavg=zero
               do i=1,nelmnen(itettyp(it))
                  i1=itet1(itetoff(it)+i)
                  xavg=xavg+xic_table(i1)
                  yavg=yavg+yic_table(i1)
                  zavg=zavg+zic_table(i1)
               enddo
               xavg=xavg/real(nelmnen(itettyp(it)))
               yavg=yavg/real(nelmnen(itettyp(it)))
               zavg=zavg/real(nelmnen(itettyp(it)))
               dist=(xavg-xa)**2+
     *              (yavg-ya)**2+
     *              (zavg-za)**2
               if(dist.lt.distmin) then
                  itnext=it
                  distmin=dist
               endif
            enddo
            iconvex=iconvex+1
            goto 200
         else
            distmin=1.0d+100
            do i=1,nnodes
               if(ireal1(i).eq.1) then
                  dist=(xic_table(i)-xa)**2+
     *                 (yic_table(i)-ya)**2+
     *                 (zic_table(i)-za)**2
                  if(dist.lt.distmin) then
                     iclose=i
                     distmin=dist
                  endif
               endif
            enddo
            ielement(ip)=-iclose
         endif
      else
         ielement(ip)=itnext
         itstart=itnext
      endif
C
      enddo
C
      if(iprintcount.gt.0) then
         write(logmess,'(a,i10)') 'Total lost paths: ',iprintcount
         call writloga('default',1,logmess,1,ierr)
      endif
C
      goto 9999
 9999 continue
      call mmrelprt(isubname,icscode)
      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
Founder and Davis resident Matt Blake is forming the Afghan Pepper Company (http://www.afghanpeppercompany.com/ website here) to give Afghani farmers a chance to grow a legal trade crop, hot sauce, which will be processed in Afghanistan and sold around the world.  Here is what he hopes will happen: Farmers grow produce other than opium poppy, demand for barrels and bottles and safe transportation rises, food processing expands, jobs are created, stability ensues, governments are supported by the governed, troops come home. All that and the world has a new, very delicious hot sauce to enjoy.  This http://www.davisenterprise.com/localnews/changingtheworldwithhotsauce/ Davis Enterprise article has further details.

Blake sought a total of $25,000 of http://www.kickstarter.com/projects/matthewpblake/theafghanpeppercompany Kickstarter money for his project and was able to secure enough funding by the deadline (Feb 1, 2:59am EST).  The total amount pledged was $27,344.  As those of you who contributed money to LocalWiki know, Kickstarter funding is an allornothing prospect.

The Afghan Pepper Company was also featured on http://gooddaysacramento.cbslocal.com/video/6686049afghanpeppercompany/ Good Day Sacramento.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 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 viewfactors(textpart,iviewfile,istep,inpc,
     &  istat,n,key,iline,ipol,inl,ipoinp,inp,jobnamec,ipoinpc,ier)
!
!     reading the input deck: *VIEWFACTOR
!
      implicit none
!
      character*1 inpc(*)
      character*132 textpart(16),jobnamec(*)
!
      integer i,iviewfile,istep,n,istat,iline,ipol,inl,ipoinp(2,*),
     &  inp(3,*),key,j,k,l,ipoinpc(0:*),ier
!
      if(istep.lt.1) then
         write(*,*) '*ERROR reading *VIEWFACTOR: *VIEWFACTOR can '
         write(*,*) '       only be used within a STEP'
         ier=1
         return
      endif
!
      do i=2,n
         if(textpart(i)(1:4).eq.'READ') then
            if(iviewfile.eq.0) then
               iviewfile=-1
            elseif(iviewfile.gt.0) then
               write(*,*) '*ERROR reading *VIEWFACTOR: READ and WRITE/'
               write(*,*) '       WRITE ONLY are mutually exclusive'
               call inputerror(inpc,ipoinpc,iline,
     &              "*VIEWFACTOR%",ier)
               return
            endif
         elseif(textpart(i)(1:8).eq.'NOCHANGE') then
            if(istep.eq.1) then
               write(*,*) '*ERROR reading *VIEWFACTOR: NO CHANGE cannot'
               write(*,*) '       be used in the first step'
               call inputwarning(inpc,ipoinpc,iline,
     &"*VIEWFACTOR%")
            elseif(iviewfile.le.0) then
               iviewfile=-2
            elseif(iviewfile.gt.0) then
               write(*,*) '*ERROR reading *VIEWFACTOR: NO CHANGE and'
               write(*,*) '       WRITE/WRITE ONLY are mutually'
               write(*,*) '       exclusive'
               call inputerror(inpc,ipoinpc,iline,
     &              "*VIEWFACTOR%",ier)
               return
            endif
         elseif(textpart(i)(1:9).eq.'WRITEONLY') then
            if(iviewfile.eq.0) then
               iviewfile=3
            elseif(iviewfile.lt.0) then
               write(*,*) '*ERROR reading *VIEWFACTOR: '
               write(*,*) '       WRITE ONLY and READ/NO CHANGE'
               write(*,*) '       are mutually exclusive'
               call inputerror(inpc,ipoinpc,iline,
     &              "*VIEWFACTOR%",ier)
               return
            endif
         elseif(textpart(i)(1:5).eq.'WRITE') then
            if(iviewfile.eq.0) then
               iviewfile=2
            elseif(iviewfile.lt.0) then
               write(*,*) '*ERROR reading *VIEWFACTOR: WRITE'
               write(*,*) '       and READ/NO CHANGE'
               write(*,*) '       are mutually exclusive'
               call inputerror(inpc,ipoinpc,iline,
     &              "*VIEWFACTOR%",ier)
               return
            endif
         elseif(textpart(i)(1:6).eq.'INPUT=') then
            jobnamec(2)(1:126)=textpart(i)(7:132)
            jobnamec(2)(127:132)='      '
            loop1: do j=1,126
               if(jobnamec(2)(j:j).eq.'"') then
                  do k=j+1,126
                     if(jobnamec(2)(k:k).eq.'"') then
                        do l=k-1,126
                           jobnamec(2)(l:l)=' '
                           exit loop1
                        enddo
                     endif
                     jobnamec(2)(k-1:k-1)=jobnamec(2)(k:k)
                  enddo
                  jobnamec(2)(126:126)=' '
               endif
            enddo loop1
         elseif(textpart(i)(1:7).eq.'OUTPUT=') then
            jobnamec(3)(1:125)=textpart(i)(8:132)
            jobnamec(3)(126:132)='      '
            loop2: do j=1,125
               if(jobnamec(3)(j:j).eq.'"') then
                  do k=j+1,125
                     if(jobnamec(3)(k:k).eq.'"') then
                        do l=k-1,125
                           jobnamec(3)(l:l)=' '
                           exit loop2
                        enddo
                     endif
                     jobnamec(3)(k-1:k-1)=jobnamec(3)(k:k)
                  enddo
                  jobnamec(3)(125:125)=' '
               endif
            enddo loop2
         else
            write(*,*)
     &        '*WARNING reading *VIEWFACTOR: parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*VIEWFACTOR%")
         endif
      enddo
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
// @author David Han
//
// @kind function
//
// @fileOverview Create a basic line chart.
//
// @example Create and render a line chart specification:
// .qp.go[600;400;] .line.basic.f03[]
//
// @returns {table} a GG specification (for example, ``.qp.point[t; `x; `y; ::]``)

{[]
  xfmt:{n:`long$last string x;$[0=n mod 2;x;`]};
  yfmt:{`$string[0.1*floor 0.5+1000*x],"%"};

  .qp.line[.line.iexMktShare;`quarter;`mktshr;]
     .qp.s.labels[`x`y!(`Quarter;`$"Market Share (%)")]
    ,.qp.s.scale[`x;.gg.scale.format[xfmt;.gg.scale.categorical[]]]
    ,.qp.s.scale[`y;.gg.scale.format[yfmt;] .gg.scale.linear]
  }

max_stars_repo_name: dpkwhan/qplot

max_issues_repo_name: dpkwhan/qplot

max_forks_repo_name: dpkwhan/qplot


text:
      subroutine zonevols
c***********************************************************************
c     calculates absolute cell volumes (cc) needed for plankian source.                           *
c***********************************************************************
c     implicit double precision (a-h,o-z)
      include 'params.inc'
      include 'geomz.inc'

      do 20 j = 1,nzones
         g_volcl(j) = 0.0d0
 20   continue

      do 21 ks = 1,4
         do 22 j = 1,nzones
            if (ng_itype(j,ks) .ne. 2 .and. ng_itype(j,ks) .ne. 4) then
               tv1 = (g_zz(j,ks+1) -
     .              g_zz(j,ks))*(g_rr(j,ks+1)*(g_rr(j,ks+1)
     .              + g_rr(j,ks)) + g_rr(j,ks)*g_rr(j,ks))
               g_volcl(j) = g_volcl(j) + 1.0471976*tv1
            endif
 22      continue
 21   continue

      return
      end

max_stars_repo_name: gmlueck/llvm-test-suite

max_issues_repo_name: gmlueck/llvm-test-suite

max_forks_repo_name: gmlueck/llvm-test-suite


text:
module forno1d

contains

subroutine calcula_perfil(Ta, T1, T2, T3, R1, R2, R3,T)
    implicit none

    real, intent(in) :: Ta, T1, T2, T3, R1, R2, R3
    integer, parameter :: n = 30   ! número de elementos finitos
    real, dimension(n), intent(out) :: T

    real :: L, sor, dx, T_media, Tn
    integer :: p_n1, p_n2, p_n3, i, j ,k
    real, dimension(n) :: fonte

    L = 2.0         ! comprimento do forno em metros
    sor = 1.85      ! parâmetro sor
    p_n1 = n/4      ! posição do primeiro sensor de Temperatura
    p_n2 = n/2      ! segundo sensor
    p_n3 = 3*n/4    ! terceiro sensor

    dx = L/(n+1)
    T_media = (2*Ta + T1 + T2 + T3)/5

    ! valores iniciais da temperatura
    do i=1, n-1
        T(i) = T_media
    enddo
    T(1) = Ta   ! tempa ambiente nas pontas (direichilet)
    T(n) = Ta   ! tempa ambiente nas pontas (direichilet)

    ! valores da fonte
    do i=1, INT(n/3)
        fonte(i) = R1
    enddo
    do i=INT(n/3)+1, INT(2*n/3)
        fonte(i) = R2
    enddo
    do i=INT(2*n/3)+1, n
        fonte(i) = R3
    enddo

    do j=1, 200
        do i=2 , n-1
            if (i == p_n1) then
                T(i) = T1
            else if (i == p_n2) then
                T(i) = T2
            else if (i == p_n3) then
                T(i) = T3
            else
                Tn = (T(i-1) + T(i+1))/2.0 +dx*dx*fonte(i)
                T(i) = T(i) + sor*( Tn - T(i) )
            end if
        enddo
    enddo

end subroutine calcula_perfil

end module forno1d

max_stars_repo_name: Atzingen/controleForno-interface

max_issues_repo_name: Atzingen/controleForno-interface

max_forks_repo_name: Atzingen/controleForno-interface


text:
module projection
      !Generate csf by projection
      use prep, only: i16b, ARRAY_START_LENGTH,ARRAY_LONG_LENGTH, ARRAY_SHORT_LENGTH, DET_MAX_LENGTH, equals0, REAL_MIN
      use, intrinsic :: iso_fortran_env, only: rk => real64
      use detgen, only: lpls, lms, eposinit, spls, sms, Szt2_det, Lz_det

      implicit none

      contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!################## OPERATORS ##################
! subroutine scalar_append(const, inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)
! subroutine Lplus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)

! subroutine Lminus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)

!subroutine Splus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)
! subroutine Sminus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)
!
! subroutine LzLzplus1_append(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)
!subroutine SzSzplus1_append(inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)

! subroutine Proj_S(Smin2, Smax2, Sdes2, &
!              &inbasis, incoefs, inepos, iniszeros, innum, &
!               & basislist, coeflist, eposlist, iszerolist, num)
! subroutine Proj_L(Lmin, Lmax, Ldes, &
!               & inbasis, incoefs, inepos, iniszeros, innum,  &
!              & basislist, coeflist, eposlist, iszerolist, num)

!################ OPERATOR HELPERS #################
!!!!!!!! Operators with single det pair as inputs !!!!!!!!!!!!!!!!!!!!
! subroutine Sminus_single(detup, detdn, eposup, eposdn, coef, iszero, &
!               & basislist, coeflist, eposlist, iszerolist, num)
! subroutine Splus_single(detup, detdn, eposup, eposdn, coef, iszero, &
!               & basislist, coeflist, eposlist, iszerolist, num)
! subroutine Lminus_single(detup, detdn, eposup, eposdn, coef, iszero, &
!               & basislist, coeflist, eposlist, iszerolist, num )
! subroutine Lplus_single(detup, detdn, eposup, eposdn, coef, iszero, &
!               & basislist, coeflist, eposlist, iszerolist, num )

!################# CSF TABLE SUBROUTINES ##################
! subroutine normalizetable(coeftable, ndets, ncsf)

! subroutine getallsigns(basislist, coeflist, eposlist, iszerolist, num)
!               getsign(detup, detdn, eposup, eposdn, coef)
! subroutine collect_csf(basislist, coeflist, iszerolist, num,&
!               & allbasis, coeftable, ndets, ncsf)

!################## OTHER SMALL HELPERS #################
! initlists(det, basis, coefs, eposes, iszeros, num) init from single det pair to basislist and other lists
!TODO  initlists_fromlists(detlist, incoeflist,  innum,coeflist, eposlist, iszerolist) not tested
! integer function  countswps(array, num)
! integer function findindetlist(det, allbasis, ndets)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !Moved to prep.f90
      !logical function equals0(coef)
      !    real(rk), intent(in) :: coef
      !    if(abs(coef).le.REAL_MIN) then
      !        equals0 = .true.
      !    else
      !        equals0 = .false.
      !    endif
      !end function equals0

      subroutine initlists(det, basis, coefs, eposes, iszeros, num)
          !init from single det pair
          integer(i16b), intent(in) :: det(:)
          integer(i16b), allocatable :: basis(:, :)
          real(rk), allocatable :: coefs(:)
          integer, allocatable :: eposes(:, :, :)
          logical, allocatable :: iszeros(:)
          integer ::  num
          if(.not.allocated(basis)) then
              allocate(basis(ARRAY_SHORT_LENGTH, 2))
          endif
          if(.not.allocated(coefs)) then
              allocate(coefs(ARRAY_SHORT_LENGTH))
          endif
          if(.not.allocated(eposes)) then
              allocate(eposes(ARRAY_SHORT_LENGTH, 2, DET_MAX_LENGTH))
          endif
          if(.not.allocated(iszeros)) then
              allocate(iszeros(ARRAY_SHORT_LENGTH))
          endif
          basis(1, 1:2) = det(1:2)
          call eposinit(eposes(1, 1, :), eposes(1, 2, :), det(1), det(2))
          coefs(1) = 1.0
          iszeros(1) = .false.
          num = 1
      end subroutine initlists

      subroutine initlists_fromlists(detlist, incoeflist, innum, coeflist, eposlist, iszerolist)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          !!
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          integer(i16b), intent(in) :: detlist(:, :)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          integer :: i
          real(rk), intent(in) :: incoeflist(:)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)

          if(.not.allocated(basislist)) then
              allocate(basislist(ARRAY_SHORT_LENGTH, 2))
          endif
          if(.not.allocated(coeflist)) then
              allocate(coeflist(ARRAY_SHORT_LENGTH))
          endif
          if(.not.allocated(eposlist)) then
              allocate(eposlist(ARRAY_SHORT_LENGTH, 2, DET_MAX_LENGTH))
          endif
          if(.not.allocated(iszerolist)) then
              allocate(iszerolist(ARRAY_SHORT_LENGTH))
          endif
          do i = 1, innum
              call eposinit(eposlist(i, 1, :), eposlist(i, 2, :), detlist(i, 1), detlist(i, 2))
              coeflist(i) = incoeflist(i)

              iszerolist(i) = .false.
          enddo
      end subroutine initlists_fromlists

      subroutine Proj_S(Smin2, Smax2, Sdes2, &
              &inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk), intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: Smax2, Smin2, Sdes2, innum
          integer(i16b), allocatable :: basislist(:, :), tpbasis(:, :), basislist1(:, :)
          real(rk), allocatable :: coeflist(:), tpcoefs(:), coeflist1(:)
          integer, allocatable :: eposlist(:, :, :), tpepos(:, :, :), eposlist1(:, :, :)
          logical, allocatable :: iszerolist(:), tpiszeros(:), iszerolist1(:)
          integer :: num, tpnum, i, Sp2, num1
          real(rk) :: tail
          ! write(*, *) '>>>>>>>>>>>>>>>>>>>>>>> Projection P_S '
          allocate(tpbasis(ARRAY_LONG_LENGTH, 2))
          allocate(tpcoefs(ARRAY_LONG_LENGTH))
          allocate(tpepos(ARRAY_LONG_LENGTH, 2, DET_MAX_LENGTH))
          allocate(tpiszeros(ARRAY_LONG_LENGTH))
          !----------Initialize tp.. as in..
          do i = 1, innum
              tpbasis(i, 1:2) = inbasis(i, 1:2)
              tpcoefs(i) = incoefs(i)
              tpepos(i, 1, :) = inepos(i, 1, :)
              tpepos(i, 2, :) = inepos(i, 2, :)
              tpiszeros(i) = iniszeros(i)
          enddo
          tpnum = innum
          !---------- Apply operators for each Sp2 -----------
          do Sp2 = Smin2, Smax2
              if(Sp2.eq.Sdes2) then
                  cycle
              endif

              !S^2 = S-S+.....( + Sz^2 + Sz later)
              ! write(*, *) ' '
              ! write(*,  '("====================== Now at 2Sp =", 1I3)') Sp2
              !-------set intermediates (basislist1...) and output ( basislist ... ) empty
              num1 = 0
              num = 0
              !----------S^2 = S-S+ + Sz^2 + Sz
              call Splus_multiple(tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist1, coeflist1, eposlist1, iszerolist1, num1)
              ! write(*, *) ' '
              ! write(*, *) 'Calling Sminus'
              call Sminus_multiple(basislist1, coeflist1, eposlist1, iszerolist1, num1, &
                  & basislist, coeflist, eposlist, iszerolist, num)
              call SzSzplus1_append(tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist, coeflist, eposlist, iszerolist, num)

              !-------------------(- Sp(Sp+1))---------------------
              ! write(*, *) ' ------------------- Append tail to basislist '
              tail = real(-Sp2/2.0 * (Sp2/2.0 + 1.0), rk)
              ! write(*, '("tail = ", 1f15.10)') tail
              call scalar_append(tail, tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist, coeflist, eposlist, iszerolist, num)
              ! write(*, '("Total number of dets in basislist", 1I10)') num

              !----------------------------tp.. <- output ( basislist, ...)
              tpnum = 0
              do i = 1, num
                  tpbasis(i, 1:2) = basislist(i, 1:2)
                  tpcoefs(i) = coeflist(i)
                  tpepos(i, 1, :) = eposlist(i, 1, :)
                  tpepos(i, 2, :) = eposlist(i, 2, :)
                  tpiszeros(i) = iszerolist(i)
              enddo
              tpnum = num
          enddo
          deallocate(tpbasis)
          deallocate(tpcoefs)
          deallocate(tpepos)
          deallocate(tpiszeros)
      end subroutine Proj_S

      subroutine Proj_L(Lmin, Lmax, Ldes, &
              & inbasis, incoefs, inepos, iniszeros, innum,  &
              & basislist, coeflist, eposlist, iszerolist, num)
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk), intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: Lmax, Lmin, Ldes, innum
          integer(i16b), allocatable :: basislist(:, :), tpbasis(:, :), basislist1(:, :)
          real(rk), allocatable :: coeflist(:), tpcoefs(:), coeflist1(:)
          integer, allocatable :: eposlist(:, :, :), tpepos(:, :, :), eposlist1(:, :, :)
          logical, allocatable :: iszerolist(:), tpiszeros(:), iszerolist1(:)
          integer :: num, tpnum, i, Lp, num1
          real(rk) :: tail
          ! write(*, *) '>>>>>>>>>>>>>>>>>>>>>>> Projection P_L '
          allocate(tpbasis(ARRAY_START_LENGTH, 2))
          allocate(tpcoefs(ARRAY_START_LENGTH))
          allocate(tpepos(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
          allocate(tpiszeros(ARRAY_START_LENGTH))
          !----------Initialize tp.. as in..
          do i = 1, innum
              tpbasis(i, 1:2) = inbasis(i, 1:2)
              tpcoefs(i) = incoefs(i)
              tpepos(i, 1, :) = inepos(i, 1, :)
              tpepos(i, 2, :) = inepos(i, 2, :)
              tpiszeros(i) = iniszeros(i)
          enddo
          tpnum = innum
          !---------- Apply operators for each Lp -----------
          do Lp = Lmin, Lmax
              if(Lp.eq.Ldes) then
                  cycle
              endif

              !L^2 = L-L+ + Lz^2 + Lz
              ! write(*, *) ' '
              ! write(*,  '("====================== Now at Lp =", 1I3)') Lp
              !-------set intermediates (basislist1...) and output ( basislist ... ) empty
              num1 = 0
              num = 0
              !----------L^2 = L-L+ + Lz^2 + Lz
              call Lplus_multiple(tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist1, coeflist1, eposlist1, iszerolist1, num1)
              call Lminus_multiple(basislist1, coeflist1, eposlist1, iszerolist1, num1, &
                  & basislist, coeflist, eposlist, iszerolist, num)
              call LzLzplus1_append(tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist, coeflist, eposlist, iszerolist, num)
              !-------------------(- Lp(Lp+1))---------------------
              ! write(*, *) ' ------------------- Append tail to basislist '
              tail = real(-Lp * (Lp + 1), rk)
              ! write(*, '("tail = ", 1f15.10)') tail
              call scalar_append(tail, tpbasis, tpcoefs, tpepos, tpiszeros, tpnum, &
                  & basislist, coeflist, eposlist, iszerolist, num)
              ! write(*, '("Total number of dets in basislist", 1I4)') num

              !----------------------------tp.. <- output ( basislist, ...)
              tpnum = 0
              do i = 1, num
                  tpbasis(i, 1:2) = basislist(i, 1:2)
                  tpcoefs(i) = coeflist(i)
                  tpepos(i, 1, :) = eposlist(i, 1, :)
                  tpepos(i, 2, :) = eposlist(i, 2, :)
                  tpiszeros(i) = iszerolist(i)
              enddo
              tpnum = num
          enddo
          deallocate(tpbasis)
          deallocate(tpcoefs)
          deallocate(tpepos)
          deallocate(tpiszeros)

      end subroutine Proj_L

      subroutine normalizetable(coeftable, ndets, ncsf)
          !call this after coeftable is complete
          real(rk) :: coeftable(:, :), sumcol
          integer, intent(in) :: ndets, ncsf
          integer :: i, j
          ! write(*, *) '-------- Normalize coef table ------------'
          do i = 1, ncsf
            sumcol = 0
            do j = 1, ndets
              sumcol = sumcol + coeftable(j, i)**2
            enddo
            sumcol = sqrt(sumcol)
            do j = 1, ndets
              coeftable(j, i) = coeftable(j, i)/sumcol
            enddo
          enddo
      end subroutine normalizetable

      subroutine getallsigns(basislist, coeflist, eposlist, iszerolist, num)
          integer(i16b), intent(in) :: basislist(:, :)
          real(rk) :: coeflist(:)
          integer, intent(in) :: eposlist(:, :, :), num
          logical, intent(in) :: iszerolist(:)
          integer :: i

          ! write(*, *) '>>>>>>>>>>>>>>>>>>>> Get all signs of basislist <<<<<<<<<<<<<<<<<<<'
          do i = 1, num
              if(.not.iszerolist(i)) then
                  call getsign(basislist(i, 1), basislist(i, 2), eposlist(i, 1, :),&
                      & eposlist(i, 2, :), coeflist(i))
              endif
          enddo
      end subroutine getallsigns

      subroutine getsign(detup, detdn, eposup, eposdn, coef)
          integer(i16b), intent(in) :: detup, detdn
          integer, intent(in) :: eposup(:), eposdn(:)
          integer, allocatable :: tparray(:)
          real(rk) :: coef
          integer(i16b) :: tpdet
          integer :: num, i, swaps
          ! write(*, *) 'Getsign()'
          ! write(*, '(2b16)') detup, detdn
          ! write(*, '(15I3)') eposup(1:15)
          ! write(*, '(15I3)') eposdn(1:15)

          allocate(tparray(2*size(eposup)))
          tpdet = detup
          num = 0
          do while (tpdet.ne.0)
              i = trailz(tpdet)
              num = num + 1
              tparray(num) = eposup(i+1)
              tpdet = ibclr(tpdet, i)
          enddo
          tpdet = detdn
          do while (tpdet.ne.0)
              i = trailz(tpdet)
              num = num + 1
              tparray(num) = eposdn(i+1)
              tpdet = ibclr(tpdet, i)
          enddo
          swaps = countswps(tparray, num)
          i = mod(swaps, 2)
          if(i.ne.0) then
              coef = - coef
          endif
          ! write(*, '("Sign:", 1I4)') i
          deallocate(tparray)
      end subroutine getsign

      integer function  countswps(array, num)
            integer, allocatable:: array(:)
            integer:: num, i, tp
            countswps = 0
            ! write(*, '("Start array:")')
            ! write(*, *) array(1:num)
            do i = 1, num
                do while (i.ne.array(i))
                    tp = array(i)
                    array(i) = array(tp)
                    array(tp) = tp
                    countswps = countswps + 1
                enddo
            enddo
            ! write(*, '("End array:, number of swps")')
            ! write(*, *) array(1:num)
            ! write(*, '("Count swps = ", 1I4)') countswps
      end function countswps

      subroutine collect_csf(basislist, coeflist, iszerolist, num,&
              & allbasis, coeftable, ndets, ncsf)
          !add a csf to coeftable
          !num = dimension of the basislist
          !ndets, ncsf are coeftable sizes
          !coeflist after calling getsign
          integer(i16b), intent(in):: basislist(:, :)
          real(rk), intent(in) :: coeflist(:)
          logical, intent(in) :: iszerolist(:)
          integer, intent(in) :: num
          integer(i16b), allocatable :: allbasis(:, :)
          real(rk), allocatable :: coeftable(:, :)!ndets, ncsf
          integer :: ndets, ncsf, i,j, tp
          ! write(*, *) '>>>>>>>>>>>>>>>>>>>>>>>>> Collect one csf <<<<<<<<<<<<<<<<<<<<<<<<'

          !Firsty entry initialize
          if(ndets.eq.0 .OR. ncsf.eq.0) then
              if(.not.allocated(allbasis)) then
                  allocate(allbasis(ARRAY_START_LENGTH, 2))
                  ! write(*, '("Allocated allbasis")')
              endif
              if(.not.allocated(coeftable)) then
                  allocate(coeftable(ARRAY_START_LENGTH, ARRAY_SHORT_LENGTH))
                  ! write(*, '("Allocated coeftable")')
              endif
          endif
          !Initialize new column to be 0, ndets 0s
          if(ndets.ne. 0) then
              do i = 1, ndets
                coeftable(i, ncsf+1) = .0
              enddo
              ! write(*, '("Initialized new column to 0")')
          endif

          do i = 1, num
              ! write(*, '("Processing det# from list:", 1I3)') i
              if(iszerolist(i)) then
                  cycle
              endif

              if(equals0(coeflist(i))) then
                  cycle
              endif

              tp = findindetlist(basislist(i, :), allbasis(1:ndets, :), ndets)
              ! write(*, '("findindetlist()=", 1I4)') tp
              !If not in the list, add new det
              !! now this is only appending, may need to order dets later
              if(tp.eq.0) then
                  allbasis(ndets + 1, :) = basislist(i, :)
                  do j = 1, ncsf
                    coeftable(ndets + 1, j) = .0
                  enddo
                  ! write(*, *) 'Initialized new row'
                  coeftable(ndets +1, ncsf+1) = coeflist(i)
                  ndets = ndets + 1
                  ! write(*, '("Added new det row", 1I5)') ndets
              else!If in the list, add coef
                  coeftable(tp, ncsf+1) = coeftable(tp, ncsf+1) + coeflist(i)
                  ! write(*, '("Combined coef with existing det")')
              endif
          enddo
          ncsf = ncsf + 1
          ! write(*, '("--------------- Collected csf #:", 1I3)') ncsf
          ! write(*, '("Basis:      detup         detdn     coefs of new csf")')
          do i = 1, ndets
             ! write(*, '(2b16, 1f18.10)') allbasis(i, 1:2), coeftable(i, ncsf)
          enddo
          !--------------------------------------
          !if all coefs in the coef list are zero
          !ncsf will not be incremented, so even if the new column is initialized
          !it will still be overwritten by the next loop, or accessible if no next loop
      end subroutine collect_csf

      integer function findindetlist(det, allbasis, ndets)
          integer(i16b), intent(in) :: det(:), allbasis(:, :)
          integer, intent(in) :: ndets
          integer :: i
          !! write(*, *) 'findindetlist(), ndets in all basis = ', ndets
          !! write(*, '(2B16)') det(1:2)
          !! write(*, '(2B16)') allbasis(1, 1:2)
          do i = 1, ndets
            if(det(1).eq.allbasis(i, 1)) then
                if(det(2).eq.allbasis(i, 2)) then
                    findindetlist = i
                    return
                endif
            endif
          enddo
          findindetlist = 0
      end function findindetlist

      subroutine scalar_append(const, inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk),  intent(in) :: incoefs(:), const
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i
          if(equals0(const)) then
              ! write(*, *) '0.0, exit append'
              return
          endif
          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif

          do i = 1, innum
            ! write(*, '("At input basis det #", 1I4)') i
            if(iniszeros(i)) then
                ! write(*, '("iszero, onto next det")')
                cycle
            endif
            !combine if same basis
            !tp = findindetlist(inbasis(i, 1:2), basislist(1:num, :), num)
            !! write(*, '("find in inbasis = ", 1I4)') tp
            !if(tp.ne.0) then
            !    !TODO check epos equal?
            !    coeflist(tp) = incoefs(tp) * const
            !    ! write(*, *) 'det in basislist, changed coef'
            !    cycle
            !endif
            num = num + 1
            basislist(num, 1:2) = inbasis(i, 1:2)
            coeflist(num) = incoefs(i)*const
            eposlist(num, 1, :) = inepos(i, 1, :)
            eposlist(num, 2, :) = inepos(i, 2, :)
            iszerolist(num) = iszerolist(i)
            ! write(*, '("Appended det #", 1I4)') num
          end do
      end subroutine scalar_append

      subroutine LzLzplus1_append(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk),  intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i
          real(rk) :: tpval

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          do i = 1, innum
            if(iniszeros(i)) then
                cycle
            endif
            ! write(*, *)
            ! write(*, '("=================== Applying Lz(Lz+1) to multiple dets, now at No.", 1I5)') i
            tpval = real(Lz_det(inbasis(i, 1:2)), rk) !Lz
            tpval = tpval * (tpval + 1.0)  !Lz(Lz + 1)
            if(equals0(tpval)) then
                cycle
            else
                num = num + 1
                basislist(num, 1:2) = inbasis(i, 1:2)
                coeflist(num) = incoefs(i) * tpval
                eposlist(num, 1, :) = inepos(i, 1, :)
                eposlist(num, 2, :) = inepos(i, 2, :)
                iszerolist(num) = .false.
                ! write(*, '(" Added Lz(Lz+1) = ", 1F10.5)') tpval
            endif
          end do
          ! write(*, '("Final size of output = ", 1I8)') num
      end subroutine LzLzplus1_append

      subroutine SzSzplus1_append(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk),  intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i
          real(rk) :: tpval

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_LONG_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_LONG_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_LONG_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_LONG_LENGTH))
              endif
          endif
          do i = 1, innum
            if(iniszeros(i)) then
                cycle
            endif
            ! write(*, *)
            ! write(*, '("=================== Applying Sz(Sz+1) to multiple dets, now at No.", 1I5)') i
            tpval = real(Szt2_det(inbasis(i, 1:2)), rk)/2.0 !Sz
            tpval = tpval * (tpval + 1.0)  !Sz(Sz + 1)
            if(equals0(tpval)) then
                cycle
            else
                num = num + 1
                basislist(num, 1:2) = inbasis(i, 1:2)
                coeflist(num) = incoefs(i) * tpval
                eposlist(num, 1, :) = inepos(i, 1, :)
                eposlist(num, 2, :) = inepos(i, 2, :)
                iszerolist(num) = .false.
                ! write(*, '(" Added Sz(Sz+1) = ", 1F10.5)') tpval
            endif
          end do
          ! write(*, '("Final size of output = ", 1I8)') num
      end subroutine SzSzplus1_append

      subroutine Sminus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk),  intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i
          ! write(*, *) 'subroutine Sminus_multiple'

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_LONG_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_LONG_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_LONG_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_LONG_LENGTH))
              endif
          endif
          ! write(*, '("innum = ", 1I8)') innum
          do i = 1, innum
            ! write(*, *)
            ! write(*, '("=================== Applying S- ta multiple dets, now at No.", 1I5)') i
            call Sminus_single(inbasis(i,1), inbasis(i, 2), inepos(i, 1, :), &
                &inepos(i, 2, :), incoefs(i), iniszeros(i), &
                & basislist, coeflist, eposlist, iszerolist, num)
          enddo
      end subroutine Sminus_multiple

      subroutine Splus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk), intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_LONG_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_LONG_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_LONG_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_LONG_LENGTH))
              endif
          endif
          do i = 1, innum
            ! write(*, *)
            ! write(*, '("--------------- Applying S+ to multiple dets, now at No.", 1I5)') i
            call Splus_single(inbasis(i,1), inbasis(i, 2), inepos(i, 1, :), &
                &inepos(i, 2, :), incoefs(i), iniszeros(i), &
                & basislist, coeflist, eposlist, iszerolist, num)
          enddo
      end subroutine Splus_multiple

      subroutine Lminus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to output list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk), intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          do i = 1, innum
            ! write(*, *)
            ! write(*, '("---------------------- Applying L- to multiple dets, now at No.", 1I5)') i
            call Lminus_single(inbasis(i,1), inbasis(i, 2), inepos(i, 1, :), &
                &inepos(i, 2, :), incoefs(i), iniszeros(i), &
                & basislist, coeflist, eposlist, iszerolist, num)
          enddo
      end subroutine Lminus_multiple

      subroutine Lplus_multiple(inbasis, incoefs, inepos, iniszeros, innum, &
              & basislist, coeflist, eposlist, iszerolist, num)
          !append to outaut list
          integer(i16b), intent(in) :: inbasis(:, :)
          real(rk),  intent(in) :: incoefs(:)
          integer, intent(in) :: inepos(:, :, :)
          logical, intent(in) :: iniszeros(:)
          integer, intent(in) :: innum
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)
          logical, allocatable :: iszerolist(:)
          integer :: num, i

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          do i = 1, innum
            ! write(*, *)
            ! write(*, '("-------------------- Applying L+ to multiple dets, now at No.", 1I5)') i
            call Lplus_single(inbasis(i,1), inbasis(i, 2), inepos(i, 1, :), &
                &inepos(i, 2, :), incoefs(i), iniszeros(i), &
                & basislist, coeflist, eposlist, iszerolist, num)
          enddo
      end subroutine Lplus_multiple

      subroutine Sminus_single(detup, detdn, eposup, eposdn, coef, iszero, &
              & basislist, coeflist, eposlist, iszerolist, num)
          integer(i16b), intent(in) :: detup, detdn
          integer, intent(in) :: eposup(:), eposdn(:)
          real(rk), intent(in) :: coef
          logical, intent(in) :: iszero
          logical, allocatable :: iszerolist(:)
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)!#, up/down, electron position
          integer :: num !number of basis determinants
          integer(i16b) :: tpdet, outdetup, outdetdn
          integer :: i
          integer, allocatable :: outeposup(:), outeposdn(:)
          logical :: outis0
          real(rk) :: outcoef
          ! write(*, '("----------------------- Apply S- on detup, detdn: ", 2B16)') detup, detdn
          ! write(*, '(15I3)') eposup(1:15)
          ! write(*, '(15I3)') eposdn(1:15)
          ! write(*, *) coef, iszero, num

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          allocate(outeposup(DET_MAX_LENGTH))
          allocate(outeposdn(DET_MAX_LENGTH))
          !apply s- on detup
          tpdet = detup
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdetup = detup
              outdetdn = detdn
              outeposup(:) = eposup(:)
              outeposdn(:) = eposdn(:)
              outcoef = coef
              outis0 = iszero.or.equals0(outcoef)
              call sms(i, outdetup, outdetdn, outeposup, outeposdn, outcoef, outis0)
              !outdet outepos outcoef outis0 are all updated here
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("s- on e at up pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/outdetup, outdetdn/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = outeposup(:)
                  eposlist(num, 2, :) = outeposdn(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
          deallocate(outeposup)
          deallocate(outeposdn)
      end subroutine Sminus_single

      subroutine Splus_single(detup, detdn, eposup, eposdn, coef, iszero, &
              & basislist, coeflist, eposlist, iszerolist, num)
          integer(i16b), intent(in) :: detup, detdn
          integer, intent(in) :: eposup(:), eposdn(:)
          real(rk), intent(in) :: coef
          logical, intent(in) :: iszero
          logical, allocatable :: iszerolist(:)
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)!#, up/down, electron position
          integer :: num !number of basis determinants
          integer(i16b) :: tpdet, outdetup, outdetdn
          integer :: i
          integer, allocatable :: outeposup(:), outeposdn(:)
          logical :: outis0
          real(rk) :: outcoef
          ! write(*, *) ' '
          ! write(*, '("--------------------- Apply S+ on detup, detdn: ", 2B16)') detup, detdn
          ! write(*, '(15I3)') eposup(1:15)
          ! write(*, '(15I3)') eposdn(1:15)
          ! write(*, *) coef, iszero, num

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          allocate(outeposup(DET_MAX_LENGTH))
          allocate(outeposdn(DET_MAX_LENGTH))
          !apply s+ on detdn
          tpdet = detdn
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdetup = detup
              outdetdn = detdn
              outeposup(:) = eposup(:)
              outeposdn(:) = eposdn(:)
              outcoef = coef
              outis0 = iszero.or.equals0(outcoef)
              call spls(i, outdetup, outdetdn, outeposup, outeposdn, outcoef, outis0)
              !outdet outepos outcoef outis0 are all updated here
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("s+ on e at up pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/outdetup, outdetdn/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = outeposup(:)
                  eposlist(num, 2, :) = outeposdn(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
          deallocate(outeposup)
          deallocate(outeposdn)
      end subroutine Splus_single

          !TODO add original det in L+single, L-single? may not need to
      subroutine Lminus_single(detup, detdn, eposup, eposdn, coef, iszero, &
              & basislist, coeflist, eposlist, iszerolist, num )
          integer(i16b), intent(in) :: detup, detdn
          integer, intent(in) :: eposup(:), eposdn(:)
          real(rk), intent(in) :: coef
          logical, intent(in) :: iszero
          logical, allocatable :: iszerolist(:)
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)!#, up/down, electron position
          integer :: num !number of basis determinants
          integer(i16b) :: tpdet, outdet
          integer :: i
          integer, allocatable :: outepos(:)
          logical :: outis0
          real(rk) :: outcoef
          ! write(*, '("---------------------- Apply L- on detup, detdn: ", 2B16)') detup, detdn
          ! write(*, '(15I3)') eposup(1:15)
          ! write(*, '(15I3)') eposdn(1:15)
          ! write(*, *) coef, iszero, num

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          allocate(outepos(DET_MAX_LENGTH))
          !apply l- on detup
          tpdet = detup
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdet = detup
              outepos(:) = eposup(:)
              outcoef = coef
              outis0 = iszero.or.equals0(outcoef)
              call lms(i, outdet, outepos, outcoef, outis0)
              !outdet outepos outcoef outis0 are all updated here
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("l- on e at up pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/outdet, detdn/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = outepos(:)
                  eposlist(num, 2, :) = eposdn(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
           !apply l+ on detdn
          tpdet = detdn
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdet = detdn
              outepos(:) = eposdn(:)
              outcoef = coef
              outis0 = iszero.or.equals0(coef)
              call lms(i, outdet, outepos, outcoef, outis0)
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("l- on e at dn pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/detup, outdet/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = eposup(:)
                  eposlist(num, 2, :) = outepos(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
         deallocate(outepos)

     end subroutine Lminus_single

      subroutine Lplus_single(detup, detdn, eposup, eposdn, coef, iszero, &
              & basislist, coeflist, eposlist, iszerolist, num )
          integer(i16b), intent(in) :: detup, detdn
          integer, intent(in) :: eposup(:), eposdn(:)
          real(rk), intent(in) :: coef
          logical, intent(in) :: iszero
          logical, allocatable :: iszerolist(:)
          integer(i16b), allocatable :: basislist(:, :)
          real(rk), allocatable :: coeflist(:)
          integer, allocatable :: eposlist(:, :, :)!#, up/down, electron position
          integer :: num !number of basis determinants
          integer(i16b) :: tpdet, outdet
          integer :: i
          integer, allocatable :: outepos(:)
          logical :: outis0
          real(rk) :: outcoef
          ! write(*, '("---------------------- Apply L+ on detup, detdn: ", 2B16)') detup, detdn
          ! write(*, '(15I3)') eposup(1:15)
          ! write(*, '(15I3)') eposdn(1:15)
          ! write(*, *) coef, iszero, num

          if(num.eq.0) then
              if(.not.allocated(basislist)) then
                  allocate(basislist(ARRAY_START_LENGTH, 2))
              endif
              if(.not.allocated(coeflist)) then
                  allocate(coeflist(ARRAY_START_LENGTH))
              endif
              if(.not.allocated(eposlist)) then
                  allocate(eposlist(ARRAY_START_LENGTH, 2, DET_MAX_LENGTH))
              endif
              if(.not.allocated(iszerolist)) then
                  allocate(iszerolist(ARRAY_START_LENGTH))
              endif
          endif
          allocate(outepos(DET_MAX_LENGTH))
          !apply l+ on detup
          tpdet = detup
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdet = detup
              outepos(:) = eposup(:)
              outcoef = coef
              outis0 = iszero.or.equals0(coef)
              call lpls(i, outdet, outepos, outcoef, outis0)
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("l+ on e at up pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/outdet, detdn/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = outepos(:)
                  eposlist(num, 2, :) = eposdn(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
           !apply l+ on detdn
          tpdet = detdn
          do while(tpdet.ne.0)
              i = trailz(tpdet) ! starting from 0
              outdet = detdn
              outepos(:) = eposdn(:)
              outcoef = coef
              outis0 = iszero.or.equals0(outcoef)
              call lpls(i, outdet, outepos, outcoef, outis0)
              !Original det is not added to the list here when outis0 = true
              ! write(*, '("l+ on e at dn pos:", 1I4)') i
              ! write(*, *) outis0
              if(.not.outis0) then
                  num = num + 1
                  basislist(num, 1:2) = (/detup, outdet/)
                  coeflist(num) = outcoef
                  eposlist(num, 1, :) = eposup(:)
                  eposlist(num, 2, :) = outepos(:)
                  iszerolist(num) = outis0
                  ! write(*, '("Output detup, detdn: ", 2B16)') basislist(num, 1:2)
                  ! write(*, '(15I3)') eposlist(num, 1, 1:15)
                  ! write(*, '(15I3)') eposlist(num, 2, 1:15)
                  ! write(*, *) coeflist(num), iszerolist(num), num
              endif

              tpdet = ibclr(tpdet, i)
          enddo
         deallocate(outepos)

     end subroutine Lplus_single

end module projection

max_stars_repo_name: MingtongH/csfcal

max_issues_repo_name: MingtongH/csfcal

max_forks_repo_name: MingtongH/csfcal


text:
!>
!! @par Copyright
!!
!! @brief Program to drive the MACv2-SP module "mo_simple_plumes"
!!
!! @author Karsten Peters, Stephanie Fiedler and Bjorn Stevens MPI-M, Hamburg (2015-12-05)
!!
!! @change-log:
!!          - 2016-12-05: beta release (KP, SF and BS, MPI-Met)
!!          - 2016-10-15: modified treatment of orography (BS)
!!          - 2016-11-09: changed lambda from  517.46442 to 550 nm 517.46442 nm used to make plots in paper (BS)
!!          - 2016-11-09: changed nlev from 48 to 80 and stretched orography, 48 used to make plots in paper (BS)
!!          - 2016-11-09: capability to output 2D fields (SF)
!!          - 2016-11-09: increasing level thickness with height, constant dz=500m used in plots for paper (KP)
!!
!! $ID: n/a$
!!
!!
!! @par Copyright
!!
!
PROGRAM MACv2SP

  USE mo_simple_plumes, ONLY: sp_aop_profile
  USE netcdf

  IMPLICIT NONE

  INTEGER, PARAMETER :: nlev = 80, nmon = 12

  REAL, PARAMETER    :: mons(nmon) = (/1., 2., 3., 4., 5., 6., 7., 8., 9.,10.,11.,12./)
  REAL, PARAMETER    :: lambda     = 550.0

  INTEGER ::           &
       nlat          , & !< number of latitudes, defined by orography file
       nlon          , & !< number of longitudes, defined by orography file
       iret          , & !< netCDF reading return variable
       ncid          , & !< netCDF file ID
       latID         , & !< pointer to latitude dimension in netCDF file
       lonID         , & !< pointer to longitude dimension in netCDF file
       VarID         , & !< pointer to generic dimension in netCDF file
       levID         , & !< pointer to level dimension in netCDF file
       monID         , & !< pointer to month dimension in netCDF file
       var_t_ID      , & !< pointer to time variable in netCDF file
       var_lat_ID    , & !< pointer to latitude variable in netCDF file
       var_lon_ID    , & !< pointer to longitude variable in netCDF file
       var_z_ID      , & !< pointer to height variable in netCDF file
       var_aod_ID    , & !< pointer to AOD variable in netCDF file
       var_asy_ID    , & !< pointer to asymmetry parameter variable in netCDF file
       var_ssa_ID    , & !< pointer to single-scattering albedo variable in netCDF file
       var_dNovrN_ID , & !< pointer to change in droplet number variable in netCDF file
       var_aod_ID_2d , & !< pointer to 2D AOD variable in netCDF file
       var_asy_ID_2d , & !< pointer to 2D asymmetry parameter variable in netCDF file
       var_ssa_ID_2d , & !< pointer to 2D single-scattering albedo variable in netCDF file
       iyear         , & !< index for year loop
       k             , & !< index for height loop
       i             , & !< index for latitude loop
       j             , & !< index for longitude loop
       m             , & !< index for month loop
       ifile             !< output file index

  REAL ::                    &
       lat_wght            , & !< latitudinal weights for global average
       global_aod          , & !< globally averaged aerosol optical depth (diagnosis)
       global_dNovrN       , & !< global averaged change in droplet number (diagnosis)
       year_fr

  REAL, ALLOCATABLE ::       &
       lat(:)              , & !< array of latitudes allocated based on nlat
       lon(:)              , & !< array of longitudes allocated based on nlon
       oro(:,:)            , & !< array of orographic heights (lon,lat)
       aod(:,:,:,:)        , & !< array for output aerosol optical depth (lon, lat, level, time)
       ssa(:,:,:,:)        , & !< array for output aerosol single scattering albedo (lon, lat, level, time)
       asy(:,:,:,:)        , & !< array for output aerosol asymmetry parameter (lon, lat, level, time)
       aod_int(:,:,:)      , & !< array for output 2D aerosol optical depth (lon, lat, time)
       ssa_int(:,:,:)      , & !< array for output 2D aerosol single scattering albedo (lon, lat, time)
       asy_int(:,:,:)      , & !< array for output 2D aerosol asymmetry parameter (lon, lat, time)
       dNovrN(:,:,:)       , & !< array for output of change in droplet number (lon, lat, level, time)
       aod_prof(:,:)       , & !< aerosol optical depth profiles by colums (ncolumns,level)
       ssa_prof(:,:)       , & !< aerosol single scattering albedo profiles by colums (ncolumns,level)
       asy_prof(:,:)       , & !< aerosol asymmetry parameter by colums (ncolumns,level)
       z(:,:)              , & !< heights by colums (ncolumns,level)
       dz(:,:)             , & !< layer thicknesses (ncolumns,level)
       col_lat(:)          , & !< latitudes of columns (ncolumns)
       col_lon(:)              !< longitude of columns (ncolumns)
  !
  ! read in the orography data and use this to define the input (lat/lon) grid
  !
  iret = nf90_open("./orography_T63.nc", NF90_NOWRITE, ncid)
  IF (iret /= NF90_NOERR) STOP 'Error in opening orography file'
  iret = nf90_inq_dimid(ncid, "lat", VarID)
  iret = nf90_inquire_dimension(ncid, VarID, len = nlat)
  iret = nf90_inq_dimid(ncid, "lon", VarID)
  iret = nf90_inquire_dimension(ncid, VarID, len = nlon)

  ALLOCATE (lat(nlat))
  ALLOCATE (lon(nlon))
  ALLOCATE (oro(nlon,nlat))
  ALLOCATE (aod(nlon,nlat,nlev,nmon))
  ALLOCATE (ssa(nlon,nlat,nlev,nmon))
  ALLOCATE (asy(nlon,nlat,nlev,nmon))
  ALLOCATE (dNovrN(nlon,nlat,nmon))
  ALLOCATE (aod_int(nlon,nlat,nmon))
  ALLOCATE (ssa_int(nlon,nlat,nmon))
  ALLOCATE (asy_int(nlon,nlat,nmon))
  ALLOCATE (z(nlon,nlev))
  ALLOCATE (dz(nlon,nlev))
  ALLOCATE (col_lat(nlon))
  ALLOCATE (col_lon(nlon))
  ALLOCATE (aod_prof(nlon,nlev))
  ALLOCATE (ssa_prof(nlon,nlev))
  ALLOCATE (asy_prof(nlon,nlev))

  iret = nf90_inq_varid(ncid, "lat", VarID)
  iret = nf90_get_var(ncid, VarID, lat(:)  , start=(/1/)  ,count=(/nlat/))
  IF (iret /= NF90_NOERR) STOP 'Error in reading latitudes'
  iret = nf90_inq_varid(ncid, "lon", VarID)
  iret = nf90_get_var(ncid, VarID, lon(:)  , start=(/1/)  ,count=(/nlon/))
  IF (iret /= NF90_NOERR) STOP 'Error in reading longitudes'
  iret = nf90_inq_varid(ncid, "asl", VarID)
  iret = nf90_get_var(ncid, VarID, oro(:,:), start=(/1,1/),count=(/nlon,nlat/))
  IF (iret /= NF90_NOERR) STOP 'Error in reading orographic height'
  iret = nf90_close  (ncid)
  !
  ! define the height array
  !
  z(:,1)  = 0.
  DO k = 2,nlev
      dz(:,k) = 50. * exp(0.03*(k-1))
      z(:,k) = z(:,k-1) + dz(:,k)
  END DO
  !
  ! loop over annual cycle for given year, or year rangeand calculate aerosol optical properties at each lat and
  ! lon point
  !
  DO iyear = 2005,2005
    DO m = 1,nmon
      year_fr = iyear + (m-0.5) / nmon

      DO i = 1,nlat
        !
        ! here we process ncolumns = nlon in one call to sp_aop_profile, output is written to 2D arrays and then
        ! transfered to 4D output arrays to avoid striding through memory in call to sp_aop_profile
        !
        col_lat(:) = lat(i)
        col_lon(:) = lon(:)
        CALL sp_aop_profile                                                               ( &
             nlev         ,nlon       ,lambda     ,oro(:,i)     ,col_lon      ,col_lat    , &
             year_fr      ,z(:,:)     ,dz(:,:)    ,dNovrN(:,i,m),aod_prof     ,ssa_prof   , &
             asy_prof     )
        DO j = 1,nlon
          DO k=1,nlev
            aod(j,i,k,m) = aod_prof(j,k)
            ssa(j,i,k,m) = ssa_prof(j,k)
            asy(j,i,k,m) = asy_prof(j,k)
          END DO
        END DO
        !
        ! 2D vertically integrated values
        !
        DO j = 1,nlon
          aod_int(j,i,m) = SUM(aod(j,i,:,m))
          ssa_int(j,i,m) = SUM(aod(j,i,:,m)*ssa(j,i,:,m))/aod_int(j,i,m)
          asy_int(j,i,m) = SUM(aod(j,i,:,m)*ssa(j,i,:,m)*asy(j,i,:,m))/(ssa_int(j,i,m)*aod_int(j,i,m))
        END DO
      END DO
    END DO
    !
    ! compute global and monthly average of selected quantities for diagnostic output
    !
    global_aod    = 0.
    global_dNovrN = 0.
    lat_wght      = 0.
    DO i=1,nlat
      lat_wght      = lat_wght      + COS(lat(i)*ASIN(1.)/90.)
      global_aod    = global_aod    + SUM(aod(:,i,:,:))*COS(lat(i)*ASIN(1.)/90.)
      global_dNovrN = global_dNovrN + SUM(dNovrN(:,i,:))*COS(lat(i)*ASIN(1.)/90.)
    END DO

    global_aod    = 0.
    global_dNovrN = 0.
    lat_wght      = 0.
    DO i=1,nlat/2
      lat_wght      = lat_wght      + COS(lat(i)*ASIN(1.)/90.)
      global_aod    = global_aod    + SUM(aod(:,i,:,:))*COS(lat(i)*ASIN(1.)/90.)
      global_dNovrN = global_dNovrN + SUM(dNovrN(:,i,:))*COS(lat(i)*ASIN(1.)/90.)
    END DO
    !
    ! write output as a netcdf file
    !
    IF (iyear == 2005) THEN
      iret = NF90_NOERR
      iret = iret + nf90_create("./MACv2-SP_3Dfields_T63_1km_2005.nc", NF90_CLOBBER, ncid)
      iret = iret + nf90_def_dim(ncid, 'lat'   ,nlat , latID)
      iret = iret + nf90_def_dim(ncid, 'lon'   ,nlon , lonID)
      iret = iret + nf90_def_dim(ncid, 'time'  ,nmon , monID)
      iret = iret + nf90_def_dim(ncid, 'z'     ,nlev , levID)
      IF (iret /= 6*NF90_NOERR) STOP 'Error in Creating File Dimensions'
      !
      iret = NF90_NOERR
      iret = iret + nf90_def_var(ncid, 'time'         , NF90_FLOAT, monID, var_t_ID)
      iret = iret + nf90_def_var(ncid, 'lat'          , NF90_FLOAT, latID, var_lat_ID)
      iret = iret + nf90_def_var(ncid, 'lon'          , NF90_FLOAT, lonID, var_lon_ID)
      iret = iret + nf90_def_var(ncid, 'z'            , NF90_FLOAT, levID, var_z_ID)
      iret = iret + nf90_def_var(ncid, 'aod'          , NF90_FLOAT, (/lonID,latID,levID,monID/), var_aod_ID)
      iret = iret + nf90_def_var(ncid, 'ssa'          , NF90_FLOAT, (/lonID,latID,levID,monID/), var_ssa_ID)
      iret = iret + nf90_def_var(ncid, 'asy'          , NF90_FLOAT, (/lonID,latID,levID,monID/), var_asy_ID)
      iret = iret + nf90_def_var(ncid, 'dNovrN'       , NF90_FLOAT, (/lonID,latID,monID/), var_dNovrN_ID)
      iret = iret + nf90_def_var(ncid, 'aod_2D'       , NF90_FLOAT, (/lonID,latID,monID/), var_aod_ID_2d)
      iret = iret + nf90_def_var(ncid, 'ssa_2D'       , NF90_FLOAT, (/lonID,latID,monID/), var_ssa_ID_2d)
      iret = iret + nf90_def_var(ncid, 'asy_2D'       , NF90_FLOAT, (/lonID,latID,monID/), var_asy_ID_2d)

      iret = iret + nf90_put_att(ncid, var_t_ID       , "long_name", "month of year")
      iret = iret + nf90_put_att(ncid, var_t_ID       , "units"    , "month")
      iret = iret + nf90_put_att(ncid, var_lat_ID     , "long_name", "latitude")
      iret = iret + nf90_put_att(ncid, var_lat_ID     , "units"    , "degrees_north")
      iret = iret + nf90_put_att(ncid, var_lon_ID     , "long_name", "longitude")
      iret = iret + nf90_put_att(ncid, var_lon_ID     , "units"    , "degrees_east")
      iret = iret + nf90_put_att(ncid, var_z_ID       , "long_name", "height above sea-level")
      iret = iret + nf90_put_att(ncid, var_z_ID       , "units"    , "m")
      iret = iret + nf90_put_att(ncid, var_aod_ID     , "long_name", "aerosol optical depth")
      iret = iret + nf90_put_att(ncid, var_aod_ID     , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_ssa_ID     , "long_name", "aerosol single scattering albedo")
      iret = iret + nf90_put_att(ncid, var_ssa_ID     , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_asy_ID     , "long_name", "aerosol asymmetry parameter")
      iret = iret + nf90_put_att(ncid, var_asy_ID     , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_dNovrN_ID  , "long_name", "normalized change in drop number")
      iret = iret + nf90_put_att(ncid, var_dNovrN_ID  , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_aod_ID_2d  , "long_name", "2D aerosol optical depth")
      iret = iret + nf90_put_att(ncid, var_aod_ID_2d  , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_ssa_ID_2d  , "long_name", "2D aerosol single scattering albedo")
      iret = iret + nf90_put_att(ncid, var_ssa_ID_2d  , "units"    ," ")
      iret = iret + nf90_put_att(ncid, var_asy_ID_2d  , "long_name", "2D aerosol asymmetry parameter")
      iret = iret + nf90_put_att(ncid, var_asy_ID_2d  , "units"    ," ")
      iret = iret + nf90_enddef(ncid)
      IF (iret /= 26*NF90_NOERR) STOP 'Error in creating file variables'
      !
      iret = NF90_NOERR
      iret = iret + nf90_put_var(ncid, var_t_ID       , values=mons)
      iret = iret + nf90_put_var(ncid, var_lat_ID     , values=lat)
      iret = iret + nf90_put_var(ncid, var_lon_ID     , values=lon)
      iret = iret + nf90_put_var(ncid, var_z_ID       , values=z(1,:))
      iret = iret + nf90_put_var(ncid, var_aod_ID     , values=aod)
      iret = iret + nf90_put_var(ncid, var_ssa_ID     , values=ssa)
      iret = iret + nf90_put_var(ncid, var_asy_ID     , values=asy)
      iret = iret + nf90_put_var(ncid, var_dNovrN_ID  , values=dNovrN)
      iret = iret + nf90_put_var(ncid, var_aod_ID_2d  , values=aod_int)
      iret = iret + nf90_put_var(ncid, var_ssa_ID_2d  , values=ssa_int)
      iret = iret + nf90_put_var(ncid, var_asy_ID_2d  , values=asy_int)
      iret = iret + nf90_close(ncid)
      IF (iret /= 10*NF90_NOERR) STOP 'error writing data or in closing file'
    END IF

  END DO

END PROGRAM MACv2SP

max_stars_repo_name: fdisante/RegCM

max_issues_repo_name: fdisante/RegCM

max_forks_repo_name: fdisante/RegCM