text:
      SUBROUTINE EWDAT_CPV(PAR,PROB)

      IMPLICIT NONE

      DOUBLE PRECISION PAR(*),PROB(*)
      DOUBLE PRECISION Pi,aux,NMB0,NMB22,NMH,CEW0,BEW1,FSF1,DEW0,DEW27
      DOUBLE PRECISION PiZMZ,PIZO,PiWMW,PiWO,Drho,rho2SM,dVB

      DOUBLE PRECISION GF,MZ,MW,g1,g2,alSMZ,S2TW,ALEMMZ
      DOUBLE PRECISION mt,mb,mtau,mmu,mel,MS,MC,MBP,MPI,MSTRANGE
      DOUBLE PRECISION tanb,cosb,sinb,vu,vd
      DOUBLE PRECISION MSL3,MSE3,MSL1,MSE1,ATAU,AMU
      DOUBLE PRECISION MSQ1,MSU1,MSD1
      DOUBLE PRECISION MSQ3,MSU3,MSD3,AT,AB
      DOUBLE PRECISION mur,M1r,M2r,msi

      COMMON/EWPAR/GF,MZ,MW,g1,g2,alSMZ,S2TW,ALEMMZ
      COMMON/SMFERM/mt,mb,mtau,mmu,mel,MS,MC,MBP,MPI,MSTRANGE
      COMMON/TBPAR/tanb,cosb,sinb,vu,vd
      COMMON/SLEPPAR/MSL3,MSE3,MSL1,MSE1,ATAU,AMU
      COMMON/SQUPAR/MSQ1,MSU1,MSD1
      COMMON/RADCOR2/MSQ3,MSU3,MSD3,AT,AB
      COMMON/GAUGINOPAR/mur,M1r,M2r,msi

      Pi=4d0*datan(1d0)

c      1) DRbar fine-structure constant

      aux=0.06887d0-1d0/(137.036d0*2d0*Pi)
     . *(-7d0*dlog(MW/MZ)+16d0/9d0*dlog(mt/MZ)+dlog(PAR(23)/MZ)/3d0
     .  +4d0/9d0*(2d0*dlog(MSQ1/MZ**2)+2d0*dlog(MSU1/MZ**2)
     .                    +dlog(MSQ3/MZ**2)+dlog(MSU3/MZ**2))/2d0
     .  +1d0/9d0*(2d0*dlog(MSQ1/MZ**2)+2d0*dlog(MSD1/MZ**2)
     .                    +dlog(MSQ3/MZ**2)+dlog(MSD3/MZ**2))/2d0
     .  +1d0/3d0*(2d0*dlog(MSL1/MZ**2)+2d0*dlog(MSE1/MZ**2)
     .                    +dlog(MSL3/MZ**2)+dlog(MSE3**2/MZ))/2d0
     .  +4d0/3d0*(dlog(M2r/MZ)+dlog(mur/MZ)))

      ALEMMZ=1d0/137.036d0/(1d0-aux)

c      2) EW gauge boson self-energies

      PiZMZ=MZ**2*NMB0(MZ**2,125d0**2,MZ**2,MZ**2)
     . -2d0*NMB22(MZ**2,125d0**2,MZ**2,MZ**2)
     . -2d0*NMB22(MZ**2,PAR(23)**2,PAR(23)**2,MZ**2)
     . -(1d0-2d0*S2TW)**2*(NMB22(MZ**2,MW**2,MW**2,MZ**2)
     .                 +NMB22(MZ**2,PAR(23)**2,PAR(23)**2,MZ**2))
     . -2d0*(1d0-S2TW)**2*(2d0*MZ**2+MW**2-MZ**2*S2TW**2/(1-S2TW))
     .                              *NMB0(MZ**2,MW**2,MW**2,MZ**2)
     . -8d0*(1d0-S2TW)**2*NMB22(MZ**2,MW**2,MW**2,MZ**2)
     . +3d0*((1d0/2d0-2d0/3d0*S2TW)**2+(2d0/3d0*S2TW)**2)
     .  *(2d0*NMH(MZ**2,0d0,0d0,MZ**2)+NMH(MZ**2,mt**2,mt**2,MZ**2))
     . +12d0*(2d0/3d0*S2TW)*(1d0/2d0-2d0/3d0*S2TW)*mt**2
     .                      *NMB0(MZ**2,mt**2,mt**2,MZ**2)
     . +3d0*((-1d0/2d0+S2TW/3d0)**2+(S2TW/3d0)**2)
     .  *(2d0*NMH(MZ**2,0d0,0d0,MZ**2)+NMH(MZ**2,mb**2,mb**2,MZ**2))
     . +12d0*(-1d0/3d0*S2TW)*(-1d0/2d0+1d0/3d0*S2TW)*mb**2
     .                      *NMB0(MZ**2,mb**2,mb**2,MZ**2)
     . +3d0*(1d0/2d0)**2*NMH(MZ**2,0d0,0d0,MZ**2)
     . +((-1d0/2d0+S2TW)**2+(S2TW)**2+(1d0/2d0)**2)
     .  *(2d0*NMH(MZ**2,0d0,0d0,MZ**2)+NMH(MZ**2,mtau**2,mtau**2,MZ**2))
     . +4d0*(-S2TW)*(-1d0/2d0+S2TW)*mtau**2
     .                      *NMB0(MZ**2,mtau**2,mtau**2,MZ**2)
     . -12d0*((1d0/2d0-2d0/3d0*S2TW)**2+(-1d0/2d0+S2TW/3d0)**2)
     .  *(2d0*NMB22(MZ**2,MSQ1,MSQ1,MZ**2)
     .                            +NMB22(MZ**2,MSQ3,MSQ3,MZ**2))
     . -12d0*(2d0/3d0*S2TW)**2
     .  *(2d0*NMB22(MZ**2,MSU1,MSU1,MZ**2)
     .                            +NMB22(MZ**2,MSU3,MSU3,MZ**2))
     . -12d0*(-1d0/3d0*S2TW)**2
     .  *(2d0*NMB22(MZ**2,MSD1,MSD1,MZ**2)
     .                            +NMB22(MZ**2,MSD3,MSD3,MZ**2))
     . -4d0*((1d0/2d0)**2+(-1d0/2d0+S2TW)**2)
     .  *(2d0*NMB22(MZ**2,MSL1,MSL1,MZ**2)
     .                            +NMB22(MZ**2,MSL3,MSL3,MZ**2))
     . -4d0*(-S2TW)**2
     .  *(2d0*NMB22(MZ**2,MSE1,MSE1,MZ**2)
     .                            +NMB22(MZ**2,MSE3,MSE3,MZ**2))
     . +(NMH(MZ**2,mur**2,mur**2,MZ**2)
     .                +2d0*mur**2*NMB0(MZ**2,mur**2,mur**2,MZ**2))/2d0
     . +2d0*(1d0-S2TW)**2*(NMH(MZ**2,M2r**2,M2r**2,MZ**2)
     .                +2d0*M2r**2*NMB0(MZ**2,M2r**2,M2r**2,MZ**2))
     . +(1d0-2d0*S2TW)**2*(NMH(MZ**2,mur**2,mur**2,MZ**2)
     .                +2d0*mur**2*NMB0(MZ**2,mur**2,mur**2,MZ**2))/2d0

      PiZMZ=(g1+g2)*PiZMZ/16d0/Pi**2

      PiZO=-MZ**2*FSF1(125d0**2,MZ**2,MZ**2)
     . -2d0*NMB22(0d0,125d0**2,MZ**2,MZ**2)
     . -2d0*NMB22(0d0,PAR(23)**2,PAR(23)**2,MZ**2)
     . -(1d0-2d0*S2TW)**2*(NMB22(0d0,MW**2,MW**2,MZ**2)
     .                 +NMB22(0d0,PAR(23)**2,PAR(23)**2,MZ**2))
     . +2d0*(1d0-S2TW)**2*(MW**2-MZ**2*S2TW**2/(1-S2TW))
     .                              *FSF1(MW**2,MW**2,MZ**2)
     . -8d0*(1d0-S2TW)**2*NMB22(0d0,MW**2,MW**2,MZ**2)
     . +3d0*((1d0/2d0-2d0/3d0*S2TW)**2+(2d0/3d0*S2TW)**2)
     .  *(2d0*0d0+NMH(0d0,mt**2,mt**2,MZ**2))
     . -12d0*(2d0/3d0*S2TW)*(1d0/2d0-2d0/3d0*S2TW)*mt**2
     .                      *FSF1(mt**2,mt**2,MZ**2)
     . +3d0*((-1d0/2d0+S2TW/3d0)**2+(S2TW/3d0)**2)
     .  *(2d0*0d0+NMH(0d0,mb**2,mb**2,MZ**2))
     . -12d0*(-1d0/3d0*S2TW)*(-1d0/2d0+1d0/3d0*S2TW)*mb**2
     .                      *FSF1(mb**2,mb**2,MZ**2)
     . +3d0*(1d0/2d0)**2*0d0
     . +((-1d0/2d0+S2TW)**2+(S2TW)**2+(1d0/2d0)**2)
     .  *(2d0*0d0+NMH(0d0,mtau**2,mtau**2,MZ**2))
     . -4d0*(-S2TW)*(-1d0/2d0+S2TW)*mtau**2
     .                      *FSF1(mtau**2,mtau**2,MZ**2)
     . -12d0*((1d0/2d0-2d0/3d0*S2TW)**2+(-1d0/2d0+S2TW/3d0)**2)
     .  *(2d0*NMB22(0d0,MSQ1,MSQ1,MZ**2)
     .                            +NMB22(0d0,MSQ3,MSQ3,MZ**2))
     . -12d0*(2d0/3d0*S2TW)**2
     .  *(2d0*NMB22(0d0,MSU1,MSU1,MZ**2)
     .                            +NMB22(0d0,MSU3,MSU3,MZ**2))
     . -12d0*(-1d0/3d0*S2TW)**2
     .  *(2d0*NMB22(0d0,MSD1,MSD1,MZ**2)
     .                            +NMB22(0d0,MSD3,MSD3,MZ**2))
     . -4d0*((1d0/2d0)**2+(-1d0/2d0+S2TW)**2)
     .  *(2d0*NMB22(0d0,MSL1,MSL1,MZ**2)
     .                            +NMB22(0d0,MSL3,MSL3,MZ**2))
     . -4d0*(-S2TW)**2
     .  *(2d0*NMB22(0d0,MSE1,MSE1,MZ**2)
     .                            +NMB22(0d0,MSE3,MSE3,MZ**2))
     . +(NMH(0d0,mur**2,mur**2,MZ**2)
     .                -2d0*mur**2*FSF1(mur**2,mur**2,MZ**2))/2d0
     . +2d0*(1d0-S2TW)**2*(NMH(0d0,M2r**2,M2r**2,MZ**2)
     .                -2d0*M2r**2*FSF1(M2r**2,M2r**2,MZ**2))
     . +(1d0-2d0*S2TW)**2*(NMH(0d0,mur**2,mur**2,MZ**2)
     .                -2d0*mur**2*FSF1(mur**2,mur**2,MZ**2))/2d0

      PiZO=(g1+g2)*PiZO/16d0/Pi**2

      PiWMW=MW**2*NMB0(MW**2,125d0**2,MW**2,MZ**2)
     . -NMB22(MW**2,125d0**2,MW**2,MZ**2)
     . -NMB22(MW**2,MZ**2,MW**2,MZ**2)
     . -2d0*NMB22(MW**2,PAR(23)**2,PAR(23)**2,MZ**2)
     . -((4d0*MW**2+MW**2+MZ**2)*(1d0-S2TW)-MZ**2*S2TW**2)
     .                        *NMB0(MW**2,MZ**2,MW**2,MZ**2)
     . -8d0*(1d0-S2TW)*NMB22(MW**2,MZ**2,MW**2,MZ**2)
     . -4d0*(2d0*NMB22(MW**2,MW**2,0d0,MZ**2)
     .                    +MW**2*NMB0(MW**2,MW**2,0d0,MZ**2))
     . +4d0*NMH(MW**2,0d0,0d0,MZ**2)
     . +3d0/2d0*NMH(MW**2,mt**2,mb**2,MZ**2)
     . +NMH(MW**2,mtau**2,0d0,MZ**2)/2d0
     . -6d0*(2d0*NMB22(MW**2,MSQ1,MSQ1,MZ**2)
     .             +NMB22(MW**2,MSQ3,MSQ3,MZ**2))
     . -2d0*(2d0*NMB22(MW**2,MSL1,MSL1,MZ**2)
     .             +NMB22(MW**2,MSL3,MSL3,MZ**2))
     . +2d0*(NMH(MW**2,M2r**2,M2r**2,MZ**2)
     .       +2d0*M2r**2*NMB0(MW**2,M2r**2,M2r**2,MZ**2))
     . +2d0*(NMH(MW**2,mur**2,mur**2,MZ**2)
     .       +mur**2*NMB0(MW**2,mur**2,mur**2,MZ**2))

      PiWMW=g2*PiWMW/16d0/Pi**2

      PiWO=-MW**2*FSF1(125d0**2,MW**2,MZ**2)
     . -NMB22(0d0,125d0**2,MW**2,MZ**2)
     . -NMB22(0d0,MZ**2,MW**2,MZ**2)
     . -2d0*NMB22(0d0,PAR(23)**2,PAR(23)**2,MZ**2)
     . +((MW**2+MZ**2)*(1d0-S2TW)-MZ**2*S2TW**2)
     .                        *FSF1(MZ**2,MW**2,MZ**2)
     . -8d0*(1d0-S2TW)*NMB22(0d0,MZ**2,MW**2,MZ**2)
     . -4d0*(2d0*NMB22(0d0,MW**2,0d0,MZ**2)
     .                    -MW**2*FSF1(MW**2,0d0,MZ**2))
     . +3d0/2d0*NMH(0d0,mt**2,mb**2,MZ**2)
     . +NMH(0d0,mtau**2,0d0,MZ**2)/2d0
     . -6d0*(2d0*NMB22(0d0,MSQ1,MSQ1,MZ**2)
     .             +NMB22(0d0,MSQ3,MSQ3,MZ**2))
     . -2d0*(2d0*NMB22(0d0,MSL1,MSL1,MZ**2)
     .             +NMB22(0d0,MSL3,MSL3,MZ**2))
     . +2d0*(NMH(0d0,M2r**2,M2r**2,MZ**2)
     .       -2d0*M2r**2*FSF1(M2r**2,M2r**2,MZ**2))
     . +2d0*(NMH(0d0,mur**2,mur**2,MZ**2)
     .       -mur**2*FSF1(mur**2,mur**2,MZ**2))

      PiWO=g2*PiWO/16d0/Pi**2

c      3) Custodial symmetry parameter

      Drho=(PiZMZ/MZ**2-PiWMW/MW**2)/(1d0+PiZMZ/MZ**2)
     . +ALEMMZ*ALSMZ/4d0/Pi**2/S2TW*(-2.145d0*(mt/MZ)**2
     .   +1.262d0*dlog(mt/MZ)-2.24d0-0.55*(MZ/mt))
     . +3*GF**2*mt**4/128d0/Pi**4*rho2SM(125d0/mt)

c      4) Vertex / Box contributions to the mu-nu_mu-e-nu_e eff. coupling

      dVB=ALEMMZ/4.d0/Pi/S2TW/(1d0-Drho)*(6d0+dlog(MW**2/MZ**2)
     . /(1d0-MW**2/MZ**2)*(1d0+5d0/2d0*MW**2/MZ**2
     . -(1d0-MW**2/MZ**2)*(5d0-3d0/2d0*MW**2/MZ**2/(1-S2TW))))

      aux=g1/2d0*BEW1(M1r**2,MSL1,MZ**2)
     . +3d0*g2/2d0*BEW1(M2r**2,MSL1,MZ**2)
     . +2d0*g2*M2r**2*CEW0(MSL1,M2r**2,M2r**2)
     . +g2*(FSF1(M2r**2,M2r**2,MZ**2)
     .          -MSL1*CEW0(MSL1,M2r**2,M2r**2)+1d0/2d0)
     . -g1/2d0/dsqrt(2d0)*(FSF1(MSL1,MSL1,MZ**2)
     .          -M1r**2*CEW0(M1r**2,MSL1,MSL1)-1d0/2d0)
     . +g2/2d0/dsqrt(2d0)*(FSF1(MSL1,MSL1,MZ**2)
     .          -M2r**2*CEW0(M2r**2,MSL1,MSL1)-1d0/2d0)

      dVB=dVB-2d0*aux/16d0/Pi**2

      aux=g1*g2/2d0*M1r*M2r*DEW0(MSL1,M1r**2,M2r**2)
     . -(g2*M2r)**2/2d0*DEW0(MSL1,M2r**2,M2r**2)
     . +g1*g2*DEW27(MSL1,M1r**2,M2r**2)
     . +g2**2*DEW27(MSL1,M2r**2,M2r**2)

      dVB=dVB-S2TW*(1d0-S2TW)/2d0/Pi/ALEMMZ*MZ**2/16d0/Pi**2*aux

c      5) EW vev's and couplings
c      print*,vu,vd,g1,g2
c      vu=vu/dsqrt(1d0-PiWO/MW**2-dVB)
c      vd=vd/dsqrt(1d0-PiWO/MW**2-dVB)

c      g2=2d0*MW**2/(vu**2+vd**2)*(1d0+PiWMW/MW**2)
c      g1=2d0*MZ**2/(vu**2+vd**2)*(1d0+PiZMZ/MZ**2)-g2
c      print*,vu,vd,g1,g2

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function NMB22(x,y,z,t)

c            ->NMB22(p^2,m1^2,m2^2,Q^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,t,aux,NMB0,NMA0,FSF1

      IF(x.ge.1d-4)THEN
      aux=y+z-x/3d0-(x-2d0*y-2d0*z)*NMB0(x,y,z,t)/2d0+NMA0(y,t)/2d0
     .      +NMA0(z,t)/2d0+(z-y)/2d0/x*(NMA0(z,t)-NMA0(y,t)
     .                                     -(z-y)*NMB0(x,y,z,t))
      ELSE
       IF(Max(y,z).ge.1d-4)THEN
      IF(min(y,z).ge.1d-4)THEN
      IF(dabs(z-y).ge.1d-4)THEN
      aux=y+z-(y+z)*FSF1(y,z,t)+NMA0(y,t)/2d0
     .      +NMA0(z,t)/2d0+y*z*dlog(z/y)/2d0/(z-y)-(y+z)/4d0
      ELSE
      aux=y+z-(y+z)*FSF1(y,z,t)+NMA0(y,t)/2d0+NMA0(z,t)/2d0
      ENDIF
      ELSE
       aux=-Max(y,z)/4d0
      ENDIF
       ELSE
      aux=0d0
       ENDIF
      ENDIF

      NMB22=aux/6d0-(NMA0(y,t)+NMA0(z,t))/4d0

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function NMH(x,y,z,t)

c            ->NMH(p^2,m1^2,m2^2,Q^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,t,aux,NMB0,NMA0,NMB22,FSF1

      IF(x.ge.1d-4)THEN
      aux=(x-y-z)*NMB0(x,y,z,t)
      ELSE
       IF(Max(y,z).ge.1d-4)THEN
      aux=(y+z)*FSF1(y,z,t)
       ELSE
      aux=0d0
       ENDIF
      ENDIF

      NMH=4d0*NMB22(x,y,z,t)+aux

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function rho2SM(x)

c            ->rho2SM(x)

      IMPLICIT NONE
      DOUBLE PRECISION x,Pi

      Pi=4d0*datan(1d0)

      rho2SM=19d0-33d0*x/2d0+43d0*x**2/12d0+7d0*x**3/120d0
     . -Pi*dsqrt(x)*(4d0-3d0/2d0*x+3d0/32d0*x**2+x**3/256d0)
     . -Pi**2*(2d0-2d0*x+x**2/2d0)-dlog(x)*(3d0*x-x**2/2d0)

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function BEW1(x,y,z)

c            ->BEW1(m1^2,m2^2,Q^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,aux

      IF(min(x,y).ge.1d-4)THEN
      IF(dabs(x-y).ge.1d-4)THEN
       aux=1d0+dlog(z/y)+(x/(x-y))**2*dlog(y/x)+(x+y)/(x-y)/2d0
      ELSE
       aux=dlog(z/x)
      ENDIF
      ELSE
       IF(x.ge.1d-4)aux=dlog(z/x)+3d0/2d0
       IF(y.ge.1d-4)aux=dlog(z/y)+1d0/2d0
      ENDIF

      BEW1=aux/2d0

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function CEW0(x,y,z)

c            ->CEW0(m1^2,m2^2,m3^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,aux

      IF(min(x,y,z).ge.1d-4)THEN
      IF(dabs(y-z).ge.1d-4)THEN
       IF(dabs(x-y).ge.1d-4)THEN
      IF(dabs(x-z).ge.1d-4)THEN
      aux=(y*dlog(y/x)/(x-y)-z*dlog(z/x)/(x-z))/(y-z)
      ELSE
      aux=(y*dlog(x/y)-x+y)/(x-y)**2
      ENDIF
       ELSE
      IF(dabs(x-z).ge.1d-4)THEN
      aux=(z*dlog(x/z)-x+z)/(x-z)**2
      ELSE
      aux=-1/2d0/x
        ENDIF
       ENDIF
      ELSE
       IF(dabs(x-y).ge.1d-4)THEN
      aux=(x*dlog(y/x)+x-y)/(x-y)**2
       ELSE
      aux=-1/2d0/x
       ENDIF
      ENDIF
      ELSE
       IF(x.ge.1d-4.and.y.ge.1d-4)THEN
      IF(dabs(y-x).ge.1d-4)THEN
      aux=dlog(y/x)/(x-y)
      ELSE
      aux=-1/x
      ENDIF
       ENDIF
       IF(x.ge.1d-4.and.y.ge.1d-4)THEN
      IF(dabs(y-x).ge.1d-4)THEN
      aux=dlog(y/x)/(x-y)
      ELSE
      aux=-1/x
      ENDIF
       ENDIF
       IF(x.ge.1d-4.and.z.ge.1d-4)THEN
      IF(dabs(z-x).ge.1d-4)THEN
      aux=dlog(z/x)/(x-z)
      ELSE
      aux=-1/x
      ENDIF
       ENDIF
       IF(y.ge.1d-4.and.z.ge.1d-4)THEN
      IF(dabs(z-y).ge.1d-4)THEN
      aux=dlog(z/y)/(y-z)
      ELSE
      aux=-1/y
      ENDIF
       ENDIF
      ENDIF

      CEW0=aux

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function DEW0(x,y,z)

c            ->DEW0(m1^2,m2^2,m3^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,aux

      IF(min(x,y,z).ge.1d-4)THEN
       IF(dabs(y-z).ge.1d-4)THEN
      IF(dabs(x-y).ge.1d-4)THEN
       IF(dabs(x-z).ge.1d-4)THEN
        aux=-(y*(x-y+x*dlog(y/x))*(x-z)**2
     . -z*(x-z+x*dlog(z/x))*(x-y)**2)/x/(y-z)/(x-y)**2/(x-z)**2
       ELSE
        aux=(x**2-y**2+2d0*x*y*dlog(y/x))/2d0/x/(x-y)**3
       ENDIF
      ELSE
       IF(dabs(x-z).ge.1d-4)THEN
        aux=(x**2-z**2+2d0*x*z*dlog(z/x))/2d0/x/(x-z)**3
       ELSE
        aux=1d0/6d0/x**2
       ENDIF
      ENDIF
       IF(dabs(x-y).ge.1d-4)THEN
        aux=-(2d0*(x-y)+(x+y)*dlog(x/y))/(x-y)**3
       ELSE
        aux=1d0/6d0/x**2
       ENDIF
       ENDIF
      ELSE
       IF(x.ge.1d-4)THEN
      IF(max(y,z).ge.1d-4)THEN
         aux=(-x+Max(y,z)-x*dlog(Max(y,z)/x))/x/(x-Max(y,z))**2
      ENDIF
       ENDIF
      ENDIF

      DEW0=aux

      RETURN
      END

**********************************************************************

      DOUBLE PRECISION function DEW27(x,y,z)

c            ->DEW27(m1^2,m2^2,m3^2)

      IMPLICIT NONE
      DOUBLE PRECISION x,y,z,aux

      IF(min(x,y,z).ge.1d-4)THEN
       IF(dabs(y-z).ge.1d-4)THEN
      IF(dabs(x-y).ge.1d-4)THEN
       IF(dabs(x-z).ge.1d-4)THEN
        aux=-(y**2*(x-z)**2*dlog(y/x)-z**2*(x-y)**2*dlog(z/x)
     .        +x*(x-y)*(x-z)*(y-z))/4d0/(y-z)/(x-y)**2/(x-z)**2
       ELSE
        aux=-(x**2-4d0*x*y+3d0*y**2-2d0*y**2*dlog(y/x))
     .          /8d0/(x-y)**3
       ENDIF
      ELSE
       IF(dabs(x-z).ge.1d-4)THEN
        aux=-(x**2-4d0*x*z+3d0*z**2-2d0*z**2*dlog(z/x))
     .          /8d0/(x-z)**3
       ELSE
        aux=-1d0/12d0/x
       ENDIF
      ENDIF
       IF(dabs(x-y).ge.1d-4)THEN
        aux=(-x**2+y**2-2d0*x*y*dlog(y/x))/4d0/(x-y)**3
       ELSE
        aux=-1d0/12d0/x
       ENDIF
       ENDIF
      ELSE
       IF(x.ge.1d-4)THEN
      IF(Max(y,z).ge.1d-4)THEN
         aux=(x-Max(y,z)+Max(y,z)*dlog(Max(y,z)/x))/4d0
     .          /(x-Max(y,z))**2
      ELSE
       aux=-1/4d0/x
      ENDIF
       ELSE
      IF(min(y,z).ge.1d-4)THEN
       IF(dabs(y-z).ge.1d-4)THEN
        aux=-dlog(y/z)/4d0/(y-z)
       ELSE
        aux=-1/4d0/y
       ENDIF
      ENDIF
       ENDIF
      ENDIF

      DEW27=aux

      RETURN
      END

ext: f


text:
	SUBROUTINE IPWTH ( icode, ix, iy, ixoff, iyoff, size, iret )
C************************************************************************
C* IPWTH								*
C*									*
C* This subroutine computes the line segments necessary to draw		*
C* past weather symbols.						*
C*									*
C* IPWTH ( ICODE, IX, IY, IXOFF, IYOFF, SIZE, IRET )			*
C*									*
C* Input parameters:							*
C*	ICODE		INTEGER		Symbol code			*
C*	IX		INTEGER		X coordinate in device units	*
C*	IY		INTEGER		Y coordinate in device units	*
C*	IXOFF		INTEGER		X offset in half characters	*
C*	IYOFF		INTEGER		Y offset in half characters	*
C*	SIZE		REAL		Symbol size			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* K. Brill		11/91						*
C* G. Krueger		 5/95	Fixed mwtwid -> mpwwid			*
C* M. Linda/GSC		 9/96	Moved most logic into ISYMB, cleaned up	*
C* M. Linda/GSC		12/96	Removed color fill			*
C************************************************************************
	INCLUDE		'DEVACT.CMN'
	INCLUDE		'ERROR.PRM'
C*
C------------------------------------------------------------------------
C
C*	Declare symbol shape table.
C
	INTEGER symtbl ( 0:32, 3:9 )
C
C------------------------------------------------------------------------
C
C*	Declare symbol matrices.
C
	INTEGER symbol3 (0:15), symbol4 (0: 6), symbol5 (0: 5)
	INTEGER symbol6 (0: 2), symbol7 (0: 6), symbol8 (0: 4)
	INTEGER symbol9 (0: 9)
C
C------------------------------------------------------------------------
C
C*	Associate symbol shape table with symbol matrices.
C
	EQUIVALENCE ( symtbl ( 0, 3 ), symbol3 )
	EQUIVALENCE ( symtbl ( 0, 4 ), symbol4 )
	EQUIVALENCE ( symtbl ( 0, 5 ), symbol5 )
	EQUIVALENCE ( symtbl ( 0, 6 ), symbol6 )
	EQUIVALENCE ( symtbl ( 0, 7 ), symbol7 )
	EQUIVALENCE ( symtbl ( 0, 8 ), symbol8 )
	EQUIVALENCE ( symtbl ( 0, 9 ), symbol9 )
C
C------------------------------------------------------------------------
C
C*	Define symbol shapes.
C@
C@	Each DATA statement defines the shape of one symbol.
C@
C@	The first element in each DATA statement is the number of
C@	elements that follow.  The subsequent elements define dots,
C@	points connected by lines, and shaded polygons.
C@
C@	Each symbol is formed out of dots, line segments, and polygons
C@	that are drawn within an 11 x 11 matrix.  Each X and Y coordinate
C@	pair is coded as a single number where values 100 to 1100
C@	represent the X coordinate, and values 1 to 11 represent the
C@	Y coordinate.  The origin of the matrix is the lower left
C@	corner, specified by 101.  Point 1111 represents the upper right
C@	corner.  The center of the symbol is point 606.
C@
C@	Negative coordinates imply moves to the specified location
C@	with the pen up.  Positive coordinates imply moves with the
C@	pen down.
C@
C@	Dots are specified as two co-located points (connected by a null
C@	line), and a multiple of 10,000 is added to the coordinate code.
C@	The dot size is the number of 10,000's.  For example, the dot
C@	specified by 40606 is twice as big as the dot specified by 20606.
c@	Coordinate sequence ...-30606, 30606... defines a single dot of
C@	size 3 that is located at point 606.
C@
C@	A zero (0) indicates a filled region defined by the shape that
C@	precedes the zero.  For example, ...-503, 504, 603, 503, 0,...
C@	is a shaded (filled in) triangle.
C@
C------------------------------------------------------------------------
C
	DATA symbol3   / 15, -808, 709, 509, 408, 407, 805, 804, 703,
     +			 503, 404, -306, 1006, 907, 1006, 905 /
C!	Blowing sand or snow
C
	DATA symbol4   / 6, -304, 904, -306, 906, -308, 908 /
C!	Fog
C
	DATA symbol5   / 5, -20607, 20607, -607, 606, 505 /
C!	Drizzle
C
	DATA symbol6   / 2, -30606, 30606 /
C!	Rain
C
	DATA symbol7   / 6, -406, 806, -505, 707, -507, 705 /
C!	Snow
C
	DATA symbol8   / 4, -307, 601, 907, 307 /
C!	Shower
C
	DATA symbol9   / 9, -407, 807, 605, 902, 903, 902, 802, -502,
     +			 507 /
C!	Thunderstorm
C
C------------------------------------------------------------------------
C
	iret = NORMAL
C
C*	Check symbol code validity.
C
	IF ( ( icode .lt. 3 ) .or. ( icode .gt. 9 ) ) RETURN
C
C*	Plot the specified symbol.
C
	iwidth = mpwwid
	CALL ISYMB ( ix, iy, ixoff, iyoff, size, iwidth, symtbl(0,icode),
     +		     iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!***********************************/
!	Name:
!	Alias:PointSwap
!	Description:
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.01.24
!	Update:
!	Other:not parallelize
!***********************************/
    subroutine MergeBesideCell(tGCC,dGCC)
    use StructVar_Mod
    use FrequentOperation
    implicit none

    !integer, intent(in) :: iLocalEdge !of Global Cell (=1,2,3,4)
    !integer, intent(in) :: iCenterX, iCenterY
    !type(Geometry), intent(in) :: Geom
    type(CutCell), intent(inout) :: tGCC, dGCC
    !type(InternalPointOfCell)
    type(InternalEdgeOfCell), pointer :: TempIEC
    integer :: iCell

        TempIEC => tGCC%IEC

        do while(associated(TempIEC%Next))
            TempIEC => TempIEC%Next
        end do

            TempIEC%Next => dGCC%IEC

        return
    end subroutine MergeBesideCell


text:
http://asucd.ucdavis.edu/media/aspapers AS PAPERs is actually an acronym for Associated Students Providing Alternative Publications Equal Resources, and served as a forum for students to create and run their own student run publications. AS PAPERs was created as an ASUCD Unit in 1993 following a student referendum. The referendum abolished then Unit, Third World Forum (TWF).

Over Winter break 2005, 30 years of http://www.californiaaggie.com/media/storage/paper981/news/2006/01/12/CampusNews/LossOf.As.Papers.Archives.Due.To.Communication.Failure1340706.shtml?norewrite200607270643&sourcedomainwww.californiaaggie.com AS PAPERs archives were thrown away by hired movers. The movers had been hired by the University to make storage space for Picnic Day. Many of the archives can still be found at Shields Library.

In 2009, the unit was cut from the ASUCD Budget. It http://theaggie.org/article/senatehearsnextyearsbudget continued to operate as a function of the ASUCD Vice President. In 2011, the unit was cut entirely from the ASUCD Budget.

 Past Publications
  La Palabra
  AlWasilah
  Third World Forum
  Vent Magazine
  the VOICE
  xPress
  Accents
  The Davis Enforcer
  Greek Times
  International Affairs Journal
  Jewish Quarterly
  The Sneeze
  Vox Feminista
  Advocacy Now!
  Awaaz
  The Business Letter
  The California Jewish Review
  Chilam Balam
  Ergo Sum
  Libertys Flame
  One Drop
  ouRevolution
  ouT LOUD

The Unit has three paid positions including the director and assistant director who are assigned duties to keep the office in working condition, and the Distributor, who distributes the publications across campus. The Director and Assistant Director also make budget plans and ensure those working in AS Papers feel welcome.

In theory, AS PAPERs is comprised of publications that operate fairly independently. But, the staff on each paper shares the office in Lower Freeborn, which houses six computers, a big printer, scanner, and two adjacent rooms; one a filing room, the other storage. The office also has two couches; the larger one is often used for students who dont get enough sleep. Although it has been rumored that students live in AS PAPERs, none actually do, although you may encounter students pulling all nighters working on their publications.

Like many bureaucratic enterprises, being part of AS PAPERs comes with its own set of rules. The publications are supposed to be official student organizations, registered with SPAC. Every publication must read (and supposedly, comply) with the AS PAPERs Admin Plan and sign a contract with AS PAPERs to become an active publication. New implementations for the 2006/2007 year include a required quarterpage AS PAPERs advertisement in any AS PAPERs publication as well as a bimonthly meeting that needs to be attended by a representative of each active publication. Publications that would like to seek additional funding, they must register and then go through the Club Finance Council. The papers usually use The Davis Enterprise or Kinkos to publish.

See also Local Publications.

20050615 20:48:08 nbsp What exactly does AS papers do? Provide a location where the designers can design the paper if they dont have the software? At first I saw them and thought they were a service that helped fund publications, but it seems thats still the creators responsibility. Hmm. I want to create a magazine/publication thing, but I dont know how to get the funding. Somehow I dont think SPAC is going to help. I did find a really neat little pamphlet in the bathroom today, though. Its called Phlegm Magazine and appears to only have 34 advertisers.... Users/MichaelGiardina
    20050615 21:02:26 nbsp It was created in 1993 as a backlash against Third World Forum.  TWF was a ASUCD funded paper that went out with the Aggie on a regular basis, and was extremely inflammatory in the White Man Bad sort of rhetoric.  It so alienated members of the ASUCD (presenate) that they took TWFs funding and created AS papers to in theory open it up to anyone, but it was more of a way of eliminating funding for TWF. I have no idea what they actually did, but it did severely curtail the publication of TWF. Users/RogerClark

20060727 04:14:18 nbsp I have been part of the AS PAPERs office for somewhere around a year now, and have had both good and bad experiences. I think its important to note that AS PAPERs was my first encounter with any sort of serious antiSemitism. Before working there, the most Id run into was a bad Jew joke. Generally speaking, the office is comprised of many different minorities, most of whom I have heard express views against Israel. Ive seen both subtle to more open hostility to white people and especially white Jews for no reason whatsoever, and I am not alone in observing this. Ive seen ASUCD turn a blind eye to this, cloaking it under the fact that AS PAPERs is a diverse environment. Why yes it is; one that fosters reverseracism. All that said, I dont mean to overly generallize. There have been quite a few people in this office that have been very kind, not to mention broadened my horizons and beliefs on some very interesting cultures. SS

    There is a Jewish publication, what do you suggest ASUCD should do?   Users/SteveOstrowski

        Quite honestly I think there is a lot ASUCD can probably do that it isnt doing. But I think for starters (since I know this is under ASUCD control), you guys should reconsider your hiring practices. I applied for Director around a year ago and Rob Roy told me I did not get the job because other candidates would be better qualified to handle the diverse environment. Im sure the ambiguity of that answer can lead to arguments, but let me just say that everyone else I saw applying for that job were women who were part of a minority culture, and for the past three directors, those are the people who have been hired. Mind you, maybe they were the most qualified ... I dont want you to think Im bringing this up because Im bitter over not getting the job. I just think more active scrutiny can be taken in the hiring process. For example: consider not hiring someone who is a member of SJP or Aggies For Israel. A political activist in a minority office that has very strong feelings about Israel would probably not serve well to also be a director.

        Im confused by your comment. Are you saying because were only one Jewish publication that it shouldnt be that big of deal? Because there are others outside of The California Jewish Review that have expressed concern. SS

            No, I was curious because there is a California Jewish Review which led me to believe there may be balance in AS Papers.  For example many conservatives say AS Papers is very liberal but I point out Libertys Flame.  But I see what you mean. Users/SteveOstrowski

20070507 01:07:37 nbsp wow. this page is interesting! ive rarely picked up one of the as paper issues, but the few times i did so i found the articles to be REALLY interesting. i liked how it had both international, national and local issues presented. thanks for existing! Users/JessicaRockwell

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20081115 15:54:31 nbsp Welcome to the Wiki Howdy Sanya and Welcome to the Wiki!     My names Evan, pleased to meet you!  Thanks for adding to the Olive Drive Market... it would be great it youd add things like specials to the entry itself.  Just click the Edit icon and you can edit the whole entry.  In fact, as an editor on the wiki, you can edit anything on the entire site.  If youre connected to the business, the wiki community wrote an informal introduction to the wiki that you can read Welcome to the Wiki/Business Owner by clicking here.  Once again, welcome to the wiki! Users/JabberWokky Evan JabberWokky Edwards


text:
! “ButterflyPACK” Copyright (c) 2018, The Regents of the University of California, through
! Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from the
! U.S. Dept. of Energy). All rights reserved.

! If you have questions about your rights to use or distribute this software, please contact
! Berkeley Lab's Intellectual Property Office at  IPO@lbl.gov.

! NOTICE.  This Software was developed under funding from the U.S. Department of Energy and the
! U.S. Government consequently retains certain rights. As such, the U.S. Government has been
! granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable
! worldwide license in the Software to reproduce, distribute copies to the public, prepare
! derivative works, and perform publicly and display publicly, and to permit other to do so.

! Developers: Yang Liu
!             (Lawrence Berkeley National Lab, Computational Research Division).

! This exmple works with double precision data
#define DAT 1

#include "ButterflyPACK_config.fi"

module APPLICATION_MODULE
use BPACK_DEFS
implicit none

	!**** define your application-related variables here
	type quant_app
		integer ntrain,ntest,diagflag ! size of training points and test points
		character(LEN=500)::trainfile_p,trainfile_tree,trainfile_l,trainfile_d,testfile_p,testfile_l !Kernel Regression: file pointers to train data, preordered tree, train labels, test data, and test labels
		real(kind=8)::lambda
		integer Nunk ! size of the matrix
		! real(kind=8),allocatable:: xyz(:,:)   ! coordinates of the points
		real(kind=8),allocatable:: matZ_glo(:,:)
		integer,allocatable:: perms(:)
		real(kind=8),allocatable:: diags(:)
	end type quant_app

contains

	!**** user-defined subroutine to sample Z_mn
	subroutine Zelem_FULL(m,n,value_e,quant)
		use BPACK_DEFS
		implicit none

		class(*),pointer :: quant
		integer, INTENT(IN):: m,n
		real(kind=8)::value_e

		real(kind=8) r_mn
		integer dimn

		select TYPE(quant)
		type is (quant_app)
			if(m<n)then
			value_e = quant%matZ_glo(m,n)
			else
			value_e = quant%matZ_glo(n,m)
			endif
			if(m==n)value_e = value_e +quant%lambda
		class default
			write(*,*)"unexpected type"
			stop
		end select
	end subroutine Zelem_FULL

end module APPLICATION_MODULE

PROGRAM ButterflyPACK_FullKRR
    use BPACK_DEFS
    use APPLICATION_MODULE

	use BPACK_structure
	use BPACK_factor
	use BPACK_constr
	use omp_lib
	use MISC_Utilities

    implicit none

	! include "mkl_vml.fi"

    real(kind=8) para
    real(kind=8) tolerance
    integer Primary_block, nn, mm,kk,mn,rank,ii,jj
    integer i,j,k, threads_num
	integer times(8)
	real(kind=8) t1,t2,x,y,z,r,theta,phi

	character(len=:),allocatable  :: string
	character(len=1024)  :: strings,strings1
	character(len=6)  :: info_env
	integer :: length
	integer :: edge_m,edge_n,flag,nargs
	integer :: ierr
	integer*8 oldmode,newmode
	type(Hoption)::option
	type(Hstat)::stats
	type(mesh)::msh
	type(kernelquant)::ker
	type(quant_app),target::quant
	type(Bmatrix)::bmat
	integer,allocatable:: groupmembers(:)
	integer nmpi
	type(proctree)::ptree
	CHARACTER (LEN=1000) DATA_DIR
	integer,allocatable::Permutation(:)
	integer Nunk_loc
	integer v_major,v_minor,v_bugfix
	real(kind=8),allocatable:: tmpline(:)

	!**** nmpi and groupmembers should be provided by the user
	call MPI_Init(ierr)
	call MPI_Comm_size(MPI_Comm_World,nmpi,ierr)
	allocate(groupmembers(nmpi))
	do ii=1,nmpi
		groupmembers(ii)=(ii-1)
	enddo

	!**** create the process tree
	call CreatePtree(nmpi,groupmembers,MPI_Comm_World,ptree)
	deallocate(groupmembers)

	if(ptree%MyID==Main_ID)then
    write(*,*) "-------------------------------Program Start----------------------------------"
    write(*,*) "ButterflyPACK_FullKRR"
    write(*,*) "   "
	endif

   !***********************************************************************
   if(ptree%MyID==Main_ID)then
   write (*,*) ''
   write (*,*) 'FULLKER computing'
   call BPACK_GetVersionNumber(v_major,v_minor,v_bugfix)
   write(*,'(A23,I1,A1,I1,A1,I1,A1)') " ButterflyPACK Version:",v_major,".",v_minor,".",v_bugfix
   write (*,*) ''
   endif
   !***********************************************************************

	!**** initialize stats and option
	call InitStat(stats)
	call SetDefaultOptions(option)

	!**** set solver parameters
	option%xyzsort=NATURAL
	option%nogeo=1
	option%Nmin_leaf=500
	option%ErrFillFull=0
	option%ErrSol=1

	!**** intialize the user-defined derived type quant
	quant%trainfile_p='../EXAMPLE/FULLMAT_DATA/FullMatKrr/kernel_train5k_test10k.txt'
	quant%trainfile_l='../EXAMPLE/FULLMAT_DATA/FullMatKrr/label_train5k_test10k.txt'
    quant%ntrain=5000
    quant%ntest=10000
    quant%lambda=0

	quant%diagflag=0
	nargs = iargc()
	ii=1
	do while(ii<=nargs)
		call getarg(ii,strings)
		if(trim(strings)=='-quant')then ! user-defined quantity parameters   !**** read data file directory, data set dimension, size of training and testing sets, RBF parameters and solver parameters
			flag=1
			do while(flag==1)
				ii=ii+1
				if(ii<=nargs)then
					call getarg(ii,strings)
					if(strings(1:2)=='--')then
						ii=ii+1
						call getarg(ii,strings1)
						if(trim(strings)=='--trainfile_p')then
							quant%trainfile_p=trim(strings1)
						else if	(trim(strings)=='--trainfile_l')then
							quant%trainfile_l=trim(strings1)
						else if	(trim(strings)=='--trainfile_d')then
							quant%trainfile_d=trim(strings1)
							quant%diagflag=1
						else if	(trim(strings)=='--ntrain')then
							read(strings1,*)quant%ntrain
						else if	(trim(strings)=='--ntest')then
							read(strings1,*)quant%ntest
						else if	(trim(strings)=='--lambda')then
							read(strings1,*)quant%lambda
						else
							if(ptree%MyID==Main_ID)write(*,*)'ignoring unknown quant: ', trim(strings)
						endif
					else
						flag=0
					endif
				else
					flag=0
				endif
			enddo
		else if(trim(strings)=='-option')then ! options of ButterflyPACK
			call ReadOption(option,ptree,ii)
		else
			if(ptree%MyID==Main_ID)write(*,*)'ignoring unknown argument: ',trim(strings)
			ii=ii+1
		endif
	enddo

	quant%Nunk = quant%ntrain

	call PrintOptions(option,ptree)

	t1 = OMP_get_wtime()
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "reading fullmatrix......"
	allocate(quant%perms(quant%ntrain+quant%ntest))
	do ii=1,quant%ntrain+quant%ntest
		quant%perms(ii)=ii
	enddo
	! call rperm(quant%ntrain+quant%ntest,quant%perms)

	allocate(tmpline(quant%ntrain))
	open (90,file=quant%trainfile_p)
	allocate (quant%matZ_glo(quant%ntrain+quant%ntest,quant%ntrain))
	do edge_m=1,quant%ntrain+quant%ntest
	read (90,*) tmpline
	do edge_n=1,quant%ntrain
		quant%matZ_glo(quant%perms(edge_m),quant%perms(edge_n)) = tmpline(edge_n)
	enddo
	enddo
	close(90)
	deallocate(tmpline)
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "reading fullmatrix finished"
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "    "
	t2 = OMP_get_wtime()

	!**** register the user-defined function and type in ker
	ker%QuantApp => quant
	ker%FuncZmn => Zelem_FULL

	!**** initialization of the construction phase
	allocate(Permutation(quant%Nunk))
	call BPACK_construction_Init(quant%Nunk,Permutation,Nunk_loc,bmat,option,stats,msh,ker,ptree)
	deallocate(Permutation) ! caller can use this permutation vector if needed

	!**** computation of the construction phase
    call BPACK_construction_Element(bmat,option,stats,msh,ker,ptree)

	!**** factorization phase
    call BPACK_Factorization(bmat,option,stats,ptree,msh)

	!**** solve phase
    call FULLKER_solve(bmat,option,msh,quant,ptree,stats)

	!**** deletion of quantities
	call delete_proctree(ptree)
	call delete_Hstat(stats)
	call delete_mesh(msh)
	call delete_kernelquant(ker)
	call BPACK_delete(bmat)

    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "-------------------------------program end-------------------------------------"

	call blacs_exit(1)
	call MPI_Finalize(ierr)

end PROGRAM ButterflyPACK_FullKRR

subroutine FULLKER_solve(bmat,option,msh,quant,ptree,stats)

    use BPACK_DEFS
	use APPLICATION_MODULE
	use omp_lib
	use BPACK_Solve_Mul
	use MISC_Utilities

    implicit none

    integer i, j, ii, jj, iii, jjj, ierr, ntest,Dimn,edge_m,edge_n,ncorrect
    integer level, blocks, edge, patch, node, group
    integer rank, index_near, m, n, length, flag, num_sample, n_iter_max, iter, N_unk, N_unk_loc
    real(kind=8) theta, phi, dphi, rcs_V, rcs_H,error
    real(kind=8) T0,T1
    real(kind=8) n1,n2,rtemp
    real(kind=8) value_Z
    real(kind=8),allocatable:: Voltage_pre(:),x(:,:),b(:,:),vout(:,:),vout_tmp(:,:),vout_test(:,:),matrixtemp1(:,:),matrixtemp2(:,:),matrixtemp(:,:)
	real(kind=8):: rel_error
	type(Hoption)::option
	type(mesh)::msh
	type(quant_app)::quant
	type(proctree)::ptree
	type(Bmatrix)::bmat
	type(Hstat)::stats
	real(kind=8),allocatable:: current(:),voltage(:)
	real(kind=8), allocatable:: labels(:)
	integer, allocatable :: ipiv(:)
	real(kind=8) r_mn
	real(kind=8) label

    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "Solve and Prediction......"

	N_unk=msh%Nunk
	N_unk_loc = msh%idxe-msh%idxs+1
	ntest=quant%ntest

	if(option%ErrSol==1)then
		call BPACK_Test_Solve_error(bmat,N_unk_loc,option,ptree,stats)
	endif

	!**** read training label as local RHS and solve for the weights

	allocate(labels(quant%ntrain+quant%ntest))
	allocate (x(N_unk_loc,1))
	x=0
	allocate (b(N_unk_loc,1))
	open (91,file=quant%trainfile_l)
	! read(91,*)N_unk,Dimn
	do ii=1,quant%ntrain+quant%ntest
		read(91,*)labels(quant%perms(ii))
	enddo
	close(91)
	do ii=1,N_unk_loc
		b(ii,1) = labels(msh%new2old(ii-1+msh%idxs))
	enddo
	allocate (vout_test(ntest,1))
	vout_test(:,1) = labels(1+quant%ntrain:quant%ntrain+quant%ntest)
	deallocate(labels)

	n1 = OMP_get_wtime()

	call BPACK_Solution(bmat,x,b,N_unk_loc,1,option,ptree,stats)

	n2 = OMP_get_wtime()
	stats%Time_Sol = stats%Time_Sol + n2-n1
	call MPI_ALLREDUCE(stats%Time_Sol,rtemp,1,MPI_DOUBLE_PRECISION,MPI_MAX,ptree%Comm,ierr)
	if(ptree%MyID==Main_ID .and. option%verbosity>=0)write (*,*) 'Solving:',rtemp,'Seconds'
	call MPI_ALLREDUCE(stats%Flop_Sol,rtemp,1,MPI_DOUBLE_PRECISION,MPI_SUM,ptree%Comm,ierr)
	if(ptree%MyID==Main_ID .and. option%verbosity>=0)write (*,'(A13Es14.2)') 'Solve flops:',rtemp

	!**** prediction on the test sets

	if(quant%diagflag==1)then
		allocate(quant%diags(quant%ntrain))
		open (92,file=quant%trainfile_d)
		! read(91,*)N_unk,Dimn
		do ii=1,quant%ntrain
			read(92,*)quant%diags(quant%perms(ii))
		enddo
		close(92)
		do edge=1, N_unk_loc
			x(edge,1) = x(edge,1)*quant%diags(msh%new2old(edge+msh%idxs-1))
		enddo
		deallocate(quant%diags)
	endif

	T0 = OMP_get_wtime()

	allocate (vout(ntest,1))
	allocate (vout_tmp(ntest,1))
	vout_tmp = 0
	do edge=1, N_unk_loc
		do edge_m=1,ntest
			value_Z = quant%matZ_glo(edge_m+quant%ntrain,msh%new2old(edge+msh%idxs-1))
			vout_tmp(edge_m,1) = vout_tmp(edge_m,1) + value_Z*x(edge,1)
		enddo
	enddo

	call MPI_REDUCE(vout_tmp, vout, ntest,MPI_double_precision, MPI_SUM, Main_ID, ptree%Comm,ierr)

		! allocate(matrixtemp2(N_unk,N_unk))
		! allocate(matrixtemp1(N_unk,N_unk))
		! matrixtemp2=quant%matZ_glo(msh%new2old,msh%new2old)
		! call GeneralInverse(N_unk,N_unk,matrixtemp2,matrixtemp1,1d-10)

		! allocate(ipiv(N_unk))
		! allocate(matrixtemp1(N_unk,N_unk))
		! matrixtemp1=quant%matZ_glo(msh%new2old,msh%new2old)
		! ipiv=0
		! call getrff90(matrixtemp1,ipiv)
		! call getrif90(matrixtemp1,ipiv)
		! deallocate(ipiv)

		! allocate(matrixtemp(ntest,N_unk))
		! matrixtemp=quant%matZ_glo(quant%ntrain+1:quant%ntrain+quant%ntest,msh%new2old)
		! call gemmf90(matrixtemp1,N_unk,b,N_unk,x,N_unk,'N','N',N_unk,1,N_unk,BPACK_cone,BPACK_czero)
		! allocate (vout(ntest,1))
		! vout=0
		! call gemmf90(matrixtemp,ntest,x,N_unk,vout,ntest,'N','N',ntest,1,N_unk,BPACK_cone,BPACK_czero)

	T1 = OMP_get_wtime()
	if (ptree%MyID==Main_ID) then
		vout_tmp = vout-vout_test
		error=0
		do ii=1,ntest
			error = error + vout_tmp(ii,1)**2d0
		enddo
		error = error/ntest
		error = sqrt(error)

		write (*,*) ''
		write (*,*) 'Prediction time:',T1-T0,'Seconds'
		write (*,*) 'Prediction error:',error
		write (*,*) ''

	endif

	deallocate (vout)
	deallocate (vout_tmp)
	deallocate (vout_test)

	deallocate(x)
	deallocate(b)

	if(allocated(quant%matZ_glo))deallocate(quant%matZ_glo)
	if(allocated(quant%perms))deallocate(quant%perms)

	call MPI_barrier(ptree%Comm,ierr)

    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "Solve and Prediction finished"
    if(ptree%MyID==Main_ID .and. option%verbosity>=0)write(*,*) "    "

    return

end subroutine FULLKER_solve

max_stars_repo_name: liuyangzhuan/hod-lr-bf

max_issues_repo_name: liuyangzhuan/hod-lr-bf

max_forks_repo_name: liuyangzhuan/hod-lr-bf


text:
      subroutine fwdcpfft2(psi,nr1,nr2,nr3)
c------------------------------------------------------
c
c Written by A. Canning (CRAY-EPFL) 25th July 94
c
c
c output = psi(x,y,z)  wavefunction in fourier space load balanced
c                      as small sphere columns in x direction
c input = psi(z,y,x) real space grid with each PE having consecutive
c                    sets of z columns.
c
c computes forward fft specifically for CP algo
c ie taking cube and going to two slices (one slice for gamma)
c then  to a sphere.  See fftprep.f for more details
c fftprep must be called once before this routine
c for setups.

      use fft_data
      use load_data
      use data

      implicit none

      include  'param.escan_real'
      include "mpif.h"

      real*8,allocatable,dimension(:) :: worknr1,worknr2,worknr3
      real*8,allocatable,dimension(:) :: psiy,combuf1,combuf2

      real*8 psi(mr_n)
      real*8 fac,s

c scalars used

      integer i,ib,ic,idum,ii,iloc,ilocadd,isc,isign,itar,itaradd,
     c        itnode,iw,ix,iy,izb,j,jcol,ngy,ngyadd,ngz,nr1,nr2,nr3,
     c        nr3u,iloc_dum,ierr
      integer nworknr1,nworknr2,nworknr3

c mpi version arrays
c
      integer mpistatus(mpi_status_size)
      integer ireq(nnodes)

      call mpi_barrier(mpi_comm_k,ierr)
c
      nworknr1 = 20000+2.28*nr1x
      nworknr2 = 20000+2.28*nr2x
      nworknr3 = 20000+2.28*nr3x

      allocate(worknr1(nworknr1))
      allocate(worknr2(nworknr2))
      allocate(worknr3(nworknr3))
      allocate(psiy(mr_n))
      allocate(combuf1(mr_n))
      allocate(combuf2(mr_n))

c
c for gamma point fft
c
      nr3u = nr3/2 + 1
c
c memory allocation for local arrays used with SHMEM routines
c with SHMEM routines address on all PE's must be the same for
c same variable
c
c
      isign = 1
c
c full FFT on the cube psi this is on the z dir
c layout is nz,ny,nx  colulmns of height nr3 complex
c reslt not put directly into psi (problem for scfft)
c due to different sizes of input and output
c

       do i = 1,ncolz2
        ilocadd = 1+(i-1)*nr3
       iloc_dum = 1+(i-1)*(2*nr3u)

        call system_scfft(0,isign,nr3,1.0d0,psi(ilocadd),
     &    psiy(iloc_dum),tabnr3lrc,worknr3,0,ntabnr3,nworknr3)
       enddo
c
c now transpose nz,ny,nx to ny,nz,nx in the two slice mode
c into psiy. Each PE will have ncoly columns in psiy
c

      idum = 1
      do i = 1,nnodes
        do j = 1,ivunpn2l(i)
          combuf1(idum) = psiy(ivunp2l(idum))
          idum = idum + 1
        enddo
      enddo

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
       call mpi_isend(combuf1(idum),ivunpn2l(i),mpi_real8,i-1,inode,
     &                mpi_comm_k,ireq(i),ierr)
       idum = idum + ivunpn2l(i)
      enddo

      idum = 1
      do i = 1,nnodes
       call mpi_recv(combuf2(idum),ivpacn2l(i),mpi_real8,i-1,i,
     &               mpi_comm_k,mpistatus,ierr)
       idum = idum + ivpacn2l(i)
      enddo

      do i = 1, nnodes
         call mpi_wait(ireq(i), mpistatus, ierr)
      end do

      idum = 1
      do i = 1,nnodes
        do j = 1,ivpacn2l(i)
          psiy(ivpac2l(idum)) = combuf2(idum)
          idum = idum + 1
        enddo
      enddo

c
c now do FFT's on the one slice
c
       do i = 1,ncoly2
        ilocadd = 1+(i-1)*(2*nr2)
       call system_ccfft(0,isign,nr2,1.0d0,psiy(ilocadd),
     c     psiy(ilocadd),tabnr2lfw,worknr2,0,ntabnr2,nworknr2)
       enddo
c
c now transpose back to format of program into psi
c ie into x columns load balanced
c
      idum = 1
      do i = 1,nnodes
        do j = 1,ivunpn1l(i)
          combuf1(idum) = psiy(ivunp1l(idum))
          idum  = idum + 1
        enddo
      enddo

      call mpi_barrier(mpi_comm_k,ierr)

      idum = 1
      do i = 1,nnodes
       call mpi_isend(combuf1(idum),ivunpn1l(i),mpi_real8,i-1,inode,
     &                mpi_comm_k,ireq(i),ierr)
       idum = idum + ivunpn1l(i)
      enddo

      idum = 1
      do i = 1,nnodes
       call mpi_recv(combuf2(idum),ivpacn1l(i),mpi_real8,i-1,i,
     &               mpi_comm_k,mpistatus,ierr)
       idum = idum + ivpacn1l(i)
      enddo

      do i = 1, nnodes
         call mpi_wait(ireq(i), mpistatus, ierr)
      end do

      idum = 1
      do i = 1,nnodes
        do j = 1,ivpacn1l(i)
          psi(ivpac1l(idum)) = combuf2(idum)
          idum = idum + 1
        enddo
      enddo

c
c now do last FFT in place in the format of the program
c for psi in the x direction
c
       do jcol = 1,ncol2(inode)
        ilocadd = 1+(jcol-1)*(2*nr1)
        call system_ccfft(0,isign,nr1,1.0d0,psi(ilocadd),
     c    psi(ilocadd),tabnr1lfw,worknr1,0,ntabnr1,nworknr1)
       enddo !i
c
c now do the scaling of the FFT
c
      fac=1.d0/dfloat(nr1*nr2*nr3)

      do i = 1,ncol2(inode)*2*nr1
        psi(i) = psi(i)*fac
      enddo

      call mpi_barrier(MPI_COMM_K,ierr)

      deallocate(worknr1)
      deallocate(worknr2)
      deallocate(worknr3)
      deallocate(psiy)
      deallocate(combuf1)
      deallocate(combuf2)

      end

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE GRFNEU (NPTS, XPTS, YPTS, PLTITL, TXLAB, TYLAB)
C=======================================================================

C   --*** WRTNEU *** (XYPLOT) Write curve to neutral file
C   --   Written by Amy Gilkey - revised 04/21/88
C   --
C   --WRTNEU writes the data for a curve to a neutral file.  The first
C   --time the routine is called, the neutral file is opened.
C   --
C   --The format of the neutral file is described in "GRAFAID Code User
C   --Manual" under Neutral File Format.
C   --
C   --Parameters:
C   --   NPTS - IN - the number of points on the curve
C   --   XPTS, YPTS - IN - the points on the curve
C   --   PLTITL - IN - the plot title describing the curve
C   --      (e.g. "TIME vs SIGXX at ELEMENT 30")
C   --   TXLAB, TYLAB - IN - the X and Y axis labels, either the
C   --      user-input labels or the plot variable descriptions
C   --
C   --Common Variables:
C   --   Uses TITLE, CREATE, MODIFY, DRAW of /DBTITL/
C   --   Uses DOQA, CAPTN of /LEGOPT/
C   --   Uses XMIN, XMAX, YMIN, YMAX of /LIMITS/
C   --   Uses NEU, NUMCRV, INCCRV, CRVNAM of /NEUTR./
C   --   Uses and sets NEUOPN of /NEUTR./

      include 'dbname.blk'
      include 'dbtitl.blk'
      include 'legopt.blk'
      include 'xylim.blk'
      include 'neutr.blk'

      REAL XPTS(NPTS), YPTS(NPTS)
      CHARACTER*(*) PLTITL
      CHARACTER*(*) TXLAB, TYLAB

      CHARACTER*2048 FILNAM, ERRMSG
      CHARACTER*4 XTYP
      CHARACTER*4 AXTYP
      CHARACTER AUX
      CHARACTER*80 CURVE

      DATA AUX /'F'/
      DATA AXTYP /'NOLO'/

      IF (.NOT. GRFOPN) THEN

C      --Open the neutral file and write the title line

        filnam = basenam(:lenstr(basenam)) // '.grf'

C      --Open the neutral file and write the title line
        open (unit=neugrf, file=filnam(:lenstr(filnam)),
     *    form='formatted', status='unknown', iostat=ierr)
        IF (IERR .NE. 0) THEN
          ERRMSG = 'Neutral file "'//FILNAM(:LENSTR(FILNAM))//
     *      '" could not be opened.'
          CALL PRTERR ('CMDERR', ERRMSG(:LENSTR(ERRMSG)))
          GOTO 170
        END IF
        GRFOPN = .TRUE.

        WRITE (NEUGRF, 10000) DRAW(1)(:LENSTR(DRAW(1))),
     &    CREATE(1)(:LENSTR(CREATE(1))), CREATE(3), CREATE(4),
     &    MODIFY(1)(:LENSTR(MODIFY(1))), MODIFY(3), MODIFY(4)
10000   FORMAT (A, ': CREATED BY ', A, ' ', A, ' ', A,
     &    ', MODIFIED BY ', A, ' ', A, ' ', A)
      END IF

C   --Figure out if the X points are monotonic or not

      XTYP = 'MONO'
      DO 100 I = 2, NPTS
         IF (XPTS(I-1) .GE. XPTS(I)) THEN
            XTYP = 'NONM'
            GOTO 110
         END IF
  100 CONTINUE
  110 CONTINUE

C   --Get the curve name

      WRITE (CURVE, '(A32, I7)', IOSTAT=IDUM) CRVNAM, NUMCRV
      CALL PCKSTR (1, CURVE)
      NUMCRV = NUMCRV + INCCRV

      WRITE (*, 10090) 'Writing ', PLTITL(:LENSTR(PLTITL))

C   --Write the begin curve record with the curve name

      WRITE (NEUGRF, 10010) 'BEGIN CURVE', CURVE
10010  FORMAT (A, ',', A15)

C   --Write the title lines

      DO 120 IEND = 3, 1, -1
         IF (CAPTN(IEND,2) .NE. ' ') GOTO 130
  120 CONTINUE
  130 CONTINUE
      IF (DOQA(2)) THEN
         WRITE (NEUGRF, 10020) 1+MAX(IEND,1), TITLE
         IF (IEND .GT. 0) THEN
            DO 140 I = 1, IEND
               WRITE (NEUGRF, 10030) CAPTN(I,2)
  140       CONTINUE
         ELSE
            WRITE (NEUGRF, 10030) PLTITL
         END IF
      ELSE IF (IEND .GT. 0) THEN
         WRITE (NEUGRF, 10020) IEND, CAPTN(1,2)
         DO 150 I = 2, IEND
            WRITE (NEUGRF, 10030) CAPTN(I,2)
  150    CONTINUE
      ELSE
         WRITE (NEUGRF, 10020) 1, PLTITL
      END IF
10020  FORMAT (I1, ',', A80)
10030  FORMAT (A80)

C   --Write the X and Y labels

      WRITE (NEUGRF, 10040) TXLAB
      WRITE (NEUGRF, 10040) TYLAB
10040  FORMAT (A40)

C   --Write the min/max, the number of points, and the auxiliary data flag

      WRITE (NEUGRF, 10050) XMIN, XMAX, YMIN, YMAX, NPTS, AUX
10050  FORMAT (4 (1PE15.7E3, ','), I5, ',', A1)
      WRITE (NEUGRF, 10060) AXTYP, XTYP, ' '
10060  FORMAT (A4, ',', A4, ',', A4)

C   --Write the data points

      DO 160 I = 1, NPTS
         WRITE (NEUGRF, 10070) XPTS(I), YPTS(I)
10070     FORMAT (2 (1PE15.7E3, :, ','))
  160 CONTINUE

C   --Write the end curve record with the curve name

      WRITE (NEUGRF, 10080) 'END CURVE', CURVE
10080  FORMAT (A, ',', A)

  170 CONTINUE
      RETURN
10090  FORMAT (1X, 5A)
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
      MODULE mpcbds_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  08:35:05  03/09/06
      SUBROUTINE mpcbds (BONDAB)
      USE vast_kind_param,ONLY: DOUBLE
      use molkst_C, only : NUMAT
      REAL(DOUBLE), DIMENSION((NUMAT*(NUMAT + 1))/2), INTENT(IN) :: BONDAB
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


module n_bodies
    use integrators
    implicit none
    contains

    function adjust_timestep(dt, err, q, t_tot)
        implicit none
        real(8), intent(in) :: dt, err, t_tot
        integer(8), intent(in) :: q
        real(8) :: adjust_timestep
        real(8) :: fac
!        real :: facmax=6., facmin=0.01
!        real :: fac
        fac = 0.38**(1. / (1. + q))
        adjust_timestep = dt * fac * err**(-1. / (1. + q))
!        print *, "Adjust timestep", adjust_timestep, dt, fac, err
!        if (adjust_timestep > facmax * dt) then
!            adjust_timestep = facmax * dt
!        else if (adjust_timestep < facmin * dt) then
!            adjust_timestep = facmin * dt
        if ((adjust_timestep / t_tot < 1e-12) .or. isnan(adjust_timestep)) then
            adjust_timestep = t_tot * 1e-12
        end if
    end function adjust_timestep

    function center_of_mass(n, masses, positions)
        implicit none
        integer(8), intent(in) :: n
        real(8), intent(in), dimension(n) :: masses
        real(8), intent(in), dimension(n, 3) :: positions
        real(8), dimension(3) :: center_of_mass
        integer(8) :: i
        center_of_mass = 0.
        do i=1, n
            center_of_mass(:) = center_of_mass(:) + masses(i) * positions(i, :)
        end do
    end function center_of_mass

    function average_vel(n, masses, velocities)
        implicit none
        integer(8), intent(in) :: n
        real(8), intent(in), dimension(n) :: masses
        real(8), intent(in), dimension(n, 3) :: velocities
        real(8), dimension(3) :: average_vel
        integer(8) :: i
        average_vel = 0.
        do i=1, n
            average_vel(:) = average_vel(:) + masses(i) * velocities(i, :)
        end do
    end function average_vel

    subroutine n_body_model(integrator, n, n_times, masses, i_pos, i_vel, &
                            t_tot, dt, g, update_com, times, positions, &
                            velocities, &
                            status)
        implicit none
        character(265), intent(in) :: integrator
        integer(8) :: n, n_times, n_times_internal
        real(8), dimension(n), intent(in) :: masses
        real(8), dimension(n, 3), intent(in) :: i_pos, i_vel
        real(8), dimension(n, 3) :: poses, vels, null_poses, null_vels
        real(8), dimension(n, 3) :: tmp_poses, tmp_vels, accels
        real(8), dimension(3) :: com, avv
        real(8), intent(in) :: t_tot, dt
        real(8), intent(in) :: g
        logical, intent(in) :: update_com
        real(8) :: t_now
        integer(8) :: t, i, j
        real(8), dimension(n_times + 1), intent(out) :: times
        ! shape n_steps, n, 3
        real(8), dimension(n_times + 1, n, 3), intent(out) :: positions, velocities
        character(256), intent(out) :: status
        status = "Ok"

        if (n_times < int(t_tot/dt)) then
            print *, "SizeError, n_times too small for t_tot/dt"
            status = "SizeError, n_times too small for t_tot/dt"
            stop
        else
            n_times_internal = floor(t_tot / dt)
        end if
!        allocate(times(n_steps))
!        allocate(positions(n_steps, n, 3))
!        allocate(velocities(n_steps, n, 3))
!
!
        ! Centers the system inputs
        if (update_com) then
            com(:) = center_of_mass(n, masses, i_pos)
            avv(:) = average_vel(n, masses, i_vel)
            do i=1, n
                poses(i, :) = i_pos(i, :) - com
                vels(i, :) = i_vel(i, :) - avv
            end do
        end if

        timeloop: do t=1, n_times_internal
!            print *, t
            times(t) = t_now
            do i=1, n
                positions(t, i, :) = poses(i, :)
!                print *, poses(i, :)
                velocities(t, i, :) = vels(i, :)
            end do
!             Run the integrator of choice for a timestep, using temporary
!             variables to avoid simultaneous IO for variables
            if (integrator(1:len("euler")) == "euler") then
                call euler_forward(n, masses, poses, vels, dt, g, tmp_poses, &
                        tmp_vels)
!                print *, poses(2, :), tmp_poses(2, :)
!                print *, vels(2, :), tmp_vels(2, :)
                do i=1, n
                    poses(i, :) = tmp_poses(i, :)
                    vels(i, :) = tmp_vels(i, :)
                end do
            else if (integrator(1:len("kdk")) == "kdk") then
                call kdk(n, masses, poses, vels, dt, g, tmp_poses, tmp_vels)
                poses = tmp_poses
                vels = tmp_vels
            else if (integrator(1:len("rk3")) == "rk3") then
                call rk3_classic(n, masses, poses, vels, dt, g, tmp_poses, tmp_vels)
                poses = tmp_poses
                vels = tmp_vels
            else if (integrator == "rk4") then
                call rkf45(n, masses, poses, vels, dt, g, tmp_poses, tmp_vels, &
                           null_poses, null_vels)
                poses = tmp_poses
                vels = tmp_vels
            else if (integrator == "rk5") then
                call rkf45(n, masses, poses, vels, dt, g, null_poses, null_vels, &
                        tmp_poses, tmp_vels)
                poses = tmp_poses
                vels = tmp_vels
            end if
            t_now = t_now + dt

!             Check validity of the integrator's calculation
            do i=1, n
                do j=1, 3
                    if (isnan(poses(i, j))) then
                        print *, "NAN encountered for position at", i, j
                        status = "NAN encountered for position"
                        exit timeloop
                    end if
                    if (isnan(vels(i, j))) then
                        print *, "NAN encountered for velocity at", i, j
                        status = "NAN encountered for velocity"
                        exit timeloop
                    end if
                end do
            end do
!
!            print *, vels(2, :)
            ! Recenters the system
            if (update_com) then
                com(:) = center_of_mass(n, masses, poses)
                avv(:) = average_vel(n, masses, vels)
                do i=1, n
                    poses(i, :) = poses(i, :) - com(:)
                    vels(i, :) = vels(i, :) - avv(:)
                end do
            end if

!            print *, vels(2, :)

        end do timeloop

    end subroutine n_body_model

    subroutine adaptive_n_body_model(integrator, n, n_times, masses, i_pos, &
                                     i_vel, t_tot, tolerance, dt_0, q, g, &
                                     update_com, times, positions, &
                                     velocities, status)
        implicit none
        character(265), intent(in) :: integrator
        integer(8) :: n, n_times, n_times_internal
        real(8), dimension(n), intent(in) :: masses
        real(8), dimension(n, 3), intent(in) :: i_pos, i_vel
        real(8), dimension(n, 3) :: poses, vels, high_poses, high_vels
        real(8), dimension(n, 3) :: low_poses, low_vels, accels
        real(8), dimension(3) :: com, avv
        real(8), intent(in) :: t_tot, tolerance, g, dt_0
        integer(8), intent(in) :: q
        logical, intent(in) :: update_com
        real(8) :: t_now, err_sum, dt
        integer(8) :: t, i, j
        real(8), dimension(n_times + 1), intent(out) :: times
        ! shape n_steps, n, 3
        real(8), dimension(n_times + 1, n, 3), intent(out) :: positions, velocities
        character(256), intent(out) :: status

        status = "Ok"
        dt = dt_0
        if (n_times < int(t_tot/dt)) then
            print *, "SizeRisk, n_times too small for t_tot/dt=", int(t_tot/dt)
            status = "SizeRisk, n_times too small for t_tot/dt="
        else
            n_times_internal = floor(t_tot / dt)
        end if
       ! allocate(times(n_steps))
       ! allocate(positions(n_steps, n, 3))
       ! allocate(velocities(n_steps, n, 3))
!
!
        ! Centers the system inputs
        if (update_com) then
            com(:) = center_of_mass(n, masses, i_pos)
            avv(:) = average_vel(n, masses, i_vel)
            do i=1, n
                poses(i, :) = i_pos(i, :) - com
                vels(i, :) = i_vel(i, :) - avv
            end do
        end if
!        print *, poses(2, :)
!        print *, "Looping for...", n_times_internal
        timeloop: do t=1, n_times_internal
       !     print *, t
            if (t_now > t_tot) then
                exit timeloop
            end if
            times(t) = t_now
            do i=1, n
                positions(t, i, :) = poses(i, :)
!                print *, poses(i, :), vels(i, :)
                velocities(t, i, :) = vels(i, :)
            end do
       !      Run the integrator of choice for a timestep, using temporary
       !      variables to avoid simultaneous IO for variables
            if (integrator(1:len("kdkrk3")) == "kdkrk3") then
                call kdkrk3(n, masses, poses, vels, dt, g, low_poses, low_vels, &
                            high_poses, high_vels)
                poses = high_poses
                vels = high_vels
            else if (integrator(1:len("rkf45")) == "rkf45") then
                call rkf45(n, masses, poses, vels, dt, g, low_poses, low_vels, &
                           high_poses, high_vels)
!                print *, low_poses(2, :), high_poses(2, :)
                do i=1, n
                    poses(i, :) = high_poses(i, :)
                    vels(i, :) = high_vels(i, :)
                end do
            end if
!
       !      Check validity of the integrator's calculation
            do i=1, n
                do j=1, 3
                    if (isnan(poses(i, j))) then
                        print *, "NAN encountered for position at", i, j
                        status = "NAN encountered for position"
                        exit timeloop
                    end if
                    if (isnan(vels(i, j))) then
                        print *, "NAN encountered for velocity at", i, j
                        status = "NAN encountered for velocity"
                        exit timeloop
                    end if
                end do
            end do

       !     print *, vels(2, :)
            ! Recenters the system
            if (update_com) then
                com(:) = center_of_mass(n, masses, poses)
                avv(:) = average_vel(n, masses, vels)
                do i=1, n
                    poses(i, :) = poses(i, :) - com(:)
                    vels(i, :) = vels(i, :) - avv(:)
                end do
            end if
!
       !     print *, vels(2, :)
            err_sum = 0.
            do i=1, n
                err_sum = err_sum + vec_mag(3, (high_poses - low_poses) &
                        / tolerance)
                err_sum = err_sum + vec_mag(3, (high_vels - low_vels) &
                        / tolerance)
            end do
            err_sum = err_sum / real(n * 3 * 2)**0.5
!            print *, tolerance
!
            dt = adjust_timestep(dt, err_sum, q, t_tot)
!            print *, dt
!
            t_now = t_now + dt
!
!
         end do timeloop

    end subroutine adaptive_n_body_model
end module n_bodies


text:
!! R845 if-stmt
!    is  IF ( scalar-logical-expr ) action-stmt
!
! Not tested here are: scalar-logical-expr and action-stmt.
if(x.eq.y) x = 1

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
      subroutine fem_vonMises_stress_update(i, j, gp_stress, gp_strain,
     & iUnload )
!***********************************************************************
! Copyright 2011 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!
! Computes the stresses from strains for the von Mises plasticity model
!
! Author : Sai Rapaka
!

      use comsi, only: iPlastic, plasticModel, modelNumber, e1, e2, e3
      use comsi, only: isPlastic, plasticParam1, plastic_strain
      use comsi, only: du, dv, dw, alp, bulk
      use comdi, only: t, phi, tini, phini
      use comai, only: iout, iptty, iad, ns, nei
      use comfem

      implicit none
      integer                      :: i, j
      real*8,  dimension(6)        :: gp_stress, gp_strain, dstrain
      real*8,  dimension(6)        :: tot_strain, dev_strain
      integer, dimension(8)        :: node
      real*8,  dimension(6, 24)    :: B
      real*8,  dimension(24)       :: disp
      logical                      :: iUnload
      real*8,  dimension(8)        :: alpha, beta, deltaT, deltaP
      real*8,  dimension(8)        :: lambda, G
      real*8,  dimension(6,6)      :: D
      real*8                       :: alphadeltaT, betadeltaP
      real*8 vol_strain, pressure
      real*8 bulk_mod, lambda_bar, G_bar, J2, Q, deltaGamma
      real*8 phi_trial, phi_tol, fac, sigma_y
      integer itmp, iModel
      integer k

      if(iPlastic.eq.0) then
        write(iout,*) '***ERROR : von Mises stress update called
     &   without plastic flag being set '
      endif

      if(.not.allocated(conv_pstrain)) then
        allocate(conv_pstrain(nei, numgausspoints))
        conv_pstrain = 0.0d0
      endif

      itmp = modelNumber(elnode(i, 1))
      iModel = plasticModel(itmp)

      do k=2,ns
        itmp = modelNumber(elnode(i, k))
        if(iModel.ne.plasticModel(itmp)) then
          write(iout, *) 'Multiple plastic models being used !
     &        Not supported at this time! '
          write(iptty, *) 'Multiple plastic models being used !
     &        Not supported at this time! '
          stop
        endif
      enddo

      phi_tol = 1.0d-6
      isPlastic(i,j) = 0

      bulk_mod = 0.0d0
      lambda_bar = 0.0d0
      G_bar = 0.0d0

      sigma_y = plasticParam1(itmp)

      alphadeltaT = 0.0d0
      betadeltaP  = 0.0d0

      do k=1, ns
        node(k) = elnode(i, k)
        lambda(k) = e2(node(k))
        G(k) = 0.5d0*(e1(k) - e2(k))
        lambda_bar = lambda_bar + lambda(k)*Psi(i,j,k)
        G_bar = G_bar + G(k)*Psi(i,j,k)
        alpha(k) = alp(node(k))
        deltaT(k) = t(node(k)) - tini(node(k))
        beta(k) = bulk(node(k))
        deltaP(k) = phi(node(k)) - phini(node(k))
        alphadeltaT = alphadeltaT +
     &                 Psi(i,j,k)*alpha(k)*deltaT(k)
        betadeltaP = betadeltaP +
     &                 Psi(i,j,k)*beta(k)*deltaP(k)

      enddo

!      if(iad.eq.0) then
!        fem_strain(i,j,:) = 0.0d0
!      endif

      ! Change gp_strain from incremental strain to total new strain
      tot_strain = conv_strain(i, j, :) + gp_strain

      ! tot_strain(1) = tot_strain(1) - alphadeltaT - betadeltaP
      ! tot_strain(2) = tot_strain(2) - alphadeltaT - betadeltaP
      ! tot_strain(3) = tot_strain(3) - alphadeltaT - betadeltaP

      bulk_mod = lambda_bar + (2.0d0/3.0d0)*G_bar
      vol_strain = tot_strain(1) + tot_strain(2) + tot_strain(3)
      ! vol_strain = vol_strain - 3.0d0*(alphadeltaT + betadeltaP)
      pressure = bulk_mod*(vol_strain -3.d0*(alphadeltaT + betadeltaP))

      dev_strain = tot_strain
      dev_strain(1) = dev_strain(1) - vol_strain/3.0d0
      dev_strain(2) = dev_strain(2) - vol_strain/3.0d0
      dev_strain(3) = dev_strain(3) - vol_strain/3.0d0
      dev_strain(4) = dev_strain(4)/2.0d0
      dev_strain(5) = dev_strain(5)/2.0d0
      dev_strain(6) = dev_strain(6)/2.0d0

      call J2Invariant(J2, 2.0d0*G_bar*dev_strain)
      Q = sqrt(3.0d0*J2)
      phi_trial = Q - sigma_y

      if((phi_trial.ge.phi_tol).and.(iUnload.eqv..false.)) then
        isPlastic(i, j) = 1
        deltaGamma = phi_trial/(3.0d0*G_bar)
        plastic_strain(i,j) = conv_pstrain(i,j) + deltaGamma

        !  phi_trial = Q - 3.0d0*G_bar*deltaGamma - sigma_y/sqrt(3.0d0)

        fac = 2.0d0*G_bar*(1.0d0 - 3.0d0*G_bar*deltaGamma/Q)
        gp_stress(1) = fac*dev_strain(1) + pressure
        gp_stress(2) = fac*dev_strain(2) + pressure
        gp_stress(3) = fac*dev_strain(3) + pressure
        gp_stress(4) = fac*dev_strain(4)
        gp_stress(5) = fac*dev_strain(5)
        gp_stress(6) = fac*dev_strain(6)

        fac = fac/(2.0d0*G_bar)
        tot_strain(1) = fac*dev_strain(1) + vol_strain/3.0d0
        tot_strain(2) = fac*dev_strain(2) + vol_strain/3.0d0
        tot_strain(3) = fac*dev_strain(3) + vol_strain/3.0d0
        tot_strain(4) = 2.0d0*fac*dev_strain(4)
        tot_strain(5) = 2.0d0*fac*dev_strain(5)
        tot_strain(6) = 2.0d0*fac*dev_strain(6)

      else
        fac = 2.0d0*G_bar
        gp_stress(1) = fac*dev_strain(1) + pressure
        gp_stress(2) = fac*dev_strain(2) + pressure
        gp_stress(3) = fac*dev_strain(3) + pressure
        gp_stress(4) = fac*dev_strain(4)
        gp_stress(5) = fac*dev_strain(5)
        gp_stress(6) = fac*dev_strain(6)
      endif

      fem_strain(i,j,:) = tot_strain
      fem_stress(i,j,:) = gp_stress

      end subroutine fem_vonMises_stress_update

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
subroutine sld_turnonGeneral(a)
   ! DESCRIPTION
   !    This routine opens the run for the elastic solid prob
   !-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_Solids
   implicit none
   class(SolidsProblem)  :: a
   integer(ip)           :: ndime,idime,count,icount,npoin,nod,sz
   real(rp), allocatable :: auxpoin(:,:)
   logical(lg)           :: pass

   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)

   !Material constants
   a%mu       = a%young/(2.0_rp + 2.0_rp*a%poisson)
   a%lambda   = a%poisson*a%young/((1.0_rp + a%poisson)*(1.0_rp - 2.0_rp*a%poisson))

   a%bulkMod  = a%young/(3.0_rp*(1.0_rp - 2.0_rp*a%poisson))

   !We modify the size of pwSource so that PETSC can recieve only
   !the nodes that will be modified
   if (allocated(a%pwSource)) then

       call a%Memor%alloc(a%ndofn+1,npoin,auxpoin,'auxpoin','sld_turnon')

       !Find non zero entries
       count=0
       do nod = 1,npoin

           pass =.false.
           if(a%pwSource(1,nod)/=0) then
               pass=.true.
           endif

           if (pass) then
               auxpoin(:,count+1) = a%pwSource(:,nod)
               count = count +1
           endif

       enddo

       a%pwSourceDim = a%ndofn + 1

       if(count > 0) then

           a%pwSourceSize=count

           call a%Memor%dealloc(a%ndofn+1,npoin,a%pwSource,'pwSource','sld_turnon')
           call a%Memor%dealloc(npoin,a%pwSourceId,'pwSourceId','sld_turnon')

           !Realloc with reduced entries
           call a%Memor%alloc(a%ndofn+1,count,a%pwSource,'pwSource','sld_turnon')
           call a%Memor%alloc(count,a%pwSourceId,'pwSourceId','sld_turnon')

           do icount = 1,count
               do idime= 1,a%ndofn+1
                   a%pwSource(idime,icount)=auxpoin(idime,icount)
               enddo
           enddo
       else
           a%pwSourceSize= npoin

       end if

       call a%Memor%dealloc(a%ndofn+1,npoin,auxpoin,'auxpoin','sld_turnon')

       call a%ModifyPointForces(a%ndofn,a%pwSourceSize,a%pwSource(2:a%pwSourceDim,1:a%pwSourceSize))

   end if

   call a%SolidSpecificTurnon

   if(a%sld_model .eq.'LINEA') then
       a%sld_type = 'LINEA'
   else

       a%sld_type = 'NONLI'
       a%kfl_linea= 2 !For now we default nonlinear sld to NR linearization

       !So we can move the mesh
       call a%Mesh%SetALE(1_ip)
       call a%Mesh%SetDisplacements(a%disp)
   endif

end subroutine sld_turnonGeneral

subroutine sld_turnonSpecific(a)
   ! DESCRIPTION
   !    This routine opens the run for the elastic solid prob
   !-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_Solids
   implicit none
   class(SolidsProblem) :: a
   integer(ip)          :: ndime,sz

   interface
      subroutine sld_initElasticTe(a,ndime,sz,young,nu,c)
         use typre
         use Mod_Solids
         implicit none
         class(SolidsProblem) :: a
         integer(ip)          :: ndime,sz
         real(rp)             :: young,nu
         real(rp)             :: c(sz,sz)
      end subroutine
   end interface

   call a%Mesh%GetNdime(ndime)
   sz = (ndime*(ndime+1))/2

   call sld_initElasticTe(a,ndime,sz,a%young,a%poisson,a%c_elastic)

end subroutine sld_turnonSpecific

subroutine sld_initElasticTe(a,ndime,sz,young,nu,c)
   use typre
   use Mod_Solids
   implicit none

   class(SolidsProblem) :: a
   integer(ip)          :: ndime,sz
   real(rp) :: c(sz,sz)
   real(rp) :: young  ,&           ! Young's modulus
                nu     ,&          ! Poisson
                lam    ,&
                mu     ,&
                aux    ,&
                aux_esf,&
                aux_dev

        !Auxiliary variables for esferic and deviatoric
        !components
        aux_esf=nu/(1.0_rp-nu)
        aux_dev=(1.0_rp-2.0_rp*nu)/(2.0_rp-2.0_rp*nu)
        aux=young*(1.0_rp-nu)/(1.0_rp-nu-2.0_rp*nu*nu)

        !Elastic tensor for isotropic materials
        if (ndime.eq.2) then

            c =  reshape([ 1.0_rp, aux_esf,  0.0_rp, &
                          aux_esf,  1.0_rp,  0.0_rp, &
                           0.0_rp,  0.0_rp, aux_dev], &
                          [3,3])

        elseif (ndime.eq.3) then

            c =  reshape([ 1.0_rp, aux_esf, aux_esf,  0.0_rp, 0.0_rp , 0.0_rp, &
                          aux_esf,  1.0_rp, aux_esf,  0.0_rp, 0.0_rp , 0.0_rp, &
                          aux_esf, aux_esf,  1.0_rp,  0.0_rp, 0.0_rp , 0.0_rp, &
                           0.0_rp,  0.0_rp,  0.0_rp, aux_dev,  0.0_rp, 0.0_rp, &
                           0.0_rp,  0.0_rp,  0.0_rp,  0.0_rp, aux_dev, 0.0_rp, &
                           0.0_rp,  0.0_rp,  0.0_rp,  0.0_rp,  0.0_rp, aux_dev], &
                          [6,6])

        endif

        c =  aux*c

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      subroutine mch_pulay(w_in,w_out,nint0,AA,nreset,islda)
******************************************
cc     Written by Lin-Wang Wang, March 30, 2001.
*************************************************************************
**  copyright (c) 2003, The Regents of the University of California,
**  through Lawrence Berkeley National Laboratory (subject to receipt of any
**  required approvals from the U.S. Dept. of Energy).  All rights reserved.
*************************************************************************

******************************************
ccccc  Now, this is in n1L,n2L,n3L

      use fft_data
      use load_data
      use data

      implicit double precision (a-h,o-z)
      include 'param.escan_real'

      real*8 w_in(mr_nL,islda),w_out(mr_nL,islda)

      real*8 AA(npulay_max,npulay_max)
      real*8, allocatable, dimension(:,:)  :: AA1
      real*8, allocatable, dimension(:)  :: B
      integer nreset

      allocate(AA1(npulay_max,npulay_max))
      allocate(B(npulay_max))

      alpha2=1.d0
cccc alpha2 controls how many recent charge densities
cccc to be used. If alpha2 > 1, then, the very old
cccc charge density is not used effectivly.
cccc We find that alpha2=1 is O.K.
******************************************************
      if(nint0.eq.1) nreset=0
      nint=nint0-nreset

      if(nint.gt.npulay_max) then
      write(6,*) "restart pulay, nint0,npulay_max",
     &    nint0,npulay_max
      nreset=nint0-1
      nint=1
      endif

      if(nint.eq.1) then

      do iislda=1,islda
      do i=1,nr_nL
      R0(i,iislda)=w_out(i,iislda)-w_in(i,iislda)
      w_in0(i,iislda)=w_in(i,iislda)
      enddo
      enddo

      endif
******************************************************
      if(nint.gt.1) then

      do iislda=1,islda
      do i=1,nr_nL
      dw(i,nint-1,iislda)=w_in(i,iislda)-w_in0(i,iislda)
      dR(i,nint-1,iislda)=w_out(i,iislda)-w_in(i,iislda)-
     &                 R0(i,iislda)
      R0(i,iislda)=w_out(i,iislda)-w_in(i,iislda)
      w_in0(i,iislda)=w_in(i,iislda)
      enddo
      enddo

      do m=1,nint-1
      s=0.d0

      do iislda=1,islda
      do i=1,nr_nL
      s=s+dR(i,m,iislda)*R0(i,iislda)
      enddo
      enddo

      call global_sumr(s)

      s=s*vol/nrL
      B(m)=-s
      enddo

******************************************************
      do m=1,nint-1
      s1=0.d0

      do iislda=1,islda
      do i=1,nr_nL
      s1=s1+dR(i,m,iislda)*dR(i,nint-1,iislda)
      enddo
      enddo

      call global_sumr(s1)

      s1=s1*vol/nrL
      AA(m,nint-1)=s1
      AA(nint-1,m)=s1
      enddo

cccccccccc pulay optimization
**********************************************************
      do m1=1,nint-1
      do m2=1,nint-1
      AA1(m1,m2)=AA(m1,m2)
      enddo
      enddo

      w=1.0d0
      do m=nint-1,1,-1
      AA1(m,m)=AA1(m,m)*w
      w=w*alpha2
      enddo

**********************************************************

      call gaussj(AA1,nint-1,npulay_max,B,1,1)

      do iislda=1,islda
      do i=1,nr_nL
      w_out(i,iislda)=R0(i,iislda)
      enddo
      enddo

      do m=1,nint-1
      do iislda=1,islda
      do i=1,nr_nL
      w_in(i,iislda)=w_in(i,iislda)+B(m)*dw(i,m,iislda)
      w_out(i,iislda)=w_out(i,iislda)+B(m)*dR(i,m,iislda)
      enddo
      enddo
      enddo

      do iislda=1,islda
      do i=1,nr_nL
      w_out(i,iislda)=w_out(i,iislda)+w_in(i,iislda)
      enddo
      enddo

      endif

      deallocate(AA1)
      deallocate(B)

      return
      end


text:
Name: Michelle
Age: 19
From: Stockton, CA
Occupation: UC Davis Student
Year: Sophomore
Miscellaneous Info: I follow the heart, not the anatomy. {{{^__}}}

About Me:

       Hey everyone! What it do? Im Michelle and its good to see that people are interested in me, or that they have nothing really to do and are just cruising around on Davis Wiki. {{{^__^}}} Well anyways, Im currently a sophomore at UC Davis majoring in Animal Biology. Hopefully that major could get me into vet school if I can manage to keep my grades up (Stupid chemistry..). Im a layedback kind of gal who enjoys having fun with friends and family. One of the best ways to do that is to play in the California Aggie Marching BandUh Cal Aggie Marching Banduh! and play rockin, kick ass tunes! Other than that I like doing the usual things of staying home and enjoying a relaxing day or going out there and play a quick game of basketball. Other than that you can check me out on Facebook and give me a message. See ya! {{{^__^}}}


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.01.09.02/06                        *
C  *    TEST TITLE : Error 408 handled by system           *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      CALL INITGL ('09.01.09.02/06')
      CALL E09408 (2)
      CALL ENDIT
      END


function velocity(x,y,z,ips)
common/flat_sf/key_flat1
vel_min=0.5

rz=6371

depth = sph_flat(key_flat1,x,y,z)

v0=vrefmod(depth,ips)

dv_apr= 0.01 * vert_anom(x,y,z,ips) * v0

dv = anom_3D_xyz_lin_v(x,y,z,ips)

velocity = v0 + dv + dv_apr

if(velocity .lt. vel_min) velocity = vel_min

return
end


text:
C
C $Id: mdloch.f,v 1.10 2008-09-18 00:42:17 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MDLOCH (RLON,CHRS,NCHR)
C
C This routine, given a longitude RLON and a character buffer CHRS,
C returns CHRS and NCHR such that CHRS(1:NCHR), when written by PLCHHQ,
C will yield the nautical representation of RLON, in degrees, minutes,
C and seconds east or west of the prime meridian at Greenwich.
C
        DOUBLE PRECISION RLON
        CHARACTER*(*)    CHRS
        INTEGER          NCHR
C
C Declare local variables.
C
        DOUBLE PRECISION ALON
        INTEGER          IDEG,IMIN,ISEC,NSEC
C
C Compute the equivalent value of the longitude in the range -180 to
C +180 degrees.
C
        ALON=MOD(MOD(RLON,360.D0)+360.D0,360.D0)
        IF (ALON.GE.180.D0) ALON=ALON-360.D0
C
C Compute the absolute value of the longitude in seconds, limiting it to
C the range from 0 to 180 degrees (0 to 648,000 seconds), and rounding
C to the nearest second.
C
        NSEC=MIN(648000,INT(3600.D0*ABS(ALON)+.5D0))
C
C Compute the degree, minute, and second quantities, as integers.
C
        IDEG=NSEC/3600
        IMIN=MOD(NSEC,3600)/60
        ISEC=MOD(NSEC,60)
C
C If the value is zero, return a string that will make PLCHHQ write a
C 0 followed by a degree sign.
C
        IF (IDEG.EQ.0.AND.IMIN.EQ.0.AND.ISEC.EQ.0) THEN
          NCHR=7
          CHRS='0:F34:0'
          RETURN
        END IF
C
C Otherwise, put blanks in the buffer and zero the current character
C count.
C
        CHRS=' '
        NCHR=0
C
C Put the number of degrees in the buffer, followed by the code for a
C degree sign.
C
        CALL MDINCH (IDEG,CHRS,NCHR)
        NCHR=NCHR+9
        CHRS(NCHR-8:NCHR)=':F34:0:F:'
C
C If the number of minutes and seconds are both zero, we're done except
C for, perhaps, a "W" or an "E".
C
        IF (IMIN.EQ.0.AND.ISEC.EQ.0) THEN
          IF (IDEG.EQ.180) RETURN
          GO TO 101
        END IF
C
C Put the number of minutes in the buffer, followed by the code for a
C single quote.
C
        CALL MDINCH (IMIN,CHRS,NCHR)
        NCHR=NCHR+15
        CHRS(NCHR-14:NCHR)=':F29H-5:'':H-4F:'
C
C If the number of seconds is zero, we're done except for, perhaps, a
C "W" or an "E".
C
        IF (ISEC.EQ.0) GO TO 101
C
C Put the number of seconds in the buffer, followed by the code for a
C double quote.
C
        CALL MDINCH (ISEC,CHRS,NCHR)
        NCHR=NCHR+15
        CHRS(NCHR-14:NCHR)=':F29H-5:":FH-4:'
C
C Finally, put either a "W" or an "E" in the buffer.
C
  101   NCHR=NCHR+1
C
        IF (ALON.LT.0.) THEN
          CHRS(NCHR:NCHR)='W'
        ELSE
          CHRS(NCHR:NCHR)='E'
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
\ essential "types"

struct: %color
    %color svar color.r
    %color svar color.g
    %color svar color.b
    %color svar color.a

struct: %rect
    %rect svar rect.x
    %rect svar rect.y
    %rect svar rect.w
    %rect svar rect.h

: x@    @ ;                       : x!    ! ;
: y@    cell+ @ ;                 : y!    cell+ ! ;
: w@    rect.w @ ;                : w!    rect.w ! ;
: h@    rect.h @ ;                : h!    rect.h ! ;
: xy@   2@ ;                      : xy!   2! ;
: wh@   rect.w 2@ ;               : wh!   rect.w 2! ;
: xywh@  4@ ;                     : xywh! 4! ;
: x2@   dup x@ swap w@ + ;        : x2!   >r r@ x@ - r> w! ;
: y2@   dup y@ swap h@ + ;        : y2!   >r r@ y@ - r> h! ;
: xy2@  dup 2@ rot wh@ 2+ ;       : xy2!  >r r@ xy@ 2- r> wh! ;

: aabb@ xywh@ 2over 2+ ;

create srcrect  0 , 0 , 0 , 0 ,

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine