C *************************************************************************
C Program CENMAG is a version of MAFIC for all datasets
C *************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*4 AR
      CHARACTER*6 FREDA
      DIMENSION LAR(3000),NAR(3000),FRAR(3000),EFRAR(3000),TCI(20),
     #     TCC(20),DEFC(20),CHI2(20),PM(20),EPM(20)
      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(20),GAMAG(20),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)
C------ parameters for p-mode calculations
      OPEN(UNIT=16,FILE='GBC',STATUS='UNKNOWN')
      OPEN(UNIT=17,FILE='EGBC',STATUS='UNKNOWN')
      OPEN(UNIT=18,FILE='CHIC',STATUS='UNKNOWN')
      WRITE(18,188) (XBOT(IV),IV=1,NIV)
188   FORMAT(10F8.5)      
      LCOW=20
      DSIM=5.0D-8
      NSE=0	
C p-mode data from Jesper 
      OPEN (UNIT=15,FILE='DATAF',STATUS='OLD')
999   CONTINUE
      NSE=NSE+1
      READ(15,*) NC,FREDA,Y,DY
      Y=1900+Y
      WRITE(*,*) FREDA, Y
      OPEN(2,FILE=FREDA,STATUS='OLD')
      IF(NSE.EQ.1) THEN
	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)
	GO TO 999
      ENDIF
      TOCI=0.
      TOCC=0.
      TII=0.
      TTW=0.
      DO IV=1,NIV
        TCI(IV)=0.
	TCC(IV)=0.
      ENDDO
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.16.and.NROR.eq.26) fres=4517.93
      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.23) FRES=4206.56
      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   
C      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
      TTW=TTW+1/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
C      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
      GAM=TOCI/TOCC
      FDF=TII-TOCI**2/TOCC
      CHIO2=FDF/(NPA-1)
      ERG=SQRT(TII/TOCC-GAM**2)
      DO IV=1,NIV       	
        FAC=1.
	DV=1-XBOT(IV)
	RDEPTH=DPHOT/DV
        IF(RDEPTH.LT.1.) FAC=FV(RDEPTH)
	TTT=TCI(IV)/TCC(IV)
        PM(IV)=100*FAC*TTT
	FDF=TII-TCI(IV)**2/TCC(IV)
	CHI2(IV)=FDF/(NPA-1)
	EPM(IV)=100*FAC*SQRT(TII/TCC(IV)-TTT**2)
      ENDDO
      CLOSE(2)
      WRITE(16,142)Y,GAM,(PM(IV),IV=1,NIV)
      WRITE(17,143)NSE,NPA,ERG,(EPM(IV),IV=1,NIV)
      WRITE(18,143)NSE,NPA,CHI0,(CHI2(IV),IV=1,NIV)	
142   FORMAT(15F8.3)
143   FORMAT(I3,I5,15F8.5)
      IF(NSE.LT.24) GO TO 999
      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),DT(5),DTP(5),UP(2,3),DUD(2,3),URK(2,5)
     #         ,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(20),GAMAG(20),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
         DO 67 J=1,3
         DO 67 I=1,2
67       URK(I,J)=UP(I,J)+FF*DUD(I,J)
         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=0
      IF(L.GT.0)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