text:
      subroutine geneqmdnode
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C**********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 To generate equations fo multiply defined nodes
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 04-04-96     G. Zyvoloski   00022   Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/geneqmdnode.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:06   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:06:04   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:09:22   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:24:06   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:02:28   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:41:40 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.3   Thu Jun 06 15:10:48 1996   robinson
CD2 Fixed implementation for dpdp
CD2
CD2    Rev 1.2   Wed May 22 08:55:58 1996   gaz
CD2 corrections and additions
CD2
C**********************************************************************
CD3
CD3 REQUIREMENTS TRACEABILITY
CD3
CD3 2.3.2 Heat- and mass-transfer equations
CD3 2.3.3 Noncondensible gas flow equations
CD3
C**********************************************************************
CD4
CD4 SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C**********************************************************************

      use comflow
      use davidi
      use comji
      use comfi
      use comgi
      use comei
      use comdi
      use comci
      use combi
      use comdti
      use comai
      implicit none

      integer i,i1,i2,icd,ii,jj,ii1,ii2,idg,iq,jmi,jml,jmia,jm,neqp1
      integer ij,ij1,ij2,iz,kb,neighc,iau,ial,kz,nmatavw,iz4m1
      real*8 sx1d,pvii,phii,dlpi,dlpkb,dvpi,dvpkb,dlei,dvei,dlekb
      real*8 dvekb,dlapi,dvapi,dlapkb,dvapkb,pvikb,phikb,axyd,axy,axyf
      real*8 pvxy,pxy,pxyh,pxyi,vxyd,vxy,vxyf,sx4d,sx4h,dpvti
      real*8 dlaei,dlaekb,dvaei,dvaekb,sx_md,fid,fid1,dili,dilpi,dilei
      real*8 dilkb,dilpkb,dilekb,divi,divpi,divei,divkb,divpkb,divekb
      real*8 heatc,ti,cnvi,dcvi,dcvei,dcvci,cnvkb,dcvkb,dcvekb,dcvckb
      real*8 dlci,dlckb,dlaci,dlackb,dvci,dvckb,dvaci,dvackb
      logical bit

      neqp1=neq+1

c     loop on nodes
c
      do i=1,n
c
c     form constants for i>neq
c
         if(i.gt.neq.and.idualp.eq.0) then
            icd=neq
         else
            icd=0
         endif
c
c     loop on mdnodes
c
         do jj=1,abs(mdnodes(i-icd))
            ii=mdnode(i-icd,jj)
            if(ii.gt.i-icd) then
c
c     storage locations for fluxes
c
               if(i.gt.neq) then
                  nmatavw=ldna
               else
                  nmatavw=0
               endif
c
c     node i value of parameters
c
               sx1d=sx1(i)
               if(i.gt.neq.and.idpdp.ne.0) then
c     means double permeability
                  if(ico2.le.0) then
                     nmat(1)=ldna*(11-1)
                     nmat(2)=ldna*(12-1)
                     nmat(3)=ldna*(15-1)
                     nmat(4)=ldna*(16-1)
                     nrhs(1)=(3-1)*neq
                     nrhs(2)=(4-1)*neq
                  else if(ico2.gt.0) then
                     nmat(1)=ldna*(22-1)
                     nmat(2)=ldna*(23-1)
                     nmat(3)=ldna*(24-1)
                     nmat(4)=ldna*(28-1)
                     nmat(5)=ldna*(29-1)
                     nmat(6)=ldna*(30-1)
                     nmat(7)=ldna*(34-1)
                     nmat(8)=ldna*(35-1)
                     nmat(9)=ldna*(36-1)
                     nrhs(1)=(4-1)*neq
                     nrhs(2)=(5-1)*neq
                     nrhs(3)=(6-1)*neq
                  endif
               else if(idpdp.ne.0) then
                  if(ico2.le.0) then
                     nmat(1)=0
                     nmat(2)=ldna
                     nmat(3)=ldna*(5-1)
                     nmat(4)=ldna*(6-1)
                     nrhs(1)=0
                     nrhs(2)=neq
                  else if(ico2.gt.0) then
                     nmat(1)=0
                     nmat(2)=ldna
                     nmat(3)=ldna*2
                     nmat(4)=ldna*(7-1)
                     nmat(5)=ldna*(8-1)
                     nmat(6)=ldna*(9-1)
                     nmat(7)=ldna*(13-1)
                     nmat(8)=ldna*(14-1)
                     nmat(9)=ldna*(15-1)
                     nrhs(1)=0
                     nrhs(2)=neq
                     nrhs(3)=neq+neq
                  endif
               endif
               iz=i-icd
c
c     calculate flow coefficient
c
               sx_md=0.0
               i1=nelm(iz)+1
               i2=nelm(iz+1)
               do jm=i1,i2
                  iw=istrw(jm-neqp1)
                  sx_md=sx_md+(sx(iw,1)+sx(iw,2)+sx(iw,3))
               enddo
               i1=nelm(ii)+1
               i2=nelm(ii+1)
               do jm=i1,i2
                  iw=istrw(jm-neqp1)
                  sx_md=sx_md+(sx(iw,1)+sx(iw,2)+sx(iw,3))
               enddo
               sx_md=sx_md*sx_mult
               iz4m1 = 4*(iz-1)+1
               ii1=nelm(iz)+1
               ii2=nelm(iz+1)
               idg=nelmdg(iz)-ii1+1
               neqp1=neq+1
               jmi=nelmdg(iz)
               jmia=jmi-neqp1

c     Take care of source/sink term by putting it into empty slot (from
c     same node to same node).
c     If this is an isothermal air-water simulation
               a_axy(jmia+nmatavw)=sk(i)
c     Take care of source/sink term
c     If this is an isothermal air-water simulation
               a_vxy(jmia+nmatavw)=qh(i)

               iq=0
               do  jm=jmi+1,ii2
                  if(nelm(jm).eq.ii) then
                     iq=iq+1
                     it8(iq)=nelm(jm)+icd
                     it9(iq)=jm-ii1+1
                     it10(iq)=istrw(jm-neqp1)
                     it11(iq)=jm-neqp1
                     ij1=nelm(nelm(jm))+1
                     ij2=nelmdg(nelm(jm))-1
                     do  ij=ij1,ij2
                        if(nelm(ij).eq.iz) then
                           it12(iq)=ij-neqp1
                        endif
                     enddo
                  endif
               enddo
c
c     if block on degree of freedom
c
               if(idof.gt.1) then
c     non heat conduction only
c
c     ifblock on equation type
c
                  pvii=phi(i)
                  phii=pvii-pcp(i)
                  dpvti=dpcef(i)
                  if(ico2.lt.0) then

c     3-d and 2-d geometry
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iw=it10(jm)
                     pvikb=phi(kb)
                     phikb=pvikb-pcp(kb)
                     pxy=sx_md*
     &               max(pnx(i),pny(i),pnz(i),
     &               pnx(kb),pny(kb),pnz(kb))
                     pvxy=pxy
                     pxyi=pxy*(phikb-phii)
                     pxyh=pvxy*(pvikb-pvii)
                     t1(neighc)=pxyi
                     t2(neighc)=pxyh
                     t3(neighc)=pxy
                     t4(neighc)=pvxy
                     t5(neighc)=sx_md*max(thx(i),thy(i),thz(i),
     &               thx(kb),thy(kb),thz(kb))
                     t6(neighc)=-grav*t3(neighc)
                     t7(neighc)=-grav*t4(neighc)
c
c     liquid phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyi=t1(neighc)
                     sx4d=t6(neighc)
                     axyd=pxyi+0.5*sx4d*(rolf(i)+rolf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=axyd
c
c     determine upwind nodes and if liquid phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        axyd=t8(neighc)
                        if(axyd.lt.0.0) fid=dnwgt
                        if(axyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_l(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_l(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     dili=s(i)
                     dilpi=0.0
                     dilei=1.0
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyi=t1(neighc)
                     pxy=t3(neighc)
                     sx4d=t6(neighc)
                     axyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     dilkb=s(kb)
                     dilpkb=0.0
                     dilekb=1.0
                     dlpi=-pxy+0.5*sx4d*dglp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlpkb=pxy+0.5*sx4d*dglp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlei=pxy*dpvti+0.5*sx4d*dgle(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlekb=-pxy*dpcef(kb)+0.5*sx4d*dgle(kb)
     *                    *(cord(kz,igrav)-cord(iz,igrav))

                     axyf=(fid*dilkb+fid1*dili)
                     axy=axyd*axyf
                     dlapi=dlpi*axyf+axyd*fid1*dilpi
                     dlapkb=dlpkb*axyf+axyd*fid*dilpkb
                     dlaei=dlei*axyf+axyd*fid1*dilei
                     dlaekb=dlekb*axyf+axyd*fid*dilekb

                     a_axy(iau+nmatavw)=axy
                     a_axy(ial+nmatavw)=-axy

                     bp(iz+nrhs(1))=bp(iz+nrhs(1))+axy
                     bp(kz+nrhs(1))=bp(kz+nrhs(1))-axy
                     a(jmia+nmat(1))=a(jmia+nmat(1))+dlapi
                     a(jmia+nmat(2))=a(jmia+nmat(2))+dlaei
                     a(ial+nmat(1))=a(ial+nmat(1))-dlapi
                     a(ial+nmat(2))=a(ial+nmat(2))-dlaei
                     a(iau+nmat(1))=a(iau+nmat(1))+dlapkb
                     a(iau+nmat(2))=a(iau+nmat(2))+dlaekb
                     a(jml+nmat(1))=a(jml+nmat(1))-dlapkb
                     a(jml+nmat(2))=a(jml+nmat(2))-dlaekb
c
c     vapour phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyh=t2(neighc)
                     sx4h=t7(neighc)
                     vxyd=pxyh+0.5*sx4h*(rovf(i)+rovf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=vxyd
c
c     determine upwind nodes and if vapour phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        vxyd=t8(neighc)
                        if(vxyd.lt.0.0) fid=dnwgt
                        if(vxyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_v(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_v(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     divi=1.0-s(i)
                     divpi=0.0
                     divei=-1.0
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyh=t2(neighc)
                     pvxy=t4(neighc)
                     sx4h=t7(neighc)
                     vxyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     divkb=1.0-s(kb)
                     divpkb=0.0
                     divekb=-1.0
                     dvpi=-pvxy+0.5*sx4h*dgvp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvpkb=pvxy+0.5*sx4h*dgvp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvei=0.5*sx4h*dgve(i)*(cord(kz,igrav)-
     &                    cord(iz,igrav))
                     dvekb=0.5*sx4h*dgve(kb)
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     vxyf=(fid*divkb+fid1*divi)
                     vxy=vxyd*vxyf
                     dvapi=dvpi*vxyf+vxyd*fid1*divpi
                     dvapkb=dvpkb*vxyf+vxyd*fid*divpkb
                     dvaei=dvei*vxyf+vxyd*fid1*divei
                     dvaekb=dvekb*vxyf+vxyd*fid*divekb

                     a_vxy(iau+nmatavw)=vxy
                     a_vxy(ial+nmatavw)=-vxy

                     bp(iz+nrhs(2))=bp(iz+nrhs(2))+vxy
                     bp(kz+nrhs(2))=bp(kz+nrhs(2))-vxy
                     a(jmia+nmat(3))=a(jmia+nmat(3))+dvapi
                     a(jmia+nmat(4))=a(jmia+nmat(4))+dvaei
                     a(ial+nmat(3))=a(ial+nmat(3))-dvapi
                     a(ial+nmat(4))=a(ial+nmat(4))-dvaei
                     a(iau+nmat(3))=a(iau+nmat(3))+dvapkb
                     a(iau+nmat(4))=a(iau+nmat(4))+dvaekb
                     a(jml+nmat(3))=a(jml+nmat(3))-dvapkb
                     a(jml+nmat(4))=a(jml+nmat(4))-dvaekb
                  else if(ico2.eq.0.and.idof.ne.1) then
c
c     ifblock for heat/water (no air) equation type
c
c     3-d and 2-d geometry
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iw=it10(jm)
                     pvikb=phi(kb)
                     phikb=pvikb-pcp(kb)
                     pxy=sx_md*
     &               max(pnx(i),pny(i),pnz(i),
     &               pnx(kb),pny(kb),pnz(kb))
                     pvxy=pxy
                     pxyi=pxy*(phikb-phii)
                     pxyh=pvxy*(pvikb-pvii)
                     t1(neighc)=pxyi
                     t2(neighc)=pxyh
                     t3(neighc)=pxy
                     t4(neighc)=pvxy
                     t5(neighc)=sx_md*max(thx(i),thy(i),thz(i),
     &               thx(kb),thy(kb),thz(kb))
                     t6(neighc)=-grav*t3(neighc)
                     t7(neighc)=-grav*t4(neighc)
c
c     liquid phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyi=t1(neighc)
                     sx4d=t6(neighc)
                     axyd=pxyi+0.5*sx4d*(rolf(i)+rolf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=axyd
c
c     determine upwind nodes and if liquid phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        axyd=t8(neighc)
                        if(axyd.lt.0.0) fid=dnwgt
                        if(axyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_l(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_l(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     dili=s(i)
                     dilpi=0.0
                     if(ieos(i).ne.2) then
                        dilei=0.0
                     else
                        dilei=1.0
                     endif
c
c     enli=0.01*t(i)
c     deli=0.01*dtpa(i)
c     delei=0.01*dtpae(i)
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyi=t1(neighc)
                     pxy=t3(neighc)
                     sx4d=t6(neighc)
                     axyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     dilkb=s(kb)
                     dilpkb=0.0
                     if(ieos(kb).ne.2) then
                        dilekb=0.0
                     else
                        dilekb=1.0
                     endif
                     dlpi=-pxy+0.5*sx4d*dglp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlpkb=pxy+0.5*sx4d*dglp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlei=pxy*dpvti+0.5*sx4d*dgle(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlekb=-pxy*dpcef(kb)+0.5*sx4d*dgle(kb)
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     axyf=(fid*dilkb+fid1*dili)
                     axy=axyd*axyf
                     dlapi=dlpi*axyf+axyd*fid1*dilpi
                     dlapkb=dlpkb*axyf+axyd*fid*dilpkb
                     dlaei=dlei*axyf+axyd*fid1*dilei
                     dlaekb=dlekb*axyf+axyd*fid*dilekb

                     a_axy(iau+nmatavw)=axy
                     a_axy(ial+nmatavw)=-axy

                     bp(iz+nrhs(1))=bp(iz+nrhs(1))+axy
                     bp(kz+nrhs(1))=bp(kz+nrhs(1))-axy
                     a(jmia+nmat(1))=a(jmia+nmat(1))+dlapi
                     a(jmia+nmat(2))=a(jmia+nmat(2))+dlaei
                     a(ial+nmat(1))=a(ial+nmat(1))-dlapi
                     a(ial+nmat(2))=a(ial+nmat(2))-dlaei
                     a(iau+nmat(1))=a(iau+nmat(1))+dlapkb
                     a(iau+nmat(2))=a(iau+nmat(2))+dlaekb
                     a(jml+nmat(1))=a(jml+nmat(1))-dlapkb
                     a(jml+nmat(2))=a(jml+nmat(2))-dlaekb
c
c
c     vapour phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyh=t2(neighc)
                     sx4h=t7(neighc)
                     vxyd=pxyh+0.5*sx4h*(rovf(i)+rovf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=vxyd
c
c     determine upwind nodes and if vapour phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        vxyd=t8(neighc)
                        if(vxyd.lt.0.0) fid=dnwgt
                        if(vxyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_v(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_v(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     divi=1.0-s(i)
                     divpi=0.0
                     if(ieos(i).ne.2) then
                        divei=0.0
                     else
                        divei=-1.0
                     endif
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyh=t2(neighc)
                     pvxy=t4(neighc)
                     sx4h=t7(neighc)
                     vxyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     divkb=1.0-s(kb)
                     divpkb=0.0
                     if(ieos(kb).ne.2) then
                        divekb=0.0
                     else
                        divekb=-1.0
                     endif
                     dvpi=-pvxy+0.5*sx4h*dgvp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvpkb=pvxy+0.5*sx4h*dgvp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvei=0.5*sx4h*dgve(i)*(cord(kz,igrav)-
     &                    cord(iz,igrav))
                     dvekb=0.5*sx4h*dgve(kb)
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     vxyf=(fid*divkb+fid1*divi)
                     vxy=vxyd*vxyf
                     dvapi=dvpi*vxyf+vxyd*fid1*divpi
                     dvapkb=dvpkb*vxyf+vxyd*fid*divpkb
                     dvaei=dvei*vxyf+vxyd*fid1*divei
                     dvaekb=dvekb*vxyf+vxyd*fid*divekb

                     a_vxy(iau+nmatavw)=vxy
                     a_vxy(ial+nmatavw)=-vxy

                     bp(iz+nrhs(1))=bp(iz+nrhs(1))+vxy
                     bp(kz+nrhs(1))=bp(kz+nrhs(1))-vxy
                     a(jmia+nmat(1))=a(jmia+nmat(1))+dvapi
                     a(jmia+nmat(2))=a(jmia+nmat(2))+dvaei
                     a(ial+nmat(1))=a(ial+nmat(1))-dvapi
                     a(ial+nmat(2))=a(ial+nmat(2))-dvaei
                     a(iau+nmat(1))=a(iau+nmat(1))+dvapkb
                     a(iau+nmat(2))=a(iau+nmat(2))+dvaekb
                     a(jml+nmat(1))=a(jml+nmat(1))-dvapkb
                     a(jml+nmat(2))=a(jml+nmat(2))-dvaekb
c
c     add heat conduction
c
                     ti=t(i)
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     heatc=t5(neighc)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     bp(iz+nrhs(2))=bp(iz+nrhs(2))+heatc*(t(kb)-ti)
                     bp(kz+nrhs(2))=bp(kz+nrhs(2))-heatc*(t(kb)-ti)
                     a(jmia+nmat(3))=a(jmia+nmat(3))-heatc*dtpa(i)
                     a(jmia+nmat(4))=a(jmia+nmat(4))-heatc*dtpae(i)
                     a(ial+nmat(3))=a(ial+nmat(3))+heatc*dtpa(i)
                     a(ial+nmat(4))=a(ial+nmat(4))+heatc*dtpae(i)
                     a(iau+nmat(3))=a(iau+nmat(3))+heatc*dtpa(kb)
                     a(iau+nmat(4))=a(iau+nmat(4))+heatc*dtpae(kb)
                     a(jml+nmat(3))=a(jml+nmat(3))-heatc*dtpa(kb)
                     a(jml+nmat(4))=a(jml+nmat(4))-heatc*dtpae(kb)
c     end ifblock for heat/water (no air) equation type
                  else if(ico2.gt.0) then
c
c     ifblock for heat/water/air equation type
c
c     3-d and 2-d geometry
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iw=it10(jm)
                     pvikb=phi(kb)
                     phikb=pvikb-pcp(kb)
                     pxy=sx_md*
     &               max(pnx(i),pny(i),pnz(i),
     &               pnx(kb),pny(kb),pnz(kb))
                     pvxy=pxy
                     pxyi=pxy*(phikb-phii)
                     pxyh=pvxy*(pvikb-pvii)
                     t1(neighc)=pxyi
                     t2(neighc)=pxyh
                     t3(neighc)=pxy
                     t4(neighc)=pvxy
                     t5(neighc)=sx_md*max(thx(i),thy(i),thz(i),
     &               thx(kb),thy(kb),thz(kb))
                     t6(neighc)=-grav*t3(neighc)
                     t7(neighc)=-grav*t4(neighc)
c
c     liquid phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyi=t1(neighc)
                     sx4d=t6(neighc)
                     axyd=pxyi+0.5*sx4d*(rolf(i)+rolf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=axyd
c
c     determine upwind nodes and if liquid phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        axyd=t8(neighc)
                        if(axyd.lt.0.0) fid=dnwgt
                        if(axyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_l(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_l(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     dili=s(i)
                     dilpi=0.0
                     if(ieos(i).ne.2) then
                        dilei=0.0
                     else
                        dilei=1.0
                     endif
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyi=t1(neighc)
                     pxy=t3(neighc)
                     sx4d=t6(neighc)
                     axyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     dilkb=s(kb)
                     dilpkb=0.0
                     if(ieos(kb).ne.2) then
                        dilekb=0.0
                     else
                        dilekb=1.0
                     endif
                     dlpi=-pxy+0.5*sx4d*dglp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlpkb=pxy+0.5*sx4d*dglp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlei=pxy*dpvti+0.5*sx4d*dgle(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlekb=-pxy*dpcef(kb)+0.5*sx4d*dgle(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlci=0.5*sx4d*dglc(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dlckb=0.5*sx4d*dglc(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     axyf=(fid*dilkb+fid1*dili)
                     axy=axyd*axyf
                     dlapi=dlpi*axyf+axyd*fid1*dilpi
                     dlapkb=dlpkb*axyf+axyd*fid*dilpkb
                     dlaei=dlei*axyf+axyd*fid1*dilei
                     dlaekb=dlekb*axyf+axyd*fid*dilekb
                     dlaci=dlci*axyf
                     dlackb=dlckb*axyf

                     a_axy(iau+nmatavw)=axy
                     a_axy(ial+nmatavw)=-axy

                     bp(iz+nrhs(1))=bp(iz+nrhs(1))+axy
                     bp(kz+nrhs(1))=bp(kz+nrhs(1))-axy
                     a(jmia+nmat(1))=a(jmia+nmat(1))+dlapi
                     a(jmia+nmat(2))=a(jmia+nmat(2))+dlaei
                     a(jmia+nmat(3))=a(jmia+nmat(3))+dlaci
                     a(ial+nmat(1))=a(ial+nmat(1))-dlapi
                     a(ial+nmat(2))=a(ial+nmat(2))-dlaei
                     a(ial+nmat(3))=a(ial+nmat(3))-dlaci
                     a(iau+nmat(1))=a(iau+nmat(1))+dlapkb
                     a(iau+nmat(2))=a(iau+nmat(2))+dlaekb
                     a(iau+nmat(3))=a(iau+nmat(3))+dlackb
                     a(jml+nmat(1))=a(jml+nmat(1))-dlapkb
                     a(jml+nmat(2))=a(jml+nmat(2))-dlaekb
                     a(jml+nmat(3))=a(jml+nmat(3))-dlackb
c
c
c     vapour phase calculations
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     pxyh=t2(neighc)
                     sx4h=t7(neighc)
                     vxyd=pxyh+0.5*sx4h*(rovf(i)+rovf(kb))
     *                    *(cord(kz,igrav)-cord(iz,igrav))
                     t8(neighc)=vxyd
c
c     determine upwind nodes and if vapour phase exists
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
c     add coding to save upwind position
                     if(iad.le.iad_up) then
                        fid=0.5
                        vxyd=t8(neighc)
                        if(vxyd.lt.0.0) fid=dnwgt
                        if(vxyd.gt.0.0) fid=upwgt
                        t9(neighc)=fid
                        call setbit(nbits,neighc,upwind_v(iz4m1),fid)
                     else
                        if(bit(nbits,neighc,upwind_v(iz4m1))) then
                           t9(neighc)=1.0
                        else
                           t9(neighc)=0.0
                        endif
                     endif
c
c     form equations
c
                     divi=1.0-s(i)
                     divpi=0.0
                     if(ieos(i).ne.2) then
                        divei=0.0
                     else
                        divei=-1.0
                     endif
c
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     pxyh=t2(neighc)
                     pvxy=t4(neighc)
                     sx4h=t7(neighc)
                     vxyd=t8(neighc)
                     fid=t9(neighc)
                     fid1=1.0-fid
                     divkb=1.0-s(kb)
                     divpkb=0.0
                     if(ieos(kb).ne.2) then
                        divekb=0.0
                     else
                        divekb=-1.0
                     endif
                     dvpi=-pvxy+0.5*sx4h*dgvp(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvpkb=pvxy+0.5*sx4h*dgvp(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvei=0.5*sx4h*dgve(i)*(cord(kz,igrav)-
     &                    cord(iz,igrav))
                     dvekb=0.5*sx4h*dgve(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvci=0.5*sx4h*dgvc(i)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     dvckb=0.5*sx4h*dgvc(kb)*
     *                    (cord(kz,igrav)-cord(iz,igrav))
                     vxyf=(fid*divkb+fid1*divi)
                     vxy=vxyd*vxyf
                     dvapi=dvpi*vxyf+vxyd*fid1*divpi
                     dvapkb=dvpkb*vxyf+vxyd*fid*divpkb
                     dvaei=dvei*vxyf+vxyd*fid1*divei
                     dvaekb=dvekb*vxyf+vxyd*fid*divekb
                     dvaci=dvci*vxyf
                     dvackb=dvckb*vxyf

                     a_vxy(iau+nmatavw)=vxy
                     a_vxy(ial+nmatavw)=-vxy

                     bp(iz+nrhs(1))=bp(iz+nrhs(1))+vxy
                     bp(kz+nrhs(1))=bp(kz+nrhs(1))-vxy
                     a(jmia+nmat(1))=a(jmia+nmat(1))+dvapi
                     a(jmia+nmat(2))=a(jmia+nmat(2))+dvaei
                     a(jmia+nmat(3))=a(jmia+nmat(3))+dvaci
                     a(ial+nmat(1))=a(ial+nmat(1))-dvapi
                     a(ial+nmat(2))=a(ial+nmat(2))-dvaei
                     a(ial+nmat(3))=a(ial+nmat(3))-dvaci
                     a(iau+nmat(1))=a(iau+nmat(1))+dvapkb
                     a(iau+nmat(2))=a(iau+nmat(2))+dvaekb
                     a(iau+nmat(3))=a(iau+nmat(3))+dvackb
                     a(jml+nmat(1))=a(jml+nmat(1))-dvapkb
                     a(jml+nmat(2))=a(jml+nmat(2))-dvaekb
                     a(jml+nmat(3))=a(jml+nmat(3))-dvackb
c
c
c     add air diffusion
c
                     cnvi=cnvf(i)+cnlf(i)
                     dcvi=dcvf(i)+dclf(i)
                     dcvei=dcvef(i)+dclef(i)
                     dcvci=dcvcf(i)+dclcf(i)
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     heatc=t5(neighc)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     heatc=t5(neighc)
                     cnvkb=cnvf(kb)+cnlf(kb)
                     dcvkb=dcvf(kb)+dclf(kb)
                     dcvekb=dcvef(kb)+dclef(kb)
                     dcvckb=dcvcf(kb)+dclcf(kb)
                     bp(iz+nrhs(3))=bp(iz+nrhs(3))+heatc*(cnvkb-cnvi)
                     bp(kz+nrhs(3))=bp(kz+nrhs(3))-heatc*(cnvkb-cnvi)
                     a(jmia+nmat(7))=a(jmia+nmat(7))-heatc*dcvi
                     a(jmia+nmat(8))=a(jmia+nmat(8))-heatc*dcvei
                     a(jmia+nmat(9))=a(jmia+nmat(9))-heatc*dcvci
                     a(ial+nmat(7))=a(ial+nmat(7))+heatc*dcvi
                     a(ial+nmat(8))=a(ial+nmat(8))+heatc*dcvei
                     a(ial+nmat(9))=a(ial+nmat(9))+heatc*dcvci
                     a(iau+nmat(7))=a(iau+nmat(7))+heatc*dcvkb
                     a(iau+nmat(8))=a(iau+nmat(8))+heatc*dcvekb
                     a(iau+nmat(9))=a(iau+nmat(9))+heatc*dcvckb
                     a(jml+nmat(7))=a(jml+nmat(7))-heatc*dcvkb
                     a(jml+nmat(8))=a(jml+nmat(8))-heatc*dcvekb
                     a(jml+nmat(9))=a(jml+nmat(9))-heatc*dcvckb
c
c     add heat conduction
c
                     ti=t(i)
                     jm=1
                     kb=it8(jm)
                     kz=kb-icd
                     neighc=it9(jm)
                     heatc=t5(neighc)
                     iau=it11(jm)
                     ial=it12(jm)
                     jml=nelmdg(kz)-neqp1
                     heatc=t5(neighc)
                     bp(iz+nrhs(2))=bp(iz+nrhs(2))+heatc*(t(kb)-ti)
                     bp(kz+nrhs(2))=bp(kz+nrhs(2))-heatc*(t(kb)-ti)
                     a(jmia+nmat(4))=a(jmia+nmat(4))-heatc*dtpa(i)
                     a(jmia+nmat(5))=a(jmia+nmat(5))-heatc*dtpae(i)
                     a(jmia+nmat(6))=a(jmia+nmat(6))-heatc*dtpac(i)
                     a(ial+nmat(4))=a(ial+nmat(4))+heatc*dtpa(i)
                     a(ial+nmat(5))=a(ial+nmat(5))+heatc*dtpae(i)
                     a(ial+nmat(6))=a(ial+nmat(6))+heatc*dtpac(i)
                     a(iau+nmat(4))=a(iau+nmat(4))+heatc*dtpa(kb)
                     a(iau+nmat(5))=a(iau+nmat(5))+heatc*dtpae(kb)
                     a(iau+nmat(6))=a(iau+nmat(6))+heatc*dtpac(kb)
                     a(jml+nmat(4))=a(jml+nmat(4))-heatc*dtpa(kb)
                     a(jml+nmat(5))=a(jml+nmat(5))-heatc*dtpae(kb)
                     a(jml+nmat(6))=a(jml+nmat(6))-heatc*dtpac(kb)
c     end ifblock for heat/water/air equation type
                  endif
               else
c     heat conduction only(idof=1)
c     start ifblock for heat/water/air equation type
c
c     3-d and 2-d geometry
c
c
c     form equations
c
                  ti=t(i)
                  jm=1
                  kb=it8(jm)
                  kz=kb-icd
                  t5(neighc)=sx_md*max(thx(i),thy(i),thz(i),
     &             thx(kb),thy(kb),thz(kb))
                  neighc=it9(jm)
                  heatc=t5(neighc)
                  iau=it11(jm)
                  ial=it12(jm)
                  jml=nelmdg(kz)-neqp1
                  heatc=t5(neighc)
                  bp(i)=bp(i)+heatc*(t(kb)-ti)
                  bp(kb)=bp(kb)-heatc*(t(kb)-ti)
                  a(jmia)=a(jmia)-heatc*dtpae(i)
                  a(ial)=a(ial)+heatc*dtpae(i)
                  a(iau)=a(iau)+heatc*dtpae(kb)
                  a(jml)=a(jml)-heatc*dtpae(kb)
c     end ifblock for heat conduction only
               endif
c     end ifblock for mdnodes
            endif
c     end do loop on neighbors of md nodes
         enddo
c     end do loop on nodes
      enddo

      if(idpdp.ne.0) then
         do i=1,idof*idof
            nmat(i)=(i-1)*ldna
         enddo
         do i=1,idof
            nrhs(i)=(i-1)*neq
         enddo
      endif

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine gen3dforc(ikforc,ilforc,nforc,nforc_,nodeforc,
     &  ndirforc,xforc,iamforc,ntrans,inotr,trab,rig,ipompc,nodempc,
     &  coefmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,labmpc,iponoel,inoel,
     &  iponoelmax,kon,ipkon,lakon,ne,iponor,xnor,knor,nam,nk,nk_,
     &  co,thicke,nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
     &  iamboun,typeboun,xboun,nmethod,iperturb,istep,vold,mi,idefforc)
!
!     connects nodes of 1-D and 2-D elements, for which
!     concentrated forces were
!     defined, to the nodes of their expanded counterparts
!
      implicit none
!
      logical add,fixed,user,quadratic
!
      character*1 type,typeboun(*)
      character*8 lakon(*)
      character*20 labmpc(*),label
!
      integer ikforc(*),ilforc(*),nodeforc(2,*),ndirforc(*),itr,idirref,
     &  iamforc(*),idim,ier,matz,nodeact,linc,lend,lstart,nnodes,
     &  nforc,nforc_,ntrans,inotr(2,*),rig(*),ipompc(*),nodempc(3,*),
     &  nmpc,nmpc_,mpcfree,ikmpc(*),ilmpc(*),iponoel(*),inoel(3,*),
     &  iponoelmax,kon(*),ipkon(*),ne,iponor(2,*),knor(*),nforcold,
     &  i,node,index,ielem,j,indexe,indexk,nam,iamplitude,idir,
     &  irotnode,nk,nk_,newnode,idof,id,mpcfreenew,k,isector,ndepnodes,
     &  idepnodes(80),l,iexpnode,indexx,irefnode,imax,isol,mpcfreeold,
     &  nod,impc,istep,nodeboun(*),ndirboun(*),ikboun(*),ilboun(*),
     &  nboun,nboun_,iamboun(*),nmethod,iperturb,nrhs,ipiv(3),info,m,
     &  mi(*),idefforc(*),nedge
!
      real*8 xforc(*),trab(7,*),coefmpc(*),xnor(*),val,co(3,*),dot,
     &  thicke(mi(3),*),pi,xboun(*),xnoref(3),dmax,d(3,3),e(3,3,3),
     &  alpha,q(3),w(3),xn(3),a(3,3),a1(3),a2(3),dd,c1,c2,c3,ww,c(3,3),
     &  vold(0:mi(2),*),e1(3),e2(3),t1(3),b(3,3),x(3),y(3),fv1(3),
     &  fv2(3),z(3,3),xi1,xi2,xi3,u(3,3),r(3,3)
!
      data d /1.,0.,0.,0.,1.,0.,0.,0.,1./
      data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
     &        0.,0.,1.,0.,0.,0.,-1.,0.,0.,
     &        0.,-1.,0.,1.,0.,0.,0.,0.,0./
!
      label='                    '
      fixed=.false.
!
      add=.false.
      user=.false.
      pi=4.d0*datan(1.d0)
      isector=0
!
      nforcold=nforc
      do i=1,nforcold
         node=nodeforc(1,i)
         if(node.gt.iponoelmax) then
            if(ndirforc(i).gt.4) then
               write(*,*) '*WARNING: in gen3dforc: node ',i,
     &              ' does not'
               write(*,*) '       belong to a beam nor shell'
               write(*,*) '       element and consequently has no'
               write(*,*) '       rotational degrees of freedom'
            endif
            cycle
         endif
         index=iponoel(node)
         if(index.eq.0) then
            if(ndirforc(i).gt.4) then
               write(*,*) '*WARNING: in gen3dforc: node ',i,
     &              ' does not'
               write(*,*) '       belong to a beam nor shell'
               write(*,*) '       element and consequently has no'
               write(*,*) '       rotational degrees of freedom'
            endif
            cycle
         endif
         ielem=inoel(1,index)
!
!        checking whether element is linear or quardratic
!
         if((lakon(ielem)(4:4).eq.'6').or.
     &      (lakon(ielem)(4:4).eq.'8')) then
            quadratic=.false.
         else
            quadratic=.true.
         endif
!
!        checking whether element is 3-sided or 4-sided
!
         if((lakon(ielem)(4:4).eq.'6').or.
     &      (lakon(ielem)(4:5).eq.'15')) then
            nedge=3
         else
            nedge=4
         endif
!
         j=inoel(2,index)
         indexe=ipkon(ielem)
         indexk=iponor(2,indexe+j)
         if(nam.gt.0) iamplitude=iamforc(i)
         idir=ndirforc(i)
         val=xforc(i)
!
         if(rig(node).ne.0) then
            if(idir.gt.4) then
               if(rig(node).lt.0) then
                  write(*,*) '*ERROR in gen3dforc: in node ',node
                  write(*,*) '       a rotational DOF is loaded;'
                  write(*,*) '       however, the elements to which'
                  write(*,*) '       this node belongs do not have'
                  write(*,*) '       rotational DOFs'
                  call exit(201)
               endif
c               val=xforc(i)
               k=idir-4
               irotnode=rig(node)
               call forcadd(irotnode,k,val,nodeforc,
     &              ndirforc,xforc,nforc,nforc_,iamforc,
     &              iamplitude,nam,ntrans,trab,inotr,co,
     &              ikforc,ilforc,isector,add,user,idefforc,
     &              ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
            endif
         else
!
!           check for moments defined in any but the first step
!           nonlinear dynamic case: creation of knots
!           knots (expandable rigid bodies) can take rotational
!           values arbitrarily exceeding 90 degrees
!
            if((idir.gt.4).and.((nmethod.eq.4).and.(iperturb.gt.1)))then
!
!              create a knot: determine the knot
!
               ndepnodes=0
               if(lakon(ielem)(7:7).eq.'L') then
                  do k=1,3
                     ndepnodes=ndepnodes+1
                     idepnodes(ndepnodes)=knor(indexk+k)
                  enddo
                  idim=1
               elseif(lakon(ielem)(7:7).eq.'B') then
                  do k=1,8
                     ndepnodes=ndepnodes+1
                     idepnodes(ndepnodes)=knor(indexk+k)
                  enddo
                  idim=3
               else
                  write(*,*)
     &           '*ERROR in gen3dboun: a rotational DOF was applied'
                  write(*,*)
     &           '*      to node',node,' without rotational DOFs'
                  call exit(201)
               endif
!
!              remove all MPC's in which the knot nodes are
!              dependent nodes
!
               do k=1,ndepnodes
                  nod=idepnodes(k)
                  do l=1,3
                     idof=8*(nod-1)+l
                     call nident(ikmpc,idof,nmpc,id)
                     if(id.gt.0) then
                        if(ikmpc(id).eq.idof) then
                           impc=ilmpc(id)
                           call mpcrem(impc,mpcfree,nodempc,nmpc,
     &                       ikmpc,ilmpc,labmpc,coefmpc,ipompc)
                        endif
                     endif
                  enddo
               enddo
!
!              generate a rigid body knot
!
               irefnode=node
               nk=nk+1
               if(nk.gt.nk_) then
                  write(*,*) '*ERROR in rigidbodies: increase nk_'
                  call exit(201)
               endif
               irotnode=nk
               rig(node)=irotnode
               nk=nk+1
               if(nk.gt.nk_) then
                  write(*,*) '*ERROR in rigidbodies: increase nk_'
                  call exit(201)
               endif
               iexpnode=nk
               do k=1,ndepnodes
                  call knotmpc(ipompc,nodempc,coefmpc,irefnode,
     &                 irotnode,iexpnode,
     &                 labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,nk,nk_,
     &                 nodeboun,ndirboun,ikboun,ilboun,nboun,nboun_,
     &                 idepnodes,typeboun,co,xboun,istep,k,ndepnodes,
     &                 idim,e1,e2,t1)
               enddo
!
!              determine the location of the center of gravity of
!              the section and its displacements
!
               do l=1,3
                  q(l)=0.d0
                  w(l)=0.d0
               enddo
               if(ndepnodes.eq.3) then
                  do k=1,ndepnodes,2
                     nod=idepnodes(k)
                     do l=1,3
                        q(l)=q(l)+co(l,nod)
                        w(l)=w(l)+vold(l,nod)
                     enddo
                  enddo
                  do l=1,3
                     q(l)=q(l)/2.d0
                     w(l)=w(l)/2.d0
                  enddo
               else
                  do k=1,ndepnodes
                     nod=idepnodes(k)
                     do l=1,3
                        q(l)=q(l)+co(l,nod)
                        w(l)=w(l)+vold(l,nod)
                     enddo
                  enddo
                  do l=1,3
                     q(l)=q(l)/ndepnodes
                     w(l)=w(l)/ndepnodes
                  enddo
               endif
!
!              determine the uniform expansion
!
               alpha=0.d0
               do k=1,ndepnodes
                  nod=idepnodes(k)
                  dd=(co(1,nod)-q(1))**2
     &              +(co(2,nod)-q(2))**2
     &              +(co(3,nod)-q(3))**2
                  if(dd.lt.1.d-20) cycle
                  alpha=alpha+dsqrt(
     &             ((co(1,nod)+vold(1,nod)-q(1)-w(1))**2
     &             +(co(2,nod)+vold(2,nod)-q(2)-w(2))**2
     &             +(co(3,nod)+vold(3,nod)-q(3)-w(3))**2)/dd)
               enddo
               alpha=alpha/ndepnodes
!
!              determine the displacements of irotnodes
!
               do l=1,3
                  do m=1,3
                     a(l,m)=0.d0
                  enddo
                  xn(l)=0.d0
               enddo
               do k=1,ndepnodes
                  nod=idepnodes(k)
                  dd=0.d0
                  do l=1,3
                     a1(l)=co(l,nod)-q(l)
                     a2(l)=vold(l,nod)-w(l)
                     dd=dd+a1(l)*a1(l)
                  enddo
                  dd=dsqrt(dd)
                  if(dd.lt.1.d-10) cycle
                  do l=1,3
                     a1(l)=a1(l)/dd
                     a2(l)=a2(l)/dd
                  enddo
                  xn(1)=xn(1)+(a1(2)*a2(3)-a1(3)*a2(2))
                  xn(2)=xn(2)+(a1(3)*a2(1)-a1(1)*a2(3))
                  xn(3)=xn(3)+(a1(1)*a2(2)-a1(2)*a2(1))
                  do l=1,3
                     do m=1,3
                        a(l,m)=a(l,m)+a1(l)*a1(m)
                     enddo
                  enddo
               enddo
!
               do l=1,3
                  do m=1,3
                     a(l,m)=a(l,m)/ndepnodes
                  enddo
                  xn(l)=xn(l)/ndepnodes
                  a(l,l)=1.d0-a(l,l)
               enddo
!
               m=3
               nrhs=1
               call dgesv(m,nrhs,a,m,ipiv,xn,m,info)
               if(info.ne.0) then
                  write(*,*) '*ERROR in gen3dforc:'
                  write(*,*) '       singular system of equations'
                  call exit(201)
               endif
!
               dd=0.d0
               do l=1,3
                  dd=dd+xn(l)*xn(l)
               enddo
               dd=dsqrt(dd)
               do l=1,3
                  xn(l)=dasin(dd/alpha)*xn(l)/dd
               enddo
!
!              determine the displacements of irefnode
!
               ww=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
!
               c1=dcos(ww)
               if(ww.gt.1.d-10) then
                  c2=dsin(ww)/ww
               else
                  c2=1.d0
               endif
               if(ww.gt.1.d-5) then
                  c3=(1.d0-c1)/ww**2
               else
                  c3=0.5d0
               endif
!
!              rotation matrix r
!
               do k=1,3
                  do l=1,3
                     r(k,l)=c1*d(k,l)+
     &                  c2*(e(k,1,l)*xn(1)+e(k,2,l)*xn(2)+
     &                  e(k,3,l)*xn(3))+c3*xn(k)*xn(l)
                  enddo
               enddo
!
!              copying the displacements
!
               do l=1,3
                  vold(l,irefnode)=w(l)
                  vold(l,irotnode)=xn(l)
               enddo
c               vold(1,iexpnode)=alpha
               vold(1,iexpnode)=alpha-1.d0
!
!              correction of the expansion values for beam sections
!
               if(idim.eq.2) then
!
!              initializing matrices b and c
!
                  do l=1,3
                     do m=1,3
                        b(l,m)=0.d0
                        c(l,m)=0.d0
                     enddo
                  enddo
!
!              solving a least squares problem to determine
!              the transpose of the deformation gradient:
!              c.F^T=b
!
                  do k=1,ndepnodes
                     nod=idepnodes(k)
                     do l=1,3
                        x(l)=co(l,nod)-q(l)
                        y(l)=x(l)+vold(l,nod)-w(l)
                     enddo
                     do l=1,3
                        do m=1,3
                           c(l,m)=c(l,m)+x(l)*x(m)
                           b(l,m)=b(l,m)+x(l)*y(m)
                        enddo
                     enddo
                  enddo
!
!              solving the linear equation system
!
                  m=3
                  nrhs=3
                  call dgesv(m,nrhs,c,m,ipiv,b,m,info)
                  if(info.ne.0) then
                     write(*,*) '*ERROR in gen3dforc:'
                     write(*,*) '       singular system of equations'
                     call exit(201)
                  endif
!
!              now b=F^T
!
!              constructing the right stretch tensor
!              U=F^T.R
!
                  do l=1,3
                     do m=l,3
                        u(l,m)=b(l,1)*r(1,m)+b(l,2)*r(2,m)+
     &                       b(l,3)*r(3,m)
                     enddo
                  enddo
                  u(2,1)=u(1,2)
                  u(3,1)=u(1,3)
                  u(3,2)=u(2,3)
!
!              determining the eigenvalues and eigenvectors of U
!
                  m=3
                  matz=1
                  ier=0
                  call rs(m,m,u,w,matz,z,fv1,fv2,ier)
                  if(ier.ne.0) then
                     write(*,*)
     &                   '*ERROR in knotmpc while calculating the'
                     write(*,*) '       eigenvalues/eigenvectors'
                     call exit(201)
                  endif
!
                  if((dabs(w(1)-1.d0).lt.dabs(w(2)-1.d0)).and.
     &                 (dabs(w(1)-1.d0).lt.dabs(w(3)-1.d0))) then
                     l=2
                     m=3
                  elseif((dabs(w(2)-1.d0).lt.dabs(w(1)-1.d0)).and.
     &                    (dabs(w(2)-1.d0).lt.dabs(w(3)-1.d0))) then
                     l=1
                     m=3
                  else
                     l=1
                     m=2
                  endif
                  xi1=datan2((z(1,l)*e2(1)+z(2,l)*e2(2)+z(3,l)*e2(2)),
     &                 (z(1,l)*e1(1)+z(2,l)*e1(2)+z(3,l)*e1(2)))
                  xi2=w(l)-1.d0
                  xi3=w(m)-1.d0
!
                  vold(1,iexpnode)=xi1
                  vold(2,iexpnode)=xi2
                  vold(3,iexpnode)=xi3
               endif
!
!              apply the moment
!
               idir=idir-4
               val=xforc(i)
               call forcadd(irotnode,idir,val,nodeforc,
     &              ndirforc,xforc,nforc,nforc_,iamforc,
     &              iamplitude,nam,ntrans,trab,inotr,co,
     &              ikforc,ilforc,isector,add,user,idefforc,
     &              ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
!
!              check for shells whether the rotation about the normal
!              on the shell has been eliminated
!
               if(lakon(ielem)(7:7).eq.'L') then
                  indexx=iponor(1,indexe+j)
                  do k=1,3
                     xnoref(k)=xnor(indexx+k)
                  enddo
                  dmax=0.d0
                  imax=0
                  do k=1,3
                     if(dabs(xnoref(k)).gt.dmax) then
                        dmax=dabs(xnoref(k))
                        imax=k
                     endif
                  enddo
!
!                 check whether a SPC suffices
!
                  if(dabs(1.d0-dmax).lt.1.d-3) then
                     val=0.d0
                     if(nam.gt.0) iamplitude=0
                     type='R'
                     call bounadd(irotnode,imax,imax,val,nodeboun,
     &                    ndirboun,xboun,nboun,nboun_,iamboun,
     &                    iamplitude,nam,ipompc,nodempc,coefmpc,
     &                    nmpc,nmpc_,mpcfree,inotr,trab,ntrans,
     &                    ikboun,ilboun,ikmpc,ilmpc,co,nk,nk_,labmpc,
     &                    type,typeboun,nmethod,iperturb,fixed,vold,
     &                    irotnode,mi,label)
                  else
!
!                    check for an unused rotational DOF
!
                     isol=0
                     do l=1,3
                        idof=8*(node-1)+4+imax
                        call nident(ikboun,idof,nboun,id)
                        if((id.gt.0).and.(ikboun(id).eq.idof)) then
                           imax=imax+1
                           if(imax.gt.3) imax=imax-3
                           cycle
                        endif
                        isol=1
                        exit
                     enddo
!
!     if one of the rotational dofs was not used so far,
!     it can be taken as dependent side for fixing the
!     rotation about the normal. If all dofs were used,
!     no additional equation is needed.
!
                     if(isol.eq.1) then
                        idof=8*(irotnode-1)+imax
                        call nident(ikmpc,idof,nmpc,id)
                        nmpc=nmpc+1
                        if(nmpc.gt.nmpc_) then
                           write(*,*)
     &                          '*ERROR in gen3dnor: increase nmpc_'
                           call exit(201)
                        endif
!
                        ipompc(nmpc)=mpcfree
                        labmpc(nmpc)='                    '
!
                        do l=nmpc,id+2,-1
                           ikmpc(l)=ikmpc(l-1)
                           ilmpc(l)=ilmpc(l-1)
                        enddo
                        ikmpc(id+1)=idof
                        ilmpc(id+1)=nmpc
!
                        nodempc(1,mpcfree)=irotnode
                        nodempc(2,mpcfree)=imax
                        coefmpc(mpcfree)=xnoref(imax)
                        mpcfree=nodempc(3,mpcfree)
                        imax=imax+1
                        if(imax.gt.3) imax=imax-3
                        nodempc(1,mpcfree)=irotnode
                        nodempc(2,mpcfree)=imax
                        coefmpc(mpcfree)=xnoref(imax)
                        mpcfree=nodempc(3,mpcfree)
                        imax=imax+1
                        if(imax.gt.3) imax=imax-3
                        nodempc(1,mpcfree)=irotnode
                        nodempc(2,mpcfree)=imax
                        coefmpc(mpcfree)=xnoref(imax)
                        mpcfreeold=mpcfree
                        mpcfree=nodempc(3,mpcfree)
                        nodempc(3,mpcfreeold)=0
                     endif
                  endif
               endif
               cycle
            endif
!
!           all cases except nonlinear dynamic case: creation
!           of mean rotation MPC's
!
            if((idir.gt.4).and.((nmethod.ne.4).or.(iperturb.le.1)))then
!
!              create a mean rotation MPC
!              advantage: more accurate since less constraining
!              disadvantage: cannot exceed 90 degrees rotation
!
!              if a mean rotation MPC has already been created
!              for idof, ilforc(id) contains the index of the
!              SPC created for the angle value
!
               idof=8*(node-1)+idir
               call nident(ikforc,idof,nforc,id)
               if(ilforc(id).ne.i) cycle

               idirref=idir-4
!
               if(lakon(ielem)(7:7).eq.'L') then
                  lstart=3
                  lend=1
                  linc=-2
!
!                 calculating the normal vector =
!                 vector along the drilling direction
!
                  do j=1,3
                     xn(j)=co(j,knor(indexk+3))-co(j,knor(indexk+1))
                  enddo
                  dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
                  do j=1,3
                     xn(j)=xn(j)/dd
                  enddo
!
               elseif(lakon(ielem)(7:7).eq.'B') then
                  lstart=4
                  lend=1
                  linc=-1
               endif
!
!              check for transformations
!
               if(ntrans.le.0) then
                  itr=0
               elseif(inotr(1,node).eq.0) then
                  itr=0
               else
                  itr=inotr(1,node)
               endif
!
!              determine a unit vector on the rotation axis
!
               if(itr.eq.0) then
                  do j=1,3
                     do k=1,3
                        a(j,k)=0.d0
                     enddo
                     a(j,j)=1.d0
                  enddo
               else
                  call transformatrix(trab(1,itr),co(1,node),a)
               endif
!
!              check whether the rotation vector does not
!              have a component along the drilling direction
!
               if(lakon(ielem)(7:7).eq.'L') then
                  dot=a(1,idirref)*xn(1)+a(2,idirref)*xn(2)+
     &                a(3,idirref)*xn(3)
                  if(dot.gt.0.05) then
                     write(*,*) '*ERROR in gen3dforc: applied'
                     write(*,*) '       moment in node ',node
                     write(*,*) '       and direction ',idir-1
                     write(*,*) '       has a significant'
                     write(*,*)
     &                    '       component along the drilling'
                     write(*,*) '       direction; this is not'
                     write(*,*) '       allowed'
                     call exit(201)
                  endif
               endif
!
!              specific label for mean rotations for beams and
!              shells
!
               label='MEANROTBS           '
               nnodes=0
               do j=lstart,lend,linc
                  nodeact=knor(indexk+j)
                  do k=1,3
                     nnodes=nnodes+1
                     call usermpc(ipompc,nodempc,coefmpc,
     &                    labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
     &                    nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
     &                    nboun,nboun_,nnodes,nodeact,co,label,
     &                    typeboun,iperturb,node,idirref,xboun)
                  enddo
               enddo
!
!              rotation value term
!
               nodeact=nk+1
               do k=1,3
                  co(k,nodeact)=a(k,idirref)
               enddo
               nnodes=nnodes+1
               call usermpc(ipompc,nodempc,coefmpc,
     &              labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
     &              nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
     &              nboun,nboun_,nnodes,nodeact,co,label,
     &              typeboun,iperturb,node,idirref,xboun)
!
!              inhomogeneous term
!
               nodeact=0
               call usermpc(ipompc,nodempc,coefmpc,
     &              labmpc,nmpc,nmpc_,mpcfree,ikmpc,ilmpc,
     &              nk,nk_,nodeboun,ndirboun,ikboun,ilboun,
     &              nboun,nboun_,nnodes,nodeact,co,label,
     &              typeboun,iperturb,node,idirref,xboun)
!
!              end meanrotationmpc
!
!              SPC angle term
!
               if(nodeact.ne.-1) then
                  idir=1
                  call forcadd(nk,idir,val,nodeforc,
     &              ndirforc,xforc,nforc,nforc_,iamforc,
     &              iamplitude,nam,ntrans,trab,inotr,co,
     &              ikforc,ilforc,isector,add,user,idefforc,
     &              ipompc,nodempc,nmpc,ikmpc,ilmpc,labmpc)
!
!                 storing the index of the SPC with the angle
!                 value in ilboun(id)
!
                  ilforc(id)=nforc
               endif
!
               cycle
            endif
!
!           2d shell element: generate MPC's
!
            if(lakon(ielem)(7:7).eq.'L') then
               newnode=knor(indexk+1)
               idof=8*(newnode-1)+idir
               call nident(ikmpc,idof,nmpc,id)
               if((id.le.0).or.(ikmpc(id).ne.idof)) then
                  nmpc=nmpc+1
                  if(nmpc.gt.nmpc_) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
                  labmpc(nmpc)='                    '
                  ipompc(nmpc)=mpcfree
                  do k=nmpc,id+2,-1
                     ikmpc(k)=ikmpc(k-1)
                     ilmpc(k)=ilmpc(k-1)
                  enddo
                  ikmpc(id+1)=idof
                  ilmpc(id+1)=nmpc
!
!                 for middle nodes: u_1+u_3-2*u_node=0
!                 for end nodes: -u_1+4*u_2-u_3-2*u_node=0
!
!                 u_1 corresponds to knor(indexk+1)....
!
                  nodempc(1,mpcfree)=newnode
                  nodempc(2,mpcfree)=idir
                  if((j.gt.nedge).or.(.not.quadratic)) then
                     coefmpc(mpcfree)=1.d0
                  else
                     coefmpc(mpcfree)=-1.d0
                  endif
                  mpcfree=nodempc(3,mpcfree)
                  if(mpcfree.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
!
                  if((j.le.nedge).and.(quadratic)) then
                     nodempc(1,mpcfree)=knor(indexk+2)
                     nodempc(2,mpcfree)=idir
                     coefmpc(mpcfree)=4.d0
                     mpcfree=nodempc(3,mpcfree)
                     if(mpcfree.eq.0) then
                        write(*,*)
     &                       '*ERROR in gen3dforc: increase nmpc_'
                        call exit(201)
                     endif
                  endif
!
                  nodempc(1,mpcfree)=knor(indexk+3)
                  nodempc(2,mpcfree)=idir
                  if((j.gt.nedge).or.(.not.quadratic)) then
                     coefmpc(mpcfree)=1.d0
                  else
                     coefmpc(mpcfree)=-1.d0
                  endif
                  mpcfree=nodempc(3,mpcfree)
                  if(mpcfree.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
!
                  nodempc(1,mpcfree)=node
                  nodempc(2,mpcfree)=idir
                  coefmpc(mpcfree)=-2.d0
                  mpcfreenew=nodempc(3,mpcfree)
                  if(mpcfreenew.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
!
                  nodempc(3,mpcfree)=0
                  mpcfree=mpcfreenew
               endif
            elseif(lakon(ielem)(7:7).eq.'B') then
!
!                       1d beam element: generate MPC's
!
               newnode=knor(indexk+1)
               idof=8*(newnode-1)+idir
               call nident(ikmpc,idof,nmpc,id)
               if((id.le.0).or.(ikmpc(id).ne.idof)) then
                  nmpc=nmpc+1
                  if(nmpc.gt.nmpc_) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
                  labmpc(nmpc)='                    '
                  ipompc(nmpc)=mpcfree
                  do k=nmpc,id+2,-1
                     ikmpc(k)=ikmpc(k-1)
                     ilmpc(k)=ilmpc(k-1)
                  enddo
                  ikmpc(id+1)=idof
                  ilmpc(id+1)=nmpc
                     do k=1,4
                        nodempc(1,mpcfree)=knor(indexk+k)
                        nodempc(2,mpcfree)=idir
                        coefmpc(mpcfree)=1.d0
                        mpcfree=nodempc(3,mpcfree)
                        if(mpcfree.eq.0) then
                           write(*,*)
     &                          '*ERROR in gen3dforc: increase nmpc_'
                           call exit(201)
                        endif
                     enddo
                  nodempc(1,mpcfree)=node
                  nodempc(2,mpcfree)=idir
                  coefmpc(mpcfree)=-4.d0
                  mpcfreenew=nodempc(3,mpcfree)
                  if(mpcfreenew.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dforc: increase nmpc_'
                     call exit(201)
                  endif
                  nodempc(3,mpcfree)=0
                  mpcfree=mpcfreenew
               endif
            else
!
!              2d plane stress, plane strain or axisymmetric
!              element: MPC in all but z-direction
!
               newnode=knor(indexk+2)
               idof=8*(newnode-1)+idir
               call nident(ikmpc,idof,nmpc,id)
               if(((id.le.0).or.(ikmpc(id).ne.idof)).and.
     &              (idir.ne.3)) then
                  nmpc=nmpc+1
                  if(nmpc.gt.nmpc_) then
                     write(*,*)
     &                    '*ERROR in gen3dmpc: increase nmpc_'
                     call exit(201)
                  endif
                  labmpc(nmpc)='                    '
                  ipompc(nmpc)=mpcfree
                  do j=nmpc,id+2,-1
                     ikmpc(j)=ikmpc(j-1)
                     ilmpc(j)=ilmpc(j-1)
                  enddo
                  ikmpc(id+1)=idof
                  ilmpc(id+1)=nmpc
                  nodempc(1,mpcfree)=newnode
                  nodempc(2,mpcfree)=idir
                  coefmpc(mpcfree)=1.d0
                  mpcfree=nodempc(3,mpcfree)
                  if(mpcfree.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dmpc: increase nmpc_'
                     call exit(201)
                  endif
                  nodempc(1,mpcfree)=node
                  nodempc(2,mpcfree)=idir
                  coefmpc(mpcfree)=-1.d0
                  mpcfreenew=nodempc(3,mpcfree)
                  if(mpcfreenew.eq.0) then
                     write(*,*)
     &                    '*ERROR in gen3dmpc: increase nmpc_'
                     call exit(201)
                  endif
                  nodempc(3,mpcfree)=0
                  mpcfree=mpcfreenew
               endif
            endif
         endif
      enddo
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
program ch0904
    ! Imperial Pints and Litres
    implicit none

    integer :: fluid
    real :: litres
    real :: pints

    print *, ' Imperial      Litre(s)'
    print *, ' pint(s)               '
    do fluid = 1, 10
        litres = fluid/1.75
        pints = fluid*1.75
        print 100, pints, fluid, litres
    end do

100 format (' ', f6.2, '  ', i3, '   ', f5.2)
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
c-----7---------------------------------------------------------------72
c
c                 CMORPH 0.25 degree 3-hourly Information
c
c                    Last Update: February 13 2007
c
c The 0.25 degree 3-hourly CMORPH can be found on the ftp.cpc.ncep.noaa.gov
c server in the precip/global_CMORPH/3-hourly_025deg directory
c
c
c The data are compressed using the standard Unix compress function (files have a
c suffix of ".Z"). Each file is composed of 16 direct access binary
c ("big_endian") records that are defined as follows:
c
c Record 1: contains the merged microwave precipitation only for 00 UTC
c Record 2: contains the "CMORPH" precipitation estimates for 00 UTC
c
c Record 3: contains the merged microwave precipitation only for 03 UTC
c Record 4: contains the "CMORPH" precipitation estimates for 03 UTC
c .
c .
c .
c Record 15: contains the merged microwave precipitation only for 21 UTC
c Record 16: contains the "CMORPH" precipitation estimates for 21 UTC
c
c All units are "mm/hr".  Missing data are denoted by values of "-9999."
c Each record contains a 1440 x 480 REAL*4 array of data which is oriented
c from 0.125E EASTward and from 59.875N SOUTHward, with a grid increment
c of 0.25 degrees of latitude and longitude.  Thus, the grid locations
c are the centers on a 0.25 degree lat/lon grid box.  Note that these
c estimates represent spatial averages, so the data are grid-centered,
c rather than lattice-centered.
c
c For example (1,1) is 0.125E, 59.875N
c             (2,2) is 0.375E, 59.625N, etc.
c
c There is an associated GrADS ctl file with this data set:
c CMORPH+MWCOMB_025deg-3hr.ctl
c
c
c Below is FORTRAN code (run on an SGI) that will read the data sets:
c
c-----7---------------------------------------------------------------72

      program cmorph_025deg

      parameter (IRMAX=1440, JRMAX=480, PMAX=2)

      real*4 ssmipropREG(IRMAX,JRMAX,PMAX)
      character*115 in1ssmipropREG

      iyr = 2011
      iyrmndy = 20110901
      irec = 0

c    Open input regional blended sensor rainfall file for 3
c      hour CPC micro rainfall file

      in1ssmipropREG = '/glade/data02/dsswork/tcram/tigge/cmorph/data/'
     &//'025deg/YYYY/YYYY0M0D_3hr-025deg_cpc+comb'

      write (in1ssmipropREG(54:57),'(i4)') iyr
      write (in1ssmipropREG(59:66),'(i8)') iyrmndy
      write(6,'('' file '',a88)') in1ssmipropREG

      open (unit=15, file=in1ssmipropREG, access='direct',
     &              recl=1440*480*2*4, form='unformatted',
     &              convert='big_endian')

c    Loop over 8 3 hourly periods for one day
      do ihr = 0, 21, 3

       irec = irec + 1

       write(6,'(//,''proc 3 hr starting '',i8,i4)') iyrmndy, ihr

       read (15, rec = irec) ssmipropREG

       write (6,'(/,'' cpc combined MW '')')
       do jreg = 271, 280
        write(6,'(15(f12.1))') (ssmipropREG(ireg,jreg,1),ireg=1161,1175)
       enddo

       write (6,'(/,'' CMORPH '')')
       do jreg = 271, 280
        write(6,'(15(f12.1))') (ssmipropREG(ireg,jreg,2),ireg=1161,1175)
       enddo

c    End 3 hour loop

      enddo

      close (15)

      stop

      end

max_stars_repo_name: NCAR/rda-cmorph

max_issues_repo_name: NCAR/rda-cmorph

max_forks_repo_name: NCAR/rda-cmorph


text:
Im a librarian at the http://www.skokielibrary.info/ Skokie Public Library.  I helped create this page with my colleague Users/tgreenwalt.

20070822 16:45:13 nbsp Hi, perhaps you came from the wikispot hub? This is the wiki for the Town of Davis. Your library appears to be in Illinois.

Regardless, welcome to the wiki(s) Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!**********************************************************************************************************************************
!
!  MODULE: Orca_Driver_Types  - This module contains types used by the OrcaFlexInterface Driver program to store arguments passed in
!
!  The types listed here are used within the OrcaFlexInterface Driver program to store the settings. These settings are read in as
!  command line arguments, then stored within these types.
!
!**********************************************************************************************************************************
!
!..................................................................................................................................
! LICENSING
! Copyright (C) 2015  National Renewable Energy Laboratory
!
!    This file is part of OrcaFlexInterface.
!
!    OrcaFlexInterface is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as
!    published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
!    of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License along with OrcaFlexInterface.
!    If not, see <http://www.gnu.org/licenses/>.
!
!**********************************************************************************************************************************
! File last committed: $Date: 2014-07-29 13:30:04 -0600 (Tue, 29 Jul 2014) $
! (File) Revision #: $Rev: 169 $
! URL: $HeadURL: https://windsvn.nrel.gov/OrcaFlexInterface/branches/modularization2/Source/Driver/OrcaDriver_Types.f90 $
!**********************************************************************************************************************************

MODULE OrcaDriver_Types

   USE NWTC_Library
   USE OrcaFlexInterface_Types

   IMPLICIT NONE

      !> This contains flags to note if the settings were made.  This same data structure is
      !! used both during the driver input file and the command line options.
      !!
      !! NOTE: The WindFileType is only set if it is given as a command line option.  Otherwise
      !!       it is handled internally by InflowWInd.
      !!
      !! NOTE: The wind direction is specified by the OrcaFlexInterface input file.
   TYPE     :: OrcaDriver_Flags
      LOGICAL                 :: DvrIptFile           = .FALSE.      !< Was an input file name given on the command line?
      LOGICAL                 :: OrcaIptFile          = .FALSE.      !< Was an OrcaFlexInterface input file requested?
      LOGICAL                 :: AddedMass            = .FALSE.      !< create an added mass table at command line?
      LOGICAL                 :: AddedMassFile        = .FALSE.      !< create an added mass file?
      LOGICAL                 :: DT                   = .FALSE.      !< specified a resolution in time
      LOGICAL                 :: DTDefault            = .FALSE.      !< specified a 'DEFAULT' for the time resolution

      LOGICAL                 :: Degrees              = .FALSE.      !< angles are specified in degrees

      LOGICAL                 :: PtfmCoord            = .FALSE.      !< (x,y,z,R1,R2,R3) coordinate specified
      LOGICAL                 :: PtfmVeloc            = .FALSE.      !< (x,y,z,R1,R2,R3) coordinate specified
      LOGICAL                 :: PtfmAccel            = .FALSE.      !< (x,y,z,R1,R2,R3) coordinate specified

      LOGICAL                 :: PointsFile           = .FALSE.      !< points filename to read in
      LOGICAL                 :: PointsDegrees        = .FALSE.      !< points in the pointsfile are specified in degrees

      LOGICAL                 :: AddedMassOutputInit  = .FALSE.      !< Is the WindGridOut file initialized
      LOGICAL                 :: PointsOutputInit     = .FALSE.      !< Is the Points output file initialized
      LOGICAL                 :: Verbose              = .FALSE.      !< Verbose error reporting
      LOGICAL                 :: VVerbose             = .FALSE.      !< Very Verbose error reporting
   END TYPE    OrcaDriver_Flags

      ! This contains all the settings (possible passed in arguments).
   TYPE     :: OrcaDriver_Settings
      CHARACTER(1024)         :: DvrIptFileName                !< Driver input file name
      CHARACTER(1024)         :: OrcaIptFileName               !< Filename of OrcaFlexInterface input file to read (if no driver input file)
      CHARACTER(1024)         :: AddedMassFileName             !< Filename for the added mass matrix output

      CHARACTER(1024)         :: PointsFileName                !< Filename of points file to read in
      CHARACTER(1024)         :: PointsOutputName              !< Filename for output from points read in from points file
      INTEGER(IntKi)          :: AddedMassOutputUnit           !< Unit number for the output file for the AddedMass matrix
      INTEGER(IntKi)          :: PointsOutputUnit              !< Unit number for the output file for the Points file output
      REAL(DbKi)              :: DT                            !< resolution of time
      REAL(ReKi)              :: TMax                          !< Maximum time (we calculate this based on the number of points and timestep)

      REAL(ReKi)              :: PtfmCoord(1:6)                !< (x,y,z,R1,R2,R3) coordinate and rotations to calculate at
      REAL(ReKi)              :: PtfmVeloc(1:6)                !< instantaneous velocities corresponding to the PtfmCoord
      REAL(ReKi)              :: PtfmAccel(1:6)                !< instantaneous velocities corresponding to the PtfmCoord

      TYPE(ProgDesc)          :: ProgInfo                      !< Program info
      TYPE(ProgDesc)          :: OrcaProgInfo                  !< Program info for OrcaFlexInterface

   END TYPE    OrcaDriver_Settings

END MODULE OrcaDriver_Types

max_stars_repo_name: abhineet-gupta/openfast

max_issues_repo_name: abhineet-gupta/openfast

max_forks_repo_name: abhineet-gupta/openfast


text:
named_block : BLOCK
  INTEGER :: b
  a = 3
END BLOCK named_block
BLOCK
  c = 13
END BLOCK
BLOCK
END BLOCK
END PROGRAM

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf