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 printoutint(prlab,ipkon,lakon,stx,eei,xstate,ener,
     &  mi,nstate_,ii,nelem,qfx,orab,ielorien,norien,co,konf,
     &  ielmat,thicke,eme,ielprop,prop,nelel,ithermal,orname)
!
!     stores integration point results for element "nelem" in the .dat file
!
!     nelem is the element number from the input deck
!     nelel is a local number, created e.g. by renumbering; thus
!           far only applicable for CFD
!
      implicit none
!
      character*6 prlab(*)
      character*8 lakon(*)
      character*80 orname(*)
!
      integer ipkon(*),mi(*),nstate_,nelem,l,ii,mint3d,j,k,nope,
     &  ielorien(mi(3),*),norien,konf(*),konl,indexe,m,iorien,iflag,
     &  ielmat(mi(3),*),nopes,mint2d,kk,ki,kl,nlayer,ilayer,
     &  null,ielprop(*),nelel,ithermal(2)
!
      real*8 stx(6,mi(1),*),eei(6,mi(1),*),xstate(nstate_,mi(1),*),
     &  ener(mi(1),*),qfx(3,mi(1),*),xi,et,ze,xl(3,20),xsj,shp(4,20),
     &  coords(3,27),weight,orab(7,*),co(3,*),a(3,3),b(3,3),c(3,3),
     &  qfxl(3),thicke(mi(3),*),xsj2(3),shp2(7,8),xl2(3,8),xs2(3,7),
     &  thickness,tlayer(4),dlayer(4),xlayer(mi(3),4),eme(6,mi(1),*),
     &  prop(*)
!
      include "gauss.f"
!
      data iflag /1/
      data null /0/
!
      if(ipkon(nelel).lt.0) then
         return
      else
         indexe=ipkon(nelel)
      endif
!
!     check whether transformation is necessary (if orientation
!     is applied and output in local system is requested)
!
      if(lakon(nelel)(7:8).ne.'LC') then
         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
            iorien=0
         else
            iorien=ielorien(1,nelel)
         endif
      elseif(lakon(nelel)(4:5).eq.'20') then
!
!     composite materials
!
         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
            iorien=0
         else
!
!           check whether at least one layer has a transformation
!
            iorien=0
            do k=1,mi(3)
               if(ielorien(k,nelel).ne.0) then
                  iorien=ielorien(k,nelel)
                  exit
               endif
            enddo
         endif
!
!     determining the number of layers
!
         mint2d=4
         nopes=8
!
         nlayer=0
         do k=1,mi(3)
            if(ielmat(k,nelel).ne.0) then
               nlayer=nlayer+1
            endif
         enddo
!
!     determining the layer thickness and global thickness
!     at the shell integration points
!
         iflag=1
         do kk=1,mint2d
            xi=gauss3d2(1,kk)
            et=gauss3d2(2,kk)
            call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
            tlayer(kk)=0.d0
            do k=1,nlayer
               thickness=0.d0
               do j=1,nopes
                  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
               enddo
               tlayer(kk)=tlayer(kk)+thickness
               xlayer(k,kk)=thickness
            enddo
         enddo
         iflag=3
!
         ilayer=0
         do k=1,4
            dlayer(k)=0.d0
         enddo
!
!
      elseif(lakon(nelel)(4:5).eq.'15') then
!
!     composite materials
!
         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
            iorien=0
         else
!
!           check whether at least one layer has a transformation
!
            iorien=0
            do k=1,mi(3)
               if(ielorien(k,nelel).ne.0) then
                  iorien=ielorien(k,nelel)
                  exit
               endif
            enddo
         endif
!
!     determining the number of layers
!
         mint2d=3
         nopes=6
!
         nlayer=0
         do k=1,mi(3)
            if(ielmat(k,nelel).ne.0) then
               nlayer=nlayer+1
            endif
         enddo
!
!     determining the layer thickness and global thickness
!     at the shell integration points
!
         iflag=1
         do kk=1,mint2d
            xi=gauss3d10(1,kk)
            et=gauss3d10(2,kk)
            call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
            tlayer(kk)=0.d0
            do k=1,nlayer
               thickness=0.d0
               do j=1,nopes
                  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
               enddo
               tlayer(kk)=tlayer(kk)+thickness
               xlayer(k,kk)=thickness
            enddo
         enddo
         iflag=3
!
         ilayer=0
         do k=1,3
            dlayer(k)=0.d0
         enddo
!
      endif
!
!     number of integration points
!
      if((lakon(nelel)(4:5).eq.'8R').or.
     &   (lakon(nelel)(1:1).eq.'F')) then
         mint3d=1
      elseif(lakon(nelel)(4:7).eq.'20RB') then
         if((lakon(nelel)(8:8).eq.'R').or.
     &      (lakon(nelel)(8:8).eq.'C')) then
            mint3d=50
         else
            call beamintscheme(lakon(nelel),mint3d,ielprop(nelel),prop,
     &           null,xi,et,ze,weight)
         endif
      elseif((lakon(nelel)(4:4).eq.'8').or.
     &       (lakon(nelel)(4:6).eq.'20R')) then
         if(lakon(nelel)(7:8).eq.'LC') then
            mint3d=8*nlayer
         else
            mint3d=8
         endif
      elseif(lakon(nelel)(4:4).eq.'2') then
         mint3d=27
      elseif(lakon(nelel)(4:5).eq.'10') then
         mint3d=4
      elseif(lakon(nelel)(4:4).eq.'4') then
         mint3d=1
      elseif(lakon(nelel)(4:5).eq.'15') then
         if(lakon(nelel)(7:8).eq.'LC') then
            mint3d=6*nlayer
         else
            mint3d=9
         endif
      elseif(lakon(nelel)(4:4).eq.'6') then
         mint3d=2
      else
         return
      endif
!
!     calculation of the integration point coordinates for
!     output in the local system (if needed)
!
      if(iorien.ne.0) then
         if(lakon(nelel)(4:4).eq.'2') then
            nope=20
         elseif(lakon(nelel)(4:4).eq.'8') then
            nope=8
         elseif(lakon(nelel)(4:5).eq.'10') then
            nope=10
         elseif(lakon(nelel)(4:4).eq.'4') then
            nope=4
         elseif(lakon(nelel)(4:5).eq.'15') then
            nope=15
         elseif(lakon(nelel)(4:4).eq.'6') then
            nope=6
         endif
!
c         indexe=ipkon(nelel)
         do j=1,nope
            konl=konf(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl)
            enddo
         enddo
!
         do j=1,mint3d
            if((lakon(nelel)(4:5).eq.'8R').or.
     &         (lakon(nelel)(1:4).eq.'F3D8')) then
               xi=gauss3d1(1,j)
               et=gauss3d1(2,j)
               ze=gauss3d1(3,j)
               weight=weight3d1(j)
            elseif(lakon(nelel)(4:8).eq.'20RB') then
               if((lakon(nelel)(8:8).eq.'B').or.
     &            (lakon(nelel)(8:8).eq.'C')) then
                  xi=gauss3d13(1,j)
                  et=gauss3d13(2,j)
                  ze=gauss3d13(3,j)
                  weight=weight3d13(j)
               else
                  call beamintscheme(lakon(nelel),mint3d,ielprop(nelel),
     &                  prop,j,xi,et,ze,weight)
               endif
            elseif((lakon(nelel)(4:4).eq.'8').or.
     &             (lakon(nelel)(4:6).eq.'20R'))
     &              then
               if(lakon(nelel)(7:8).ne.'LC') then
                  xi=gauss3d2(1,j)
                  et=gauss3d2(2,j)
                  ze=gauss3d2(3,j)
                  weight=weight3d2(j)
               else
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
!
                  xi=gauss3d2(1,kl)
                  et=gauss3d2(2,kl)
                  ze=gauss3d2(3,kl)
                  weight=weight3d2(kl)
!
                  ki=mod(j,4)
                  if(ki.eq.0) ki=4
!
                  if(kl.eq.1) then
                     ilayer=ilayer+1
                     if(ilayer.gt.1) then
                        do k=1,4
                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
                        enddo
                     endif
                  endif
                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
     &                 tlayer(ki)-1.d0
                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
               endif
            elseif(lakon(nelel)(4:4).eq.'2') then
               xi=gauss3d3(1,j)
               et=gauss3d3(2,j)
               ze=gauss3d3(3,j)
               weight=weight3d3(j)
            elseif(lakon(nelel)(4:5).eq.'10') then
               xi=gauss3d5(1,j)
               et=gauss3d5(2,j)
               ze=gauss3d5(3,j)
               weight=weight3d5(j)
            elseif(lakon(nelel)(4:4).eq.'4') then
               xi=gauss3d4(1,j)
               et=gauss3d4(2,j)
               ze=gauss3d4(3,j)
               weight=weight3d4(j)
            elseif(lakon(nelel)(4:5).eq.'15') then
               if(lakon(nelel)(7:8).ne.'LC') then
                  xi=gauss3d8(1,j)
                  et=gauss3d8(2,j)
                  ze=gauss3d8(3,j)
                  weight=weight3d8(j)
               else
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
!
                  xi=gauss3d10(1,kl)
                  et=gauss3d10(2,kl)
                  ze=gauss3d10(3,kl)
                  weight=weight3d10(kl)
!
                  ki=mod(j,3)
                  if(ki.eq.0) ki=3
!
                  if(kl.eq.1) then
                     ilayer=ilayer+1
                     if(ilayer.gt.1) then
                        do k=1,3
                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
                        enddo
                     endif
                  endif
                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
     &                 tlayer(ki)-1.d0
                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
               endif
            elseif(lakon(nelel)(1:4).eq.'C3D6') then
               xi=gauss3d7(1,j)
               et=gauss3d7(2,j)
               ze=gauss3d7(3,j)
               weight=weight3d7(j)
            elseif(lakon(nelel)(1:4).eq.'F3D6') then
               xi=gauss3d14(1,j)
               et=gauss3d14(2,j)
               ze=gauss3d14(3,j)
               weight=weight3d14(j)
            endif
!
            if(nope.eq.20) then
               call shape20h(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
            else
               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
            endif
!
            do k=1,3
               coords(k,j)=0.d0
               do l=1,nope
                  coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l)
               enddo
            enddo
         enddo
      endif
!
      if((prlab(ii)(1:4).eq.'S   ').or.(prlab(ii)(1:4).eq.'SVF ')) then
         do j=1,mint3d
!
!           composite materials
!
            if(lakon(nelel)(7:8).eq.'LC') then
               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
                  iorien=0
               elseif(lakon(nelel)(4:5).eq.'20') then
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
                  ilayer=(j-kl)/8+1
                  iorien=ielorien(ilayer,nelel)
               elseif(lakon(nelel)(4:5).eq.'15') then
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
                  ilayer=(j-kl)/6+1
                  iorien=ielorien(ilayer,nelel)
               endif
            endif
!
            if(iorien.eq.0) then
               write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
     &              (stx(k,j,nelel),k=1,6)
            else
               call transformatrix(orab(1,iorien),coords(1,j),a)
               b(1,1)=stx(1,j,nelel)
               b(2,2)=stx(2,j,nelel)
               b(3,3)=stx(3,j,nelel)
               b(1,2)=stx(4,j,nelel)
               b(1,3)=stx(5,j,nelel)
               b(2,3)=stx(6,j,nelel)
               b(2,1)=b(1,2)
               b(3,1)=b(1,3)
               b(3,2)=b(2,3)
               do k=1,3
                  do l=1,3
                     c(k,l)=0.d0
                     do m=1,3
                        c(k,l)=c(k,l)+b(k,m)*a(m,l)
                     enddo
                  enddo
               enddo
               do k=1,3
                  do l=k,3
                     b(k,l)=0.d0
                     do m=1,3
                        b(k,l)=b(k,l)+a(m,k)*c(m,l)
                     enddo
                  enddo
               enddo
               write(5,'(i10,1x,i3,1p,6(1x,e13.6),1x,a20)') nelem,j,
     &              b(1,1),b(2,2),b(3,3),b(1,2),b(1,3),b(2,3),
     &              orname(iorien)(1:20)
            endif
         enddo
      elseif(prlab(ii)(1:4).eq.'E   ') then
         do j=1,mint3d
!
!           composite materials
!
            if(lakon(nelel)(7:8).eq.'LC') then
               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
                  iorien=0
               elseif(lakon(nelel)(4:5).eq.'20') then
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
                  ilayer=(j-kl)/8+1
                  iorien=ielorien(ilayer,nelel)
               elseif(lakon(nelel)(4:5).eq.'15') then
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
                  ilayer=(j-kl)/6+1
                  iorien=ielorien(ilayer,nelel)
               endif
            endif
!
            if(iorien.eq.0) then
               write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
     &              (eei(k,j,nelel),k=1,6)
            else
               call transformatrix(orab(1,iorien),coords(1,j),a)
               b(1,1)=eei(1,j,nelel)
               b(2,2)=eei(2,j,nelel)
               b(3,3)=eei(3,j,nelel)
               b(1,2)=eei(4,j,nelel)
               b(1,3)=eei(5,j,nelel)
               b(2,3)=eei(6,j,nelel)
               b(2,1)=b(1,2)
               b(3,1)=b(1,3)
               b(3,2)=b(2,3)
               do k=1,3
                  do l=1,3
                     c(k,l)=0.d0
                     do m=1,3
                        c(k,l)=b(k,m)*a(m,l)
                     enddo
                  enddo
               enddo
               do k=1,3
                  do l=k,3
                     b(k,l)=0.d0
                     do m=1,3
                        b(k,l)=a(m,k)*c(m,l)
                     enddo
                  enddo
               enddo
               write(5,'(i10,1x,i3,1p,6(1x,e13.6),1x,a20)') nelem,j,
     &              b(1,1),b(2,2),b(3,3),b(1,2),b(1,3),b(2,3),
     &              orname(iorien)(1:20)
            endif
         enddo
      elseif(prlab(ii)(1:4).eq.'ME  ') then
         do j=1,mint3d
!
!           composite materials
!
            if(lakon(nelel)(7:8).eq.'LC') then
               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
                  iorien=0
               elseif(lakon(nelel)(4:5).eq.'20') then
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
                  ilayer=(j-kl)/8+1
                  iorien=ielorien(ilayer,nelel)
               elseif(lakon(nelel)(4:5).eq.'15') then
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
                  ilayer=(j-kl)/6+1
                  iorien=ielorien(ilayer,nelel)
               endif
            endif
!
            if(iorien.eq.0) then
               write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
     &              (eme(k,j,nelel),k=1,6)
            else
               call transformatrix(orab(1,iorien),coords(1,j),a)
               b(1,1)=eme(1,j,nelel)
               b(2,2)=eme(2,j,nelel)
               b(3,3)=eme(3,j,nelel)
               b(1,2)=eme(4,j,nelel)
               b(1,3)=eme(5,j,nelel)
               b(2,3)=eme(6,j,nelel)
               b(2,1)=b(1,2)
               b(3,1)=b(1,3)
               b(3,2)=b(2,3)
               do k=1,3
                  do l=1,3
                     do m=1,3
                        c(k,l)=b(k,m)*a(m,l)
                     enddo
                  enddo
               enddo
               do k=1,3
                  do l=k,3
                     do m=1,3
                        b(k,l)=a(m,k)*c(m,l)
                     enddo
                  enddo
               enddo
               write(5,'(i10,1x,i3,1p,6(1x,e13.6),1x,a20)') nelem,j,
     &              b(1,1),b(2,2),b(3,3),b(1,2),b(1,3),b(2,3),
     &              orname(iorien)(1:20)
            endif
         enddo
      elseif(prlab(ii)(1:4).eq.'PEEQ') then
         do j=1,mint3d
            write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
     &           xstate(1,j,nelel)
         enddo
      elseif(prlab(ii)(1:4).eq.'ENER') then
         do j=1,mint3d
            write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
     &           ener(j,nelel)
         enddo
      elseif(prlab(ii)(1:4).eq.'SDV ') then
         do j=1,mint3d
!
!           composite materials
!
            if(lakon(nelel)(7:8).eq.'LC') then
               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
                  iorien=0
               elseif(lakon(nelel)(4:5).eq.'20') then
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
                  ilayer=(j-kl)/8+1
                  iorien=ielorien(ilayer,nelel)
               elseif(lakon(nelel)(4:5).eq.'15') then
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
                  ilayer=(j-kl)/6+1
                  iorien=ielorien(ilayer,nelel)
               endif
            endif
!
            if(iorien.ne.0) then
               write(*,*) '*WARNING in printoutint: SDV cannot be'
               write(*,*) '         printed in the local system'
               write(*,*) '         results are in the global system'
            endif
            write(5,'(i10,1x,i3,1p,99(1x,e13.6))') nelem,j,
     &           (xstate(k,j,nelel),k=1,nstate_)
         enddo
      elseif(((prlab(ii)(1:4).eq.'HFL ').or.(prlab(ii)(1:4).eq.'HFLF'))
     &        .and.(ithermal(1).gt.1)) then
         do j=1,mint3d
!
!           composite materials
!
            if(lakon(nelel)(7:8).eq.'LC') then
               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
                  iorien=0
               elseif(lakon(nelel)(4:5).eq.'20') then
                  kl=mod(j,8)
                  if(kl.eq.0) kl=8
                  ilayer=(j-kl)/8+1
                  iorien=ielorien(ilayer,nelel)
               elseif(lakon(nelel)(4:5).eq.'15') then
                  kl=mod(j,6)
                  if(kl.eq.0) kl=6
                  ilayer=(j-kl)/6+1
                  iorien=ielorien(ilayer,nelel)
               endif
            endif
!
            if(iorien.eq.0) then
               write(5,'(i10,1x,i3,1p,3(1x,e13.6))') nelem,j,
     &              (qfx(k,j,nelel),k=1,3)
            else
               do k=1,3
                  qfxl(k)=qfx(k,j,nelel)
               enddo
               call transformatrix(orab(1,iorien),coords(1,j),a)
               write(5,'(i10,1x,i3,1p,3(1x,e13.6),1x,a20)') nelem,j,
     &              qfxl(1)*a(1,1)+qfxl(2)*a(2,1)+qfxl(3)*a(3,1),
     &              qfxl(1)*a(1,2)+qfxl(2)*a(2,2)+qfxl(3)*a(3,2),
     &              qfxl(1)*a(1,3)+qfxl(2)*a(2,3)+qfxl(3)*a(3,3),
     &              orname(iorien)(1:20)
            endif
         enddo
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
 -lint

+incdir+$OVM_HOME/src
$OVM_HOME/src/ovm_pkg.sv
my_env_pkg.sv
top.sv

ext: f


text:
module multiply_two
  implicit none
  integer, parameter, public :: two = 2
contains
  function multiply_two_numbers(a, b) result(c)
    real, intent(in) :: a, b
    real :: c
    c = a * b
    return
  end function multiply_two_numbers
end module multiply_two

max_stars_repo_name: DWesl/setup-cygwin-ci

max_issues_repo_name: DWesl/setup-cygwin-ci

max_forks_repo_name: DWesl/setup-cygwin-ci


text:
*------------------------------------------------------------------------*
      subroutine dia4op(me_dia,ecore,xdia,str_info,orb_info)
*------------------------------------------------------------------------*
*     set up diagonal hamiltonian
*     ecore contains the core energy
*     xdia contains the diagonal one-electron Hamiltonian (Fock matrix)
*     andreas, spring 2007
*------------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'
      include 'def_filinf.h'
      include 'def_operator.h'
      include 'def_me_list.h'
      include 'ifc_memman.h'
      include 'ifc_baserout.h'
      include 'hpvxseq.h'
      include 'multd2h.h'

      integer, parameter ::
     &     ntest = 00

      real(8), intent(in) ::
     &     ecore, xdia(*)
      type(me_list), intent(in) ::
     &     me_dia
      type(strinf), intent(in), target ::
     &     str_info
      type(orbinf), intent(in), target ::
     &     orb_info

      logical ::
     &     first, first_str, open_close_ffdia
      integer ::
     &     iblk, msamax, mscmax, ica, ihpvdx, ihpv, nsym, ngas,
     &     msa, msc, igama, igamc, igamstr, ms_str,
     &     idxms, igrph, ndis, idis, did, njoined,
     &     ncidx, naidx, maxidx, nidx, idx, ifree,
     &     nblk, nblkmax, idum, maxbuff,
     &     maxstrbuf, nstrbuf, nloop,
     &     istr, nouter, n_inner1, n_inner2, lenblk,
     &     iouter, ioffbuf, ioffbuf0, ioff_inner1, ioff_inner2,
     &     idxbuf, idx_inner1, len_blk

      integer ::
     &     msdst(ngastp,2), igamdst(ngastp,2),
     &     ioff_xsum(2*ngastp), nstr(2*ngastp), nincr(2*ngastp)
      real(8) ::
     &     xsum_outer, xsum_i1, fac,
     &     cpu, sys, wall, cpu0, sys0, wall0

      ! some pointers for easy typing and efficency
      integer, pointer ::
     &     iocc(:,:), idx_graph(:,:),
     &     igas_restr(:,:,:,:,:), mostnd(:,:,:),
     &     igamorb(:), ngas_hpv(:), idx_gas(:)

      ! allocatable stuff:
      integer, pointer ::
     &     idxorb(:),idxspn(:),idxdss(:)
      real(8), pointer ::
     &     buffer(:), xsum(:)

      type(filinf), pointer ::
     &     ffdia
      type(operator), pointer ::
     &     op

      real(8), external ::
     &     dnrm2
      logical, external ::
     &     next_msgamdist, next_string

      ifree = mem_setmark('dia4op')

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'this is dia4op')
      end if

      op => me_dia%op
      ffdia => me_dia%fhand

      njoined = op%njoined

      if (njoined.ne.1)
     &     call quit(1,'dia4op','njoined!=1 occurred!')

      call atim_csw(cpu0,sys0,wall0)

      igas_restr => str_info%igas_restr
      mostnd => orb_info%mostnd
      igamorb => orb_info%igamorb
      ngas_hpv => orb_info%ngas_hpv
      idx_gas => orb_info%idx_gas
      nsym = orb_info%nsym
      ngas = orb_info%ngas

      open_close_ffdia = ffdia%unit.le.0
      if (open_close_ffdia) call file_open(ffdia)

      maxbuff=0
      do iblk = 1, op%n_occ_cls
        iocc => op%ihpvca_occ(1:ngastp,1:2,iblk)
c        ! currently only tested for pure excitations, so:
c        if (imltlist(0,iocc(1:,1),ngastp,1).lt.3.or.
c     &      imltlist(0,iocc(1:,2),ngastp,1).lt.3.or.
c     &      iocc(ihole,1).gt.0 .or.
c     &      iocc(ipart,2).gt.0 ) then
c          call wrt_occ(lulog,op%ihpvca_occ(1,1,iblk))
c          call quit(1,'dia4op',
c     &         'routine not tested for this kind of occupation')
c        end if
        maxbuff = max(maxbuff,me_dia%len_op_occ(iblk))
      end do

      ! we should of course use a better memory management
      ! when going to large calculations
      if (maxbuff.gt.ifree) then
        write(lulog,*) 'maxbuff, ifree: ',maxbuff,ifree
        call mem_map(.true.)
        call quit(1,'dia4op','inadequate memory management')
      end if

      if (ntest.ge.100)
     &     write(lulog,*)'allocating result buffer of size: ',maxbuff

      ifree = mem_alloc_real(buffer,maxbuff,'buffer')
      ! loop over operator elements
      occ_cls: do iblk = 1, op%n_occ_cls

        if(op%formal_blk(iblk))cycle

        ! buffer: start from new
        idxbuf = 0
        buffer(1:maxbuff) = 0d0

        iocc => op%ihpvca_occ(1:ngastp,1:2,iblk)
        idx_graph => me_dia%idx_graph(1:ngastp,1:2,iblk)

        if (ntest.ge.100) then
          write(lulog,*) 'current occupation class:'
          call wrt_occ(lulog,op%ihpvca_occ(1,1,iblk))
        end if

        ncidx = op%ica_occ(1,iblk)
        naidx = op%ica_occ(2,iblk)
        maxidx = max(ncidx,naidx)

        mscmax = ncidx
        msamax = naidx

        ifree = mem_setmark('dia4op_str')
        ifree = mem_alloc_int(idxorb,maxidx,'idxorb')
        ifree = mem_alloc_int(idxspn,maxidx,'idxspn')
        ifree = mem_alloc_int(idxdss,maxidx,'idxdss')

        ioffbuf0 = 0
        ! inquire maximum number of strings:
        maxstrbuf = 0
        do msa = msamax, -msamax, -2
          msc = msa + me_dia%mst
          if (abs(msc).gt.mscmax) cycle
          do igama = 1, nsym
            igamc = multd2h(igama,me_dia%gamt)
            first = .true.
            distr_loop0: do
              if (.not.next_msgamdist(first,
     &           msa,msc,igama,igamc,iocc,nsym,
     &           msdst,igamdst)) exit distr_loop0
              first = .false.

              nstrbuf = 0
              do ica = 1, 2
                do ihpv = 1, ngastp
                  if (iocc(ihpv,ica).eq.0) cycle
                  igrph = idx_graph(ihpv,ica)
                  igamstr = igamdst(ihpv,ica)
                  idxms = (iocc(ihpv,ica)-msdst(ihpv,ica))/2+1
                  nstrbuf = nstrbuf +
     &                 str_info%g(igrph)%lenstr_gm(igamstr,idxms)
                end do
              end do
              maxstrbuf = max(maxstrbuf,nstrbuf)
            end do distr_loop0
          end do
        end do

        if (ntest.ge.100)
     &       write(lulog,*) 'allocating xsum-buffer of size ',maxstrbuf
        ! buffer for orbital energy sums
        ifree = mem_alloc_real(xsum,maxstrbuf,'xsum')

        idxms = 0
        msa_loop: do msa = msamax, -msamax, -2

          ! C <-> A  means alpha <-> beta !!
          msc = msa + me_dia%mst
          if (abs(msc).gt.mscmax) cycle msa_loop
          idxms = idxms+1

          igama_loop: do igama = 1, nsym
            igamc = multd2h(igama,me_dia%gamt)

            if (me_dia%len_op_gmo(iblk)%gam_ms(igama,idxms).eq.0)
     &           cycle igama_loop

c            idis = 0
c            first = .true.
            ! for pure P/H excitations we will not need this:
            ! separate loop / procedure ?
            ndis = me_dia%off_op_gmox(iblk)%ndis(igama,idxms)
            distr_loop: do idis = 1, ndis

              did = me_dia%off_op_gmox(iblk)%did(idis,igama,idxms)

              call did2msgm(msdst,igamdst,did,
     &               iocc,nsym,njoined)

c              if (.not.next_msgamdist(first,
c     &           msa,msc,igama,igamc,iocc,
c     &                   orb_info%nsym,
c     &           msdst,igamdst)) exit distr_loop
c              first = .false.
c              idis = idis+1

c              if (me_dia%len_op_gmox(iblk)%
c     &             d_gam_ms(idis,igama,idxms).eq.0)
c     &             cycle distr_loop

              if (ntest.ge.100) then
                write(lulog,*) 'msc,msa,igamc,igama: ',
     &               msc,msa,igamc,igama
                write(lulog,*) 'current distribution (MS,GAMMA):'
                call wrt_occ(lulog,msdst)
                call wrt_occ(lulog,igamdst)
              end if

              ! strings for HPV/CA
              nloop = 0
              istr = 0
              ! sequence: from outer to inner loop
              do ihpvdx = ngastp, 1, -1
                ihpv = hpvxseq(ihpvdx)
                do ica = 1, 2
                  if (iocc(ihpv,ica).eq.0) cycle
                  fac = 1d0
                  if (ica.eq.2) fac = -1d0
                  nloop = nloop+1
                  first_str = .true.
                  ioff_xsum(nloop) = istr
                  nidx = iocc(ihpv,ica)
                  igrph = idx_graph(ihpv,ica)
                  ms_str = msdst(ihpv,ica)
                  igamstr = igamdst(ihpv,ica)
c dbg
c                  print *,'ihpv,ica,iocc: ',ihpv,ica,nidx,igrph,ms_str,
c     &                 igamstr
c dbg
                  str_loop: do
                    if (.not.next_string(idxorb,idxspn,idxdss,
     &                 nidx,ms_str,igamstr,first_str,
     &                 igas_restr(1,1,1,1,igrph),
     &                 mostnd(1,1,idx_gas(ihpv)),igamorb,
     &                 nsym,ngas_hpv(ihpv))
     &                 ) exit str_loop

                    first_str = .false.

                    istr = istr+1
                    xsum(istr) = 0d0
                    do idx = 1, nidx
                      ! need to patch for UHF:
                      xsum(istr) = xsum(istr) + fac*xdia(idxorb(idx))
                    end do
                  end do str_loop
                  nstr(nloop) = istr-ioff_xsum(nloop)
                  ! do not consider empty strings
                  if (nstr(nloop).eq.0) nloop = nloop-1
                end do
              end do
c dbg
c              do idx = 1, nloop
c                print *,'idx, off, len: ',idx, ioff_xsum(idx), nstr(idx)
c                print *,'xsum = ',xsum(ioff_xsum(idx)+1:
c     &                                 ioff_xsum(idx)+nstr(idx))
c              end do
c dbg
              ! no contribution at all?
              if (nloop.eq.0) cycle

              ! only one loop?
              if (nloop.eq.1) then
                buffer(ioffbuf0+1:ioffbuf0+nstr(1)) =
     &               xsum(1:nstr(1))
                ioffbuf0 = ioffbuf0+nstr(1)
                cycle
              end if

              ! more than one loop:
              ! setup block-length information
              nouter = nloop-2
              n_inner2 = nstr(nloop)
              n_inner1 = nstr(nloop-1)
              lenblk = n_inner1*n_inner2
              do iouter = nouter, 1, -1
                nincr(iouter) = lenblk
                lenblk = lenblk*nstr(iouter)
              end do

              ! assemble the actual diagonal:
              ioffbuf = ioffbuf0
              ioff_inner1 = ioff_xsum(nouter+1)
              ioff_inner2 = ioff_xsum(nouter+2)
              do while(ioffbuf.lt.ioffbuf0+lenblk)
                ! collect contributions from formal outer loops
                ! (over more than two strings):
                idxbuf = ioffbuf - ioffbuf0 !+1
                xsum_outer = 0d0
                do iouter = 1, nouter
                  idx = ioff_xsum(iouter) + idxbuf/nincr(iouter)+1
                  xsum_outer = xsum_outer + xsum(idx)
c dbg
c                  if (idx.gt.nstr(iouter).or.idx.lt.1) then
c                    write(6,*) 'range(O) iouter,idx,nstr: ',
c     &                 iouter,idx,nstr(iouter)
c                    write(6,*) 'range(O) '//
c     &                   'idxbuf,nincr(iouter),ioff_xsum(iouter): ',
c     &                   idxbuf,nincr(iouter),
c     &                 ioff_xsum(iouter)
c                    write(6,*) 'off,idx/nicr: ',
c     &                   ioff_xsum(iouter),idxbuf/nincr(iouter),
c     &                   ioff_xsum(iouter) + idxbuf/nincr(iouter)+1
c                  end if
c dbg
                  idxbuf = mod(idxbuf,nincr(iouter))
                end do
                ! explicit loops for the two innermost strings
                do idx_inner1 = ioff_inner1+1, ioff_inner1+n_inner1
                  xsum_i1 = xsum(idx_inner1) + xsum_outer
c dbg
c                  if (ioffbuf+1.lt.1) write(6,*)'range1:',ioffbuf+1
c                  if (ioffbuf+n_inner2.gt.maxbuff)
c     &                 write(6,*)'range2:',ioffbuf+n_inner2
c dbg
                  buffer(ioffbuf+1:ioffbuf+n_inner2) =
     &                xsum(ioff_inner2+1:ioff_inner2+n_inner2) + xsum_i1
                  ioffbuf = ioffbuf+n_inner2
                end do
              end do
              ioffbuf0 = ioffbuf

            end do distr_loop

          end do igama_loop
        end do msa_loop

        ! add core energy
        len_blk = me_dia%len_op_occ(iblk)
        buffer(1:len_blk) = buffer(1:len_blk) + ecore
c dbg
c        if (ntest.ge.100)
c     &       print *,'final buffer: ',
c     &       buffer(1:me_dia%len_op_occ(iblk))
c dbg
        ! put buffer to disc
        call put_vec(ffdia,buffer,me_dia%off_op_occ(iblk)+1
     &       + ffdia%length_of_record*(ffdia%current_record-1),
     &                            me_dia%off_op_occ(iblk)+len_blk,
     &       + ffdia%length_of_record*(ffdia%current_record-1))

        ifree = mem_flushmark()

      end do occ_cls

      if (open_close_ffdia) call file_close_keep(ffdia)

      ifree = mem_flushmark()

      call atim_csw(cpu,sys,wall)

      if (iprlvl.ge.5)
     &     call prtim(lulog,'time in dia4op ',
     &                cpu-cpu0,sys-sys0,wall-wall0)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      subroutine splie2(x2a,ya,m,n,y2a)

      parameter (nn=100)
      integer m,n,j,k
      real x2a(n),ya(m,n),y2a(m,n),ytmp(nn),y2tmp(nn)
      do 13 j=1,m
        do 11 k=1,n
          ytmp(k)=ya(j,k)
11      continue
        call spline(x2a,ytmp,n,1.e30,1.e30,y2tmp)
        do 12 k=1,n
          y2a(j,k)=y2tmp(k)
12      continue
13    continue
      return
      end

max_stars_repo_name: xcorail/OTB

max_issues_repo_name: xcorail/OTB

max_forks_repo_name: xcorail/OTB


program test_ncdres_fclose_crash_1
    use nc_diag_res, only: nc_diag_load_resource_file, &
        nc_diag_close_resource_file

    ! Nothing is open, so this should fail!
    call nc_diag_close_resource_file
end program test_ncdres_fclose_crash_1


text:
Alisa Palma is a candidate for ASUCD ASUCD Senate Senate in the Fall 2013 ASUCD Election running on the NOW slate.

 Candidate Statement

Hey Aggies!

My name is Alisa Palma. I am running for ASUCD Senate with the NOW Slate. I am a 2nd year Economics major with a minor in management. I am from San Ramon, California, and currently a member of the Economics and Business Student Association, as well as the Swim Club at UC Davis. I am also on the Club Finance Council board and the CSI Student Advisory Council where I help allocate a portion of ASUCD’s funds to student organizations.

I want to serve as your next ASUCD Senator so students will better understand where their fees paid to ASUCD are being spent. I will ensure that our student government is responsible with all of its financial expenses.

ADD ACCREDITED PROGRAMS TO EXPERIMENTAL COLLEGE
I wish to add accredited programs to the ASUCD Experimental College, such as nursing and kinesiology. UC Davis offers many elaborate graduate programs, and yet there are no undergraduate courses for students to gain an interest in pursuing a higher degree in those fields. More importantly, Sac State offers these courses, and I believe implementing them into the Experimental College will help make up for the difference between schools, especially when credits are nontransferable. It is vital that our campus is actively involved with helping students get more engaged with their future.

EXPANDING ON AND OFF CAMPUS COMPOSTING AVAILABILITY
I would like to expand composting availability within the Davis community by placing compost bins around the quad, off campus housing such as West Village, and participating Greek houses. Our school’s compost goes straight back into the campus for the agriculture department, as well as other uses. Many people waste food without even realizing it, and by providing more compost bins, students can help the campus go green by giving back to the environment and the campus. Smog has been a prevalent issue arising in California, and compost is a lowcost solution to reducing landfill emissions. Inserting more compost bins will make the campus more sustainable and therefore students would be utilizing their leftovers in helping reduce harmful pollutants.

I am excited to represent you by making changes to this campus in order for every student to gain a better outlook of their future. Please VOTE ALISA PALMA #1 FOR ASUCD SENATE AND NOW#26!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module m
implicit none
type :: vec
   real, allocatable :: x(:)
end type vec
contains
function ran(n) result(xran)
integer, intent(in) :: n
real :: xran(n)
call random_number(xran)
end function ran
!
elemental function moment(v,m) result(xmom)
! return the nth moment of the array component of v
type(vec), intent(in) :: v
integer  , intent(in) :: m
real                  :: xmom
real                  :: xmean
integer               :: n
n = max(1,size(v%x))
xmean = sum(v%x)/n
xmom  = sum((v%x-xmean)**m)/n
end function moment
end module m
!
program main
use m, only: moment,vec,ran
implicit none
integer            :: i, nmom=4
type(vec)          :: v
character (len=20) :: fmt = "(a10,*(f8.4))"
v = vec(x = ran(10**4))
! central moments of uniform distribution
print "(a10,*(i8))","moment",(i,i=1,nmom)
print fmt,"sim",moment(v,[(i,i=1,nmom)])
print fmt,"true",0.0,1/12.0,0.0,1/80.0
end program main
! sample output:
!     moment       1       2       3       4
!        sim  0.0000  0.0836  0.0000  0.0125
!       true  0.0000  0.0833  0.0000  0.0125

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
!! https://github.com/flang-compiler/flang/issues/313
!!
!! $ flang -c flang-20171123a.f90
!! F90-S-0142-size is not a component of this OBJECT (flang-20171123a.f90: 17)
!! F90-S-0075-Subscript, substring, or argument illegal in this context for size (flang-20171123a.f90: 17)
!! F90-S-0155-FUNCTION may not be declared character*(*) when in an INTERFACE - foo_func (flang-20171123a.f90: 14)
!!   0 inform,   0 warnings,   3 severes, 0 fatal for example
!!
!! The first error message concerning the declaration of STRING is erroneous.

module example

  type, abstract :: foo
  contains
    procedure(foo_size), deferred :: size
    procedure(foo_func), deferred :: func
  end type

  abstract interface
    pure integer function foo_size (this)
      import foo
      class(foo), intent(in) :: this
    end function
    function foo_func (this) result (string)
      import :: foo
      class(foo) :: this
      character(this%size()) :: string
    end function
  end interface

end module

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
c     DRDLESUM
c>> 1995-05-28 DRDLESUM Krogh  Changes to use M77CON
c>> 1994-08-09 DRDLESUM WVS  Set up for CHGTYP
c>> 1994-07-14 DRDLESUM CLL
c>> 1992-04-29 DRDLESUM CAO Replaced '1' in format.
c>> 1991-11-19 DRDLESUM CLL
c>> 1987-12-09 DRDLESUM Lawson  Initial Code.
c--D replaces "?": ?LESUM, DR?LESUM
c
c     Demonstration driver for evaluation of a Legendre series.
c     ------------------------------------------------------------------
      integer j
      double precision x,a(0:5),y,w,z
      data a/0.07d0, 0.27d0, 0.20d0, 0.28d0, 0.08d0, 0.08d0/
c     ------------------------------------------------------------------
      print '(1x,3x,a1,14x,a1,17x,a1/)','x','y','z'
      do 20 j = -10,10,2
        x = dble(j) /10.d0
        call dlesum (x,5,a,y)
        w = 0.35d0 * (x**4) + 0.63d0 * (x**5)
        z = y - w
        print '(1x,f5.2,5x,g15.7,g15.2)',x,y,z
   20 continue
      end


text:
!***************************************************************************************************************************
! This is the code from the section (DYNAMIC == INTEGR) lines 6419 - 6541 of the original CSCAS code. The names of the
! dummy arguments are the same as in the original CSCAS code and the call statement and are declared here. The variables
! that are not arguments are declared in module YCA_First_Trans_m. Unless identified as by MF, all comments are those of
! the original CSCAS.FOR code.
!
! Subroutine YCA_Integ_WthrSum calculates summaries of weather and soil variables and par utilization.
!***************************************************************************************************************************

    SUBROUTINE YCA_Integ_WthrSum ( &
        CO2         , DOY         , DRAIN       , IRRAMT      , RAIN        , RUNOFF      , SRAD        , TLCHD       , &
        TMAX        , TMIN        , TNIMBSOM    , TNOXD       , TOMINFOM    , TOMINSOM    , TOMINSOM1   , TOMINSOM2   , &
        TOMINSOM3   , YEAR        , ISWWAT      &
        )

        USE YCA_First_Trans_m
        USE YCA_Control_Plant

        IMPLICIT NONE

        INTEGER DOY         , YEAR

        REAL    CO2         , DRAIN       , IRRAMT      , RAIN        , RUNOFF      , SRAD        , TLCHD       , TMAX
        REAL    TMIN        , TNIMBSOM    , TNOXD       , TOMINFOM    , TOMINSOM    , TOMINSOM1   , TOMINSOM2   , TOMINSOM3

        CHARACTER(LEN=1) ISWWAT
        !-----------------------------------------------------------------------
        !         Calculate weather and soil summary variables
        !-----------------------------------------------------------------------

        ! Cumulatives
        !LPM 24MAR2016 TTCUM estimated after germination, before is DAGERM and GEUCUM
        IF (GERMFR > 0.0) THEN
            IF (DAE > 0.0) THEN
                TTCUM = TTCUM + TT
            ELSE
                TTCUM = TTCUM + TTGEM
            ENDIF
        ELSE
            TTCUM = 0.0
        ENDIF

        IF (ISWWAT == 'Y') THEN
            DAWWP = DAWWP + (TT*WFG) !LPM 31JUL2015 Added to have a new clock with water stress
        ELSE
            DAWWP = TTCUM
        ENDIF
        IF (DAWWP > 900.0) THEN
            TTCUMLS = TTCUMLS + TTlfsize   ! LPM 12JUL2015 added to consider a different optimum temperature for potential leaf size
            IF (ISWWAT == 'Y') THEN
                DALS = DALS + (TTlfsize*WFG) !LPM 24APR2015 Added to have a new clock with water stress
            ELSE
                DALS = TTCUMLS
            ENDIF
        ENDIF
        RAINC = RAINC + RAIN
        DRAINC = DRAINC + DRAIN
        RUNOFFC = RUNOFFC + RUNOFF
        IRRAMTC = IRRAMTC + IRRAMT
        SRADC = SRADC + SRAD
        PARMJC = PARMJC + PARMJFAC*SRAD
        PARMJIC = PARMJIC + PARMJFAC*SRAD*PARI + PARMJIADJ
        TOMINC = TOMINC + TOMIN
        TOFIXC = TOFIXC + TNIMBSOM
        TOMINFOMC = TOMINFOMC + TOMINFOM
        TOMINSOMC = TOMINSOMC + TOMINSOM
        IF (TOMINSOM1 >= 0.0) THEN
            TOMINSOM1C = TOMINSOM1C + TOMINSOM1
            TOMINSOM2C = TOMINSOM2C + TOMINSOM2
            TOMINSOM3C = TOMINSOM3C + TOMINSOM3
        ELSE
            TOMINSOM1C = -99.0
            TOMINSOM2C = -99.0
            TOMINSOM3C = -99.0
        ENDIF
        TLCHC = TLCHC + TLCHD
        TNOXC = TNOXC + TNOXD
        ! Extremes
        TMAXX = AMAX1(TMAXX,TMAX)
        TMINN = AMIN1(TMINN,TMIN)
        CO2MAX = AMAX1(CO2MAX,CO2)

        ! Growing season means

        TMEANNUM = TMEANNUM + 1
        TMEANSUM = TMEANSUM + TMEAN

        ! 20-day means
        SRAD20S = 0.0
        TMEAN20S = 0.0
        STRESS20S = 0.0
        STRESS20NS = 0.0
        STRESS20WS = 0.0
        TT20S = 0.0
        DO L = 20,2,-1
            SRADD(L) = SRADD(L-1)
            SRAD20S = SRAD20S + SRADD(L)
            TMEAND(L) = TMEAND(L-1)
            TMEAN20S = TMEAN20S + TMEAND(L)
            STRESS(L) = STRESS(L-1)
            STRESSN(L) = STRESSN(L-1)
            STRESSW(L) = STRESSW(L-1)
            STRESS20S = STRESS20S + STRESS(L)
            STRESS20NS = STRESS20NS + STRESSN(L)
            STRESS20WS = STRESS20WS + STRESSW(L)
            TTD(L) = TTD(L-1)
            TT20S = TT20S + TTD(L)
            WUPRD(L) = WUPRD(L-1)
        ENDDO
        SRADD(1) = SRAD
        SRAD20S = SRAD20S + SRAD
        TMEAND(1) = TMEAN
        TMEAN20S = TMEAN20S + TMEAND(1)
        STRESS(1) = AMIN1(WFG,NFG)
        STRESSN(1) = NFG
        STRESSW(1) = WFG
        STRESS20S = STRESS20S + STRESS(1)
        STRESS20NS = STRESS20NS + STRESSN(1)
        STRESS20WS = STRESS20WS + STRESSW(1)
        TTD(1) = TT
        TT20S = TT20S + TTD(1)
        WUPRD(1) = AMAX1(0.0,AMIN1(10.0,WUPR))
        IF (TMEANNUM >= 20) THEN
            IF (TMEANNUM <= 20) TMEAN20P = TMEAN20S/20.0
            SRAD20 = SRAD20S/20.0
            TMEAN20 = TMEAN20S/20.0
            TT20 = TT20S/20.0
            STRESS20 = STRESS20S/20.0
            STRESS20N = STRESS20NS/20.0
            STRESS20W = STRESS20WS/20.0
        ELSE
            SRAD20 = 0.0
            IF (TMEANNUM == 1) THEN       !CHP
              TT20 = TT20S                !CHP
            ELSE                          !CHP
              TT20 = TT20S/(TMEANNUM-1)   !LPM 28FEB15 to have a value different than 0
            ENDIF                         !CHP
            TMEAN20 = 0.0
            STRESS20 = 0.0
            STRESS20N = 0.0
            STRESS20N = 0.0
        ENDIF

        ! Monthly means
        CALL Calendar (year,doy,dom,month)
        IF (DOM > 1) THEN
            TMAXSUM = TMAXSUM + TMAX
            TMINSUM = TMINSUM + TMIN
            DAYSUM = DAYSUM + 1.0
        ELSE
            IF (DAYSUM > 0) THEN
                IF (TMAXM < TMAXSUM/DAYSUM) TMAXM=TMAXSUM/DAYSUM
                IF (TMINM > TMINSUM/DAYSUM) TMINM=TMINSUM/DAYSUM
            ENDIF
            TMAXSUM = TMAX
            TMINSUM = TMIN
            DAYSUM =  1
        ENDIF

        !-----------------------------------------------------------------------
        !         Calculate PAR utilization efficiencies
        !-----------------------------------------------------------------------

        IF (PARMJC > 0.0) THEN
            PARUEC = AMAX1(0.0,(totalWeight()+SENTOPLITTER+SENROOT-SEEDUSE)* PLTPOP / PARMJC)
        ENDIF
        IF (PARMJIC > 0.0) THEN
            PARIUED = AMAX1(0.0,(totalWeight()+SENTOPLITTER+SENROOT-SEEDUSE)* PLTPOP / PARMJIC)
        ENDIF
        IF (CARBOBEG > 0.0) THEN
            PARIUE = (CARBOBEG*PLTPOP)/(PARMJFAC*SRAD*PARI)
        ENDIF

    END SUBROUTINE YCA_Integ_WthrSum

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
asset: %font
    %font svar font.fnt
    %font svar font.size
    %font svar font.flags
: >fnt ( font - ALLEGRO_FONT ) font.fnt @ ;

: reload-font  ( font - )
    >r  r@ srcfile count findfile zstring  r@ font.size @ 1i  r@ font.flags @  al_load_font  r> font.fnt ! ;

: unload-font  ( font - )
    font.fnt @ al_destroy_font ;

: init-font  ( path c size flags font - )
    >r  r@ font.flags !  r@ font.size !
    r@ srcfile place  ['] reload-font ['] unload-font r@ register
    r> reload-font ;

: font:  ( path c size flags - <name> )
    create  %font *struct  init-font ;

create (chr)  0 c, 0 c,
: chrw    ( font chr - n ) (chr) c!  >fnt (chr) al_get_text_width 1p ;
: chrh    ( font - n ) >fnt al_get_font_line_height 1p ;

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine


text:
! stand alone driver to read in lattice vectors and atomic basis
! and return the spacegroup rotations and shift.
PROGRAM spacegroup
use num_types
use symmetry
implicit none

real(dp), dimension(3,3) :: LV
real(dp), allocatable    :: d(:,:)
integer, allocatable     :: typ(:)
character(120)           :: LVfile, dFile
integer iD, nD
integer nArg

nArg = command_argument_count()
if (nArg)/= 2
stop "ERROR: requires two argments. File names for lattice vectors and for atomic basis vectors"

! Read the two file names for LV and d from the command line

! Cound the number of atoms, allocate atom positions and typ

! Call the spacegroup routine

! Write out the

END PROGRAM spacegroup

max_stars_repo_name: braydenbekker/enumlib

max_issues_repo_name: braydenbekker/enumlib

max_forks_repo_name: braydenbekker/enumlib


text:
      subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest,eps,
     *  nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
      implicit none
c  subroutine sphere determines a smooth bicubic spherical spline
c  approximation s(teta,phi), 0 <= teta <= pi ; 0 <= phi <= 2*pi
c  to a given set of data points (teta(i),phi(i),r(i)),i=1,2,...,m.
c  such a spline has the following specific properties
c
c    (1) s(0,phi)  = constant   0 <=phi<= 2*pi.
c
c    (2) s(pi,phi) = constant   0 <=phi<= 2*pi
c
c         j             j
c        d s(teta,0)   d s(teta,2*pi)
c    (3) ----------- = ------------   0 <=teta<=pi, j=0,1,2
c             j             j
c        d phi         d phi
c
c        d s(0,phi)    d s(0,0)             d s(0,pi/2)
c    (4) ----------  = -------- *cos(phi) + ----------- *sin(phi)
c        d teta        d teta               d teta
c
c        d s(pi,phi)   d s(pi,0)            d s(pi,pi/2)
c    (5) ----------- = ---------*cos(phi) + ------------*sin(phi)
c        d teta        d teta               d teta
c
c  if iopt =-1 sphere calculates a weighted least-squares spherical
c  spline according to a given set of knots in teta- and phi- direction.
c  if iopt >=0, the number of knots in each direction and their position
c  tt(j),j=1,2,...,nt ; tp(j),j=1,2,...,np are chosen automatically by
c  the routine. the smoothness of s(teta,phi) is then achieved by mini-
c  malizing the discontinuity jumps of the derivatives of the spline
c  at the knots. the amount of smoothness of s(teta,phi) is determined
c  by the condition that fp = sum((w(i)*(r(i)-s(teta(i),phi(i))))**2)
c  be <= s, with s a given non-negative constant.
c  the spherical spline is given in the standard b-spline representation
c  of bicubic splines and can be evaluated by means of subroutine bispev
c
c calling sequence:
c     call sphere(iopt,m,teta,phi,r,w,s,ntest,npest,eps,
c    *  nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
c
c parameters:
c  iopt  : integer flag. on entry iopt must specify whether a weighted
c          least-squares spherical spline (iopt=-1) or a smoothing
c          spherical spline (iopt=0 or 1) must be determined.
c          if iopt=0 the routine will start with an initial set of knots
c          tt(i)=0,tt(i+4)=pi,i=1,...,4;tp(i)=0,tp(i+4)=2*pi,i=1,...,4.
c          if iopt=1 the routine will continue with the set of knots
c          found at the last call of the routine.
c          attention: a call with iopt=1 must always be immediately pre-
c                     ceded by another call with iopt=1 or iopt=0.
c          unchanged on exit.
c  m     : integer. on entry m must specify the number of data points.
c          m >= 2. unchanged on exit.
c  teta  : real array of dimension at least (m).
c  phi   : real array of dimension at least (m).
c  r     : real array of dimension at least (m).
c          before entry,teta(i),phi(i),r(i) must be set to the spherical
c          co-ordinates of the i-th data point, for i=1,...,m.the order
c          of the data points is immaterial. unchanged on exit.
c  w     : real array of dimension at least (m). before entry, w(i) must
c          be set to the i-th value in the set of weights. the w(i) must
c          be strictly positive. unchanged on exit.
c  s     : real. on entry (in case iopt>=0) s must specify the smoothing
c          factor. s >=0. unchanged on exit.
c          for advice on the choice of s see further comments
c  ntest : integer. unchanged on exit.
c  npest : integer. unchanged on exit.
c          on entry, ntest and npest must specify an upper bound for the
c          number of knots required in the teta- and phi-directions.
c          these numbers will also determine the storage space needed by
c          the routine. ntest >= 8, npest >= 8.
c          in most practical situation ntest = npest = 8+sqrt(m/2) will
c          be sufficient. see also further comments.
c  eps   : real.
c          on entry, eps must specify a threshold for determining the
c          effective rank of an over-determined linear system of equat-
c          ions. 0 < eps < 1.  if the number of decimal digits in the
c          computer representation of a real number is q, then 10**(-q)
c          is a suitable value for eps in most practical applications.
c          unchanged on exit.
c  nt    : integer.
c          unless ier=10 (in case iopt >=0), nt will contain the total
c          number of knots with respect to the teta-variable, of the
c          spline approximation returned. if the computation mode iopt=1
c          is used, the value of nt should be left unchanged between
c          subsequent calls.
c          in case iopt=-1, the value of nt should be specified on entry
c  tt    : real array of dimension at least ntest.
c          on succesful exit, this array will contain the knots of the
c          spline with respect to the teta-variable, i.e. the position
c          of the interior knots tt(5),...,tt(nt-4) as well as the
c          position of the additional knots tt(1)=...=tt(4)=0 and
c          tt(nt-3)=...=tt(nt)=pi needed for the b-spline representation
c          if the computation mode iopt=1 is used, the values of tt(1),
c          ...,tt(nt) should be left unchanged between subsequent calls.
c          if the computation mode iopt=-1 is used, the values tt(5),
c          ...tt(nt-4) must be supplied by the user, before entry.
c          see also the restrictions (ier=10).
c  np    : integer.
c          unless ier=10 (in case iopt >=0), np will contain the total
c          number of knots with respect to the phi-variable, of the
c          spline approximation returned. if the computation mode iopt=1
c          is used, the value of np should be left unchanged between
c          subsequent calls.
c          in case iopt=-1, the value of np (>=9) should be specified
c          on entry.
c  tp    : real array of dimension at least npest.
c          on succesful exit, this array will contain the knots of the
c          spline with respect to the phi-variable, i.e. the position of
c          the interior knots tp(5),...,tp(np-4) as well as the position
c          of the additional knots tp(1),...,tp(4) and tp(np-3),...,
c          tp(np) needed for the b-spline representation.
c          if the computation mode iopt=1 is used, the values of tp(1),
c          ...,tp(np) should be left unchanged between subsequent calls.
c          if the computation mode iopt=-1 is used, the values tp(5),
c          ...tp(np-4) must be supplied by the user, before entry.
c          see also the restrictions (ier=10).
c  c     : real array of dimension at least (ntest-4)*(npest-4).
c          on succesful exit, c contains the coefficients of the spline
c          approximation s(teta,phi).
c  fp    : real. unless ier=10, fp contains the weighted sum of
c          squared residuals of the spline approximation returned.
c  wrk1  : real array of dimension (lwrk1). used as workspace.
c          if the computation mode iopt=1 is used the value of wrk1(1)
c          should be left unchanged between subsequent calls.
c          on exit wrk1(2),wrk1(3),...,wrk1(1+ncof) will contain the
c          values d(i)/max(d(i)),i=1,...,ncof=6+(np-7)*(nt-8)
c          with d(i) the i-th diagonal element of the reduced triangular
c          matrix for calculating the b-spline coefficients. it includes
c          those elements whose square is less than eps,which are treat-
c          ed as 0 in the case of presumed rank deficiency (ier<-2).
c  lwrk1 : integer. on entry lwrk1 must specify the actual dimension of
c          the array wrk1 as declared in the calling (sub)program.
c          lwrk1 must not be too small. let
c            u = ntest-7, v = npest-7, then
c          lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m
c  wrk2  : real array of dimension (lwrk2). used as workspace, but
c          only in the case a rank deficient system is encountered.
c  lwrk2 : integer. on entry lwrk2 must specify the actual dimension of
c          the array wrk2 as declared in the calling (sub)program.
c          lwrk2 > 0 . a save upper bound  for lwrk2 = 48+21*v+7*u*v+
c          4*(u-1)*v**2 where u,v are as above. if there are enough data
c          points, scattered uniformly over the approximation domain
c          and if the smoothing factor s is not too small, there is a
c          good chance that this extra workspace is not needed. a lot
c          of memory might therefore be saved by setting lwrk2=1.
c          (see also ier > 10)
c  iwrk  : integer array of dimension (kwrk). used as workspace.
c  kwrk  : integer. on entry kwrk must specify the actual dimension of
c          the array iwrk as declared in the calling (sub)program.
c          kwrk >= m+(ntest-7)*(npest-7).
c  ier   : integer. unless the routine detects an error, ier contains a
c          non-positive value on exit, i.e.
c   ier=0  : normal return. the spline returned has a residual sum of
c            squares fp such that abs(fp-s)/s <= tol with tol a relat-
c            ive tolerance set to 0.001 by the program.
c   ier=-1 : normal return. the spline returned is a spherical
c            interpolating spline (fp=0).
c   ier=-2 : normal return. the spline returned is the weighted least-
c            squares constrained polynomial . in this extreme case
c            fp gives the upper bound for the smoothing factor s.
c   ier<-2 : warning. the coefficients of the spline returned have been
c            computed as the minimal norm least-squares solution of a
c            (numerically) rank deficient system. (-ier) gives the rank.
c            especially if the rank deficiency which can be computed as
c            6+(nt-8)*(np-7)+ier, is large the results may be inaccurate
c            they could also seriously depend on the value of eps.
c   ier=1  : error. the required storage space exceeds the available
c            storage space, as specified by the parameters ntest and
c            npest.
c            probably causes : ntest or npest too small. if these param-
c            eters are already large, it may also indicate that s is
c            too small
c            the approximation returned is the weighted least-squares
c            spherical spline according to the current set of knots.
c            the parameter fp gives the corresponding weighted sum of
c            squared residuals (fp>s).
c   ier=2  : error. a theoretically impossible result was found during
c            the iteration proces for finding a smoothing spline with
c            fp = s. probably causes : s too small or badly chosen eps.
c            there is an approximation returned but the corresponding
c            weighted sum of squared residuals does not satisfy the
c            condition abs(fp-s)/s < tol.
c   ier=3  : error. the maximal number of iterations maxit (set to 20
c            by the program) allowed for finding a smoothing spline
c            with fp=s has been reached. probably causes : s too small
c            there is an approximation returned but the corresponding
c            weighted sum of squared residuals does not satisfy the
c            condition abs(fp-s)/s < tol.
c   ier=4  : error. no more knots can be added because the dimension
c            of the spherical spline 6+(nt-8)*(np-7) already exceeds
c            the number of data points m.
c            probably causes : either s or m too small.
c            the approximation returned is the weighted least-squares
c            spherical spline according to the current set of knots.
c            the parameter fp gives the corresponding weighted sum of
c            squared residuals (fp>s).
c   ier=5  : error. no more knots can be added because the additional
c            knot would (quasi) coincide with an old one.
c            probably causes : s too small or too large a weight to an
c            inaccurate data point.
c            the approximation returned is the weighted least-squares
c            spherical spline according to the current set of knots.
c            the parameter fp gives the corresponding weighted sum of
c            squared residuals (fp>s).
c   ier=10 : error. on entry, the input data are controlled on validity
c            the following restrictions must be satisfied.
c            -1<=iopt<=1,  m>=2, ntest>=8 ,npest >=8, 0<eps<1,
c            0<=teta(i)<=pi, 0<=phi(i)<=2*pi, w(i)>0, i=1,...,m
c            lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m
c            kwrk >= m+(ntest-7)*(npest-7)
c            if iopt=-1: 8<=nt<=ntest , 9<=np<=npest
c                        0<tt(5)<tt(6)<...<tt(nt-4)<pi
c                        0<tp(5)<tp(6)<...<tp(np-4)<2*pi
c            if iopt>=0: s>=0
c            if one of these conditions is found to be violated,control
c            is immediately repassed to the calling program. in that
c            case there is no approximation returned.
c   ier>10 : error. lwrk2 is too small, i.e. there is not enough work-
c            space for computing the minimal least-squares solution of
c            a rank deficient system of linear equations. ier gives the
c            requested value for lwrk2. there is no approximation re-
c            turned but, having saved the information contained in nt,
c            np,tt,tp,wrk1, and having adjusted the value of lwrk2 and
c            the dimension of the array wrk2 accordingly, the user can
c            continue at the point the program was left, by calling
c            sphere with iopt=1.
c
c further comments:
c  by means of the parameter s, the user can control the tradeoff
c   between closeness of fit and smoothness of fit of the approximation.
c   if s is too large, the spline will be too smooth and signal will be
c   lost ; if s is too small the spline will pick up too much noise. in
c   the extreme cases the program will return an interpolating spline if
c   s=0 and the constrained weighted least-squares polynomial if s is
c   very large. between these extremes, a properly chosen s will result
c   in a good compromise between closeness of fit and smoothness of fit.
c   to decide whether an approximation, corresponding to a certain s is
c   satisfactory the user is highly recommended to inspect the fits
c   graphically.
c   recommended values for s depend on the weights w(i). if these are
c   taken as 1/d(i) with d(i) an estimate of the standard deviation of
c   r(i), a good s-value should be found in the range (m-sqrt(2*m),m+
c   sqrt(2*m)). if nothing is known about the statistical error in r(i)
c   each w(i) can be set equal to one and s determined by trial and
c   error, taking account of the comments above. the best is then to
c   start with a very large value of s ( to determine the least-squares
c   polynomial and the corresponding upper bound fp0 for s) and then to
c   progressively decrease the value of s ( say by a factor 10 in the
c   beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
c   approximation shows more detail) to obtain closer fits.
c   to choose s very small is strongly discouraged. this considerably
c   increases computation time and memory requirements. it may also
c   cause rank-deficiency (ier<-2) and endager numerical stability.
c   to economize the search for a good s-value the program provides with
c   different modes of computation. at the first call of the routine, or
c   whenever he wants to restart with the initial set of knots the user
c   must set iopt=0.
c   if iopt=1 the program will continue with the set of knots found at
c   the last call of the routine. this will save a lot of computation
c   time if sphere is called repeatedly for different values of s.
c   the number of knots of the spline returned and their location will
c   depend on the value of s and on the complexity of the shape of the
c   function underlying the data. if the computation mode iopt=1
c   is used, the knots returned may also depend on the s-values at
c   previous calls (if these were smaller). therefore, if after a number
c   of trials with different s-values and iopt=1, the user can finally
c   accept a fit as satisfactory, it may be worthwhile for him to call
c   sphere once more with the selected value for s but now with iopt=0.
c   indeed, sphere may then return an approximation of the same quality
c   of fit but with fewer knots and therefore better if data reduction
c   is also an important objective for the user.
c   the number of knots may also depend on the upper bounds ntest and
c   npest. indeed, if at a certain stage in sphere the number of knots
c   in one direction (say nt) has reached the value of its upper bound
c   (ntest), then from that moment on all subsequent knots are added
c   in the other (phi) direction. this may indicate that the value of
c   ntest is too small. on the other hand, it gives the user the option
c   of limiting the number of knots the routine locates in any direction
c   for example, by setting ntest=8 (the lowest allowable value for
c   ntest), the user can indicate that he wants an approximation which
c   is a cubic polynomial in the variable teta.
c
c  other subroutines required:
c    fpback,fpbspl,fpsphe,fpdisc,fpgivs,fprank,fprati,fprota,fporde,
c    fprpsp
c
c  references:
c   dierckx p. : algorithms for smoothing data on the sphere with tensor
c                product splines, computing 32 (1984) 319-342.
c   dierckx p. : algorithms for smoothing data on the sphere with tensor
c                product splines, report tw62, dept. computer science,
c                k.u.leuven, 1983.
c   dierckx p. : curve and surface fitting with splines, monographs on
c                numerical analysis, oxford university press, 1993.
c
c  author:
c    p.dierckx
c    dept. computer science, k.u. leuven
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c  creation date : july 1983
c  latest update : march 1989
c
c  ..
c  ..scalar arguments..
      real*8 s,eps,fp
      integer iopt,m,ntest,npest,nt,np,lwrk1,lwrk2,kwrk,ier
c  ..array arguments..
      real*8 teta(m),phi(m),r(m),w(m),tt(ntest),tp(npest),
     * c((ntest-4)*(npest-4)),wrk1(lwrk1),wrk2(lwrk2)
      integer iwrk(kwrk)
c  ..local scalars..
      real*8 tol,pi,pi2,one
      integer i,ib1,ib3,ki,kn,kwest,la,lbt,lcc,lcs,lro,j,
     * lbp,lco,lf,lff,lfp,lh,lq,lst,lsp,lwest,maxit,ncest,ncc,ntt,
     * npp,nreg,nrint,ncof,nt4,np4
c  ..function references..
      real*8 atan
c  ..subroutine references..
c    fpsphe
c  ..
c  set constants
      one = 0.1e+01
c  we set up the parameters tol and maxit.
      maxit = 20
      tol = 0.1e-02
c  before starting computations a data check is made. if the input data
c  are invalid,control is immediately repassed to the calling program.
      ier = 10
      if(eps.le.0. .or. eps.ge.1.) go to 80
      if(iopt.lt.(-1) .or. iopt.gt.1) go to 80
      if(m.lt.2) go to 80
      if(ntest.lt.8 .or. npest.lt.8) go to 80
      nt4 = ntest-4
      np4 = npest-4
      ncest = nt4*np4
      ntt = ntest-7
      npp = npest-7
      ncc = 6+npp*(ntt-1)
      nrint = ntt+npp
      nreg = ntt*npp
      ncof = 6+3*npp
      ib1 = 4*npp
      ib3 = ib1+3
      if(ncof.gt.ib1) ib1 = ncof
      if(ncof.gt.ib3) ib3 = ncof
      lwest = 185+52*npp+10*ntt+14*ntt*npp+8*(m+(ntt-1)*npp**2)
      kwest = m+nreg
      if(lwrk1.lt.lwest .or. kwrk.lt.kwest) go to 80
      if(iopt.gt.0) go to 60
      pi = atan(one)*4
      pi2 = pi+pi
      do 20 i=1,m
        if(w(i).le.0.) go to 80
        if(teta(i).lt.0. .or. teta(i).gt.pi) go to 80
        if(phi(i) .lt.0. .or. phi(i).gt.pi2) go to 80
  20  continue
      if(iopt.eq.0) go to 60
      ntt = nt-8
      if(ntt.lt.0 .or. nt.gt.ntest) go to 80
      if(ntt.eq.0) go to 40
      tt(4) = 0.
      do 30 i=1,ntt
         j = i+4
         if(tt(j).le.tt(j-1) .or. tt(j).ge.pi) go to 80
  30  continue
  40  npp = np-8
      if(npp.lt.1 .or. np.gt.npest) go to 80
      tp(4) = 0.
      do 50 i=1,npp
         j = i+4
         if(tp(j).le.tp(j-1) .or. tp(j).ge.pi2) go to 80
  50  continue
      go to 70
  60  if(s.lt.0.) go to 80
  70  ier = 0
c  we partition the working space and determine the spline approximation
      kn = 1
      ki = kn+m
      lq = 2
      la = lq+ncc*ib3
      lf = la+ncc*ib1
      lff = lf+ncc
      lfp = lff+ncest
      lco = lfp+nrint
      lh = lco+nrint
      lbt = lh+ib3
      lbp = lbt+5*ntest
      lro = lbp+5*npest
      lcc = lro+npest
      lcs = lcc+npest
      lst = lcs+npest
      lsp = lst+m*4
      call fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eps,tol,maxit,
     * ib1,ib3,ncest,ncc,nrint,nreg,nt,tt,np,tp,c,fp,wrk1(1),wrk1(lfp),
     * wrk1(lco),wrk1(lf),wrk1(lff),wrk1(lro),wrk1(lcc),wrk1(lcs),
     * wrk1(la),wrk1(lq),wrk1(lbt),wrk1(lbp),wrk1(lst),wrk1(lsp),
     * wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier)
  80  return
      end

max_stars_repo_name: FRidh/scipy

max_issues_repo_name: FRidh/scipy

max_forks_repo_name: FRidh/scipy


!This script is used to coupling the NWM streams. It returns the featureID
!of the intersected streams and the elementID of the intersected
!land boundary of the SCHISM grid.
!It assumes the land boundary is the 1st land boundary and there is only
!one land boundary.
!Inputs files:
!      NWM_shp_ll.nc : which is converted from the shape file using the NWM streams.
!      hgrid.ll: the gr3 format SCHISM grid; please only keep land
!      bnd segments needed for coupling (islands are OK). The # of land
!      bnd nodes is only used for dimensioning so give it a large number
!      if you want.
!             Note that the input files are in the same projection.
!             Here we use lat/lon projection
!      Epsilon: a small distance used to nudge the intersection point towards the downstream vetices
!               suggested value is 1.e-3, which is about 100m when using
!               lat/lon projection.
!      dd,mm,yyyy: are day, month, and year, respectively provided by user.
!
!Outputs files:
!      source_sink.in   : contains the element ID for each
!                         intersection of the NWM streams and the land boundary.
!      msource.th       : contains the salinity and temprature
!                         of the source element along the land boundary.
!                         Salinity is set to be 0, temp = -9999.
!      vsource.th       : input of the stream flows of source elements.
!      vsind.th         : input of the stream flows of sink elements.
!      fort.95 fort.96  : output the maximum and average stream flow for
!      sourch and sink nodes.
!starting dates are set by users when asking for day, month, and year
!
!Author: Wei Huang
!
!
!
!ifort -O2 -CB -g -traceback -o coupling_nwm  ~/git/schism/src/Utility/UtilLib/julian_date.f90 ~/git/schism/src/Utility/UtilLib/schism_geometry.f90 ~/git/schism/src/Utility/UtilLib/pt_in_poly_test.f90 coupling_nwm.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -L$NETCDF/lib -lnetcdf -lnetcdff

       program coupling_nwm
       use netcdf
       use schism_geometry_mod
       use pt_in_poly_test
       use ifport
       implicit real*8(a-h,o-z)

       character(len=*),parameter::REPODir='/sciclone/home20/whuang07/&
                                             data10/NWM/CHRTOUT/'

       character(len=*),parameter::FILE_NAME='NWM_shp_ll.nc'
       integer :: ncid
       character(len =*), parameter :: lat_NAME = 'lat'
       character(len =*), parameter :: lon_NAME = 'lon'
       character(len =*), parameter :: featureID_NAME = 'featureID'
       character(len =*), parameter :: origID_NAME = 'ORIG_FID'
       character(len=154)::NWMfile
       character(len=2)::dd1,mm1
       character(len=1)::hh1,hh2
       character(len=4)::yy
       integer :: lat_varid, lon_varid, featureID_varid, origID_varid
       integer :: ne, np, id, id1, ele
       integer :: nID,recorddimid,latdimid,nwmdimid,n_nwm
       integer :: nwmID_varid
       integer,parameter::max_len=100000
       integer,allocatable :: nwmID(:)
       integer,allocatable ::featureID(:),origID(:),isbnd(:),nlnd(:),lndid(:,:)
       integer,allocatable::segn(:),seg_sink(:),seg_source(:),source_bnd(:),source_seg(:)
       integer,allocatable::bnd_sink(:),bnd_source(:),sink_bnd(:),sink_seg(:)
       integer,allocatable::eid(:),i34(:),elnode(:,:),ic3(:,:),elside(:,:),isdel(:,:),isidenode(:,:), nne(:),indel(:,:)
       integer,allocatable::uniso_ele(:),uniso_nwm(:),eleso_uni(:),nwmso_uni(:)
       integer,allocatable::unisi_ele(:),unisi_nwm(:),elesi_uni(:),nwmsi_uni(:)
       real*8,allocatable :: lats(:),lons(:),gx(:),gy(:),xcj(:),ycj(:)
       real*8,allocatable :: tstep(:)
       integer :: sf_varid,t_varid,n
       integer::ntime,nday,imm,idd,iyy,stime,timearray(9)
       integer,allocatable :: tm1(:),ty1(:),td1(:),td2(:),th1(:),th2(:)
       real*8::x11,x12,y11,y12,x21,x22,y21,y22,x3,y3
       real*8::arco(3),x4,y4,dup
       integer:: inside,nodel(3)
       real*8,allocatable :: salt(:),temp(:),streamflow(:)
       real*8,allocatable :: SF_so(:,:),SF_si(:,:),sflow(:)
       integer,parameter::nodata=-999900 !nodata value in nwm stream flows

       !Read in epsilon (nudging ratio from side center to downstream pt
       !of NWM)
       print*, 'Input nudging ratio (suggest 1.e-3):'
       read*, epsilon

       print*, 'Input rain_rate:'
       rain_rate=0.03 !m/hour

       print*, 'Input number of days'
       read*,nday
       ntime=nday*24

       print*, 'Enter start time - dd,mm,yyyy (e.g. 1 1 1992)'
       read(*,*) idd,imm,iyy

       open(14,file='hgrid.ll',status='old')!lambert projection, NWM has shorter precision
       read(14,*)
       read(14,*)ne,np
       write(*,*)'# of elements',ne,'# of nodes',np
       allocate(gx(np),gy(np),isbnd(np),nne(np))
       allocate(eid(ne),i34(ne),elnode(4,ne))
       do i=1,np
         read(14,*)j,gx(i),gy(i),dp
       enddo
       do i=1,ne
         read(14,*)eid(i),i34(i),elnode(1:i34(i),i)
       enddo
        !     Bnd
        isbnd=0
       read(14,*) nope
       read(14,*) neta
       ntot=0
       do k=1,nope
         read(14,*) nond
         do i=1,nond
           read(14,*) nd !iond(k,i)
           isbnd(nd)=k
         enddo !i
       enddo !k
       !write(*,*)gx(1),gy(1)
       !Land bnds
       read(14,*) nland !total # of land segments
       read(14,*) nvel !total # of nodes for all land boundaries
       if(nland<=0.or.nvel<=0) stop 'nland<=0.or.nvel<=0'
       allocate(nlnd(nland),lndid(nvel,nland))

       do k=1,nland
         read(14,*)nlnd(k) !total # of nodes on the land boundary
         if(nlnd(k)>nvel) stop 'check # of land nodes'
         do i=1,nlnd(k)
           read(14,*)lndid(i,k)
         enddo !i
       enddo !k
       close(14)

       write(*,*)'# of nodes on land boundary=',nlnd(:),nland

!     We also need elem ball
       nne=0
       do i=1,ne
         do j=1,i34(i)
           nd=elnode(j,i)
           nne(nd)=nne(nd)+1
         enddo
       enddo
       mnei=maxval(nne)

       allocate(indel(mnei,np),stat=istat)
       if(istat/=0) stop 'Failed to alloc. indel'
       nne=0
       do i=1,ne
         do j=1,i34(i)
           nd=elnode(j,i)
           nne(nd)=nne(nd)+1
           if(nne(nd)>mnei) then
             write(*,*)'Too many neighbors',nd
             stop
           endif
           indel(nne(nd),nd)=i
         enddo
       enddo !i

       ! Compute geometry
       call compute_nside(np,ne,i34,elnode(1:4,1:ne),ns2)
       allocate(ic3(4,ne),elside(4,ne),isdel(2,ns2),isidenode(2,ns2),xcj(ns2),ycj(ns2),stat=istat)
       if(istat/=0) stop 'Allocation error: side(0)'
       call schism_geometry_double(np,ne,ns2,gx,gy,i34,elnode(1:4,1:ne),ic3(1:4,1:ne), &
     &elside(1:4,1:ne),isdel,isidenode,xcj,ycj)

       call check(nf90_open(FILE_NAME, nf90_nowrite,ncid))
       !get the size of the featureID
       call check(nf90_inq_dimid(ncid,'RecordID', latdimid))
       call check(nf90_inquire_dimension(ncid, latdimid, len = nID))
       write(*,*)'# of NWM points',nID
       allocate(featureID(nID))
       allocate(origID(nID))
       allocate(lats(nID))
       allocate(lons(nID))
       call check(nf90_inq_varid(ncid,lat_NAME,lat_varid))
       call check(nf90_inq_varid(ncid,lon_NAME,lon_varid))
       call check(nf90_inq_varid(ncid,featureID_NAME,featureID_varid))
       call check(nf90_inq_varid(ncid,origID_NAME,origID_varid))

       call check(nf90_inq_varid(ncid,featureID_NAME,featureID_varid))
       call check(nf90_get_var(ncid,lat_varid,lats))
       call check(nf90_get_var(ncid,lon_varid,lons))
       call check(nf90_get_var(ncid,origID_varid,origID))
       call check(nf90_get_var(ncid,featureID_varid,featureID))
       call check(nf90_close(ncid))
       id=nID !size(lats) !# of vertices
       id1=maxval(origID)+1 !# of seg's
       write(*,*)'# NWM line segments=',id1
       !write(*,*)lons(1),lats(1)

       allocate(segn(id1)) !list of vertices along a seg
!       do i=1,id1+1
!         segn(i)=0
!         do j=1,id
!          if (origID(j).eq.i-1) then
!            segn(i)=segn(i)+1 ! get # of points on each line segments
!
!          endif
!         enddo
!       enddo

       segn=0
       do i=1,nID
         itmp=origID(i)+1 !seg ID
         if(itmp>id1) then
           print*, 'Overflow:',itmp,id1,i
           stop
         endif
         segn(itmp)=segn(itmp)+1
       enddo !i

       !p1 represents point pairs on NWM streams
       !p2 represents point pairs on grid land boudnary
       ninseg=0
       n_sink=0
       n_source=0
       allocate(seg_sink(max_len))
       allocate(bnd_sink(max_len))
       allocate(seg_source(max_len))
       allocate(bnd_source(max_len))

       do k1=1,nland !loop over all land segments of hgrid
       do k2=1,nlnd(k1)-1 !land seg
         nd1=lndid(k2,k1); nd2=lndid(k2+1,k1)
         x21=gx(nd1);  y21=gy(nd1)
         x22=gx(nd2); y22=gy(nd2)

         !Bounding box
         xmin2=min(x21,x22); xmax2=max(x21,x22)
         ymin2=min(y21,y22); ymax2=max(y21,y22)

        !write(*,*)p2(1)%x,p2(2)%x,p2(1)%y,p2(2)%y
        !Find bnd elem
         ele=0 !init as flag
         isd0=0
         loop1: do m=1,nne(nd1) !ball
                  ie=indel(m,nd1)
                  do mm=1,i34(ie)
                    isd=elside(mm,ie)
                    if((isidenode(1,isd)==nd1.and.isidenode(2,isd)==nd2).or.(isidenode(1,isd)==nd2.and.isidenode(2,isd)==nd1)) then
!                     nd3=sum(elnode(1:i34(ie),ie))-nd1-nd2
                      isd0=isd
                      ele=ie
                      exit loop1
                    endif
                  enddo !mm
                enddo loop1 !m
         if(ele==0.or.isd0==0) then
           write(*,*)'failed to find an elem'
           stop
         endif

         write(99,*) 'Found bnd elem:',k1,k2,nd1,nd2,ele,isd0

         dx2 = x21-x22
         dy2 = y21-y22
         if(dx2==0.) then
           a2=0.
           b2=y21
         else
           a2=dy2/dx2
           b2=y21-a2*x21
         endif
       !write(*,*)'a2,b2',a2,b2
         n=1
         do j=1,id1 !NWM seg

           do i=1,segn(j)-1 !vertices

             x11=lons(n+i-1);   y11=lats(n+i-1);!upstream point
             x12=lons(n+i); y12=lats(n+i);!downstream point
             !Bounding box
             xmin1=min(x11,x12); xmax1=max(x11,x12)
             ymin1=min(y11,y12); ymax1=max(y11,y12)

             !Skip bounding boxes that r outside each other
             if(xmin1>xmax2.or.xmin2>xmax1.or.ymin1>ymax2.or.ymin2>ymax1) cycle

             dx1 = x11-x12
             dy1 = y11-y12
             !write(*,*)p1(1)%x,p1(1)%y,dx1,dy1
             if(dx1==0.) then
               a1=0.
               b1=y11
             else
               a1=dy1/dx1 !dacter(len=*),parameterx1
               b1=y11-a1*x11
             endif
             x_inter=(b2-b1)/(a1-a2)
             y_inter=a1*x_inter+b1
!            write(*,*)x_inter,y_inter
             if(a1.ne.a2.and.x_inter.ge.min(x21,x22).and.x_inter.le.max(x21,x22).and.&
                         y_inter.ge.min(y21,y22).and.y_inter.le.max(y21,y22).and.&
                         x_inter.ge.min(x11,x12).and.x_inter.le.max(x11,x12).and.&
                         y_inter.ge.min(y11,y12).and.y_inter.le.max(y11,y12)) then
               ninseg=ninseg+1

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !--determine the source or sink by point in a poly-------!
        !--select the fourth point which is calculated by adding
        ! a short distance from side center towards the second point of NWM stream
              x4=xcj(isd0)*(1-epsilon)+epsilon*x12
              y4=ycj(isd0)*(1-epsilon)+epsilon*y12

          !    write(*,*)nd1,nd2,ele
              call pt_in_poly_double(i34(ele),gx(elnode(1:i34(ele),ele)),gy(elnode(1:i34(ele),ele)),x4,y4,inside,arco,nodel)

              if(inside.eq.1) then
                n_source=n_source+1
                if(n_source.gt.max_len) then
                  stop 'n_source exceeds max_len'
                endif
                seg_source(n_source)=featureID(n)
                bnd_source(n_source)=ele
              else
                n_sink=n_sink+1
                if(n_sink.gt.max_len) then
                  stop 'n_sink exceeds max_len'
                endif
                seg_sink(n_sink)=featureID(n)
                bnd_sink(n_sink)=ele
              endif

              write(98,*)featureID(n),ele,ninseg
            endif !inside bounding box

          enddo !i: vertices
       !write(*,*)n
          n=n+segn(j)
        enddo !j !NWM seg

      enddo !k2: each land seg
      enddo !k1: SCHISM land bnd seg

      write(98,*)'intersection#',ninseg,'sink#',n_sink,'source#',n_source

      allocate(sink_seg(max(1,n_sink)))
      allocate(sink_bnd(max(1,n_sink)))

      do i=1,n_sink
        sink_seg(i)=seg_sink(i)
        sink_bnd(i)=bnd_sink(i)
      enddo

      allocate(source_seg(max(1,n_source)))
      allocate(source_bnd(max(1,n_source)))

      do i=1,n_source
        source_seg(i)=seg_source(i)
        source_bnd(i)=bnd_source(i)
      enddo

      deallocate(seg_sink,bnd_sink,seg_source,bnd_source)
!count unique number of sink and source elements

      if (n_sink.ne.0) then
        allocate(unisi_ele(max_len))
        allocate(unisi_nwm(max_len))
        j=0;nsi=0;
        do i=1,n_sink-1
          if (sink_bnd(i).eq.sink_bnd(i+1)) then !remove duplicate
            j=j+1
          else
            nsi=nsi+1
            unisi_ele(nsi)=sink_bnd(i) !unique hgrid elem ID
            unisi_nwm(nsi)=sink_seg(i) !NWM feature ID
          endif
        enddo
        allocate(elesi_uni(nsi),nwmsi_uni(nsi))
        do i=1,nsi
          elesi_uni(i)=unisi_ele(i) !element id is unique now
          nwmsi_uni(i)=unisi_nwm(i) !nwm featureID is not necessarily unique
        enddo
      else
        allocate(unisi_ele(1))
        allocate(unisi_nwm(1))
        unisi_ele(1)=0
        unisi_nwm(1)=0
        nsi=0
      endif
      if (n_source.ne.0) then
        allocate(uniso_ele(max_len))
        allocate(uniso_nwm(max_len))
        j=0;nso=0;
        do i=1,n_source-1
          if (source_bnd(i).eq.source_bnd(i+1)) then !remove duplicate
            j=j+1
          else
            nso=nso+1
            uniso_ele(nso)=source_bnd(i)
            uniso_nwm(nso)=source_seg(i)
          endif
        enddo
        allocate(eleso_uni(nso),nwmso_uni(nso))
        do i=1,nso
          eleso_uni(i)=uniso_ele(i) !element id is unique now
          nwmso_uni(i)=uniso_nwm(i) !nwm featureID is not necessarily unique
        enddo
      else
        nso=0
        allocate(uniso_ele(1))
        allocate(uniso_nwm(1))
        uniso_ele(1)=0
        uniso_nwm(1)=0
      endif
      deallocate(uniso_ele,uniso_nwm)

!define time
      allocate(tstep(ntime))
      tstep(1)=0
      do i=2,ntime
        tstep(i)=3600*(i-1)
      enddo
      !write(*,*)tstep(1),tstep(2)
!write msource.th file
      allocate(salt(nso),temp(nso))
      do i=1,nso
         salt(i)=0
         temp(i)=-9999
      enddo
      open(15,file='msource.th')
      do i=1,ntime
      write(15,'(50000F15.3)')tstep(i),(temp(k),k=1,nso),(salt(k),k=1,nso)

      enddo

!following is to read in the NWM output
      allocate(tm1(ntime),td1(ntime),ty1(ntime),&
               th1(ntime),th2(ntime))
!convert time from julian date to calendar date:
! for function gmtime: Numeric time data to be formatted. Number of
! seconds since 00:00:00 Greenwich mean time, January 1, 1970
! The output of gmtime is One-dimensional array with 9 elements used to
! contain numeric time data.
! tarray(1)Seconds (0-61, where 60-61 can be returned for leap seconds)
! tarray(2)Minutes (0-59)
! tarray(3)Hours (0-23)
! tarray(4)Day of month (1-31)
! tarray(5)Month (0-11)
! tarray(6) Number of years since 1900
! tarray(7)Day of week (0-6, where 0 is Sunday)
! tarray(8)Day of year (0-365)
! tarray(9)Daylight saving flag (0 if standard time, 1 if daylight
! saving time)

      do i=1,ntime
         stime=julian_date(iyy,imm,idd)*24*3600
         stime=stime-julian_date(1970,1,1)*24*3600
         stime=stime+(i-1)*3600
         call gmtime(stime, timearray)
         tm1(i)=timearray(5)+1
         td1(i)=timearray(4)
         ty1(i)=timearray(6)+1900
      enddo

      do i=1,nday
         do j=24*(i-1)+1,24*i
           if((j-24*(i-1))<=10) then
             th1(j)=0
             th2(j)=int(j-24*(i-1))-1
           else
             th=int((j-24*(i-1)))-1
             th1(j)=int(th/10)
             th2(j)=th-10*th1(j)
           endif
         enddo
       enddo

!      allocate(SF_so(ntime,max(1,nso)))
!      allocate(SF_si(ntime,max(1,nsi)))
      if(n_source.ne.0) then
         allocate(SF_so(ntime,max(1,nso)))
      else
         allocate(SF_so(1,1))
      endif
      if(n_sink.ne.0) then
         allocate(SF_si(ntime,max(1,nsi)))
      else
         allocate(SF_si(1,1))
      endif

      do i=1,ntime
        write(yy,'(i4)')ty1(i)
        write(dd1,'(i2.2)')td1(i)
        write(mm1,'(i2.2)')tm1(i)
        write(hh1,'(i1)')th1(i)
        write(hh2,'(i1)')th2(i)

       ! print*, yy,dd1,mm1,hh1,hh2

        NWMfile=REPODir//yy//mm1//dd1&
                     //hh1//hh2//'00.CHRTOUT_DOMAIN1.comp'
        write(98,*)'Inside time iteration:',i,NWMfile

       call check(nf90_open(NWMfile, nf90_nowrite,ncid))
       call check(nf90_inq_dimid(ncid,'feature_id', nwmdimid))
       call check(nf90_inquire_dimension(ncid, nwmdimid, len = n_nwm))
       !write(*,*)'# of NWM output streams',n_nwm
       allocate(streamflow(n_nwm))
       allocate(nwmID(n_nwm))
       call check(nf90_inq_varid(ncid,'streamflow',sf_varid))
       call check(nf90_inq_varid(ncid,'time',t_varid))
       call check(nf90_inq_varid(ncid,'feature_id',nwmID_varid))
       call check(nf90_get_var(ncid,sf_varid,streamflow))
       call check(nf90_get_var(ncid,t_varid,time1))
       call check(nf90_get_var(ncid,nwmID_varid,nwmID))
       !write(*,*)nso,nsi,nwmso_uni(1)
       ilow=minval(nwmID); ihigh=maxval(nwmID)
       allocate(sflow(ilow:ihigh))
       write(98,*)'ilow,ihigh:', ilow,ihigh

       do k=1,n_nwm
         itmp=nwmID(k)
         if(itmp<ilow.or.itmp>ihigh) then
           print*, 'Overflow(2):',itmp,ilow,ihigh
           stop
         endif
         if (streamflow(k)==nodata)then
           streamflow(k)=0
         endif
         sflow(itmp)=0.01*streamflow(k)
         !if(nwmID(k)==2590185)then
         !  write(96,*)streamflow(k)
         !endif
       enddo !k
       !write(96,*)2590185,860,sflow(2590185)

!unique the element id during following process
!combine flux if there are muptiple streams going through one element

       if(n_source.ne.0) then
         dup=0;nso=0;
         do j=1,n_source-1
           if(source_seg(j)<ilow.or.source_seg(j)>ihigh) then
             write(97,*)'Overflow(3):',source_seg(j),ilow,ihigh
             SF_so(i,nso)=0
           !sum streams from same source elements
           else if (source_bnd(j).eq.source_bnd(j+1)) then
             dup=dup+sflow(source_seg(j))
           else
             nso=nso+1
             SF_so(i,nso)=dup+sflow(source_seg(j))
             dup=0
           endif
         enddo !j
       else
         SF_so(1,1)=0
       endif
       !write(*,*)'next step 1'
       if(n_sink.ne.0) then
         dup=0;nsi=0;
         do j=1,n_sink-1
           if(sink_seg(j)<ilow.or.sink_seg(j)>ihigh) then
             write(97,*)'Overflow(4):',sink_seg(j),ilow,ihigh
             SF_si(i,nsi)=0
           ! sum the flux from same sink elements
           else if (sink_bnd(j).eq.sink_bnd(j+1)) then
             dup=dup+sflow(sink_seg(j))
           else
             nsi=nsi+1
             SF_si(i,nsi)=dup+sflow(sink_seg(j))
             SF_si(i,nsi)=-SF_si(i,nsi)
             dup=0
           endif
         enddo !j
       else
         SF_si(1,1)=0
       endif
       !write(*,*)'next step 2'
       deallocate(streamflow,nwmID,sflow)
      enddo !i=1,ntime

      !print out statistics of streamflows
      if(nso.ne.0) then
        do i=1,nso
          write(96,*)maxval(SF_so(1:ntime,i),dim=1),sum(SF_so(1:ntime,i),dim=1)/ntime
        enddo
      endif
      if(nsi.ne.0) then
        do i=1,nsi
          write(95,*)maxval(SF_si(1:ntime,i),dim=1),sum(SF_si(1:ntime,i),dim=1)/ntime
        enddo
      endif

!write source_sink.in file
      open(13,file='source_sink.in')
      if (n_source.ne.0) then
        write(13,*)nso
        do i=1,nso
           write(13,*)eleso_uni(i)
        enddo
      else
        write(13,*)0
      endif
      write(13,*)
      if (n_sink.ne.0) then
        write(13,*)nsi
        do i=1,nsi
           write(13,*)elesi_uni(i)
        enddo
      else
        write(13,*)0
      endif

!write vsource.th file
      if (n_source.ne.0) then
        open(16,file='vsource.th',status='replace')
        do i=1,ntime
          write(16,'(10000F15.3)')tstep(i),(SF_so(i,k),k=1,nso)
        end do
      endif
!      write(*,*)'next step 4'
!write vsink.th
      if (n_sink.ne.0) then
        open(17,file='vsink.th')
        do i=1,ntime
          write(17,'(10000F15.3)')tstep(i),(SF_si(i,k),k=1,nsi)
        end do
      endif
!      write(*,*)'next step 5'

      call check(nf90_close(ncid))

      contains
         subroutine check(status)
          integer, intent ( in) :: status

          if(status /= nf90_noerr) then
          print *, trim(nf90_strerror(status))
          stop "Stopped"
          end if
         end subroutine check
      end


text:
      SUBROUTINE DISKS(RHO,TH,PH,EPS,FV)
      DOUBLE PRECISION RP(3),ALI(3,3),PHII(3,3),AI(3,3),PLI(3,3),RI(3),
     1ALPT(3,3),FV(3),CO,    RM2,AI1,AKI,SKI,AT,H,V,FM(3,3),VF(3)
      DOUBLE PRECISION TH, PH, EPS, U, RHO, RC, COM(1), THRAD, PHRAD,
     1 RCOS, DLCO, PLCO, RMAI,UI,A,PHI,DL,RH
      COMMON /PRAM/ NMAS, NHARM, NDISK
      COMMON RC,U,COM
      DATA CO /.1745329251994430D-1/
COMMENTS ABOUT CHANGES TO STATEMENT FUNCTIONS:
C... I1 - REPRESENTS THE COUNTER VARIABLE I IN LOOP 100
C... I2 - REPRESENTS THE NUMBER OF DISKS(NDISK) TO PROCESS
C... I3 - REPRESENTS THE ACTUAL ADDRESS IN THE COM ARRAY
C
       UI(I1,I3)=COM(I1+I3)
       A(I1,I2,I3)=COM(I1+I2+I3)
       PHI(I1,I2,I3)=COM(I1+2*I2+I3)
       DL(I1,I2,I3)=COM(I1+3*I2+I3)
       RH(I1,I2,I3)=COM(I1+4*I2+I3)
      IA = 4*NMAS
      ID = IA + 2*NHARM**2 + NHARM
      THRAD = TH*CO
      PHRAD = PH*CO
      RCOS = RHO*DCOS(THRAD)
      RP(1) = RCOS*DCOS(PHRAD)
      RP(2) = RCOS*DSIN(PHRAD)
      RP(3) = RHO*DSIN(THRAD)
      DO 1 I=1,3
      FV(I)=0.D0
      DO 1 J=1,3
      ALI(I,J)=0.D0
      PHII(I,J)=0.D0
    1 AI(I,J)=0.D0
      PHII(2,2)=1.D0
      ALI(3,3)=1.D0
      DO 100 I=1,NDISK
      DLCO = DL(I,NDISK,ID)*CO
      ALI(1,1) = DCOS(DLCO)
      ALI(2,2)=ALI(1,1)
      ALI(1,2) = DSIN(DLCO)
      ALI(2,1)=-ALI(1,2)
      PLCO = PHI(I,NDISK,ID)*CO
      PHII(1,1) = DSIN(PLCO)
      PHII(3,3)=PHII(1,1)
      PHII(3,1) = DCOS(PLCO)
      PHII(1,3)=-PHII(3,1)
      CALL MATMUL(PHII,ALI,3,3,3,PLI)
      CALL MATVEC(PLI,RP,3,3,RI)
      RI(3)=RI(3)-RH(I,NDISK,ID)
      RM2=RI(1)*RI(1)+RI(2)*RI(2)+RI(3)*RI(3)
      IF(RM2/UI(I,ID).GT.EPS)GO TO 100
      AI1 = 1.0D0/A(I,NDISK,ID)**2
      RMAI = RM2*AI1 - 1.0D0
      AKI = (RMAI + DSQRT(RMAI**2 + 4.0D0*RI(3)**2*AI1))*0.5D0
      SKI=DSQRT(AKI)
      AT=DATAN(1.D0/SKI)
      H=AT-SKI/(1.D0+AKI)
      V=2.D0/SKI-2.D0*AT
      AI(1,1)=H
      AI(2,2)=H
      AI(3,3)=V
      DO 2 J=1,3
      DO 2 K=1,3
    2 ALPT(J,K)=PLI(K,J)
      CALL MATMUL(ALPT,AI,3,3,3,FM)
      CALL MATVEC(FM,RI,3,3,VF)
      DO 3 J=1,3
    3 FV(J)=VF(J)*3.D0*U*UI(I,ID)*AI1/(2.D0*A(I,NDISK,ID))+FV(J)
  100 CONTINUE
      RETURN
      END

ext: for