text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1995  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine estrtor2  --  atomwise stretch-torsion Hessian  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "estrtor2" calculates the stretch-torsion potential energy
c     second derivatives with respect to Cartesian coordinates
c
c
      subroutine estrtor2 (i)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bond.i'
      include 'bound.i'
      include 'group.i'
      include 'hessn.i'
      include 'strtor.i'
      include 'torpot.i'
      include 'tors.i'
      integer i,j,k,istrtor
      integer ia,ib,ic,id
      real*8 dedphi,d2edphi2,fgrp
      real*8 rt2,ru2,rtru,rcb
      real*8 xt,yt,zt,xu,yu,zu
      real*8 xtu,ytu,ztu
      real*8 v1,v2,v3
      real*8 c1,c2,c3
      real*8 s1,s2,s3
      real*8 sine,cosine
      real*8 sine2,cosine2
      real*8 sine3,cosine3
      real*8 xia,yia,zia
      real*8 xib,yib,zib
      real*8 xic,yic,zic
      real*8 xid,yid,zid
      real*8 xba,yba,zba
      real*8 xcb,ycb,zcb
      real*8 xdc,ydc,zdc
      real*8 xca,yca,zca
      real*8 xdb,ydb,zdb
      real*8 phi1,phi2,phi3
      real*8 dphi1,dphi2,dphi3
      real*8 d2phi1,d2phi2,d2phi3
      real*8 dr,ddr,d2dr
      real*8 ddrdx,ddrdy,ddrdz
      real*8 d2drdxx,d2drdyy,d2drdzz
      real*8 d2drdxy,d2drdxz,d2drdyz
      real*8 dphidxt,dphidyt,dphidzt
      real*8 dphidxu,dphidyu,dphidzu
      real*8 dphidxia,dphidyia,dphidzia
      real*8 dphidxib,dphidyib,dphidzib
      real*8 dphidxic,dphidyic,dphidzic
      real*8 dphidxid,dphidyid,dphidzid
      real*8 xycb2,xzcb2,yzcb2
      real*8 rcbxt,rcbyt,rcbzt,rcbt2
      real*8 rcbxu,rcbyu,rcbzu,rcbu2
      real*8 dphidxibt,dphidyibt,dphidzibt
      real*8 dphidxibu,dphidyibu,dphidzibu
      real*8 dphidxict,dphidyict,dphidzict
      real*8 dphidxicu,dphidyicu,dphidzicu
      real*8 dxia,dyia,dzia
      real*8 dxib,dyib,dzib
      real*8 dxic,dyic,dzic
      real*8 dxid,dyid,dzid
      real*8 dxiaxia,dyiayia,dziazia
      real*8 dxibxib,dyibyib,dzibzib
      real*8 dxicxic,dyicyic,dziczic
      real*8 dxidxid,dyidyid,dzidzid
      real*8 dxiayia,dxiazia,dyiazia
      real*8 dxibyib,dxibzib,dyibzib
      real*8 dxicyic,dxiczic,dyiczic
      real*8 dxidyid,dxidzid,dyidzid
      real*8 dxiaxib,dxiayib,dxiazib
      real*8 dyiaxib,dyiayib,dyiazib
      real*8 dziaxib,dziayib,dziazib
      real*8 dxiaxic,dxiayic,dxiazic
      real*8 dyiaxic,dyiayic,dyiazic
      real*8 dziaxic,dziayic,dziazic
      real*8 dxiaxid,dxiayid,dxiazid
      real*8 dyiaxid,dyiayid,dyiazid
      real*8 dziaxid,dziayid,dziazid
      real*8 dxibxic,dxibyic,dxibzic
      real*8 dyibxic,dyibyic,dyibzic
      real*8 dzibxic,dzibyic,dzibzic
      real*8 dxibxid,dxibyid,dxibzid
      real*8 dyibxid,dyibyid,dyibzid
      real*8 dzibxid,dzibyid,dzibzid
      real*8 dxicxid,dxicyid,dxiczid
      real*8 dyicxid,dyicyid,dyiczid
      real*8 dzicxid,dzicyid,dziczid
      logical proceed
c
c
c     compute the Hessian elements of the stretch-torsions
c
      do istrtor = 1, nstrtor
         j = ist(1,istrtor)
         ia = itors(1,j)
         ib = itors(2,j)
         ic = itors(3,j)
         id = itors(4,j)
c
c     decide whether to compute the current interaction
c
         proceed = (i.eq.ia .or. i.eq.ib .or. i.eq.ic .or. i.eq.id)
         if (proceed .and. use_group)
     &      call groups (proceed,fgrp,ia,ib,ic,id,0,0)
c
c     compute the value of the torsional angle
c
         if (proceed) then
            xia = x(ia)
            yia = y(ia)
            zia = z(ia)
            xib = x(ib)
            yib = y(ib)
            zib = z(ib)
            xic = x(ic)
            yic = y(ic)
            zic = z(ic)
            xid = x(id)
            yid = y(id)
            zid = z(id)
            xba = xib - xia
            yba = yib - yia
            zba = zib - zia
            xcb = xic - xib
            ycb = yic - yib
            zcb = zic - zib
            xdc = xid - xic
            ydc = yid - yic
            zdc = zid - zic
            if (use_polymer) then
               call image (xba,yba,zba)
               call image (xcb,ycb,zcb)
               call image (xdc,ydc,zdc)
            end if
            xt = yba*zcb - ycb*zba
            yt = zba*xcb - zcb*xba
            zt = xba*ycb - xcb*yba
            xu = ycb*zdc - ydc*zcb
            yu = zcb*xdc - zdc*xcb
            zu = xcb*ydc - xdc*ycb
            xtu = yt*zu - yu*zt
            ytu = zt*xu - zu*xt
            ztu = xt*yu - xu*yt
            rt2 = xt*xt + yt*yt + zt*zt
            ru2 = xu*xu + yu*yu + zu*zu
            rtru = sqrt(rt2 * ru2)
            if (rtru .ne. 0.0d0) then
               rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb)
               cosine = (xt*xu + yt*yu + zt*zu) / rtru
               sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru)
c
c     set the stretch-torsional parameters for this angle
c
               v1 = kst(1,istrtor)
               c1 = tors1(3,i)
               s1 = tors1(4,i)
               v2 = kst(2,istrtor)
               c2 = tors2(3,i)
               s2 = tors2(4,i)
               v3 = kst(3,istrtor)
               c3 = tors3(3,i)
               s3 = tors3(4,i)
c
c     compute the multiple angle trigonometry and the phase terms
c
               cosine2 = cosine*cosine - sine*sine
               sine2 = 2.0d0 * cosine * sine
               cosine3 = cosine*cosine2 - sine*sine2
               sine3 = cosine*sine2 + sine*cosine2
               phi1 = 1.0d0 + (cosine*c1 + sine*s1)
               phi2 = 1.0d0 + (cosine2*c2 + sine2*s2)
               phi3 = 1.0d0 + (cosine3*c3 + sine3*s3)
               dphi1 = (cosine*s1 - sine*c1)
               dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2)
               dphi3 = 3.0d0 * (cosine3*s3 - sine3*c3)
               d2phi1 = -(cosine*c1 + sine*s1)
               d2phi2 = -4.0d0 * (cosine2*c2 + sine2*s2)
               d2phi3 = -9.0d0 * (cosine3*c3 + sine3*s3)
c
c     calculate the bond-stretch for the central bond
c
               k = ist(2,istrtor)
               dr = rcb - bl(k)
c
c     calculate the stretch-torsion master chain rule terms
c
               dedphi = storunit * (v1*dphi1 + v2*dphi2 + v3*dphi3)
               d2edphi2 = storunit * dr
     &                       * (v1*d2phi1 + v2*d2phi2 + v3*d2phi3)
               ddr = 1.0d0 / rcb
               d2dr = -storunit * (v1*phi1 + v2*phi2 + v3*phi3) / rcb**3
c
c     scale the interaction based on its group membership
c
               if (use_group) then
                  dedphi = dedphi * fgrp
                  d2edphi2 = d2edphi2 * fgrp
                  d2dr = d2dr * fgrp
               end if
c
c     first and second derivative components for the bond stretch
c
               ddrdx = xcb * ddr
               ddrdy = ycb * ddr
               ddrdz = zcb * ddr
               d2drdxx = (xcb**2-rcb**2) * d2dr
               d2drdyy = (ycb**2-rcb**2) * d2dr
               d2drdzz = (zcb**2-rcb**2) * d2dr
               d2drdxy = xcb * ycb * d2dr
               d2drdxz = xcb * zcb * d2dr
               d2drdyz = ycb * zcb * d2dr
c
c     abbreviations for first derivative chain rule terms
c
               xca = xic - xia
               yca = yic - yia
               zca = zic - zia
               xdb = xid - xib
               ydb = yid - yib
               zdb = zid - zib
               if (use_polymer) then
                  call image (xca,yca,zca)
                  call image (xdb,ydb,zdb)
               end if
               dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb)
               dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb)
               dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb)
               dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb)
               dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb)
               dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb)
c
c     abbreviations for second derivative chain rule terms
c
               xycb2 = xcb*xcb + ycb*ycb
               xzcb2 = xcb*xcb + zcb*zcb
               yzcb2 = ycb*ycb + zcb*zcb
               rcbxt = -2.0d0 * rcb * dphidxt
               rcbyt = -2.0d0 * rcb * dphidyt
               rcbzt = -2.0d0 * rcb * dphidzt
               rcbt2 = rcb * rt2
               rcbxu = 2.0d0 * rcb * dphidxu
               rcbyu = 2.0d0 * rcb * dphidyu
               rcbzu = 2.0d0 * rcb * dphidzu
               rcbu2 = rcb * ru2
               dphidxibt = yca*dphidzt - zca*dphidyt
               dphidxibu = zdc*dphidyu - ydc*dphidzu
               dphidyibt = zca*dphidxt - xca*dphidzt
               dphidyibu = xdc*dphidzu - zdc*dphidxu
               dphidzibt = xca*dphidyt - yca*dphidxt
               dphidzibu = ydc*dphidxu - xdc*dphidyu
               dphidxict = zba*dphidyt - yba*dphidzt
               dphidxicu = ydb*dphidzu - zdb*dphidyu
               dphidyict = xba*dphidzt - zba*dphidxt
               dphidyicu = zdb*dphidxu - xdb*dphidzu
               dphidzict = yba*dphidxt - xba*dphidyt
               dphidzicu = xdb*dphidyu - ydb*dphidxu
c
c     intermediate terms for first derivative components
c
               dphidxia = zcb*dphidyt - ycb*dphidzt
               dphidyia = xcb*dphidzt - zcb*dphidxt
               dphidzia = ycb*dphidxt - xcb*dphidyt
               dphidxib = dphidxibt + dphidxibu
               dphidyib = dphidyibt + dphidyibu
               dphidzib = dphidzibt + dphidzibu
               dphidxic = dphidxict + dphidxicu
               dphidyic = dphidyict + dphidyicu
               dphidzic = dphidzict + dphidzicu
               dphidxid = zcb*dphidyu - ycb*dphidzu
               dphidyid = xcb*dphidzu - zcb*dphidxu
               dphidzid = ycb*dphidxu - xcb*dphidyu
c
c     chain rule terms for first derivative components
c
               dxia = dedphi * dphidxia
               dyia = dedphi * dphidyia
               dzia = dedphi * dphidzia
               dxib = dedphi * dphidxib
               dyib = dedphi * dphidyib
               dzib = dedphi * dphidzib
               dxic = dedphi * dphidxic
               dyic = dedphi * dphidyic
               dzic = dedphi * dphidzic
               dxid = dedphi * dphidxid
               dyid = dedphi * dphidyid
               dzid = dedphi * dphidzid
               dedphi = dedphi * dr
c
c     chain rule terms for second derivative components
c
               dxiaxia = rcbxt*dphidxia
               dxiayia = rcbxt*dphidyia - zcb*rcb/rt2
               dxiazia = rcbxt*dphidzia + ycb*rcb/rt2
               dxiaxic = rcbxt*dphidxict + xcb*xt/rcbt2
               dxiayic = rcbxt*dphidyict - dphidzt
     &                      - (xba*zcb*xcb+zba*yzcb2)/rcbt2
               dxiazic = rcbxt*dphidzict + dphidyt
     &                      + (xba*ycb*xcb+yba*yzcb2)/rcbt2
               dxiaxid = 0.0d0
               dxiayid = 0.0d0
               dxiazid = 0.0d0
               dyiayia = rcbyt*dphidyia
               dyiazia = rcbyt*dphidzia - xcb*rcb/rt2
               dyiaxib = rcbyt*dphidxibt - dphidzt
     &                      - (yca*zcb*ycb+zca*xzcb2)/rcbt2
               dyiaxic = rcbyt*dphidxict + dphidzt
     &                      + (yba*zcb*ycb+zba*xzcb2)/rcbt2
               dyiayic = rcbyt*dphidyict + ycb*yt/rcbt2
               dyiazic = rcbyt*dphidzict - dphidxt
     &                      - (yba*xcb*ycb+xba*xzcb2)/rcbt2
               dyiaxid = 0.0d0
               dyiayid = 0.0d0
               dyiazid = 0.0d0
               dziazia = rcbzt*dphidzia
               dziaxib = rcbzt*dphidxibt + dphidyt
     &                      + (zca*ycb*zcb+yca*xycb2)/rcbt2
               dziayib = rcbzt*dphidyibt - dphidxt
     &                      - (zca*xcb*zcb+xca*xycb2)/rcbt2
               dziaxic = rcbzt*dphidxict - dphidyt
     &                      - (zba*ycb*zcb+yba*xycb2)/rcbt2
               dziayic = rcbzt*dphidyict + dphidxt
     &                      + (zba*xcb*zcb+xba*xycb2)/rcbt2
               dziazic = rcbzt*dphidzict + zcb*zt/rcbt2
               dziaxid = 0.0d0
               dziayid = 0.0d0
               dziazid = 0.0d0
               dxibxic = -xcb*dphidxib/(rcb*rcb)
     &             - (yca*(zba*xcb+yt)-zca*(yba*xcb-zt))/rcbt2
     &             - 2.0d0*(yt*zba-yba*zt)*dphidxibt/rt2
     &             - (zdc*(ydb*xcb+zu)-ydc*(zdb*xcb-yu))/rcbu2
     &             + 2.0d0*(yu*zdb-ydb*zu)*dphidxibu/ru2
               dxibyic = -ycb*dphidxib/(rcb*rcb) + dphidzt + dphidzu
     &             - (yca*(zba*ycb-xt)+zca*(xba*xcb+zcb*zba))/rcbt2
     &             - 2.0d0*(zt*xba-zba*xt)*dphidxibt/rt2
     &             + (zdc*(xdb*xcb+zcb*zdb)+ydc*(zdb*ycb+xu))/rcbu2
     &             + 2.0d0*(zu*xdb-zdb*xu)*dphidxibu/ru2
               dxibxid = rcbxu*dphidxibu + xcb*xu/rcbu2
               dxibyid = rcbyu*dphidxibu - dphidzu
     &                      - (ydc*zcb*ycb+zdc*xzcb2)/rcbu2
               dxibzid = rcbzu*dphidxibu + dphidyu
     &                      + (zdc*ycb*zcb+ydc*xycb2)/rcbu2
               dyibzib = ycb*dphidzib/(rcb*rcb)
     &             - (xca*(xca*xcb+zcb*zca)+yca*(ycb*xca+zt))/rcbt2
     &             - 2.0d0*(xt*zca-xca*zt)*dphidzibt/rt2
     &             + (ydc*(xdc*ycb-zu)+xdc*(xdc*xcb+zcb*zdc))/rcbu2
     &             + 2.0d0*(xu*zdc-xdc*zu)*dphidzibu/ru2
               dyibxic = -xcb*dphidyib/(rcb*rcb) - dphidzt - dphidzu
     &             + (xca*(zba*xcb+yt)+zca*(zba*zcb+ycb*yba))/rcbt2
     &             - 2.0d0*(yt*zba-yba*zt)*dphidyibt/rt2
     &             - (zdc*(zdb*zcb+ycb*ydb)+xdc*(zdb*xcb-yu))/rcbu2
     &             + 2.0d0*(yu*zdb-ydb*zu)*dphidyibu/ru2
               dyibyic = -ycb*dphidyib/(rcb*rcb)
     &             - (zca*(xba*ycb+zt)-xca*(zba*ycb-xt))/rcbt2
     &             - 2.0d0*(zt*xba-zba*xt)*dphidyibt/rt2
     &             - (xdc*(zdb*ycb+xu)-zdc*(xdb*ycb-zu))/rcbu2
     &             + 2.0d0*(zu*xdb-zdb*xu)*dphidyibu/ru2
               dyibxid = rcbxu*dphidyibu + dphidzu
     &                      + (xdc*zcb*xcb+zdc*yzcb2)/rcbu2
               dyibyid = rcbyu*dphidyibu + ycb*yu/rcbu2
               dyibzid = rcbzu*dphidyibu - dphidxu
     &                      - (zdc*xcb*zcb+xdc*xycb2)/rcbu2
               dzibxic = -xcb*dphidzib/(rcb*rcb) + dphidyt + dphidyu
     &             - (xca*(yba*xcb-zt)+yca*(zba*zcb+ycb*yba))/rcbt2
     &             - 2.0d0*(yt*zba-yba*zt)*dphidzibt/rt2
     &             + (ydc*(zdb*zcb+ycb*ydb)+xdc*(ydb*xcb+zu))/rcbu2
     &             + 2.0d0*(yu*zdb-ydb*zu)*dphidzibu/ru2
               dzibzic = -zcb*dphidzib/(rcb*rcb)
     &             - (xca*(yba*zcb+xt)-yca*(xba*zcb-yt))/rcbt2
     &             - 2.0d0*(xt*yba-xba*yt)*dphidzibt/rt2
     &             - (ydc*(xdb*zcb+yu)-xdc*(ydb*zcb-xu))/rcbu2
     &             + 2.0d0*(xu*ydb-xdb*yu)*dphidzibu/ru2
               dzibxid = rcbxu*dphidzibu - dphidyu
     &                      - (xdc*ycb*xcb+ydc*yzcb2)/rcbu2
               dzibyid = rcbyu*dphidzibu + dphidxu
     &                      + (ydc*xcb*ycb+xdc*xzcb2)/rcbu2
               dzibzid = rcbzu*dphidzibu + zcb*zu/rcbu2
               dxicxid = rcbxu*dphidxicu - xcb*(zdb*ycb-ydb*zcb)/rcbu2
               dxicyid = rcbyu*dphidxicu + dphidzu
     &                      + (ydb*zcb*ycb+zdb*xzcb2)/rcbu2
               dxiczid = rcbzu*dphidxicu - dphidyu
     &                      - (zdb*ycb*zcb+ydb*xycb2)/rcbu2
               dyicxid = rcbxu*dphidyicu - dphidzu
     &                      - (xdb*zcb*xcb+zdb*yzcb2)/rcbu2
               dyicyid = rcbyu*dphidyicu - ycb*(xdb*zcb-zdb*xcb)/rcbu2
               dyiczid = rcbzu*dphidyicu + dphidxu
     &                      + (zdb*xcb*zcb+xdb*xycb2)/rcbu2
               dzicxid = rcbxu*dphidzicu + dphidyu
     &                      + (xdb*ycb*xcb+ydb*yzcb2)/rcbu2
               dzicyid = rcbyu*dphidzicu - dphidxu
     &                      - (ydb*xcb*ycb+xdb*xzcb2)/rcbu2
               dziczid = rcbzu*dphidzicu - zcb*(ydb*xcb-xdb*ycb)/rcbu2
               dxidxid = rcbxu*dphidxid
               dxidyid = rcbxu*dphidyid + zcb*rcb/ru2
               dxidzid = rcbxu*dphidzid - ycb*rcb/ru2
               dyidyid = rcbyu*dphidyid
               dyidzid = rcbyu*dphidzid + xcb*rcb/ru2
               dzidzid = rcbzu*dphidzid
c
c     get some second derivative chain rule terms by difference
c
               dxiaxib = -dxiaxia - dxiaxic - dxiaxid
               dxiayib = -dxiayia - dxiayic - dxiayid
               dxiazib = -dxiazia - dxiazic - dxiazid
               dyiayib = -dyiayia - dyiayic - dyiayid
               dyiazib = -dyiazia - dyiazic - dyiazid
               dziazib = -dziazia - dziazic - dziazid
               dxibxib = -dxiaxib - dxibxic - dxibxid
               dxibyib = -dyiaxib - dxibyic - dxibyid
               dxibzib = -dxiazib - dzibxic - dzibxid
               dxibzic = -dziaxib - dxibzib - dxibzid
               dyibyib = -dyiayib - dyibyic - dyibyid
               dyibzic = -dziayib - dyibzib - dyibzid
               dzibzib = -dziazib - dzibzic - dzibzid
               dzibyic = -dyiazib - dyibzib - dzibyid
               dxicxic = -dxiaxic - dxibxic - dxicxid
               dxicyic = -dyiaxic - dyibxic - dxicyid
               dxiczic = -dziaxic - dzibxic - dxiczid
               dyicyic = -dyiayic - dyibyic - dyicyid
               dyiczic = -dziayic - dzibyic - dyiczid
               dziczic = -dziazic - dzibzic - dziczid
c
c     increment diagonal and off-diagonal Hessian elements
c
               if (i .eq. ia) then
                  hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxia
     &                             + d2edphi2*dphidxia*dphidxia
                  hessy(1,ia) = hessy(1,ia) + dedphi*dxiayia
     &                             + d2edphi2*dphidxia*dphidyia
                  hessz(1,ia) = hessz(1,ia) + dedphi*dxiazia
     &                             + d2edphi2*dphidxia*dphidzia
                  hessx(2,ia) = hessx(2,ia) + dedphi*dxiayia
     &                             + d2edphi2*dphidxia*dphidyia
                  hessy(2,ia) = hessy(2,ia) + dedphi*dyiayia
     &                             + d2edphi2*dphidyia*dphidyia
                  hessz(2,ia) = hessz(2,ia) + dedphi*dyiazia
     &                             + d2edphi2*dphidyia*dphidzia
                  hessx(3,ia) = hessx(3,ia) + dedphi*dxiazia
     &                             + d2edphi2*dphidxia*dphidzia
                  hessy(3,ia) = hessy(3,ia) + dedphi*dyiazia
     &                             + d2edphi2*dphidyia*dphidzia
                  hessz(3,ia) = hessz(3,ia) + dedphi*dziazia
     &                             + d2edphi2*dphidzia*dphidzia
                  hessx(1,ib) = hessx(1,ib) + dedphi*dxiaxib
     &                             + d2edphi2*dphidxia*dphidxib
     &                             - dxia*ddrdx
                  hessy(1,ib) = hessy(1,ib) + dedphi*dyiaxib
     &                             + d2edphi2*dphidyia*dphidxib
     &                             - dyia*ddrdx
                  hessz(1,ib) = hessz(1,ib) + dedphi*dziaxib
     &                             + d2edphi2*dphidzia*dphidxib
     &                             - dzia*ddrdx
                  hessx(2,ib) = hessx(2,ib) + dedphi*dxiayib
     &                             + d2edphi2*dphidxia*dphidyib
     &                             - dxia*ddrdy
                  hessy(2,ib) = hessy(2,ib) + dedphi*dyiayib
     &                             + d2edphi2*dphidyia*dphidyib
     &                             - dyia*ddrdy
                  hessz(2,ib) = hessz(2,ib) + dedphi*dziayib
     &                             + d2edphi2*dphidzia*dphidyib
     &                             - dzia*ddrdy
                  hessx(3,ib) = hessx(3,ib) + dedphi*dxiazib
     &                             + d2edphi2*dphidxia*dphidzib
     &                             - dxia*ddrdz
                  hessy(3,ib) = hessy(3,ib) + dedphi*dyiazib
     &                             + d2edphi2*dphidyia*dphidzib
     &                             - dyia*ddrdz
                  hessz(3,ib) = hessz(3,ib) + dedphi*dziazib
     &                             + d2edphi2*dphidzia*dphidzib
     &                             - dzia*ddrdz
                  hessx(1,ic) = hessx(1,ic) + dedphi*dxiaxic
     &                             + d2edphi2*dphidxia*dphidxic
     &                             + dxia*ddrdx
                  hessy(1,ic) = hessy(1,ic) + dedphi*dyiaxic
     &                             + d2edphi2*dphidyia*dphidxic
     &                             + dyia*ddrdx
                  hessz(1,ic) = hessz(1,ic) + dedphi*dziaxic
     &                             + d2edphi2*dphidzia*dphidxic
     &                             + dzia*ddrdx
                  hessx(2,ic) = hessx(2,ic) + dedphi*dxiayic
     &                             + d2edphi2*dphidxia*dphidyic
     &                             + dxia*ddrdy
                  hessy(2,ic) = hessy(2,ic) + dedphi*dyiayic
     &                             + d2edphi2*dphidyia*dphidyic
     &                             + dyia*ddrdy
                  hessz(2,ic) = hessz(2,ic) + dedphi*dziayic
     &                             + d2edphi2*dphidzia*dphidyic
     &                             + dzia*ddrdy
                  hessx(3,ic) = hessx(3,ic) + dedphi*dxiazic
     &                             + d2edphi2*dphidxia*dphidzic
     &                             + dxia*ddrdz
                  hessy(3,ic) = hessy(3,ic) + dedphi*dyiazic
     &                             + d2edphi2*dphidyia*dphidzic
     &                             + dyia*ddrdz
                  hessz(3,ic) = hessz(3,ic) + dedphi*dziazic
     &                             + d2edphi2*dphidzia*dphidzic
     &                             + dzia*ddrdz
                  hessx(1,id) = hessx(1,id) + dedphi*dxiaxid
     &                             + d2edphi2*dphidxia*dphidxid
                  hessy(1,id) = hessy(1,id) + dedphi*dyiaxid
     &                             + d2edphi2*dphidyia*dphidxid
                  hessz(1,id) = hessz(1,id) + dedphi*dziaxid
     &                             + d2edphi2*dphidzia*dphidxid
                  hessx(2,id) = hessx(2,id) + dedphi*dxiayid
     &                             + d2edphi2*dphidxia*dphidyid
                  hessy(2,id) = hessy(2,id) + dedphi*dyiayid
     &                             + d2edphi2*dphidyia*dphidyid
                  hessz(2,id) = hessz(2,id) + dedphi*dziayid
     &                             + d2edphi2*dphidzia*dphidyid
                  hessx(3,id) = hessx(3,id) + dedphi*dxiazid
     &                             + d2edphi2*dphidxia*dphidzid
                  hessy(3,id) = hessy(3,id) + dedphi*dyiazid
     &                             + d2edphi2*dphidyia*dphidzid
                  hessz(3,id) = hessz(3,id) + dedphi*dziazid
     &                             + d2edphi2*dphidzia*dphidzid
               else if (i .eq. ib) then
                  hessx(1,ib) = hessx(1,ib) + dedphi*dxibxib
     &                             + d2edphi2*dphidxib*dphidxib
     &                             - 2.0d0*dxib*ddrdx + d2drdxx
                  hessy(1,ib) = hessy(1,ib) + dedphi*dxibyib
     &                             + d2edphi2*dphidxib*dphidyib
     &                             - dyib*ddrdx - dxib*ddrdy + d2drdxy
                  hessz(1,ib) = hessz(1,ib) + dedphi*dxibzib
     &                             + d2edphi2*dphidxib*dphidzib
     &                             - dzib*ddrdx - dxib*ddrdz + d2drdxz
                  hessx(2,ib) = hessx(2,ib) + dedphi*dxibyib
     &                             + d2edphi2*dphidxib*dphidyib
     &                             - dxib*ddrdy - dyib*ddrdx + d2drdxy
                  hessy(2,ib) = hessy(2,ib) + dedphi*dyibyib
     &                             + d2edphi2*dphidyib*dphidyib
     &                             - 2.0d0*dyib*ddrdy + d2drdyy
                  hessz(2,ib) = hessz(2,ib) + dedphi*dyibzib
     &                             + d2edphi2*dphidyib*dphidzib
     &                             - dzib*ddrdy - dyib*ddrdz + d2drdyz
                  hessx(3,ib) = hessx(3,ib) + dedphi*dxibzib
     &                             + d2edphi2*dphidxib*dphidzib
     &                             - dxib*ddrdz - dzib*ddrdx + d2drdxz
                  hessy(3,ib) = hessy(3,ib) + dedphi*dyibzib
     &                             + d2edphi2*dphidyib*dphidzib
     &                             - dyib*ddrdz - dzib*ddrdy + d2drdyz
                  hessz(3,ib) = hessz(3,ib) + dedphi*dzibzib
     &                             + d2edphi2*dphidzib*dphidzib
     &                             - 2.0d0*dzib*ddrdz + d2drdzz
                  hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxib
     &                             + d2edphi2*dphidxib*dphidxia
     &                             - dxia*ddrdx
                  hessy(1,ia) = hessy(1,ia) + dedphi*dxiayib
     &                             + d2edphi2*dphidyib*dphidxia
     &                             - dxia*ddrdy
                  hessz(1,ia) = hessz(1,ia) + dedphi*dxiazib
     &                             + d2edphi2*dphidzib*dphidxia
     &                             - dxia*ddrdz
                  hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxib
     &                             + d2edphi2*dphidxib*dphidyia
     &                             - dyia*ddrdx
                  hessy(2,ia) = hessy(2,ia) + dedphi*dyiayib
     &                             + d2edphi2*dphidyib*dphidyia
     &                             - dyia*ddrdy
                  hessz(2,ia) = hessz(2,ia) + dedphi*dyiazib
     &                             + d2edphi2*dphidzib*dphidyia
     &                             - dyia*ddrdz
                  hessx(3,ia) = hessx(3,ia) + dedphi*dziaxib
     &                             + d2edphi2*dphidxib*dphidzia
     &                             - dzia*ddrdx
                  hessy(3,ia) = hessy(3,ia) + dedphi*dziayib
     &                             + d2edphi2*dphidyib*dphidzia
     &                             - dzia*ddrdy
                  hessz(3,ia) = hessz(3,ia) + dedphi*dziazib
     &                             + d2edphi2*dphidzib*dphidzia
     &                             - dzia*ddrdz
                  hessx(1,ic) = hessx(1,ic) + dedphi*dxibxic
     &                             + d2edphi2*dphidxib*dphidxic
     &                             + (dxib-dxic)*ddrdx - d2drdxx
                  hessy(1,ic) = hessy(1,ic) + dedphi*dyibxic
     &                             + d2edphi2*dphidyib*dphidxic
     &                             + dyib*ddrdx - dxic*ddrdy - d2drdxy
                  hessz(1,ic) = hessz(1,ic) + dedphi*dzibxic
     &                             + d2edphi2*dphidzib*dphidxic
     &                             + dzib*ddrdx - dxic*ddrdz - d2drdxz
                  hessx(2,ic) = hessx(2,ic) + dedphi*dxibyic
     &                             + d2edphi2*dphidxib*dphidyic
     &                             + dxib*ddrdy - dyic*ddrdx - d2drdxy
                  hessy(2,ic) = hessy(2,ic) + dedphi*dyibyic
     &                             + d2edphi2*dphidyib*dphidyic
     &                             + (dyib-dyic)*ddrdy - d2drdyy
                  hessz(2,ic) = hessz(2,ic) + dedphi*dzibyic
     &                             + d2edphi2*dphidzib*dphidyic
     &                             + dzib*ddrdy - dyic*ddrdz - d2drdyz
                  hessx(3,ic) = hessx(3,ic) + dedphi*dxibzic
     &                             + d2edphi2*dphidxib*dphidzic
     &                             + dxib*ddrdz - dzic*ddrdx - d2drdxz
                  hessy(3,ic) = hessy(3,ic) + dedphi*dyibzic
     &                             + d2edphi2*dphidyib*dphidzic
     &                             + dyib*ddrdz - dzic*ddrdy - d2drdyz
                  hessz(3,ic) = hessz(3,ic) + dedphi*dzibzic
     &                             + d2edphi2*dphidzib*dphidzic
     &                             + (dzib-dzic)*ddrdz - d2drdzz
                  hessx(1,id) = hessx(1,id) + dedphi*dxibxid
     &                             + d2edphi2*dphidxib*dphidxid
     &                             - dxid*ddrdx
                  hessy(1,id) = hessy(1,id) + dedphi*dyibxid
     &                             + d2edphi2*dphidyib*dphidxid
     &                             - dxid*ddrdy
                  hessz(1,id) = hessz(1,id) + dedphi*dzibxid
     &                             + d2edphi2*dphidzib*dphidxid
     &                             - dxid*ddrdz
                  hessx(2,id) = hessx(2,id) + dedphi*dxibyid
     &                             + d2edphi2*dphidxib*dphidyid
     &                             - dyid*ddrdx
                  hessy(2,id) = hessy(2,id) + dedphi*dyibyid
     &                             + d2edphi2*dphidyib*dphidyid
     &                             - dyid*ddrdy
                  hessz(2,id) = hessz(2,id) + dedphi*dzibyid
     &                             + d2edphi2*dphidzib*dphidyid
     &                             - dyid*ddrdz
                  hessx(3,id) = hessx(3,id) + dedphi*dxibzid
     &                             + d2edphi2*dphidxib*dphidzid
     &                             - dzid*ddrdx
                  hessy(3,id) = hessy(3,id) + dedphi*dyibzid
     &                             + d2edphi2*dphidyib*dphidzid
     &                             - dzid*ddrdy
                  hessz(3,id) = hessz(3,id) + dedphi*dzibzid
     &                             + d2edphi2*dphidzib*dphidzid
     &                             - dzid*ddrdz
               else if (i .eq. ic) then
                  hessx(1,ic) = hessx(1,ic) + dedphi*dxicxic
     &                             + d2edphi2*dphidxic*dphidxic
     &                             + 2.0d0*dxic*ddrdx + d2drdxx
                  hessy(1,ic) = hessy(1,ic) + dedphi*dxicyic
     &                             + d2edphi2*dphidxic*dphidyic
     &                             + dyic*ddrdx + dxic*ddrdy + d2drdxy
                  hessz(1,ic) = hessz(1,ic) + dedphi*dxiczic
     &                             + d2edphi2*dphidxic*dphidzic
     &                             + dzic*ddrdx + dxic*ddrdz + d2drdxz
                  hessx(2,ic) = hessx(2,ic) + dedphi*dxicyic
     &                             + d2edphi2*dphidxic*dphidyic
     &                             + dxic*ddrdy + dyic*ddrdx + d2drdxy
                  hessy(2,ic) = hessy(2,ic) + dedphi*dyicyic
     &                             + d2edphi2*dphidyic*dphidyic
     &                             + 2.0d0*dyic*ddrdy + d2drdyy
                  hessz(2,ic) = hessz(2,ic) + dedphi*dyiczic
     &                             + d2edphi2*dphidyic*dphidzic
     &                             + dzic*ddrdy + dyic*ddrdz + d2drdyz
                  hessx(3,ic) = hessx(3,ic) + dedphi*dxiczic
     &                             + d2edphi2*dphidxic*dphidzic
     &                             + dxic*ddrdz + dzic*ddrdx + d2drdxz
                  hessy(3,ic) = hessy(3,ic) + dedphi*dyiczic
     &                             + d2edphi2*dphidyic*dphidzic
     &                             + dyic*ddrdz + dzic*ddrdy + d2drdyz
                  hessz(3,ic) = hessz(3,ic) + dedphi*dziczic
     &                             + d2edphi2*dphidzic*dphidzic
     &                             + 2.0d0*dzic*ddrdz + d2drdzz
                  hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxic
     &                             + d2edphi2*dphidxic*dphidxia
     &                             + dxia*ddrdx
                  hessy(1,ia) = hessy(1,ia) + dedphi*dxiayic
     &                             + d2edphi2*dphidyic*dphidxia
     &                             + dxia*ddrdy
                  hessz(1,ia) = hessz(1,ia) + dedphi*dxiazic
     &                             + d2edphi2*dphidzic*dphidxia
     &                             + dxia*ddrdz
                  hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxic
     &                             + d2edphi2*dphidxic*dphidyia
     &                             + dyia*ddrdx
                  hessy(2,ia) = hessy(2,ia) + dedphi*dyiayic
     &                             + d2edphi2*dphidyic*dphidyia
     &                             + dyia*ddrdy
                  hessz(2,ia) = hessz(2,ia) + dedphi*dyiazic
     &                             + d2edphi2*dphidzic*dphidyia
     &                             + dyia*ddrdz
                  hessx(3,ia) = hessx(3,ia) + dedphi*dziaxic
     &                             + d2edphi2*dphidxic*dphidzia
     &                             + dzia*ddrdx
                  hessy(3,ia) = hessy(3,ia) + dedphi*dziayic
     &                             + d2edphi2*dphidyic*dphidzia
     &                             + dzia*ddrdy
                  hessz(3,ia) = hessz(3,ia) + dedphi*dziazic
     &                             + d2edphi2*dphidzic*dphidzia
     &                             + dzia*ddrdz
                  hessx(1,ib) = hessx(1,ib) + dedphi*dxibxic
     &                             + d2edphi2*dphidxic*dphidxib
     &                             - (dxic-dxib)*ddrdx - d2drdxx
                  hessy(1,ib) = hessy(1,ib) + dedphi*dxibyic
     &                             + d2edphi2*dphidyic*dphidxib
     &                             - dyic*ddrdx + dxib*ddrdy - d2drdxy
                  hessz(1,ib) = hessz(1,ib) + dedphi*dxibzic
     &                             + d2edphi2*dphidzic*dphidxib
     &                             - dzic*ddrdx + dxib*ddrdz - d2drdxz
                  hessx(2,ib) = hessx(2,ib) + dedphi*dyibxic
     &                             + d2edphi2*dphidxic*dphidyib
     &                             - dxic*ddrdy + dyib*ddrdx - d2drdxy
                  hessy(2,ib) = hessy(2,ib) + dedphi*dyibyic
     &                             + d2edphi2*dphidyic*dphidyib
     &                             - (dyic-dyib)*ddrdy - d2drdyy
                  hessz(2,ib) = hessz(2,ib) + dedphi*dyibzic
     &                             + d2edphi2*dphidzic*dphidyib
     &                             - dzic*ddrdy + dyib*ddrdz - d2drdyz
                  hessx(3,ib) = hessx(3,ib) + dedphi*dzibxic
     &                             + d2edphi2*dphidxic*dphidzib
     &                             - dxic*ddrdz + dzib*ddrdx - d2drdxz
                  hessy(3,ib) = hessy(3,ib) + dedphi*dzibyic
     &                             + d2edphi2*dphidyic*dphidzib
     &                             - dyic*ddrdz + dzib*ddrdy - d2drdyz
                  hessz(3,ib) = hessz(3,ib) + dedphi*dzibzic
     &                             + d2edphi2*dphidzic*dphidzib
     &                             - (dzic-dzib)*ddrdz - d2drdzz
                  hessx(1,id) = hessx(1,id) + dedphi*dxicxid
     &                             + d2edphi2*dphidxic*dphidxid
     &                             + dxid*ddrdx
                  hessy(1,id) = hessy(1,id) + dedphi*dyicxid
     &                             + d2edphi2*dphidyic*dphidxid
     &                             + dxid*ddrdy
                  hessz(1,id) = hessz(1,id) + dedphi*dzicxid
     &                             + d2edphi2*dphidzic*dphidxid
     &                             + dxid*ddrdz
                  hessx(2,id) = hessx(2,id) + dedphi*dxicyid
     &                             + d2edphi2*dphidxic*dphidyid
     &                             + dyid*ddrdx
                  hessy(2,id) = hessy(2,id) + dedphi*dyicyid
     &                             + d2edphi2*dphidyic*dphidyid
     &                             + dyid*ddrdy
                  hessz(2,id) = hessz(2,id) + dedphi*dzicyid
     &                             + d2edphi2*dphidzic*dphidyid
     &                             + dyid*ddrdz
                  hessx(3,id) = hessx(3,id) + dedphi*dxiczid
     &                             + d2edphi2*dphidxic*dphidzid
     &                             + dzid*ddrdx
                  hessy(3,id) = hessy(3,id) + dedphi*dyiczid
     &                             + d2edphi2*dphidyic*dphidzid
     &                             + dzid*ddrdy
                  hessz(3,id) = hessz(3,id) + dedphi*dziczid
     &                             + d2edphi2*dphidzic*dphidzid
     &                             + dzid*ddrdz
               else if (i .eq. id) then
                  hessx(1,id) = hessx(1,id) + dedphi*dxidxid
     &                             + d2edphi2*dphidxid*dphidxid
                  hessy(1,id) = hessy(1,id) + dedphi*dxidyid
     &                             + d2edphi2*dphidxid*dphidyid
                  hessz(1,id) = hessz(1,id) + dedphi*dxidzid
     &                             + d2edphi2*dphidxid*dphidzid
                  hessx(2,id) = hessx(2,id) + dedphi*dxidyid
     &                             + d2edphi2*dphidxid*dphidyid
                  hessy(2,id) = hessy(2,id) + dedphi*dyidyid
     &                             + d2edphi2*dphidyid*dphidyid
                  hessz(2,id) = hessz(2,id) + dedphi*dyidzid
     &                             + d2edphi2*dphidyid*dphidzid
                  hessx(3,id) = hessx(3,id) + dedphi*dxidzid
     &                             + d2edphi2*dphidxid*dphidzid
                  hessy(3,id) = hessy(3,id) + dedphi*dyidzid
     &                             + d2edphi2*dphidyid*dphidzid
                  hessz(3,id) = hessz(3,id) + dedphi*dzidzid
     &                             + d2edphi2*dphidzid*dphidzid
                  hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxid
     &                             + d2edphi2*dphidxid*dphidxia
                  hessy(1,ia) = hessy(1,ia) + dedphi*dxiayid
     &                             + d2edphi2*dphidyid*dphidxia
                  hessz(1,ia) = hessz(1,ia) + dedphi*dxiazid
     &                             + d2edphi2*dphidzid*dphidxia
                  hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxid
     &                             + d2edphi2*dphidxid*dphidyia
                  hessy(2,ia) = hessy(2,ia) + dedphi*dyiayid
     &                             + d2edphi2*dphidyid*dphidyia
                  hessz(2,ia) = hessz(2,ia) + dedphi*dyiazid
     &                             + d2edphi2*dphidzid*dphidyia
                  hessx(3,ia) = hessx(3,ia) + dedphi*dziaxid
     &                             + d2edphi2*dphidxid*dphidzia
                  hessy(3,ia) = hessy(3,ia) + dedphi*dziayid
     &                             + d2edphi2*dphidyid*dphidzia
                  hessz(3,ia) = hessz(3,ia) + dedphi*dziazid
     &                             + d2edphi2*dphidzid*dphidzia
                  hessx(1,ib) = hessx(1,ib) + dedphi*dxibxid
     &                             + d2edphi2*dphidxid*dphidxib
     &                             - dxid*ddrdx
                  hessy(1,ib) = hessy(1,ib) + dedphi*dxibyid
     &                             + d2edphi2*dphidyid*dphidxib
     &                             - dyid*ddrdx
                  hessz(1,ib) = hessz(1,ib) + dedphi*dxibzid
     &                             + d2edphi2*dphidzid*dphidxib
     &                             - dzid*ddrdx
                  hessx(2,ib) = hessx(2,ib) + dedphi*dyibxid
     &                             + d2edphi2*dphidxid*dphidyib
     &                             - dxid*ddrdy
                  hessy(2,ib) = hessy(2,ib) + dedphi*dyibyid
     &                             + d2edphi2*dphidyid*dphidyib
     &                             - dyid*ddrdy
                  hessz(2,ib) = hessz(2,ib) + dedphi*dyibzid
     &                             + d2edphi2*dphidzid*dphidyib
     &                             - dzid*ddrdy
                  hessx(3,ib) = hessx(3,ib) + dedphi*dzibxid
     &                             + d2edphi2*dphidxid*dphidzib
     &                             - dxid*ddrdz
                  hessy(3,ib) = hessy(3,ib) + dedphi*dzibyid
     &                             + d2edphi2*dphidyid*dphidzib
     &                             - dyid*ddrdz
                  hessz(3,ib) = hessz(3,ib) + dedphi*dzibzid
     &                             + d2edphi2*dphidzid*dphidzib
     &                             - dzid*ddrdz
                  hessx(1,ic) = hessx(1,ic) + dedphi*dxicxid
     &                             + d2edphi2*dphidxid*dphidxic
     &                             + dxid*ddrdx
                  hessy(1,ic) = hessy(1,ic) + dedphi*dxicyid
     &                             + d2edphi2*dphidyid*dphidxic
     &                             + dyid*ddrdx
                  hessz(1,ic) = hessz(1,ic) + dedphi*dxiczid
     &                             + d2edphi2*dphidzid*dphidxic
     &                             + dzid*ddrdx
                  hessx(2,ic) = hessx(2,ic) + dedphi*dyicxid
     &                             + d2edphi2*dphidxid*dphidyic
     &                             + dxid*ddrdy
                  hessy(2,ic) = hessy(2,ic) + dedphi*dyicyid
     &                             + d2edphi2*dphidyid*dphidyic
     &                             + dyid*ddrdy
                  hessz(2,ic) = hessz(2,ic) + dedphi*dyiczid
     &                             + d2edphi2*dphidzid*dphidyic
     &                             + dzid*ddrdy
                  hessx(3,ic) = hessx(3,ic) + dedphi*dzicxid
     &                             + d2edphi2*dphidxid*dphidzic
     &                             + dxid*ddrdz
                  hessy(3,ic) = hessy(3,ic) + dedphi*dzicyid
     &                             + d2edphi2*dphidyid*dphidzic
     &                             + dyid*ddrdz
                  hessz(3,ic) = hessz(3,ic) + dedphi*dziczid
     &                             + d2edphi2*dphidzid*dphidzic
     &                             + dzid*ddrdz
               end if
            end if
         end if
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
c$Header: /data/petsun4/data1/src_solaris/mprtir/RCS/powsub.f,v 1.5 2010/11/29 06:35:32 avi Exp $
c$Log: powsub.f,v $
c Revision 1.5  2010/11/29  06:35:32  avi
c correct termination of powstr
c
c Revision 1.4  2008/01/02  01:24:01  avi
c gcc v4 compliant
c
c Revision 1.3  1999/12/08  22:59:51  avi
c right justify powstring
c
c Revision 1.2  1999/02/18  01:24:48  avi
c improve speed
c
c Revision 1.1  1999/02/15  10:07:28  avi
c Initial revision
      subroutine evalpow_test
      parameter(n=2)
      parameter(nterm=((n+3)*(n+2)*(n+1))/6)
      real*4 a(nterm),b(nterm),t(nterm,nterm)
      pointer (pa,a),(pb,b),(pt,t)
      character*256 string
      data x,y,z/1.,2.,3./

      pa=malloc(4*nterm)
      pb=malloc(4*nterm)
      pt=malloc(4*nterm**2)

      write(*,"('nterm=',i5)")nterm
      if(8*nterm.lt.255)then
        call powstring(n,8,string)
        write(*,"(a)")string(1:8*nterm)
      endif
      call evalpo1(x,y,z,n,a)
      call evalpow(x,y,z,n,b)
      do i=1,nterm
        write(*,"(i10,2f10.4)")i,a(i),b(i)
      enddo

      call free(pa)
      call free(pb)
      call free(pt)
      call exit(0)
      end

      subroutine evalpow(x,y,z,norder,p)
      real*4 p(0:((norder+3)*(norder+2)*(norder+1))/6-1)

      p(0)=1.
      if(norder.lt.1)return
      p(1)=z
      p(2)=y
      p(3)=x
      if(norder.lt.2)return
      p(4)=z*z
      p(5)=y*z
      p(6)=y*y
      p(7)=x*z
      p(8)=x*y
      p(9)=x*x
      if(norder.lt.3)return
      j=10
      do 7 i=3,norder
      fx=1.
      do 4 k=0,norder
      fy=fx
      do 5 l=0,norder
      fz=fy
      do 6 m=0,norder
      if(k+l+m.eq.i)then
        p(j)=fz
        j=j+1
      endif
    6 fz=fz*z
    5 fy=fy*y
    4 fx=fx*x
    7 continue

      return
      end

      subroutine poly3ev(x,y,z,norder,a,p,g)
      real*4 a(0:((norder+3)*(norder+2)*(norder+1))/6-1)
      real*4 p(0:((norder+3)*(norder+2)*(norder+1))/6-1)

      p(0)=1.
      g=a(0)
      if(norder.lt.1)return
      p(1)=z
      p(2)=y
      p(3)=x
      g=g+z*a(1)+y*a(2)+x*a(3)
      if(norder.lt.2)return
      p(4)=z*z
      p(5)=y*z
      p(6)=y*y
      p(7)=x*z
      p(8)=x*y
      p(9)=x*x
      g=g+a(4)*p(4)+a(5)*p(5)+a(6)*p(6)+a(7)*p(7)+a(8)*p(8)+a(9)*p(9)
      if(norder.lt.3)return
      j=10
      do 7 i=3,norder
      fx=1.
      do 4 k=0,norder
      fy=fx
      do 5 l=0,norder
      fz=fy
      do 6 m=0,norder
      if(k+l+m.eq.i)then
        p(j)=fz
        g=g+a(j)*fz
        j=j+1
      endif
    6 fz=fz*z
    5 fy=fy*y
    4 fx=fx*x
    7 continue

      return
      end

      subroutine evalpo1(x,y,z,norder,a)
      real*4 a(((norder+3)*(norder+2)*(norder+1))/6)

      j=0
      do 7 i=0,norder
      do 7 k=0,norder
      do 7 l=0,norder
      do 7 m=0,norder
      if(k+l+m.ne.i)goto 7
      j=j+1
      a(j)=x**k*y**l*z**m
    7 continue

      return
      end

      subroutine powstring(norder,lenfield,string)
      character*255 string
      logical*4 ldebug/.false./

      nterm=((norder+3)*(norder+2)*(norder+1))/6
      string=' '
      j=lenfield-5
      do 7 i=0,norder
      do 7 k=0,norder
      do 7 l=0,norder
      do 7 m=0,norder
      if(k+l+m.ne.i)goto 7
      write(string(j:j+lenfield-1),"('x',i1,'y',i1,'z',i1)")k,l,m
      j=j+lenfield
    7 continue
      if(ldebug)write(*,"(a)")string(1:lenfield*nterm)

      string(lenfield*nterm+1:)=char(0)
      return
      end

ext: f

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
      double precision function pdg2pdf(ih,ipdg,x,xmu)
c***************************************************************************
c     Based on pdf.f, wrapper for calling the pdf of MCFM
c***************************************************************************
      implicit none
c
c     Arguments
c
      DOUBLE  PRECISION x,xmu
      INTEGER IH,ipdg
C
C     Include
C
      include 'pdf.inc'
C
      double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
      integer mode,Irt,i,j
      double precision xlast(2),xmulast(2),pdflast(-7:7,2),q2max
      character*7 pdlabellast(2)
      double precision epa_electron,epa_proton
      integer ipart,ireuse,iporg,ihlast(2), iisr
      save xlast,xmulast,pdflast,pdlabellast,ihlast
      data xlast/2*-99d9/
      data xmulast/2*-99d9/
      data pdflast/30*-99d9/
      data pdlabellast/2*'abcdefg'/
      data ihlast/2*-99/

c     Make sure we have a reasonable Bjorken x. Note that even though
c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
c     instead of stopping the code, as this might accidentally happen.
      if (x.eq.0d0) then
         pdg2pdf=0d0
         return
      elseif (x.lt.0d0 .or. x.gt.1d0) then
         write (*,*) 'PDF not supported for Bjorken x ', x
         open(unit=26,file='../../../error',status='unknown')
         write(26,*) 'Error: PDF not supported for Bjorken x ',x
         stop 1
      endif

      ipart=ipdg
      iisr=0
      if(iabs(ipart).eq.21) ipart=0
      if(iabs(ipart).eq.22) ipart=7
      if(iabs(ipart).eq.11) iisr=1
      if(iabs(ipart).eq.11) ipart=7
      iporg=ipart

c     This will be called for any PDG code, but we only support up to 7
      if(iabs(ipart).gt.7)then
         write(*,*) 'PDF not supported for pdg ',ipdg
         write(*,*) 'For lepton colliders, please set the lpp* '//
     $    'variables to 0 in the run_card'
         open(unit=26,file='../../../error',status='unknown')
         write(26,*) 'Error: PDF not supported for pdg ',ipdg
         stop 1
      endif

      ireuse = 0
      do i=1,2
c     Check if result can be reused since any of last two calls
         if (x.eq.xlast(i) .and. xmu.eq.xmulast(i) .and.
     $        pdlabel.eq.pdlabellast(i) .and. ih.eq.ihlast(i)) then
            ireuse = i
         endif
      enddo

c     Reuse previous result, if possible
      if (ireuse.gt.0.)then
         if (pdflast(iporg,ireuse).ne.-99d9) then
            pdg2pdf=pdflast(iporg,ireuse)
            return
         endif
      endif

c     Bjorken x and/or facrorization scale and/or PDF set are not
c     identical to the saved values: this means a new event and we
c     should reset everything to compute new PDF values. Also, determine
c     if we should fill ireuse=1 or ireuse=2.
      if (ireuse.eq.0.and.xlast(1).ne.-99d9.and.xlast(2).ne.-99d9)then
         do i=1,2
            xlast(i)=-99d9
            xmulast(i)=-99d9
            do j=-7,7
               pdflast(j,i)=-99d9
            enddo
            pdlabellast(i)='abcdefg'
            ihlast(i)=-99
         enddo
c     everything has been reset. Now set ireuse=1 to fill the first
c     arrays of saved values below
         ireuse=1
      else if(ireuse.eq.0.and.xlast(1).ne.-99d9)then
c     This is first call after everything has been reset, so the first
c     arrays are already filled with the saved values (hence
c     xlast(1).ne.-99d9). Fill the second arrays of saved values (done
c     below) by setting ireuse=2
         ireuse=2
      else if(ireuse.eq.0)then
c     Special: only used for the very first call to this function:
c     xlast(i) are initialized as data statements to be equal to -99d9
         ireuse=1
      endif

c     Give the current values to the arrays that should be
c     saved. 'pdflast' is filled below.
      xlast(ireuse)=x
      xmulast(ireuse)=xmu
      pdlabellast(ireuse)=pdlabel
      ihlast(ireuse)=ih

      if(iabs(ipart).eq.7.and.ih.gt.1) then
         q2max=xmu*xmu
         if(ih.eq.3 .and. iisr.ne.1 ) then       !from the electron
            pdg2pdf=epa_electron(x,q2max)
         elseif(ih.eq.3 .and. iisr.eq.1 ) then       !from the electron
            pdg2pdf=1.0d0
         elseif(ih .eq. 2) then !from a proton without breaking
            pdg2pdf=epa_proton(x,q2max)
         endif
         pdflast(iporg,ireuse)=pdg2pdf
         return
      endif

      if (pdlabel(1:5) .eq. 'cteq3') then
C
         if (pdlabel .eq. 'cteq3_m') then
            mode=1
         elseif (pdlabel .eq. 'cteq3_l') then
            mode=2
         elseif (pdlabel .eq. 'cteq3_d') then
            mode=3
         endif

         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
     $      ipart=sign(3-iabs(ipart),ipart)

         pdg2pdf=Ctq3df(mode,ipart,x,xmu,Irt)/x

         if(ipdg.ge.1.and.ipdg.le.2)
     $      pdg2pdf=pdg2pdf+Ctq3df(mode,-ipart,x,xmu,Irt)/x

C
      elseif (pdlabel(1:5) .eq. 'cteq4') then
C
         if (pdlabel .eq. 'cteq4_m') then
            mode=1
         elseif (pdlabel .eq. 'cteq4_d') then
            mode=2
         elseif (pdlabel .eq. 'cteq4_l') then
            mode=3
         elseif (pdlabel .eq. 'cteq4a1') then
            mode=4
         elseif (pdlabel .eq. 'cteq4a2') then
            mode=5
         elseif (pdlabel .eq. 'cteq4a3') then
            mode=6
         elseif (pdlabel .eq. 'cteq4a4') then
            mode=7
         elseif (pdlabel .eq. 'cteq4a5') then
            mode=8
         elseif (pdlabel .eq. 'cteq4hj') then
            mode=9
         elseif (pdlabel .eq. 'cteq4lq') then
            mode=10
         endif

         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
     $      ipart=sign(3-iabs(ipart),ipart)

         pdg2pdf=Ctq4Fn(mode,ipart,x,xmu)
C
      elseif (pdlabel .eq. 'cteq5l1') then
C
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
     $      ipart=sign(3-iabs(ipart),ipart)

         pdg2pdf=Ctq5L(ipart,x,xmu)
C
      elseif ((pdlabel(1:5) .eq. 'cteq5') .or.
     .        (pdlabel(1:4) .eq. 'ctq5')) then
C
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
     $      ipart=sign(3-iabs(ipart),ipart)

         pdg2pdf=Ctq5Pdf(ipart,x,xmu)
C
      elseif (pdlabel(1:5) .eq. 'cteq6') then
C
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
     $      ipart=sign(3-iabs(ipart),ipart)

         pdg2pdf=Ctq6Pdf(ipart,x,xmu)
      else
         call pftopdg(ih,x,xmu,pdflast(-7,ireuse))
         pdg2pdf=pdflast(iporg,ireuse)
      endif

      pdflast(iporg,ireuse)=pdg2pdf
      return
      end

max_stars_repo_name: oliviermattelaer/MGISR

max_issues_repo_name: oliviermattelaer/MGISR

max_forks_repo_name: oliviermattelaer/MGISR


text:
less_toxic(bb1,cc1).
less_toxic(ii1,f1).
less_toxic(x1,p1).
less_toxic(ee1,d1).
less_toxic(hh1,i1).
less_toxic(t1,f1).
less_toxic(k1,h1).
less_toxic(ii1,t1).
less_toxic(hh1,cc1).
less_toxic(ee1,l1).
less_toxic(k1,c1).
less_toxic(x1,w1).
less_toxic(o1,w1).
less_toxic(s1,e1).
less_toxic(s1,t1).
less_toxic(kk1,jj1).
less_toxic(x1,dd1).
less_toxic(bb1,w1).
less_toxic(x1,d1).
less_toxic(j1,t1).
less_toxic(n1,jj1).
less_toxic(hh1,a1).
less_toxic(w1,t1).
less_toxic(cc1,u1).
less_toxic(ll1,jj1).
less_toxic(bb1,c1).
less_toxic(e1,f1).
less_toxic(ee1,p1).
less_toxic(ll1,cc1).
less_toxic(ll1,l1).
less_toxic(s1,z1).
less_toxic(y1,f1).
less_toxic(ii1,u1).
less_toxic(ii1,jj1).
less_toxic(n1,l1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!decision control using if
program name
    implicit none
    INTEGER :: a, b
    ! Syntax:
    ! if(condition) then
    ! statement
    ! statement
    ! statement
    !.
    !.
    !.
    !else
    !statement
    !statement
    !statement
    !.
    !.
    !.
    !endif

    if (.true.) then
        print *, "It is true statement"
    end if

    a = 5
    b = 4

    if (a .gt. b) then
        print *, "a is greater than b"
    else
        print *, "a is less than b"
    end if

end program name

max_stars_repo_name: SusheelThapa/Code-With-FORTRAN

max_issues_repo_name: SusheelThapa/Code-With-FORTRAN

max_forks_repo_name: SusheelThapa/Code-With-FORTRAN


text:
program example3
  ! array dependencies. All of the first loop must complete before
  ! starting the second loop
  implicit none
  integer :: i
  real :: a(4)
  a = 0.0
  do i=1,4
     !$omp task out(a)
     call my_work1(i,a)
     !$omp end task
  end do

  do i=1,4
     !$omp task in(a)
     call my_work2(i,a)
     !$omp end task
  end do
end program example3

subroutine my_work1(i,a)
  implicit none
  real :: a(4)
  integer :: i
  print *, "hello", i
  a(i)=i
end subroutine my_work1

subroutine my_work2(i,a)
  implicit none
  real :: a(i)
  integer :: i
  print *, "world", i,a(i)
end subroutine my_work2

max_stars_repo_name: stfc/PSycloneBench

max_issues_repo_name: stfc/PSycloneBench

max_forks_repo_name: stfc/PSycloneBench


text:
      SUBROUTINE CKINPT(NSUBX)
C
C        'CKINPT' CALCULATES THE TRANSFORMATION MATRIX BETWEEN THE
C        ELEMENT FRAME AND THE BODY FRAME
C
        IMPLICIT REAL*8 (A-H,O-Z)
      REAL * 8 NSUBX(3)
      INTEGER*4 ACNTRL
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/DEBUG2/ IOUT,JOUT,KLUGE
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/MOMENT/ ACNTRL,IVISCS,IATTDE,IMGMTS,IWHEEL,NPULSE
C
      COMMON/RPOOL5/ CKMAT(3,3,10),FM2(3,3)
C
      COMMON/XIN3  / ALFAEK(10),BETAEK(10),GAMAEK(10),ETTAD,ZETTAD
C
      DIMENSION S(3,3),CK(3,3,10)
C
C
C
      DO 10 I=1,3
      DO 10 J=1,3
      FM2(I,J)=0.D0
   10 FM2(I,I)=1.D0
C
C
      IF(IVISCS.EQ.0) GO TO 15
C
C                    CALCULATE NSUBX ARRAY
      CETA=DCOS(ETTAD*RADIAN)
      SETA=DSIN(ETTAD*RADIAN)
      SZETA=DSIN(ZETTAD*RADIAN)
      CZETA=DCOS(ZETTAD*RADIAN)
C
      NSUBX(1)=SETA*CZETA
      NSUBX(2)=SETA*SZETA
      NSUBX(3)=CETA
C
C
   15 DO 30 K=1,IK
C
C
      SALFAE=DSIN(ALFAEK(K)*RADIAN)
      CALFAE=DCOS(ALFAEK(K)*RADIAN)
      SBETAE=DSIN(BETAEK(K)*RADIAN)
      CBETAE=DCOS(BETAEK(K)*RADIAN)
      SGAMAE=DSIN(GAMAEK(K)*RADIAN)
      CGAMAE=DCOS(GAMAEK(K)*RADIAN)
C
      S(1,1)=CGAMAE*CBETAE + SALFAE*SGAMAE*SBETAE
      S(1,2)=CGAMAE*SALFAE*SBETAE - SGAMAE*CBETAE
      S(1,3)=CALFAE*SBETAE
C
      S(2,1)=SGAMAE*CALFAE
      S(2,2)=CALFAE*CGAMAE
      S(2,3)=-SALFAE
C
      S(3,1)=SGAMAE*SALFAE*CBETAE - SBETAE*CGAMAE
      S(3,2)=CGAMAE*SALFAE*CBETAE + SBETAE*SGAMAE
      S(3,3)=CALFAE*CBETAE
C
C
      DO 20 I=1,3
      DO 20 J=1,3
   20 CK(I,J,K)=S(I,J)
C
   30 CONTINUE
C
      N=1
      DO 70 K=1,IK
      L=IK-K1+N
      IF((K-K1).GT.0) L=K-K1
C
      DO 40 I=1,3
      DO 40 J=1,3
   40 CKMAT(I,J,K)=CK(I,J,L)
      N=N + 1
C
C
   70 CONTINUE
      IF(IOUT.EQ.1) RETURN
      WRITE(6,10000)
      DO 80 K=1,IK
   80 WRITE(6,10001) ((CKMAT(I,J,K),J=1,3),I=1,3)
C
      RETURN
10000 FORMAT('0',3X,'CKMAT')
10001 FORMAT('0',3X,3(G15.8,2X)/3X,3(G15.8,2X)/3X,3(G15.8,2X))
      END

ext: for


text:
C
C	$Id: gray.f,v 1.4 2008-07-27 00:17:14 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GRAY
C
C SUBROUTINE GRAY COLORS HALF-TONE CELL (I,J) WITH INTENSITY INTEN.
C THE ROUTINE ASSUMES 8X8 CELL SIZE ON A VIRTUAL SCREEN 1024X1024.
C
      DIMENSION       IFOT(16)   ,JFOT(16)
      DIMENSION       WIND2(4)   ,VIEW2(4)
CCC      DIMENSION       MX(16)  ,MY(16)
      COMMON /HAFTO1/ I          ,J          ,INTEN
      COMMON /HAFTO4/ NPTMAX     ,NPOINT     ,XPNT(50)  ,YPNT(50)
      SAVE
C
      DATA
     1 IFOT(1),IFOT(2),IFOT(3),IFOT(4),IFOT(5),IFOT(6),IFOT(7),IFOT(8)/
     2   1,      5,      1,      5,      3,      7,      3,      7   /
      DATA
     1 IFOT(9),IFOT(10),IFOT(11),IFOT(12),IFOT(13),IFOT(14),IFOT(15)/
     2   3,      7,       3,       7,       1,       5,       1/,
     3 IFOT(16)/
     4   5    /
C
      DATA
     1 JFOT(1),JFOT(2),JFOT(3),JFOT(4),JFOT(5),JFOT(6),JFOT(7),JFOT(8)/
     2   1,      5,      5,      1,      3,     7,      7,      3   /
      DATA
     1 JFOT(9),JFOT(10),JFOT(11),JFOT(12),JFOT(13),JFOT(14),JFOT(15)/
     2   1,      5,       5,       1,       3,       7,       7/,
     3 JFOT(16)/
     4   3    /
      IF (INTEN) 103,103,101
  101 I1 = I*8
      J1 = J*8
      IF ((NPOINT+INTEN) .LE.NPTMAX)  GO TO 1015
      CALL GETSET (VIEW2(1),VIEW2(2),VIEW2(3),VIEW2(4),
     +             WIND2(1),WIND2(2),WIND2(3),WIND2(4),IOLLS2)
      CALL SET (0.,1.,0.,1.,0.,1023.,0.,1023.,1)
      CALL POINTS (XPNT,YPNT,NPOINT,0,0)
      CALL SET (VIEW2(1),VIEW2(2),VIEW2(3),VIEW2(4),
     +          WIND2(1),WIND2(2),WIND2(3),WIND2(4),IOLLS2)
      NPOINT = 0
 1015 DO 102 I2=1,INTEN
      NPOINT = NPOINT + 1
      XPNT(NPOINT) = I1+IFOT(I2)
      YPNT(NPOINT) = J1+JFOT(I2)
  102 CONTINUE
  103 RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


    Program zcgesv_example

!     ZCGESV Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_interfaces, Only: zcgesv
      Use lapack_precision, Only: sp, dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, info, iter, j, lda, ldb, ldx, n, r
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:, :), work(:), x(:, :)
      Complex (Kind=sp), Allocatable :: swork(:)
      Real (Kind=dp), Allocatable :: rwork(:)
      Integer, Allocatable :: ipiv(:)
!     .. Executable Statements ..
      Write (nout, *) 'ZCGESV Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n, r
      lda = n
      ldb = n
      ldx = n
      Allocate (a(lda,n), b(ldb,r), work(n*r), x(ldx,r), swork(n*(n+ &
        r)), ipiv(n), rwork(n))

!     Read A and B from data file

      Read (nin, *)(a(i,1:n), i=1, n)
      Read (nin, *)(b(i,1:r), i=1, n)

!     Solve the equations Ax = b for x

      Call zcgesv(n, r, a, lda, ipiv, b, ldb, x, ldx, work, swork, rwork, &
        iter, info)

      If (info==0) Then

!       Print solution

        Write (nout, *) 'Solution'
        Write (nout, 100)((x(i,r),j=1,r), i=1, n)

!       Print pivot indices

        Write (nout, *)
        Write (nout, *) 'Pivot indices'
        Write (nout, 110) ipiv(1:n)

      Else
        Write (nout, 120) 'The (', info, ',', info, ')', &
          ' element of the factor U is zero'
      End If

100   Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
110   Format (1X, 7I11)
120   Format (1X, A, I3, A, I3, A, A)
    End Program


text:
Address(Emerald Bay Drive) is a residential street in South Davis.

 Intersecting Streets

  Danbury Street
  Chateau Lane
  Ranier Street
  Chateau Lane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! ers_euler_ideal_1d.f90
!
! Exact Riemann Solver for ideal 1D euler equations, producing a
! converged exact solution.
!
! https://github.com/JOThurgood/SimpleCFD/wiki/ers_euler_1d.f90

module shared_data ! a common block essentially

  implicit none

  integer, parameter :: num=selected_real_kind(p=15) !p=6 or 15
  integer :: niter !debug

  real(num) :: rhol, ul, pl, rhor, ur, pr !W_l and W_r
  real(num) :: gamma, gm, gp, g1, g2, g3, g4, g5, g6, g7 !gamma + const
  real(num) :: al,ar, bl, br !more constants
  real(num) :: cl, cr !sound speeds
  real(num) :: ps, us, rhosl, rhosr !star region vars

  integer :: nx !sampling stuff
  real(num) :: x_min, x_max
  real(num) :: t, x0
  real(num), allocatable, dimension(:) :: x, rho_x, u_x, p_x, en_x

  public

  contains

  subroutine constants
    !gamma
    gm = gamma - 1.0_num
    gp = gamma + 1.0_num
    g1 = gm / 2.0_num / gamma
    g2 = gp / 2.0_num / gamma
    g3 = 2.0_num * gamma / gm
    g4 = 2.0_num / gm
    g5 = 2.0_num / gp
    g6 = gm / gp
    g7 = gm / 2.0_num
    !data dependant constants for finding the root of the "f" function
    al = g5 / rhol !only used in pressure function - move out of shared?
    ar = g5 / rhor
    bl = g6 * pl
    br = g6 * pr
    !sound speeds
    cl = sqrt(gamma * pl /rhol)
    cr = sqrt(gamma * pr / rhor)
  end subroutine constants

end module shared_data

module user

  use shared_data

  implicit none

  public

  contains

  subroutine initial_conditions

    ! set these for your problem - default is Sod
    rhol = 1.0_num
    ul = 0.0_num
    pl = 1.0_num
    rhor = 0.125_num
    ur = 0.0_num
    pr = 0.10_num
    gamma = 1.4_num

  end subroutine initial_conditions

  subroutine control !sampling and output options
    t = 0.25_num
    x_min = 0.0_num
    x_max = 1.0_num
    x0 = 0.5_num
    nx = 1000
  end subroutine

end module user

module riemann !subroutines related to calculating star states

  use shared_data

  implicit none

  private

  public :: check_positivity, pstar, ustar, rhostar, sample, &
    newton_raphson

  !if when all is said and done there are any vars only used at
  ! this stage it would be sensible to move here from shared_data

  contains

  subroutine check_positivity
    real(num) :: du, du_crit
    du = ul - ur
    du_crit = cl * g4 + cr * g4
    if (du_crit < du) then
      print *, 'warning - pressure positivity condition violated'
      print *, 'your ICs will create vacuum'
      print *, 'STOP'
      STOP
    else
      print *, 'initial conditions pass p. positivity test OK'
    endif
  end subroutine check_positivity

  subroutine pstar
    call guess_ps
    call newton_raphson
  end subroutine pstar

  subroutine ustar
    us = 0.5_num * (ul + ur) + &
      & 0.5_num * (f(ps,rightonly=1) - f(ps,leftonly=1))
  end subroutine ustar

  subroutine rhostar
    integer :: i
    real(num) :: ck, pk, prat, rhok
    real(num) :: rhos !temp => rhosl or rhosr

    rhosl = 1e15_num !make errors obvious
    rhosr = 1e15_num
    do i = 0,1
      rhos = 1e15_num
      if (i == 0) then
        ck = cl
        pk =  pl
        rhok = rhol
        prat = ps / pk
      else
        ck = cr
        pk = pr
        rhok = rhor
        prat = ps / pk
      endif
      if (ps > pk) then
        rhos = rhok * ( (prat + g6) / (g6 * prat + 1.0_num) )
      else
        rhos = rhok * prat**(1.0_num/gamma)
      endif
      if (i == 0) then
        rhosl = rhos
      else
        rhosr = rhos
      endif
    enddo

  end subroutine rhostar

  subroutine guess_ps !initial guess on pstar
    ! in the interests of speed you, can base your guess on the value of
    ! pstar as the  given by an appropriate choice of approximate
    ! R solver. This should mean that the Newton-Raphson iterates in as
    ! few steps as possible.

    ps = 0.5_num * (pl + pr) ! For now just always use arithmetic mean
  end subroutine guess_ps

  subroutine newton_raphson

    integer :: i
    real(num) :: pold, delta
    real(num) :: tol = 1.0e-6_num
    real(num) :: rpc

    i = 0 !need to define these instead of header
    rpc = 1e15_num ! for sake of test module

    do
      if (rpc <= tol) exit !condition on tolerance
      pold = ps
      delta = f(pold)/fprime(pold)
      ps = pold - delta !f(pold)/fprime(pold)
      rpc = 2.0_num * abs((ps - pold) / (ps + pold))
      i = i + 1
      if (ps < 0.0_num) ps = tol !to correct -ve p guesses
      if (i > 20) print *,'nonconvergence of pstar iteration'
      if (i > 20) exit
    enddo

    niter = i !for test module
  end subroutine newton_raphson

  real(num) function f(p,leftonly,rightonly) !root function
    real(num), intent(in) :: p
    integer, intent(in), optional :: leftonly, rightonly
    real(num) :: ak,bk, pk, prat, ck
    integer :: i, i0, i1

    f = 0.0_num
    i0 = 0
    i1 = 1

    if (present(leftonly)) i1 = 0
    if (present(rightonly)) i0 = 1
    if (present(leftonly) .and. present(rightonly)) then
      PRINT *, 'what are you doing!?'
      STOP
    endif

    do i = i0,i1

      if (i == 0) then
        ak = al
        bk = bl
        ck = cl
        pk =  pl
        prat = p / pk
      else
        ak = ar
        bk = br
        ck = cr
        pk = pr
        prat = p / pk
      endif

      if (p > pk) then
        f = f + (p - pk) * sqrt(ak / (p + bk))
      else
        f = f + g4 * ck * (prat**g1 - 1.0_num)
      endif

    enddo

    f = f + (ur - ul)

    return
  end function

  real(num) function fprime(p) ! 1st derivative
    real(num), intent(in) :: p
    real(num) :: ak,bk,ck, pk, rhok, prat
    integer :: i

    fprime = 0.0_num

    do i = 0,1

      if (i == 0) then
        ak = al
        bk = bl
        ck = cl
        pk = pl
        rhok = rhol
        prat = p / pk
      else
        ak = ar
        bk = br
        ck = cr
        pk = pr
        rhok = rhor
        prat = p / pk
      endif

      if (p > pk) then
        fprime = fprime + &
          & sqrt(ak / (bk + p)) * (1.0_num - 0.5_num * (p-pk)/(bk+p))
      else
        fprime = fprime + prat**(-g2) / (rhok * ck)
      endif

    enddo

    return
  end function

  subroutine sample

    real(num) :: S
    real(num) :: Sk,Shk, Stk
    real(num) :: Sl,Shl, Stl
    real(num) :: Sr,Shr, Str
    real(num) :: ck, pk, rhok, rhosk, uk
    integer :: ix
    real(num) :: dx

    allocate(x(1:nx))
    allocate(rho_x(1:nx))
    allocate(p_x(1:nx))
    allocate(u_x(1:nx))
    allocate(en_x(1:nx))

    dx = (x_max - x_min) / REAL(nx,num)
    do ix = 1, nx
     x(ix) = x_min + REAL(ix,num) * dx
    enddo

    !sampling constants
    Sl = ul - cl * sqrt(g2*ps/pl + g1)
    Sr = ur + cr * sqrt(g2*ps/pr + g1)
    Shl = ul - cl
    Shr = ur + cr
    Stl = us - cl * (ps/pl)**(g1)
    Str = us + cr * (ps/pr)**(g1)

    do ix = 1, nx
      S = (x(ix)-x0) / t

      if (S < us) then !left of contact
!        print *, x(ix),'left of contact'
        Sk = Sl !set stuff like S_k to s_l
        Shk = Shl
        Stk = Stl
        pk = pl
        rhok = rhol
        rhosk = rhosl
        ck = cl
        uk = ul
        if (ps > pk) then !left shock
          if (S < Sk) then !W=Wk
!            print *, x(ix), 'region undisturbed (shock moving into it)'
            rho_x(ix) = rhok
            p_x(ix) = pk
            u_x(ix) = uk
          else !W=W*k
!            print *, x(ix), 'region shocked'
            rho_x(ix) = rhosk
            p_x(ix) = ps
            u_x(ix) = us
          endif
        else  !left rarefaction
          if (S < Shk) then !W=Wk
 !           print *, x(ix), 'region undisturbed (rarefaction moving into it)'
            rho_x(ix) = rhok
            p_x(ix) = pk
            u_x(ix) = uk
          else
            if (S>Stk) then !W = W*k
!            print *, x(ix), 'region between rarefaction tail and CD'
              rho_x(ix) = rhosk
              p_x(ix) = ps
              u_x(ix) = us
            else  !W = Wfan
!            print *, x(ix), 'rarefaction fan'
              rho_x(ix) = rhok * ( g5 + g6 / ck * (uk-S))**g4
              p_x(ix) = pk *  ( g5  + g6 / ck * (uk-S))**g3
              u_x(ix) = g5 * (ck + g7*uk + S)
            endif
          endif
        endif

      else !right of contact
!        print *,x(ix), 'right of contact'
        Sk = Sr
        Shk = Shr
        Stk = Str
        pk = pr
        rhok = rhor
        rhosk = rhosr
        ck = cr
        uk = ur

        if (ps > pk) then !right shock
          if (S > Sk) then !W=Wk
 !           print *, x(ix), 'region undisturbed (shock moving into it)'
            rho_x(ix) = rhok
            p_x(ix) = pk
            u_x(ix) = uk
          else !W=W*k
 !           print *, x(ix), 'region shocked'
            rho_x(ix) = rhosk
            p_x(ix) = ps
            u_x(ix) = us
          endif
        else  !right rarefaction
          if (S > Shk) then !W=Wk
 !           print *, x(ix), 'region undisturbed (rarefaction moving into it)'
            rho_x(ix) = rhok
            p_x(ix) = pk
            u_x(ix) = uk
          else
            if (S<Stk) then !W = W*k
  !          print *, x(ix), 'region between rarefaction tail and CD'
              rho_x(ix) = rhosk
              p_x(ix) = ps
              u_x(ix) = us
            else  !W = Wfan
  !          print *, x(ix), 'rarefaction fan'
              rho_x(ix) = rhok * ( g5 - g6 / ck * (uk-S))**g4
              p_x(ix) = pk *  ( g5 - g6 / ck * (uk-S))**g3
              u_x(ix) = g5 * (-ck + g7 * uk + S)
            endif
          endif
        endif

      endif !left or right
    end do !ix

    en_x = p_x / gm / rho_x

  end subroutine sample

end module riemann

module tests !subroutines for automatic testing

  use shared_data
  use riemann

  implicit none

  private

  public :: test_starvals, test_1,test_2,test_3,test_4,test_5

  logical :: test_star = .true.
  logical :: verbose = .false.

  contains

  subroutine test_starvals
    call test_pstar
    call test_ustar
    call test_rhostar
    if (test_star) then
      print *, 'passed all comparison tests for star vals'
    else
      print *, 'failed a star test - run with verbose to debug'
    endif
  end subroutine test_starvals
  subroutine test_pstar !check pstar and iteration info
    logical :: test1 = .true.
    real(num) :: tol = 1e-6_num
    real(num), dimension(5) :: p_toro
    real(num), dimension(4,5) :: p0_toro
    integer, dimension(4,5) :: i_toro
    integer :: i

    !toro
    p_toro = (/0.30313_num, 0.00189_num, 460.894_num, 46.0950_num, &
      & 1691.64_num /)

    i_toro = reshape( (/3,5,3,5, &
                      1, 8, 8 ,9, &
                      5, 4, 3, 4, &
                      5, 4, 3, 4, &
                      4, 5, 4, 6/), shape(i_toro))!, order = (/2,1/) )

    p0_toro = reshape((/0.30677_num, 0.55_num, 0.31527_num, 0.55_num, &
      & 0.00189387344_num, tol, tol, 0.4_num, &
      & 912.449_num, 500.005_num, 464.108_num, 500.005_num, &
      & 82.9831_num, 50.005_num, 46.4162_num, 50.005_num, &
      & 2322.65_num, 781.353_num, 1241.21_num, 253.494_num/), &
      & shape(p0_toro) )

!    do i=1,5
!      print *,'i_toro',i_toro(:,i)
!      print *,'p0_toro',p0_toro(:,i)
!    enddo

    call test_1 ! initialise  Sods' test
    do i = 1, 4
      ps = p0_toro(i,1)
      call newton_raphson
      if (verbose) &
        & print *, 'test 1: ','p0', p0_toro(i,1), 'ps',ps,'niter',niter
      if (niter .ne. i_toro(i,1)) then
        print *,'fail'
        test1 = .false.
      endif
    enddo

    call test_2
    do i = 1, 4
      ps = p0_toro(i,2)
      call newton_raphson
      if (verbose) &
        & print *, 'test 2: ','p0', p0_toro(i,2), 'ps',ps,'niter',niter
      if (niter .ne. i_toro(i,2)) print *,'fail'
    enddo

    call test_3
    do i = 1, 4
      ps = p0_toro(i,3)
      call newton_raphson
      if (verbose) &
        & print *, 'test 3: ','p0', p0_toro(i,3), 'ps',ps,'niter',niter
      if (niter .ne. i_toro(i,3)) print *,'fail'
    enddo

    call test_4
    do i = 1, 4
      ps = p0_toro(i,4)
      call newton_raphson
      if (verbose) &
        & print *, 'test 4: ','p0', p0_toro(i,4), 'ps',ps,'niter',niter
      if (niter .ne. i_toro(i,4)) print *,'fail'
    enddo

    call test_5
    do i = 1, 4
      ps = p0_toro(i,5)
      call newton_raphson
      if (verbose) &
        & print *, 'test 5: ','p0', p0_toro(i,5), 'ps',ps,'niter',niter
      if (niter .ne. i_toro(i,5)) print *,'fail'
    enddo

    if (test1) then
      print *,'passed pstar test'
    else
      print *,'failed pstar test'
      test_star = .false.
    endif

  end subroutine test_pstar

  subroutine test_ustar
    real (num) :: diff, utoro
    logical :: test2 = .true.

    call test_1
    call pstar
    call ustar
    utoro = 0.92745_num
    diff = 2.0_num * (us - utoro) / (us + utoro)
    if (verbose) print *, 'ustar-test1 and diff',us,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'ustar fail on test1',diff
      test2 = .false.
    endif

    call test_2
    call pstar
    call ustar
    utoro =  0.00000_num
    diff = 2.0_num * (us - utoro) / (us + utoro)
    if (verbose) print *, 'ustar-test2 and diff',us,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *, 'ustar fail on test2',diff
      test2 = .false.
    endif

    call test_3
    call pstar
    call ustar
    utoro = 19.5975_num
    diff = 2.0_num * (us - utoro) / (us + utoro)
    if (verbose) print *, 'ustar-test3 and diff',us,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *, 'ustar fail on test3',diff
      test2 = .false.
    endif

    call test_4
    call pstar
    call ustar
    utoro = -6.19633_num
    diff = 2.0_num * (us - utoro) / (us + utoro)
    if (verbose) print *, 'ustar-test4 and diff',us,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *, 'ustar fail on test4',diff
      test2 = .false.
    endif

    call test_5
    call pstar
    call ustar
    utoro = 8.68975_num
    diff = 2.0_num * (us - utoro) / (us + utoro)
    if (verbose) print *, 'ustar-test5 and diff',us,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *, 'ustar fail on test4',diff
      test2 = .false.
    endif

    if (test2) then
      print *,'passed ustar test'
    else
      print *,'failed ustar test'
      test_star = .false.
    endif
  end subroutine test_ustar

  subroutine test_rhostar

    real (num) :: diff, ltoro, rtoro
    logical :: test3 = .true.

    call test_1
    call pstar
    call ustar
    call rhostar

    ltoro = 0.42632_num
    diff = rhosl-ltoro
!    diff = 2.0_num * (rhosl - ltoro) / (rhosl + ltoro)
    if (verbose) print *, 'rhostar left-test1 and diff',rhosl,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test1 left'
      test3 = .false.
    endif

    rtoro = 0.26557_num
    diff = rhosr-rtoro
!   diff = 2.0_num * (rhosr - rtoro) / (rhosr + rtoro)
    if (verbose) print *, 'rhostar right -test1 and diff',rhosr,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test1 right'
      test3 = .false.
    endif

    call test_2
    call pstar
    call ustar
    call rhostar

    ltoro = 0.02185_num
    diff = rhosl-ltoro
    !diff = 2.0_num * (rhosl - ltoro) / (rhosl + ltoro)
    if (verbose) print *, 'rhostar left-test2 and diff',rhosl,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test2 left'
      test3 = .false.
    endif

    rtoro = ltoro
    diff = rhosr-rtoro
!    diff = 2.0_num * (rhosr - rtoro) / (rhosr + rtoro)
    if (verbose) print *, 'rhostar right -test2 and diff',rhosr,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test2 right'
      test3 = .false.
    endif

    call test_3
    call pstar
    call ustar
    call rhostar

    ltoro = 0.57506_num
    diff = rhosl-ltoro
    !diff = 2.0_num * (rhosl - ltoro) / (rhosl + ltoro)
    if (verbose) print *, 'rhostar left-test3 and diff',rhosl,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test3 left'
      test3 = .false.
    endif

    rtoro = 5.99924_num
    diff = rhosr-rtoro
!    diff = 2.0_num * (rhosr - rtoro) / (rhosr + rtoro)
    if (verbose) print *, 'rhostar right -test3 and diff',rhosr,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test3 right'
      test3 = .false.
    endif

    call test_4
    call pstar
    call ustar
    call rhostar

    ltoro = 5.99242_num
    diff = rhosl-ltoro
    !diff = 2.0_num * (rhosl - ltoro) / (rhosl + ltoro)
    if (verbose) print *, 'rhostar left-test4 and diff',rhosl,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test4 left'
      test3 = .false.
    endif

    rtoro = 0.57511_num
    diff = rhosr-rtoro
!    diff = 2.0_num * (rhosr - rtoro) / (rhosr + rtoro)
    if (verbose) print *, 'rhostar right -test4 and diff',rhosr,diff
    if (abs(diff) > 5e-6_num) then
      if (verbose) print *,  'rhostar fail on test4 rightt'
      test3 = .false.
    endif

    call test_5
    call pstar
    call ustar
    call rhostar

    ltoro = 14.2823_num
    diff = rhosl-ltoro
    !diff = 2.0_num * (rhosl - ltoro) / (rhosl + ltoro)
    if (verbose) print *, 'rhostar left-test5 and diff',rhosl,diff
    if (abs(diff) > 5e-5_num) then
      if (verbose) print *,  'rhostar fail on test5 left'
      test3 = .false.
    endif

    rtoro = 31.0426_num
    diff = rhosr-rtoro
!    diff = 2.0_num * (rhosr - rtoro) / (rhosr + rtoro)
    if (verbose) print *, 'rhostar right -test5 and diff',rhosr,diff
    if (abs(diff) > 5e-5_num) then !bit higher tol
      !Toro gives only to certain sf's - by inspection this is
      !acceptable within roundoff
      if (verbose) print *,  'rhostar fail on test5 right'
      test3 = .false.
    endif
    if (test3) then
      print *,'passed rhostar test'
    else
      print *,'failed rhostar test'
      test_star = .false.
    endif
  endsubroutine test_rhostar

  ! predefined initial conditions

  subroutine test_1 !Sod's test
    rhol = 1.0_num
    ul = 0.0_num
    pl = 1.0_num
    rhor = 0.125_num
    ur = 0.0_num
    pr = 0.10_num
    gamma = 1.4_num
    call constants !reset constants
  end subroutine test_1

  subroutine test_1_r !Reverse Sod's test for symmetry test in p*
    rhor = 1.0_num
    ur = 0.0_num
    pr = 1.0_num
    rhol = 0.125_num !1.0_num / 8.0_num
    ul = 0.0_num
    pl = 0.10_num
    gamma = 1.4_num
    call constants !reset constants
    t = 0.25_num
  end subroutine test_1_r

  subroutine test_2 !1,2,3 problem
    rhol = 1.0_num
    ul = -2.0_num
    pl = 0.4_num
    rhor = 1.0_num
    ur = 2.0_num
    pr = 0.4_num
    gamma = 1.4_num
    call constants !reset constants
    t = 0.15_num
  end subroutine test_2

  subroutine test_3 !left half blast
    rhol = 1.0_num
    ul = 0.0_num
    pl = 1000.0_num
    rhor = 1.0_num
    ur = 0.0_num
    pr = 0.01_num
    gamma = 1.4_num
    call constants !reset constants
    t = 0.012_num
  end subroutine test_3

  subroutine test_4 !right half
    rhol = 1.0_num
    ul = 0.0_num
    pl = 0.01_num
    rhor = 1.0_num
    ur = 0.0_num
    pr = 100.0_num
    gamma = 1.4_num
    call constants !reset constants
    t = 0.035_num
  end subroutine test_4

  subroutine test_5 !both halfs of the W+C blast
    rhol = 5.99924_num
    ul = 19.5975_num
    pl = 460.894_num
    rhor= 5.99242_num
    ur = -6.19633_num
    pr = 46.0950_num
    gamma = 1.4_num
    call constants !reset constants
    t = 0.035_num
  end subroutine

end module tests

module diagnostics

  use shared_data

  implicit none

  contains

  subroutine do_io
    integer :: out_unit =10
    call execute_command_line("rm -rf x.dat rho.dat p.dat u.dat en.dat")
       !^ dont stream to existing
    open(out_unit, file="x.dat", access="stream")
    write(out_unit) x
    close(out_unit)
    open(out_unit, file="rho.dat", access="stream")
    write(out_unit) rho_x
    close(out_unit)
    open(out_unit, file="p.dat", access="stream")
    write(out_unit) p_x
    close(out_unit)
    open(out_unit, file="u.dat", access="stream")
    write(out_unit) u_x
    close(out_unit)
    open(out_unit, file="en.dat", access="stream")
    write(out_unit) en_x
    close(out_unit)
    call execute_command_line("python plot1.py")
  end subroutine do_io

end module diagnostics

program ers_euler_ideal_1d

  use shared_data
  use riemann
  use user
  use tests
  use diagnostics

  implicit none

  call test_starvals !check passes all numerical tests

  call initial_conditions
  call constants
  call control

!  print *, 'warning: overwrite user initial conditions with test case'
!  call test_5 !_1, _2, _3, _4, _5
!  uncomment above to qualitatively test
!  sampling routines against Figs 4.7- 4.11 in Toro

  call check_positivity

  call pstar
  call ustar
  call rhostar

  call sample

  call do_io

end program ers_euler_ideal_1d

max_stars_repo_name: JOThurgood/SimpleCFD

max_issues_repo_name: JOThurgood/SimpleCFD

max_forks_repo_name: JOThurgood/SimpleCFD


text:
program decl_double
double precision :: i
end program decl_double

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


!==============================================================================
module ModHeatEq
  ! This is a "module" subprogram.
  ! Modules are ways to share values and subroutines across many parts
  ! of a larger F90 code.
  ! This contains methods and variables to solve the heat equation for a simple
  ! initial condition.

  ! Begin the initialization/declaration part of the code.

  ! This ensures that NO values are implicitly declared.
  ! ALWAYS include implicit none!!!
  implicit none
  save  ! "save" can be used to ensure that variables retain their values
  ! between calls to this module.  If you don't "save", variables are
  ! reset between calls to this subprogram.

  ! Variable declarations.
  ! Statically typed- variable meanings can't change within any part
  ! of the program!
  logical :: DoTest = .true.
  integer :: nT, nX
  real    :: dt, dx, tLimit=0.2, xLim(2) = (/0,1/)

  ! If you don't know how big an array is upon declaration, we need to
  ! make them allocatable. This means we allocate their size at run time.
  real, allocatable :: Domain_II(:,:)
  real, allocatable :: xGrid(:), tGrid(:)

  ! Constants: the parameter attribute means they can't change their values.
  real, parameter :: cDiffusion = 1.0

  ! The "contains" statement says "everything below this point has access to
  ! the module variables and all functions/subroutines below can be shared
  ! through the module".
contains

  !============================================================================
  subroutine run_simulation
    ! Advance our simulation from the initial condition to the end of the time
    ! range.  This should only be called AFTER the simulation domain has
    ! been initialized using "init_sim".

    ! Note our local variable declarations that are not recognized elsewhere
    ! within the module.
    integer :: j
    real    :: r
    !------------------------------------------------------------------------
    write(*,*) 'Beginning Simulation.'

    ! Check for stability:
    if( .not.IsStable(dT, dx, cDiffusion) ) then
       ! If we're not stable, we need to tell the user and
       ! end the program.  Create a well formatted message:
       write(*,"(a,f10.8,' >',f10.8)")'ERROR: Stability criterion not met! ', &
            dt, (dx**2.0/(2*cDiffusion**2.0))
       stop  ! The stop statement ends the whole program!
    end if

    ! Integrate!
    write(*,*) '...integrating...'
    ! Calculate a coefficient...
    r= cDiffusion**2.0 * dt / dx**2.0

    ! For each point from now until the second to last time,
    ! calculate the next state.  See the full PDF for details.
    do j=1, nT-1
       ! Note the ampersand is the line continuation marker.
       Domain_II(2:nX-1,j+1) = (1.0-2.0*r) * Domain_II(2:nX-1, j) + &
            r * (Domain_II(1:nx-2,j) + Domain_II(3:nx,j))
    end do

    ! Write results to file.  We do this by calling an external subroutine
    ! (see write2d.f90).  Note the ampersand is the line continuation marker.
    call write2d('results.txt', nX, nT, Domain_II, &
         xGrid, tGrid, 'NMUM Chapter 10 Example 3')

  end subroutine run_simulation

  !============================================================================
  logical function IsStable(dt, dx, c)
    ! Given a dt, dx, and a diffusion coefficient, c, determine if the
    ! simulation is stable.  Return a logical indicating the result.

    ! Functions have 0 or more inputs but only ONE output.  The output
    ! is the variable declared with the function name (i.e., "IsStable").
    ! IsStable is sent back to caller, e.g., x = IsStable(dt, dx, c).

    real, intent(in) :: dt, dx, c
    !------------------------------------------------------------------------

    ! Note how we use the name of the function as a variable inside the function
    IsStable = .true.
    if ( dt > (dx**2.0/(2.0*c**2.0)) ) IsStable=.false.

    ! We don't need a "return" statement.  Functions return only one value
    ! every time.

  end function IsStable

  !============================================================================
  subroutine init_sim(DxIn, DtIn)

    ! Subroutines have 0 or more inputs and 0 or more outputs.  We must use
    ! variable attributes to set the "intent" of the variables.  Are they
    ! input, output, or both?  Inputs CAN NOT be changed within the
    ! subroutine!  This is a good way to prevent mistakes and is a good
    ! habit!

    ! Set up the simulation.
    ! Note that for all inputs and outputs, we need to set their
    ! intent. intent(in) are values that cannot be edited/changed inside
    ! the subroutine.  intent(out) are returned to the caller.
    ! intent(inout) are inputs that *can* be changed.  Pick carefully!
    real, intent(in) :: DxIn, DtIn

    ! We also need just one integer.
    integer:: i
    !------------------------------------------------------------------------
    write(*,*) '...initializing simulation...'

    ! Set dx, dt based on inputs.  Because we cannot change our inputs, we
    ! copy them to new memory locations.
    dx = DxIn
    dt = DtIn

    ! Determine size of arrays.  "ceiling" rounds to the next highest integer.
    nX = ceiling( (xLim(2)-xLim(1))/dx + 1.0 )
    nT = ceiling( tLimit/dt +  1.0 )

    if(DoTest) write(*,*) '   Grid size (nX, nT) = ', nX, nT

    ! Now that we know size of arrays, allocate our arrays.
    allocate(Domain_II(nX, nT))
    allocate(xGrid(nX))
    allocate(tGrid(nT))

    ! It's usually a good idea to fill matrices with zeros.
    ! Make a habit of doing this.
    Domain_II = 0
    xGrid     = 0

    ! Set up grid and initial conditions.
    ! See PDF for details.
    do i=1, nX
       xGrid(i)       = (i-1) * dx
       Domain_II(i,1) = 4.0*xGrid(i) - 4.0*xGrid(i)**2.0
    end do

    ! Create time grid.
    do i=1, nT
       tGrid(i) = real(i-1) * dT
    end do

    ! We are ready to go!

  end subroutine init_sim

  !============================================================================
  subroutine finalize_sim
    ! Finalize the simulation by deallocating arrays.  This is a good habit!

    !------------------------------------------------------------------------
    deallocate(Domain_II)
    deallocate(xGrid)
    deallocate(tGrid)

    ! ALWAYS deallocate your arrays!  Otherwise, you will get memory leaks.

    write(*,*)'Finished simulation.  Produced ', nX*nT, ' points.'
  end subroutine finalize_sim

  !============================================================================

end module ModHeatEq
!==============================================================================


text:
!##############################################################################
!# ****************************************************************************
!# <name> AllenCahn_basic </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains basic problem structures for the heat conduction
!# problem.
!# </purpose>
!##############################################################################

MODULE AllenCahn_basic

  use fsystem
  use linearsolver
  use boundary
  use matrixfilters
  use vectorfilters
  use triangulation
  use spatialdiscretisation
  use sortstrategy
  use timestepping
  use discretebc
  use linearsystemblock
  use multilevelprojection
  use filtersupport

  use collection
  use paramlist

  IMPLICIT NONE

!<types>

!<typeblock description="Type block defining all information about one level">

  TYPE t_ACproblem_lvl

    ! An object for saving the triangulation on the domain
    TYPE(t_triangulation) :: rtriangulation

    ! An object specifying the discretisation (trial/test functions,...)
    TYPE(t_blockDiscretisation), POINTER :: p_rdiscretisation

!~~~~~~~MCai
    ! An object specifying the discretisation (trial/test functions,...)
    TYPE(t_blockDiscretisation)  :: rdiscretisation
!~~~~~
    ! The static matrix (containing Laplace, convection,...) which does not
    ! change in time.
    TYPE(t_matrixBlock) :: rmatrixStatic

    TYPE(t_matrixBlock) :: rmatrixPoly

    TYPE(t_matrixBlock) :: rmatrixJacobian

    TYPE(t_matrixBlock) :: rmatrixConv

    ! The mass matrix
    TYPE(t_matrixBlock) :: rmatrixMass

    ! System matrix. May change during the time iteration
    TYPE(t_matrixBlock) :: rmatrix

    ! A variable describing the discrete boundary conditions.
    ! The variable points to NULL until the first boundary conditions
    ! are assembled.
    TYPE(t_discreteBC), POINTER :: p_rdiscreteBC => NULL()

  END TYPE

!</typeblock>

!<typeblock description="Configuration block for the time stepping.">

  TYPE t_ACproblem_nonst

    ! Configuration block of the time stepping scheme.
    TYPE(t_explicitTimeStepping)        :: rtimestepping

    ! Number of current time step
    INTEGER                             :: iiteration

    ! Maximum number of time steps
    INTEGER                             :: niterations

    ! Start time
    REAL(DP)                            :: dtimemin

    ! Current time
    REAL(DP)                            :: dtime

    ! Maximum time
    REAL(DP)                            :: dtimemax

  END TYPE

!</typeblock>

!<typeblock description="Application-specific type block for AllenCahn problem">

  TYPE t_ACproblem

    ! Minimum refinement level; = Level i in RlevelInfo
    INTEGER :: NLMIN

    ! Maximum refinement level
    INTEGER :: NLMAX

    ! An object for saving the domain:
    TYPE(t_boundary) :: rboundary

    ! A RHS vector on the finest level used for solving linear systems
    TYPE(t_vectorBlock) :: rrhs

    ! A solver node that accepts parameters for the linear solver
    TYPE(t_linsolNode), POINTER :: p_rsolverNode

    ! An interlevel projection structure for changing levels
    TYPE(t_interlevelProjectionBlock) :: rprojection

    ! A filter chain to filter vectors during the solution process
    TYPE(t_filterChain), DIMENSION(1) :: RfilterChain

    ! A parameter block for everything that controls the time dependence.
    TYPE(t_ACproblem_nonst) :: rtimedependence

    ! An array of t_ACproblem_lvl structures, each corresponding
    ! to one level of the discretisation. There is currently
    ! only one level supported, identified by NLMAX!
    TYPE(t_ACproblem_lvl), DIMENSION(:), POINTER :: RlevelInfo

    ! A collection object that saves structural data and some
    ! problem-dependent information which is e.g. passed to
    ! callback routines.
    TYPE(t_collection) :: rcollection

  END TYPE

!</typeblock>

!</types>

END MODULE

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
module Cmnnext
  integer, save :: Cmnnext_c4
  real, save :: Cmnnext_c5
end module Cmnnext

max_stars_repo_name: raehik/camfort

max_issues_repo_name: camfort/camfort

max_forks_repo_name: raehik/camfort


text:
integer :: i, j = -1, k = 42
logical :: a

i = bit_size(j)       ! returns the number of bits in the given INTEGER variable

! bitwise boolean operations on integers
i = iand(k, j)        ! returns bitwise AND of K and J
i = ior(k, j)         ! returns bitwise OR of K and J
i = ieor(k, j)        ! returns bitwise EXCLUSIVE OR of K and J
i = not(j)            ! returns bitwise NOT of J

! single-bit integer/logical operations (bit positions are zero-based)
a = btest(i, 4)       ! returns logical .TRUE. if bit position 4 of I is 1, .FALSE. if 0
i = ibclr(k, 8)       ! returns value of K with 8th bit position "cleared" (set to 0)
i = ibset(k, 13)      ! returns value of K with 13th bit position "set" (set to 1)

! multi-bit integer operations
i = ishft(k, j)       ! returns value of K shifted by J bit positions, with ZERO fill
                      !    (right shift if J < 0 and left shift if J > 0).
i = ishftc(k, j)      ! returns value of K shifted CIRCULARLY by J bit positions
                      !    (right circular shift if J < 0 and left if J > 0)
i = ishftc(k, j, 20)  ! returns value as before except that ONLY the 20 lowest order
                      !    (rightmost) bits are circularly shifted
i = ibits(k, 7, 8)    ! extracts 8 contiguous bits from K starting at position 7 and
                      !    returns them as the rightmost bits of an otherwise
                      !    zero-filled integer. For non-negative K this is
                      !    arithmetically equivalent to:   MOD((K / 2**7), 2**8)

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


program main
   use m_npy
   real    :: x(1000)
   real    :: step
   integer :: n
   real, allocatable :: P(:,:)

   step = 2.0 / (size(x) - 1)

   x(1) = -1.0
   do n = 2,size(x)
      x(n) = x(n-1) + step
   enddo

   call legendre_poly(10, x, P)

   call save_npy("leg.npy", P)

contains
   subroutine legendre_poly(n_max, x, P)
      implicit none
      integer, intent(in)              :: n_max
      real, intent(in)                 :: x(:)
      real, intent(inout), allocatable :: P(:,:)

      integer :: n

      if(allocated(P)) then
         if(size(P,1) /= size(x) .or. size(P,2) /= n_max+1) then
            deallocate(P)
         endif
      endif

      if(.not. allocated(P)) then
         allocate(P(size(x), 0:n_max), source=0.0)
      endif

      if (n_max >= 0) then
         P(:, 0) = 1.0
      endif

      if (n_max >= 1) then
         P(:, 1) = x
      endif

      do n = 1, n_max - 1
         P(:, n + 1) = ((2*n + 1)*x*P(:, n) - n*P(:, n - 1))/(n + 1.0)
      enddo
   end subroutine legendre_poly
end program main