C Last change: CA 21 Aug 2002 4:28 pm C ************************************************************************* C Date: 30.06.91. C KERNELS FOR RANDOM MAGNETIC (1kG) AND VELOCITY (km/s) 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='MAGN',STATUS='UNKNOWN') OPEN (UNIT=15,FILE='TURB',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 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) WRITE(16,141) L,SIGT,FRE*1.0E3,PE/60.,XTUR 141 FORMAT(I4,4F12.4) ENORM=4.0E8/DDD COB=1.06E-13*ENORM/FRE COV=3.76E-8*ENORM/FRE XPHOT=EXP(AR(1,NPHOT)) NFT=NOST DO KK=1,8 SUKER(KK)=0. ENDDO 111 NFT=NFT-1 X=FORK(1,NFT) DX=FORK(1,NFT+1)-X CO=COV DO KK=1,8 SUKER(KK)=SUKER(KK)+0.5*CO*(FORK(KK+1,NFT)+FORK(KK+1,NFT+1))*DX IF(KK.EQ.4) CO=COB ENDDO D=696*(XPHOT-X) WRITE(15,142)X,D,(COV*FORK(K,NFT),K=2,5) WRITE(14,142)X,D,(COB*FORK(K,NFT),K=6,9) WRITE(16,142)X,D,(SUKER(K),K=1,4) IF(X.GT.0.96) GO TO 111 142 FORMAT(F8.6,F8.2,1P4E12.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*CCC AA=0.0 IF (L.GT.0) AA=BH**2/PL DEK=(DRM**2+AA)*DEK ER=ER+DER*HI EK=EK+HI*DEK FORK(1,NOS)=EXP(ALX) XS=EXP(2*ALX) BH=BE*DHM Z=BH/PL HI=2.*H ALX=AD(1)+FF*D(1) 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 IF(EXP(ALX).GT.0.95) THEN c KERNELS FLAM=DDM FLAMS=FLAM**2 EC=DDEK Z=DHM/BC ZLA=Z*FLAM ZZ=PL*Z**2 DCS0=-GAMMA*(1+GAP+(1+GARO)/CHIR)*FLAMS+2*PL/CHIR*((1-CHIR)*EC- # (GAM+CHIR)*ZLA-ZZ) DCT=-GAMMA*(1+GARO)*FLAMS+PL*(EC-GAMMA*ZLA-ZZ) DCTC=DCT*CHIT/CHIR GDCDP=GAMMA*(1+GAP)-(VA/VG+1)*(1-GARO) GARD=0 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) c turbulent velocity DVV=-FLAMS-PL*(2*ZLA+EC) DVH=-PL*EC/2 RCV=DCS0+DVV RCH=DVH COV=RHOR*XSQ FORK(2,NOS)=COV*RCV FORK(3,NOS)=COV*RCH RCV=DCSV+DVV RCH=-DCSH/2+DVH FORK(4,NOS)=COV*RCV FORK(5,NOS)=COV*RCH C turbulent magnetic field DMV=2*PL*(ZLA+EC) DMH=2*FLAMS+PL*(4*ZLA+3*EC) RCV=DMV-DCS0 RCH=DMH+DCS0 FORK(6,NOS)=XSQ*RCV FORK(7,NOS)=XSQ*RCH RCV=DMV-DCSV+DCSH RCH=DMH+DCSV FORK(8,NOS)=XSQ*RCV FORK(9,NOS)=XSQ*RCH ENDIF 10 CONTINUE IF(N.EQ.NPHOT.AND.IC.EQ.1) DDD=DR**2 IF (N.LT.NT) GO TO 5 NOST=NOS RETURN END