C ************************************************************************* C Program MAGFIS-interprets even-a coefficients in terms of small-scale C magnetic field. Current version is for vertical(radial) field change C for horizontal see line 269 and replace BKR with BKH in line 530 C ************************************************************************* IMPLICIT REAL*8(A-H,O-Z) REAL*4 AR CHARACTER*6 SPLIDA, OUT1,OUT2 DIMENSION NMODE(18),AEI(18),EAE(18),TII(18),AEC(18,20),TCI(18,20), # TOCC(18),TOCI(18),TCC(18,20),PAR(3),PARPH(3),CHI2(3),AEROT(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) COMMON/ROT/NBIN,XCO,ET(5,30),OMEGA(3,30),TI(5,5),UR(2,3,3000) # ,FMO(3),AP(5),AQ(5),AO(3) C CONSTANTS PI=3.14159265D0 G=6.672D-8 SOM=1.989D33 SOL=3.86D33 SOR=6.96D10 DO 52 J=1,5 TT=1.0/(2.0*J+1.0) DO 51 K=1,J TT=TT*(J-K+1.0)/(J+K+0.5) TI(J,K)=(4.0*K+1.)*TT 51 CONTINUE 52 CONTINUE 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=3,FILE='OM',STATUS='OLD',FORM='FORMATTED') 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 centrifugal force READ(3,*)NOM,XCO READ(3,*) ((OMEGA(I,N),I=1,3),N=1,NOM) OMS=OMEGA(1,NOM) OMSS=OMS**2 EPS=3.0*(OMS*2.0E-9*PI*DTS)**2 DO 42 N=1,NOM ET(1,N)=OMEGA(1,N)**2/OMSS ET(2,N)=2.0*OMEGA(2,N)*OMEGA(1,N)/OMSS ET(3,N)=(OMEGA(2,N)**2+2.0*OMEGA(3,N)*OMEGA(1,N))/OMSS ET(4,N)=2.0*OMEGA(3,N)*OMEGA(2,N)/OMSS 42 ET(5,N)=OMEGA(3,N)**2/OMSS NBIN=NOM-1 CALL CLAR OPEN (UNIT=15,FILE='MAGFIS.DET',STATUS='UNKNOWN') WRITE(*,434) (EPS*FMO(K),K=1,3) WRITE(15,434) (EPS*FMO(K),K=1,3) WRITE(*,435) (-EPS*UR(1,K,NT),K=1,3) WRITE(15,435) (-EPS*UR(1,K,NT),K=1,3) 434 FORMAT(2X,' Js : ',1P5E13.3) 435 FORMAT(1X,'DR/R : ',1P5E13.3) 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 DO K=1,18 NMODE(K)=0 TII(K)=0. TOCI(K)=0 TOCC(K)=0 DO IV=1,NIV TCI(K,IV)=0. TCC(K,IV)=0. ENDDO ENDDO NUMOD=0 C p-mode data from Jesper OPEN (UNIT=17,FILE='DATMS',STATUS='OLD') WRITE(*,*) 'TYPE DATA SET NUMBER ' READ(*,*) NSE DO N=1,NSE READ(17,*) NC,SPLIDA,Y,DY,OUT1,OUT2 ENDDO Y=1900+Y WRITE(*,*) SPLIDA, Y OPEN(2,FILE=SPLIDA,STATUS='OLD') OPEN(UNIT=16,FILE=OUT1,STATUS='UNKNOWN') OPEN(UNIT=18,FILE=OUT2,STATUS='UNKNOWN') OPEN(UNIT=19,FILE='FORPLOT',STATUS='UNKNOWN') C..................................................... C begin p-mode loop NUMOD=NUMOD+1 100 continue READ(2,*,end=900)l,nror,fres,(aei(k),k=1,18),(eae(k),k=1,18) 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 if(l.gt.200) go to 900 if(nror.eq.0) go to 100 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 calculate a_2k COB=7.4E-8/(EK*FREMIC) COV=2.6E-2/(EK*FREMIC) DPHOT=1.0-EXP(AR(1,NPHOT)) KLAST=18 IF(L.LT.KLAST) KLAST=L FAMP=1.0/L DO K=1,KLAST NMODE(K)=NMODE(K)+1 AEI(K)=AEI(K)*1.0E-3 EE=EAE(K)*1.0E-3 IF(EE.LT.2.0E-5) EE=1.0E10 FAMP=-FAMP*(2.0-1.0/K)*(L+1.0-K)/(2.0*(L+K)+1.0) IF(K.LE.3) THEN AEROT(K)=AP(K)*FAMP*OMSS*1.0E-6/(2*FREMIC*EK) AEI(K)=AEI(K)-AEROT(K) ENDIF TII(K)=TII(K)+(AEI(K)/EE)**2 BBB=FAMP*DDD/(EK*4.0E8) TOCC(K)=TOCC(K)+(BBB/EE)**2 TOCI(K)=TOCI(K)+BBB*AEI(K)/EE**2 DO IV=1,NIV BBB=COB*GAMAG(IV) AEC(K,IV)=FAMP*BBB TCC(K,IV)=TCC(K,IV)+(AEC(K,IV)/EE)**2 TCI(K,IV)=TCI(K,IV)+AEC(K,IV)*AEI(K)/EE**2 ENDDO ENDDO WRITE(15,141) L,NROR,IT,FREMIC,XTUR,EK/DDD,(AEROT(K),K=1,3) 141 FORMAT(3I3,f10.3,F10.5,1p4e12.3) C IF(NUMOD.EQ.50) STOP GO TO 100 900 CONTINUE ZC=1.48E-11 RHO=AR(5,NPHOT)/AR(2,NPHOT) V=AR(3,NPHOT) VRHO=V*AR(4,NPHOT)+AR(6,NPHOT) CHIT=CHITAR(NPHOT) CHIR=CHIRAR(NPHOT) GAMMA1=TOCI(1)/TOCC(1) CHIO2=(TII(1)-TOCI(1)**2/TOCC(1))/(NMODE(1)-1) WRITE(19,144)Y,GAMMA1,CHIO2 DO IV=1,NIV FAC=1. DV=1-XBOT(IV) RDEPTH=DPHOT/DV DFAC=0. IF(RDEPTH.LT.1.) THEN FAC=FV(RDEPTH) DFAC=-4*(1-RDEPTH)/(RDEPTH*(3-2*RDEPTH)*DV) ENDIF DO K=1,3 PAR(K)=TCI(K,IV)/TCC(K,IV) CHI2(K)=(TII(K)-TCI(K,IV)**2/TCC(K,IV))/(NMODE(K)-1) ENDDO WRITE(16,142)XBOT(IV),(PAR(K),FAC*PAR(K),CHI2(K),K=1,3) Z=ZC*FAC*PAR(1)/RHO DLP=V*Z DLRHO=2*Z*(DFAC+VRHO) C for horizontal field use C DLP=-V*Z C DLRHO=-Z*(DFAC+VRHO) DLT=(DLP-CHIR*DLRHO)/CHIT WRITE(18,143)XBOT(IV),DLP,DLRHO,DLT DBG=-0.5*FAC*100*PAR(1) DTG=-0.5*DLT*5873 DBMM=696*(DV-DPHOT) WRITE(19,144)DBMM,CHI2(1),DBG,DTG ENDDO 142 FORMAT(F7.5,9F8.3) 143 FORMAT(f7.5,1P6E12.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),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) COMMON/ROT/NBIN,XCO,ET(5,30),OMEGA(3,30),TI(5,5),UR(2,3,3000) # ,FMO(3),AP(5),AQ(5),AO(3) XSTEP=(1.0-XCO)/NBIN 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) THEN F=F/YA1(1) DO 60 J=1,5 AP(J)=0.0 AQ(J)=0.0 IF(J.LT.4) THEN AO(J)=0. OMEG(J)=OMEGA(J,1) ENDIF ETA(J)=ET(J,1) DETA(J)=0.0 60 CONTINUE ENDIF N=NIC 5 K=N N=N+1 NAK=KD(N) F=1.D0/NAK DO 6 J=1,8 AD(J)=AR(J,K) 6 D(J)=AR(J,N)-AD(J) IF(IC.EQ.1) THEN DO 65 J=1,3 DO 65 I=1,2 UP(I,J)=UR(I,J,K) 65 DUD(I,J)=UR(I,J,N)-UP(I,J) ENDIF 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) 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 DG1=1+GAP DG2=1+GARO DGDLP=GAMMA*(GAP+GARO*(VA+VG)/V) DG3=GARD+GARO*(3-GAMMA+VA/VG-DGDLP) C DG1=1 C DG2=1 C DG3=0 CHV=2*BC*(DRM*DDM*DG2+F4F)-F1F*DG3-2*F3F*GARO CH=F1F*DG2+2*(F2F+F3F) BKR=4*F2F+2*F3F+2*CHV-CH BKH=2*F1F+3*F2F+4*F3F-CHV VKR=(F1F+F2F+2*F3F+CHV)*RHOUN VKH=0.5*(F2F+CHV+CH)*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 CALCULATE QUADRATIC EFFECT OF ROTATION FX=(X-XCO)/XSTEP N1=FX+0.5 IF(N1.EQ.0) THEN DO 25 J=1,5 CET=ET(J,2)-ET(J,1) ETA(J)=ET(J,1)+CET*FX**2 IF(J.LT.4) OMEG(J)=OMEGA(J,1)+(OMEGA(J,2)-OMEGA(J,1))*FX**2 DETA(J)=2.*FX*X/XSTEP*CET 25 DDETA(J)=2.*(X/XSTEP)**2*CET ENDIF IF(N1.GT.0) THEN IF(N1.GT.NBIN-1) N1=NBIN-1 N2=N1+1 N3=N1+2 DX=FX+0.5-N1 DO 27 J=1,5 IF(J.LT.4) OMEG(J)=0.125*(3.*OMEGA(J,N1)+6.*OMEGA(J,N2)- # OMEGA(J,N3))+DX*(OMEGA(J,N2)-OMEGA(J,N1))+ # 0.5*DX**2*(OMEGA(J,N1)-2.0*OMEGA(J,N2)+OMEGA(J,N3)) CET=ET(J,N2)-ET(J,N1) CETA=ET(J,N1)-2.*ET(J,N2)+ET(J,N3) ETA(J)=0.125*(3.*ET(J,N1)+6.*ET(J,N2)-ET(J,N3))+DX*CET+ # 0.5*DX**2*CETA DETA(J)=X/XSTEP*(CET+DX*CETA) DDETA(J)=(X/XSTEP)**2*CETA 27 CONTINUE ENDIF A8=AR(8,N-1)+FF*(AR(8,N)-AR(8,N-1)) ZS=Z**2 YZ=DRM*Z YS=DRM**2 F10=2.0*BC*(YS+(2.-U)*YZ)-DRM*(3.0*F1M+(4.-2.*U)*F0M) F11=BC*ZS/2.0 F3=DRM*(BC*Z-F0M) F20=2.0*(BC*(YS*(2.-U)+PL*ZS)+F3*((2.-U)**2-2.)+ # U*(DDM*DRM-YS*(VA+VG))+DRM*(F1M*(1.5*U-2.)-PL*F0M)+PL*Z*F1M) F21=ZS*(1.-U/2.)+F3 F40=-BC*(YS+PL*ZS)+DDM**2*A8/VG+YS*(U-3.)+2.*(DRM*F1M+PL*Z*F0M) F41=0.5*(BC*ZS-YZ+2.*Z*F0M) HH=U*X**5 HC=HH*SI/BC DO 31 IS=1,5 LAS=2.*IS*(2.*IS+1) DOK=0. BOK=0. DDOK=0. UH=0. DUH=0. DO 30 J=IS,5 TU=TI(J,IS) DU=0.5*TU/J IF(IS.GT.3) THEN UH=UH+DU*ETA(J)*XS DUH=DUH+DU*(2.0*ETA(J)+DETA(J))*XS ENDIF DOK=DOK+TU*(0.5*DETA(J)-(1.-1.0/J)*ETA(J)) DDOK=DDOK+TU*(0.5*DETA(J)+0.5*DDETA(J)-(1.-1.0/J)*DETA(J)) IF(J.LT.5) THEN DOK=DOK+TU*ETA(J+1) DDOK=DDOK+TU*DETA(J+1) ENDIF BOK=BOK+DU*(DDETA(J)+6.0*DETA(J)+(6.-LAS)*ETA(J)) 30 CONTINUE BOK=BOK+U*DOK IF(IS.GT.3) THEN URK(1,IS)=UH URK(2,IS)=DUH ENDIF F1=F10+LAS*F11 F2=F20+LAS*F21 F4=F40+LAS*F41 DBD=F3*BOK+F4*DOK-YS*DDOK DT(IS)=((URK(2,IS)*F1+URK(1,IS)*F2)/XS+DBD)*HC AP(IS)=AP(IS)+DT(IS)*H*2.0 AQ(IS)=AQ(IS)+ETA(IS)*DEK IF(IS.LT.4) AO(IS)=AO(IS)+OMEG(IS)*DEK 31 CONTINUE 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.............................................................. C SUBROUTINE CLAR IMPLICIT REAL*8(A-H,O-Z) REAL*4 AR DIMENSION F(3),FD(3),ETA(5),DETA(5) 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/ROT/NBIN,XCO,ET(5,30),OMEGA(3,30),TI(5,5),UR(2,3,3000) # ,FMO(3),AP(5),AQ(5),AO(3) XSTEP=(1.0-XCO)/NBIN DO 20 J=1,3 FN=2.*J*(2*J+1) IN=0 5 IN=IN+1 IC=3-IN UU=AR(5,1) UAVP=UU*(AR(3,1)*AR(4,1)+AR(6,1)) AOK=0.0 DOK=0.0 DO 60 JK=J,5 ETA(JK)=ET(JK,1) DETA(JK)=0.0 ETN=0.0 IF(JK.LT.5) ETN=ET(JK+1,1) AOK=AOK+TI(JK,J)*ETA(JK)/JK/2. 60 DOK=DOK+TI(JK,J)*(ETN-(1.0-1.0/JK)*ETA(JK)) XS=EXP(2.0*AR(1,1)) RHP=XS*(DOK*UU-AOK*UAVP) IF(IN.EQ.1)THEN F(1)=0. F(2)=1 FD(1)=0. FD(2)=2.*J ELSE F(1)=CIN FD(1)=2.*J*CIN UR(1,J,1)=CIN UR(2,J,1)=FD(1) ENDIF DO 10 N=2,NT DX=(AR(1,N)-AR(1,N-1))/2. UU=AR(5,N) UAV=UU*(AR(3,N)*AR(4,N)+AR(6,N)) X=EXP(AR(1,N)) FX=(X-XCO)/XSTEP N1=FX+0.5 DO 62 JK=J,5 IF(N1.EQ.0) THEN CET=ET(JK,2)-ET(JK,1) ETA(JK)=ET(JK,1)+CET*FX**2 DETA(JK)=2.*FX*X/XSTEP*CET ENDIF IF(N1.GT.0) THEN IF(N1.GT.NBIN-1) N1=NBIN-1 N2=N1+1 N3=N1+2 DXS=FX+0.5-N1 CET=ET(JK,N2)-ET(JK,N1) CETA=ET(JK,N1)-2.*ET(JK,N2)+ET(JK,N3) ETA(JK)=0.125*(3.*ET(JK,N1)+6.*ET(JK,N2)-ET(JK,N3))+DXS*CET+ # 0.5*DXS**2*CETA DETA(JK)=X/XSTEP*(CET+DXS*CETA) ENDIF 62 CONTINUE AOK=0.0 DOK=0.0 DO 64 JK=J,5 ETN=0.0 IF(JK.LT.5) ETN=ETA(JK+1) AOK=AOK+TI(JK,J)*ETA(JK)/JK/2. 64 DOK=DOK+TI(JK,J)*(ETN-(1.0-1.0/JK)*ETA(JK)+0.5*DETA(JK)) XS=X**2 RH=XS*(DOK*UU-AOK*UAV) AB=(UAV+UAVP)/2.0 AC=DX*(RH+RHP) AA=DX*(FN-AB) W=1.+DX-DX*AA DO 8 I=1,IC A1=F(I)+DX*FD(I) A2=AA*F(I)+(1.-DX)*FD(I) IF(I.EQ.1) A2=A2+AC F(I)=(A1*(1.+DX)+A2*DX)/W 8 FD(I)=(A1*AA+A2)/W IF(IC.EQ.1) THEN AA=AOK*XS UR(1,J,N)=F(1)+AA UR(2,J,N)=FD(1)+2.*AA+DOK*XS ENDIF UAVP=UAV RHP=RH 10 CONTINUE IF(IN.EQ.1) THEN CIN=-(FD(1)+(2*J+1)*F(1))/(FD(2)+(2*J+1)*F(2)) FMO(J)=F(1)+F(2)*CIN GO TO 5 ENDIF 20 CONTINUE RETURN END C******************************************************************** DOUBLE PRECISION FUNCTION FV(X) IMPLICIT REAL*8(A-H,O-Z) FV=X**2*(3.-2.*X) RETURN END