cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Program
c Microscopic Global Optical Potential for nucleon-nucleus system
c with MPa interactrion based on RMF and HF+BCS
c (mgop2.f)
c
c  2017/02/07   version 1.00
c  2017/11/02   version 1.01
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit real*8(a-h,o-z)
      character*40 input,output
      parameter(mp=100,ig=12)
      dimension amp1(ig),amp2(ig),amp3(ig),amp4(ig),rgauss(ig)
      dimension yamp1(0:mp,ig),yamp2(0:mp,ig)
      dimension yamp3(0:mp,ig),yamp4(0:mp,ig)
c input
      open(1,file='mgopN.in')
      read(1,'(a40)')input
      read(1,'(a40)')output
      read(1,'(i6,f10.3)')inch,ee
      read(1,'(2f6.2)')amass,zmass
      read(1,'(2f6.2)')rmax,dr
      close(1)
c read input file
      open(1,file=input)
      read(1,*)iden,zmas2,minmas,maxmas
      idis=(maxmas-minmas)/2
      idnu=(nint(amass)-minmas)/2
      num=13
      ir=ig
      do 3001 m=1,idnu*num*ir
      read(1,*)
 3001 continue
      do 3002 m=0,num-1
      do 3002 i=1,ir
      if(inch.eq.1) then
      read(1,*)yamp1(m,i),yamp2(m,i),yamp3(m,i),yamp4(m,i)
      else
      read(1,*)d1,d2,d3,d4,yamp1(m,i),yamp2(m,i),yamp3(m,i),yamp4(m,i)
      end if
 3002 continue
      close(1)
c check input parameter
      if(inch.ne.1.and.inch.ne.0) then
      write(*,*)'Error! inch=',inch
      stop
      end if
      ene=ee
      if(ene.lt.50.d0.or.ene.gt.400.d0) then
      write(*,*)'Error! ee=',ene,'MeV/u'
      stop
      end if
      if(nint(zmass).ne.nint(zmas2)) then
      write(*,*)'Error! zmass=',zmass
      stop
      end if
      if(idnu.lt.0.or.idnu.gt.idis) then
      write(*,*)'Error! amass=',amass
      stop
      end if
c search parameter
      call search(num,ir,ene,amp1,mp,yamp1)
      call search(num,ir,ene,amp2,mp,yamp2)
      call search(num,ir,ene,amp3,mp,yamp3)
      call search(num,ir,ene,amp4,mp,yamp4)
      rgauss(12)=0.5d0*amass**(1./3.)+1.5d0
      do 3012 i=1,ir-1
      rgauss(i)=rgauss(12)/4.d0**(dble(12-i)/11.d0)
 3012 continue
c make potential & write
 1000 format('# Microscopic global optical potential based on RMF')
 1001 format('# Microscopic global optical potential based on HF+BCS')
 1002 format('# Incident proton energy  = ',f12.5,'(MeV)')
 1003 format('# Incident neutron energy = ',f12.5,'(MeV)')
 1004 format('# target        (A,Z,N) = (',f5.0,',',f5.0,',',f5.0,')')
 1005 format('# ')
 1006 format('# R(fm)    V_real(MeV)     W_imag(MeV)    Vso_real(MeV)   
     &  Wso_imag(MeV)')
 1007 format(f8.4,1p4d16.6)
      open(1,file=output)
      if(iden.eq.1) write(1,1000)
      if(iden.eq.2) write(1,1001)
      if(inch.eq.1) write(1,1002)ene
      if(inch.eq.0) write(1,1003)ene
      write(1,1004)amass,zmass,amass-zmass
      write(1,1005)
      write(1,1006)
      nrrmax=rmax/dr+1.0001
      do 100 i=1,nrrmax
      r=(i-1)*dr
      v=0.d0
      w=0.d0
      vls=0.d0
      wls=0.d0
      do 200 n=1,ir
      v  =v  +amp1(n)*exp(-r**2/rgauss(n)**2)
      w  =w  +amp2(n)*exp(-r**2/rgauss(n)**2)
      vls=vls+amp3(n)*exp(-r**2/rgauss(n)**2)
      wls=wls+amp4(n)*exp(-r**2/rgauss(n)**2)
 200  continue
      write(1,1007)r,v,w,vls,wls
 100  continue
      close(1)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine search(num,ir,ene,amp,mp,yamp)
      implicit real*8(a-h,o-z)
      dimension amp(ir)
      dimension X(0:MP),Y(0:MP),yamp(0:mp,ir)
      x( 0)= 50.d0
      x( 1)= 60.d0
      x( 2)= 70.d0
      x( 3)= 80.d0
      x( 4)= 90.d0
      x( 5)=100.d0
      x( 6)=120.d0
      x( 7)=140.d0
      x( 8)=160.d0
      x( 9)=180.d0
      x(10)=200.d0
      x(11)=300.d0
      x(12)=400.d0
c energy
      do 1003 n=1,ir
      do 1004 i=0,num-1
      y(i)=yamp(i,n)
 1004 continue
      call spline(num,ene,yy,x,y,mp)
      amp(n)=yy
 1003 continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine spline(N2,xx,yy,x,y,mp)
      implicit real*8(a-h,o-z)
c
c     N2:number of points(input)
c     xx:your wish point (input)
c     yy:your wish point (output)
c     x :x component(input)
c     y :y component(input)
c     mp:dimension of x and y components(input)
C
C     SPLINE METHOD
C
      DIMENSION X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
C
      N=N2-2
      CALL SPL21(X,Y,G,H,N,MP)
C     DO 40 I=0,N+1
C  40 WRITE(6,*) X(I),Y(I),G(I),H(I)
      YY=SPL22(G,H,X,Y,XX,N,MP)
      return
      END
C
C     CALCULATION OF F(X) WITH CUBIC SPLINE METHOD
C
      FUNCTION SPL22(G,H,X,Y,XX,N,MP)
      implicit real*8(a-h,o-z)
      DIMENSION X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
      DO 10 I=1,N+1
      IF(X(I-1).LE.XX.AND.XX.LE.X(I)) GO TO 20
   10 CONTINUE
   20 XI=XX-X(I-1)
      XL=X(I)-X(I-1)
      SPL22=Y(I-1)+XI*(G(I-1)
     1     +0.5*XI*(H(I-1)+XI*(H(I)-H(I-1))/(3.0*XL)))
      RETURN
      END
C
C     CALCULATION OF THE COEFFICIENTS
C
      SUBROUTINE SPL21(X,Y,G,H,N,MP)
      implicit real*8(a-h,o-z)
      PARAMETER (NP=100)
      DIMENSION A(NP,NP+1),X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
      N1=N+1
      DO 10 I=1,N
      DO 10 J=1,N
   10 A(J,I)=0.0
      DO 11 I=1,N
      XL=X(I)-X(I-1)
      XR=X(I+1)-X(I)
   11 A(I,N1)=6.0*(Y(I-1)/XL-(XL+XR)/XL/XR*Y(I)+Y(I+1)/XR)
      A(1,1)=2.0*(X(2)-X(0))
      A(1,2)=X(2)-X(1)
      A(N,N)=2.0*(X(N1)-X(N-1))
      A(N,N-1)=X(N)-X(N-1)
      DO 12 I=2,N-1
      A(I,I-1)=X(I)-X(I-1)
      A(I,I)=2.0*(X(I+1)-X(I-1))
   12 A(I,I+1)=X(I+1)-X(I)
      CALL GASJO(A,DET,N,NP)
      DO 13 I=1,N
   13 H(I)=A(I,N1)
      H(0)=0.0
      H(N1)=0.0
      XL=X(1)-X(0)
      G(0)=(Y(1)-Y(0))/XL-XL*H(1)/6.0
      DO 14 I=1,N1
   14 G(I)=G(I-1)+0.5*(H(I)+H(I-1))*(X(I)-X(I-1))
      RETURN
      END
C
C     GAUSS-JORDAN METHOD
C
      SUBROUTINE GASJO(A,DET,N,NP)
      implicit real*8(a-h,o-z)
C
      PARAMETER (MP=100)
      DIMENSION A(NP,NP+1),IP(MP),IN(MP,2)
C
      IF(MP.LT.NP) THEN
        WRITE(6,*) 'MP should be, MP.GE.NP'
        STOP
        ENDIF
C
      N1=N+1
      DET=1.0
      DO 5 J=1,N
    5 IP(J)=0
      DO 200 I=1,N
      T=0.0
      DO 60 J=1,N
      IF(IP(J).EQ.1) GO TO 60
      DO 59 K=1,N
      IF(IP(K).EQ.1) GO TO 59
      TTX=ABS(A(J,K))
      IF(T.GT.TTX) GO TO 59
      IR=J
      IC=K
      T=TTX
   59 CONTINUE
   60 CONTINUE
      IF(I.EQ.1) TMAX=T
      IP(IC)=1
      IF(IR.EQ.IC) GO TO 7
      DET=-DET
      DO 8 L=1,N1
      T=A(IR,L)
      A(IR,L)=A(IC,L)
    8 A(IC,L)=T
    7 IN(I,1)=IR
      IN(I,2)=IC
      PI=A(IC,IC)
      DET=DET*PI
      IF(ABS(PI)/TMAX.LT.1.E-6) GO TO 111
      A(IC,IC)=1.0
      DO 9 L=1,N1
    9 A(IC,L)=A(IC,L)/PI
      DO 100 LI=1,N
      IF(LI.EQ.IC) GO TO 100
      T=A(LI,IC)
      A(LI,IC)=0.0
      DO 11 L=1,N1
   11 A(LI,L)=A(LI,L)-A(IC,L)*T
  100 CONTINUE
  200 CONTINUE
      DO 13 I=1,N
      L=N-I+1
      IF(IN(L,1).EQ.IN(L,2))  GO TO 13
      JR=IN(L,1)
      JC=IN(L,2)
      DO 12 K=1,N
      T=A(K,JR)
      A(K,JR)=A(K,JC)
      A(K,JC)=T
   12 CONTINUE
   13 CONTINUE
      GO TO 9999
  111 WRITE(6,2000)
 2000 FORMAT(1H,'*** MATRIX A IS SINGULAR ***')
 9999 RETURN
      END      
