C *************************************************************************
C Program MAGFIC-interprets frequency changes in terms of small-scale magnetic
C field changes. Current version is for vertical(radial) field change
C for horizontal replace BKR with BKH in line 473
C *************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*4 AR
      CHARACTER*6 SPLIDA,SPLIDAF,OUT1,OUT2 
      DIMENSION LAR(3000),NAR(3000),FRAR(3000),EFRAR(3000),TCI(30),
     #     TCC(30),DEFC(30)
      COMMON /GENER/ SM,SL,RT,AR(9,3000),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(3000),
     #       NPHOT,NOST,NT,L,NIC,LCOW,IC,NUMOD,NTOP
      COMMON/KER/NIV,XBOT(30),GAMAG(30),CHITAR(3000),CHIRAR(3000)
C    CONSTANTS
      PI=3.14159265D0
      G=6.672D-8
      SOM=1.989D33
      SOL=3.86D33
      SOR=6.96D10
C ....................................................
C Reading data on the reference model
C     the following two numbers are for SUNJCD      
      NPHOT=2401
      IKNI=1602
      OPEN (UNIT=1,FILE='SUNJCD.B',STATUS='OLD',FORM='UNFORMATTED')
      OPEN (UNIT=11,FILE='CHIJCD',STATUS='OLD',FORM='FORMATTED') 
       OPEN (UNIT=15,FILE='MAGFIC.DET',STATUS='UNKNOWN')
      READ (1) SM,RT,SL,XP,ZP,XC,NT
      DO 503 IK=1,NT
      READ (1) (AR(JK,IK),JK=1,9)
      IF(IK.GE.IKNI) READ(11,*) NIK,dep,x,ar1,CHITAR(IK),CHIRAR(IK)
  503 CONTINUE
      CLOSE (UNIT=1)
C
      WRITE (*,299)
      WRITE (*,402) SM,RT,SL,XP,ZP,XC,NT
      WRITE (*,299)
      WRITE (*,592) (AR(LK,1),LK=1,9)
      WRITE (*,592) (AR(LK,NT),LK=1,9)
      WRITE (*,299)
      ROM=0.75D0*SM*SOM/(PI*(RT*SOR)**3)
      DTS=DSQRT(1./(4.0*PI*G*ROM))
      SIM=30.0D0
      FN=48.
      NFI=0
      DO 1 N=2,NT
      A3N=AR(3,N)
      IF(AR(5,N).GT.1.0E-4) NTOP=N
C      IF(AR(6,N-1).LT.0.0.AND.AR(6,N).GT.0.0)NPHOT=N
      KD(N)=IDINT((AR(1,N)-AR(1,N-1))*SIM*DSQRT(AR(2,N)*A3N*
     #      AR(4,N))/PI*FN)+1
      NFI=NFI+KD(N)
    1 CONTINUE
C end model ......................................................

C-------magnetic field
      OPEN (UNIT=12,FILE='VERSIONS',STATUS='OLD',FORM='FORMATTED')      
      READ(12,*) NIV 
      READ(12,*) (XBOT(IV),IV=1,NIV)
      CLOSE(12)
C------ parameters for p-mode calculations
      LCOW=20
      DSIM=5.0D-8
      TOCI=0.
      TOCC=0.
      TII=0.
      DO IV=1,NIV
        TCI(IV)=0.
	TCC(IV)=0.
      ENDDO
      NUMOD=0	
C p-mode data from Jesper 
      OPEN (UNIT=17,FILE='DATAF',STATUS='OLD')
      WRITE(*,*) 'TYPE DATA SET NUMBERS: NF NL    '
      READ(*,*) NF, NL
      DO N=1,NL
        READ(17,*) NC,SPLIDA,Y,DY,OUT1,OUT2
        IF(N.EQ.NF) THEN
          OPEN(2,FILE=SPLIDA,STATUS='OLD')
	  KF=0
41        KF=KF+1
          READ(2,*)LAR(KF),NAR(KF),FRAR(KF),EFRAR(KF)
	  YF=1900+Y
	  IF(LAR(KF).LE.200) GO TO 41
	  KFT=KF
	  CLOSE(2)
	ENDIF
      ENDDO
      Y=1900+Y
      WRITE(*,*) SPLIDA, Y
      OPEN(2,FILE=SPLIDA,STATUS='OLD')
      OPEN(UNIT=16,FILE=OUT1,STATUS='UNKNOWN')
      OPEN(UNIT=19,FILE='FORPLOT',STATUS='UNKNOWN')
      OPEN(UNIT=18,FILE=OUT2,STATUS='UNKNOWN')
C.....................................................
C begin p-mode loop
      NUMOD=NUMOD+1
      NPA=0  
100   continue
      ID=0
      READ(2,*,END=900)L,NROR,FRES,EFRE
      IF(L.GT.200) GO TO 900	
      IF(NROR.EQ.0) GO TO 100
      DO KF=1,KFT
	IF(NROR.EQ.NAR(KF).AND.L.EQ.LAR(KF)) ID=KF
      ENDDO
      IF(ID.EQ.0) GO TO 100
      NPA=NPA+1
      DEF=FRES-FRAR(ID)
      EDEF=SQRT(EFRAR(ID)**2+EFRE**2)
      IF(L.EQ.14.and.NROR.eq.21) FRES=3723.16
      IF(L.EQ.14.and.NROR.eq.24) FRES=4150.57
      IF(L.EQ.17.and.NROR.eq.26) FRES=4552.25
      IF(L.EQ.18.and.NROR.eq.23) FRES=4160.92
      IF(L.EQ.19.and.NROR.eq.25) FRES=4494.7
      IF(L.EQ.23.and.NROR.eq.24) FRES=4489.08
876   CONTINUE      	
      SIGT=PI*DTS*FRES/5.0D5
      DSIS=0.1D0
      PL=L*(L+1)
      NIC=1
      IT=0
      IC=2
      SI=SIGT**2
      DSI=DSIS
C                   --------------------
      IF (L.GE.LCOW) CALL SETNIC
C                   --------------------
c      WRITE (*,*) ' NIC = ', NIC
   20 IT=IT+1
C-------------------------
      CALL CADPIN
C-------------------------
      W=AR(2,NT)*SI
      V=AR(3,NT)
      VG=V*AR(4,NT)
      RF1=PL/W
      W=VG*(4.+W)-3.*V
      DELT=(V-3)**2-4.*W
      DELTA=DELT-4.*(RF1*(V-VG-1.)-PL)
      RF2=0.5D0*(V-3.-DSQRT(DELTA))
      FLA=V-2.0*RF2
      RF2=(RF2-RF1+3.)/(VG-RF1)
      W=YA3(1)-YA1(1)*(1.0-RF2)
      IF(DABS(DSI).GT.DSIM) THEN
        IF(IT.GT.1) DSI=DSI/(W1/W-1.0)
	SI=SI+DSI
      SIGT=DSQRT(SI)
        W1=W
        IF(IT.LE.25) GO TO   20
        WRITE (*,111) L,FRES
  111 FORMAT (1X,'L,FRES :',I4,F10.5,' NEW FRES (<0 TO STOP)')
        READ(*,*) FRES
        IF(FRES.GT.0.)GO TO 876
        IF(FRES.LT.0.) STOP	
      ENDIF
      PE=2.0*PI*DTS/SIGT
      FRE=1.0D3/PE
C      WRITE (*,*) L,SIGT,FRE
      IC=1
      EK=0.
      ER=0.
      FRC=0.
      DO IV=1,NIV
	GAMAG(IV)=0
      ENDDO
C-----------------------------
      CALL CADPIN
C-----------------------------
      EKS=(1.0+PL*(YA3(1)/3.0D0/SI)**2)*AR(5,NT)/3.0D0/FLA
      EK0=EK
      EK=EK+EKS
      IF (L.GE.LCOW) FRC=-FRE*FRC/(4.D0*L+2.D0)/EK0*SI
      DD=YA3(1)
      ER=ER+AR(5,NT)*DD/9.0D0
      FRE=FRE+FRC
      FREQ=DSQRT(ER/EK0)*5.0D2/PI/DTS+FRC
      TEST=FREQ/FRE-1
      ACT=L*(L+1)/SIGT**2
      ACN=0.
      XTUR=0.
      N=2
      X=0.
877   N=N+1
      XP=X
      X=EXP(AR(1,N))
      ACP=ACN
      ACN=AR(2,N)*AR(3,N)*AR(4,N)
      IF(ACP.LT.ACN.AND.ACN.GE.ACT) THEN
        XTUR=XP+(X-XP)*(ACT-ACP)/(ACN-ACP)
	GO TO 880
      ENDIF
      GO TO 877
880   continue   
      WRITE (*,888) L,NROR, IT,XTUR
888   FORMAT (1X,'l, nror, it, XTUR  ',3i5,F8.5)
      FREMIC=FRE*1.0E3
c calculated DEF
      COB=7.4E-8/(EK*FREMIC)
      COV=2.6E-2/(EK*FREMIC)
      DPHOT=1.0-EXP(AR(1,NPHOT))
      TII=TII+(DEF/EDEF)**2
      BBB=DDD/(EK*4.0E8)
      TOCC=TOCC+(BBB/EDEF)**2
      TOCI=TOCI+BBB*DEF/EDEF**2
      DO IV=1,NIV
        BBB=COB*GAMAG(IV)
	DEFC(IV)=BBB
        TCC(IV)=TCC(IV)+(DEFC(IV)/EDEF)**2
	TCI(IV)=TCI(IV)+DEFC(IV)*DEF/EDEF**2
      ENDDO
      WRITE(15,141) L,NROR,NPA,FREMIC,XTUR,DEF,EDEF
141   FORMAT(2I3,I5,f10.3,F10.5,1p4e12.3)
      GO TO 100
900   CONTINUE
      PAR=TOCI/TOCC
      CHI2=(TII-TOCI**2/TOCC)/(NPA-1)
      WRITE(16,143)NPA,PAR,CHI2,TII/NPA 
143   FORMAT(I7,1P4E12.3)
      WRITE(19,144)Y,PAR,CHI2
      DO IV=1,NIV       	
        FAC=1.
	DV=1-XBOT(IV)
	RDEPTH=DPHOT/DV
        IF(RDEPTH.LT.1.) FAC=FV(RDEPTH)
        PAR=TCI(IV)/TCC(IV)
	CHI2=(TII-TCI(IV)**2/TCC(IV))/(NPA-1)
	WRITE(16,142)XBOT(IV),PAR,FAC*PAR,CHI2
	DBG=FAC*100*PAR
	DBMM=696*(DV-DPHOT)
	WRITE(19,144)DBMM,CHI2,DBG
      ENDDO	  
142   FORMAT(F7.5,1P3E12.3)
144   FORMAT(4f10.3)
      STOP
  299 FORMAT (1X)
  402 FORMAT (1X,6F12.7,I5)
  502 FORMAT (1X,1P9E13.6,0PF7.1)
  521 FORMAT (1X,2F10.6,1P1D24.15,I6)
  522 FORMAT (I5,1P6D24.15)
  523 FORMAT (1X,I4,1P6E12.4)
  592 FORMAT (1X,1P10E13.6)
  666 FORMAT (1X,2I5,0P1F12.4,I5,0P1F12.4)
      END
C *******************************************************************
      SUBROUTINE SETNIC
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      REAL*4 AR
      COMMON /GENER/ SM,SL,RT,AR(9,3000),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(3000),
     #       NPHOT,NOST,NT,L,NIC,LCOW,IC,NUMOD,NTOP
      RC1=10.0
      NA0=1 
      NA1=NT
      NA2=NT
      NIC=NT
      NS=0
      N=NT+1
      TAU=0.0
   10 N=N-1
      SSP=SS
      A3N=AR(3,N)
      SS=(L+0.5D0)**2-AR(2,N)*A3N*AR(4,N)*SI
      DX=AR(1,N+1)-AR(1,N)
      IF(SS.LT.0.) GO TO 10
      IF(TAU.GT.1.0D-8) GO TO 15
      D=DX/(1.-SSP/SS)
      XT=AR(1,N)+D
      RS=DSQRT(SS)
      TAU=D*RS/2.D0
      NA2=N
      GO TO 10
15    RSP=RS
      RS=DSQRT(SS)
      TAU=TAU+DX*(RS+RSP)/2.0D0
      A3N=AR(3,N)
      SC=RC1*0.5*DABS((RSP-RS)/DX/RS+AR(6,N)+1.0-A3N*AR(4,N))
      IF(TAU.LT.1.0.OR.RS.LT.SC) GO TO 10
      NIC=N
      RETURN
      END
c *************************************************************
      SUBROUTINE CADPIN
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*4 AR
      DIMENSION AD(9),D(9)
     #         ,APM(5),ETA(5),DETA(5),DDETA(5),OMEG(3)
      COMMON /GENER/ SM,SL,RT,AR(9,3000),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(3000),
     #       NPHOT,NOST,NT,L,NIC,LCOW,IC,NUMOD,NTOP
      COMMON/KER/NIV,XBOT(30),GAMAG(30),CHITAR(3000),CHIRAR(3000)
      ITYP=2
      F=1.
      IF (IC.EQ.1) F=1.0/YA1(1)
      IF (L.EQ.0.OR.L.GE.LCOW) ITYP=1
      IF (ITYP.EQ.1) THEN
         VG=AR(3,NIC)*AR(4,NIC)
         BC=AR(2,NIC)*SI
         IF (L.EQ.0) THEN
            GAM=1.D0/AR(4,NIC)
            YA1(1)=F*(1.0+VG*(3.*GAM-4.-BC)/10.D0)
            YA3(1)=-F*(3.0/VG-3.+1.5*GAM+0.5*BC)
         ELSE
            YA1(1)=F
            YA3(1)=F*BC/DSQRT(PL-VG*BC)
         ENDIF
      ELSE
         BB=4*L+6
         BA=AR(6,NIC)/BB
         BG=AR(3,NIC)*AR(4,NIC)/BB
         BC=AR(2,NIC)*SI
         BE=PL/BC
         YA1(2)=BG*(L+2)-BA*BE
         YA2(2)=1.0+0.6*(BA+BG)*(2*L+3)-3.0*BG
         YA3(2)=BG*BC-BA*(L+3)
         YA4(2)=L*YA2(2)-6.0*BG
         BB=BC/L
         YA1(1)=1.0+(1.0-BB)*YA1(2)
         YA2(1)=3.0*(BA+BB*BG)
         YA3(1)=BB+(1.0-BB)*YA3(2)
         YA4(1)=(L+2)*YA2(1)
         IF (IC.EQ.1) THEN
            RF1=F
            RF2=FSAVE*F
            YA1(1)=YA1(1)*RF1+YA1(2)*RF2
            YA2(1)=YA2(1)*RF1+YA2(2)*RF2
            YA3(1)=YA3(1)*RF1+YA3(2)*RF2
            YA4(1)=YA4(1)*RF1+YA4(2)*RF2
         ENDIF
      ENDIF
      N=NIC
    5 NM=N
      N=N+1
      NAK=KD(N)
      F=1.D0/NAK
      DO 6 J=1,8
      AD(J)=AR(J,NM)
6     D(J)=AR(J,N)-AD(J)
      H=0.5*F*D(1)
      DO 10 K=1,NAK
      FF=F*(K-0.5)
      BC=SI*(AD(2)+FF*D(2))
      VG=(AD(3)+FF*D(3))*(AD(4)+FF*D(4))
      U=AD(5)+FF*D(5)
      VA=AD(6)+FF*D(6)
      GAP=AD(7)+FF*D(7)
      GARO=AD(8)+FF*D(8)
      GAMMA=1.0/(AD(4)+FF*D(4))
      V=GAMMA*VG
      GARD=-D(8)/D(1)/V
      BE=PL/BC
      IF(IC.EQ.1) THEN
         CHIT=1.
         CHIR=1.
       IF(NM.GT.1602) THEN 
         CHIT=CHITAR(NM)+(CHITAR(N)-CHITAR(NM))*FF 
         CHIR=CHIRAR(NM)+(CHIRAR(N)-CHIRAR(NM))*FF
       ENDIF	
      ENDIF
      IF (ITYP.EQ.2.AND.N.EQ.NTOP) THEN
         ITYP=1
         IF (IC.EQ.2) THEN
            FSAVE=-(YA4(1)+(L+1)*YA2(1))/(YA4(2)+(L+1)*YA2(2))
            YA1(1)=YA1(1)+FSAVE*YA1(2)
            YA3(1)=YA3(1)+FSAVE*YA3(2)
         ENDIF
      ENDIF
      IF (ITYP.EQ.1) THEN
         AE=BC-VA
         IF (L.EQ.0) AE=AE+U
         S11=1.D0+H*(U-VA-1.D0) 
         S12=H*(BE-VG)
         S21=H*AE
         S22=1.+H*(3.-VG)
         DET=S11*S22-S12*S21
         DR=((2.D0*S11-DET)*YA1(1)+2.D0*S12*YA3(1))/DET
         DH=(2.D0*S21*YA1(1)+(2.D0*S22-DET)*YA3(1))/DET
         IF (IC.EQ.2) GO TO 69
         DRM=(DR+YA1(1))/2.D0
         DHM=(DH+YA3(1))/2.D0
	 F0M=0
	 F1M=0
         DDM=VG*(DRM-DHM)
   69    CONTINUE
         YA1(1)=DR
         YA3(1)=DH
      ELSE
         AL1=1.0D0/(1.0D0+H*(3.0D0-VG))
         AL2=1.0D0/(1.0D0+H*(U-1.0D0))
         P31=BC-VA
         P41=U*VA
         P32=-VA
         P42=PL-U*VG
         H11=AL1*(BE-VG)*H
         H12=VG*AL1*AL2*H**2
         H22=AL2*H
         S11=1.0D0+H*(U-P41*H12-P42*H22)
         S12=H*(P31*H12+P32*H22)
         S21=H*(U*VG+P41*H11)
         S22=1.0D0+H*(U-VA-1.0-P31*H11)
         DET=S11*S22-S12*S21
         BA=P31*VG*AL1*H+P32
         BB=P41*VG*AL1*H+P42
         H=2.0D0*H
         G11=H*AL1*(P31*S11+P41*S12)
         G12=H*AL2*(BA*S11+BB*S12)
         G21=H*AL1*(P31*S21+P41*S22)
         G22=H*AL2*(BA*S21+BB*S22)
         H=0.5D0*H
         DO 9 I=1,IC
         DH=((2.0D0*S11-DET)*YA3(I)+2.0D0*S12*YA4(I)+G11*YA1(I)+G12*
     #   YA2(I))/DET
         F1=((2.0D0*S22-DET)*YA4(I)+2.0D0*S21*YA3(I)+G21*YA1(I)+G22*
     #   YA2(I))/DET
         BA=DH+YA3(I)
         BB=F1+YA4(I)
         DR=(2.0D0*AL1-1.0D0)*YA1(I)+2.0D0*VG*AL1*AL2*YA2(I)*H+H11*BA+
     #   H12*BB
         F0=(2.0D0*AL2-1.0D0)*YA2(I)+H22*BB
         IF (IC.EQ.2) GO TO 89
         DRM=(DR+YA1(1))/2.0D0
         DHM=BA/2.D0
         F0M=(F0+YA2(1))/2.D0
         F1M=(F1+YA4(1))/2.D0
         DDM=VG*(DRM-DHM+F0M)
   89    CONTINUE
         YA1(I)=DR
         YA3(I)=DH
         YA2(I)=F0
         YA4(I)=F1
    9    CONTINUE
      ENDIF
    8 IF (IC.EQ.2) GO TO 10
      BH=BE*DHM
      Z=BH/PL
      HI=2.*H
      RHOUN=SI*U/BC
      ALX=AD(1)+FF*D(1)
      DEK=RHOUN*DEXP(5.0*ALX)
      DER=DDM**2/VG-2.*DDM*DRM+DRM**2*(VA+VG)
      IF (L.EQ.0) DER=DER-DRM**2*U
      IF (ITYP.EQ.2) DER=DER+DRM*F1M+BH*F0M
      DER=DER*DEK*SI/BC
      AA=0.0
      IF (L.GT.0) AA=BH**2/PL
      DDEK=DRM**2+AA
      DEK=DDEK*DEK
      ER=ER+DER*HI
      EK=EK+HI*DEK
      XS=EXP(2*ALX)
      X=EXP(ALX)
      XQ=EXP(3*ALX)
      IF(X.GT.0.71) THEN
      F1F=DDM**2
      F2F=PL*DDEK
      F3F=PL*DDM*Z
      F4F=2*PL*Z*DRM
c      F2F=0
c      F3F=0
c      F4F=0
      RHOP=1.0/CHIR
      CV=F1F*GAMMA*(1+RHOP*(1+GARO)+GAP)+2*F2F*(1-RHOP)+
     #   2*F3F*(RHOP*GAMMA)+2*PL**2*Z**2*RHOP
      BKR=4*F2F+2*F3F+CV
      BKH=2*F1F+3*F2F+4*F3F-CV
      VKR=(F1F+F2F+2*F3F+CV)*RHOUN
      VKH=0.5*F2F*RHOUN
      DEPTH=1-X
      DO IV=1,NIV
        FAC=1.
	DV=1-XBOT(IV)
	RDEPTH=DEPTH/DV
        IF(RDEPTH.LT.1.) FAC=FV(RDEPTH)
        GAMAG(IV)=GAMAG(IV)+FAC*BKR*XQ*HI
      ENDDO
      ENDIF
C        ----------------------------------------------------	
   10 CONTINUE
      IF(N.EQ.NPHOT.AND.IC.EQ.1) DDD=DR**2
C      IF(IC.EQ.1.AND.NUMOD.EQ.50)WRITE(15,144)N1,X,(ETA(K),K=1,5)
     
144   FORMAT(I3,1P7E11.2)      
      IF (N.LT.NT) GO TO 5
C       IF(IC.EQ.1)WRITE(15,144)NUMOD,F3,BOK,F4,DOK,YS,DDOK,DBD
      RETURN
      END       
C..........................................................
      DOUBLE PRECISION FUNCTION FV(X)
      IMPLICIT REAL*8(A-H,O-Z)
      FV=X**2*(3.-2.*X)
      RETURN
      END