text:
c --------------------------------------------------------------------
      real*8 function LQ_QBB(s,t1,ms1)

      implicit none

      real*8 Pi,Nc,Cf,alphas,fac
      real*8 s,t1,u1,ms1

      Pi = 4.D0*atan(1.D0)
      Nc = 3.D0
      Cf = 4.D0/3.D0
      alphas = 1.D0

      u1 = - s - t1

      fac = 1.D0/(64*Pi*Nc**2)

      LQ_QBB =
     +  + Nc*Cf*Pi**2*alphas**2 * (  - 64*ms1**2*s**(-3) + 64*s**(-4)*
     +    t1*u1 )

      LQ_QBB = fac * LQ_QBB

      end

c --------------------------------------------------------------------
      real*8 function LQ_QB1(s,t1,ms1,mu)

      implicit none

      real*8 Pi,Nc,Cf,Nf,alphas
      real*8 s,t1,u1,ms1,sca,mu
      real*8 conalp,consof,LQ_QBB

      Pi = 4.D0*atan(1.D0)
      Nc = 3.D0
      Cf = 4.D0/3.D0
      Nf = 6.D0
      alphas = 1.D0

      sca = (mu/ms1)**2

      u1 = - s - t1

c renormalization of coupling with light flavours
      conalp = alphas/2.D0/Pi * LQ_QBB(s,t1,ms1)
     &       * ( 11.D0/3.D0*Nc - 2.D0/3.D0 * (Nf-1.D0)) * log(sca)

      consof = -alphas*Cf/2.D0/Pi * LQ_QBB(s,t1,ms1)
     &       * ( - 2.D0 * log(-t1/ms1**2)
     &           - 2.D0 * log(-u1/ms1**2)
     &           + 3.D0                  ) * log(sca)

      LQ_QB1 = conalp + consof

      end

c --------------------------------------------------------------------
      real*8 function LQ_QB2(s,t1,s4,ms1,del,s4max,mu)

      implicit none

      real*8 Pi,Cf,alphas
      real*8 s,t1,s4,del,s4max,ms1,sca,mu
      real*8 dldel1,LQ_QBB

      Pi = 4.D0*atan(1.D0)
      Cf = 4.D0/3.D0
      alphas = 1.D0

      sca = (mu/ms1)**2

      dldel1 = log(s4max/ms1**2)
     &       - ( s4max - del )/s4

      LQ_QB2 = -alphas*Cf/2.D0/Pi * LQ_QBB(s,t1,ms1)
     &      * 4.D0*dldel1 * log(sca)

      end

c --------------------------------------------------------------------
      real*8 function LQ_QB3(s,t1,s4,ms1,mu)

      implicit none

      real*8 Pi,Nc,Cf,fac,alphas
      real*8 s,t1,u1,s4,ms1,sca,mu
      real*8 COLO1(9)

      Pi = 4.D0*atan(1.D0)
      Cf = 4.D0/3.D0
      Nc = 3.D0
      alphas = 1.D0

      sca = (mu/ms1)**2

      u1 = s4 - s - t1

      fac = s4 /(s4+ms1**2) /8.D0 /(16*Pi**2)**2 /Nc**2

      COLO1(9) = -log(sca)

      LQ_QB3 = 0.D0
      LQ_QB3 = LQ_QB3 + COLO1(9)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-4)*t1**(-1)*u1 + 1024*ms1**2*s**(-4)*t1*u1**(-1) + 4096
     +    *ms1**2*s**(-4)*t1*s4**(-1) - 4096*ms1**2*s**(-4)*t1**2*
     +    s4**(-2) - 2048*ms1**2*s**(-4) - 1024*ms1**2*s**(-3)*t1**(-2)
     +    *u1 + 1024*ms1**2*s**(-3)*t1**(-1) - 1024*ms1**2*s**(-3)*t1*
     +    u1**(-2) - 4096*ms1**2*s**(-3)*t1*s4**(-2) + 1024*ms1**2*
     +    s**(-3)*u1**(-1) - 2048*ms1**2*s**(-3)*s4**(-1) - 1024*ms1**2
     +    *s**(-2)*t1**(-2) - 1024*ms1**2*s**(-2)*u1**(-2) - 1024*
     +    ms1**4*s**(-3)*t1**(-2) + 2048*ms1**4*s**(-3)*t1**(-1)*
     +    s4**(-1) - 1024*ms1**4*s**(-3)*u1**(-2) + 2048*ms1**4*s**(-3)
     +    *u1**(-1)*s4**(-1) - 4096*ms1**4*s**(-3)*s4**(-2) + 1024*
     +    s**(-4)*t1**(-1)*u1**2 + 3072*s**(-4)*t1 + 1024*s**(-4)*t1**2
     +    *u1**(-1) - 4096*s**(-4)*t1**2*s4**(-1) - 1024*s**(-4)*u1 +
     +    1024*s**(-3)*t1**(-1)*u1 + 1024*s**(-3)*t1*u1**(-1) - 4096*
     +    s**(-3)*t1*s4**(-1) )

      LQ_QB3 = fac * LQ_QB3

      end

c --------------------------------------------------------------------
      real*8 function LQ_QBD(s,t1,s4,ms1,del,s4max)

      implicit none

      real*8 Pi,Nc,Cf,alphas,fac
      real*8 s,t1,u1,ms1,s4,del,s4max
      real*8 xw,COLO1(9),dldel1,dldel2
      real*8 m2(5)

      Pi = 4.D0*atan(1.D0)
      Nc = 3.D0
      Cf = 4.D0/3.D0
      alphas = 1.D0

      u1 = - s - t1
      xw = sqrt(1.D0-4.D0*ms1**2/s)

      m2(1) = s
      m2(2) = t1
      m2(3) = u1
      m2(4) = s4
      m2(5) = ms1**2

      call LOGAS_LQ(m2,COLO1)

      dldel1 = log(s4max/ms1**2)
     &       - ( s4max - del )/s4
      dldel2 = log(s4max/ms1**2)**2
     &       - 2.D0* ( s4max - del )/s4 * log(s4/ms1**2)

      fac = 1.D0/(64*Pi*Nc**2)

      LQ_QBD =
     +  + Nc*Cf**2*Pi*alphas**3*dldel1 * ( 128*ms1**2*s**(-3) + 128*
     +    s**(-4)*t1**2 + 128*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + Nc*Cf**2*Pi*alphas**3*dldel2 * (  - 128*ms1**2*
     +    s**(-3) - 128*s**(-4)*t1**2 - 128*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(1)*Nc*Cf**2*Pi*alphas**3*dldel1 * ( -256*
     +    ms1**2*s**(-5)*t1**2*xw**(-1) - 256*ms1**2*s**(-4)*t1*
     +    xw**(-1) + 128*ms1**2*s**(-3)*xw**(-1) - 256*ms1**4*s**(-4)*
     +    xw**(-1) + 128*s**(-4)*t1**2*xw**(-1) + 128*s**(-3)*t1*
     +    xw**(-1) )
      LQ_QBD = LQ_QBD + COLO1(1)*Nc**2*Cf*Pi*alphas**3*dldel1 * ( 128*
     +    ms1**2*s**(-5)*t1**2*xw**(-1) + 128*ms1**2*s**(-4)*t1*
     +    xw**(-1) - 64*ms1**2*s**(-3)*xw**(-1) + 128*ms1**4*s**(-4)*
     +    xw**(-1) - 64*s**(-4)*t1**2*xw**(-1) - 64*s**(-3)*t1*xw**(-1)
     +     )
      LQ_QBD = LQ_QBD + COLO1(3)*Nc*Cf**2*Pi*alphas**3*dldel1 * ( -128*
     +    ms1**2*s**(-3) - 128*s**(-4)*t1**2 - 128*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(3)*Nc**2*Cf*Pi*alphas**3*dldel1 * ( 64*
     +    ms1**2*s**(-3) + 64*s**(-4)*t1**2 + 64*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(4)*Nc*Cf**2*Pi*alphas**3*dldel1 * ( -384*
     +    ms1**2*s**(-3) - 384*s**(-4)*t1**2 - 384*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(4)*Nc**2*Cf*Pi*alphas**3*dldel1 * ( 128*
     +    ms1**2*s**(-3) + 128*s**(-4)*t1**2 + 128*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(8)*Nc*Cf**2*Pi*alphas**3*dldel1 * ( 640*
     +    ms1**2*s**(-3) + 640*s**(-4)*t1**2 + 640*s**(-3)*t1 )
      LQ_QBD = LQ_QBD + COLO1(8)*Nc**2*Cf*Pi*alphas**3*dldel1 * ( -256*
     +    ms1**2*s**(-3) - 256*s**(-4)*t1**2 - 256*s**(-3)*t1 )

      LQ_QBD = fac * LQ_QBD

      end

c --------------------------------------------------------------------
      real*8 function LQ_QBV(s,t1,ms1,mt)

      implicit none

      real*8 Pi,Nc,Cf,Nf,alphas,fac,Nl
      real*8 s,t1,u1,ms1,mt
      real*8 ms12
      real*8 m2(2,6),softin(4)
      real*8 SK3B0A(5),SK3B0B(6),SK3B0C(2),SK3B0D(2,2),SK3BP(3)
      real*8 SK3C0A(8),SK3C0B(4),SK3C0C(5,2),SK3D0(7,2)
      real*8 SOF1(8)
      real*8 mur,dec,LQ_QBB

      Pi=4.D0*atan(1.D0)
      Nc = 3.D0
      Cf = 4.D0/3.D0
      Nf = 6.D0
      Nl = 1.D0
      alphas = 1.D0

      u1 = - s - t1
      fac = 1.D0/(64*Pi*Nc**2)

      m2(1,1) = s
      m2(1,2) = t1 + ms1**2
      m2(1,3) = u1 + ms1**2
      m2(1,4) = t1
      m2(1,5) = u1
      m2(2,1) = ms1**2
      m2(2,2) = 1.D0
      m2(2,3) = 1.D0
      m2(2,4) = 1.D0
      m2(2,5) = mt**2
      m2(2,6) = ms1**2

      call SCALAR_ARRAY_LQ_B(m2,SK3B0A,SK3B0B,SK3B0C,SK3B0D,SK3BP)
      call SCALAR_ARRAY_LQ_C(m2,SK3C0A,SK3C0B,SK3C0C,SK3D0)

c      write(99,*)'test1 ',m2

      softin(1) = s
      softin(2) = t1
      softin(3) = u1
      softin(4) = ms1**2

      call SOFT_ARRAY_LQ(softin,SOF1)

      mur = ms1

c decoupling of the heavy flavours from alpha_s
      dec = 2.D0*alphas/Pi * LQ_QBB(s,t1,ms1)*
     &      (  Nl/24.D0        * log(mur**2/ms1**2)
     &       + 1.D0/6.D0       * log(mur**2/mt**2)  )

      ms12 = ms1**2

      LQ_QBV =
     +  + Nc*Cf*Pi*alphas**3 * (  - 64./9.*Nf*ms1**2*s**(-3) - 64./9.*
     +    Nf*s**(-4)*t1**2 - 64./9.*Nf*s**(-3)*t1 + 32./9.*Nl*ms1**2*
     +    s**(-3) + 32./9.*Nl*s**(-4)*t1**2 + 32./9.*Nl*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + Nc*Cf**2*Pi*alphas**3 * ( 64*ms1**2*s**(-3)+64*
     +    s**(-4)*t1**2 + 64*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + Nc**2*Cf*Pi*alphas**3 * (  - 320./9.*ms1**2*
     +    s**(-3) - 320./9.*s**(-4)*t1**2 - 320./9.*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + SK3B0A(2)*Nc*Cf*Pi*alphas**3 * ( -128./3.*mt**2
     +    *ms1**2*s**(-4) - 128./3.*mt**2*s**(-5)*t1**2 - 128./3.*mt**2
     +    *s**(-4)*t1 )
      LQ_QBV = LQ_QBV + SK3B0B(1)*Nc*Cf*Pi*alphas**3*( 64./3.*Nf*ms1**2
     +    *s**(-3) + 64./3.*Nf*s**(-4)*t1**2 + 64./3.*Nf*s**(-3)*t1 -
     +    64./3.*ms1**2*s**(-3) - 64./3.*s**(-4)*t1**2 - 64./3.*s**(-3)
     +    *t1 )
      LQ_QBV = LQ_QBV + SK3B0B(1)*Nc*Cf**2*Pi*alphas**3 * ( 96*ms1**2*
     +    s**(-3) + 96*s**(-4)*t1**2 + 96*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + SK3B0B(1)*Nc**2*Cf*Pi*alphas**3 * ( 32*
     +    (s-4*ms12)**(-1)*s**(-3)*t1**2 + 32*(s-4*ms12)**(-1)*s**(-2)*
     +    t1 + 8*(s-4*ms12)**(-1)*s**(-1) - 256./3.*ms1**2*s**(-3) -
     +    256./3.*s**(-4)*t1**2 - 256./3.*s**(-3)*t1 - 8*s**(-2) )
      LQ_QBV = LQ_QBV + SK3B0B(2)*Nc*Cf*Pi*alphas**3 * (  - 64./3.*Nl*
     +    ms1**2*s**(-5)*t1**2 - 64./3.*Nl*ms1**2*s**(-4)*t1 + 16./3.*
     +    Nl*ms1**2*s**(-3) - 64./3.*Nl*ms1**4*s**(-4) + 16./3.*Nl*
     +    s**(-4)*t1**2 + 16./3.*Nl*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + SK3B0B(2)*Nc*Cf**2*Pi*alphas**3 * ( 64*
     +    (s-4*ms12)**(-1)*s**(-3)*t1**2 + 64*(s-4*ms12)**(-1)*s**(-2)*
     +    t1 + 16*(s-4*ms12)**(-1)*s**(-1) + 64*ms1**2*s**(-3) + 64*
     +    s**(-4)*t1**2 + 64*s**(-3)*t1 - 16*s**(-2) )
      LQ_QBV = LQ_QBV + SK3B0B(2)*Nc**2*Cf*Pi*alphas**3 * (  - 32*
     +    (s-4*ms12)**(-1)*s**(-3)*t1**2 - 32*(s-4*ms12)**(-1)*s**(-2)*
     +    t1 - 8*(s-4*ms12)**(-1)*s**(-1) - 32*ms1**2*s**(-3) - 32*
     +    s**(-4)*t1**2 - 32*s**(-3)*t1 + 8*s**(-2) )
      LQ_QBV = LQ_QBV + SK3B0B(4)*Nc*Cf*Pi*alphas**3 * ( 128./3.*mt**2*
     +    ms1**2*s**(-4) + 128./3.*mt**2*s**(-5)*t1**2 + 128./3.*mt**2*
     +    s**(-4)*t1 + 64./3.*ms1**2*s**(-3) + 64./3.*s**(-4)*t1**2 +
     +    64./3.*s**(-3)*t1 )
      LQ_QBV = LQ_QBV + SK3B0C(1)*Nc*Cf**2*Pi*alphas**3 * (  - 64*
     +    (s-4*ms12)**(-1)*s**(-3)*t1**2 - 64*(s-4*ms12)**(-1)*s**(-2)*
     +    t1 - 16*(s-4*ms12)**(-1)*s**(-1) - 64*ms1**2*s**(-3) - 64*
     +    s**(-4)*t1**2 - 64*s**(-3)*t1 + 16*s**(-2) )
      LQ_QBV = LQ_QBV + SK3BP(1)*Nc*Cf**2*Pi*alphas**3 * ( -128*ms1**2*
     +    s**(-4)*t1**2 - 128*ms1**2*s**(-3)*t1 - 128*ms1**4*s**(-3) )
      LQ_QBV = LQ_QBV + SK3C0A(1)*Nc*Cf**2*Pi*alphas**3 * ( 64*ms1**2*
     +    s**(-2) + 64*s**(-3)*t1**2 + 192*s**(-2)*t1 + 64*s**(-1) )
      LQ_QBV = LQ_QBV + SK3C0A(1)*Nc**2*Cf*Pi*alphas**3 * ( -48*ms1**2*
     +    s**(-2) - 64*s**(-3)*t1**2 - 112*s**(-2)*t1 - 32*s**(-1) )
      LQ_QBV = LQ_QBV + SK3C0B(1)*Nc*Cf**2*Pi*alphas**3 * ( -128*ms1**2
     +    *s**(-4)*t1**2 - 128*ms1**2*s**(-3)*t1 + 64*ms1**2*s**(-2) -
     +    128*ms1**4*s**(-3) + 64*s**(-3)*t1**2 + 64*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SK3C0B(1)*Nc**2*Cf*Pi*alphas**3 * ( 64*ms1**2*
     +    s**(-4)*t1**2 + 64*ms1**2*s**(-3)*t1 - 32*ms1**2*s**(-2) + 64
     +    *ms1**4*s**(-3) - 32*s**(-3)*t1**2 - 32*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SK3C0B(2)*Nc*Cf**2*Pi*alphas**3 * ( -256*ms1**2
     +    *s**(-3)*t1 - 128*ms1**2*s**(-2) + 128*s**(-2)*t1 + 64*
     +    s**(-1) )
      LQ_QBV = LQ_QBV + SK3C0B(2)*Nc**2*Cf*Pi*alphas**3 * ( 16*
     +    (s-4*ms12)**(-1)*s**(-2)*t1**2 + 16*(s-4*ms12)**(-1)*s**(-1)*
     +    t1 + 4*(s-4*ms12)**(-1) + 96*ms1**2*s**(-3)*t1 + 32*ms1**2*
     +    s**(-2) - 48*s**(-3)*t1**2 - 96*s**(-2)*t1 - 36*s**(-1) )
      LQ_QBV = LQ_QBV + SK3C0C(1,1)*Nc*Cf**2*Pi*alphas**3 *(128*ms1**2*
     +    s**(-3)*t1 + 128*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SK3C0C(1,1)*Nc**2*Cf*Pi*alphas**3 * (  - 32*
     +    ms1**2*s**(-3)*t1 - 32*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SK3C0C(1,2)*Nc*Cf**2*Pi*alphas**3 *(128*ms1**2*
     +    s**(-3)*t1 + 128*ms1**2*s**(-2) - 128*s**(-3)*t1**2 - 256*
     +    s**(-2)*t1 - 128*s**(-1) )
      LQ_QBV = LQ_QBV + SK3C0C(1,2)*Nc**2*Cf*Pi*alphas**3 * (  - 64*
     +    ms1**2*s**(-3)*t1 - 64*ms1**2*s**(-2) + 64*s**(-3)*t1**2 +
     +    128*s**(-2)*t1 + 64*s**(-1) )
      LQ_QBV = LQ_QBV + SK3D0(1,1)*Nc*Cf**2*Pi*alphas**3 * ( 64*ms1**2*
     +    s**(-2)*t1 + 128*s**(-3)*t1**3 + 64*s**(-2)*t1**2 )
      LQ_QBV = LQ_QBV + SK3D0(1,1)*Nc**2*Cf*Pi*alphas**3 * ( -16*ms1**2
     +    *s**(-2)*t1 - 32*s**(-3)*t1**3 - 16*s**(-2)*t1**2 )
      LQ_QBV = LQ_QBV + SK3D0(1,2)*Nc*Cf**2*Pi*alphas**3 * ( 64 + 64*
     +    ms1**2*s**(-2)*t1 + 64*ms1**2*s**(-1) + 128*s**(-3)*t1**3 +
     +    320*s**(-2)*t1**2 + 256*s**(-1)*t1 )
      LQ_QBV = LQ_QBV + SK3D0(1,2)*Nc**2*Cf*Pi*alphas**3 * ( -32 - 32*
     +    ms1**2*s**(-2)*t1 - 32*ms1**2*s**(-1) - 64*s**(-3)*t1**3 -
     +    160*s**(-2)*t1**2 - 128*s**(-1)*t1 )
      LQ_QBV = LQ_QBV + SOF1(1)*Nc*Cf**2*Pi*alphas**3 * ( 512*ms1**2*
     +    s**(-4)*t1**2 + 512*ms1**2*s**(-3)*t1 - 256*ms1**2*s**(-2) +
     +    512*ms1**4*s**(-3) - 256*s**(-3)*t1**2 - 256*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(1)*Nc**2*Cf*Pi*alphas**3 * ( -256*ms1**2*
     +    s**(-4)*t1**2 - 256*ms1**2*s**(-3)*t1 + 128*ms1**2*s**(-2) -
     +    256*ms1**4*s**(-3) + 128*s**(-3)*t1**2 + 128*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(2)*Nc*Cf**2*Pi*alphas**3 * ( 256*ms1**2*
     +    s**(-4)*t1**2 + 256*ms1**2*s**(-3)*t1 + 256*ms1**4*s**(-3) )
      LQ_QBV = LQ_QBV + SOF1(3)*Nc*Cf**2*Pi*alphas**3 * ( 256*ms1**2*
     +    s**(-4)*t1**2 + 256*ms1**2*s**(-3)*t1 + 256*ms1**4*s**(-3) )
      LQ_QBV = LQ_QBV + SOF1(4)*Nc*Cf**2*Pi*alphas**3 * ( -512*ms1**2*
     +    s**(-3)*t1 - 512*ms1**2*s**(-2) - 512*s**(-4)*t1**3 - 1024*
     +    s**(-3)*t1**2 - 512*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(4)*Nc**2*Cf*Pi*alphas**3 * ( 256*ms1**2*
     +    s**(-3)*t1 + 256*ms1**2*s**(-2) + 256*s**(-4)*t1**3 + 512*
     +    s**(-3)*t1**2 + 256*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(5)*Nc*Cf**2*Pi*alphas**3 * ( -512*ms1**2*
     +    s**(-3)*t1 - 512*s**(-4)*t1**3 - 512*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SOF1(5)*Nc**2*Cf*Pi*alphas**3 * ( 128*ms1**2*
     +    s**(-3)*t1 + 128*s**(-4)*t1**3 + 128*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SOF1(6)*Nc*Cf**2*Pi*alphas**3 * ( -512*ms1**2*
     +    s**(-3)*t1 - 512*s**(-4)*t1**3 - 512*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SOF1(6)*Nc**2*Cf*Pi*alphas**3 * ( 128*ms1**2*
     +    s**(-3)*t1 + 128*s**(-4)*t1**3 + 128*s**(-3)*t1**2 )
      LQ_QBV = LQ_QBV + SOF1(7)*Nc*Cf**2*Pi*alphas**3 * ( -512*ms1**2*
     +    s**(-3)*t1 - 512*ms1**2*s**(-2) - 512*s**(-4)*t1**3 - 1024*
     +    s**(-3)*t1**2 - 512*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(7)*Nc**2*Cf*Pi*alphas**3 * ( 256*ms1**2*
     +    s**(-3)*t1 + 256*ms1**2*s**(-2) + 256*s**(-4)*t1**3 + 512*
     +    s**(-3)*t1**2 + 256*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(8)*Nc*Cf**2*Pi*alphas**3 * ( -256*ms1**2*
     +    s**(-2) - 256*s**(-3)*t1**2 - 256*s**(-2)*t1 )
      LQ_QBV = LQ_QBV + SOF1(8)*Nc**2*Cf*Pi*alphas**3 * ( 128*ms1**2*
     +    s**(-2) + 128*s**(-3)*t1**2 + 128*s**(-2)*t1 )

      LQ_QBV = fac * LQ_QBV  + dec

      end

c --------------------------------------------------------------------
      real*8 function LQ_QBH(s,t1,s4,ms1)

      implicit none

      real*8 Pi,Nc,Cf,alphas,fac
      real*8 s,t1,u1,ms1,s4
      real*8 hardin(5)
      real*8 ANG4(71),COLO1(9)

      Pi = 4.D0*atan(1.D0)
      Nc = 3.D0
      Cf = 4.D0/3.D0
      alphas = 1.D0

      u1 = s4 - s - t1
      fac = s4 /(s4+ms1**2) /8.D0 /(16*Pi**2)**2 /Nc**2

      hardin(1) = s
      hardin(2) = t1
      hardin(3) = u1
      hardin(4) = s4
      hardin(5) = ms1**2

      call LOGAS_LQ(hardin,COLO1)
      call ANGULAR_ARRAY_LQ(hardin,ANG4)

      LQ_QBH =
     +  + Nc*Cf**2*Pi**4*alphas**3 * ( 1024*ms1**2*s**(-4)*t1**(-1)*u1
     +     + 1024*ms1**2*s**(-4)*t1*u1**(-1) - 4096*ms1**2*s**(-4)*t1*
     +    s4**(-1) + 2048*ms1**2*s**(-4)*t1**2*s4**(-2) + 2048*ms1**2*
     +    s**(-4) - 1024*ms1**2*s**(-3)*t1**(-2)*u1 - 1024*ms1**2*
     +    s**(-3)*t1**(-1) - 1024*ms1**2*s**(-3)*t1*u1**(-2) + 2048*
     +    ms1**2*s**(-3)*t1*s4**(-2) - 1024*ms1**2*s**(-3)*u1**(-1) -
     +    1024*ms1**2*s**(-2)*t1**(-2) - 1024*ms1**2*s**(-2)*u1**(-2)
     +     - 1024*ms1**4*s**(-3)*t1**(-2) - 1024*ms1**4*s**(-3)*
     +    u1**(-2) + 2048*ms1**4*s**(-3)*s4**(-2) + 1024*s**(-4)*
     +    t1**(-1)*u1**2 + 1024*s**(-4)*t1 + 1024*s**(-4)*t1**2*
     +    u1**(-1) + 1024*s**(-4)*u1 + 1024*s**(-3)*t1**(-1)*u1 + 1024*
     +    s**(-3)*t1*u1**(-1) + 1024*s**(-3)*t1*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(9)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2*
     +    s**(-4) )
      LQ_QBH = LQ_QBH + ANG4(10)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3) )
      LQ_QBH = LQ_QBH + ANG4(11)*Nc**2*Cf*Pi**4*alphas**3 * ( -1024*
     +    ms1**2*s**(-5)*s4*(t1+u1)**(-1) + 1024*ms1**2*s**(-5) )
      LQ_QBH = LQ_QBH + ANG4(12)*Nc*Cf**2*Pi**4*alphas**3 *(2048*ms1**2
     +    *s**(-4)*t1*s4**(-1) - 4096*ms1**2*s**(-4)*s4*(t1+u1)**(-1)
     +     + 2048*ms1**2*s**(-4) + 1024*ms1**2*s**(-3)*s4**(-1) - 1024*
     +    s**(-3)*t1*s4**(-1) + 512*s**(-3) - 512*s**(-2)*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(12)*Nc**2*Cf*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-5)*t1*s4*(t1+u1)**(-1) - 1024*ms1**2*s**(-5)*t1 - 1024*
     +    ms1**2*s**(-4)*t1*s4**(-1) + 512*ms1**2*s**(-4)*s4*
     +    (t1+u1)**(-1) - 512*ms1**2*s**(-4) - 512*ms1**2*s**(-3)*
     +    s4**(-1) + 512*s**(-3)*t1*s4**(-1) - 256*s**(-3) + 256*
     +    s**(-2)*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(13)*Nc**2*Cf*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-4) )
      LQ_QBH = LQ_QBH + ANG4(14)*Nc**2*Cf*Pi**4*alphas**3 *(2048*ms1**2
     +    *s**(-4)*u1 + 1024*ms1**2*s**(-3) )
      LQ_QBH = LQ_QBH + ANG4(15)*Nc**2*Cf*Pi**4*alphas**3 * (  - 1024*
     +    ms1**2*s**(-5)*s4*(t1+u1)**(-1) + 1024*ms1**2*s**(-5) )
      LQ_QBH = LQ_QBH + ANG4(16)*Nc*Cf**2*Pi**4*alphas**3 * (  - 4096*
     +    ms1**2*s**(-4)*s4*(t1+u1)**(-1) + 4096*ms1**2*s**(-4) )
      LQ_QBH = LQ_QBH + ANG4(16)*Nc**2*Cf*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-5)*t1*s4*(t1+u1)**(-1) - 1024*ms1**2*s**(-5)*t1 - 1024*
     +    ms1**2*s**(-4)*t1*s4**(-1) + 512*ms1**2*s**(-4)*s4*
     +    (t1+u1)**(-1) - 512*ms1**2*s**(-4) - 512*ms1**2*s**(-3)*
     +    s4**(-1) + 512*s**(-3)*t1*s4**(-1) - 256*s**(-3) + 256*
     +    s**(-2)*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(36)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**4
     +    *s**(-3) )
      LQ_QBH = LQ_QBH + ANG4(37)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3)*t1*s4**(-1) - 2048*ms1**2*s**(-3)*s4*(t1+u1)**(-1)
     +     + 2048*ms1**2*s**(-3) - 1024*ms1**2*s**(-2)*s4**(-1) + 2048*
     +    ms1**4*s**(-3)*s4**(-1) + 1024*s**(-3)*t1 - 1024*s**(-3)*u1
     +     - 512*s**(-2)*t1*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(37)*Nc**2*Cf*Pi**4*alphas**3 *(512*ms1**2*
     +    s**(-4)*t1*s4*(t1+u1)**(-1) - 512*ms1**2*s**(-4)*t1 - 512*
     +    ms1**2*s**(-3)*t1*s4**(-1) + 512*ms1**2*s**(-3)*s4*
     +    (t1+u1)**(-1) - 768*ms1**2*s**(-3) + 512*ms1**2*s**(-2)*
     +    s4**(-1) - 1024*ms1**4*s**(-3)*s4**(-1) + 512*s**(-4)*t1**2*
     +    s4*(t1+u1)**(-1) - 512*s**(-4)*t1**2 - 768*s**(-3)*t1 + 512*
     +    s**(-3)*u1 + 256*s**(-2)*t1*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(38)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-2) )
      LQ_QBH = LQ_QBH + ANG4(38)*Nc**2*Cf*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-4)*u1**2 + 1024*ms1**2*s**(-3)*u1 + 512*ms1**2*s**(-2)
     +     )
      LQ_QBH = LQ_QBH + ANG4(39)*Nc*Cf**2*Pi**4*alphas**3 *(4096*ms1**2
     +    *s**(-3)*t1*s4**(-1) - 2048*ms1**2*s**(-3)*s4*(t1+u1)**(-1)
     +     - 2048*ms1**2*s**(-3) + 2048*ms1**2*s**(-2)*s4**(-1) - 2048*
     +    s**(-2)*t1*s4**(-1) + 1024*s**(-2) - 1024*s**(-1)*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(39)*Nc**2*Cf*Pi**4*alphas**3 *(512*ms1**2*
     +    s**(-4)*t1*s4*(t1+u1)**(-1) - 1536*ms1**2*s**(-4)*t1 + 1024*
     +    ms1**2*s**(-4)*t1**2*s4**(-1) + 1024*ms1**2*s**(-4)*u1 - 512*
     +    ms1**2*s**(-3)*t1*s4**(-1) + 512*ms1**2*s**(-3)*s4*
     +    (t1+u1)**(-1) + 1024*ms1**2*s**(-3) - 256*ms1**2*s**(-2)*
     +    s4**(-1) + 512*s**(-4)*t1**2*s4*(t1+u1)**(-1) - 512*s**(-4)*
     +    t1**2 - 512*s**(-3)*t1 + 512*s**(-3)*t1**2*s4**(-1) + 256*
     +    s**(-3)*u1 + 1280*s**(-2)*t1*s4**(-1) - 512*s**(-2) + 512*
     +    s**(-1)*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(48)*Nc*Cf**2*Pi**4*alphas**3 * (  - 1024*
     +    ms1**2*s**(-2) - 512*s**(-3)*t1**2 - 512*s**(-3)*u1**2 - 512*
     +    s**(-2)*t1 - 512*s**(-2)*u1 )
      LQ_QBH = LQ_QBH + ANG4(48)*Nc**2*Cf*Pi**4*alphas**3 *(512*ms1**2*
     +    s**(-2) + 256*s**(-3)*t1**2 + 256*s**(-3)*u1**2 + 256*s**(-2)
     +    *t1 + 256*s**(-2)*u1 )
      LQ_QBH = LQ_QBH + ANG4(49)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3)*t1*s4*(t1+u1)**(-1) - 2048*ms1**2*s**(-3)*t1 + 2048*
     +    s**(-4)*t1**3*s4*(t1+u1)**(-1) - 2048*s**(-4)*t1**3 + 1024*
     +    s**(-3)*t1**2*s4*(t1+u1)**(-1) - 2048*s**(-3)*t1**2 )
      LQ_QBH = LQ_QBH + ANG4(49)*Nc**2*Cf*Pi**4*alphas**3 * (  - 256*
     +    ms1**2*s**(-3)*t1*s4*(t1+u1)**(-1) + 512*ms1**2*s**(-3)*t1 -
     +    512*s**(-4)*t1**3*s4*(t1+u1)**(-1) + 512*s**(-4)*t1**3 - 256*
     +    s**(-3)*t1**2*s4*(t1+u1)**(-1) + 512*s**(-3)*t1**2 )
      LQ_QBH = LQ_QBH + ANG4(50)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3)*t1*s4*(t1+u1)**(-1) - 1024*ms1**2*s**(-3)*t1 + 1024*
     +    ms1**2*s**(-3)*u1 - 1024*ms1**2*s**(-2) + 2048*s**(-4)*t1**3*
     +    s4*(t1+u1)**(-1) - 2048*s**(-4)*t1**3 + 2048*s**(-3)*t1*u1 -
     +    1024*s**(-3)*t1**2*s4*(t1+u1)**(-1) - 1024*s**(-3)*t1**2 -
     +    1024*s**(-3)*u1**2 + 1024*s**(-2)*t1 - 1024*s**(-2)*u1 )
      LQ_QBH = LQ_QBH + ANG4(50)*Nc**2*Cf*Pi**4*alphas**3 * (  - 512*
     +    ms1**2*s**(-3)*t1*s4*(t1+u1)**(-1) + 512*ms1**2*s**(-3)*t1 -
     +    512*ms1**2*s**(-3)*u1 + 512*ms1**2*s**(-2) - 1024*s**(-4)*
     +    t1**3*s4*(t1+u1)**(-1) + 1024*s**(-4)*t1**3 - 1024*s**(-3)*t1
     +    *u1 + 512*s**(-3)*t1**2*s4*(t1+u1)**(-1) + 512*s**(-3)*t1**2
     +     + 512*s**(-3)*u1**2 - 512*s**(-2)*t1 + 512*s**(-2)*u1 )
      LQ_QBH = LQ_QBH + ANG4(52)*Nc*Cf**2*Pi**4*alphas**3 *(512*ms1**2*
     +    s**(-1) )
      LQ_QBH = LQ_QBH + ANG4(54)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3)*t1*s4*(t1+u1)**(-1) - 1024*ms1**2*s**(-3)*t1 - 1024*
     +    ms1**2*s**(-2)*t1*s4**(-1) + 512*ms1**2*s**(-2) - 1024*ms1**2
     +    *s**(-1)*s4**(-1) + 2048*s**(-4)*t1**3*s4*(t1+u1)**(-1) -
     +    2048*s**(-4)*t1**3 + 1024*s**(-3)*t1**2*s4*(t1+u1)**(-1) -
     +    512*s**(-3)*t1**2 - 2048*s**(-3)*t1**3*s4**(-1) + 512*s**(-3)
     +    *u1**2 + 2560*s**(-2)*t1 - 5120*s**(-2)*t1**2*s4**(-1) - 1024
     +    *s**(-2)*u1 - 4096*s**(-1)*t1*s4**(-1) + 1024*s**(-1) - 1024*
     +    s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(54)*Nc**2*Cf*Pi**4*alphas**3 * (  - 256*
     +    ms1**2*s**(-3)*t1*s4*(t1+u1)**(-1) + 256*ms1**2*s**(-3)*t1 +
     +    512*ms1**2*s**(-2)*t1*s4**(-1) - 256*ms1**2*s**(-2) + 512*
     +    ms1**2*s**(-1)*s4**(-1) - 512*s**(-4)*t1**3*s4*(t1+u1)**(-1)
     +     + 512*s**(-4)*t1**3 + 256*s**(-3)*t1*u1 - 256*s**(-3)*t1**2*
     +    s4*(t1+u1)**(-1) - 256*s**(-3)*t1**2 + 1024*s**(-3)*t1**3*
     +    s4**(-1) - 256*s**(-3)*u1**2 - 1280*s**(-2)*t1 + 2560*s**(-2)
     +    *t1**2*s4**(-1) + 512*s**(-2)*u1 + 2048*s**(-1)*t1*s4**(-1)
     +     - 512*s**(-1) + 512*s4**(-1) )
      LQ_QBH = LQ_QBH + ANG4(56)*Nc*Cf**2*Pi**4*alphas**3 *(512*ms1**2*
     +    s**(-1) )
      LQ_QBH = LQ_QBH + ANG4(58)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-3)*t1*s4*(t1+u1)**(-1) - 1024*ms1**2*s**(-3)*t1 - 1024*
     +    ms1**2*s**(-2)*t1*s4**(-1) + 512*ms1**2*s**(-2) + 2048*
     +    s**(-4)*t1**3*s4*(t1+u1)**(-1) - 2048*s**(-4)*t1**3 - 1024*
     +    s**(-3)*t1**2*s4*(t1+u1)**(-1) + 1536*s**(-3)*t1**2 - 2048*
     +    s**(-3)*t1**3*s4**(-1) + 512*s**(-3)*u1**2 + 1024*s**(-2)*t1
     +     - 1024*s**(-2)*t1**2*s4**(-1) + 512*s**(-2)*u1 )
      LQ_QBH = LQ_QBH + ANG4(58)*Nc**2*Cf*Pi**4*alphas**3 * (  - 512*
     +    ms1**2*s**(-3)*t1*s4*(t1+u1)**(-1) + 512*ms1**2*s**(-3)*t1 +
     +    256*ms1**2*s**(-2)*t1*s4**(-1) - 1024*s**(-4)*t1**3*s4*
     +    (t1+u1)**(-1) + 1024*s**(-4)*t1**3 - 256*s**(-3)*t1*u1 + 512*
     +    s**(-3)*t1**2*s4*(t1+u1)**(-1) - 256*s**(-3)*t1**2 + 512*
     +    s**(-3)*t1**3*s4**(-1) - 512*s**(-2)*t1 + 256*s**(-2)*t1**2*
     +    s4**(-1) )
      LQ_QBH = LQ_QBH + COLO1(9)*Nc*Cf**2*Pi**4*alphas**3 *(1024*ms1**2
     +    *s**(-4)*t1**(-1)*u1 + 1024*ms1**2*s**(-4)*t1*u1**(-1) + 4096
     +    *ms1**2*s**(-4)*t1*s4**(-1) - 4096*ms1**2*s**(-4)*t1**2*
     +    s4**(-2) - 2048*ms1**2*s**(-4) - 1024*ms1**2*s**(-3)*t1**(-2)
     +    *u1 + 1024*ms1**2*s**(-3)*t1**(-1) - 1024*ms1**2*s**(-3)*t1*
     +    u1**(-2) - 4096*ms1**2*s**(-3)*t1*s4**(-2) + 1024*ms1**2*
     +    s**(-3)*u1**(-1) - 2048*ms1**2*s**(-3)*s4**(-1) - 1024*ms1**2
     +    *s**(-2)*t1**(-2) - 1024*ms1**2*s**(-2)*u1**(-2) - 1024*
     +    ms1**4*s**(-3)*t1**(-2) + 2048*ms1**4*s**(-3)*t1**(-1)*
     +    s4**(-1) - 1024*ms1**4*s**(-3)*u1**(-2) + 2048*ms1**4*s**(-3)
     +    *u1**(-1)*s4**(-1) - 4096*ms1**4*s**(-3)*s4**(-2) + 1024*
     +    s**(-4)*t1**(-1)*u1**2 + 3072*s**(-4)*t1 + 1024*s**(-4)*t1**2
     +    *u1**(-1) - 4096*s**(-4)*t1**2*s4**(-1) - 1024*s**(-4)*u1 +
     +    1024*s**(-3)*t1**(-1)*u1 + 1024*s**(-3)*t1*u1**(-1) - 4096*
     +    s**(-3)*t1*s4**(-1) )

      LQ_QBH = fac * LQ_QBH

      end

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text:
C/MEMBR ADD NAME=COSHIN,SSI=0
c
c fonction:  coshin
c fonction cosinus hyperbolique inverse de x
c en double precision
c acheve le 05/12/85
c ecrit par philippe touron
c
c
c sous programmes appeles: aucun
c
c
c     Copyright INRIA
      double precision function coshin(x)
      double precision x
      coshin=log(x+sqrt(x*x-1.0d+0))
      return
      end

ext: f

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
subroutine interp_tracer(t0,t1,IPP)
    use global,only: Nx,Ny,dtp,xyz,dic,djc,dkc,pi2c,pj2c,pk2c,Npts,theta,salt,gam,tsg,NPP,tt,phihyd,fn_PHIHYD
    implicit none
    integer*8,intent(in) :: t0,t1,IPP
    integer*8 :: i,j,k,ip
    real*8 :: tmp0,tmp1
    real*4 :: tmp(2,2,2)

   if (trim(fn_PHIHYD) .ne. '') then
   call load_PHIHYD(tt)
   endif

!$OMP PARALLEL DO PRIVATE(ip,i,j,k,tmp0,tmp1) SHARED(IPP,t0,t1,dic,djc,dkc,theta,salt,gam,tsg,dtp,phihyd) SCHEDULE(dynamic)
        do ip=1,npts

            if (xyz(ip,2,IPP)<Ny-1) then
                tmp0=0d0
                tmp1=0d0
                i=floor(mod(xyz(ip,1,IPP),real(Nx-1))-0.5d0)
                j=floor(xyz(ip,2,IPP)-0.5d0)
                k=floor(xyz(ip,3,IPP)-0.5d0)

                tmp=theta(i:i+1,j:j+1,k:k+1,t0)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp0)
                tmp=theta(i:i+1,j:j+1,k:k+1,t1)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp1)
                tsg(ip,1,IPP) = (tmp1-tmp0)*dtp+ tmp0

                tmp=salt(i:i+1,j:j+1,k:k+1,t0)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp0)
                tmp=salt(i:i+1,j:j+1,k:k+1,t1)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp1)
                tsg(ip,2,IPP) = (tmp1-tmp0)*dtp + tmp0

                tmp=gam(i:i+1,j:j+1,k:k+1,t0)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp0)
                tmp=gam(i:i+1,j:j+1,k:k+1,t1)
                where(tmp==0)tmp=maxval(tmp)
                call interp_trilinear(dic(ip,IPP),djc(ip,IPP),dkc(ip,IPP),tmp,tmp1)
                tsg(ip,3,IPP) = (tmp1-tmp0)*dtp + tmp0

                if (trim(fn_PHIHYD) .ne. '') then
                    call interp_bilinear(dic(ip,IPP),djc(ip,IPP),phihyd(i:i+1,j:j+1),tmp0)
                    tsg(ip,4,IPP) = tmp0
                endif

            endif
    enddo
!$OMP END PARALLEL DO

end subroutine interp_tracer

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
program test
  use, intrinsic :: iso_fortran_env, only : error_unit

  print *, "Error unit ", error_unit

end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
! Copyright (c) 2018 Etienne Descamps
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
!    this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation and/or
!    other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be
!    used to endorse or promote products derived from this software without specific prior
!    written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
! INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
! EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Module mlf_ode_class
  Use ieee_arithmetic
  Use iso_c_binding
  Use iso_fortran_env
  Use mlf_intf
  Use mlf_rsc_array
  Use mlf_models
  Use mlf_step_algo
  Use mlf_fun_intf
  Use mlf_utils
  Use mlf_supervised_model
  IMPLICIT NONE
  PRIVATE

  Public :: mlf_ode_algo_init, mlf_ode_algo_reinit, mlf_init_ode_model
  Type, Public, Abstract, Extends(mlf_step_obj) :: mlf_ode_algo
    real(c_double), pointer :: X0(:), X(:), K(:,:)
    real(c_double), pointer :: t0, t, tMax, hMax
    real(c_double), allocatable :: Cont(:,:)
    real(c_double) :: lastT
    integer(c_int64_t), pointer :: nFun
    class(mlf_ode_fun), pointer :: fun
    integer :: nK
  Contains
    procedure :: setFun => mlf_ode_algo_setFun
    procedure :: reinit => mlf_ode_algo_reinit
    procedure :: initODE => mlf_ode_algo_initODE
    procedure :: updateODE => mlf_ode_algo_updateODE
    procedure :: denseEvaluation => mlf_ode_algo_denseEvaluation
    procedure(mlf_ode_algo_initHandler), deferred :: initHandler
  End Type mlf_ode_algo

  Type, Public, Abstract, Extends(mlf_experience_model) :: mlf_ode_model
    class(mlf_ode_algo), pointer :: ode
  End Type mlf_ode_model

  Abstract Interface
    Integer Function mlf_ode_algo_initHandler(this, data_handler)
      import :: mlf_ode_algo, mlf_data_handler
      class(mlf_ode_algo), intent(inout), target :: this
      class(mlf_data_handler), intent(inout) :: data_handler
    End Function mlf_ode_algo_initHandler
  End Interface

Contains
  Integer Function mlf_ode_algo_updateODE(this, X0, t0, tMax, hMax) Result(info)
    class(mlf_ode_algo), intent(inout), target :: this
    real(c_double), intent(in), optional, target :: X0(:)
    real(c_double), intent(in), optional :: t0, tMax, hMax
    this%lastT = -HUGE(this%lastT)
    If(PRESENT(X0)) Then
      this%X0 = X0
      this%X = this%X0
    Endif
    If(PRESENT(t0)) Then
      this%t0 = t0
      this%t = this%t0
    Endif
    If(PRESENT(tMax)) this%tMax = tMax
    If(PRESENT(hMax)) this%hMax = hMax
  End Function mlf_ode_algo_updateODE

  Integer Function mlf_ode_algo_initODE(this, X0, t0, tMax, hMax) Result(info)
    class(mlf_ode_algo), intent(inout), target :: this
    real(c_double), intent(in), optional, target :: X0(:)
    real(c_double), intent(in), optional :: t0, tMax, hMax
    this%lastT = -HUGE(this%lastT)
    If(PRESENT(X0)) this%X0 = X0
    this%X = this%X0
    If(PRESENT(t0)) Then
      this%t0 = t0
    Else
      this%t0 = 0d0
    Endif
    this%t = this%t0
    If(PRESENT(tMax)) Then
      this%tMax = tMax
    Else
      this%tMax = HUGE(this%tMax)
    Endif
    If(PRESENT(hMax)) Then
      this%hMax = hMax
    Else
      this%hMax = HUGE(this%hMax)
    Endif
    info = 0
  End Function mlf_ode_algo_initODE

  ! Generic function for dense evaluation that use only X, X0 and their derivatives
  ! More advanced methods (such as ODE45) may rewrite this function more appropriatly
  Subroutine DefaultUpdateDense(this)
    class(mlf_ode_algo), intent(inout), target :: this
    real(c_double) :: dt
    dt = this%t-this%t0
    If(.NOT. ALLOCATED(this%Cont)) ALLOCATE(this%Cont(SIZE(this%X), 3))
    ASSOCIATE(A => this%Cont, KX0 => this%K(:,1), KX => this%K(:,this%nK))
      A(:,1) = this%X-this%X0
      A(:,2) = dt*KX0(:)-A(:,1)
      A(:,3) = -dt*KX(:)+A(:,1)-A(:,2)
    END ASSOCIATE
    this%lastT = this%t
  End Subroutine DefaultUpdateDense

  Subroutine mlf_ode_algo_denseEvaluation(this, t, Y, K)
    class(mlf_ode_algo), intent(inout), target :: this
    real(c_double), intent(in) :: t
    real(c_double), intent(out), target :: Y(:)
    real(c_double), intent(out), optional, target :: K(:)
    real(c_double) :: th, th1
    If(this%lastT < this%t) CALL DefaultUpdateDense(this)
    th = (t-this%t0)/(this%lastT-this%t0)
    th1 = th*(1d0-th)
    ASSOCIATE(X0 => this%X0, A => this%Cont)
      Y = X0+th*A(:,1)+th1*A(:,2)+th*th1*A(:,3)
      If(PRESENT(K)) Then
        K = A(:,1)+(1-2*th)*A(:,2)+th*(2-3*th)*A(:,3)
      Endif
    END ASSOCIATE
  End Subroutine mlf_ode_algo_denseEvaluation

  Integer Function mlf_ode_algo_reinit(this) result(info)
    class(mlf_ode_algo), intent(inout), target :: this
    info = mlf_step_obj_reinit(this)
    this%nFun = 0
    this%hMax = HUGE(this%hMax); this%tMax = HUGE(this%tMax)
    this%X0 = 0; this%X = 0; this%t0 = 0; this%t = 0
    this%lastT = -HUGE(this%lastT)
  End Function mlf_ode_algo_reinit

  Subroutine mlf_ode_algo_setFun(this, fun)
    class(mlf_ode_algo), intent(inout), target :: this
    class(mlf_ode_fun), intent(inout), target :: fun
    this%fun => fun
  End Subroutine mlf_ode_algo_setFun

  Integer Function mlf_ode_algo_init(this, numFields, nK, nX, data_handler) &
      Result(info)
    class(mlf_ode_algo), intent(inout), target :: this
    class(mlf_step_numFields), intent(inout) :: numFields
    class(mlf_data_handler), intent(inout), optional :: data_handler
    integer, intent(in), optional :: nK
    integer(c_int64_t), intent(in), optional :: nX
    integer(c_int64_t) :: N, nK2(2)
    CALL numFields%addFields(nRPar = 2, nRsc = 3, nIVar = 1, nRVar = 2)
    info = mlf_step_obj_init(this, numFields, data_handler)
    If(info < 0) RETURN
    N = -1; nK2 = -1
    If(PRESENT(nX)) N = nX
    info = this%add_rarray(numFields, N, this%X0, C_CHAR_"X0", data_handler = data_handler)
    If(info < 0) RETURN
    info = this%add_rarray(numFields, N, this%X, C_CHAR_"X", data_handler = data_handler)
    If(info < 0) RETURN
    If(PRESENT(nK)) nK2 = [N, INT(nK, KIND = 8)]
    info = this%add_rmatrix(numFields, nK2, this%K, C_CHAR_"K", data_handler = data_handler)
    If(info < 0) RETURN
    this%nK = INT(nK2(2), KIND = 4)
    ! Integer variables
    CALL this%addIVar(numFields, this%nFun, "nFun")
    ! Real parameters
    CALL this%addRPar(numFields, this%hMax, "hMax")
    CALL this%addRPar(numFields, this%tMax, "tMax")
    ! Real variables
    CALL this%addRVar(numFields, this%t, "t")
    CALL this%addRVar(numFields, this%t0, "t0")
  End Function mlf_ode_algo_init

  Integer Function mlf_init_ode_model(this, numFields, ode, data_handler) Result(info)
    class(mlf_ode_model), intent(inout), target :: this
    class(mlf_step_numFields), intent(inout) :: numFields
    class(mlf_ode_algo), pointer, intent(inout) :: ode
    class(mlf_data_handler), intent(inout), optional :: data_handler
    class(mlf_obj), pointer :: obj
    info = -1
    If(.NOT. ASSOCIATED(ode)) RETURN
    info = mlf_step_obj_init(this, numFields, data_handler)
    If(info<0) RETURN
    obj => ode
    If(PRESENT(data_handler)) Then
      info = ode%initHandler(data_handler%getSubObject(C_CHAR_"ode"))
      If(info < 0) RETURN
    Endif
    CALL this%add_subobject(C_CHAR_"ode", obj)
    this%ode => ode
    ode => NULL()
  End Function mlf_init_ode_model

End Module mlf_ode_class

max_stars_repo_name: Etdescamps/MlFortran

max_issues_repo_name: Etdescamps/MlFortran

max_forks_repo_name: Etdescamps/MlFortran


text:
	SUBROUTINE GPNPLY ( ictyp, line, loci, info, iret )
C************************************************************************
C* GPNPLY								*
C*									*
C*									*
C**									*
C* Log:									*
C* S. Jacobs/SSAI	12/91						*
C* C. Bailey/HPC	10/06		Changed call to IN_LINE		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	line, loci, info
C*
	CHARACTER	sys*1, plypts(20)*24
	REAL		tmpt(2), xplypt(20), yplypt(20)
	REAL		smooth, linfltr
	LOGICAL		scflag
C------------------------------------------------------------------------
C*	Get plot bounds in Normal coordinates.
C
        CALL GQBND ( 'N', xleft, ybot, xrght, ytop, ier )
C
C*	Parse the fill input.
C
	CALL IN_LINE ( line, 1, 1, icolor, itype, iwidth, ilabel,
     +		       smooth, linfltr, scflag, ier )
C
	IF  ( icolor .ne. 0 )  THEN
C
C*	    Set the coordinate system to be used.
C
	    sys = 'N'
	    IF  ( loci(1:1) .eq. '@' )  THEN
		sys = 'G'
		loci = loci(2:)
	    ELSE IF  ( loci(1:1) .eq. '#' )  THEN
		sys = 'M'
		loci = loci(2:)
	    END IF
C
C*	    Get the point(s) for the plot.
C
	    CALL ST_CLST ( loci, '/', ' ', 20, plypts, num, ier )
C
C*	    Get the line type.
C
	    CALL ST_NUMB ( info, ltype, ier )
C
	    DO  i = 1, num
		CALL ST_RLST ( plypts(i), ';', 0., 2, tmpt, num2,
     +			       ier )
		xplypt(i) = tmpt(1)
		yplypt(i) = tmpt(2)
C
C*		If using Normal coords, correct for PS device.
C
		IF  ( sys .eq. 'N' )  THEN
		    xplypt(i) = (xrght-xleft) * xplypt(i)
		    yplypt(i) = (ytop -ybot ) * yplypt(i)
		END IF
C
	    END DO
C
C*	    Set attributes.
C
	    CALL GSCOLR  ( icolor, ier )
	    CALL GSLINE  ( ltype, 0, iwidth, 0, ier )
C
C*	    Plot polygon.
C
            IF ( ictyp .eq. 1 ) THEN
                num = num + 1
                xplypt (num) = xplypt (1)
                yplypt (num) = yplypt (1)
                CALL GLINE ( sys, num, xplypt, yplypt, ier )
            ELSE
                CALL GFILL ( sys, num, xplypt, yplypt, ier )
            END IF
C
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Some of Daviss elections are considered mailonly elections, although strictly speaking, that is not the only way that one can vote in such elections.  The May 2011 Election was the first mailonly election, followed by the March 2012 Election and the March 2013 Election.

Here are the ways you can vote:

 1. By mail  before election day
 Your ballot may be dropped into a Mail Drop Boxes U.S. Mail Box.  It must be received in the Yolo County Elections Office no later than 8:00 PM on election day.  Postmarks do not count.

 1. In person  on or before election day
 Your ballot may be returned by you in person to the Yolo County Elections Office on or before election day or to any polling place in Yolo County, by 8:00 PM on election day.

 1. On or before election day
 If you are ill or disabled your ballot may be returned by anyone from the specified list to the Yolo County Elections Office on or before election day or to any polling place in Yolo County on election day.  In this event, you must complete the authorization shown on the Voter Return Envelope.

In the March 2013 Election election you may also vote in one of the following ways:
  Drop your ballot in the secure ballot box in the Yolo County Library  Mary L. Stephens Branch Stephens Library on 14th Street by 8:00 PM, March 5, 2013.
  Vote in person at the Veterans Memorial Center, 203 14th St. in Davis, from 7 a.m. to 8 p.m. March 5, 2013.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine creeps(inpc,textpart,nelcon,nmat,ntmat_,npmat_,
     &        plicon,nplicon,elcon,iplas,iperturb,nstate_,ncmat_,
     &        matname,irstrt,istep,istat,n,iline,ipol,inl,ipoinp,inp,
     &        ipoinpc)
!
!     reading the input deck: *CREEP
!
      implicit none
!
      logical iso
!
      character*1 inpc(*)
      character*80 matname(*)
      character*132 textpart(16)
!
      integer nelcon(2,*),nmat,ntmat_,ntmat,istep,npmat_,nstate_,
     &  n,key,i,j,iplas,iperturb(*),istat,nplicon(0:ntmat_,*),ncmat_,
     &  k,id,irstrt,iline,ipol,inl,ipoinp(2,*),inp(3,*),ipoinpc(0:*)
!
      real*8 temperature,elcon(0:ncmat_,ntmat_,*),t1l,
     &  plicon(0:2*npmat_,ntmat_,*)
!
      iso=.true.
      ntmat=0
!
      if((istep.gt.0).and.(irstrt.ge.0)) then
         write(*,*) '*ERROR reading *CREEP: *CREEP should be placed'
         write(*,*) '  before all step definitions'
         call exit(201)
      endif
!
      if(nmat.eq.0) then
         write(*,*) '*ERROR reading *CREEP: *CREEP should be preceded'
         write(*,*) '  by a *MATERIAL card'
         call exit(201)
      endif
!
!     check for anisotropic creep: assumes a ucreep routine
!
      if((nelcon(1,nmat).ne.2).and.(nelcon(1,nmat).ne.-51)) then
         write(*,*) '*ERROR reading *CREEP: *CREEP should be'
         write(*,*) '       preceded by an *ELASTIC,TYPE=ISO card'
         call exit(201)
      endif
!
!     if the *CREEP card is not preceded by a *PLASTIC card, a zero
!     yield surface is assumed
!
!        elastic isotropic
!        plasticity or no plasticity
!        creep (Norton or user): -52
!
         if(nelcon(1,nmat).ne.-51) then
!
!           elastic isotropic
!           no plasticity -> zero yield plasticity
!           creep (Norton or user)
!
            nplicon(0,nmat)=1
            nplicon(1,nmat)=2
            plicon(0,1,nmat)=0.d0
            plicon(1,1,nmat)=0.d0
            plicon(2,1,nmat)=0.d0
            plicon(3,1,nmat)=0.d0
            plicon(4,1,nmat)=10.d10
         endif
!
         iperturb(1)=3
         iplas=1
         nelcon(1,nmat)=-52
         nstate_=max(nstate_,13)
!
         do i=2,n
            if(textpart(i)(1:8).eq.'LAW=USER') then
               do j=1,nelcon(2,nmat)
                  elcon(3,j,nmat)=-1.d0
               enddo
               call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &              ipoinp,inp,ipoinpc)
               return
            else
               write(*,*)
     &              '*WARNING reading *CREEP: parameter not recognized:'
               write(*,*) '         ',
     &              textpart(i)(1:index(textpart(i),' ')-1)
               call inputwarning(inpc,ipoinpc,iline,
     &"*CREEP%")
            endif
         enddo
!
!        before interpolation: data are stored in positions 6-9:
!        A,n,m,temperature
!        after interpolation: data are stored in positions 3-5:
!        A,n,m
!
         do
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) exit
            ntmat=ntmat+1
            if(ntmat.gt.ntmat_) then
               write(*,*) '*ERROR reading *CREEP: increase ntmat_'
               call exit(201)
            endif
            do i=1,3
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &               elcon(i+5,ntmat,nmat)
               if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*CREEP%")
            enddo
            if(elcon(6,ntmat,nmat).le.0.d0) then
               write(*,*) '*ERROR reading *CREEP: parameter A'
               write(*,*) '       in the Norton law is nonpositive'
               call exit(201)
            endif
            if(elcon(7,ntmat,nmat).le.0.d0) then
               write(*,*) '*ERROR reading *CREEP: parameter n'
               write(*,*) '       in the Norton law is nonpositive'
               call exit(201)
            endif
            if(textpart(4)(1:1).ne.' ') then
               read(textpart(4)(1:20),'(f20.0)',iostat=istat)
     &               temperature
               if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*CREEP%")
            else
               temperature=0.d0
            endif
            elcon(9,ntmat,nmat)=temperature
         enddo
!
         if(ntmat.eq.0) then
            write(*,*) '*ERROR reading *CREEP: Norton law assumed,'
            write(*,*) '       yet no constants given'
            call exit(201)
         endif
!
!        interpolating the creep data at the elastic temperature
!        data points
!
         do i=1,nelcon(2,nmat)
            t1l=elcon(0,i,nmat)
            call ident2(elcon(9,1,nmat),t1l,ntmat,ncmat_+1,id)
            if(ntmat.eq.0) then
               continue
            elseif((ntmat.eq.1).or.(id.eq.0)) then
               elcon(3,i,nmat)=elcon(6,1,nmat)
               elcon(4,i,nmat)=elcon(7,1,nmat)
               elcon(5,i,nmat)=elcon(8,1,nmat)
            elseif(id.eq.ntmat) then
               elcon(3,i,nmat)=elcon(6,id,nmat)
               elcon(4,i,nmat)=elcon(7,id,nmat)
               elcon(5,i,nmat)=elcon(8,id,nmat)
            else
               do k=3,5
                  elcon(k,i,nmat)=elcon(k+3,id,nmat)+
     &               (elcon(k+3,id+1,nmat)-elcon(k+3,id,nmat))*
     &               (t1l-elcon(9,id,nmat))/
     &               (elcon(9,id+1,nmat)-elcon(9,id,nmat))
               enddo
            endif
         enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


    Subroutine x04ebfn(matrix, diag, m, n, a, lda, form, title, labrow, rlabs, &
      labcol, clabs, ncols, indent, errbuf, ifail)

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

!     Prints a general integer matrix.

!     .. Use Statements ..
      Use lapack_example_aux, Only: p01abf, x04abfn, x04bafn, x04cbzn
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: one = 1, zero = 0
      Character (6), Parameter :: srname = 'X04EBF'
!     .. Scalar Arguments ..
      Integer :: ifail, indent, lda, m, n, ncols
      Character (1) :: diag, labcol, labrow, matrix
      Character (200) :: errbuf
      Character (*) :: form, title
!     .. Array Arguments ..
      Integer :: a(lda, *)
      Character (*) :: clabs(*), rlabs(*)
!     .. Local Scalars ..
      Integer :: aa, cl, clbwid, cleft, clshft, cright, diagel, finis2, &
        finish, i, ierr, incols, indnt, j, k, lwid, maxel, minel, nd, nelems, &
        nels, nout, nrec, nsets, ntitle, numwid, offset, rlbwid, start, start2
      Logical :: generl, lower, prdiag
      Character (132) :: blanks, tilter
      Character (87) :: fmat
      Character (81) :: formt
      Character (10) :: icform, irform
      Character (266) :: infil2
      Character (133) :: infile
!     .. Local Arrays ..
      Character (80) :: rec(2)
!     .. Intrinsic Procedures ..
      Intrinsic :: char, len, max, min
!     .. Executable Statements ..

!-INT INTENT(IN) :: MATRIX, DIAG, M, N, A, LDA, FORM, TITLE, LABROW
!-INT INTENT(IN) :: RLABS, LABCOL, CLABS, NCOLS, INDENT
!-INT INTENT(INOUT) :: IFAIL
!-INT INTENT(OUT) :: ERRBUF
!     .. Executable Statements ..
      Continue
      errbuf(1:1) = char(0)
!
      ierr = 0
      generl = matrix == 'G' .Or. matrix == 'g'
      lower = matrix == 'L' .Or. matrix == 'l'
      If (ncols<=0 .Or. ncols>132) Then
        incols = 80
      Else
        incols = ncols
      End If
      If (indent<0 .Or. indent>=incols) Then
        indnt = 0
      Else
        indnt = indent
      End If
      incols = incols - indnt
      blanks = ' '
!
!     Check for incorrect arguments.
      If (.Not. generl .And. .Not. lower .And. matrix/='U' .And. matrix/='u') &
        Then
        ierr = 1
        nrec = 1
        Write (rec, Fmt=130) matrix
        Write (errbuf, Fmt=140) srname, 1, 99999, 1, matrix
      Else If ((labrow/='N' .And. labrow/='n' .And. labrow/='I' .And. &
          labrow/='i' .And. labrow/='C' .And. labrow/='c') .Or. &
          (labcol/='N' .And. labcol/='n' .And. labcol/='I' .And. &
          labcol/='i' .And. labcol/='C' .And. labcol/='c')) Then
        ierr = 6
        nrec = 2
        Write (rec, Fmt=190) labrow, labcol
        Write (errbuf, Fmt=200) srname, 6, 99996, 2, labrow, labcol
      Else If (m>lda) Then
        ierr = 3
        nrec = 1
        Write (rec, Fmt=150) m, lda
        Write (errbuf, Fmt=160) srname, 3, 99998, 2, m, lda
      Else If (.Not. generl) Then
        If (diag/='U' .And. diag/='u' .And. diag/='N' .And. diag/='n' .And. &
          diag/='B' .And. diag/='b') Then
          ierr = 2
          nrec = 1
          Write (rec, Fmt=170) matrix, diag
          Write (errbuf, Fmt=180) srname, 2, 99997, 2, matrix, diag
        End If
      End If
      If (ierr/=0) Go To 120
!
!     End of argument checking.
!
!     Get the advisory message unit number.
      Call x04abfn(0, nout)
!
      formt = form
      If (formt==' ') Then
!       Construct the FORMAT code of smallest field width that is
!       large enough to print all matrix elements.
!       First find the largest and smallest elements to be printed.
        maxel = zero
        minel = zero
        If (generl) Then
          Do j = 1, n
            Do i = 1, m
              aa = a(i, j)
              If (aa>maxel) maxel = aa
              If (aa<minel) minel = aa
            End Do
          End Do
        Else
          If (diag=='N' .Or. diag=='n') Then
            Do j = 1, min(m, n)
              aa = a(j, j)
              If (aa>maxel) maxel = aa
              If (aa<minel) minel = aa
            End Do
          End If
          Do j = 1, n
            If (lower) Then
              Do i = j + 1, m
                aa = a(i, j)
                If (aa>maxel) maxel = aa
                If (aa<minel) minel = aa
              End Do
            Else
              Do i = 1, j - 1
                aa = a(i, j)
                If (aa>maxel) maxel = aa
                If (aa<minel) minel = aa
              End Do
            End If
          End Do
        End If
!
        Write (infile, Fmt='(I40)') maxel
        Call x04cbzn(infile, start, finish)
        Write (infile, Fmt='(I40)') minel
        Call x04cbzn(infile, start2, finis2)
        k = max(finish-start, finis2-start2) + 2
        formt = 'I  '
        Write (formt(2:), Fmt='(I2)') k
      End If
!
!     Construct the format statement to be used internally.
      Call x04cbzn(formt, start, finish)
      If (finish-start+1>80) Then
!       The length of FORM is too great.
        ierr = 4
        nrec = 2
        Write (rec, Fmt=210) form(start:start+74)
        Write (errbuf, Fmt=220) srname, 4, 99995, 1, form(start:start+74)
        Go To 120
      End If
!
      fmat = '(999(' // formt(start:finish) // '))'
!
!     Decide how wide each column of numbers is going to be,
!     by writing the number 0 to an internal file and measuring
!     the width needed. Since the width may include trailing blanks,
!     we also write 0 twice with the same format, and compute
!     the required field width from the two widths. Note that if
!     FORM has more than one edit descriptor in it, with different
!     field widths, for example FORM = 'I12,I13', then it is
!     not possible to compute a sensible value, in which case
!     the columns of the output matrix will be a little skew.
      Write (infile, Fmt=fmat, Err=100) zero
      Call x04cbzn(infile, start, finish)
      Write (infil2, Fmt=fmat, Err=100) zero, zero
      Call x04cbzn(infil2, start, finis2)
      numwid = finis2 - finish
!     NUMWID is the width of a number as printed using FORMT.
      Go To 110
100   Continue
!     The format in FMAT caused an error when used to print a number.
      ierr = 5
      nrec = 2
      Write (rec, Fmt=230) form(1:min(75,len(form)))
      Write (errbuf, Fmt=240) srname, 5, 99994, 2, form(1:min(75,len(form)))
      Go To 120
110   Continue
!
!     What kind of row labelling is required?
      If (labrow=='N' .Or. labrow=='n') Then
!       No row labelling.
        rlbwid = 1
      Else If (labrow=='I' .Or. labrow=='i') Then
!       Numeric row labelling.
        Write (infile, Fmt='(I16)') m
        Call x04cbzn(infile, start, finish)
        rlbwid = finish - start + 2
        irform = '(I    )'
        Write (irform(3:6), Fmt='(I4)') rlbwid
      Else
!       User supplied row labelling.
        rlbwid = 1
        Do i = 1, m
          Call x04cbzn(rlabs(i), start, finish)
          rlbwid = max(rlbwid, finish-start+2)
        End Do
      End If
!
!     What kind of column labelling is required?
      If (labcol=='I' .Or. labcol=='i') Then
!       Numeric column labelling.
        Write (infile, Fmt='(I16)') n
        Call x04cbzn(infile, start, finish)
        clbwid = finish - start + 2
        icform = '(999I    )'
        Write (icform(6:9), Fmt='(I4)') numwid
      Else If (labcol/='N' .And. labcol/='n') Then
!       User supplied column labelling.
        clbwid = len(clabs(1))
      End If
!
      nelems = (incols-1-rlbwid)/numwid
      If (nelems<1) Then
        ierr = 7
        nrec = 2
        Write (rec, Fmt=250) incols + indnt, indnt
        Write (errbuf, Fmt=260) srname, 7, 99993, 2, incols + indnt, indnt
        Go To 120
      End If
!     NELEMS is the number of elements that can fit into INCOLS columns.
!
      nsets = (n-1)/nelems + 1
!     NSETS is the number of pieces that the matrix must be split into.
!
!     Print the title, splitting it up if more than INCOLS-1 characters.
      Call x04cbzn(title, start, finish)
      If (finish/=0) Then
        ntitle = (finish-1)/(incols-1) + 1
        Do i = 1, ntitle - 1
          tilter = blanks(1:indnt+1) // title((i-1)*(incols-1)+1:i*(incols-1))
          Call x04bafn(nout, tilter)
        End Do
        tilter = blanks(1:indnt+1) // title((ntitle-1)*(incols-1)+1:finish)
        Call x04bafn(nout, tilter)
      End If
!
!     Exit after printing the title if M or N is less than 1.
      If (m<1 .Or. n<1) Go To 120
!
!     Print the matrix, with row and column labels if requested.
      clshft = 0
!     CLSHFT is the offset into the current set of columns, when
!     the matrix cannot be printed in one go but has to be split.
      Do i = 1, nsets
        If (i==nsets) Then
          nels = n - (nsets-1)*nelems
        Else
          nels = nelems
        End If
        If (labcol=='I' .Or. labcol=='i') Then
!         Construct the numeric column labels.
          infile = ' '
          Write (infile(rlbwid+2+indnt:), Fmt=icform)(cl, cl=clshft+1, &
            clshft+nels)
        Else If (labcol/='N' .And. labcol/='n') Then
!         Process the user-supplied column labels.
          infile = ' '
          lwid = min(clbwid, numwid)
          offset = rlbwid + 1 + indnt
          Do k = clshft + 1, clshft + nels
            Call x04cbzn(clabs(k)(1:lwid), start, finish)
            If (start==0) Then
              start = 1
              finish = 1
            End If
            infile(offset+numwid-finish+start:offset+numwid+finish-start) &
              = clabs(k)(start:finish)
            offset = offset + numwid
          End Do
        End If
!       Output the column labels.
        If (labcol/='N' .And. labcol/='n') Then
          Call x04bafn(nout, infile(1:ncols))
        End If
!
!       Now print each row in turn.
        Do j = 1, m
          infile = ' '
!
!         Insert the row label.
          If (labrow=='I' .Or. labrow=='i') Then
            Write (infile(indnt+1:indnt+rlbwid), Fmt=irform) j
          Else If (labrow/='N' .And. labrow/='n') Then
            Call x04cbzn(rlabs(j), start, finish)
            If (start==0) Then
              start = 1
              finish = 1
            End If
            infile(indnt+rlbwid-finish+start:indnt+rlbwid) = rlabs(j) &
              (start:finish)
          End If
!
          If (generl) Then
!           General rectangular matrix.
            Write (infile(indnt+rlbwid+2:), Fmt=fmat)(a(j,cl), cl=clshft+1, &
              clshft+nels)
          Else
!           Upper or lower triangular matrix.
            nd = max(0, j-(i-1)*nelems)
!           ND is the position of the Jth row diagonal element.
            If (lower) Then
              cleft = clshft + 1
              cright = clshft + min(nd-1, nels)
            Else
              cleft = clshft + nd + 1
              cright = clshft + nels
            End If
!           CLEFT and CRIGHT are the leftmost and rightmost elements
!           of the current row to be printed, excluding the diagonal.
            prdiag = diag /= 'B' .And. diag /= 'b' .And. nd > 0 .And. &
              nd <= nels
!           PRDIAG is true if a diagonal element appears in the
!           current matrix row section, and it is to be printed.
            If (prdiag) Then
              If (diag=='U' .Or. diag=='u') Then
                diagel = one
              Else
                diagel = a(j, j)
              End If
            End If
!
            If (lower) Then
!             Lower triangular matrix.
              If (prdiag) Then
                If (cleft<=cright) Then
!                 Optimization bug on Alpha forces us to check
!                 CLEFT.LE.CRIGHT; DIAGEL didn't get printed
!                 out when CLEFT.GT.CRIGHT.
                  Write (infile(indnt+rlbwid+2:), Fmt=fmat)(a(j,cl), cl=cleft, &
                    cright), diagel
                Else
                  Write (infile(indnt+rlbwid+2:), Fmt=fmat) diagel
                End If
              Else
                Write (infile(indnt+rlbwid+2:), Fmt=fmat)(a(j,cl), cl=cleft, &
                  cright)
              End If
            Else
!             Upper triangular matrix.
              If (prdiag) Then
                Write (infile(indnt+rlbwid+2+numwid*(nd- &
                  1):), Fmt=fmat) diagel, (a(j,cl), cl=cleft, cright)
              Else
                If (cleft<=cright) Then
!                 Have to do the check on CLEFT and CRIGHT to
!                 avoid INDNT+RLBWID+2+NUMWID*ND possibly being
!                 out of range.
                  Write (infile(indnt+rlbwid+2+numwid*nd:), Fmt=fmat)(a(j,cl), &
                    cl=cleft, cright)
                End If
              End If
            End If
          End If
!
!         Output the (partial) matrix row.
          Call x04bafn(nout, infile(1:ncols))
!
        End Do
!
        clshft = clshft + nelems
        If (i/=nsets) Then
          Call x04bafn(nout, ' ')
        End If
      End Do
!
120   Continue
      errbuf(200:200) = char(0)
      ifail = p01abf(ifail, ierr, srname, nrec, rec)
!
      Return
!
130   Format (' ** On entry, MATRIX is not valid: MATRIX = ''', A, '''.')
140   Format (A, 3(1X,I5), 1X, A)
150   Format (' ** On entry, M .lt. LDA: M = ', I16, ', LDA = ', I16, '.')
160   Format (A, 3(1X,I5), 1X, I16, 1X, I16)
170   Format (' ** On entry, MATRIX = ''', A, ''', but DIAG is not val', &
        'id: DIAG = ''', A, '''.')
180   Format (A, 3(1X,I5), 1X, A, 1X, A)
190   Format (' ** On entry, either LABROW or LABCOL is not valid', /, 4X, &
        'LABROW = ''', A, ''', LABCOL = ''', A, '''.')
200   Format (A, 3(1X,I5), 1X, A, 1X, A)
210   Format (' ** On entry, FORM has more than 80 characters: the', &
        ' first 75 are', /, 4X, A, '.')
220   Format (A, 3(1X,I5), 1X, A)
230   Format (' ** The format specifier in FORM cannot be used to ', &
        'print a number: FORM =', /, 4X, A, '.')
240   Format (A, 3(1X,I5), 1X, A)
250   Format (' ** On entry, NCOLS-INDENT is not wide enough to hold', &
        ' at least one matrix', /, 4X, 'column: NCOLS = ', I16, ', INDE', &
        'NT =', I16, '.')
260   Format (A, 3(1X,I5), 1X, I16, 1X, I16)
    End Subroutine


text:
SUBROUTINE MULTIMODETIMEEVOLUTINOPERATOR(D,NM,MODES_NUM,U_F_MODES,E_MULTIFLOQUET,D_BARE,FIELD,T1,T2,U,INFO)

  ! TIME EVOLUTION OPERATOR OF A MULTIMODE DRESSED SYSTEM. THE EVOLUTION OPERATOR IS WRITEN IN THE BASIS USED TO EXPRESS THE
  ! MULTIMODE FLOQUET HAMILTONIAN
  ! U : MATRIX OF AMPLITUED OF PROBABILITIES FOR TRANSITIONS BETWEEN T1 TO T2
!!$  D              (IN)   : DIMENSION OF THE EXTENDED HILBERT SPACE (SIZE OF THE MULTIMODE FLOQUET MATRIX)
!!$  NM             (IN)   : NUMBER OF MODES
!!$  MODES_NUM      (IN)   : VECTOR (NM) INDICATING THE NUMBER OF HARMONICS OF EACH MODE
!!$  U_F_MODES      (IN)   : TRANSFORMATION, DIMENSOON (D,D)
!!$  E_MULTIFLOQUET (IN)   : MULTIMODE FLOQUET SPECTRUM
!!$  D_BARE         (IN)   : DIMENSION OF THE BARE HILBERT SPACE
!!$  FIELD          (IN)   : STRUCTURE DESCRIBING THE COUPLINGS
!!$  T1             (IN)   : INITIAL TIME
!!$  T2             (IN)   : FINAL TIME
!!$  U              (OUT)  : TRANFORMATION BETWEEN THE EXTENDED BARE BASIS AND THE FLOQUET STATES, DIMENSION (D_BARE,D)
!!$  INFO           (INOUT): (POSSIBLE) ERROR FLAG

  USE TYPES
  USE SUBINTERFACE_LAPACK

  IMPLICIT NONE
  INTEGER,                                    INTENT(IN)    :: D,D_BARE,NM ! DIMENSION OF THE MULTIMODE FLOQUET SPACE AND THE BARE BASIS
  INTEGER,                                    INTENT(INOUT) :: INFO
  INTEGER,          DIMENSION(NM),            INTENT(IN)    :: MODES_NUM
  TYPE(MODE),       DIMENSION(NM),            INTENT(IN)    :: FIELD  ! FIELDS PROPERTIES: FREQUENCY, AMPLITUDES AND PHASES
  DOUBLE PRECISION,                           INTENT(IN)    :: T1,T2  ! IN SECONDS
  DOUBLE PRECISION, DIMENSION(D),             INTENT(IN)    :: E_MULTIFLOQUET ! SET OF MULTIMODE FLOQUET ENERGIES, IN Hz, TO AVOID HBAR FACTORS
  COMPLEX*16,       DIMENSION(D,D),           INTENT(IN)    :: U_F_MODES   ! TRANFORMATION MATRIX BETWEEN DRESSED AND BARE BASIS
  COMPLEX*16,       DIMENSION(D_BARE,D_BARE), INTENT(OUT)   :: U           ! EVOLUTION OPERATOR U(T2,T1)

  COMPLEX*16, DIMENSION(D,D)      :: U_DIAGONAL
  COMPLEX*16, DIMENSION(D_BARE,D) :: U_AUX

  INTEGER :: MULTIMODE_HARMONICS, n,i,j,m,index0,index1,FIELD_INDEX
  DOUBLE PRECISION, DIMENSION(NM) :: OMEGA_VEC
  INTEGER, DIMENSION(NM) :: N_FLOQUET
  TYPE(HARMONIC_FACTORS), DIMENSION(:),ALLOCATABLE:: U_MODES_n
  DOUBLE PRECISION :: PHASE

  !write(*,*) size(U_F_MODES)
  MULTIMODE_HARMONICS   = D/D_BARE

  N_FLOQUET = 0
  OMEGA_VEC = 0

  DO n=2,NM
     FIELD_INDEX  = 2+SUM(MODES_NUM(2:n-1))
     N_FLOQUET(n) = FIELD(FIELD_INDEX)%N_Floquet
     OMEGA_VEC(n) = FIELD(FIELD_INDEX)%OMEGA
     !write(*,*) n,N_FLOQUET(n),field_index,modes_num(n)
     IF(modes_num(n).GT.N_FLOQUET(n)+1) THEN
        WRITE(*,*) "TO BUILD THE EXTENDED HAMILTONIAN THE NUMBER OF FLOQUET MODES MUST BE DEFINED"
        WRITE(*,*) "LARGER THAN THE NUMBER OF FIELD MODES"
        INFO = -10
     END IF
  END DO

  ALLOCATE(U_MODES_n(MULTIMODE_HARMONICS))

  DO n=1,MULTIMODE_HARMONICS
     ALLOCATE(U_MODES_n(n)%U(D_BARE,D))
     U_MODES_n(n)%U = 0.0
     ALLOCATE(U_MODES_n(n)%n(NM))
     U_MODES_n(n)%n = 0.0
  END DO

  DO i=0,MULTIMODE_HARMONICS-1
     U_MODES_n(i+1)%n = 0
     index0 = i
     DO j=2,NM
        U_MODES_n(i+1)%n(j)= -N_FLOQUET(j) + MOD(index0,(2*N_FLOQUET(j)+1))
        index0 = INT(index0/(2*N_FLOQUET(j)+1))
     END DO
  END DO

  i  = 1
  DO index1=1,MULTIMODE_HARMONICS
     DO n=1,D
        j = n
        i = 1 + (index1-1)*D_BARE
        DO m=1,D_BARE
           U_MODES_n(index1)%U(m,n) = U_F_MODES(i,j)
           i = i + 1
        END DO
     END DO
!     write(*,*) u_modes_n(index1)%U
  END DO

  index0 = (MULTIMODE_HARMONICS -1)/2 + 1
  U_AUX  =  U_MODES_n(index0)%U

  U_DIAGONAL = 0.0
  DO i=1,D
     U_DIAGONAL(i,i) = EXP(-DCMPLX(0.0,1.0)*E_MULTIFLOQUET(i)*(T2-T1))
!     write(*,*) D,i,U_DIAGONAL(i,i),E_MULTIFLOQUET(i)!*(T2-T1),t2-t1
  END DO

  DO index1=1,MULTIMODE_HARMONICS
     U_MODES_n(index1)%U = MATMUL(U_MODES_n(index1)%U,U_DIAGONAL)
  END DO

  U = 0.0
  DO index1=1,MULTIMODE_HARMONICS
     PHASE = DOT_PRODUCT(U_MODES_n(index1)%n,omega_vec)*T2
     U = U +  MATMUL(U_MODES_n(index1)%U,TRANSPOSE(CONJG(U_AUX)))*EXP(DCMPLX(0.0,1.0)*PHASE)
  END DO
  !write(*,*) abs(U)**2
END SUBROUTINE MULTIMODETIMEEVOLUTINOPERATOR

SUBROUTINE MULTIMODETIMEEVOLUTINOPERATOR_RESTRICTED(D,NF,U_F_MODES,E_MULTIFLOQUET,D_BARE,FIELD,T1,T2,U,INFO)

  ! TIME EVOLUTION OPERATOR OF A MULTIMODE DRESSED SYSTEM. THE EVOLUTION OPERATOR IS WRITEN IN THE BASIS USED TO EXPRESS THE
  ! MULTIMODE FLOQUET HAMILTONIAN
  ! U : MATRIX OF AMPLITUED OF PROBABILITIES FOR TRANSITIONS BETWEEN T1 TO T2
  !
  ! BUT THE SUM OVER FLOQUET EIGENVALUES IS RESTRICTED TO THOSE ACCURATELY DEFINED, FAR FROM THE EDGES OF THE MATRICES
  ! THE SUBROUTINE HAS BEEN DESIGNED FOR 87Rb driven by RF+MW fields, AND IT IS NOT EXPECTED TO WORK FOR OTHER CONFIGURATIONS
  ! D: THE DIMENSOIN OF THE MULTIMODE FLOQUET IS NOW DEFINED AS THE NUMBER OF FLOQUET MODES USED TO CALCULATE THE EVOLUTION OPERATOR
  ! the point of this is to make a comparison with MULTIMODETIMEEVOLUTIONOPERATOR AND CHECK WHETHER.
  ! MULTIMODETIMEEVOLUTIONOPERATOR SHOULD WORK FOR ANY ATOM/FIELD CONFIGURATION.
  ! MULTIMODETIMEEVOLUTIONOPEATOR_RESTRICTED IS DIFFICULT TO RESTRICT IN GENERAL SINCE WE NEED TO IDENTIFY THE "CENTRAL" MANIFOLD OF FLOQUET
  ! STATES

!!$  D              (IN)   : DIMENSION OF THE EXTENDED HILBERT SPACE (SIZE OF THE MULTIMODE FLOQUET MATRIX)
!!$  NM             (IN)   : NUMBER OF MODES
!!$  MODES_NUM      (IN)   : VECTOR (NM) INDICATING THE NUMBER OF HARMONICS OF EACH MODE
!!$  U_F_MODES      (IN)   : TRANSFORMATION, DIMENSOON (D,D)
!!$  E_MULTIFLOQUET (IN)   : MULTIMODE FLOQUET SPECTRUM
!!$  D_BARE         (IN)   : DIMENSION OF THE BARE HILBERT SPACE
!!$  FIELD          (IN)   : STRUCTURE DESCRIBING THE COUPLINGS
!!$  T1             (IN)   : INITIAL TIME
!!$  T2             (IN)   : FINAL TIME
!!$  U              (OUT)  : TRANFORMATION BETWEEN THE EXTENDED BARE BASIS AND THE FLOQUET STATES, DIMENSION (D_BARE,D)
!!$  INFO           (INOUT): (POSSIBLE) ERROR FLAG

  USE TYPES
  USE SUBINTERFACE_LAPACK

  IMPLICIT NONE
  INTEGER,                                    INTENT(IN)    :: D,D_BARE,NF ! DIMENSION OF THE MULTIMODE FLOQUET SPACE AND THE BARE BASIS
  INTEGER,                                    INTENT(INOUT) :: INFO
  TYPE(MODE),       DIMENSION(NF),     INTENT(IN)    :: FIELD  ! FIELDS PROPERTIES: FREQUENCY, AMPLITUDES AND PHASES
  DOUBLE PRECISION,                           INTENT(IN)    :: T1,T2  ! IN SECONDS
  DOUBLE PRECISION, DIMENSION(D),             INTENT(IN)    :: E_MULTIFLOQUET ! SET OF MULTIMODE FLOQUET ENERGIES, IN Hz, TO AVOID HBAR FACTORS
  COMPLEX*16,       DIMENSION(D,D),           INTENT(IN)    :: U_F_MODES   ! TRANFORMATION MATRIX BETWEEN DRESSED AND BARE BASIS
  COMPLEX*16,       DIMENSION(D_BARE,D_BARE), INTENT(OUT)   :: U           ! EVOLUTION OPERATOR U(T2,T1)

  COMPLEX*16, DIMENSION(:,:),ALLOCATABLE :: U_DIAGONAL
  COMPLEX*16, DIMENSION(:,:),ALLOCATABLE :: U_AUX

  INTEGER :: MULTIMODE_HARMONICS, n,i,j,m,index0,index1,D_r,i_0,i_r,Fup,Fdown
  DOUBLE PRECISION, DIMENSION(NF) :: OMEGA_VEC
  TYPE(HARMONIC_FACTORS), DIMENSION(:),ALLOCATABLE:: U_MODES_n
  DOUBLE PRECISION :: PHASE
  DOUBLE PRECISION :: A,HBAR

  MULTIMODE_HARMONICS   = D/D_BARE
  IF(NF.GT.2) THEN
     D_r                   = (2*FIELD(3)%N_FLOQUET - 2)*(2*(FIELD(2)%N_FLOQUET-2)-1)*D_bare
  END IF
  Fup                   = 2
  Fdown                 = 1
  hbar                  = 1.054E-34
  A                     = 2.0*4.0*ATAN(1.0)*hbar*3.417341E9

  DO i=1,NF
     OMEGA_VEC(i) = FIELD(i)%OMEGA
  END DO

  ALLOCATE(U_MODES_n(MULTIMODE_HARMONICS))

  DO n=1,MULTIMODE_HARMONICS
     ALLOCATE(U_MODES_n(n)%U(D_BARE,D))
     U_MODES_n(n)%U = 0.0
     ALLOCATE(U_MODES_n(n)%U_r(D_BARE,D_r))
     U_MODES_n(n)%U_r = 0.0
     ALLOCATE(U_MODES_n(n)%n(NF))
     U_MODES_n(n)%n = 0.0
  END DO

  DO i=0,MULTIMODE_HARMONICS-1
     U_MODES_n(i+1)%n = 0
     index0 = i
     DO j=2,NF
        U_MODES_n(i+1)%n(j)= -FIELD(j)%N_FLOQUET + MOD(index0,(2*FIELD(j)%N_FLOQUET+1))
        index0 = INT(index0/(2*FIELD(j)%N_FLOQUET+1))
     END DO
  END DO

  i  = 1
  DO index1=1,MULTIMODE_HARMONICS
     DO n=1,D
        j = n
        i = 1 + (index1-1)*D_BARE
        DO m=1,D_BARE
           U_MODES_n(index1)%U(m,n) = U_F_MODES(i,j)
           i = i + 1
        END DO
     END DO
  END DO

  index0 = (MULTIMODE_HARMONICS -1)/2 + 1

  !HOW MANY FLOQUET MODES AND WHICH ONES ARE WE GOING TO USE?
  !
  ! D_r  = (2*FIELD(3)%N_FLOQUET - 2)*(2*(FIELD(2)%N_FLOQUET-2)-1)*RB87%D_bare

  ! (2*FIELD(3)%N_FLOQUET - 2)       : number of mw manifolds, we neglect one up and one down
  ! (2*(FIELD(2)%N_FLOQUET - 2) - 1) : number of complete (i.e. having Rb87%D_bare) RF manifolds (not in general, it depends on Bdc and omega_rf)
  ! Rb87%D_bare                : number of bare states

  !WITHIN EACH MW MANIFOLD, WE KEEP FLOQUETS MODES WITH INDICES
  ! I \IN  I_0+1, I_0 +  (2*(FIELD(2)%N_FLOQUET - 2) - 1)*Rb87%D_bare
  ! WITH
  ! I_0 = 3*(2*Fdown+2*Fup+2)+1
  ! and the starting index for each manifold is
  ! (2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+1)+(m)*(2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+2*Fup+2)
  ! with m the manifold number, but we keep manifolds m=1,(2*FIELD(3)%N_FLOQUET - 2)

  ALLOCATE(U_DIAGONAL(D_r,D_r))
  ALLOCATE(U_AUX(D_bare,D_r))

  n  = 1
  U_DIAGONAL = 0.0
  DO i=1,(2*FIELD(3)%N_FLOQUET - 2)
     I_0 = (2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+1) + i*(2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+2*Fup+2)
     DO j=1,(2*(FIELD(2)%N_FLOQUET-2)-1)*D_bare

        i_r = I_0 + 3*(2*Fdown+2*Fup+2)+1 + j
        U_DIAGONAL(n,n) = EXP(-DCMPLX(0.0,1.0)*E_MULTIFLOQUET(i_r)*(T2-T1))
        n = n + 1
     END DO
  END DO

  DO index1=1,MULTIMODE_HARMONICS
     n = 1
     DO i=1,(2*FIELD(3)%N_FLOQUET - 2)
        I_0 = (2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+1) + i*(2*FIELD(2)%N_FLOQUET+1)*(2*Fdown+2*Fup+2)
        DO j=1,(2*(FIELD(2)%N_FLOQUET-2)-1)*D_bare
           i_r = I_0 + 3*(2*Fdown+2*Fup+2)+1 + j
           U_MODES_n(index1)%U_r(:,n) = U_MODES_n(index1)%U(:,i_r)
           n = n+1
        END DO
     END DO
  END DO

  index0 = (MULTIMODE_HARMONICS -1)/2 + 1
  U_AUX  =  U_MODES_n(index0)%U_r

  DO index1=1,MULTIMODE_HARMONICS
     U_MODES_n(index1)%U_r = MATMUL(U_MODES_n(index1)%U_r,U_DIAGONAL)
  END DO

  U = 0.0
  DO index1=1,MULTIMODE_HARMONICS
     PHASE = DOT_PRODUCT(U_MODES_n(index1)%n,omega_vec)*T2
     U = U +  MATMUL(U_MODES_n(index1)%U_r,TRANSPOSE(CONJG(U_AUX)))*EXP(DCMPLX(0.0,1.0)*PHASE)
  END DO

END SUBROUTINE MULTIMODETIMEEVOLUTINOPERATOR_RESTRICTED

max_stars_repo_name: gsinuco/OPENMMF

max_issues_repo_name: gsinuco/OPENMMF

max_forks_repo_name: gsinuco/OPENMMF


program newtonmd

implicit none
real:: error
integer::istep,n,i,z,pointer,j, steps, q, points,outputfile,num, kappa
real :: mu1, mu2, stepsize, mu
real:: func
real, dimension (:), allocatable :: w0, muarray
real, dimension (:,:), allocatable :: jarray,c,fx,x0, x1,xf, bin0
character (LEN = 60) :: inputval, filename1, filename2
!CHARACTER*29 :: path ='C:\Users\17029\Desktop\gfortran\datfiles\'

error=1.0e-06
inputval = 'S'
kappa = 0

! initialising n
print*,"Enter the dimension"
READ(*,*)  n

allocate ( x0(n,1) )
allocate ( bin0(n,1) )
allocate ( w0(n) )
allocate ( fx(n,1))
allocate ( x1(n,1))
allocate ( xf(n,1))
allocate ( jarray(n,n) )
allocate ( c(n,n))

print*,"Enter mu1 value"
READ(*,*)  mu1

print*,"Enter mu2 value"
READ(*,*)  mu2

print*,"Enter points value"
READ(*,*)  points

!print*,"Enter the Guess Vector"
!DO i = 1, n
!   READ(*,*)  x0(i,1)
!END DO

print*,"Enter the aeta Vector"
DO i = 1, n
   READ(*,*)  w0(i)
END DO

!print*,"Enter kappa value, Enter -1 for all"
!READ(*,*)  kappa

stepsize = (mu2-mu1)/(points+1)
allocate ( muarray(points + 2))

muarray(1) = mu1
do i = 2,(points+1)
muarray(i) = muarray(i-1) + stepsize
end do
muarray(points+2) = mu2
outputfile = 1

do num = 0, (2**(n) -1)

write(filename1, '("Output",I0,".dat")') outputfile
filename2=  filename1

print*, "this work 2"
call dectobin_gen(n,bin0,num)

print*, 'this works'
print*, bin0

!if ( sum(bin0) /=  kappa) CYCLE

x0 = bin0
print*, "this works 3"
open(90, file=filename2,action="write",POSITION="APPEND")

q = 1
do while ( q /= (points + 3))

mu = muarray(q)

pointer= 1

xf = x0

!do while(   ((abs(NORM2(xf)-NORM2(x1))/abs(NORM2(xf))) > 0.00001) )

!do while(   (pointer /= 70) .or. ((abs(NORM2(xf)-NORM2(x0))/abs(NORM2(x0))) == 0 )   )

do while( pointer /= 70 )

z = 1

print*, x0

!! generating jacobian

call Jac(x0, mu, n, w0, jarray)

do i = 1,n
do j = 1,n
   print*, i,j,jarray(i,j)
end do
end do

!! generating it's inverse
print*, "generating inverse"

call inverse(jarray,c,n)

print*,'c is', c

!! generating fx(x0)

do while (z /= n+1  )
call eqn(x0, w0, z, n, mu, func)
fx(z,1) = func
z = z+1
end do

print*, fx

!! computing x1

x1 = x0 - matmul(c,fx)

!!do i = 1,n
!!do j = 1,n
!!   print*, i,j,c(i,j)
!!end do
!!end do

!!print*,'fx is', fx

!!print*,'c*fx is', matmul(c, fx)

!!print*, size(c), size(fx)

!!print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
!!print*, x1
!!print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

x0 = x1
pointer = pointer + 1
end do

   write(90, '(F8.5,X)', advance='no') muarray(q)

   do j=1,n
      write(90, '(F8.5,X)', advance='no') x1(j,1)
   end do

   write(90, *) ''  ! this gives you the line break

q = q+1

end do

outputfile = outputfile + 1
close(90)
end do

end program

!!!!!!!!!! SUBROUTINES !!!!!!!!!!!!!

subroutine eqn(x0, w0, i, n, mu, f)

integer, intent(in):: n
real :: mu
real, DIMENSION(1:n) :: x0
real, DIMENSION(1:n) :: w0
real, intent(out):: f
integer:: i
real :: sumn

call sum(x0, w0, i, n, sumn)

print*,'sum = ', sumn
print*, 'mu =', mu

f = ((x0(i))**2) - (x0(i)) - (mu*sumn)

return
END

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

subroutine sum(x0, w0, i, n, sumn)

integer, intent(in):: n
real, DIMENSION(n,1) :: x0
real, DIMENSION(1:n) :: w0
real, intent(out):: sumn
integer:: i, p,k

k = i

p = k

sumn = 0

do while (k /= n)
sumn = sumn + (x0(p,1) - x0((k + 1),1))/(w0(p) - w0((k+1)) )
k = k + 1
end do

k = p

do while (k /= 1)
sumn = sumn + (x0(p,1) - x0((k - 1),1))/(w0(p) - w0((k-1)) )
k = k -1
end do

return
end

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

subroutine Jac(x0, mu, n, w0, jind)

integer :: n, j , i, k
real, DIMENSION(n,1) :: x0
real, DIMENSION(1:n) :: w0
real, DIMENSION(1:n, 1:n) :: jind
real :: mu, sumnaeta

print*, 'n value is ', n

j = 1

i = 1

do while( i /= n+1 )

do j = 1,n

print*, i, j

if (i == j) then
k = i
call sumaeta(w0, k, n, sumnaeta)
print*, 'once'
jind(i,j) =  2*(x0(i,1)) - 1 - (mu*sumnaeta)
end if

if(i /= j) then
print*, 'twice'
jind(i,j) =  mu/(w0(i)-w0(j))
end if

end do

i = i +1
end do

end

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

subroutine dectobin_gen(n,bin0,num)

integer :: n,num,i,temp
real, DIMENSION(n,1) :: bin0

do i = 1,n

bin0(i,1) = 0

end do

temp= num

do i=1,n

if (mod(temp,2)==1) bin0(i,1)=1

temp   =   temp/2

if (temp==0) exit

end do

end

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

subroutine sumaeta( w0, i, n, sumnaeta)

integer :: n
real, DIMENSION(1:n) :: w0
real :: sumnaeta, sumn
integer:: i, p

p = i

sumn = 0

do while (i /= n)
sumn = sumn + (1/(w0(p)-w0(i+1)))
i = i + 1
end do

i = p

do while (i /= 1)
sumn = sumn + (1/(w0(p) - w0(i-1)))
i = i -1
end do

sumnaeta = sumn

return
end

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

 subroutine inverse(a,c,n)
!============================================================
! Inverse matrix
! Method: Based on Doolittle LU factorization for Ax=b
! Alex G. December 2009
!-----------------------------------------------------------
! input ...
! a(n,n) - array of coefficients for matrix A
! n      - dimension
! output ...
! c(n,n) - inverse matrix of A
! comments ...
! the original matrix a(n,n) will be destroyed
! during the calculation
!===========================================================
implicit none
integer n
real :: a(n,n), c(n,n)
real :: L(n,n), U(n,n), b(n), d(n), x(n)
real :: coeff
integer i, j, k

! step 0: initialization for matrices L and U and b
! Fortran 90/95 aloows such operations on matrices
L=0.0
U=0.0
b=0.0

! step 1: forward elimination
do k=1, n-1
   do i=k+1,n
      coeff=a(i,k)/a(k,k)
      L(i,k) = coeff
      do j=k+1,n
         a(i,j) = a(i,j)-coeff*a(k,j)
      end do
   end do
end do

! Step 2: prepare L and U matrices
! L matrix is a matrix of the elimination coefficient
! + the diagonal elements are 1.0
do i=1,n
  L(i,i) = 1.0
end do
! U matrix is the upper triangular part of A
do j=1,n
  do i=1,j
    U(i,j) = a(i,j)
  end do
end do

! Step 3: compute columns of the inverse matrix C
do k=1,n
  b(k)=1.0
  d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution
  do i=2,n
    d(i)=b(i)
    do j=1,i-1
      d(i) = d(i) - L(i,j)*d(j)
    end do
  end do
! Step 3b: Solve Ux=d using the back substitution
  x(n)=d(n)/U(n,n)
  do i = n-1,1,-1
    x(i) = d(i)
    do j=n,i+1,-1
      x(i)=x(i)-U(i,j)*x(j)
    end do
    x(i) = x(i)/u(i,i)
  end do
! Step 3c: fill the solutions x(n) into column k of C
  do i=1,n
    c(i,k) = x(i)
  end do
  b(k)=0.0
end do
end subroutine inverse

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


text:
The Yolo Bypass Wildlife Area is located within the Yolo Bypass in Yolo County, California. The Wildlife Area is managed by the California Department of Fish and Wildlife with the intent of restoring and managing a variety of wildlife habitats in the Yolo Basin, a natural basin in the north part of the SacramentoSan Joaquin River Delta. The creation of the Wildlife Area was spearheaded by the Yolo Basin Foundation. The California Department of Fish and Wildlife and Yolo Basin Foundation are the core partners in the operation of this unique community resource.

History

The Wildlife Area was open to the public in 1997 after extensive restoration efforts completed by Ducks Unlimited with federal funds appropriated through the United States Army Corps of Engineers. In 1999, this 3,700acre (15 km²) restoration project was named the Vic Fazio Yolo Wildlife Area in honor of congressman Vic Fazio who lobbied hard for the funds needed to build the project. In 2001 the Wildlife Area expanded to over 16,000 acres (65 km²) with the acquisition of the Glide and Los Rios properties. These acquisitions included the 10,000acre (40 km²) Tule Ranch, a working cattle ranch with extensive vernal pool areas.

Flood Control

The Yolo Bypass is a 59,000acre (240 km²) flood control channel that protects Sacramento and other cities from flooding. The Wildlife Area was created with the understanding that it would remain completely compatible with this primary flood control function. For this reason, there are restrictions on the density of emergent vegetation and riparian trees within the Wildlife Area. These standards are determined through the use of hydrologic models.

Located at the north end of the Yolo Basin where Putah Creek enters the Yolo Bypass, this part of the Delta is known as the Putah Sinks and hosts a diverse assemblage of wildlife species inhabiting seasonal wetlands, permanent wetlands, riparian forest, uplands, vernal pools and agricultural habitats. This area remains as a critical link in the complex chain of
habitats that make up the Central Valley wetlands, a critical part of the Pacific Flyway.

Wildlife Habitat Types

Seasonal Wetlands
Most wetlands on the Wildlife Area are managed as seasonal wetlands. They go through an extensive dry period during the spring and summer months. Typically these ponds are drained April 1st to stimulate the germination of Swamp Timothy. They may receive a brief summer irrigation and then are flooded in September to provide wetland habitat for migratory waterfowl and shorebirds. Because of its low stature, Swamp Timothy provides little resistance to flood flows, thus preserving the primary flood protection purpose of the Yolo Bypass. Additionally, this plant produces small seeds for waterfowl and a suitable substrate for the production of invertebrates, an important food source for shorebirds.

Permanent Wetlands

Permanent wetlands are flooded year round and tend to be deeper and have more emergent vegetation. These ponds provide important brood water for resident aquatic birds such as mallards, common moorhens and piedbilled grebes. They also provide drinking water for wildlife as well as relief from intense summer heat. This habitat is used extensively for roosting by blackcrowned night herons, egrets, and whitefaced ibis. The Giant Garter Snake also occurs in this habitat.

Uplands

Upland habitats are used extensively by ground nesting birds such Northern Harriers, Western Meadowlarks, Mallards and Ringnecked Pheasant. Vegetation is typically dominated by annual rye grass, curly dock and wild sunflower. Typically, legumes such as vetch or various clovers are planted to increase the value of this habitat for nesting birds. Late spring irrigations through parts of the uplands can substantially increase the chick survival of Ringnecked Pheasant while increasing the occurrence of passerines such as Blue Grosbeak and Western Meadowlark.  Rodent populations in these areas provide prey for large numbers of wintering birds of prey. As rodent populations increase, so do the number of nesting raptors including Northern Harriers and Shorteared Owls.

Riparian Forest

Riparian vegetation consists of willows, cottonwoods, black walnut and other tree species. Nesting species include Swainsons hawks, greathorned owls, wood ducks, tree swallows, and black phoebe. This habitat is very important for neotropical migrants such as blue grosbeak, ashthroated flycatchers, and a variety of warbler species. Most riparian vegetation is located along Putah Creek in the central part of the Wildlife Area.

Vernal Pools
Approximately 1,800 acres (7.3 km²) of natural uplands occur in the southwest portion of the Wildlife Area. There are several vernal pools in this area that are home to such invertebrate species as tadpole shrimp, clam shrimp and the endangered Conservancy fairy shrimp. Rare and endangered plants include Heckerts pepperweed and Ferris alkali milk vetch. Grassland bird species in this area include grasshopper sparrow, savannah sparrow and burrowing owl.

Agriculture
An interesting feature of the Wildlife Area is the extensive use of agriculture to achieve its wildlife habitat goals, while providing important operating income. Cattle are used to control unwanted vegetation and are the primary management tool in the vernal pool habitat area. Their removal of thatch from the thick stands of annual rye grass help facilitate the germination of native forbs in this area, resulting in spectacular blooms of wildflowers each spring. In the northern portions of the Wildlife Area, rice is grown, which is then flooded after harvest, attracting thousands of waterfowl and shorebirds in full view of the thousands of automobiles on Interstate 80. Wildlife Managers have instituted a unique rotation within the rice fields that allows for a fallow stage every three years. During this fallow stage, the field is managed for migratory shorebirds by flooding during the mid summer months.

Public Use

Image(auto tour map.jpg, 300, left, thumbnail, Auto Tour Map.)Image(tour route.jpg, 300, right, thumbnail, Large numbers of waterfowl on the auto tour route.)The Wildlife Area is open every day (except Christmas) for wildlife viewing and fishing. October through June, there are monthly tours on the second Saturday of each month conducted by volunteers from the Yolo Basin Foundation. Other public use opportunities are provided by the Yolo Basin Foundation and Fish and Wildlife including bat tours, open houses, and a fall/winter speaker series. Information regarding these activities can be found on the Yolo Basin Foundation events page http://www.yolobasin.org/events.cfm
A six mile auto tour loop is open every day. This road travels through seasonal wetlands, permanent wetlands, uplands and agricultural areas. Especially during winter months, this route offers excellent wildlife viewing opportunities. Please park in turn outs and marked parking lots to enjoy hiking through this extensive wetland complex.

Birding
The Yolo Bypass Wildlife Area has become one of the premier birding areas in Northern California. Large numbers of waterfowl can be viewed in the winter months in both wetlands and flooded agricultural areas. They are often in the vicinity of shorebirds, many of which have traveled from the Arctic Coast to spend the winter in the Central Valley. The winter bird population also includes a wide variety of raptors. These birds of prey include Redtailed Hawks, Roughlegged Hawks, Peregrine Falcon, Praire Falcon, Merlin, Coopers Hawk, Sharpshinned Hawk and Whitetailed Kite.  Bird numbers peak from December to March, when the northbound migration is in full swing. Many ducks are moving north in March, with shorebird migration peaking in mid April. At this time, birds from subtropical areas begin moving north to spend the summer in our area. Large numbers of warblers, tanagers, flycatchers and other neotropical migrants can be found in the few areas with willow riparian forests.

Resident birds begin nesting in the late Winter and Spring. Ground nesting birds such as Ringnecked Pheasants, Mallards and Northern Harriers use the mix of rank upland growth from previous years and the fresh sprouting grasses and forbs  for nesting cover. Blacknecked Stilts and American Avocets begin nesting on the receding mud flats of the seasonal wetlands in April. The young of these species are soon observed teetering on fresh legs as they probe the mud for their invertebrate prey.

By the last week in June, shorebirds begin arriving from the far north. Specially managed fallow rice fields are flooded in anticipation of their arrival and large numbers of these birds are soon feeding and roosting by the thousands. By late July, shorebirds are joined by Mallards and Gadwall who have nested in other parts of California, and gravitated to the few remaining wetlands available to them. Soon, the Pintail and Greenwinged Teal arrive and the wintering flocks once again darken the skies over the Yolo Bypass.

Environmental EducationImage(pond study.JPG, right, thumbnail, School children learning about pond critters.)
The successful Discover the Flyway environmental education program annually brings approximately 4,000 students per year to learn about wetlands and visit the Wildlife Area. This program served school children from throughout the Sacramento area, serving over 100 schools from 15 different school districts.

The program utilizes a 3 acre wetland demonstration area at the Wildlife Area headquarters as an outdoor classroom. Educational activities are set up to learn about animal tracks, aquatic invertebrates, plant identification as well as social studies lessons in Native American plant use and several other activities. The  students then go on a tour of the Wildlife Area to see first hand, our local contribution to the Pacific Flyway chain of wetland areas.

Volunteering
There are opportunities for volunteers to contribute to the environmental education and public use programs. Docents for public tours are needed as well. There is great need for volunteers to help with special events such as Duck Days, the Tule Ranch Open House and other activities. Additionally, there are special field trips and training activities available for volunteers. For more information and to learn about volunteer opportunities at the Wildlife Area contact the Yolo Basin Foundation 3.Image(safety talk.jpg, right, thumbnail, Safety orientation given by game warden prior to the junior hunt.)

Hunting
The Department of Fish and Wildlife operates public hunting programs for waterfowl, pheasant and mourning dove  in specific areas of the Wildlife Area during the fall and winter months. Dove hunting season begins on September 1st and lasts fifteen days. Waterfowl season begins in mid to late October and usually runs through the last weekend in January. Pheasant season traditionally starts on the second Saturday in November and runs until mid to late December. Please consult the appropriate California hunting regulations 2 for more specific information. Junior hunting programs are conducted an an annual basis at the Wildlife Area. The junior pheasant hunt takes place on the second Saturday in November. Fifty birds are planted for 25 lucky hunters. The program includes safety discussions with Fish and Wildlife Game Wardens. The weekend following the end of duck season, the Wildlife Area is open for junior waterfowl hunting.

More information can be found at the Yolo Basin Foundation http://www.yolobasin.org/
and Department of Fish and Game Wildlife Area website http://www.dfg.ca.gov/lands/wa/region3/yolo/index.html

A gallery of photographs by the former Area Manager (Dave Feliz), is available for viewing:http://picasaweb.google.com/yolodave/YoloBypassWildlifeArea

20080309 22:43:09 nbsp I went to the Wildlife Area today for the first time and drove the sixmile route, stopping a dozen times to get out and walk and just listen and watch birds. It was so relaxing and quiet once we got away from the freeway traffic. The only drawback was the one or two drivers that drove faster than the speed limit (15 mph), leaving behind clouds of dust that covered us as we walked or drove. Ill go again to do more extensive walking. Users/robinlaughlin

20090207 08:04:10 nbsp I went to the Wildlife Area yesterday to take some photos and a sign said it was affected by the biweekly furlough thats hitting all state agencies.  So it will be closed Feb 20th, etc. Users/RonPluth

20100106 18:07:40 nbsp Does anyone know about dog regulations here? Can you take your dogs walking with you in the Bypass? Users/CParker
  I dont remember for sure, but I believe theyre only allowed to the north of the first parking lot.  All the areas to the south are closed.  I suspect one of the signs in the area will talk about it though.   Users/WesHardaker
  According to the Yolo Basin website: At this time bicycles and dogs are allowed only in the Causeway which is between Interstate 80 and the Railroad tracks.

20130212 17:18:36 nbsp I spent some time here for the first time last Sunday.  Great trip.  I took my parents, who are avid bird watchers, out along the driving loop.  Lots of great places to stop and look around.  Incredible variety of birds.  There are a few picnic tables and restrooms scattered around, so on a nice day its easy to settle in and enjoy yourself. Users/TomGarberson

20130313 14:57:48 nbsp We just went on the monthly (second Saturday) YBWA tour, which was led by volunteer docents Rob Floerke and Ken Ealing  and it was WONDERFUL!  Ive been on many a beginning birding tour, but this was the best even.  Their information was clear, their stories were great.  No scriptreading for these men!  Whatever was seen or found was immediately explained clearly.  When that little girl came up with a duck skeleton (feathered wings still attached), I was amazed at the info that fountained out.  Id happily go on a tour with either Rob or Ken again.  Hmmmm... maybe we will!  See you at the next second Saturday?    Users/LoisRichter


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
SUBROUTINE simple_mixing(c1,c2)
use scf_data
implicit none

integer :: ix,iy,iz
real(8),intent(IN) :: c1,c2

if(ilsda == 0)then
!$OMP parallel do private(iz,iy,ix)
  do iz=ng_sta(3),ng_end(3)
  do iy=ng_sta(2),ng_end(2)
  do ix=ng_sta(1),ng_end(1)
    rho_out(ix,iy,iz,num_rho_stock)=rho(ix,iy,iz)
  end do
  end do
  end do
elseif( ilsda==1 )then
!$OMP parallel do private(iz,iy,ix)
  do iz=ng_sta(3),ng_end(3)
  do iy=ng_sta(2),ng_end(2)
  do ix=ng_sta(1),ng_end(1)
    rho_s_out(ix,iy,iz,num_rho_stock,1:2)=rho_s(ix,iy,iz,1:2)
  end do
  end do
  end do
end if

!rho = c1*rho + c2*matmul( psi**2, occ )
if(ilsda == 0)then
  do iz=ng_sta(3),ng_end(3)
  do iy=ng_sta(2),ng_end(2)
  do ix=ng_sta(1),ng_end(1)
    rho(ix,iy,iz) = c1*rho_in(ix,iy,iz,num_rho_stock) + c2*rho_out(ix,iy,iz,num_rho_stock)
    rho_in(ix,iy,iz,num_rho_stock+1) = rho(ix,iy,iz)
  end do
  end do
  end do
else if(ilsda == 1)then
  do iz=ng_sta(3),ng_end(3)
  do iy=ng_sta(2),ng_end(2)
  do ix=ng_sta(1),ng_end(1)
    rho_s(ix,iy,iz,1) = c1*rho_s_in(ix,iy,iz,num_rho_stock,1) + c2*rho_s_out(ix,iy,iz,num_rho_stock,1)
    rho_s(ix,iy,iz,2) = c1*rho_s_in(ix,iy,iz,num_rho_stock,2) + c2*rho_s_out(ix,iy,iz,num_rho_stock,2)
    rho(ix,iy,iz) = rho_s(ix,iy,iz,1)+rho_s(ix,iy,iz,2)
    rho_s_in(ix,iy,iz,num_rho_stock+1,1) = rho_s(ix,iy,iz,1)
    rho_s_in(ix,iy,iz,num_rho_stock+1,2) = rho_s(ix,iy,iz,2)
  end do
  end do
  end do
end if

return

END SUBROUTINE simple_mixing

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON