C     Last change:  CA   31 Aug 2002    6:47 pm
C *************************************************************************
C                                                      Date: 30.06.91.
C KERNELS FOR TOROIDAL MAGNETIC (1MG)  FIELDS
C *************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (MGRID=2500, MFINE=6000)
      REAL*4 AR
      DIMENSION ARQ(10),SUKER(4)
      COMMON /GENER/ SM,SL,RT,AR(9,MGRID),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(MGRID),
     #       NPHOT,NOST,NT,L,NIC,LCOW,IC,NUMOD,NTOP
      COMMON/KER/FORK(9,MFINE),CHITAR(3000),CHIRAR(3000)
C ....................................................
C Reading array AR(1:9,1:NT)
C
C     the following two numbers are for SUNJCD      
      NPHOT=2401
      IKNI=1602
C      WRITE (*,*) ' Type name of binary file for the reference model:'
C      READ (*,*) FIN1
      OPEN (UNIT=1,FILE='SUNJCD.B',STATUS='OLD',FORM='UNFORMATTED')
      OPEN (UNIT=11,FILE='CHIJCD',STATUS='OLD',FORM='FORMATTED') 
      OPEN (UNIT=14,FILE='KERTOR',STATUS='UNKNOWN')
      OPEN (UNIT=16,FILE='CONTROL',STATUS='UNKNOWN')
      REWIND (1)
      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)
C end of reading array AR
C .......................................................
      PI=3.14159265D0
      G=6.672D-8
      SOM=1.989D33
      SOL=3.86D33
      SOR=6.96D10
      ROM=0.75D0*SM*SOM/(PI*(RT*SOR)**3)
      DTS=DSQRT(1./(4.0*PI*G*ROM))
C
       LCOW=20
      DSIM=1.0D-5
C.....................................................
        WRITE (*,*) 'L,FRE [micHz]'
        READ (*,*) L,FRES
C
      SIGT=PI*DTS*FRES/5.0D5
      NTOP=NT
      FN=50.
      NTT=0
      ACT=L*(L+1)/SIGT**2
      ACN=0.
      XTUR=0.
      JAC=0
      DO 1 N=2,NT
      XP=X
      X=EXP(AR(1,N))
      IF(AR(5,N).GT.1.0E-4) NTOP=N
      ACP=ACN
      ACN=AR(2,N)*AR(3,N)*AR(4,N)
      IF(ACP.LT.ACN.AND.ACN.GE.ACT.AND.JAC.EQ.0.AND.L.GT.0) THEN
        XTUR=XP+(X-XP)*(ACT-ACP)/(ACN-ACP)
	JAC=1
      ENDIF
      KAC=IDINT((AR(1,N)-AR(1,N-1))*SIGT*DSQRT(ACN)/PI*FN)+1
      KGR=1
      BV=AR(6,N)/AR(2,N)
      IF(BV.GT.0.7*SIGT)
     # KGR=IDINT((AR(1,N)-AR(1,N-1))*DSQRT(L*(L+1)*BV)/SIGT/X/PI*FN)+1
      NAD=MAX0(KAC,KGR)
      KD(N)=NAD
      NTT=NTT+NAD
    1 CONTINUE
      WRITE(*,*) '   NTT  =',NTT
      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                   --------------------
      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 (IT.EQ.1) GO TO         40
      DSI=DSI/(W1/W-1.0)
      IF (DABS(DSI).GT.DSIM) GO TO            40
C
      PE=2.0*PI*DTS/SIGT
      FRE=1.0D3/PE
      WRITE (*,*) L,SIGT,FRE
      IC=1
      EK=0.
      ER=0.
      FRC=0.
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
      WRITE (*,888) FREQ,ER,FRC,EK,FRE
  888 FORMAT (1X,'FREQ,ER,FRC,EK,FRE: ',F8.5,1P3E11.3,0PF8.5)
      FREMI=FRE*1.0E3
      WRITE(16,141) L,SIGT,FREMI,PE/60.,XTUR
141   FORMAT(I4,4F12.4)
      ENORM=4.0E8/DDD
      COT=500/FREMI*ENORM
      COM=0.148/FREMI*ENORM
      XPHOT=EXP(AR(1,NPHOT))
      DO  N=3,NT
      XD=EXP(AR(1,N))/xphot
      WRITE(16,142)XD,COT*FORK(2,N),(COM*FORK(K,N),K=3,5),FORK(6,N)/EK
      IF(xd.gt.0.65.and.xd.lt.0.75)
     #WRITE(14,142)XD,COT*FORK(2,N),(COM*FORK(K,N),K=3,5)
      ENDDO
142   FORMAT(F8.6,1P5E12.2)
      STOP
C
   40 SI=SI+DSI
      W1=W
      SIGT=DSQRT(SI)
C--------------------------
      IF(IT.LE.30) GO TO   20
      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)
  777 FORMAT (2I4,F10.3,F8.3)
      END
C *******************************************************************
      SUBROUTINE SETNIC
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MGRID=2500)
      REAL*4 AR
      COMMON /GENER/ SM,SL,RT,AR(9,MGRID),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(MGRID),
     #       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)
      PARAMETER (MGRID=2500, MFINE=6000)
      REAL*4 AR
      DIMENSION AD(9),D(9)
      COMMON /GENER/ SM,SL,RT,AR(9,MGRID),SI,SIGT,RF1,RF2,PL,EK,ER,E9,
     #       DDD,YA1(2),YA2(2),YA3(2),YA4(2),KD(MGRID),
     #       NPHOT,NOST,NT,L,NIC,LCOW,IC,NUMOD,NTOP
      COMMON/KER/FORK(9,MFINE),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
      IF (IC.EQ.1) NOS=NIC
      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
         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
      NOS=NOS+1
      BH=BE*DHM
      HI=2.*H
      ALX=AD(1)+FF*D(1)
      CCC=SI/BC
      RHOR=U/CCC
      DEK=RHOR*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
C   KERNELS
      FLAM=DDM
      FLAMS=FLAM**2
      EC=DDEK
      Z=DHM/BC
      ZLA=Z*FLAM
      ZZ=PL*Z**2
      DCS0A=-(GAMMA*(1+GAP)+1+GARO)*FLAMS-2*PL*((1-1/GAMMA)*EC+
     #    2*ZLA+ZZ/GAMMA)
      DCT=-GAMMA*(1+GARO)*FLAMS+PL*(EC-GAMMA*ZLA-ZZ)
      DCTC=DCT*CHIT/CHIR
      GDCDP=GAMMA*(1+GAP)-(VA/VG+1)*(1-GARO)
      PSIGAR=GARO*((2+GDCDP)*FLAMS+2*GAMMA*FLAM*(PL*Z-BC*DRM))+GARD*
     #       FLAMS
      DCSV=-2*BC*DRM*(FLAM+2*PL*Z)+PSIGAR
      DCSH=-DCSV-GAMMA*(1+GAP)*FLAMS-2*PL*(EC+ZLA)
      X=EXP(ALX)
      XSQ=EXP(2*ALX)
      RCT0=(DCS0A+2*FLAMS+3*PL*EC)/4.
      RCT1=(DCSV-2*DCSH+2*FLAMS-3*EC+(6*PL+6)*ZLA)*5/28.
      RCT2=-(6*(2*DCSV+3*DCSH)+24*FLAMS+(56*PL-120)*EC-40*(PL-6)*ZLA)
     #    /28.
      FORK(2,NOS)=3*DCT*XSQ**4*RHOR/CCC/V
      FORK(3,NOS)=RCT0*XSQ
      FORK(4,NOS)=RCT1*XSQ
      FORK(5,NOS)=RCT2*XSQ
      FORK(6,NOS)=DEK*X
   10 CONTINUE
      IF(N.EQ.NPHOT.AND.IC.EQ.1) DDD=DR**2
      IF (N.LT.NT) GO TO 5
      NOST=NOS
      RETURN
      END