C C*********************************************************************** C BLOCK DATA BELEMD INCLUDE 'PREC' CHARACTER*2 EIDENT COMMON/IDMETL/ EIDENT(16) COMMON/ELMDAT/ AWI(16),ABNDI(16) DATA EIDENT/'H ','HE','SI','MG','AL','FE','C ','NA','S ','K ', * 'CA','NI','CR','N ','O ','NE'/ DATA AWI/1.008,4.0026,28.086,24.305,26.985,55.847,12.0111,22.99, * 32.06,39.102,40.08,58.71,52.0,14.007,16.0,20.179/ END C C*********************************************************************** C SUBROUTINE ATMOS C C READS A MODEL ATMOSPHERE FROM FILE ATMOS C THE FILE CONTAINS: C ATMOSPHERIC IDENTIFICATION LINE, FORMAT A C DEPTH SCALE TYPE, STARTING WITH T OR M, C T FOR TAU5000 SCALE, M FOR MASS SCALE FORMAT A C LOG GRAVITY FREE FORMAT C NUMBER OF DEPTH POINTS (NDEP), FREE FORMAT C (LOG DEPTH(K), T(K), NE(K), V(K), VTURB(K),K=1,NDEP), FREE FORMAT C ((NHI(I,K),I=1,5),NHII(K),K=1,NDEP) FREE FORMAT C C V AND VTURB IN KM PER SECOND, ELSE CGS UNITS C IF THERE ARE NO HYDROGEN POPULATIONS IN THE INPUT FILE C THEY ARE SET TO THE LTE VALUES C C A DEPTH-SCALE IS READ FROM FILE DSCALE. C THE DEPTH SCALE CAN BE GIVEN EITHER AS LOG COLUMN MASS OR C AS LOG TAU5000. THE UNIT IS GIVEN BY THE FIRST LETTER OF C THE DEPTH SCALE TYPE STRING C AND SHOULD OF COURSE BE THE SAME IN THE ATMOS AND THE DSCALE C FILES. THE DEPTH SCALE CAN EITHER BE GIVEN EXPLICITLY OR AS C THE LG DEPTH VALUE IN THE FIRST AND LAST DEPTH POINT. THIS IS C INDICATED BY THE SIGN OF NDEP IN THE DSCALE FILE. C THE FILE DSCALE CONTAINS : C C DEPTH SCALE IDENTIFICATION LINE FORMAT A C DEPTH SCALE TYPE, FIRST NON-BLANK CHARACTER T OR M, C T FOR TAU5000 SCALE, M FOR MASS SCALE FORMAT A C NDEP, DPCON1 FREE FORMAT C C THE REST OF THE FILE DEPENDS ON THE SIGN OF NDEP. C NDEP.LT.0: C C LG DEPTH(1), LG DEPTH(NDEP) FREE FORMAT C THE DEPTH POINTS ARE DISTRIBUTED EVENLY IN LG DEPTH C C NDEP.GT.0: C C (LG DEPTH(K), K=1,NDEP) FREE FORMAT C C DPTYPE IS SET EQUAL TO THE FIRST NON-BLANK LETTER OF C THE DEPTH SCALE TYPE STRING. C THE DEPTH SCALE UNIT IS DETERMINED FROM THE VALUE OF DPTYPE C C DPTYPE C M LG DEPTH IS LG COLUMN MASS, DPCON1= LG TAU5000(1) C T LG DEPTH IS LG TAU5000, DPCON1= LG COLUMN MASS(1) C H LG DEPTH IS HEIGHT (KM) DPCON1= LG TAU5000(1) C C IF DPCON1 IS 0 A START VALUE IS CALCULATED IN ROUTINE DPCONV C C A WORD OF WARNING: C THE PROGRAM OPERATES ON A TAU5000 SCALE. IF THE INPUT DEPTH SCALE C IS A COLUMN MASS SCALE IT IS CONVERTED TO A TAU5000 SCALE USING C TRAPEZOIDAL INTEGRATION IN ROUTINE DPCONV. THIS CAN LEAD TO C ERRORS FOR SPARSE DEPTH SCALES ESPECIALLY IN THE DEEPER PARTS C OF THE ATMOSPHERE. A CONVERSION TO A TAU5000 SCALE SHOULD BE C MADE USING A DENSER DEPTH SCALE OR THE TRAPEZOIDAL INTEGRATION C IN DPCONV SHOULD BE REPLACED. C C IF HYDROSTATIC EQUILIBRIUM (HSE) INTEGRATIONS ARE TO BE DONE IN C HSEINT, THE DEPTH SCALE HAS TO BE A COLUMN MASS SCALE C C INDEP IS THE MAXIMUM NUMBER OF POINTS IN THE INPUT ATMOSPHERE C IND IS SET TO 1 FOR STATIC ATMOSPHERES AND TO 2 FOR NON-STATIC ONES C C: C: ATMOS 86-12-15 MODIFICATIONS: (MATS CARLSSON) C: NH(I,K)=0 IN INPUT ATMOSPHERIC FILE IS ALLOWED C: THIS GAVE LOG OF 0.0 IN VERSION 1.0 15-AUG-86 C: C: 90-07-25 MODIFICATIONS: (MATS CARLSSON) C: TEST NDEP.EQ.MDEP1 TAKEN AWAY. TO BE USED WITH LOCAL C: OPERATOR VERSION C: C: 97-09-24 MODIFICATIONS: (MATS CARLSSON) C: DPTYPE='H' ADDED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CSLINE' INCLUDE 'CLU' C PARAMETER (INDEP=500) DIMENSION DPIN(INDEP),TIN(INDEP),ANEIN(INDEP),VIN(INDEP), * VTURIN(INDEP),ANHIN(6,INDEP) DIMENSION DP(MDEP) CHARACTER*72 TEXTA,TEXTD LOGICAL HLTE C C INPUT OF ATMOSPHERIC MODEL C CALL CSTRIP(LATMOS,'ATMOS') READ(LDUMS,100) ATMOID 100 FORMAT(A) CALL LJUST(ATMOID) READ(LDUMS,100) TEXTA CALL LJUST(TEXTA) DPTYPE=TEXTA(1:1) IF(DPTYPE.NE.'M' .AND. DPTYPE.NE.'T' .AND. DPTYPE.NE.'H') THEN CALL STOP('ATMOS: DEPTH SCALE TYPE UNKNOWN') ENDIF READ(LDUMS,*) GRAVLG GRAV=10.0**GRAVLG READ(LDUMS,*) KDEP IF(KDEP.GT.INDEP) CALL STOP('ATMOS: KDEP.GT.INDEP') READ(LDUMS,*) (DPIN(K),TIN(K),ANEIN(K),VIN(K),VTURIN(K),K=1,KDEP) READ(LDUMS,*,END=150) ((ANHIN(I,K),I=1,6),K=1,KDEP) HLTE=.FALSE. GOTO 200 150 CONTINUE C C NO HYDROGEN POPULATIONS, SET HLTE=.TRUE. C WHEN HLTE=.TRUE., THE HYDROGEN POPULATIONS NH WILL C NOT BE INTERPOLATED FROM ANHIN BUT SET EQUAL TO 0.0 C THIS WILL CAUSE THE OPACITY PACKAGE TO SET NH TO THE C LTE POPULATIONS C HLTE=.TRUE. C C INPUT OF DEPTH-SCALE C 200 CONTINUE CALL CSTRIP(LDSCAL,'DSCALE') READ(LDUMS,100) DPID CALL LJUST(DPID) READ(LDUMS,100) TEXTD CALL LJUST(TEXTD) IF(TEXTA(1:1).NE.TEXTD(1:1)) THEN CALL STOP('ATMOS: DIFFERENT DEPTH SCALES IN ATMOS AND DSCALE') ENDIF READ(LDUMS,*) NDEP,DPCON1 IF(ABS(NDEP).GT.MDEP) CALL STOP('ATMOS: NDEP.GT.MDEP') IF(NDEP.LT.0) THEN NDEP=-NDEP READ(LDUMS,*) DP(1),DP(NDEP) DLG=(DP(NDEP)-DP(1))/(NDEP-1) DO 300 K=2,NDEP-1 DP(K)=DP(K-1)+DLG 300 CONTINUE ELSE READ(LDUMS,*) (DP(K),K=1,NDEP) ENDIF C C STORE DP IN TAU OR CMASS OR HEIGHT C STORE DPCON1 IN CMASS(1) OR TAU(1) C IF DPTYPE.EQ.'H' CHANGE SIGN TO GET INCREASING SCALE C IF(DPTYPE.EQ.'M') THEN IF(DPCON1.NE.0.0) THEN TAU(1)=10.0**DPCON1 ELSE TAU(1)=0.0 ENDIF DO 360 K=1,NDEP CMASS(K)=10.0**DP(K) 360 CONTINUE ELSE IF(DPTYPE.EQ.'T') THEN IF(DPCON1.NE.0.0) THEN CMASS(1)=10.0**DPCON1 ELSE CMASS(1)=0.0 ENDIF DO 370 K=1,NDEP TAU(K)=10.0**DP(K) 370 CONTINUE ELSE IF(DPTYPE.EQ.'H') THEN IF(DPCON1.NE.0.0) THEN TAU(1)=10.0**DPCON1 ELSE TAU(1)=0.0 ENDIF CMASS(1)=0.0 DO 375 K=1,NDEP HEIGHT(K)=DP(K) DP(K)=-DP(K) 375 CONTINUE DO 377 K=1,KDEP DPIN(K)=-DPIN(K) 377 CONTINUE ENDIF C C INTERPOLATION OF ATMOSPHERIC MODEL, LINEAR INTERPOLATION C CONVERT V AND VTURB FROM KM PER SEC TO DOPPLER UNITS C INTERPOLATE IN THE LOG FOR NE AND NH C C STORE THE LOG OF ANEIN IN ANEIN C DO 380 K=1,KDEP ANEIN(K)=LOG(ANEIN(K)) 380 CONTINUE C C STORE THE LOG OF ANHIN IN ANHIN C C MODIFIED 86-12-15 BY MATS CARLSSON C IF ANHIN(I,1)=0.0 THEN ANHIN(I,K) IS SUPPOSED TO BE 0.0 C INSTEAD OF THE LOG A VALUE OF -999. IS STORED. THIS WILL C FLAG THE INTERPOLATION TO GIVE 0.0 C IF(.NOT.HLTE) THEN DO 400 I=1,6 IF(ANHIN(I,1).NE.0.0) THEN DO 390 K=1,KDEP ANHIN(I,K)=LOG(ANHIN(I,K)) 390 CONTINUE ELSE ANHIN(I,1)=-999. ENDIF 400 CONTINUE ENDIF C C SET IVEL=1 IF VEL IS ABOVE 1.E-6 ANYWHERE IN ATMOSPHERE C IF(DP(1).LT.DPIN(1) .OR. DP(NDEP).GT.DPIN(KDEP)) * CALL STOP('ATMOS: DEPTH SCALE IS OUTSIDE RANGE OF ATMOSPHERE') IVEL=0 L=1 DO 500 K=1,NDEP 410 L=L+1 IF(DP(K).GT.DPIN(L)) GOTO 410 APOL=(DP(K)-DPIN(L-1))/(DPIN(L)-DPIN(L-1)) TEMP(K)=TIN(L-1)+APOL*(TIN(L)-TIN(L-1)) NE(K)=EXP(ANEIN(L-1)+APOL*(ANEIN(L)-ANEIN(L-1))) VEL(K)=(VIN(L-1)+APOL*(VIN(L)-VIN(L-1)))/QNORM VTURB(K)=(VTURIN(L-1)+APOL*(VTURIN(L)-VTURIN(L-1)))/QNORM IF(ABS(VEL(K)).GT.1.E-6) IVEL=1 IF(.NOT.HLTE) THEN DO 420 I=1,6 IF(ANHIN(I,1).NE.-999.) THEN NH(I,K)=EXP(ANHIN(I,L-1)+APOL*(ANHIN(I,L)-ANHIN(I,L-1))) ELSE NH(I,K)=0.0 ENDIF 420 CONTINUE ELSE DO 430 I=1,6 NH(I,K)=0.0 430 CONTINUE ENDIF L=L-1 500 CONTINUE C DO 550 KR=1,MRAD IND(KR)=1+IVEL 550 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE ATOM(ICALL) C C READS AN ATOMIC MODEL AND CALCULATES EINSTEIN A,B COEFFICIENTS. C C A : EINSTEIN A COEFFICIENT C B : B C C NK : NUMBER OF LEVELS, CONTINUUM LEVELS INCLUDED C NRAD : NUMBER OF RADIATIVE TRANSITIONS C NLINE : NUMBER OF BOUND-BOUND TRANSITIONS C NCONT : NUMBER OF BOUND-FREE TRANSITIONS I DETAIL C IWIDE : =.TRUE. FOR FREQUENCY DEPENDENT TRANSITIONS C (CONTINUA AND WIDE LINES) C KTRANS(K) TRANSITION NUMBER K IS FREQUENCY DEPENDENT C TRANSITION NUMBER KTRANS C NRFIX : NUMBER OF TRANSITIONS WITH FIXED RATES C FOR FIXED RATES VARIABLES SEE ROUTINE FIXRAD C IRAD(K): LOWER LEVEL FOR RADIATIVE TRANSITION K C JRAD(K): UPPER '' C ALAMB : VACUUM WAVELENGTH IN ANGSTROM C IN PRINTOUT ROUTINES IT IS PRINTED AS EITHER VACUUM OR AIR C DEPENDING ON THE VALUE (.LT.2000 VACUUM, .GT.2000 AIR) C HNY4P : H*NY/4PI, NY IN UNITS OF A TYPICAL DOPPLER WIDTH C G : STATISTICAL WEIGHT C F : FOR LINES : OSCILLATOR STRENGTH C FOR CONTINUA: CROSS-SECTION AT LIMIT C THE WAVELENGTH DEPENDENCE OF THE CROSS-SECTION IS C ASSUMED TO BE A=A0*(NY0/NY)**3 IF NOT GIVEN EXPLICITLY C EV : ENERGY IN EV INPUT IN CM-1 C C NSTAR : BOLTZMANN POPULATION C BP : PLANCK FUNCTION C C IND IS SET TO 2 FOR WIDE LINES AND 1 FOR CONTINUA C C ICALL DETERMINES WHICH ATOM FILE TO READ C ICALL=1 FILE=ATOM C ICALL=2 FILE=ATOM2 A REDISTRIBUTION OF POPULATIONS IS MADE: C THE POPULATIONS OF THE LEVELS IN ATOM ARE REDISTRIBUTED C ON THE LEVELS IN ATOM2 THAT HAVE THE SAME LABEL. C THIS REDISTRIBUTION IS DONE IN ROUTINE REDIST AND IS C DONE ACCORDING TO THE STATISTICAL WEIGHTS OF THE LEVELS. C C: C: ATOM 88-01-19 MODIFICATIONS: (MATS CARLSSON) C: NQ IS TESTED AGAINST MQ TO AVOID ILLEGAL ADDRESS C: C: 88-05-03 MODIFICATIONS: (MATS CARLSSON, PHILIP JUDGE) C: EXTRA OPTION ADDED FOR INPUT OF FIXED TRANSITIONS WITH C: RATES CALCULATED USING RADIATION FIELD FROM FILE JFIX C: C: 90-01-26 MODIFICATIONS: (MATS CARLSSON) C: LABEL1 CHANGED TO CHARACTER*20 C: C: 90-07-25 MODIFICATIONS: (MATS CARLSSON) C: TEST NK.EQ.MK1 TAKEN AWAY. TO BE USED WITH LOCAL C: OPERATOR VERSION C: C: 93-01-18 MODIFICATIONS: (MATS CARLSSON) C: ADDED OPTION OF HAVING SEVERAL LINES BETWEEN TERMS C: FOR DOCUMENTATION OF VARIABLES, SEE ROUTINE ATOM C: C: 95-08-21 MODIFICATIONS: (MATS CARLSSON) C: CHECKS THAT ENERGY LEVELS ARE NOT TOO CLOSE IN ENERGY C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CCONST' INCLUDE 'CLU' INCLUDE 'CFIX' INCLUDE 'CTERM' C DIMENSION G1(MK) CHARACTER*20 LABEL1(MK) CHARACTER TEXT*80 C IF(ICALL.EQ.1) THEN C C CALL FROM START. CALCULATE CONSTANTS AND READ FILE ATOM C HCE=HH*CC/EE*1.E8 HC2=2.*HH*CC *1.E24 HCK=HH*CC/BK*1.E8 EK=EE/BK HNY4P=HH*CC/QNORM/4./PI*1.E-5 CALL CSTRIP(LATOM,'ATOM') ELSE IF(ICALL.EQ.2) THEN C C CALL FROM ATOM2. READ FILE ATOM2 C SAVE OLD VALUES OF NK, G AND LABEL TO ENABLE REDISTRIBUTION IN REDIST C NK1=NK DO 30 I=1,NK1 G1(I)=G(I) LABEL1(I)=LABEL(I) 30 CONTINUE CALL CSTRIP(LATOM2,'ATOM2') ELSE CALL STOP('ATOM: ICALL.GT.2') ENDIF C C READ INPUT VALUES AND LEFT JUSTIFY ATOMID C READ(LDUMS,40) ATOMID 40 FORMAT(A) CALL LJUST(ATOMID) C READ(LDUMS,*) ABND,AWGT AWGT=AWGT*UU READ(LDUMS,*) NK,NLINE,NCONT,NRFIX NRAD=NLINE+NCONT IF(ICALL.EQ.1) THEN IF(NK.GT.MK) CALL STOP('ATOM(1): NK.GT.MK') IF(NRAD.GT.MRAD) CALL STOP('ATOM(1): NRAD.GT.MRAD') IF(NLINE.GT.MLINE) CALL STOP('ATOM(1): NLINE.GT.MLINE') IF(NRFIX.GT.MRFIX) CALL STOP('ATOM(1): NRFIX.GT.MRFIX') ELSE IF(NK.GT.MK) CALL STOP('ATOM(2): NK.GT.MK') IF(NRAD.GT.MRAD) CALL STOP('ATOM(2): NRAD.GT.MRAD') IF(NLINE.GT.MLINE) CALL STOP('ATOM(2): NLINE.GT.MLINE') IF(NRFIX.GT.MRFIX) CALL STOP('ATOM(2): NRFIX.GT.MRFIX') ENDIF DO 120 I=1,NK READ(LDUMS,*) EV(I),G(I),LABEL(I),ION(I) EV(I)=EV(I)*CC *HH/EE IF(I.GT.1) THEN IF(EV(I).LE.EV(I-1)) THEN WRITE(LJOBLO,100) LABEL(I),LABEL(I-1) 100 FORMAT('ATOM: LEVELS ARE:',A,2X,A/ * ' ENERGY NOT MONOTONICALLY INCREASING'/ * ' (MAY BE DUE TO NUMERICAL ROUNDOFF', * ' FOR VERY CLOSE LEVELS)') CALL STOP(' ') ENDIF ENDIF DO 110 J=1,NK KRAD(I,J)=0. 110 CONTINUE 120 CONTINUE C IF(ICALL.EQ.2) CALL REDIST(NK1,G1,LABEL1) C C BOUND-BOUND TRANSITIONS IN DETAIL C CALCULATE LAMBDA, A AND B C IF QMAX OR Q0.LT.0 FREQUENCY POINTS IN DOPPLER UNITS ARE READ C C TERM OPTION: C A LINE MAY BE GIVEN AS A LINEAR COMBINATION OF PROFILES C TYPICAL USE IS THE FINE-STRUCTURE IN A TERM WHERE THE UPPER C AND LOWER LEVELS ARE ASSUMED TO BE POPULATED ACCORDING TO THE C STATISTICAL WEIGHTS. THE FREQUENCY QUADRATURE IS GIVEN AS BEFORE C WITH TERM OPTION FLAGGED BY THE IWIDE INPUT VARIABLE IO. C IO=0 NORMAL LINE AS BEFORE C IO=1 NORMAL LINE AS BEFORE WITH FREQUENCY DEPENDENT BACKGROUND C IO.GE.2 NUMBER OF COMPONENTS IN TERM C THEREAFTER FOLLOWS ONE EXTRA INPUT LINE PER COMPONENT: C DETERM WTERM FTERM GATERM GWTERM GQTERM C EXTRA VARIABLES: C MTERM MAXIMUM NUMBER OF COMPONENTS IN ANY ONE TRANSITION C MTRM TOTAL MAXIMUM NUMBER OF COMPONENTS C NTERM(KR) NUMBER OF COMPONENTS OF TERM, SET TO IO C KTERM(I,KR) TERM INDEX FOR COMPONENT I OF TRANSITION KR C DETERM(KTRM) ENERGY DIFFERENCE FOR COMPONENT I OF TRANSITION KR C IS ENERGY DIFFERENCE FOR TRANSITION KR + DETERM(KTRM) C WHERE KTRM IS KTERM(I,KR). UNIT CM-1. C WTERM(KTRM) WEIGHT OF COMPONENT KTRM IN WEIGHTED AVERAGE C FTERM(KTRM) OSCILLATOR STRENGTH C GATERM(KTRM) DAMPING CONSTANTS C GWTERM(KTRM) C GQTERM(KTRM) C ADTERM(K,KTRM) ADAMP FOR COMPONENT, CALCULATED IN DAMP C C MODIFIED ROUTINES ARE ATOM, DAMP AND PROFIL C KT=0 KTRM=0 DO 200 KR=1,NLINE READ(LDUMS,*) J,I,F(KR),NQ(KR),QMAX(KR),Q0(KR),IO, * GA(KR),GW(KR),GQ(KR) NTERM(KR)=IO IF(NQ(KR).GT.MQ) CALL STOP('ATOM: NQ.GT.MQ') IF(QMAX(KR).LT.0.0 .OR. Q0(KR).LT.0.0) THEN READ(LDUMS,*) (Q(NY,KR),NY=1,NQ(KR)) IF(ABS(Q(2,KR)).LT.ABS(Q(1,KR)) .AND. * ABS(Q(NQ(KR),KR)).GT.ABS(Q(NQ(KR)-1,KR))) IND(KR)=2 DO 150 NY=2,NQ(KR) IF(Q(NY,KR).LT.Q(NY-1,KR)) THEN CALL STOP('ATOM: Q NOT MONOTONICALLY INCREASING') ENDIF 150 CONTINUE ENDIF IF(NTERM(KR).GE.2) THEN IF(NTERM(KR).GT.MTERM) CALL STOP('ATOM: NTERM.GT.MTERM') DO 160 ITRM=1,NTERM(KR) KTRM=KTRM+1 IF(KTRM.GT.MTRM) CALL STOP('ATOM: KTRM.GT.MTRM') KTERM(ITRM,KR)=KTRM READ(LDUMS,*) DETERM(KTRM),WTERM(KTRM),FTERM(KTRM), * GATERM(KTRM),GWTERM(KTRM),GQTERM(KTRM) 160 CONTINUE IND(KR)=2 ENDIF IWIDE(KR)=IO.EQ.1 IF(IWIDE(KR)) THEN KT=KT+1 KTRANS(KR)=KT IND(KR)=2 ENDIF KRAD(I,J)=KR KRAD(J,I)=KR IRAD(KR)=I JRAD(KR)=J ALAMB(KR)=HCE/(EV(J)-EV(I)) A(KR)=F(KR)*6.671E15*G(I)/G(J)/ALAMB(KR)/ALAMB(KR) B(J,I)=ALAMB(KR)**3/HC2*A(KR) B(I,J)=G(J)/G(I)*B(J,I) 200 CONTINUE IF(KT+NRAD-NLINE.GT.MWIDE) CALL STOP('ATOM: NWIDE.GT.MWIDE') C C BOUND-FREE TRANSITIONS IN DETAIL C IF QMAX.LT.0.0 FREQUENCY POINTS IN ANGSTROM (STARTING WITH C THRESHOLD AND DECREASING) AND CROSSECTIONS IN CM2 ARE READ. C UNIT CONVERSION IN ROUTINE FREQC C DO 250 KR=NLINE+1,NRAD READ(LDUMS,*) J,I,F(KR),NQ(KR),QMAX(KR) IF(NQ(KR).GT.MQ) CALL STOP('ATOM: NQ.GT.MQ') IF(QMAX(KR).LT.0.0) THEN READ(LDUMS,*) (Q(NY,KR),ALFAC(NY,KR-NLINE),NY=1,NQ(KR)) DO 230 NY=2,NQ(KR) IF(Q(NY,KR).GT.Q(NY-1,KR)) THEN WRITE(LJOBLO,220) KR-NLINE 220 FORMAT(' ATOM: CONTINUUM TRANSITION NR',I3/ * ' WAVELENGTHS NOT DECREASING') CALL STOP(' ') ENDIF 230 CONTINUE ENDIF KTRANS(KR)=KR-NLINE+KT IRAD(KR)=I JRAD(KR)=J KRAD(I,J)=KR KRAD(J,I)=KR GA(KR)=0. GW(KR)=0. GQ(KR)=0. IWIDE(KR)=.TRUE. IND(KR)=1 ALAMB(KR)=HCE/(EV(J)-EV(I)) 250 CONTINUE C C CALCULATE PLANCK FUNCTION BP C DO 300 KR=1,NRAD U=CC/ALAMB(KR)*1.E8 DO 290 K=1,NDEP BP(K,KR)=PLANCK(U,TEMP(K)) 290 CONTINUE 300 CONTINUE C C INPUT FIXED TRANSITIONS PARAMETERS C IF(ICALL.EQ.1) THEN DO 400 KF=1,NRFIX READ(LDUMS,*) JFX(KF),IFX(KF),IPHO(KF),A0(KF),TRAD(KF), * ITRAD(KF) IF(ITRAD(KF).EQ.4) THEN IF(IPHO(KF).EQ.0) THEN READ(LDUMS,*) J,I,FF(KF),NQF(KF),QMAXF(KF),Q0F(KF),IO, * GAF(KF),GWF(KF),GQF(KF) IF(NQF(KF).GT.MQ) CALL STOP('ATOM: NQF.GT.MQ') IF(QMAXF(KF).LT.0.0 .OR. Q0F(KF).LT.0.0) THEN READ(LDUMS,*) (QF(NY,KF),NY=1,NQF(KF)) IF(QF(1,KF).LT.0.0) THEN CALL STOP('ATOM: QF MUST BE GIVEN ONE SIDED') ENDIF DO 350 NY=2,NQF(KF) IF(QF(NY,KF).LT.QF(NY-1,KF)) THEN CALL STOP('ATOM: QF NOT MONOTONICALLY INCREASING') ENDIF 350 CONTINUE ENDIF ELSE READ(LDUMS,*) J,I,FF(KF),NQF(KF),QMAXF(KF) IF(NQF(KF).GT.MQ) CALL STOP('ATOM: NQF.GT.MQ') IF(QMAXF(KF).LT.0.0) THEN READ(LDUMS,*) (QF(NY,KF),ALFACF(NY,KF),NY=1,NQF(KF)) DO 360 NY=2,NQF(KF) IF(QF(NY,KF).GT.QF(NY-1,KF)) THEN CALL STOP('ATOM: QF WAVELENGTHS NOT DECREASING') ENDIF 360 CONTINUE ENDIF ENDIF ENDIF 400 CONTINUE C C READ CHOICE OF COLLISIONAL SUBROUTINE C READ(LDUMS,40) TEXT CALL LJUST(TEXT) CROUT=TEXT(1:6) ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE WCHANG(E,NK1,NDEP1,ICALL) C C WRITES PLOT OF LOGARITHMIC CHANGE OF EACH LEVEL/TRANSITION TO OUTPUT C IF IWCHAN IS GREATER THAN 0 ONLY CHANGES OF EACH LEVEL IS PRINTED C IF IWCHAN IS LESS THAN 0 BOTH CHANGES OF EACH LEVEL AND CHANGES C IN THE RADIATION FIELD IS PRINTED C EM IS MAXIMUM CHANGE AT A CERTAIN DEPTH C C: WCHANG 90-07-31 MODIFICATIONS: (MATS CARLSSON) C: CHANGED TO CONTAIN DECLARATION OF E C: E PASSED IN ARGUMENT LIST INSTEAD OF IN COMMON TO C: AVOID NK=MK1 AND NDEP=MDEP1 REQUIREMENT C: C: 93-01-27 MODIFICATIONS: (JO BRULS, MATS CARLSSON) C: LOOP LIMIT MDEP CHANGED TO NDEP TO AVOID POSSIBLE C: OVERINDEXING C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION E(NK1,NDEP1) C DIMENSION EM(MDEP) INTEGER IM(MDEP) CHARACTER BLANK,VERT,MINUS,CHAR(0:15),LINE*112 DATA BLANK,VERT,MINUS,CHAR/' ','I','-','0','1','2','3','4','5', * '6','7','8','9','A','B','C','D','E','F'/ C IF(IWCHAN.EQ.0) RETURN C CALL CPUTIME('WCHANG',0,0,1) SMALL=1.E-10 IF(IWCHAN.LT.0) THEN KSTEP=-IWCHAN ELSE IF(ICALL.EQ.2) RETURN KSTEP=IWCHAN ENDIF C IF(ICALL.EQ.1) THEN LMAX=NK ELSE LMAX=NRAD ENDIF C C IF NUMBER OF LEVELS/TRANSITIONS IS LESS THAN 16, USE 16 AS BASE, C ELSE USE 10 AND STORE MAX ERROR. C IF(LMAX.LE.15) THEN LBASE=16 ELSE LBASE=10 C C STORE MAXIMUM ERROR ON EACH DEPTH IN EM C STORE CORRESPONDING INDEX IN IM C DO 60 K=1,NDEP EM(K)=0.0 IF(ICALL.EQ.1) THEN DO 30 I=1,NK EM(K)=MAX(EM(K),E(I,K)) 30 CONTINUE ELSE DO 40 I=1,NRAD EM(K)=MAX(EM(K),ABS(EJ(I,K))) 40 CONTINUE ENDIF IM(K)=10.*LOG10(MAX(EM(K),SMALL))+102.5 IM(K)=MIN(IM(K),112) 60 CONTINUE ENDIF C DO 800 L=0,LMAX,LBASE IMAX=MIN(L+LBASE-1,LMAX) IMIN=MAX(1,L) IF(ICALL.EQ.1) THEN WRITE(LOUT,100) IMIN,IMAX,IT 100 FORMAT('1 PLOT OF LOGARITHMIC CHANGE OF LEVEL', * I3,'-',I3,' FOR ITERATION',I3) ELSE WRITE(LOUT,110) IMIN,IMAX,IT-1 110 FORMAT('1 PLOT OF LOGARITHMIC CHANGE OF TRANSITION', * I3,'-',I3,' FOR ITERATION',I3) ENDIF WRITE(LOUT,200) (I,I=-10,1) 200 FORMAT(3X,'K',1X,12(I3,7X)) DO 700 K=1,NDEP,KSTEP LINE=BLANK DO 400 I=2,112,10 LINE(I:I)=VERT 400 CONTINUE DO 500 I=IMIN,IMAX IF(ICALL.EQ.1) THEN CHANG=MAX(ABS(E(I,K)),SMALL) SIGN=E(I,K) ELSE CHANG=MAX(ABS(EJ(I,K)),SMALL) SIGN=EJ(I,K) ENDIF IC=10.*LOG10(CHANG)+102.5 IF(IC.GT.112) IC=112 LINE(IC :IC)=CHAR(MOD(I,LBASE)) IF(SIGN.LT.0.) LINE(IC-1:IC-1)=MINUS 500 CONTINUE IF(LMAX.GT.15) THEN IF(LINE(IM(K):IM(K)).EQ.' ' .OR. * LINE(IM(K):IM(K)).EQ.VERT) LINE(IM(K):IM(K))='*' ENDIF WRITE(LOUT,600) K,LINE 600 FORMAT(1X,I3,2X,A) 700 CONTINUE 800 CONTINUE CALL CPUTIME('WCHANG',0,1,1) C RETURN END C C*********************************************************************** C SUBROUTINE WEVEC(E,NK1,NDEP1) C C WRITES ERROR VECTOR TO OUTPUT LISTING C C: WEVEC 90-07-31 MODIFICATIONS: (MATS CARLSSON) C: CHANGED TO CONTAIN DECLARATION OF E C: E PASSED IN ARGUMENT LIST INSTEAD OF IN COMMON TO C: AVOID NK=MK1 AND NDEP=MDEP1 REQUIREMENT C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION E(NK1,NDEP1) C IF(IWEVEC.EQ.0) RETURN C C OUTPUT OF ERROR VECTOR C IF(IWTEST.EQ.0) WRITE(LOUT,50) IT 50 FORMAT('1 ITERATION',I3) WRITE(LOUT,70) 70 FORMAT(//' ERROR VECTOR') DO 400 M=1,NK,12 IMAX=MIN(NK,M+11) WRITE(LOUT,100) (I,I=M,IMAX) 100 FORMAT(//' ERROR VECTOR'/' K',12(I6,4X)/) DO 300 K=1,NDEP,IWEVEC WRITE(LOUT,200) K,(E(I,K),I=M,IMAX) 200 FORMAT(1X,I3,12(1P,E10.2)) 300 CONTINUE 400 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE ATOM2 C C INITIALIZES FOR FORMAL SOLUTION IF THERE IS AN ATOMIC MODEL ATOM2 C CHECKS IF THE FILE ATOM2 EXISTS. IF NOT, ATM2 IS SET .FALSE. AND C NOTHING IS DONE. C IF THE FILE EXISTS, ATM2 IS SET .TRUE. AND THE ATOMIC MODEL IS READ. C THE CONVERGED POPULATIONS ARE REDISTRIBUTED TO THE NEW LEVELS IN C ROUTINE REDIST CALLED FROM ATOM. NEW FREQUENCY QUADRATURE, C OPACITIES, DAMPING PARAMETERS AND PROFILES ARE CALCULATED TO BE C USED IN THE FORMAL SOLUTION. THE FORMAL SOLUTION IS ITERATED TO C GET THE SCATTERING CONTRIBUTION TO THE SOURCE FUNCTION CONSISTENT C WITH THE RADIATION FIELD C C: C: ATOM2 88-06-22 MODIFICATIONS: (MATS CARLSSON) C: ICONV IS SET TO 0 DURING SCATTERING LAMBDA ITERATIONS TO C: AVOID PRINTOUT IN WTEST CALLED FROM TRPT C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: SCATTERING VERSION OF TRANSP USED INSTEAD OF C: SCATTERING LAMBDA ITERATIONS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C CHARACTER TEXT DATA LMAX/20/ C C CHECK IF THERE IS AN ATOM2 FILE C CALL CPUTIME('ATOM2 ',0,0,2) CALL OPEN(LATOM2,'ATOM2',1,'UNKNOWN') READ(LATOM2,100,END=900) TEXT 100 FORMAT(A) CALL CLOSE(LATOM2) C C THERE IS AN ATOM2 FILE C CALL ATOM(2) CALL WHEAD CALL WATOM(2) CALL FREQ CALL CLOSE(LJNY) CALL CLOSE(LOPC) CALL OPAC(0) CALL WLTE CALL DAMP CALL PROFIL CALL LTEEQW C C GET CONSISTENT JBAR IN SCATTERING INTEGRAL C ISCAT0=ISCAT ICONV0=ICONV ISCAT=1 ICONV=0 CALL TRPT ISCAT=ISCAT0 ICONV=ICONV0 CALL CPUTIME('ATOM2 ',0,1,2) RETURN C C THERE IS NO ATOM2 FILE C 900 CONTINUE CALL STOP('ATOM2: ATOM2 NOT EXISTING') RETURN END C C*********************************************************************** C SUBROUTINE CA2COL C C COLLISIONAL RATES FOR CALCIUM II C THE ROUTINE WRITTEN FOR BOUND LEVELS PLUS ONE CONTINUUM LEVEL C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CCONST' INCLUDE 'CLU' C DIMENSION CE(MK-1),CI(MK-1) C C CALCULATE C BOUND-BOUND C DO 300 J=2,NK-1 IMAX=J-1 READ(LDUMS,*) (CE(I),I=1,IMAX) DO 200 I=1,J-1 DO 100 K=1,NDEP C(J,I,K)=CE(I)*SQRT(5000./TEMP(K))*NE(K) C(I,J,K)=C(J,I,K)*NSTAR(J,K)/NSTAR(I,K) 100 CONTINUE 200 CONTINUE 300 CONTINUE C C COLLISIONS WITH HYDROGEN C IF(NK-1.EQ.5) THEN I2=2 J3=3 I4=4 J5=5 DO 350 K=1,NDEP DC=1.5E-9*(TEMP(K)/1000.)**0.33*NH(1,K) C(I4,J5,K)=C(I4,J5,K)+DC C(I2,J3,K)=C(I2,J3,K)+DC C(J5,I4,K)=C(J5,I4,K)+DC*NSTAR(I4,K)/NSTAR(J5,K) C(J3,I2,K)=C(J3,I2,K)+DC*NSTAR(I2,K)/NSTAR(J3,K) 350 CONTINUE ENDIF C C CALCULATE C BOUND-FREE C READ(LDUMS,*) (CI(I),I=1,NK-1) DO 500 K=1,NDEP DO 400 I=1,NK-1 C(I,NK,K)=CI(I)*NE(K)*SQRT(TEMP(K))*EXP(-EK/TEMP(K)* * (EV(NK)-EV(I))) C(NK,I,K)=C(I,NK,K)*NSTAR(I,K)/NSTAR(NK,K) 400 CONTINUE 500 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE CNTRB(NY,KR) C C CALCULATES INTENSITY CONTRIBUTION FUNCTIONS AND TAUQ C FOR MU=NMU C: C: CNTRB 88-01-21 NEW ROUTINE: (MATS CARLSSON) C: CALCULATES INTENSITY CONTRIBUTION FUNCTIONS AND TAUQ C: FOR MU=NMU C: C: 88-14-05 MODIFICATIONS: (PHILIP JUDGE) C: CALCULATIONS FLUX CONTRIBUTION FUNCTIONS C: C: 92-11-27 MODIFICATIONS: (MATS CARLSSON) C: SAVE IREC ADDED C: C: 93-03-25 MODIFICATIONS: (GIUSEPPE SEVERINO, MATS CARLSSON) C: CALCULATES CNTRBI AND CNTRBR FOR MU=XMU(NMU) C: C: 95-10-20 MODIFICATIONS: (MATS CARLSSON) C: ERROR CAUSING CNTRB TO BE EXECUTED ONLY FOR IWLINE.NE.0 AND C: NOT WHEN IDLCNT.NE.0 CORRECTED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLGMX' INCLUDE 'CTAUQQ' INCLUDE 'CCNTRB' INCLUDE 'CLU' C LOGICAL CONT SAVE IREC C IF(IWLINE.EQ.0 .AND. IDLCNT.EQ.0) RETURN C IF(KR.EQ.1 .AND. NY.EQ.1) THEN CALL CPUTIME('CNTRB ',0,0,1) IREC=0 CALL REWIND(LPHI) ENDIF IREC=IREC+1 TAU(0)=-TAU(1) TAUQ(0)=0.0 X(0)=0.0 I=IRAD(KR) J=JRAD(KR) C C CALCULATE SOME CONSTANTS C IF(KR.GT.NLINE) THEN CONT=.TRUE. ELSE CONT=.FALSE. GIJK=G(I)/G(J) HN3C2=A(KR)/B(J,I) DO 280 K=1,NDEP Z(K)=N(I,K)-GIJK*N(J,K) SL(K,KR)=HN3C2*N(J,K)*GIJK/Z(K) GIJ(K)=GIJK 280 CONTINUE ENDIF IF(IWIDE(KR)) KT=KTRANS(KR) C CALL READJ(IREC) DO 285 K=1,NDEP SBCK(K)=SC(K)+SCAT(K)*JNY(K) 285 CONTINUE IF(CONT) THEN HN3C2=2.*HH*FRQ(NY,KT)/CC *FRQ(NY,KT)/CC *FRQ(NY,KT) DO 290 K=1,NDEP GIJ(K)=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FRQ(NY,KT)/BK/TEMP(K)) ALFA(K)=ALFAC(NY,KR-NLINE) Z(K)=N(I,K)-GIJ(K)*N(J,K) SL(K,KR)=HN3C2*N(J,K)*GIJ(K)/Z(K) 290 CONTINUE ENDIF C DO 295 MU=1,NMU IF(.NOT.CONT .AND. (IND(KR).EQ.2 .OR. MU.EQ.1)) THEN CALL READP ENDIF 295 CONTINUE MU=NMU IF(.NOT.CONT) THEN DO 300 K=1,NDEP ALFA(K)=B(I,J)*PHI(K)*HNY4P 300 CONTINUE ENDIF DO 310 K=1,NDEP X(K)=Z(K)*ALFA(K)/XNORM(K)+XCONT(K) RNY(K)=XCONT(K)/X(K) S(K)=(1.-RNY(K))*SL(K,KR)+RNY(K)*SBCK(K) C C CALCULATE TAUQ C TAUQ(K)=TAUQ(K-1)+0.5*(X(K-1)+X(K))*(TAU(K)-TAU(K-1)) TAUQQ(NY,K)=TAUQ(K) C C CALCULATE CNTRBI AND PATHL C CNTRBI=LOG(10)*TAU0*KAPPA/KAPPA0*S*EXP(-TAUQ) C IF(S(K).GT.0.0) THEN CNTRLG=LOG(TAU(K))+LOG(X(K))+LOG(S(K))-TAUQ(K)/XMU(MU) ELSE CNTRLG=-99. ENDIF IF(CNTRLG.GT.-69.) THEN CNTRBI(NY,K)=LOG(10.)/XMU(MU)*EXP(CNTRLG) C* C* 14-05-88 PGJ MODIFICATIONS START C* FIRST IS CONTRIBUTION PER UNIT LOG TAU(5000) C* (EQUIVALENT TO MAGAIN FUNCTION CI) C* SECOND IS CONTRIBUTION PER UNIT HEIGHT C* THIRD IS CONTRIBUTION PER UNIT LOG (TAU0) (LINE CENTRE OPTICAL DEPTH) C* FOURTH IS CONTRIBUTION PER UNIT LOG10 (COLUMN MASS) C*C(1). CNTRBF(NY,K) = 2.*PI*LOG(10.) * *TAU(K)*X(K)*S(K)*EXPINT(2,TAUQ(K),EX2) C* C*C(2). C* IF(K .GT. 1)THEN C DTAUDH = -(TAUQ(K) - TAUQ(K-1))/( HEIGHT(K)-HEIGHT(K-1)) C* CNTRBF(NY,K) = 2.*PI*DTAUDH*S(K)*EXPINT(2,TAUQ(K),EX2) C* ELSE C* CNTRBF(NY,K)=0. C* ENDIF C* C*C(3). C* IF(K .GT. 1)THEN C* DTDT0 = (TAUQ(K) - TAUQ(K-1)) C* * / LOG10 ( TAUQQ(1,K)/TAUQQ(1,K-1) ) C* CNTRBF(NY,K) = 2.*PI*LOG(10.)*DTDT0 C* * *S(K)*EXPINT(2,TAUQ(K),EX2) C* END IF C* C*C(4) C* IF(K .GT. 1)THEN C* DTAUDM=(TAUQ(K)-TAUQ(K-1)) / (CMASS(K) - CMASS(K-1) ) C* CNTRBF(NY,K) = 2.*PI*LOG(10.)*DTAUDM C* * *CMASS(K)*S(K)*EXPINT(2,TAUQ(K),EX2) C* ELSE C* CNTRBF(NY,K)=0. C* ENDIF C* C* 14-05-88 MODIFICATIONS END C* ELSE CNTRBI(NY,K)=0.0 C* 14-05-88 MODIFICATIONS START CNTRBF(NY,K) = 0.0 C* 14-05-88 MODIFICATIONS END ENDIF PATHL(NY,K)=1./X(K)/XNORM(K) C C CALCULATE CNTRBR C CNTRBR=LOG(10)*TAU0*KAPPAL/KAPPA0*(1-SL/IC)*EXP(-TAUR) C DTAUR=(KAPPAL+KAPPAC *SC/IC)DZ C IF(.NOT.CONT) THEN XL=Z(K)*ALFA(K)/XNORM(K) IF(K.EQ.1) THEN TAUR=TAUQ(1) XRK=XL+XCONT(K)*SBCK(K)/ICONT(KR,K) ELSE XRM=XRK XRK=XL+XCONT(K)*SBCK(K)/ICONT(KR,K) TAUR=TAUR+(XRK+XRM)/(X(K-1)+X(K))*(TAUQ(K)-TAUQ(K-1)) ENDIF IF(TAUR.GE.0.0) THEN CNTRBR(NY,K)=LOG(10.)*TAU(K)*XL*(1.-SL(K,KR)/ICONT(KR,K))* * EXP(-TAUR/XMU(MU))/XMU(MU) ELSE CNTRBR(NY,K)=0.0 ENDIF ENDIF 310 CONTINUE IF(KR.EQ.NRAD .AND. NY.EQ.NQ(KR)) THEN CALL CPUTIME('CNTRB ',0,1,1) ENDIF C RETURN END C C*********************************************************************** C FUNCTION CONVL(ALAMB) C C CONVERTS VACUUM LAMBDA TO AIR (.GT.2000 A) OR VACUUM C (.LT.2000 A) LAMBDA C CONSTANTS FROM STARLINK PROGRAM IUEDR C INCLUDE 'PREC' IF(ALAMB.LT.2000.) THEN CONVL=ALAMB ELSE CONVL=ALAMB/(1.0+2.735182E-4+131.4182/ALAMB/ALAMB+ * 2.76249E8/ALAMB/ALAMB/ALAMB/ALAMB) ENDIF C RETURN END C C********************************************************************** C SUBROUTINE CSTRIP(IUNIT,FILE) C C STRIPS COMMENT LINES (LINES WITH * IN POSITION 1) FROM INPUT FILE C UNIT IUNIT. THE FILE IS OPENED AND THE NON-COMMENT LINES ARE C WRITTEN TO UNIT LDUMS C INCLUDE 'PREC' INCLUDE 'CLU' C CHARACTER*132 TEXT CHARACTER*(*) FILE C CALL OPEN(IUNIT,FILE,1,'OLD') CALL REWIND(LDUMS) 100 CONTINUE READ(IUNIT,200,END=900) TEXT IF(TEXT(1:1).NE.'*') WRITE(LDUMS,200) TEXT 200 FORMAT(A) GOTO 100 C 900 CONTINUE CALL REWIND(LDUMS) CALL CLOSE(IUNIT) RETURN END C C********************************************************************** C SUBROUTINE CPUTIME(TEXT,ISTART,IPRI,LEVEL) C C GIVES THE CPU TIME FROM ONE CALL TO ANOTHER C IT CAN BE ACTIVATED ON UP TO 5 DIFFERENT LEVELS. C THE MEANINGS OF ISTART AND IPRI: C ISTART = 0 ZERO CLOCK ON EXIT C IPRI .GT. 1 WRITE TIME, IPRI-1 BLANK LINES BEFORE AND C AFTER PRINTOUT C INCLUDE 'PREC' INCLUDE 'CLU' C CHARACTER TEXT*6,FMT*29 DIMENSION B(5) SAVE B DATA B/0.,0.,0.,0.,0./ C A=SECOND() CPT=A-B(LEVEL) IF(ISTART.EQ.0) B(LEVEL)=A IF(IPRI.GT.0) THEN DO 50 I=1,IPRI-1 WRITE(LTIME,40) 40 FORMAT(1X) 50 CONTINUE WRITE(FMT,100) 6-LEVEL,LEVEL+2 100 FORMAT('(',I2,'X,A6,',I2,'X,''CPU='',F8.3,'' S'')') WRITE(LTIME,FMT) TEXT,CPT DO 150 I=1,IPRI-1 WRITE(LTIME,40) 150 CONTINUE ENDIF RETURN END C C************************************************************************ C FUNCTION DLAMB(QQ,KR) C C CONVERTS FROM DELTA LAMBDA IN Q UNITS TO DELTA LAMBDA IN ANGSTROM C QQ IS DELTA LAMBDA IN TYPICAL DOPPLER WIDTH UNITS C KR IS RADIATIVE TRANSITION INDEX C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CCONST' C QN=QNORM*1.E5/CC DLAMB=-ALAMB(KR)*QQ*QN/(1.0+QQ*QN) C RETURN END C C************************************************************************ C SUBROUTINE DPCONV C C CONVERTS BETWEEN TAU(5000) SCALE AND MASS-SCALE OR VICE VERSA C CALCULATES TAU(5000) FROM KAPPA AND MASS-SCALE (DPTYPE=M) OR C MASS-SCALE FROM KAPPA AND TAU(5000) (DPTYPE=T) C CALCULATES HEIGHT SCALE C TRAPEZ FORMULA IS USED C C: DPCONV 97-09-24 MODIFICATIONS: (MATS CARLSSON) C: DPTYPE='H' ADDED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' C C CALCULATE TAU SCALE C IF(DPTYPE.EQ.'M') THEN IF(TAU(1).EQ.0.) THEN TAU(1)=XNORM(1)/RHO(1)*CMASS(1) ENDIF DO 100 K=2,NDEP TAU(K)=TAU(K-1)+(XNORM(K)+XNORM(K-1))/(RHO(K)+RHO(K-1))* * (CMASS(K)-CMASS(K-1)) 100 CONTINUE ELSE IF(DPTYPE.EQ.'H') THEN IF(TAU(1).EQ.0.0) THEN DO 130 K=2,3 TAU(K)=TAU(K-1)+0.5*(XNORM(K)+XNORM(K-1))* * (HEIGHT(K-1)-HEIGHT(K))*1.E5 130 CONTINUE TAU(1)=EXP(2*LOG(TAU(2))-LOG(TAU(3))) ENDIF DO 150 K=2,NDEP TAU(K)=TAU(K-1)+0.5*(XNORM(K)+XNORM(K-1))* * (HEIGHT(K-1)-HEIGHT(K))*1.E5 150 CONTINUE ENDIF C C CALCULATE COLUMN MASS SCALE C IF(DPTYPE.NE.'M') THEN IF(CMASS(1).EQ.0.) THEN CMASS(1)=RHO(1)/XNORM(1)*TAU(1) ENDIF DO 200 K=2,NDEP CMASS(K)=CMASS(K-1)+(RHO(K)+RHO(K-1))/(XNORM(K)+XNORM(K-1))* * (TAU(K)-TAU(K-1)) 200 CONTINUE ENDIF C C CALCULATE HEIGHT SCALE. SET H=0 CLOSEST TO TAU=1, INCREASING OUTWARDS. C IF(DPTYPE.NE.'H') THEN HEIGHT(1)=0.0 DO 300 K=2,NDEP HEIGHT(K)=HEIGHT(K-1)-2.0*(CMASS(K)-CMASS(K-1))/ * (RHO(K)+RHO(K-1)) 300 CONTINUE DO 400 K=NDEP,1,-1 IF(TAU(K).LT.1.0) GOTO 500 400 CONTINUE K=1 500 K0=K IF(K0.NE.1 .AND. K0.NE.NDEP) THEN H0=HEIGHT(K0)+(HEIGHT(K0+1)-HEIGHT(K0))/(TAU(K0+1)-TAU(K0))* * (1.0-TAU(K0)) ELSE H0=HEIGHT(K0) ENDIF DO 600 K=1,NDEP HEIGHT(K)=(HEIGHT(K)-H0)*1.E-5 600 CONTINUE ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE DAMP C C CALCULATES DAMPING PARAMETERS C INPUT: GA, GW, GS IN ROUTINE ATOM WHERE C GAMMA = GR + GV + GS AND C GR=GA C GS=GQ*NE C GV=GW*GVW(H) C C GVW(H) = 8.08 * (1+0.41*NHE/NH) * VH**0.6 * NH * C6**0.4 C VH**0.6 = (8*BK*T/PI*(1/MH+1/AWGT))**0.3 C C6 = 1.01E-32 * (13.6*ION)**2 * C (1/(EV(CONT)-EV(UPPER))**2 - 1/(EV(CONT)-EV(LOWER))**2) C REFERENCE MIHALAS 1978: STELLAR ATMOSPHERES, 2ND ED C NHE/NH IS SET TO 0.1 C NB: C6 SHOULD BE DIVIDED BY 2*PI WHEN COMPARED WITH GREY AND UNSOLD C C C FOR H-ALPHA GS=4.737E-7*NH(1,K), REF. LINEAR-B C C: C: DAMP 88-01-19 MODIFICATIONS: (MATS CARLSSON) C: IF GW=0 LEVELS WITHOUT OVERLYING CONTINUUM ARE ALLOWED C: C: 88-05-20 MODIFICATIONS: (MATS CARLSSON) C: IF GQ.LT.0 IT IS ASSUMED TO BE LG(C4) C: C: 93-01-18 MODIFICATIONS: (MATS CARLSSON) C: ADDED OPTION OF HAVING SEVERAL LINES BETWEEN TERMS C: FOR DOCUMENTATION OF VARIABLES, SEE ROUTINE ATOM C: C 00-11-27 MODIFICATIONS: (PAUL BARKLEM) C ABO THEORY X-SECTIONS IMPLEMENTED C IF GW > 20. ASSUMED X-SECTION DATA C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CTERM' C ZERO=0. IF(IWDAMP.GT.0) CALL WDAMP(0,0,ZERO,ZERO,ZERO,ZERO,ZERO) C DO 400 K=1,NDEP T=TEMP(K) PE=NE(K)*BK*T TOTHI=0.0 DO 100 J=1,5 TOTHI=TOTHI+NH(J,K) 100 CONTINUE C C CALCULATE DOPPLER WIDTH C DNYD(K)=SQRT(2.*BK*T/AWGT)*1.E-5/QNORM C C INCLUDE MICROTURBULENCE C DNYD(K)=SQRT(DNYD(K)**2+VTURB(K)**2) C C CALCULATE GAMMA, ADAMP C C ADD TM C IF(K==NDEP)OPEN(58,FILE='GV_BARKLEM.dat') C END ADD TM DO 330 KR=1,NLINE I=IRAD(KR) J=JRAD(KR) C C THIS SECTION BELOW ALTERED PB, TO USE X-SECTIONS NOV 2000 C IF (GW(KR).GE.20.0) THEN GV=GW(KR) SIGMA=INT(GV)*2.80028E-21 XALPHA=GV-INT(GV) GX=2.-XALPHA*.5 GX=GX-1.0 GAMMAF=1+(-.5748646+(.9512363+(-.6998588+ * (.4245549-.1010678*GX)*GX)*GX)*GX)*GX GV=(4./PI)**(XALPHA*0.5)*GAMMAF*1.E4*SIGMA VBAR=SQRT(21173.*T*(1./1.008+1./(AWGT/UU))) GV=GV*((VBAR/1.E4)**(1.-XALPHA)) GV=GV*TOTHI*1.E6*2. C ADD TM C IF(K==NDEP) WRITE(58,'(F10.2,F8.3,ES11.2)') ALAMB(KR), GW(KR), GV C IF(K==NDEP) WRITE(58,'(A,I5,A,ES10.2)')'GV_BARKLEM(',KR,') = ',GV C IF(K==NDEP) WRITE(*,'(A,I6,A,ES10.2)')'GV_BARKLEM(',KR,') = ',GV C C625 = GV/(8.411*(8.*BK*T/PI*(1./(1.008*UU)+1./AWGT))**0.3* C * TOTHI) C IF(K==NDEP) THEN C WRITE(58,'(A,F7.2)')'Log(C6_BARKLM) = ',2.5*LOG10(C625/2.) C WRITE(*,'(A,F7.2)')'Log(C6_BARKLM) = ',2.5*LOG10(C625/2.) C ENDIF C END ADD TM C C FIND CONTINUUM LEVEL C ELSEIF(GW(KR).NE.0.0) THEN DO 200 IC=J+1,NK IF(ION(IC).EQ.ION(J)+1) GOTO 300 200 CONTINUE CALL STOP(' DAMP: NO OVERLYING CONTINUUM') 300 CONTINUE ZZ=ION(I) C625=1.283984E-12*ZZ**0.8*(1./(EV(IC)-EV(J))**2- * 1./(EV(IC)-EV(I))**2)**.4 GV=GW(KR)*8.411*(8.*BK*T/PI*(1./(1.008*UU)+1./AWGT))**0.3* * TOTHI*C625 C ADD TM C IF(K==NDEP) WRITE(58,'(A,I5,A,ES10.2)')'GV_UNSOLD(',KR,') = ',GV C IF(K==NDEP)WRITE(58,'(A,F7.2)')'Log(C6_UNSOLD) = ',2.5*LOG10(C625) C IF(K==NDEP) WRITE(*,'(A,I6,A,ES10.2)')'GV_UNSOLD(',KR,') = ',GV C IF(K==NDEP) WRITE(*,'(A,F7.2)')'Log(C6_UNSOLD) = ',2.5*LOG10(C625) C END ADD TM ELSE GV=0.0 ENDIF GR=GA(KR) IF(GQ(KR).GE.0.0) THEN GS=GQ(KR)*NE(K) ELSE C FORMULA FROM GRAY, P.237 GSLG=19.4+2./3.*GQ(KR)+LOG10(NE(K)*BK)+LOG10(TEMP(K))/6. GS=10.**GSLG ENDIF IF(ATOMID.EQ.'H ' .AND. J.EQ.3 .AND. I.EQ.2) * GS=4.737E-7*NH(1,K) GAMMA=GR+GV+GS DOP=DNYD(K)*QNORM/ALAMB(KR)*1.E13 ADAMP(K,KR)=GAMMA/(4.*PI*DOP) DDP=DNYD(K)*QNORM*1.E8/CC *ALAMB(KR) IF(IWDAMP.GT.0) CALL WDAMP(K,KR,DDP,GAMMA,GR,GV,GS) IF(NTERM(KR).GE.2) THEN DO 320 ITRM=1,NTERM(KR) KTRM=KTERM(ITRM,KR) C C THIS SECTION BELOW ALTERED PB, TO USE X-SECTIONS NOV 2000 C IF(GWTERM(KTRM).GE.20.0) THEN GV=GWTERM(KTRM) SIGMA=INT(GV)*2.80028E-21 XALPHA=GV-INT(GV) GX=2.-XALPHA*.5 GX=GX-1.0 GAMMAF=1+(-.5748646+(.9512363+(-.6998588+ * (.4245549-.1010678*GX)*GX)*GX)*GX)*GX GV=(4./PI)**(XALPHA*0.5)*GAMMAF*1.E4*SIGMA VBAR=SQRT(21173.*T*(1./1.008+1./(AWGT/UU))) GV=GV*((VBAR/1.E4)**(1.-XALPHA)) GV=GV*TOTHI*1.E6*2. C ADD TM C IF(K==NDEP) WRITE(58,'(A,I5,A,ES10.2)')'GV_BARKLEM(',KR,') = ',GV C IF(K==NDEP) WRITE(*,'(A,I6,A,ES10.2)')'GV_BARKLEM(',KR,') = ',GV C C625 = GV/(8.411*(8.*BK*T/PI*(1./(1.008*UU)+1./AWGT))**0.3* C * TOTHI) C IF(K==NDEP) THEN C WRITE(58,'(A,F7.2)')'Log(C6_BARKLM) = ',2.5*LOG10(C625/2.) C WRITE(*,'(A,F7.2)')'Log(C6_BARKLM) = ',2.5*LOG10(C625/2.) C ENDIF C END ADD TM ELSEIF(GWTERM(KTRM).NE.0.0) THEN GV=GWTERM(KTRM)*8.411*(8.*BK*T/PI* * (1./(1.008*UU)+1./AWGT))**0.3*TOTHI*C625 C ADD TM C IF(K==NDEP) WRITE(58,'(A,I5,A,ES10.2)')'GV_UNSOLD(',KR,') = ',GV C IF(K==NDEP)WRITE(58,'(A,F7.2)')'Log(C6_UNSOLD) = ',2.5*LOG10(C625) C IF(K==NDEP) WRITE(*,'(A,I6,A,ES10.2)')'GV_UNSOLD(',KR,') = ',GV C IF(K==NDEP) WRITE(*,'(A,F7.2)')'Log(C6_UNSOLD) = ',2.5*LOG10(C625) C END ADD TM ELSE GV=0.0 ENDIF GR=GATERM(KTRM) IF(GQTERM(KTRM).GE.0.0) THEN GS=GQTERM(KTRM)*NE(K) ELSE C FORMULA FROM GRAY, P.237 GSLG=19.4+2./3.*GQTERM(KTRM)+LOG10(NE(K)*BK)+ * LOG10(TEMP(K))/6. GS=10.**GSLG ENDIF IF(ATOMID.EQ.'H ' .AND. J.EQ.3 .AND. I.EQ.2) * GS=4.737E-7*NH(1,K) GAMMA=GR+GV+GS DOP=DNYD(K)*QNORM/ALAMB(KR)*1.E13 ADTERM(K,KTRM)=GAMMA/(4.*PI*DOP) 320 CONTINUE ENDIF 330 CONTINUE 400 CONTINUE C ADD TM CLOSE(58) WRITE(*,*)'DAMP ROUTINE EXECUTED.' C END ADD TM RETURN END C C*********************************************************************** C SUBROUTINE DSCAL2 C C INTERPOLATES A NEW DEPTH SCALE. C THE IDEA IS TO GET THE SAME MAX D LG TAUNY THROUGHOUT THE MODEL. C C CODED BY M.CARLSSON (JULY 1983) C: C: DSCAL2 87-12-15 MODIFICATIONS: (PHILIP JUDGE, MATS CARLSSON) C: IONIZATION RATIO CHANGES TAKEN INTO ACCOUNT FOR H RUNS C: C: 88-06-01 MODIFICATIONS: (PHILIP JUDGE, MATS CARLSSON) C: COLUMN MASS TAKEN INTO ACCOUNT C: C: 90-12-12 MODIFICATIONS: (PHILIP JUDGE, MATS CARLSSON) C: OUTPUT FORMAT CHANGED TO F15.9 FROM F12.6 C: C: 97-09-24 MODIFICATIONS: (MATS CARLSSON) C: DPTYPE='H' ADDED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CLGMX' INCLUDE 'CLU' INCLUDE 'CATOM' C DIMENSION TMX(MDEP),DP(MDEP) C C CALCULATE A ''MAX D LG TAU'' TAU SCALE C IF DPTYPE='M' INCLUDE ALSO A CHECK ON COLUMN MASS SCALE C C* ALTERATION BY M. CARLSSON/ P. JUDGE JAN 88: C* IF HYDROGEN, INCLUDE A CHECK ON THE IONIZATION FRACTION C* STEP IS NUMBER OF DEPTHPOINTS OVER IONIZATION ZONE CORRESPONDING C* TO MAX D LG TAU = 0.3 C C* ALTERATION BY P. JUDGE 22-2-1988 C* ATTEMPT TO MAKE CHANGE IN LOGM SMALLER THAN A SENSIBLE VALUE IN ORDER C* FOR THE COLUMN MASS INTEGRATIONS (E.G. FOR HEIGHT, HSEINT ETC) SMOOTH C* ENOUGH FOR REASONABLE ACCURACY C* DMASSM IS THE MEAN CHANGE IN COLUMN MASS THROUGHOUT ATMOSPHERE C* DMASS IS THE CHANGE IN COLUMN MASS AT EACH DEPTH POINT C* TMX(1) = 0.0 C* C* HYDROGEN- SPECIAL CASE C* IF (ATOMID.EQ.'H') THEN STEP = FLOAT(NDEP)/5.0 DMASSM = LOG10(CMASS(NDEP)/CMASS(1))/FLOAT(NDEP) DO 50 K = 2,NDEP DION = ABS(N(1,K-1)/TOTN(K-1)-N(1,K)/TOTN(K)) DMASS = LOG10(CMASS(K)/CMASS(K-1)) C* C* IF CHANGE IN COLUMN MASS IS GREATER THAN THE MEAN THEN PUT SOME C* POINTS IN C* IF (0.3*STEP*DION.GT.DMASS) THEN TMX(K) = TMX(K-1) + MAX(STEP*0.3*DION,DLGTMX(K)) ELSE IF (DMASS.GT.DMASSM) THEN TMX(K) = TMX(K-1) + MAX(DMASS,DLGTMX(K)) ELSE TMX(K) = TMX(K-1) + DLGTMX(K) END IF C* 50 CONTINUE ELSE DO 100 K=2,NDEP IF(DPTYPE.EQ.'M') THEN DLGCM=LOG10(CMASS(K))-LOG10(CMASS(K-1)) ELSE DLGCM=0.0 ENDIF TMX(K)=TMX(K-1)+MAX(DLGTMX(K),DLGCM) 100 CONTINUE ENDIF DTMX=TMX(NDEP)/(NDEP-1) C C PUT CURRENT DEPTH SCALE IN DP C IF(DPTYPE.EQ.'M') THEN DPCON1=LOG10(TAU(1)) DO 110 K=1,NDEP DP(K)=LOG10(CMASS(K)) 110 CONTINUE ELSE IF(DPTYPE.EQ.'T') THEN DPCON1=LOG10(CMASS(1)) DO 120 K=1,NDEP DP(K)=LOG10(TAU(K)) 120 CONTINUE ELSE IF(DPTYPE.EQ.'H') THEN DPCON1=LOG10(TAU(1)) DO 130 K=1,NDEP DP(K)=HEIGHT(K) 130 CONTINUE ENDIF C C INTERPOLATE A NEW DEPTH SCALE AND WRITE TO DSCAL2 C CALL OPEN(LDSCA2,'DSCAL2',1,'NEW') WRITE(LDSCA2,150) DPID(1:61) 150 FORMAT('* DEPTH SCALE FROM DSCAL2'/' DSCAL2 ON ',A) IF(DPTYPE.EQ.'M') THEN WRITE(LDSCA2,160) 160 FORMAT(' MASS SCALE') ELSE WRITE(LDSCA2,170) 170 FORMAT(' TAU5000 SCALE') ENDIF WRITE(LDSCA2,210) NDEP,DPCON1 210 FORMAT(I4,F15.6) WRITE(LDSCA2,220) DP(1)+0.5E-6 220 FORMAT(F15.9) T=0. L=1 DO 300 K=2,NDEP-1 T=T+DTMX 250 L=L+1 IF(T.GT.TMX(L)) GOTO 250 CM2=DP(L-1)+(T-TMX(L-1))/(TMX(L)-TMX(L-1))*(DP(L)-DP(L-1)) WRITE(LDSCA2,220) CM2 L=L-1 300 CONTINUE WRITE(LDSCA2,220) DP(NDEP)-0.5E-6 CALL CLOSE(LDSCA2) C RETURN END C C*********************************************************************** C SUBROUTINE EQSYST(N,NDIM,A,B,NEWMAT) C C SOLVES THE EQUATION SYSTEM A*X=B. C WHEN NEWMAT=TRUE, THE SYSTEM IS REARRANGED INTO U*X=L*B, WHERE U C IS UPPER AND L IS LOWER TRIANGULAR. THESE ARE THEN REUSED IN LATER C CALLS WITH NEWMAT=FALSE AND NEW RIGHT HAND SIDES B. THE SOLUTION C VECTOR IS RETURNED IN B. NO PIVOTING, I.E. THE MATRIX A IS ASSUMED C TO HAVE NONZERO DIAGONAL ELEMENTS. C C CODED BY: A. NORDLUND (FEB-1979) C C THIS IS A MODIFIED VERSION OF EQSYST WHICH TESTS FOR ZERO ELEMENTS C BELOW THE DIAGONAL AND ALSO STOPS AT THE LAST NON-ZERO ELEMENT ABOVE C THE DIAGONAL. CONSIDERABLE SAVINGS ARE OBTAINED FOR LOOSE MATRICES. C C THIS IS A COLUMN ORIENTED VERSION (M. CARLSSON JAN-1986) C TEMPORARY SCALARS BL, ALL, ALM AND BK ARE USED TO SHOW THE COMPILER C THAT THERE IS NO VECTOR DEPENDENCY IN THE INNERMOST DO-LOOP C INCLUDE 'PREC' PARAMETER (MDIM=10000) DIMENSION A(NDIM,NDIM),B(NDIM) INTEGER LASTN(MDIM) LOGICAL NEWMAT SAVE LASTN C IF(NEWMAT) THEN C C FIND THE LAST NON-ZERO ELEMENT IN EACH COLUMN C IF(N.GT.MDIM) CALL STOP('EQSYST: N.GT.MDIM') DO 120 L=1,N DO 100 K=N,L+1,-1 IF(A(K,L).NE.0.0) GOTO 110 100 CONTINUE K=L 110 LASTN(L)=K 120 CONTINUE C C COLUMN LOOP: ELIMINATE ELEMENTS BELOW THE DIAGONAL IN COLUMN L. C DO 500 L=1,N-1 C C STORE -A(K,L)/A(L,L) IN ELEMENT A(K,L) C MULTIPLY RIGHT HAND SIDE WITH -A(K,L)/A(L,L) C ALL=A(L,L) BL=B(L) DO 200 K=L+1,LASTN(L) A(K,L)=-A(K,L)/ALL B(K)=B(K)+A(K,L)*BL 200 CONTINUE C C ADD FRACTION -A(K,L)/A(L,L) OF ROW L TO ROW K. C IN EACH COLUMN GO THROUGH ALL ROWS C DO 400 M=L+1,N IF(A(L,M).NE.0.0) THEN ALM=A(L,M) LASTN(M)=MAX(LASTN(L),LASTN(M)) DO 300 K=L+1,LASTN(L) A(K,M)=A(K,M)+A(K,L)*ALM 300 CONTINUE END IF 400 CONTINUE 500 CONTINUE ELSE DO 700 L=1,N-1 BL=B(L) DO 600 K=L+1,LASTN(L) B(K)=B(K)+A(K,L)*BL 600 CONTINUE 700 CONTINUE END IF C C BACKSUBSTITUTE C DO 900 K=N,1,-1 BK=B(K) DO 800 L=K+1,N BK=BK-A(K,L)*B(L) 800 CONTINUE B(K)=BK/A(K,K) 900 CONTINUE C RETURN END C C********************************************************************** C SUBROUTINE EQWDTH(WW,KR) C C COMPUTES EQUIVALENT WIDTHS C ALGORITH FROM A.NORDLUND, IN METHODS OF RADIATIVE TRANSFER, 1984, C ED. W. KALKOFEN. C CODED BY M.SAXNER 1984 C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' C DIMENSION XL(MQ),RFLX(MQ),DFLX(MQ) C DO 100 NY=1,NQ(KR) XL(NY)=DLAMB(Q(NY,KR),KR) RFLX(NY)=FLUX(NY,KR)/FLUX(0,KR) 100 CONTINUE C IF(IND(KR).EQ.1) THEN DFLX(1)=0. ELSE DFLX(1)=(RFLX(2)-RFLX(1))/(XL(2)-XL(1)) ENDIF DFLX(NQ(KR))=(RFLX(NQ(KR))-RFLX(NQ(KR)-1))/ * (XL(NQ(KR))-XL(NQ(KR)-1)) DO 200 NY=2,NQ(KR)-1 DFLX(NY)=(RFLX(NY+1)-RFLX(NY-1))/(XL(NY+1)-XL(NY-1)) 200 CONTINUE C WSUM=0. DO 300 NY=2,NQ(KR) H=XL(NY-1)-XL(NY) WSUM=WSUM+0.5*H*(RFLX(NY-1)+RFLX(NY)-0.166667*H* * (DFLX(NY-1)-DFLX(NY))) 300 CONTINUE WW=XL(1)-XL(NQ(KR))-WSUM IF(IND(KR).EQ.1) WW=2.0*WW C RETURN END C C*********************************************************************** C SUBROUTINE ESCAPE C C APPROXIMATE SOLVER OF TRANSFER EQUATION FOR GIVEN POPULATIONS C IN 2ND ORDER ESCAPE PROBABILITY APPROXN. FOR LINES C DETAILED SOLUTION FOR CONTINUA C C SIMILAR TO ROUTINE TRPT C C APPROXIMATION IS AS FOLLOWS: C C FROM RYBICKI (IN 'METHODS IN RADIATIVE TRANSFER', ED KALKOFEN (1984), C PAGE 21), WE REPLACE THE NET RADIATIVE BRACKET, ORDNARILY COMPUTED C FROM THE TRANSFER AND STATISTICAL EQULIBRIUM EQUATIONS, WITH THE C ESCAPE PROBABILITY APPROXIMATIONS. C C THE NRB IS DEFINED AS C C NRB = 1 - (JBAR/SL) C C WHERE JBAR = MEAN INTENSITY C SL = LINE SOURCE FUNCTION C C 1. FIRST ORDER APPROXIMATION (LOCAL ESCAPE PROBABILITY) C C NRB(TAU) = PESC(TAU) C C PESC IS THE PROBABILITY OF ESCAPE FROM POINT WITH C LINE CENTRE OPTICAL DEPTH TAU C (CF. RYBICKI EQN 4.3) C C 2. SECOND ORDER APPROXIMATION C C NRB(TAU) = PESC(TAU) + INTEGRAL [ DT*PESC(T)* DS(T)/DT ] C C WHERE THE INTEGRAL EXTENDS FROM TAU TO INFINITY C C (CF. RYBICKI EQN 6.4) C C C: C: ESCAPE 87-12-11 NEW ROUTINE: (PHILIP JUDGE) C: LOCAL ESCAPE PROBABILITY TRANSFER SOLVER C: C: 88-07-15 MODIFICATIONS: (PHILIP JUDGE) C: 2ND ORDER FORM C: C: 89-03-06 MODIFICATIONS: (PHILIP JUDGE) C: FIRST ORDER COMMENT LINE ADDED: C: TO OBTAIN FIRST ORDER APPROXIMATION SIMPLY ADD THE LINE C: MARKED C*1 WHICH IS COMMENTED OUT. C: C: 89-03-07 MODIFICATIONS: (PHILIP JUDGE) C: FIRST ORDER VERSION ONLY- 2ND ORDER DIVERGES NEAR TMIN C: IN MANY STARS AND NEAR THE TOP OF THE ATMOSPHERE. THIS IS C: PRESUMABLY OWING TO THE FAST VARIATION OF S(TAU) IN THESE C: REGIONS WHICH IS CONTRARY TO THE ASSUMPTIONS OF THE METHOD C: (SEE E.G. CANFIELD ET AL., AP J 248, 82 (1981) ). C: C: 96-03-19 MODIFICATIONS: (SALVADOR ALBA VILLEGAS) C: RIJ WAS NOT SET TO ZERO IN SPITE OF COMMENT LINE SAYING SO. C: LEAD TO INCORRECT RESULTS WHEN IWEQW.NE.0 SINCE RIJ WAS C: THEN SET IN LTEEQW C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLGMX' INCLUDE 'CLU' C DIMENSION SESC(MDEP) C LOGICAL CONT DATA ZERO /0.0/ C CALL CPUTIME('ESCAP2',0,0,1) C C BOUND-BOUND TRANSITIONS C CALL REWIND(LOPC) IREC=0 DO 200 KR=1,NLINE DO 80 NY=1,NQ(KR) IREC=IREC+1 CALL READX 80 CONTINUE I=IRAD(KR) J=JRAD(KR) GIJK=G(I)/G(J) DO 110 K=1,NDEP Z(K)=N(I,K)-GIJK*N(J,K) IF(Z(K).LE.0.0) THEN IF(ICONV.EQ.1) THEN WRITE(LJOBLO,100) KR,K 100 FORMAT(' WARNING IN ESCAPE: NEGATIVE OPACITIES', * ' IN LINE',I3,' DEPTH',I4) ENDIF ENDIF C C SESC ARRAY SAVED FOR THE 2ND ORDER APPROXIMATION C C SESC = 1 / ( N(I,K)/N(J,K)*G(J)/G(I) - 1 ) C => SL = 2*HH*NU**3/CC**2 * SESC C SESC(K) = N(J,K)*GIJK/Z(K) C PROFILE FUNCTION PHI: CALL VOIGT(ADAMP(K,KR),ZERO,H) PJPHI = H / (DNYD(K)*1.772453851) ALFA(K)=B(I,J)*PJPHI*HNY4P X(K)=Z(K)*ALFA(K)/XNORM(K) 110 CONTINUE C C COMPUTE OPTICAL DEPTHS C TAUQ(1)=TAU(1)*X(1) DO 120 K=2,NDEP TAUQ(K)=TAUQ(K-1)+0.5*(X(K)+X(K-1))*(TAU(K)-TAU(K-1)) 120 CONTINUE C C OUTPUT OF THE NET RADIATIVE RATES IN TERMS OF THE ESCAPE C PROBABILITIES. C REPLACE THE OLD RADIATIVE RATES (EINSTEIN-A VALUES) C WITH ESCAPE PROBABILITY VALUES. C REPLACE THE LOWER-UPPER RATES WITH ZERO C DO 130 K=1,NDEP C C FIRST ORDER PART: C FIRST = PESC(TAUQ(K),ADAMP(K,KR),KR-NLINE) C C SECOND ORDER PART: USE THE LOCAL VARIABLE SESC PROPORTIONAL TO THE C LINE SOURCE FUNCTION VALUES. C SECORD = 0.0 DO 125 KSUM = K + 1,NDEP SECORD = SECORD + PESC((TAUQ(KSUM)+TAUQ(KSUM-1))/2.0, * (ADAMP(KSUM,KR)+ADAMP(KSUM-1,KR))/2.0,KR-NLINE)* * (SESC(KSUM)-SESC(KSUM-1)) 125 CONTINUE SECORD=SECORD/SESC(K) C* C* (INCLUDE THE FOLLOWING LINE FOR FIRST ORDER APPROXIMATION) SECORD = 0.0 C* C* 6-3-89 PGJ ALTERATION: INTEGRAL OVER SESC CHANGED SIGN C* (ERROR IN ORIGINAL VERSION) C* RJI(K,KR) = A(KR)* (FIRST-SECORD) RIJ(K,KR) = 0.0 130 CONTINUE 200 CONTINUE C C BOUND-FREE TRANSITIONS: TREATED IN DETAIL C DO 400 KR=NLINE+1,NRAD I=IRAD(KR) J=JRAD(KR) C C CALCULATE SOME CONSTANTS C DO 270 K=1,NDEP RIJ(K,KR)=0. RJI(K,KR)=0. 270 CONTINUE CONT=.TRUE. DO 275 K=1,NDEP WPHI(K,KR)=1.0 275 CONTINUE KT=KTRANS(KR) C DO 380 NY=1,NQ(KR) IREC=IREC+1 CALL READX CALL READJ(IREC) DO 285 K=1,NDEP SBCK(K)=SC(K)+SCAT(K)*JNY(K) 285 CONTINUE HN3C2=2.*HH*FRQ(NY,KT)/CC *FRQ(NY,KT)/CC *FRQ(NY,KT) DO 290 K=1,NDEP GIJ(K)=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FRQ(NY,KT)/BK/TEMP(K)) ALFA(K)=ALFAC(NY,KR-NLINE) Z(K)=N(I,K)-GIJ(K)*N(J,K) SL(K,KR)=HN3C2*N(J,K)*GIJ(K)/Z(K) 290 CONTINUE C DO 370 MU=1,NMU XMY=XMU(MU) WQMU=WQ(NY,KR)*WMU(MU) IF(IND(KR).EQ.2 .OR. MU.EQ.1) THEN DO 310 K=1,NDEP X(K)=Z(K)*ALFA(K)/XNORM(K)+XCONT(K) RNY(K)=XCONT(K)/X(K) S(K)=(1.-RNY(K))*SL(K,KR)+RNY(K)*SBCK(K) 310 CONTINUE ENDIF C CALL TRANSP C WQMUH=WQMU/HNY4P/FRQ(NY,KT)*FRQ(0,KT) DO 360 K=1,NDEP RIJ(K,KR)=RIJ(K,KR)+WQMUH*WPHI(K,KR)*ALFA(K)*P(K) RJI(K,KR)=RJI(K,KR)+WQMUH*WPHI(K,KR)*ALFA(K)* * GIJ(K)*(P(K)+HN3C2) 360 CONTINUE 370 CONTINUE 380 CONTINUE 400 CONTINUE C CALL CPUTIME('ESCAP2',0,1,1) RETURN END C C*********************************************************************** C FUNCTION EXINT1(X) C C CALCULATES THE EXPONENTIAL INTEGRAL FUNCTION E1(X) C C REFERENCE: HANDBOOK OF MATHEMATICAL FUNCTIONS, C NBS APPLIED MATHEMATICS SERIES 55 (1964) (ED. ABRAMOVITZ C AND STEGUN), EQUATIONS 5.1.53 AND 5.1.54. C INCLUDE 'PREC' DIMENSION A(6),B(4) DATA (A(I),I=1,6) /-0.57721566,0.99999193,-0.24991055,0.05519968, * -0.00976004,0.00107857/ DATA(B(I),I=1,4) /2.334733,0.250621,3.330657,1.681534/ C IF(X.GT.1.0) GO TO 10 S=0.0 Y=1.0 DO 2 I=1,6 S=S+A(I)*Y Y=Y*X 2 CONTINUE EXINT1=S-LOG(X) RETURN C 10 CONTINUE IF(X.GT.80.0) GO TO 21 S=X*X+B(1)*X+B(2) S=S/(X*X+B(3)*X+B(4)) EXINT1=S*EXP(-X)/X RETURN C 21 EXINT1=0.0 RETURN END C C********************************************************************* C FUNCTION EXPINT(N,X,EX) C C COMPUTES THE N'TH EXPONENTIAL INTEGRAL OF X C INPUT - X, INDEPENDENT VARIABLE (-100. .LE. X .LE. +100.) C N, ORDER OF DESIRED EXPONENTIAL INTEGRAL (1 .LE. N .LE. 8) C OUTPUT - EXPINT, THE DESIRED RESULT C EX, EXPF(-X) C NOTE RETURNS WITH E1(0)=0, (NOT INFINITY). C BASED ON THE SHARE ROUTINE NUEXPI, WRITTEN BY J. W. COOLEY, C COURANT INSTITUTE OF MATHEMATICAL SCIENCES, NEW YORK UNIVERSITY C OBTAINED FROM RUDOLF LOESER C-----GENERAL COMPILATION OF 1 AUGUST 1967. C INCLUDE 'PREC' DIMENSION TAB(20),XINT(7) DATA XINT/1.,2.,3.,4.,5.,6.,7./ DATA TAB /.2707662555,.2131473101,.1746297218,.1477309984, 1.1280843565,.1131470205,.1014028126,.0919145454,.0840790292, 1.0774922515,.0718735405,.0670215610,.0627878642,.0590604044, 1.0557529077,.0527977953,.0501413386,.0477402600,.0455592945, 1.0435694088/ DATA XSAVE /0./ C U=X IF(U)603,602,603 602 EX=1. IF(N-1)800,800,801 800 EXPINT=0. GOTO 777 801 EXPINT=1./XINT(N-1) GOTO 777 603 IF(U-XSAVE)604,503,604 604 XSAVE=U XM=-U EMX=EXP(XM) C C SELECT METHOD FOR COMPUTING EI(XM) C IF(XM-24.5)501,400,400 501 IF(XM-5.)502,300,300 502 IF(XM+1.)100,200,200 503 EISAVE=-ARG EXSAVE=EMX C C NOW RECURSE TO HIGHER ORDERS C IF(N-1)507,507,505 505 DO 506 I=2,N EISAVE=(U*EISAVE-EXSAVE)/(-XINT(I-1)) 506 CONTINUE 507 EXPINT=EISAVE EX=EXSAVE 777 RETURN C C EI(XM) FOR XM .LT. -1.0 C HASTINGS POLYNOMIAL APPROXIMATION C 100 ARG=((((((U+8.573328740 )*U+18.05901697 )*U+8.634760893 )*U *+.2677737343)/XM)*EMX)/((((U+9.573322345 )*U+25.63295615 )*U *+21.09965308 )*U+3.958496923 ) GOTO 503 C EI(XM) FOR -1. .LE. XM .LT. 5.0 C POWER SERIES EXPANSION ABOUT ZERO 200 ARG=LOG(ABS(XM)) ARG=((((((((((((((((.41159050E-14*XM+.71745406E-13)*XM+.76404637E- *12)*XM+.11395905E-10)*XM+.17540077E-9)*XM+.23002666E-8)*XM+.275360 *18E-7)*XM+.30588626E-6)*XM+.31003842E-5)*XM+.28346991E-4)*XM+.2314 *8057E-3)*XM+.0016666574)*XM+.010416668)*XM+.055555572)*XM+.25)*XM+ *.99999999)*XM+.57721566)+ARG GOTO 503 C C EI(XM) FOR 5.0 .LE. XM .LT. 24.5 C TABLE LOOK-UP AND INTERPOLATION C 300 I=XM+.5 XZERO=I DELTA=XZERO-XM ARG=TAB(I-4) IF(DELTA)303,305,303 303 Y=ARG DELTAX=DELTA/XZERO POWER=1./DELTAX DO 304 I=1,7 POWER=POWER*DELTAX Y=((Y-POWER/XZERO)*DELTA)/XINT(I) ARG=ARG+Y IF(ABS(Y/ARG)-1.E-8)305,304,304 304 CONTINUE 305 ARG=EMX*ARG GOTO 503 C EI(XM) FOR 24.5 .LE. XM C TRUNCATED CONTINUED FRACTION 400 ARG=((((XM-15.)*XM+58.)*XM-50.)*EMX)/((((XM-16.)*XM+72.)*XM-96.) **XM+24.) GOTO 503 END C C*************************************************************** C SUBROUTINE FIXRAD C C CALCULATES FIXED RADIATIVE RATES C FOLLOWS AUER, HEASLEY, MILKEY, KPNO 555,18 C C NRFIX : NUMBER OF FIXED TRANSITIONS C JFX : UPPER LEVEL C IFX : LOWER LEVEL C IPHO =1 CONTINUUM C =0 LINE C A0 : CROSS-SECTION AT LIMIT C TRAD : BRIGHTNESS TEMPERATURE FOR CONTINUUM C ITRAD RADIATION TEMPERATURE OPTION C =1 TR=T C =2 PHOTOSPHERIC OPTION C =3 CHROMOSPHERIC OPTION C C : COLLISIONAL TRANSITION RATES + OTHER FIXED RATES C C: C: FIXRAD 88-05-04 MODIFICATIONS: (MATS CARLSSON) C: ITRAD=4 OPTION ADDED. RATE IS THEN CALCULATED USING C: NORMAL ATOMIC PARAMETERS READ BY ATOM INTO COMMON BLOCK C: CFIX AND RADIATION FIELD READ FROM FILE JFIX C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CCONST' C DO 500 KF=1,NRFIX IF(ITRAD(KF).EQ.4) THEN IF(IPHO(KF).EQ.0) THEN CALL FREQLF(KF) CALL DAMPF(KF) ELSE CALL FREQCF(KF) ENDIF IPR=IWATOM CALL JFIX(KF,IPR) GOTO 500 ENDIF J=JFX(KF) I=IFX(KF) XXI=EV(J)-EV(I) CLAM=HCE/XXI*1.E-8 DO 400 K=1,NDEP T=TEMP(K) DELT=-1.0 IF(K.NE.NDEP) DELT=TEMP(K)-TEMP(K+1) TR=T IF(ITRAD(KF).EQ.2) THEN IF(DELT.GE.0.) TR=TRAD(KF) IF(DELT.LT.0. .AND. TRAD(KF).GT.T) TR=TRAD(KF) ELSE IF(ITRAD(KF).EQ.3 .AND. TRAD(KF).LT.T * .AND. DELT.GE.0.) THEN TR=TRAD(KF) ENDIF IF(IPHO(KF).EQ.1) THEN CA0=8.*PI*CC/CLAM/CLAM/CLAM*A0(KF) RIJK=CA0*RINT(EK/TR*XXI) RJIK=RIJK IF(TR.NE.T) RJIK=CA0*RINT(EK/T*XXI) HKTT=EK*XXI/(T*TR) STIMD=RINT1(T,TR,HKTT)*CA0 RJIK=(RJIK+STIMD)*NSTAR(I,K)/NSTAR(J,K) ELSE EX=EXP(-EK/TR*XXI) X2=(EE*XXI/HH)**2 RIJK=7.421E-22*A0(KF)*X2*EX/(1.-EX) RJIK=RIJK IF(TR.NE.T) THEN EX=EXP(-EK/T*XXI) RJIK=7.421E-22*A0(KF)*X2*EX/(1.-EX) ENDIF RJIK=RJIK*G(I)/(EX*G(J)) ENDIF C(I,J,K)=C(I,J,K)+RIJK C(J,I,K)=C(J,I,K)+RJIK 400 CONTINUE 500 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE FREQC(KR) C C GIVES QUADRATURE POINTS AND WEIGHTS FOR CONTINUA C THE POINTS ARE EQUIDISTANT IN FREQUENCY, THE WEIGHTS ARE C TRAPEZOIDAL C CALCULATES CROSSECTIONS FROM NY**-3 DEPENDENCE IF NOT GIVEN C EXPLICITLY C C QMAX .LT. 0 : Q AND ALFAC GIVEN (FROM ATOM INPUT). ONLY C WEIGHTS CALCULATED AND Q TRANSFORMED FROM ANGSTROM TO DOPPLER C UNITS. C C NQ TOTAL NUMBER OF QUADRATURE POINTS (IN) C QMAX MINIMUM VALUE OF WAVELENGTH (ANGSTROM) (IN) C FRQ QUADRATURE POINTS (S-1) (OUT) C Q QUADRATURE POINTS (DOPPLER UNITS) (OUT) C WQ QUADRATURE WEIGHTS (OUT) C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CSLINE' INCLUDE 'CCONST' C KT=KTRANS(KR) FRQ(0,KT)=CC/ALAMB(KR)*1.E8 IF(QMAX(KR).GE.0.0) THEN FRQ(NQ(KR),KT)=CC/QMAX(KR)*1.E8 DFRQ=(FRQ(NQ(KR),KT)-FRQ(0,KT))/(NQ(KR)-1) DO 100 NY=1,NQ(KR) FRQ(NY,KT)=FRQ(0,KT)+(NY-1)*DFRQ Q(NY,KR)=CC*1.E-5/QNORM*(FRQ(NY,KT)/FRQ(0,KT)-1.0) ALFAC(NY,KR-NLINE)=F(KR)*(FRQ(0,KT)/FRQ(NY,KT))**3 100 CONTINUE ELSE DO 200 NY=1,NQ(KR) FRQ(NY,KT)=CC/Q(NY,KR)*1.0E8 Q(NY,KR)=CC*1.0E-5/QNORM*(FRQ(NY,KT)/FRQ(0,KT)-1.0) 200 CONTINUE ENDIF C WQ(1,KR)=0.5*(Q(2,KR)-Q(1,KR)) DO 300 NY=2,NQ(KR)-1 WQ(NY,KR)=0.5*(Q(NY+1,KR)-Q(NY-1,KR)) 300 CONTINUE WQ(NQ(KR),KR)=0.5*(Q(NQ(KR),KR)-Q(NQ(KR)-1,KR)) C RETURN END C C****************************************************************** C SUBROUTINE FREQL(KR) C C GIVES FREQUENCY QUADRATURE POINTS AND CORRESPONDING WEIGHTS. C ONE- OR TWO-SIDED CASES. THE POINTS AND WEIGHTS ARE BASED ON C THE TRAPEZOIDAL RULE APPLIED TO A MAPPING FUNCTION THAT MAPS C EQUIDISTANT POINTS IN X TO EQUIDISTANT POINTS IN Q IN THE C CORE (Q .LT. Q0), AND EQUIDISTANT POINTS IN LOG(Q) IN THE WINGS C (Q .GT. Q0). C C NQ : TOTAL NUMBER OF QUADRATURE POINTS (IN) C QMAX : MAXIMUM VALUE OF FREQUENCY (DOPPLER UNITS) (IN) C IND : =1 FOR ONE-SIDED AND =2 FOR TWO-SIDED CASES (IN) C Q0 : TRANSITION FROM CORE TO WING (DOPPLER UNITS) (IN) C Q : QUADRATURE POINTS (OUT) C WQ : QUADRATURE WEIGHTS (OUT) C C CODED BY: AKE NORDLUND (OCT-1981). C C REVISED BY MATS CARLSSON (JAN-1984): NEW OPTIONS C QMAX=Q0 Q LINEAR IN FREQUENCY C QMAX OR Q0.LT.0 Q GIVEN AND ONLY WEIGHTS ARE CALCULATED C C: C: FREQL 88-05-03 MODIFICATIONS: (MATS CARLSSON) C: ADDED OPTIONS GAVE INCORRECT WEIGHTS FOR IND=1 C: THIS AFFECTED COOLING RATES FOR LINES WITH FREQUENCY C: QUADRATURE GIVEN WITH QMAX=Q0 OR QMAX.LT.0 OR Q0.LT.0 C: THE CORRECTED WEIGHTS ARE TWICE THE OLD ONES C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CSLINE' C IF(QMAX(KR).EQ.Q0(KR)) THEN DQ=IND(KR)*QMAX(KR)/(NQ(KR)-1) Q(1,KR)=-QMAX(KR)*(IND(KR)-1) DO 50 NY=2,NQ(KR) Q(NY,KR)=Q(NY-1,KR)+DQ WQ(NY,KR)=DQ*2./IND(KR) 50 CONTINUE WQ(1,KR)=0.5*DQ*2./IND(KR) WQ(NQ(KR),KR)=0.5*DQ*2./IND(KR) ELSE IF(QMAX(KR).GE.0.0 .AND. Q0(KR).GE.0.0) THEN C C CONSTANTS C AL10=LOG(10.) A=10.**(Q0(KR)+.5) HALF=0.5 XMAX=LOG10(A*MAX(HALF,QMAX(KR)-Q0(KR)-HALF)) C C ONE-SIDED CASE C IF (IND(KR).EQ.1) THEN DX=XMAX/(NQ(KR)-1) DO 100 J=1,NQ(KR) X=(J-1)*DX X10=10.**X Q(J,KR)=X+(X10-1.)/A WQ(J,KR)=2.*DX*(1.+X10*AL10/A) 100 CONTINUE WQ(1,KR)=WQ(1,KR)*0.5 C C TWO-SIDED CASE C ELSE DX=2.*XMAX/(NQ(KR)-1) DO 110 J=1,NQ(KR) X=-XMAX+(J-1)*DX X10=10.**X Q(J,KR)=X+(X10-1./X10)/A WQ(J,KR)=DX*(1.+(X10+1./X10)*AL10/A) 110 CONTINUE END IF ELSE DO 200 NY=2,NQ(KR)-1 WQ(NY,KR)=0.5*(Q(NY+1,KR)-Q(NY-1,KR)) 200 CONTINUE WQ(1,KR)=0.5*(Q(2,KR)-Q(1,KR)) WQ(NQ(KR),KR)=0.5*(Q(NQ(KR),KR)-Q(NQ(KR)-1,KR)) C C FOR ONE-SIDED CASE, MULTIPLY WEIGHTS BY 2.0 C IF(Q(1,KR).GE.0.0) THEN DO 300 NY=1,NQ(KR) WQ(NY,KR)=WQ(NY,KR)*2.0 300 CONTINUE ENDIF ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE FREQ C C CALCULATES THE FREQUENCY QUADRATURE C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CSLINE' INCLUDE 'CCONST' C C FREQUENCY QUADRATURE IN LINES C DO 200 KR=1,NLINE CALL FREQL(KR) IF(IWIDE(KR)) THEN KT=KTRANS(KR) FRQ(0,KT)=CC/ALAMB(KR)*1.E8 DO 100 NY=1,NQ(KR) FRQ(NY,KT)=FRQ(0,KT)*(1.0+Q(NY,KR)*QNORM/CC*1.E5) 100 CONTINUE ENDIF 200 CONTINUE C C FREQUENCY QUADRATURE IN CONTINUA C DO 300 KR=NLINE+1,NRAD CALL FREQC(KR) 300 CONTINUE C RETURN END C C************************************************************************ C SUBROUTINE GAUSI(K,A,B,AI,XI) C C SUPPLIES QUADRATURE WEIGHTS AND POINTS FOR GAUSSIAN QUADRATURE C BETWEEN A AND B. SOURCE OF DATA: LOWAN, DAVIDS, C LEVENSON, BULL. AMER. MATH. SOC. 48, P 739 (1942). C C K : INTEGRATION ORDER, K BETWEEN 1 AND 10 (IN) C A : LOWER QUADRATURE LIMIT (IN) C B : UPPER QUADRATURE LIMIT (IN) C AI : QUADRATURE WEIGHTS (OUT) C XI : QUADRATURE POINTS (OUT) C C CODED BY : B. GUSTAFSSON (1974). C INCLUDE 'PREC' INCLUDE 'PARAM' C DIMENSION AI(MMU),XI(MMU) DOUBLE PRECISION AP(29),XP(29) INTEGER INDOV(9) DATA AP/1.0D0 *,0.55555555555555D0,.88888888888888D0,.347854845137D0 *,0.65214515486254D0,0.23692688505618D0,0.47862867049936D0, * 0.56888888888888D0,0.17132449237917D0,0.36076157304813D0, * 0.46791393457269D0,0.12948496616887D0,0.27970539148927D0, * 0.38183005050511D0,0.41795918367346D0,0.10122853629037D0, * 0.22238103445337D0,0.31370664587788D0,0.36268378337836D0, * 0.08127438836157D0,0.18064816069485D0,0.26061069640293D0, * 0.31234707704000D0,0.33023935500126D0,0.06667134430868D0, * 0.14945134915058D0,0.21908636251598D0,0.26926671930999D0, * 0.29552422471475D0/ DATA XP/ *0.57735026918962D0,.77459666924148D0,.0D0,0.86113631159405D0, *0.33998104358485D0,.90617984593866D0,.53846931010568D0,.0D0, *0.93246951420315D0,.66120938646626D0,.23861918608319D0, *0.94910791234275D0,.74153118559939D0,.40584515137739D0,.0D0, *0.96028985649753D0,.79666647741362D0,.52553240991632D0, *0.18343464249565D0,.96816023950762D0,.83603110732663D0, *0.61337143270059D0,.32425342340380D0,.0D0,0.97390652851717D0, *0.86506336668898D0,.67940956829902D0,.43339539412924D0, *0.14887433898163D0/ DATA INDOV/1,3,5,8,11,15,19,24,29/ C IF (K.EQ.1) GO TO 7 C KUD=0 FLK=FLOAT(K)/2. K2=K/2 FK=FLOAT(K2) IF (ABS(FLK-FK)-1.E-7) 2,1,1 1 K2=K2+1 KUD=1 2 IOEV=INDOV(K-1) INED=IOEV-K2 C DO 3 I=1,K2 IP=INED+I XI(I)=-XP(IP)*(B-A)*0.5+(B+A)*0.5 3 AI(I)=(B-A)*0.5*AP(IP) C K2=K2+1 DO 4 I=K2,K IP=IOEV+K2-I IF (KUD) 6,6,5 5 IP=IP-1 6 CONTINUE XI(I)=XP(IP)*(B-A)*0.5+(B+A)*0.5 4 AI(I)=(B-A)*0.5*AP(IP) RETURN C 7 XI(1)=(B+A)*0.5 AI(1)=B-A RETURN C END C C********************************************************************* C SUBROUTINE GETWRD(TEXT,K0,K1,K2) C C FINDS NEXT WORD IN TEXT FROM INDEX K0. NEXT WORD IS TEXT(K1:K2) C THE NEXT WORD STARTS AT THE FIRST ALPHANUMERIC CHARACTER AT K0 C OR AFTER. IT ENDS WITH THE LAST ALPHANUMERIC CHARACTER IN A ROW C FROM THE START C INCLUDE 'PREC' INTEGER MSEPAR PARAMETER (MSEPAR=7) CHARACTER*(*) TEXT CHARACTER SEPAR(MSEPAR) INTEGER K0,K1,K2,I,J DATA SEPAR/' ','(',')','=','*','/',','/ C K1=0 DO 400 I=K0,LEN(TEXT) IF(K1.EQ.0) THEN DO 100 J=1,MSEPAR IF(TEXT(I:I).EQ.SEPAR(J)) GOTO 200 100 CONTINUE K1=I C C NOT START OF WORD C 200 CONTINUE ELSE DO 300 J=1,MSEPAR IF(TEXT(I:I).EQ.SEPAR(J)) GOTO 500 300 CONTINUE ENDIF 400 CONTINUE C C NO NEW WORD. RETURN K1=K2=0 C K1=0 K2=0 GOTO 999 C C NEW WORD IN TEXT(K1:I-1) C 500 CONTINUE K2=I-1 C 999 CONTINUE RETURN END C C*********************************************************************** C FUNCTION H1BB (I,J,T) C C DETERMINES H COLLISIONAL RATE COEFFICIENTS UP TO N=50 C C CONVERSION FACTORS: C C CM**-1 TO EV HH * CC / EE C CM**-1 TO K HH * CC / BK C C K TO HZ BK / HH C CM**-1 TO HZ CC C C CONSTANTS: C C EE ELECTRON CHARGE 1.602189E-12 C HH PLANCK CONSTANT (ERG S) C CC VELOCITY OF LIGHT (CM S**-1) C EM ELECTRON REST MASS (G) C UU ATOMIC MASS UNIT (G) C BK BOLTZMANN CONSTANT (ERG K**-1) C PI 3.14159265359 C C C INPUT: C C I LOWER LEVEL C J UPPER LEVEL C T TEMPERATURE (K) C C OUTPUT: C C H1BB UPWARDS COLLISIONAL RATE I TO J C C: H1BB 90-02-28 MODIFICATIONS: (MARTIN J. STIFT) C: REWRITTEN. UP TO 50 BOUND LEVELS + ONE CONTINUUM LEVEL C: INCLUDE 'PREC' C INTEGER I,J,K DOUBLE PRECISION * T,S,X, * AI,AMN,AJ,BM,BMN,DELMN,DERG,ENG, * DION,FIJ,GAMMN,GMX,G0,G1,G2, * EION,EV(50),EE,HH,CC,EM,UU,BK,PI,TEV, * B1S2P(9),B1S2S(9),RH1,RH2,H1BB C C POLYNOMIAL FITS FOR COLLISIONAL RATE COEFFICIENTS (LYMAN ALPHA) C C H(1S) TO H (2P) C DATA B1S2P / -2.814949375869D+1, * 1.009828023274D+1, * -4.771961915818D+0, * 1.467805963618D+0, * -2.979799374553D-1, * 3.861631407174D-2, * -3.051685780771D-3, * 1.335472720988D-4, * -2.476088392502D-6/ C C H(1S) TO H(2S) C DATA B1S2S / -2.833259375256D+1, * 9.587356325603D+0, * -4.833579851041D+0, * 1.415863373520D+0, * -2.537887918825D-1, * 2.800713977946D-2, * -1.871408172571D-3, * 6.986668318407D-5, * -1.123758504195D-6/ C C HYDROGEN MEAN ENERGY LEVELS TAKEN FROM BASHKIN & STONER (1975), PP. 2-3 C DATA EV / 0.000, * 82259.047, * 97492.279, * 102823.882, * 105291.646, * 106632.157, * 107440.441, * 107965.049, * 108324.718, * 41 * 0.000/ C C HYDROGEN IONIZATION LEVEL TAKEN FROM BASHKIN & STONER (1975), P. 3 C DATA EION / 109678.764/ C C MISCELLANOUS CONSTANTS C DATA EE / 1.602189D-12/, * HH / 6.626176D-27/, * CC / 2.99792458D10/, * EM / 9.109534D-28/, * UU / 1.6605655D-24/, * BK / 1.380662D-16/, * PI / 3.14159265359/ C C J MUST BE GREATER THAN I C IF (J.LE.I) CALL STOP('H1BB: ERROR IN ORDER OF ENERGY LEVELS') C AI = FLOAT(I) AJ = FLOAT(J) C C DETERMINE ENERGY LEVELS UP TO N=50 C DO 1 K = 10,50 EV(K) = EION * (1. - 1. / FLOAT(K*K)) 1 CONTINUE C C EV, EION, ENG IN CM**-1 C C ENG ENERGY ASSOCIATED TO TEMPERATURE C ENG = T * BK / HH / CC C C ENERGY DIFFERENCE AND IONISATION ENERGY FROM LOWER LEVEL C DERG = EV(J) - EV(I) DION = EION - EV(I) C C JOHNSON'S (1972) FORMULA FOR OSCILLATOR STRENGTH C SEE JANEV ET AL. (1987), PP. 315-316 C X = 1. - (AI/AJ) * (AI/AJ) C IF (I.EQ.1) THEN G0 = 1.1330 G1 = -0.4059 G2 = 0.07014 ELSE IF (I.EQ.2) THEN G0 = 1.0785 G1 = -0.2319 G2 = 0.02947 ELSE IF (I.GE.3) THEN G0 = 0.9935 + (0.2328 - 0.1296 / AI) / AI G1 = -(0.6282 - (0.5598 - 0.5299 / AI) / AI ) / AI G2 = (0.3887 - (1.1810 - 1.4700 / AI) / AI ) / (AI * AI) ENDIF C GMX = G0 + (G1 + G2 / X) / X C C OSCILLATOR STRENGTH - AGREES WITH VALUES IN JANEV ET AL. (1987), P. 315 C FIJ = 1.96028052 * GMX * (AI/AJ) / (AJ * AJ * X * X * X) C C APN, BP AND BPN OF VRIENS AND SMEETS (1980), EQNS. (11), (13) AND (12) C AMN = 2. * EION * FIJ / DERG C BM = ( 1.4 * LOG(AI) - 0.7 - (0.51 - (1.16 - 0.55 / AI) / * AI ) / AI ) / AI C BMN = 4. * (EION/DERG) * (EION/DERG) / (AJ * AJ * AJ) * * (1. + (4./3.) * (DION/DERG) + * BM * (DION/DERG) * (DION/DERG) ) C C DELTA_PN AND GAMMA_PN OF VRIENS AND SMEETS (1980), EQNS. (18) AND (19) C S = ABS(AI - AJ) DELMN = EXP(-BMN/AMN) + 0.06 * (S/AI) * (S/AI) / AJ GAMMN = EION * (3. + 11. * (S/AI)*(S/AI)) * * LOG(1. + AI * AI * AI * ENG / EION) / * (6. + 1.6 * AJ * S + 0.3 / (S * S) + * 0.8 * AJ * SQRT(AJ/S) * ABS(S-0.6)) C C UPWARD COLLISIONAL RATE (I TO J) : C VRIENS AND SMEETS (1980), EQN. (17), EXCEPT LYMAN ALPHA C LYMAN ALPHA ACCORDING TO JANEV ET AL. (1987) PP. 18-21 AND 257 C C TEV IS IN EV IN POLYNOMIAL EXPRESSIONS C IF (I.EQ.1 .AND. J.EQ.2) THEN C TEV = LOG(T * BK / EE) C RH1 = B1S2P(1) + TEV * (B1S2P(2) + TEV * (B1S2P(3) + * TEV * (B1S2P(4) + TEV * (B1S2P(5) + * TEV * (B1S2P(6) + TEV * (B1S2P(7) + * TEV * (B1S2P(8) + TEV * B1S2P(9) ))))))) RH2 = B1S2S(1) + TEV * (B1S2S(2) + TEV * (B1S2S(3) + * TEV * (B1S2S(4) + TEV * (B1S2S(5) + * TEV * (B1S2S(6) + TEV * (B1S2S(7) + * TEV * (B1S2S(8) + TEV * B1S2S(9) ))))))) C H1BB = EXP(RH1) + EXP(RH2) C ELSE C H1BB = 1.6D-7 * SQRT (EE / HH / CC) * * (AMN * LOG(0.3D0 * ENG / EION + DELMN) + BMN) * * SQRT(ENG) * EXP(-DERG/ENG) / (ENG + GAMMN) C ENDIF C RETURN END C C************************************************************************* C FUNCTION H1BF(I,T) C C DETERMINES H COLLISIONAL IONIZATION COEFFICIENTS UP TO N=50 C C CONVERSION FACTORS: C C CM**-1 TO EV HH * CC / EE C CM**-1 TO K HH * CC / BK C C K TO HZ BK / HH C CM**-1 TO HZ CC C C CONSTANTS: C C EE ELECTRON CHARGE 1.602189E-12 C HH PLANCK CONSTANT (ERG S) C CC VELOCITY OF LIGHT (CM S**-1) C EM ELECTRON REST MASS (G) C UU ATOMIC MASS UNIT (G) C BK BOLTZMANN CONSTANT (ERG K**-1) C PI 3.14159265359 C C C INPUT: C C I LOWER LEVEL C T TEMPERATURE (K) C C OUTPUT: C C H1BF COLLISIONAL IONIZATION RATE I TO K C C................................................................... C C: H1BF 90-02-28 MODIFICATIONS: (MARTIN J. STIFT) C: REWRITTEN. UP TO 50 BOUND LEVELS + ONE CONTINUUM LEVEL C: INCLUDE 'PREC' C INTEGER I,K DOUBLE PRECISION * T,RATIK,AI,ENG,DION, * EION,EV(50),EE,HH,CC,EM,UU,BK,PI,TEV, * B1SI(9),B2SI(9),H1BF C C H1S TO H+ C DATA B1SI / -3.271396786375D+1, * 1.353655609057D+1, * -5.739328757388D+0, * 1.563154982022D+0, * -2.877056004391D-1, * 3.482559773737D-2, * -2.631976175590D-3, * 1.119543953861D-4, * -2.039149852002D-6/ C C H2S TO H+ C DATA B2SI / -1.973476726029D+1, * 3.992702671457D+0, * -1.773436308973D+0, * 5.331949621358D-1, * -1.181042453190D-1, * 1.763136575032D-2, * -1.616005335321D-3, * 8.093908992682D-5, * -1.686664454913D-6/ C C HYDROGEN MEAN ENERGY LEVELS TAKEN FROM BASHKIN & STONER (1975), PP. 2-3 C DATA EV / 0.000, * 82259.047, * 97492.279, * 102823.882, * 105291.646, * 106632.157, * 107440.441, * 107965.049, * 108324.718, * 41 * 0.000/ C C HYDROGEN IONIZATION LEVEL TAKEN FROM BASHKIN & STONER (1975), P. 3 C DATA EION / 109678.764/ C C MISCELLANOUS CONSTANTS C DATA EE / 1.602189E-12/, * HH / 6.626176E-27/, * CC / 2.99792458E10/, * EM / 9.109534E-28/, * UU / 1.6605655E-24/, * BK / 1.380662E-16/, * PI / 3.14159265359/ C AI = FLOAT(I) C C DETERMINE ENERGY LEVELS UP TO N=50 C DO 1 K = 10,50 EV(K) = EION * (1. - 1. / FLOAT(K*K)) 1 CONTINUE C C DION AND ENG IN CM**-1 C C ENG ENERGY ASSOCIATED TO TEMPERATURE C DION IONISATION ENERGY FROM LOWER LEVEL C ENG = T * BK / HH / CC DION = EION - EV(I) C C COLLISIONAL IONIZATION RATE (I TO K) : VRIENS AND SMEETS (198), EQN. (8) C IF (I.EQ.1) THEN C TEV = LOG(T * BK / EE) RATIK = B1SI(1) + TEV * (B1SI(2) + TEV * (B1SI(3) + * TEV * (B1SI(4) + TEV * (B1SI(5) + * TEV * (B1SI(6) + TEV * (B1SI(7) + * TEV * (B1SI(8) + TEV * B1SI(9) ))))))) H1BF = EXP(RATIK) C ELSE IF (I.EQ.2) THEN C TEV = LOG(T * BK / EE) RATIK = B2SI(1) + TEV * (B2SI(2) + TEV * (B2SI(3) + * TEV * (B2SI(4) + TEV * (B2SI(5) + * TEV * (B2SI(6) + TEV * (B2SI(7) + * TEV * (B2SI(8) + TEV * B2SI(9) ))))))) H1BF = EXP(RATIK) C ELSE C H1BF = 9.56D-6 * EXP(-DION/ENG) / ENG / SQRT(ENG) * * SQRT(EE / HH / CC) * (EE / HH / CC) / * ( (DION/ENG)**2.33 + 4.38*(DION/ENG)**1.72 + 1.32*(DION/ENG) ) C ENDIF C RETURN END C C************************************************************************* C SUBROUTINE HCOL C C COLLISIONAL RATES FOR HYDROGEN. C: C: HCOL 90-02-28 MODIFCATIONS: (MARTIN J. STIFT) C: REWRITTEN. UP TO 50 BOUND LEVELS + ONE CONTINUUM LEVEL C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CCONST' C DOUBLE PRECISION T,H1BB,H1BF C C CALCULATE C BOUND-BOUND AND BOUND-FREE C REFERENCE: VRIENS AND SMEETS, PHYS. REV. A, 22, 940 (1980) C JANEV ET AL. (1987) C ELEMENTARY PROCESSES IN HYDDROGEN-HELIUM PLASMAS C DO 4 K = 1,NDEP C T = TEMP(K) C DO 2 J = 2,NK-1 DO 1 I = 1,J-1 C(I,J,K) = NE(K) * H1BB (I,J,T) C(J,I,K) = C(I,J,K) * NSTAR(I,K) / NSTAR(J,K) 1 CONTINUE 2 CONTINUE C DO 3 I = 1,NK-1 C(I,NK,K) = NE(K) * H1BF (I,T) C(NK,I,K) = C(I,NK,K) * NSTAR(I,K) / NSTAR(NK,K) 3 CONTINUE C 4 CONTINUE C RETURN END C C************************************************************************* C SUBROUTINE INCDNT C C READS INCIDENT RADIATION FIELD FROM FILE IMINUS C FORMAT OF FILE IS C SCALEI GLOBAL SCALING FACTOR C NMUI NUMBER OF ANGLE POINTS C NXL NUMBER OF WAVELENGTH POINTS C XL(I) I-(I,1)...I-(I,NMUI) WAVELENGTH (ANGSTROM), I- C ... C ONLY NMUI=1 WITH I- INDEPENDENT OF MU HAS BEEN IMPLEMENTED C ONLY B-F TRANSITIONS (PHOTOIONIZATION) ARE CONSIDERED C "EXACT" MATCH WITH WAVELENGTH IS NEEDED FOR IMINUS TO BE SET C THIS IS BECAUSE IMINUS FOR MANY APPLICATIONS CONSIST OF C CORONAL RADIATION IN DISCRETE LINES INTEGRATED OVER C WAVELENGTH BINS. IT IS THEN NECESSARY TO REMAP THE BINS TO C THE WAVELENGTH BINS USED FOR THE TRANSITION C C: INCDNT 92-10-08 NEW ROUTINE: (MATS CARLSSON) C: READS INCIDENT RADIATION FIELD FROM FILE IMINUS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLU' INCLUDE 'CIMIN' C PARAMETER (MXL=1000) DIMENSION XL0(MXL),XIMIN0(MXL) LOGICAL INCI C C READ IMINUS FILE C CALL CSTRIP(LIMIN,'IMINUS') C READ(LDUMS,*) SCALEI READ(LDUMS,*) NMUI IF(NMUI.NE.1) CALL STOP('INCDNT: NMUI.NE.1') READ(LDUMS,*) NXL IF(NXL.GT.MXL) CALL STOP('INCDNT: NXL.GT.MXL') DO 100 I=1,NXL READ(LDUMS,*) XL0(I),XIMIN0(I) 100 CONTINUE C C INITIALIZE XIMIN TO ZERO C DO 280 KR=1,NRAD DO 240 NY=1,NQ(KR) DO 200 MU=1,NMU XIMIN(NY,MU,KR)=0.0 200 CONTINUE 240 CONTINUE 280 CONTINUE C C FIND WAVELENGTH COINCIDENCES AND LOAD XIMIN C FIRST WAVELENGTH FIND: WRITE TRANSITION INFORMATION TO OUT C ALL FOLLOWING WAVELENGTHS: WRITE WAVELENGTH C DO 480 KR=NLINE+1,NRAD INCI=.FALSE. KT=KTRANS(KR) DO 440 NY=1,NQ(KR) XL=CC*1.E8/FRQ(NY,KT) DO 400 IXL=1,NXL IF(ABS(XL-XL0(IXL)).LT.0.1) THEN DO 300 MU=1,NMU XIMIN(NY,MU,KR)=XIMIN0(IXL)*SCALEI 300 CONTINUE IF(IWATOM.NE.0) THEN IF(.NOT.INCI) THEN I=IRAD(KR) J=JRAD(KR) WRITE(LOUT,310) KR,LABEL(I),LABEL(J) 310 FORMAT(' INCIDENT RADIATION FIELD IN TRANSITION',I3/ * 1X,A,' -> ',A/' WAVELENGTHS:') WRITE(LOUT,320) XL 320 FORMAT(1X,F20.3) ELSE WRITE(LOUT,320) XL ENDIF ENDIF INCI=.TRUE. GOTO 440 ENDIF 400 CONTINUE 440 CONTINUE 480 CONTINUE C RETURN END SUBROUTINE INITIA C C CALCULATES A STARTING SOLUTION C ISTART=-1 START FROM FILE RSTRT C 0 I=0 C 1 N=NSTAR C .GT.1 ISTART-1 ESCAPE PROBABILITY ITERATIONS FROM N=NSTAR C ILAMBD NUMBER OF LAMBDA ITERATIONS FROM START SOLUTION C: C: INITIA 87-05-21 MODIFICATIONS: (MATS CARLSSON, PHILIP JUDGE) C: ESCAPE PROBABILITY OPTION ADDED C: C: 88-06-22 MODIFICATIONS: (MATS CARLSSON) C: ICONV IS SET TO 0 DURING SCATTERING LAMBDA ITERATIONS TO C: AVOID PRINTOUT IN WTEST CALLED FROM TRPT C: C: 92-10-08 MODIFICATIONS: (MATS CARLSSON) C: CALLS INCDNT TO GET INCIDENT RADIATION FIELD IF C: ITRAN.GE.10 C: C: 95-03-15 MODIFICATIONS: (MATS CARLSSON) C: IF ISCAT=1 THERE ARE NO SCATTERING LAMBDA ITERATIONS C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION SWITCHED ON BY INCRAD.NE.0 INSTEAD OF C: ITRAN=10-14 C: SCATTERING VERSION OF TRANSP USED INSTEAD OF C: SCATTERING LAMBDA ITERATIONS C: C: 97-09-24 MODIFICATIONS: (MATS CARLSSON) C: DPTYPE='H' ADDED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLU' C PARAMETER (INDEP=400) DIMENSION DPIN(INDEP),ANIN(MK,INDEP),DP(MDEP) CHARACTER*72 TEXT DATA LMAX/20/ C CALL CPUTIME('INITIA',0,0,1) C C UP TO 9 ESCAPE PROBABILITY ITERATIONS C IF(ABS(ISTART) .GT. 10) THEN CALL STOP('INITIA: ISTART NOT KNOWN') ENDIF C C READ INCIDENT RADIATION FIELD FROM FILE IMINUS IF INCRAD.NE.0 C IF(INCRAD.NE.0) CALL INCDNT C IF(ISTART.EQ.-1) THEN C C READ START APPROXIMATION FROM FILE RSTRT C CALL CSTRIP(LRSTRT,'RSTRT') READ(LDUMS,100,END=900) ATMOID 100 FORMAT(A) READ(LDUMS,100) TEXT CALL LJUST(TEXT) IF(DPTYPE.NE.TEXT(1:1)) THEN CALL STOP('INITIA: DIFFERENT SCALE TYPES IN DSCALE AND RSTRT') ENDIF READ(LDUMS,*) GDUM READ(LDUMS,*) KDEP IF(KDEP.GT.INDEP) CALL STOP('INITIA:KDEP.GT.INDEP') READ(LDUMS,*) (DPIN(K),DUM,DUM,DUM,DUM,K=1,KDEP) READ(LDUMS,*) ((ANIN(I,K),I=1,NK),K=1,KDEP) C C STORE DEPTH-GRID USED IN DP C IF(DPTYPE.EQ.'M') THEN DO 110 K=1,NDEP DP(K)=LOG10(CMASS(K)) 110 CONTINUE ELSE IF(DPTYPE.EQ.'T') THEN DO 120 K=1,NDEP DP(K)=LOG10(TAU(K)) 120 CONTINUE ELSE IF(DPTYPE.EQ.'H') THEN DO 125 K=1,NDEP DP(K)=-HEIGHT(K) 125 CONTINUE DO 127 K=1,KDEP DPIN(K)=-DPIN(K) 127 CONTINUE ENDIF C C INTERPOLATE START APPROXIMATION TO DEPTH-GRID USED, C LINEAR INTERPOLATION IN THE LOG OF THE POPULATION NUMBERS C C TAKE THE LOG OF THE INPUT POPULATION NUMBERS C DO 140 K=1,KDEP DO 130 I=1,NK ANIN(I,K)=LOG(ANIN(I,K)) 130 CONTINUE 140 CONTINUE C IF(DPIN(1).GT.DP(1)) WRITE(LOUT,150) DPIN(1),DP(1) 150 FORMAT(' EXTRAPOLATION IN MIN DEPTH, DPIN=',F13.7,' DP=',F13.7) IF(DPIN(KDEP).LT.DP(NDEP)) WRITE(LOUT,160) DPIN(KDEP),DP(NDEP) 160 FORMAT(' EXTRAPOLATION IN MAX DEPTH, DPIN=',F13.7,' DP=',F13.7) L=1 DO 200 K=1,NDEP 170 L=L+1 IF(DP(K).GT.DPIN(L) .AND. L.LT.KDEP) GOTO 170 APOL=(DP(K)-DPIN(L-1))/(DPIN(L)-DPIN(L-1)) DO 180 I=1,NK N(I,K)=EXP(ANIN(I,L-1)+APOL*(ANIN(I,L)-ANIN(I,L-1))) 180 CONTINUE L=L-1 200 CONTINUE C C GET THE SCATTERING CONTRIBUTION TO THE SOURCE FUNCTION CONSISTENT C WITH THE RADIATION FIELD C ISCAT0=ISCAT ICONV0=ICONV ISCAT=1 ICONV=0 CALL TRPT ISCAT=ISCAT0 ICONV=ICONV0 C ELSE IF(ISTART.EQ.0 .OR. ISTART.LT.-1) THEN DO 340 KR=1,NRAD IF(KR.LE.NLINE) THEN DO 310 K=1,NDEP RJI(K,KR)=A(KR) RIJ(K,KR)=0. 310 CONTINUE ELSE DO 315 K=1,NDEP RJI(K,KR)=0.0 RIJ(K,KR)=0.0 315 CONTINUE KT=KTRANS(KR) I=IRAD(KR) J=JRAD(KR) DO 330 NY=1,NQ(KR) HN3C2=2.*HH*FRQ(NY,KT)/CC *FRQ(NY,KT)/CC *FRQ(NY,KT) WQMUH=WQ(NY,KR)/HNY4P/FRQ(NY,KT)*FRQ(0,KT) DO 320 K=1,NDEP GIJ(K)=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FRQ(NY,KT)/BK/TEMP(K)) ALFA(K)=ALFAC(NY,KR-NLINE) RJI(K,KR)=RJI(K,KR)+WQMUH*ALFA(K)*GIJ(K)*HN3C2 320 CONTINUE 330 CONTINUE ENDIF 340 CONTINUE CALL STATEQ(ISUM,1) C ELSE IF(ISTART.GE.1) THEN DO 410 K=1,NDEP DO 400 I=1,NK N(I,K)=NSTAR(I,K) 400 CONTINUE 410 CONTINUE ENDIF C C ADDED BY P.G. JUDGE, MAY 1987 C ISTART-1 ESCAPE PROBABILITY ITERATIONS. INSTEAD OF A FULL SOLN C TO THE RADIATIVE TRANSFER EQUATIONS THE FIRST-ORDER ESCAPE PROBABILITY C APPROXIMATION IS USED. C DO 600 IT=1,ABS(ISTART)-1 CALL ESCAPE CALL STATEQ(ISUM,1) 600 CONTINUE C* C* 06-03-89 PHILIP JUDGE ALTERATION: C* OUTPUT ESCAPE PROBABILITIES C* C* IF (ABS(ISTART)-1 .GE. 1)THEN C* CALL OPEN(LESCA,'ESCAPE',1,'NEW') C* DO 2001 KR=1,NLINE C* DO 201 K=1,NDEP C* WRITE(LESCA,2002)RJI(K,KR)/A(KR) C* 201 CONTINUE C* 2001 CONTINUE C* 2002 FORMAT(1PE11.4) C* CALL CLOSE(LESCA) C* ENDIF C* C* 06-03-89 END C* C C ILAMBD LAMBDA ITERATIONS C DO 500 IT=1,ILAMBD CALL TRPT CALL STATEQ(ISUM,1) 500 CONTINUE CALL CPUTIME('INITIA',0,1,1) RETURN 900 CALL STOP('INITIA: RSTRT NOT EXISTING') END C C*********************************************************************** C SUBROUTINE JFIX(KF,IPR) C C COMPUTES FIXED RADIATIVE RATES FOR A PHOTO-EXCITED TRANSITION C FROM SAVED MEAN INTENSITIES C READS JNU DATA WRITTEN PREVIOUSLY BY ROUTINE WRJFIX C AND COMPUTES FIXED RATES USING THESE C VALUES OF JNU: C C SECOND VERSION FOLLOWING COMMENTS OF PGJ AND MC C THIS VERSION READS ALL THE JNY DATA ONCE, SORTS AND STORES THEM C -THIS ALLOWS MULTIPLE CALLS TO BE EFFICIENT AT THE EXPENSE OF C -MEMORY REQUIREMENTS C C PHOTOEXCITATION RATE IS THUS: C R( L - U ) = INTEGRAL (JNU * B( L - U ) * PHI(NU)) DNU FOR A LINE C = INTEGRAL (JNU * ALPHA/HNY4P) DNU FOR A CONTINUUM C C INPUT: C KF = TRANSITION USED C IPR = PRINT OPTION: .GT.0 OUTPUT FREQUENCIES USED C C OUTPUT C CJI AND CIJ WITH RATES ADDED. C C CODED BY P.G. JUDGE, SEPTEMBER 7TH, 1987; DECEMBER 3RD, 1987 C C COMPLETELY REWRITTEN BY P.G. JUDGE AND M. CARLSSON APRIL-MAY 1988 C: C: JFIX 88-05-06 NEW ROUTINE: (MATS CARLSSON, PHILIP JUDGE) C: CALCULATES FIXED RATES FROM ATOMIC PARAMETERS READ BY C: ATOM INTO COMMON BLOCK CFIX AND A FIXED RADIATION FIELD C: READ FROM FILE JFIX C: C: 90-08-17 MODIFICATIONS: (MATS CARLSSON) C: JFIX FILE IS FIRST OPENED TO READ HEADER WITH INFORMATION C: ON RECORD LENGTH, THEN REOPENED WITH CORRECT RECORD C: LENGTH C: C: 97-09-24 MODIFICATIONS: (MATS CARLSSON) C: DPTYPE='H' ADDED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' INCLUDE 'CFIX' C C LOCALS: C DIMENSION WKQ(MQJ),WKK(MDEPJ),DP(MDEP),APOLK(MDEP),APOLQ(2*MQJ-1), * FREQF(2*MQ+1),QQ(2*MQJ+1) INTEGER INDX(MQJ),ICALL,NDEPJ,NQJ,KPOLK(MDEP),KPOLQ(2*MQJ-1) LOGICAL SORT DATA ICALL/0/,SORT /.FALSE./ SAVE ICALL C C INITIALIZE ATOMIC PARAMETERS. C TRANSITION I-J HAS LINE-CENTRE (OR CONTINUUM LIMIT) FREQUENCY Q00: C CGS THROUGHOUT: C ATOMIC PARAMETERS: C J = JFX(KF) I = IFX(KF) C C WAVEJ IS LAMBDA IN CM C Q00 IS FREQUENCY AT LINE CENTRE (OR CONTINUUM LIMIT) C WAVEJ = HCE/(EV(J)-EV(I))/1.E8 Q00 = CC/WAVEJ C C LINE CASE: C IF(IPHO(KF) .EQ. 0) THEN AEINST = FF(KF)* 6.67E-01 * G(I)/G(J)/WAVEJ/WAVEJ BDN = (WAVEJ*1.E8)*(WAVEJ*1.E8)*(WAVEJ*1.E8)/HC2*AEINST BUP = G(J)/G(I)*BDN DQ0=QNORM*1.E5/CC *QF(NQF(KF),KF)*Q00 QMINJ=Q00 - DQ0 QMAXJ=Q00 + DQ0 GIJK=G(I)/G(J) HN3C2=AEINST/BDN C C CONTINUUM CASE C ELSE IF (IPHO(KF) .EQ. 1) THEN QMINJ=Q00 QMAXJ=Q00*(1.+QNORM*1.E5/CC *QF(NQF(KF),KF)) ELSE CALL STOP(' JFIX: IPHO(KF) .NE. 0 OR 1') ENDIF C C FIRST CALL: READ THE JNYJ AND QJ DATA IN FILE JFIX. C IF (ICALL .EQ. 0) THEN ICALL=ICALL+1 C C CALL ROUTINE OPEN WHICH FINDS RECORD LENGTH NDEPJ C DIRECT ACCESS FILE JFIX CONTAINS C 1. A HEADER RECORD WITH NUMBER OF FREQUENCIES AND DEPTH POINTS C 2. RECORD/RECORDS WITH FREQUENCIES IN HZ C 3. RECORD WITH LG(TAU 5000) C 4. RECORD WITH LG(COL MASS) C 5. RECORDS WITH LOG(JNY(1:NDEPF)/BNU(1:NDEPF)) C ONE FOR EACH FREQUENCY C C NOTE THAT ALL INTERPOLATION IS PERFORMED IN FREQUENCY. IF IT C IS NECESSARY TO STORE JNY AT TWO POINTS VERY CLOSE IN FREQUENCY C (A RELATIVE DELTA NU OR LAMBDA ON THE ORDER OF 1.E-6) IT MIGHT C BE NECESSARY TO USE DOUBLE PRECISION FOR THE FREQUENCIES. SUCH C A CHANGE WOULD AFFECT SEVERAL ROUTINES. C IREC=1 C C READ HEADER RECORD TO GET RECORD LENGTH C NDEPJ=2 CALL OPEN(LJFIX,'JFIX',NDEPJ,'OLD') write(*,*) ndepj READ(LJFIX,REC=IREC) NQJ,NDEPJ CALL CLOSE(LJFIX) CALL OPEN(LJFIX,'JFIX',NDEPJ,'OLD') write(*,*) ndepj IF(NQJ.GT.MQJ) CALL STOP('JFIX: NQJ.GT.MQJ') IF(NDEPJ.GT.MDEPJ) CALL STOP('JFIX: NDEPJ.GT.MDEPJ') C C READ FREQUENCIES AND LOG DEPTH SCALES C DO 100 IQ0=1,NQJ,NDEPJ IREC = IREC + 1 READ(LJFIX,REC=IREC) * (QJ(IQ),IQ=IQ0,MIN(IQ0+NDEPJ-1,NQJ)) 100 CONTINUE IREC=IREC+1 READ(LJFIX,REC=IREC) (TAU5J(K),K=1,NDEPJ) IREC=IREC+1 READ(LJFIX,REC=IREC) (CMASSJ(K),K=1,NDEPJ) IRECQ0=IREC C C CHECK TO SEE IF THE ARRAYS NEED SORTING... C DO 200 IQ=2,NQJ IF(QJ(IQ).LT.QJ(IQ-1)) THEN SORT=.TRUE. GOTO 205 ENDIF 200 CONTINUE C C SORT THE QJ VALUES: C THE INDEXING USED FOR SORTING ARE RETURNED IN INDX C INDX SHOWS THE RELATION BETWEEN THE QJ VALUES BEFORE SORTING: C QJ(INDX(1)) .LT. QJ(INDX(2)) .LT. QJ(INDX(3)) ... C USES HEAPSORT ALGORITHM (NUMERICAL RECIPIES P231) C 205 CONTINUE IF(SORT) THEN CALL SORTJ(NQJ,QJ,WKQ,INDX) ELSE DO 210 IQ=1,NQJ INDX(IQ)=IQ 210 CONTINUE ENDIF C C STORE INTERPOLATION ABSCISSAS C FIND LINEAR INTERPOLATION WEIGHTS AND INDICES C K0 IS INDEX OF FIRST DEPTH-POINT WITHIN JFIX GRID C K1 IS INDEX OF LAST DEPTH-POINT WITHIN JFIX GRID C IF(DPTYPE.EQ.'T' .OR. DPTYPE.EQ.'H') THEN DO 250 K=1,NDEP DP(K)=LOG10(TAU(K)) 250 CONTINUE CALL IPOLW(TAU5J,NDEPJ,DP,NDEP,APOLK,KPOLK,K0,K1) ELSE DO 270 K=1,NDEP DP(K)=LOG10(CMASS(K)) 270 CONTINUE CALL IPOLW(CMASSJ,NDEPJ,DP,NDEP,APOLK,KPOLK,K0,K1) ENDIF C C READ JNY DATA SAVED FROM FILE. (ENTIRE SET) C READ IN ORDER OF INCREASING QJ INTO TEMPORARY VARIABLE WKK C INTERPOLATE ONTO DEPTH-GRID USED AND STORE IN JNYJ C IF DEPTH-GRID USED EXTENDS ABOVE JNY GRID, USE TOPMOST VALUE OF JNY C IF GRID EXTENDS BELOW, USE PLANCK FUNCTION C INTERPOLATION IN JNY OR LOG(JNY) WOULD LEAD TO SERIOUS ERRORS AT C THE BOTTOM OF THE ATMOSPHERE. INTERPOLATION IS THEREFORE PERFORMED C IN LOG(JNY/BNU). THIS WILL GIVE GOOD ACCURACY AT THE BOTTOM BUT C MAY LEAD TO INTERPOLATION ERRORS AT THE TOP WHERE JNY IS MORE OR C LESS CONSTANT. C DO 350 IQ=1,NQJ IREC=IRECQ0+INDX(IQ) READ(LJFIX,REC=IREC) (WKK(K),K=1,NDEPJ) IF(K0.GE.2) THEN BNUK0=PLANCK(QJ(IQ),TEMP(K0)) ENDIF DO 300 K=1,K0-1 BNUK=PLANCK(QJ(IQ),TEMP(K)) JNYJ(K,IQ)=WKK(KPOLK(K0)) + APOLK(K0)* * ( WKK(KPOLK(K0)+1)-WKK(KPOLK(K0)) ) + LOG( BNUK0/BNUK ) 300 CONTINUE DO 310 K=K0,K1 JNYJ(K,IQ)=WKK(KPOLK(K)) + APOLK(K)* * ( WKK(KPOLK(K)+1)-WKK(KPOLK(K)) ) 310 CONTINUE DO 320 K=K1+1,NDEP JNYJ(K,IQ)=0.0 320 CONTINUE 350 CONTINUE CALL CLOSE(LJFIX) ENDIF C C INTERPOLATE (LINEARLY) C AND COMPUTE THE CONTRIBUTIONS TO THE J-INTEGRAL USING WEIGHTS C FROM FREQLF AND FREQCF C FOLLOW THE ROUTINE TRPT WHERE THE DETAILED TRANSITIONS ARE COMPUTED C C LINE CASE C COMPUTE PROFILE: (FOLLOW ROUTINE PROFIL- IGNORE MU-DEPENDENCE), C USE FULL PROFILE INTEGRATION C IF(IPHO(KF) .EQ. 0) THEN C C LINE CASE C FIND INTERPOLATION WEIGHTS AND INDICES C PRINT FREQUENCIES TO OUTPUT FILE C DO 600 NY=1,2*NQF(KF)-1 IF(NY.LE.NQF(KF)-1) THEN QQ(NY)=-QF(NQF(KF)-NY+1,KF) ELSE QQ(NY)=QF(NY-NQF(KF)+1,KF) ENDIF FREQF(NY)=Q00*(1.+QNORM*1.E5/CC *QQ(NY)) 600 CONTINUE CALL IPOLW(QJ,NQJ,FREQF,2*NQF(KF)-1,APOLQ,KPOLQ,IQ0,IQ1) IF(IPR.GT.0) THEN WRITE(LOUT,800) KF,'PHOTOEXCITATION' WRITE(LOUT,810) 'LINE:' WRITE(LOUT,820) (NY,FREQF(NY),CONVL(CC/FREQF(NY)*1.E8), * NY=1,2*NQF(KF)-1) WRITE(LOUT,810) 'FROM FILE JFIX:' WRITE(LOUT,820) (IQ,QJ(IQ),CONVL(CC/QJ(IQ)*1.E8), * IQ=KPOLQ(1),KPOLQ(2*NQF(KF)-1)+1) ENDIF C DO 690 K=1,NDEP RUPF=0.0 RDNF=0.0 WPHIF=0.0 DO 620 NY=1,2*NQF(KF)-1 IF(NY.LE.NQF(KF)-1) THEN WQFF=WQF(NQF(KF)-NY+1,KF) ELSE WQFF=WQF(NY-NQF(KF)+1,KF) ENDIF V = (QQ(NY) - VEL(K))/DNYD(K) CALL VOIGT(ADAMPF(K,KF),V,H) PHIF = H/DNYD(K)/1.772453851 WPHIF=WPHIF+WQFF*PHIF ALFAF=BUP*PHIF C C INTERPOLATION/EXTRAPOLATION IN FREQUENCY C RESU = JNYJ(K,KPOLQ(NY)) + APOLQ(NY)* * ( JNYJ(K,KPOLQ(NY)+1) - JNYJ(K,KPOLQ(NY)) ) RESU = EXP(RESU)*PLANCK(FREQF(NY),TEMP(K)) C C NOW COMPUTE THE FIXED RATE EQUATIONS: C RUPF = RUPF + WQFF*ALFAF*RESU RDNF = RDNF + WQFF*ALFAF*(RESU + HN3C2)*GIJK 620 CONTINUE C(I,J,K) = C(I,J,K) + RUPF/WPHIF C(J,I,K) = C(J,I,K) + RDNF/WPHIF 690 CONTINUE ELSE C C CONTINUUM CASE C FIND INTERPOLATION WEIGHTS AND INDICES C PRINT FREQUENCIES TO OUTPUT FILE C DO 700 NY=1,NQF(KF) FREQF(NY)=Q00*(1.+QNORM*1.E5/CC *QF(NY,KF)) 700 CONTINUE CALL IPOLW(QJ,NQJ,FREQF(1),NQF(KF),APOLQ(1),KPOLQ(1),IQ0,IQ1) IF(IPR.GT.0) THEN WRITE(LOUT,800) KF,'PHOTOIONIZATION' WRITE(LOUT,810) 'CONTINUUM:' WRITE(LOUT,820) (NY,FREQF(NY),CONVL(CC/FREQF(NY)*1.E8), * NY=1,NQF(KF)) WRITE(LOUT,810) 'FROM FILE JFIX:' WRITE(LOUT,820) (IQ,QJ(IQ),CONVL(CC/QJ(IQ)*1.E8), * IQ=KPOLQ(1),KPOLQ(NQF(KF))+1) ENDIF C DO 790 K=1,NDEP RUPF=0.0 RDNF=0.0 DO 720 NY=1,NQF(KF) HN3C2=2.*HH*FREQF(NY)/CC *FREQF(NY)/CC *FREQF(NY) GIJK=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FREQF(NY)/BK/TEMP(K)) ALFAF=ALFACF(NY,KF) WQFF=WQF(NY,KF)/HNY4P/FREQF(NY)*Q00 C C INTERPOLATION/EXTRAPOLATION IN FREQUENCY C RESU = JNYJ(K,KPOLQ(NY)) + APOLQ(NY)* * ( JNYJ(K,KPOLQ(NY)+1) - JNYJ(K,KPOLQ(NY)) ) RESU = EXP(RESU)*PLANCK(FREQF(NY),TEMP(K)) C C NOW COMPUTE THE FIXED RATE EQUATIONS: C RUPF = RUPF + WQFF*ALFAF*RESU RDNF = RDNF + WQFF*ALFAF*(RESU + HN3C2)*GIJK 720 CONTINUE C(I,J,K) = C(I,J,K) + RUPF C(J,I,K) = C(J,I,K) + RDNF 790 CONTINUE ENDIF C 800 FORMAT(//' JFIX: FIXED TRANSITION NR ',I4,4X,A) 810 FORMAT(/1X,A/' NY FREQUENCY LAMBDA(A)') 820 FORMAT(I4,1P,E15.7,0P,F15.3) C RETURN END C C*********************************************************************** C SUBROUTINE FREQLF(KF) C C FREQUENCY POINTS AND CORRESPONDING WEIGHTS, FIXED TRANSITIONS C BASED ON FREQL. THE POINTS ARE CALCULATED ONLY FOR ONE SIDE, C QF(1,KF)=0.0. THE WEIGHTS ARE GIVEN FOR INTEGRATION OVER C NY=-NQF(KF),NQF(KF) WITH QF(-NY,KF)=-QF(NY,KF) C SEE FREQL FOR DESCRIPTION OF ALGORITHM C C NQF : TOTAL NUMBER OF QUADRATURE POINTS (IN) C QMAXF : MAXIMUM VALUE OF FREQUENCY (DOPPLER UNITS) (IN) C Q0F : TRANSITION FROM CORE TO WING (DOPPLER UNITS) (IN) C QF : QUADRATURE POINTS (OUT) C WQF : QUADRATURE WEIGHTS (OUT) C C: C: FREQLF 88-05-04 NEW ROUTINE: (MATS CARLSSON) C: CALCULATES FREQUENCIES FOR FIXED B-B TRANSITIONS (ITRAD=4) C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CFIX' C IF(QMAXF(KF).EQ.Q0F(KF)) THEN DQ=QMAXF(KF)/(NQF(KF)-1) QF(1,KF)=0.0 DO 50 NY=2,NQF(KF) QF(NY,KF)=QF(NY-1,KF)+DQ WQF(NY,KF)=DQ 50 CONTINUE WQF(NQF(KF),KF)=0.5*DQ ELSE IF(QMAXF(KF).GE.0.0 .AND. Q0F(KF).GE.0.0) THEN C C CONSTANTS C AL10=LOG(10.) A=10.**(Q0F(KF)+.5) HALF=0.5 XMAX=LOG10(A*MAX(HALF,QMAXF(KF)-Q0F(KF)-HALF)) C DX=XMAX/(NQF(KF)-1) DO 100 J=1,NQF(KF) X=(J-1)*DX X10=10.**X QF(J,KF)=X+(X10-1.)/A WQF(J,KF)=DX*(1.+X10*AL10/A) 100 CONTINUE C ELSE DO 200 NY=2,NQF(KF)-1 WQF(NY,KF)=0.5*(QF(NY+1,KF)-QF(NY-1,KF)) 200 CONTINUE WQF(1,KF)=(QF(2,KF)-QF(1,KF)) WQF(NQF(KF),KF)=0.5*(QF(NQF(KF),KF)-QF(NQF(KF)-1,KF)) ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE FREQCF(KF) C C GIVES QUADRATURE POINTS AND WEIGHTS FOR CONTINUA, FIXED TRANSITIONS C THE POINTS ARE EQUIDISTANT IN FREQUENCY, THE WEIGHTS ARE C TRAPEZOIDAL C CALCULATES CROSSECTIONS FROM NY**-3 DEPENDENCE IF NOT GIVEN C EXPLICITLY C C QMAX .LT. 0 : Q AND ALFAC GIVEN (FROM ATOM INPUT). ONLY C WEIGHTS CALCULATED AND Q TRANSFORMED FROM ANGSTROM TO DOPPLER C UNITS. C C NQF TOTAL NUMBER OF QUADRATURE POINTS (IN) C QMAXF MINIMUM VALUE OF WAVELENGTH (ANGSTROM) (IN) C QF QUADRATURE POINTS (DOPPLER UNITS) (OUT) C WQF QUADRATURE WEIGHTS (OUT) C C: C: FREQCF 88-05-04 NEW ROUTINE: (MATS CARLSSON) C: CALCULATES FREQUENCIES FOR FIXED B-F TRANSITIONS (ITRAD=4) C: C: 90-03-28 MODIFICATIONS: (MATS CARLSSON) C: CORRECTED MISSPELLING OF VARIABLE NAME FRQ00 C: AND INDEXING OF ALFACF IN HYDROGENIC OPTION C: (ERRORS GAVE DIVISION BY ZERO AND ERRONEOUS INDEXING) C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CCONST' INCLUDE 'CFIX' C WAVE00=HCE/(EV(JFX(KF))-EV(IFX(KF))) FRQ00=CC/WAVE00*1.E8 IF(QMAXF(KF).GE.0.0) THEN FRQMAX=CC/QMAXF(KF)*1.E8 DFRQ=(FRQMAX-FRQ00)/(NQF(KF)-1) DO 100 NY=1,NQF(KF) FRQF=FRQ00+(NY-1)*DFRQ QF(NY,KF)=CC *1.E-5/QNORM*(FRQF/FRQ00-1.0) ALFACF(NY,KF)=FF(KF)*(FRQ00/FRQF)**3 100 CONTINUE ELSE DO 200 NY=1,NQF(KF) FRQF=CC/QF(NY,KF)*1.0E8 QF(NY,KF)=CC *1.0E-5/QNORM*(FRQF/FRQ00-1.0) 200 CONTINUE ENDIF C WQF(1,KF)=0.5*(QF(2,KF)-QF(1,KF)) DO 300 NY=2,NQF(KF)-1 WQF(NY,KF)=0.5*(QF(NY+1,KF)-QF(NY-1,KF)) 300 CONTINUE WQF(NQF(KF),KF)=0.5*(QF(NQF(KF),KF)-QF(NQF(KF)-1,KF)) C RETURN END C C*********************************************************************** C SUBROUTINE SORTJ(NIN,X1,W1,INDX) C C SORTS THE INPUT FREQUENCIES IN INCREASING ORDER C FOR SUBSEQUENT INTERPOLATION C C INPUT: NIN,X1 C ARRAY X1 IS SORTED INTO INCREASING ORDER C C: C: SORTJ 88-05-04 NEW ROUTINE: (PHILIP JUDGE) C: SORTS ARRAY IN INCREASING ORDER AND RETURNS INDEXING C: USED IN (ITRAD=4) OPTION C: INCLUDE 'PREC' DIMENSION X1(NIN) DIMENSION W1(NIN) INTEGER INDX(NIN) C CALL INDEXX(NIN,X1,INDX) C DO 11 J=1,NIN W1(J)=X1(J) 11 CONTINUE DO 12 J=1,NIN X1(J)=W1(INDX(J)) 12 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE INDEXX(N,ARRIN,INDX) C C INDEXES AND ARRAY OF LENGTH N, KEY IN ARRAY INDX. HEAPSORT C ALGORITHM C C: C: INDEXX 88-05-04 NEW ROUTINE: (NUMERICAL RECIPES) C: HEAPSORT ALGORITHM C: USED IN (ITRAD=4) OPTION C: INCLUDE 'PREC' DIMENSION ARRIN(N) INTEGER INDX(N) DO 11 J=1,N INDX(J)=J 11 CONTINUE L=N/2+1 IR=N 10 CONTINUE IF (L.GT.1) THEN L=L-1 INDXT=INDX(L) Q=ARRIN(INDXT) ELSE INDXT=INDX(IR) Q=ARRIN(INDXT) INDX(IR)=INDX(1) IR=IR-1 IF(IR .EQ. 1) THEN INDX(1)=INDXT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J .LE. IR) THEN IF(J.LT.IR) THEN IF(ARRIN(INDX(J)) .LT. ARRIN(INDX(J+1)))J=J+1 ENDIF IF(Q.LT.ARRIN(INDX(J))) THEN INDX(I)=INDX(J) I=J J=J+J ELSE J=IR+1 ENDIF GOTO 20 ENDIF INDX(I)=INDXT GOTO 10 END C C*********************************************************************** C SUBROUTINE IPOLW(XA,NA,XP,NP,APOL,KPOL,K0,K1) C C FINDS LINEAR INTERPOLATION WEIGHTS AND INDICES C C INPUT: C XA INTERPOLATION ARRAY C NA DIMENSION OF XA C XP ARRAY WITH INTERPOLATION ABSCISSAS C NP DIMENSION OF XP C OUTPUT: C APOL INTERPOLATION WEIGHTS C KPOL INTERPOLATION INDICES C K0 FIRST POINT OF XP THAT IS WITHIN XA RANGE C K1 LAST POINT OF XP THAT IS WITHIN XP RANGE C C THE INTERPOLATED VALUE YP IS TO BE CALCULATED AS: C YP(K) = YA(KPOL(K)) + APOL(K)*( YA(KPOL(K)+1)-YA(KPOL(K)) ) C C: C: IPOLW 88-05-04 NEW ROUTINE: (MATS CARLSSON) C: FINDS LINEAR INTERPOLATION WEIGHTS AND INDICES C: USED IN (ITRAD=4) OPTION C: INCLUDE 'PREC' DIMENSION XA(NA),XP(NP),APOL(NP) INTEGER KPOL(NP) C L=0 K0=1 K1=NP DO 200 K=1,NP 100 CONTINUE L=L+1 IF(XP(K).GE.XA(L) .AND. L.LT.NA) GOTO 100 C C CHECK FOR EXTRAPOLATION C IF(L.EQ.1) THEN K0=K+1 KPOL(K)=1 ELSE IF(L.EQ.NA .AND. XP(K).GT.XA(L)) THEN IF(K1.EQ.NP) K1=K-1 KPOL(K)=NA-1 ELSE KPOL(K)=L-1 ENDIF APOL(K)=( XP(K)-XA(KPOL(K)) )/( XA(KPOL(K)+1)-XA(KPOL(K)) ) L=L-1 200 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE DAMPF(KF) C C CALCULATES DAMPING PARAMETERS FOR FIXED TRANSITIONS C BASED ON ROUTINE DAMP C: C: DAMPF 88-05-04 NEW ROUTINE: (MATS CARLSSON) C: CALCULATES DAMPING PARAMETERS FOR FIXED B-B TRANSITIONS C: (ITRAD=4) C: C 00-11-27 MODIFICATIONS: (PAUL BARKLEM) C ABO THEORY X-SECTIONS IMPLEMENTED C IF GW > 20. ASSUMED X-SECTION DATA C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CFIX' C SAVE ICALL DATA ICALL/0/ C ICALL=ICALL+1 I=IFX(KF) J=JFX(KF) WAVE=HCE/(EV(J)-EV(I)) C IF(GWF(KF).NE.0.0) THEN C C FIND CONTINUUM LEVEL C DO 200 IC=J+1,NK IF(ION(IC).EQ.ION(J)+1) GOTO 300 200 CONTINUE CALL STOP(' DAMPF: NO OVERLYING CONTINUUM') 300 CONTINUE ZZ=ION(I) C625=1.283984E-12*ZZ**0.8*(1./(EV(IC)-EV(J))**2- * 1./(EV(IC)-EV(I))**2)**.4 ENDIF GR=GAF(KF) DO 400 K=1,NDEP T=TEMP(K) PE=NE(K)*BK*T TOTHI=0.0 DO 100 J=1,5 TOTHI=TOTHI+NH(J,K) 100 CONTINUE IF(ICALL.EQ.1) THEN C C CALCULATE DOPPLER WIDTH C DNYD(K)=SQRT(2.*BK*T/AWGT)*1.E-5/QNORM C C INCLUDE MICROTURBULENCE C DNYD(K)=SQRT(DNYD(K)**2+VTURB(K)**2) ENDIF C C CALCULATE GAMMA, ADAMP C C THIS SECTION BELOW ALTERED PB, TO USE X-SECTIONS NOV 2000 C IF(GWF(KF).GE.20.0) THEN GV=GWF(KF) SIGMA=INT(GV)*2.80028E-21 XALPHA=GV-INT(GV) GX=2.-XALPHA*.5 GX=GX-1.0 GAMMAF=1+(-.5748646+(.9512363+(-.6998588+ * (.4245549-.1010678*GX)*GX)*GX)*GX)*GX GV=(4./PI)**(XALPHA*0.5)*GAMMAF*1.E4*SIGMA VBAR=SQRT(21173.*T*(1./1.008+1./(AWGT/UU))) GV=GV*((VBAR/1.E4)**(1.-XALPHA)) GV=GV*TOTHI*1.E6*2. ELSEIF(GWF(KF).NE.0.0) THEN GV=GWF(KF)*8.411*(8.*BK*T/PI*(1./(1.008*UU)+1./AWGT))**0.3* * TOTHI*C625 ELSE GV=0.0 ENDIF GS=GQF(KF)*NE(K) IF(ATOMID.EQ.'H ' .AND. J.EQ.3 .AND. I.EQ.2) * GS=4.737E-7*NH(1,K) GAMMA=GR+GV+GS DOP=DNYD(K)*QNORM/WAVE*1.E13 ADAMPF(K,KF)=GAMMA/(4.*PI*DOP) 400 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE WRJFIX C C WRITES JFIX2 FILE FROM JNY DATA. SWITCHED ON BY IWJFIX.NE.0 C: C: WRJFIX 88-05-05 NEW ROUTINE: (MATS CARLSSON) C: WRITES JFIX FILE FROM JNY DATA C: USED IN (ITRAD=4) OPTION C: C: 89-03-26 MODIFICATIONS: (MATS CARLSSON) C: JFIX FILE NAME CHANGED TO JFIX2 TO ENABLE FIRST INPUT OF C: FILE JFIX AND THEN OUTPUT TO JFIX2. C: C: 95-08-25 MODIFICATIONS: (MATS CARLSSON) C: JFIX WRITING SWITCHED ON BY IWJFIX.NE.0 INSTEAD OF IWRAD.LT.0 C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLGMX' INCLUDE 'CLU' C DIMENSION QJ(2*MQ*MRAD) INTEGER IRC(2*MQ*MRAD) C IF(IWJFIX.EQ.0) RETURN C C OPEN DIRECT ACCESS FILE JFIX2 FOR WRITING C CALL OPEN(LJFIX,'JFIX2',NDEP,'NEW') C C FIND FREQUENCIES FOR ALL JNY RECORDS C FOR ONE-SIDED LINES, WRITE WHOLE PROFILE C STORE JNY FILE RECORD NUMBER IN IRC C ICUM IS LAST RECORD NUMBER FOR PREVIOUS TRANSITION C IQ=0 ICUM=0 DO 200 KR=1,NRAD Q00=CC/ALAMB(KR)*1.E8 IF(IND(KR).EQ.1 .AND. KR.LE.NLINE) THEN DO 50 NY=NQ(KR),2,-1 IQ=IQ+1 QJ(IQ)=Q00*(1.-Q(NY,KR)*QNORM*1.E5/CC) IRC(IQ)=ICUM+NY 50 CONTINUE ENDIF DO 100 NY=1,NQ(KR) IQ=IQ+1 QJ(IQ)=Q00*(1.+Q(NY,KR)*QNORM*1.E5/CC) IRC(IQ)=ICUM+NY 100 CONTINUE ICUM=ICUM+NQ(KR) 200 CONTINUE NQJ=IQ C C WRITE HEADER RECORDS TO FILE JFIX2 C IREC=1 WRITE(LJFIX,REC=IREC) NQJ,NDEP DO 300 IQ0=1,NQJ,NDEP IREC=IREC+1 WRITE(LJFIX,REC=IREC) * (QJ(IQ),IQ=IQ0,MIN(IQ0+NDEP-1,NQJ)) 300 CONTINUE IREC=IREC+1 WRITE(LJFIX,REC=IREC) (LOG10(TAU(K)),K=1,NDEP) IREC=IREC+1 WRITE(LJFIX,REC=IREC) (LOG10(CMASS(K)),K=1,NDEP) C C WRITE JNY/BNU TO FILE C DO 400 IQ=1,NQJ IREC=IREC+1 CALL READJ(IRC(IQ)) WRITE(LJFIX,REC=IREC) (LOG(JNY(K)/PLANCK(QJ(IQ),TEMP(K))), * K=1,NDEP) 400 CONTINUE CALL CLOSE(LJFIX) C RETURN END C C*********************************************************************** C SUBROUTINE LJUST(TEXT) C C LEFT JUSTIFIES THE STRING TEXT C INCLUDE 'PREC' CHARACTER*(*) TEXT C L=LEN(TEXT) DO 100 J=1,L IF(TEXT(J:J).NE.' ') GOTO 200 100 CONTINUE 200 CONTINUE DO 300 I=1,L IF(J.LE.L) THEN TEXT(I:I)=TEXT(J:J) ELSE TEXT(I:I)=' ' ENDIF J=J+1 300 CONTINUE C RETURN END C C************************************************************************ C SUBROUTINE LTEEQW C C COMPUTES LTE EQUIVALENT WIDTHS C C: C: LTEEQW 88-02-02 MODIFICATIONS: (MATS CARLSSON) C: MEAN INTENSITY IS LAMBDA ITERATED TO CORRECTIONS SMALLER C: THAN ELIM2 C: C: 88-06-22 MODIFICATIONS: (MATS CARLSSON) C: ICONV IS SET TO 0 DURING SCATTERING LAMBDA ITERATIONS TO C: AVOID PRINTOUT IN WTEST CALLED FROM TRPT C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: SCATTERING VERSION OF TRANSP USED INSTEAD OF C: SCATTERING LAMBDA ITERATIONS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION OLDN(MK,MDEP) DATA LMAX/100/ C IF(IWEQW.EQ.0) RETURN C CALL CPUTIME('LTEEQW',0,0,2) DO 200 K=1,NDEP DO 100 I=1,NK OLDN(I,K)=N(I,K) N(I,K)=NSTAR(I,K) 100 CONTINUE 200 CONTINUE C C GET CONSISTENT JNY IN SCATTERING PART OF C THE SOURCE FUNCTION C ISCAT0=ISCAT ICONV0=ICONV ISCAT=1 ICONV=0 CALL TRPT CALL TRCONT ISCAT=ISCAT0 ICONV=ICONV0 DO 300 KR=1,NLINE CALL EQWDTH(WW,KR) WEQLTE(KR)=WW 300 CONTINUE C C RESET JNY TO ZERO C RESTORE POPULATIONS C DO 350 K=1,NDEP DO 330 I=1,NK N(I,K)=OLDN(I,K) 330 CONTINUE JNY(K)=0.0 350 CONTINUE IREC=0 DO 500 KR=1,NRAD DO 400 NY=1,NQ(KR) IREC=IREC+1 CALL WRITEJ(IREC) 400 CONTINUE 500 CONTINUE C CALL CPUTIME('LTEEQW',0,2,2) RETURN END C C*********************************************************************** C SUBROUTINE LTEPOP C C CALCULATES LTE POPULATIONS C: C: LTEPOP 87-04-07 MODIFICATIONS: (PHILIP JUDGE) C: A DANGER SIGN IS OUTPUT WHEN ZERO POPULATIONS ARE FOUND. C: C: 88-02-01 MODIFICATIONS: (MATS CARLSSON) C: PROVISION FOR MOLECULES ADDED C: C: 88-07-01 MODIFICATIONS: (MATS CARLSSON) C: DEBYE LOWERING OF IONIZATION POTENTIAL TAKEN INTO ACCOUNT C: C: 92-09-10 MODIFICATIONS: (MATS CARLSSON) C: DANGER SIGN LIMIT CHANGED FROM 1.E-37 TO 0.0 C: C: 98-09-04 MODIFICATIONS: (MATS CARLSSON) C: QCALC REPLACED BY QPART FOR CO, CH, CN C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CCONST' C C PGJ ADDED FOLLOWING LINE INCLUDE 'CLU' C CHARACTER*4 ELEMID C DIMENSION TNS(MK),NDXI(20) LOGICAL MOLEC C* C* 88-07-01 MATS CARLSSON C* THE ENERGIES OF IONIZATION ARE REDUCED BY ION*DXI, FOLLOWING BASCHEK ET C* AL., ABH. HAMB. VIII, 26 EQ. (10). C* SINCE ALL EV REFERS TO THE LOWEST LEVEL, IT IS NECESSARY TO SUBTRACT C* FROM THE ENERGY OF A LEVEL ALL DXI EXPERIENCED BY LOWER IONIZATION C* STAGES. THIS NUMBER IS STORED IN NDXI(ION) C* NDXI(1)=0 DO 50 I=2,20 NDXI(I)=NDXI(I-1)+(I-1) 50 CONTINUE DO 60 I=1,20 NDXI(I)=NDXI(I)-NDXI(ION(1)) 60 CONTINUE C CALL GETWRD(ATOMID,1,K1,K2) ELEMID=ATOMID(K1:K2) MOLEC=ELEMID.EQ.'CH' .OR. ELEMID.EQ.'CO' .OR. ELEMID.EQ.'CN' CCON=0.5*(HH/SQRT(2.*PI*EM)/SQRT(BK))**3 DO 500 K=1,NDEP T=TEMP(K) CONL=LOG(CCON*NE(K))-1.5*LOG(T) DXI=4.98E-4*5040./T*SQRT(NE(K)*BK*T) SUMN=1. DO 300 I=2,NK EVI=EV(I)-DXI*NDXI(ION(I)) TNSL=LOG(G(I))-LOG(G(1))-EK/T*EVI IF(ION(I).LE.ION(1)) GOTO 100 L=ION(I)-ION(1) TNSL=TNSL-FLOAT(L)*CONL 100 TNS(I)=EXP(TNSL) SUMN=SUMN+TNS(I) 300 CONTINUE IF(.NOT.MOLEC) THEN NSTAR(1,K)=TOTN(K)/SUMN ELSE NSTAR(1,K)=TOTN(K)*G(1)*EXP(-EK/T*EV(1))/QPART(ELEMID,K) ENDIF DO 400 I=2,NK NSTAR(I,K)=TNS(I)*NSTAR(1,K) C C PGJ ADDITION: IF(NSTAR(I,K) .LE. 0.0)WRITE(LJOBLO,1000)I,K 1000 FORMAT(' LTEPOP: LEVEL ',I3,' AT DEPTH ',I3,' HAS ZERO *LTE POPULATION') C END PGJ ADDITION C 400 CONTINUE C C SCALE TOTAL ABUNDANCES IF TREATED MODEL IS A MOLECULE C IF(MOLEC) THEN TOTN(K)=SUMN*NSTAR(1,K) ENDIF 500 CONTINUE RETURN END C C*********************************************************************** C FUNCTION QPART(ELEMID,K) C C CALCULATES PARTITION FUNCTIONS FOR MOLECULESN CO, CN AND CH C DATA FOR CO AND CN FROM C SAUVAL & TATUM APJ SUPP SERIES 56:193,1984, TABLE 5 C FOR PARTITION FUNCTIONS OF CH, FUNCTION QJORG IS CALLED C C: C: QPART 98-08-13 NEW ROUTINE: (ROUPPE VAN DER VOORT) C: CALCULATES PARTITION FUNCTIONS FOR CO, CH, CN C: REPLACES OLD QCALC C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CCONST' C CHARACTER*(*) ELEMID C IF(ELEMID.EQ.'CH') THEN QPART = QJORG(K) ELSE IF(ELEMID.EQ.'CO' .OR. ELEMID .EQ. 'CN') THEN IF(ELEMID.EQ.'CN') THEN A0 = 4.0078 A1 = -2.1514 A2 = 0.9226 A3 = -0.1671 ELSE IF(ELEMID .EQ. 'CO') THEN A0 = 3.6076 A1 = -1.7608 A2 = 0.4172 A3 = 0.0 ENDIF THETA = 5040./TEMP(K) C C USE A POLYNOMIAL EXPANSION OF LOG10(Q) QPART = A0 + A1*(LOG10(THETA)) + A2*(LOG10(THETA))**2 + * A3*(LOG10(THETA))**3 C CONVERT LOG10(Q) -> Q QPART = 10**QPART ELSE CALL STOP('QPART: ONLY HANDLES CO, CH AND CN') ENDIF C RETURN END C C*********************************************************************** C FUNCTION QJORG(K) C C CALCULATES PARTITION FUNCTION FOR CH C DATA FROM U.G. JORGENSEN ET AL. A&A 1996, 315, 204 C C: QJORG 98-08-13 NEW ROUTINE: (ROUPPE VAN DER VOORT) C: CALCULATES PARTITION FUNCTIONS FOR CH C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CCONST' C C DATA FROM JORGENSEN, TABLE 1. C TEMPERATURE RANGE T={1000, 8000} C NOTE THAT FUNCTION EXTAPOLATES FOR TEMP>8000 AND TEMP<1000 C FOR EXTRAPOLATION FOR TEMP>8500 OR TEMP<1000 THE RESULTS C DIFFER SIGNIFICANT FROM A POLYNOMIAL EXTRAPOLATION (POLINT) C DIMENSION PFDATA(15), TDATA(15) DATA PFDATA/203.5, 324.7, 471.2, 646.2, 850.8, 1085.1, 1348.1, * 1638.7, 1955.3, 2296.1, 2659.3, 3043.2, 3446.2, 3866.8, 4303.7/, * TDATA/1000., 1500., 2000., 2500., 3000., 3500., 4000., 4500., * 5000., 5500., 6000., 6500., 7000., 7500., 8000./ C SIMPLE LINEAR INTERPOLATION BETWEEN TWO NEAREST POINTS DO 10 I = 1, 15 IF (TEMP(K) .EQ. TDATA(I)) THEN QJORG = PFDATA(I) RETURN ELSEIF (TEMP(K) .GT. TDATA(I) .AND. TEMP(K) .LT. TDATA(I+1)) * THEN X0 = TDATA(I) X1 = TDATA(I+1) Y0 = PFDATA(I) Y1 = PFDATA(I+1) ELSEIF (TEMP(K) .LT. TDATA(1)) THEN X0 = TDATA(1) X1 = TDATA(2) Y0 = PFDATA(1) Y1 = PFDATA(2) ELSEIF (TEMP(K) .GT. TDATA(15)) THEN X0 = TDATA(14) X1 = TDATA(15) Y0 = PFDATA(14) Y1 = PFDATA(15) ENDIF 10 CONTINUE B = Y0 A = (Y1 - Y0)/(X1 - X0) QJORG = A*(TEMP(K) - X0) + B END C C*********************************************************************** C SUBROUTINE MXDLG(INIT) C C EVALUATES MAX DELTA LG TAUNY, MAX DELTA LG S AND MAX STEP IN P-S. C THIS IS DONE ONLY IN THE TAUNY RANGE TMIN-TMAX. TMIN AND TMAX C GET THEIR VALUES IN A DATA STATEMENT. C MAX DELTA LG TAUNY IS USED TO INTERPOLATE A NEW DEPTH SCALE IN DSCAL2 C: C: MXDLG 88-01-26 MODIFICATIONS: (MATS CARLSSON) C: INSTEAD OF MAX DELTA LG S AND MAX STEP IN P-S C: THE MAXIMUM AND MINIMUM LG TAUNY(MU=NMU) AT EACH DEPTHPOINT IS C: STORED C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CLGMX' INCLUDE 'CINPUT' C DATA TMIN,TMAX/1.E-3,1.E+2/ C IF(INIT.EQ.0) THEN DO 50 K=1,NDEP DLGTMX(K)=0.0 TAUMIN(K)=1.E37 TAUMAX(K)=1.E-37 50 CONTINUE ELSE KMIN=2 C C STEP SIZE IS ONLY CHECKED FOR TAUQ BETWEEN TMIN AND TMAX C DO 70 K=2,NDEP-2 IF(TAUQ(K).LT.TMIN) KMIN=K IF(TAUQ(K).GT.TMAX) GOTO 80 70 CONTINUE 80 KMAX=K+1 DO 100 K=KMIN,KMAX IF(TAUQ(K).GT.0.0 .AND. TAUQ(K-1).GT.0.0) THEN DTNY=LOG10(TAUQ(K))-LOG10(TAUQ(K-1)) ELSE DTNY=97. ENDIF DLGTMX(K)=MAX(DLGTMX(K),DTNY) 100 CONTINUE IF(XMY.EQ.XMU(NMU)) THEN DO 200 K=1,NDEP TAUMIN(K)=MIN(TAUMIN(K),TAUQ(K)) TAUMAX(K)=MAX(TAUMAX(K),TAUQ(K)) 200 CONTINUE ENDIF ENDIF RETURN END C C*********************************************************************** C C C*************************************************************** C FUNCTION PESC(TAU0,DAMP,INDEX) C C ESCAPE PROBABILITITES FOR VOIGT LINES AND CONTINUA C WHEN INDEX .GT. 0 THE CONTINUUM ESCAPE PROBABILITY IS USED C OTHERWISE LINES ARE ASSUMED C C THIS VERSION HAS CORRECTION (FUDGE) FOR NEGATIVE OPACITIES C INCLUDE 'PREC' INCLUDE 'CLU' C C CONTINUA- ON THE SPOT APPROXIMATION C IF (INDEX.GT.0) THEN PESC = 1.0 IF (TAU0.GE.1.0) PESC = 0.0 ELSE C C LINES- SEVERAL APPROXIMATIONS REPLACED BY APPROXIMATE PESC= 1/TAU C C OPTICALLY THIN LINES: PESC=1.0: C IF (ABS(TAU0).LT.1.01) THEN PESC=1.0 RETURN ENDIF TAUZ = TAU0 C C NEGATIVE OPACITIES, SET TAUZ=-TAU0 AND PESC=-1./TAUZ : C IF (TAU0.LT.0) THEN WRITE (LJOBLO,FMT=1000) 1000 FORMAT (' PESC: NEGATIVE OPACITY CHANGED TO POSITIVE') TAUZ = -TAU0 END IF PESC = 1.0/ (TAUZ) ENDIF C RETURN END FUNCTION PGJRT(FNU,ALAMB) C C A FUNCTION TO RETURN A BLACK-BODY RADIATION TEMPERATURE C FROM A GIVEN FLUX C C PHILIP JUDGE, OXFORD, OCTOBER, 1986 C C INPUT: C FNU - FLUX ( ERG/CM2/S/HZ ) C ALAMB - WAVELENGTH (ANGSTROM) C C OUTPUT: C PGJRT - RADIATION TEMPERATURE IN K C C: PGJRT 93-08-25 MODIFICATIONS: (MATS CARLSSON) C; EXPRESSION CHANGED FOR SMALL LAMBDA TO AVOID OVERFLOW C: INCLUDE 'PREC' INCLUDE 'CCONST' C A3FNU=ALAMB*ALAMB*ALAMB*FNU IF(A3FNU.GT.1.0) THEN PGJRT = (HH*CC*1.E8/BK) / ( ALAMB * * LOG(PI*HH*CC*2.E24/A3FNU + 1.0) ) ELSE IF(A3FNU.GT.0.0) THEN PGJRT = (HH*CC*1.E8/BK) / ( ALAMB * * (LOG(PI*HH*CC*2.E24)-LOG(A3FNU)) ) ELSE PGJRT = 0.0 ENDIF RETURN END C C ***************************************************************** C FUNCTION PHII(X) C C H1 COLLISIONAL ROUTINE C SOURCE: LINEAR-B C INCLUDE 'PREC' INCLUDE 'CLU' C DIMENSION XY(105), XQ(105) DATA XY/ .005, .01, .015, .02, .03, .04, .05, .06, .07, .1, .15, * .2, .25, .3, .4, .5, .6, .7, .8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.2, * 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8, 4., 4.2, 4.4, 4.6, 4.8, * 5., 5.2, 5.4, 5.6, 5.8, 6., 6.4, 6.8, 7.2, 7.6, 8., 8.4, 8.8, * 9.2, 9.6, 10., 10.4, 10.8, 11.2, 11.6, 12., 12.4, 12.8, 13.2, * 13.6, 14., 14.4, 14.8, 15.2, 15.6, 16., 16.4, 16.8, 17.2, 17.6, * 18., 19., 20., 21., 22., 23., 24., 25., 26., 28., 30., 32., * 34., 36., 38., 40., 44., 48., 52., 56., 60., 64., 68., 72., * 76., 80., 84., 88., 92., 96., 100./ DATA XQ/ 4.68662, 3.98908, 3.58170, 3.29421, 2.89337, 2.61378, * 2.40078, 2.22986, 2.08793, 1.77036, 1.43206, 1.21000, 1.04966, * 9.27071E-1, 7.50166E-1, 6.27566E-1, 5.37143E-1, 4.67571E-1, * 4.12367E-1, 3.30382E-1, 2.72591E-1, 2.29842E-1, 1.97082E-1, * 1.71284E-1, 1.50521E-1, 1.33120E-1, 1.19370E-1, 1.07461E-1, * 9.73233E-2, 8.86113E-2, 8.10617E-2, 7.44706E-2, 6.86781E-2, * 6.35567E-2, 5.90043E-2, 5.49376E-2, 5.12883E-2, 4.80000E-2, * 4.50257E-2, 4.23258E-2, 3.98670E-2, 3.76207E-2, 3.55629E-2, * 3.36727E-2, 3.19321E-2, 2.88391E-2, 2.61813E-2, 2.38798E-2, * 2.18729E-2, 2.01118E-2, 1.85576E-2, 1.71788E-2, 1.59497E-2, * 1.48492E-2, 1.38599E-2, 1.29670E-2, 1.21585E-2, 1.14239E-2, * 1.07544E-2, 1.01425E-2, 9.58173E-3, 9.06651E-3, 8.59202E-3, * 8.15405E-3, 7.74893E-3, 7.37343E-3, 7.02472E-3, 6.70030E-3, * 6.39795E-3, 6.11572E-3, 5.85184E-3, 5.60475E-3, 5.37305E-3, * 5.15549E-3, 4.95093E-3, 4.49000E-3, 4.09070E-3, 3.74260E-3, * 3.43720E-3, 3.16780E-3, 2.92900E-3, 2.71630E-3, 2.52590E-3, * 2.20080E-3, 1.93470E-3, 1.71420E-3, 1.52940E-3, 1.37300E-3, * 1.23940E-3, 1.12450E-3, 9.37820E-4, 7.94100E-4, 6.81070E-4, * 5.9057E-4,5.1699E-4,4.5636E-4,4.0580E-4,3.6321E-4,3.2699E-4, * 2.9593E-4,2.6909E-4,2.4575E-4,2.2532E-4,2.0733E-4,1.9142E-4 / C DO 100 I=2,105 IF(X.LT.XY(I)) GO TO 300 100 CONTINUE PHII=0. WRITE(LJOBLO,200) X 200 FORMAT(' PHI WAS SET TO ZERO BECAUSE X = ',E16.6) RETURN C 300 CONTINUE XZ= XQ(I-1)+(X-XY(I-1))/(XY(I)-XY(I-1))*(XQ(I)-XQ(I-1)) PHII= XZ/EXP(X) C RETURN END C C********************************************************************** C FUNCTION PLANCK(U,T) C C CALCULATES PLANCK FUNCTION BNY AT FREQUENCY U, TEMP T C INCLUDE 'PREC' INCLUDE 'CCONST' C X=HH*U/BK/T IF(X.LT.80.) THEN PLANCK=2.0*HH*U/CC*U/CC*U/(EXP(X)-1.0) ELSE PLANCK=2.0*HH*U/CC*U/CC*U*EXP(-X) ENDIF C RETURN END C C********************************************************************** C SUBROUTINE PROFIL C C CALCULATES VOIGT PROFILE AND WRITES PROFILE TO FILE PHI C C: PROFIL 90-12-28 MODIFICATIONS: (MATS CARLSSON) C: INITIALIZING PHI(NDEP+1:MDEP) C: TO ZERO TO ENABLE I/O OF WHOLE ARRAY C: C: 93-01-18 MODIFICATIONS: (MATS CARLSSON) C: ADDED OPTION OF HAVING SEVERAL LINES BETWEEN TERMS C: FOR DOCUMENTATION OF VARIABLES, SEE ROUTINE ATOM C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CLU' INCLUDE 'CTERM' C CALL CPUTIME('VOIGT ',0,0,1) CALL REWIND(LPHI) C C INITIALIZE PHI(NDEP+1:MDEP) TO ZERO C DO 5 K=NDEP+1,MDEP PHI(K)=0.0 5 CONTINUE C DO 200 KR=1,NLINE DO 100 K=1,NDEP WPHI(K,KR)=0. 100 CONTINUE 200 CONTINUE DO 700 KR=1,NLINE DO 500 NY=1,NQ(KR) DO 400 MU=1,NMU WQMU=WQ(NY,KR)*WMU(MU) DO 300 K=1,NDEP IF(IND(KR).EQ.2 .OR. MU.EQ.1) THEN IF(NTERM(KR).LE.1) THEN V=(Q(NY,KR)-XMU(MU)*VEL(K))/DNYD(K) CALL VOIGT(ADAMP(K,KR),V,H) PHI(K)=H/(DNYD(K)*1.772453851) ELSE PHI(K)=0.0 DO 250 ITRM=1,NTERM(KR) KTRM=KTERM(ITRM,KR) DQTERM=DETERM(KTRM)*ALAMB(KR)*1.E-13*CC/QNORM V=(Q(NY,KR)-DQTERM-XMU(MU)*VEL(K))/DNYD(K) CALL VOIGT(ADTERM(K,KTRM),V,H) PHI(K)=PHI(K)+WTERM(KTRM)*H 250 CONTINUE PHI(K)=PHI(K)/(DNYD(K)*1.772453851) ENDIF ENDIF WPHI(K,KR)=WPHI(K,KR)+WQMU*PHI(K) 300 CONTINUE IF(IND(KR).EQ.2 .OR. MU.EQ.1) CALL WRITEP 400 CONTINUE 500 CONTINUE DO 600 K=1,NDEP WPHI(K,KR)=1.0/WPHI(K,KR) 600 CONTINUE 700 CONTINUE CALL REWIND(LPHI) CALL CPUTIME('VOIGT ',0,1,1) RETURN END C C***************************************************************** C FUNCTION PSI(I,X,EN) C C H1 COLLISIONAL ROUTINE C SOURCE: LINEAR-B C INCLUDE 'PREC' DIMENSION H(8), A(8), R(8), XI(8) DATA H/ 1.48, 3.64, 5.93, 8.33, 10.75, 12.9, 15.05, 17.2 /, * A/ 1.3, .59, .38, .286, .229, .192, .164, .141 /, * R/ 1.83, 1.6, 1.53, 1.495, 1.475, 1.46, 1.45, 1.45 /, * XI/ 109679., 27420., 12187., 6855., 4387., 3047., 2239., 1714. / C XG= EN/XI(I) PSI= PHII(X) + 0.19*H(I)*XG**2.5*X*EXPINT(3,X,DUM)+H(I)*(XG**R(I) * +A(I)*(XG-1.))*EXPINT(2,X,DUM) RETURN END C C********************************************************************** C SUBROUTINE REDIST(NK1,G1,LABEL1) C C REDISTRIBUTES THE POPULATIONS CALCULATED USING ATOMIC MODEL ATOM C TO THE LEVELS IN ATOMIC MODEL ATOM2 ACCORDING TO THE STATISTICAL C WEIGHTS. THE LEVELS IN ATOM2 ARE CONSIDERED TO BE PART OF LEVELS C IN ATOM WITH THE SAME LABEL. ALL LABELS IN ATOM2 MUST EXIST IN C ATOM. C C: REDIST 91-06-28 MODIFICATIONS: (MATS CARLSSON) C: LENGTH OF LABEL1 CHANGED TO 20 CHARACTERS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CLU' C DIMENSION G1(MK),AN1(MK,MDEP),ANSTR1(MK,MDEP) CHARACTER*20 LABEL1(MK) C C IF THERE ARE SEVERAL LEVELS IN ATOM WITH THE SAME LABEL: C ADD UP THEIR POPULATIONS AND STATISTICAL WEIGHTS C INTO THE LOWEST LEVEL C DO 90 I=1,NK1-1 DO 80 J=I+1,NK1 IF(LABEL1(J).EQ.LABEL1(I)) THEN WRITE(LJOBLO,50) I,J 50 FORMAT(' REDIST: WARNING: LEVELS',2I3,' ADDED BEFORE', * ' SPLIT UP') G1(I)=G1(I)+G1(J) DO 70 K=1,NDEP N(I,K)=N(I,K)+N(J,K) NSTAR(I,K)=NSTAR(I,K)+NSTAR(J,K) 70 CONTINUE ENDIF 80 CONTINUE 90 CONTINUE C C STORE OLD POPULATIONS IN AN1 AND ANSTR1 C DO 200 K=1,NDEP DO 100 I=1,NK1 AN1(I,K)=N(I,K) ANSTR1(I,K)=NSTAR(I,K) 100 CONTINUE 200 CONTINUE C C REDISTRIBUTE C DO 600 I=1,NK DO 300 J=1,NK1 IF(LABEL1(J).EQ.LABEL(I)) GOTO 400 300 CONTINUE CALL STOP('REDIST: LABEL IN ATOM2 MISSING IN ATOM') 400 CONTINUE DO 500 K=1,NDEP N(I,K)=G(I)/G1(J)*AN1(J,K) NSTAR(I,K)=G(I)/G1(J)*ANSTR1(J,K) 500 CONTINUE 600 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE RINPUT C C READS AN INPUT FILE CONSISTING OF VARIABLE NAMES AND VALUES C THE VARIABLE NAMES AND VALUES MUST OCCUR IN PAIRS: C NAME, VALUE, NAME, VALUE, ETC. C IERR IS INITIALIZED TO 0. C IF A VARIABLE IN THE LIST IS NOT GIVE A VALUE IN THE INPUT FILE, C A WARNING IS ISSUED AND IERR IS SET TO 2. C IF A VARIABLE IN THE INPUT FILE IS NOT IN THE LIST C A WARNING IS ISSUED AND IERR IS SET TO 1. C IF IERR IS 2 AFTER THE READ OPERATION, STOP IS CALLED C C LIST OF LOCAL VARIABLES: C MVAR NUMBER OF VARIABLES IN THE VARIABLE LIST C VNAME ARRAY CONTAINING THE VARIABLE NAMES C CVALUE CHARACTER ARRAY WHERE THE VARIABLE VALUES ARE PUT AT INPUT C INDV = 1 FOR VARIABLE NAMES C =-1 FOR VARIABLE VALUES C C IF THE INPUT LIST IS CHANGED, CHANGES HAVE TO BE MADE THREE PLACES: C DATA VNAME/.../ C READ(LDUMI,*) C CVALUE(I)= FOR DEFAULT VALUES C C THE ORDER OF VARIABLES MUST BE THE SAME IN DATA VNAME AND IN C READ(LDUMI,*). C C IN ADDITION THE OUTPUT LIST SHOULD BE CHANGED IN WINPUT C AND RELEVANT COMMON BLOCKS SHOULD BE CHANGED C C: RINPUT 90-10-17 MODIFICATIONS: (MATS CARLSSON) C: ICONV.GE.10 SIGNALS THAT NG ACCELERATION IS NOT TO BE C: PERFORMED. NGACC IS THEN SET .FALSE. AND 10 IS SUBTRACTED C: FROM ICONV. NGACC IS USED IN ITER AND TRANSFERRED IN C: COMMON BLOCK CNGACC C: C: 95-02-27 MODIFICATIONS: (MATS CARLSSON) C: NEW VARIABLES WITH DEFAULT VALUES IF NOT GIVEN IN INPUT FILE: C: ADD DEFAULT VALUE TO CVALUE C: C: 95-08-04 MODIFICATIONS: (MATS CARLSSON) C: PARAMETER STATEMENT NO LONGER DISTINGUESHES BETWEEN C: REAL AND INTEGER TYPE VARIABLES, NEW MVAR IS THE SUM OF C: THE OLD PARAMETERS MAVAR AND MIVAR. MVAR IS ONLY UPPER C: LIMIT, THE ACTUAL NUMBER OF INPUT VARIABLES IS FOUND BY C: NAMING THE LAST+1 'END ' C: NEW VARIABLES INTRODUCED C: IDL1 WRITE IDL1 FILE C: IDLNY WRITE IDLNY FILE C: IDLCNT WRITE IDLCNT FILE C: INGACC SWITCHES ON NG ACCELERATION C: IOPACL OPACITIES FROM LINES C: ISCAT SCATTERING VERSION OF TRANSP C: INCRAD INCOMING RADIATION FIELD READ C: ICRSW REGULATES COLLISIONAL RADIATIVE SWITCHING C: IWARN REGULATES CLASSES OF WARNING MESSAGES TO WRITE TO JOBLOG C: IWJFIX WRITE JFIX2 FILE C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CLU' COMMON /CNGACC/ NGACC LOGICAL NGACC C PARAMETER (MVAR=60) CHARACTER*20 CVALUE(MVAR) CHARACTER*80 TEXT CHARACTER VNAME(MVAR)*6 DATA (VNAME(I),I=1,48) * /'DIFF ','ELIM1 ','ELIM2 ','QNORM ','THIN ', * 'IATOM2','ICONV ','IHSE ','ILAMBD','IOPAC ','IOPACL', * 'ISTART','ISUM ','ITMAX ','ITRAN ','ISCAT ','INCRAD', * 'INGACC','ICRSW ','NMU ', * 'IWABND','IWATMS','IWATOM','IWCHAN','IWDAMP','IWEMAX', * 'IWEQW ','IWEVEC','IWHEAD','IWHSE ','IWLGMX','IWLINE', * 'IWLTE ','IWN ','IWNIIT','IWOPAC','IWRAD ','IWRATE', * 'IWSTRT','IWTAUQ','IWTEST','IWWMAT','IWJFIX','IWARN ', * 'IDL1 ','IDLNY ','IDLCNT','END '/ C C FIND NUMBER OF INPUT VARIABLES BY FINDING STRING 'END ' C DO 10 I=1,MVAR IF(VNAME(I).EQ.'END ') GOTO 20 10 CONTINUE CALL STOP('RINPUT: END SHOULD BE LAST ID IN ARRAY VNAME') 20 CONTINUE NVAR=I-1 C IERR=0 INDV=-1 DO 50 I=1,NVAR CVALUE(I)=' ' 50 CONTINUE C C DEFAULT VALUES C CALL DEFVAL(VNAME,CVALUE,MVAR,'IWARN ','2') CALL DEFVAL(VNAME,CVALUE,MVAR,'IWJFIX','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'IOPACL','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'ISCAT ','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'INCRAD','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'INGACC','1') CALL DEFVAL(VNAME,CVALUE,MVAR,'ICRSW ','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'IDL1 ','1') CALL DEFVAL(VNAME,CVALUE,MVAR,'IDLNY ','0') CALL DEFVAL(VNAME,CVALUE,MVAR,'IDLCNT','0') C C GO THROUGH THE INPUT FILE LINE BY LINE C CALL CSTRIP(LINPUT,'INPUT') 100 CONTINUE K0=1 READ(LDUMS,200,END=700) TEXT 200 FORMAT(A) C C SEARCH LINE WORD BY WORD. EVERY SECOND WORD IS CONSIDERED TO BE C A VARIABLE NAME, EVERY SECOND ITS VALUE. C = IS IGNORED. K1=0 IN GETWRD IS A FLAG THAT THERE ARE NO MORE WORDS C ON THE LINE C 300 CONTINUE CALL GETWRD(TEXT,K0,K1,K2) IF(K1.EQ.0) GOTO 600 K0=K2+2 IF(TEXT(K1:K2).EQ.'=') GOTO 300 INDV=-INDV IF(INDV.EQ.1) THEN C C INDV=1 : VARIABLE NAME C DO 400 I=1,NVAR IF(TEXT(K1:K2).EQ.VNAME(I)(1:K2-K1+1)) GOTO 500 400 CONTINUE C C INPUT VARIABLE IS NOT IN LIST, OUTPUT WARNING AND SKIP THE C FOLLOWING STRING (SUPPOSED TO BE THE VALUE OF THE VARIABLE) C IERR=1 WRITE(LJOBLO,420) TEXT(K1:K2) 420 FORMAT(4X,A,' IS NOT IN LIST') CALL GETWRD(TEXT,K0,K1,K2) IF(K1.EQ.0) THEN K0=1 450 READ(LDUMS,200,END=700) TEXT CALL GETWRD(TEXT,K0,K1,K2) IF(K1.EQ.0) GOTO 450 ENDIF K0=K2+2 INDV=-INDV GOTO 300 C C INPUT VARIABLE IS IN LIST C 500 CONTINUE IVAR=I ELSE C C INDV = -1 : VARIABLE VALUE C CVALUE(IVAR)=' ' CVALUE(IVAR)(20-K2+K1:20)=TEXT(K1:K2) ENDIF GOTO 300 600 CONTINUE GOTO 100 700 CONTINUE C C ASSOCIATE VALUE TO VARIABLE C FIRST WRITE VALUES TO FILE, THEN READ VALUES IN FREE FORMAT C OUTPUT WARNING IF NO VALUE EXISTS C CALL OPEN(LDUMI,'DUMI',1,'UNKNOWN') DO 740 I=1,NVAR IF(CVALUE(I).EQ.' ') THEN IERR=2 WRITE(LJOBLO,710) VNAME(I) 710 FORMAT(4X,A,' HAS NOT A VALUE') ENDIF WRITE(LDUMI,720) CVALUE(I) 720 FORMAT(1X,A) 740 CONTINUE IF(IERR.LE.1) THEN CALL REWIND(LDUMI) READ(LDUMI,*) DIFF,ELIM1,ELIM2,QNORM,THIN, * IATOM2,ICONV,IHSE,ILAMBD,IOPAC,IOPACL, * ISTART,ISUM,ITMAX,ITRAN,ISCAT,INCRAD,INGACC,ICRSW,NMU, * IWABND,IWATMS,IWATOM,IWCHAN,IWDAMP,IWEMAX, * IWEQW,IWEVEC,IWHEAD,IWHSE,IWLGMX,IWLINE, * IWLTE,IWN,IWNIIT,IWOPAC,IWRAD,IWRATE, * IWSTRT,IWTAUQ,IWTEST,IWWMAT,IWJFIX,IWARN, * IDL1,IDLNY,IDLCNT CALL CLOSE(LDUMI) ELSE CALL STOP('RINPUT: IERR=2') ENDIF IF(ICONV.GE.10) THEN ICONV=ICONV-10 NGACC=.FALSE. ENDIF IF(INGACC.NE.0) NGACC=.TRUE. C RETURN END C C********************************************************************** C SUBROUTINE DEFVAL(VNAME,CVALUE,MVAR,CNAME,CVAL) C C SETS DEFAULT VALUE FOR INPUT VARIABLE C TYPICAL CALL C CALL DEFVAL(VNAME,CVALUE,MVAR,'ISCAT ','0') C INCLUDE 'CLU' C CHARACTER*(*) VNAME(MVAR),CVALUE(MVAR),CNAME,CVAL C DO 100 I=1,MVAR IF(CNAME.EQ.VNAME(I)) GOTO 200 100 CONTINUE WRITE(LJOBLO,150) CNAME 150 FORMAT(' CNAME=',A) CALL STOP('DEFVAL: CNAME NOT FOUND') C 200 CONTINUE C C STORE RIGHT JUSTIFIED C L2=LEN(CVALUE(I)) L1=LEN(CVAL) CVALUE(I)(L2-L1+1:L2)=CVAL C END FUNCTION RINT1(TE,TR,HKTT) C C CALCULATES STIMULATED EMISSION INTEGRAL FOR FIXED BOUND-FREE RATES C INCLUDE 'PREC' SUM=EXINT1((TE+TR)*HKTT) IF(SUM.EQ.0.0) GOTO 102 DO 101 I=2,999 XX=(FLOAT(I)*TE+TR)*HKTT DSUM=EXINT1(XX) SUM=SUM+DSUM IF(DSUM/SUM.LE.1.0E-5) GO TO 102 101 CONTINUE 102 RINT1=SUM RETURN END C C*************************************************************** C FUNCTION RINT(X) C C EVALUATES THE INTEGRAL NEEDED IN CALCULATING THE FIXED RAD. RATES C BY FINDING THE SUM OF THE EXPONENTIAL INTEGRALS. C IT LOOPS THROUGH THE SUMMATION UNTIL THE ADDED TERM IS LT C 1.0E-5 TIMES THE SUM. C INCLUDE 'PREC' SUM=EXINT1(X) IF(SUM.EQ.0.0) GOTO 999 DO 100 I=2,999 XX=FLOAT(I)*X DSUM=EXINT1(XX) SUM=SUM+DSUM IF(DSUM/SUM.LE.1.E-5) GOTO 999 100 CONTINUE 999 RINT=SUM RETURN END C C************************************************************************* C SUBROUTINE RWING(KR) C C GIVES RNY=XCONT/XTOT AT FAR WING AT MU=1 C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' C I=IRAD(KR) J=JRAD(KR) DO 100 K=1,NDEP V=(Q(NQ(KR),KR)-VEL(K))/DNYD(K) CALL VOIGT(ADAMP(K,KR),V,H) PH=H/(DNYD(K)*1.772453851) Z(K)=N(I,K)*B(I,J)-N(J,K)*B(J,I) XTOT=Z(K)*HNY4P*PH/XNORM(K)+XCONT(K) RNY(K)=XCONT(K)/XTOT 100 CONTINUE RETURN END C C************************************************************************* C FUNCTION SECOND() INCLUDE 'PREC' REAL XXX(2),ETIME DUMMY=ETIME(XXX) SECOND=XXX(1) c SECOND=0.0 RETURN END C C************************************************************************* C SUBROUTINE START C C ADMINISTERS INITIALIZATIONS C C VARIABLES: C C ATOMID 4 CHARACTER IDENTIFICATION OF ATOM. C ABND ATOMIC ABUNDANCE, LOG SCALE WITH HYDROGEN=12 C AWGT ATOMIC WEIGHT. INPUT IN ATOMIC UNITS, CONVERTED TO CGS C CROUT COLLISIONAL ROUTINE NAME C C NK NUMBER OF LEVELS INCLUDING CONTINUUM LEVELS C NRAD NUMBER OF RADIATIVE TRANSITIONS TREATED IN DETAIL C NLINE NUMBER OF RADIATIVE BOUND-BOUND TRANSITIONS C NRFIX NUMBER OF TRANSITIONS WITH FIXED RATES C NQ NUMBER OF FREQUENCIES C NMU NUMBER OF ANGLES C NDEP NUMBER OF DEPTH POINTS C C EV ENERGY ABOVE GROUND STATE. INPUT IN CM-1, CONVERTED TO EV C G STATISTICAL WEIGHT OF LEVEL C LABEL 4 CHARACTER IDENTIFICATION OF LEVEL C ION IONIZATION STAGE OF LEVEL C C JRAD JRAD(KR) IS UPPER LEVEL OF RADIATIVE TRANSITION KR C IRAD IRAD(KR) IS LOWER LEVEL OF RADIATIVE TRANSITION KR C KRAD KRAD(I,J)=KRAD(J,I) IS THE NUMBER OF THE TRANSITION I-J C F OSCILLATOR STRENGTH C IWIDE =0 FOR LINES WHERE THE BACKGROUND SOURCE FUNCTION IS C TREATED AS CONSTANT, =1 FOR WIDE LINES C GA RADIATIVE DAMPING PARAMETER C GW VAN DER WAALS DAMPING PARAMETER C GQ STARK DAMPING PARAMETER C ALAMB VACUUM WAVELENGTH IN ANGSTROM C KTRANS KTRANS(KR) IS WIDE (=WIDE LINES AND CONTINUA) TRANSITION C NUMBER FOR TRANSITION KR C C A EINSTEIN A COEFFICIENT C B EINSTEIN B COEFFICIENT C C COLLISIONAL TRANSITION RATE C C JFX JFX(KR) IS UPPER LEVEL OF FIXED TRANSITION KR C IFX IFX(KR) IS LOWER LEVEL OF FIXED TRANSITION KR C IPHO =1 CONTINUUM, =0 LINE C A0 CROSSECTION AT LIMIT C TRAD BRIGHTNESS TEMPERATURE FOR CONTINUA C ITRAD RADIATION TEMPERATURE OPTION. =1 RAD.TEMP=TEMP, C =2 PHOTOSPHERIC OPTION, RAD.TEMP=TEMP OUT TO TEMP.LT.TRAD C THEN TEMP=TRAD OUTWARDS C =3 CHROMOSPHERIC OPTION, RAD.TEMP=TEMP EXCEPT WHEN C TEMP .GT. TRAD AND TEMP INCREASING OUTWARDS. C THEN RAD.TEMP=TRAD C C DNYD DOPPLER WIDTH IN UNITS OF A TYPICAL DOPPLER WIDTH C ADAMP VOIGT DAMPING PARAMETER C C N POPULATION DENSITY IN CM-3 C NSTAR LTE POPULATION DENSITY C TOTN TOTAL POPULATION DENSITY OF ATOM C C ALFA ATOMIC CROSSECTION C ALFAC PHOTOIONIZATION CROSSECTION C Z NI - GIJ*NJ C GIJ GI/GJ C C QNORM UNIT TYPICAL DOPPLER WIDTH IN KM PER SECOND AT C LINE CENTER C HNY4P H*NY/4PI NY IN UNITS OF A TYPICAL DOPPLER WIDTH C HN3C2 H*NY**3/C**2 C C Q FREQUENCY VARIABLE, IN UNITS OF A TYPICAL DOPPLER WIDTH C POSITIVE Q FOR INCREASED FREQUENCY C QMAX MAXIMUM FREQUENCY, SAME UNITS AS Q C Q0 FREQUENCY WITHIN WHICH QUADRATURE POINTS ARE DISTRIBUTED C LINEARLY INSTEAD OF LOGARITHMICALLY C IND =1 FOR ONE SIDED QUADRATURE (SYMMETRIC PROFILE) C =2 FOR TWO SIDED QUADRATURE (ASYMMETRIC PROFILE) C WQ QUADRATURE WEIGHTS FOR Q C WQMU WQ*WMU C FRQ FREQUENCY IN HZ FOR WIDE LINES AND CONTINUA C C PHI PROFILE C WPHI WEIGHTS FOR PROFILE C C XMU COSINE OF ANGLE, ARRAY C XMY COSINE OF ANGLE, USED FOR THE COMMUNICATION WITH TRANSP C WMU WEIGHTS FOR XMU C C C ATMOID 72 CHARACTER IDENTIFICATION OF ATMOSPHERE USED C DPID 72 CHARACTER IDENTIFICATION OF DEPTH-SCALE C DPTYPE =T DEPTH SCALE IS TAUSCALE, =M DEPTH SCALE IS MASS SCALE, C SEE ROUTINE ATMOS C C GRAV GRAVITATION ACCELERATION C CMASS COLUMN MASS C TEMP TEMPERATURE C NE ELECTRON DENSITY C VEL MACROSCOPIC VELOCITY INPUT IN KM PER SEC, C CONVERTED TO DOPPLER UNITS C VTURB MICROTURBULENCE VELOCITY INPUT IN KM PER SEC, C CONVERTED TO DOPPLER UNITS C HEIGHT HEIGHT ABOVE TAU5000=1 IN KILOMETERS C C BH DEPARTURE COEFFICIENT FOR HYDROGEN C NH POPULATION DENSITY FOR HYDROGEN C RHO DENSITY C C TAU STANDARD OPTICAL DEPTH C XNORM STANDARD OPACITY IN CM**2 PER CM**3 C XCONT BACKGROUND OPACITY RELATIVE TO STANDARD OPACITY C X TOTAL OPACITY RELATIVE TO STANDARD OPACITY C RNY XCONT/X C SCAT SCATTERING RELATIVE TO BACKGROUND OPACITY C C TAUQ MONOCHROMATIC OPTICAL DEPTH C TAUQQ MONOCHROMATIC OPTICAL DEPTH. USED FOR PRINTOUT IN WTAUQ C DTAUQ DTAUQ(K)=TAUQ(K)-TAUQ(K-1) C A1 DIFFUSION PARAMETER FROM TRANSP C C1 DIFFUSION PARAMETER FROM TRANSP C C S MONOCHROMATIC SOURCE FUNCTION C SL LINE SOURCE FUNCTION C SC BACKGROUND SOURCE FUNCTION, EXCLUDING SCATTERING TERM C SBCK BACKGROUND SOURCE FUNCTION, INCLUDING SCATTERING TERM C P FEAUTRIER MEAN INTENSITY IF VEL=0, ELSE MODIFIED MEAN,SEE C ROUTINE TRANSP C PMS P - S C IPLUS I+ C IMINUS I- C JNY MEAN INTENSITY C JNYOLD MEAN INTENSITY FROM PREVIOUS ITERATION C BP PLANCK FUNCTION C C EJ MAXIMUM CORRECTION IN JNY FOR EACH TRANSITION AND EACH DEPTH C EMAXJ MAXIMUM CORRECTION IN JNY OVER ALL DEPTHS AND TRANSITIONS C C DLGTMX MAXIMUM DELTA LG TAUQ FOR EACH DEPTH C DLGSMX MAXIMUM DELTA LG S FOR EACH DEPTH C DLGPMX MAXIMUM CHANGE IN P-S FOR EACH DEPTH C C IT ITERATION NUMBER C C RIJ RATE FROM I TO J PER NI ATOM C RJI RATE FROM J TO I PER NJ ATOM C C FLUX MONOCHROMATIC PHYSICAL FLUX IN CGS UNITS C OUTINT MONOCHROMATIC SURFACE INTENSITY IN CGS UNITS C COOL COOLING FUNCTION IN ERG/CM**3/S C WEQ EQUIVALENT WIDTH C WEQLTE EQUIVALENT WIDTH IN LTE C C E ERROR TERM IN EQUATION FOR POPULATIONS C W W(I,K,J,L) MATRIX ELEMENT DESCRIBING CHANGES IN C E(I,K) DUE TO RELATIVE CHANGES IN POPULATIONS, C DN(J,L)/N(J,L) C C PHYSICAL CONSTANTS: C C EE ELECTRON CHARGE C EM ELECTRON MASS C HH PLANCK CONSTANT C CC VELOCITY OF LIGHT C BK BOLTZMANN CONSTANT C UU UNIVERSAL MASS CONSTANT C HCE HH*CC/EE*1.E8 LAMBDA(ANGSTROM)=HCE/ENERGY(EV) C HC2 2*HH*CC*1.E24 2*H*NY**3/C**2=HC2/LAMBDA(ANGSTROM)**3 C HCK HH*CC/BK*1.E8 H*NY/KT=HCK/LAMBDA(ANGSTROM)/T C EK EE/BK C PI PI C C INPUT PARAMETERS: C C ITMAX MAXIMUM NUMBER OF ITERATIONS C ELIM1 WHEN EMAX.LT.ELIM1 NO MORE RECALCULATIONS OF MATRIX W C ARE DONE C ELIM2 WHEN EMAX.LT.ELIM2 NO MORE ITERATIONS ARE DONE C DIFF WHEN DTAUQ .GT. DIFF THE DIFFUSION APPROXIMATION IS USED C ISUM THE STATISTICAL EQUILIBRIUM EQUATION FOR LEVEL ISUM IS C REPLACED BY THE EQUATION OF CONSERVATION OF NUMBERS C ISTART STARTING APPROXIMATION =-1 FROM FILE =0 I=0, =1 I=B C ILAMBD NUMBER OF LAMBDA ITERATIONS AFTER THE START APPROXIMATION C ITRAN METHOD USED FOR FORMAL SOLUTION: C =0 STANDARD FEAUTRIER C =1 FEAUTRIER, SPLINE FORMULATION C =2 FEAUTRIER, HERMITE FORMULATION C ICONV =0 ONLY CONVERGED SOLUTION WILL GIVE OUTPRINT C IOPAC =0 OPACITIES READ FROM FILE OPC C THIN RADIATION IS LAMBDA ITERATED IF TAUQ.LT.THIN C IHSE =0 NO HYDROSTATIC EQUILIBRIUM EQUATION (HSE) INTEGRATION C =1 DO HSE INTEGRATION IF ATOMID='H ' C IATOM2 DO A SECOND FORMAL SOLUTION WITH ATOMIC MODEL FROM FILE C ATOM2 WITH POPULATIONS FROM ATOM REDISTRIBUTED ACCORDING C TO THE STATISTICAL WEIGHTS C C PRINTOUT OPTION PARAMETERS: C C IWXXXX =0 INHIBITS PRINTOUT IN ROUTINE XXXX C .GT.0 WILL GIVE PRINTOUT, WITH DEPTH STEP = IWXXXX EXCEPT C IWHEAD,IWATOM,IWABND,IWEMAX,IWNIIT,IWHSE,IWWMAT C IWTEST .LT.0 WILL GIVE TEST PRINTOUTS AFTER EVERY ITERATION (LOTS) C .GT.0 WILL GIVE TEST PRINTOUTS ONLY AFTER CONVERGENCE C IWWMAT .LT.0 WILL GIVE MATRIX PRINTOUT AFTER EVERY ITERATION C .GT.0 WILL GIVE ONLY ONE MATRIX PRINTOUT C IWCHAN .LT.0 WILL GIVE CONVERGENCE MAP FOR BOTH POPULATIONS AND C MEAN INTENSITIES C .GT.0 WILL GIVE CONVERGENCE MAP ONLY FOR POPULATIONS C C C CODED BY: M.CARLSSON (SEP.1983 - JAN.1984) C MODIFIED BY M.CARLSSON JAN 1985 C C FILE STRUCTURE: C C (INPUT) MEANS THAT THIS INPUT FILE IS ONLY READ IF THE CORRESPONDING C INPUT PARAMETER IS SET .GT. 0. C C UNIT C LINPUT INPUT - INPUT : INPUT PARAMETERS C LATOM ATOM - INPUT : ATOMIC MODEL C LATMOS ATMOS - INPUT : ATMOSPHERIC STRUCTURE C LDSCAL DSCALE - INPUT : LOGARITHMIC COLUMN MASS SCALE OR C LOG TAU SCALE C LABUND ABUND - INPUT : ABUNDANCES C LATOM2 ATOM2 -(INPUT): ATOMIC MODEL 2 FOR FORMAL SOLUTION C LRSTRT RSTRT -(INPUT): EXPLICIT START APPROXIMATION C LJOBLO JOBLOG - OUTPUT: JOBLOG MESSAGES C LOUT OUT - OUTPUT: INPUT PARAMETERS, ATOMIC PARAMETERS, C ABUNDANCES, OPACITIES, ATMOSPHERIC MODEL, C LTE POPULATIONS, PROFILE PARAMETERS, C START APPROXIMATION, ITERATION PLOT, C POPULATIONS, IONIZATION FRACTIONS, C LINE PARAMETERS, FLUXES, INTENSITIES, C JBAR,SL,RATES C LTIME TIME - OUTPUT: TIMING OF ROUTINES C LRSTRT RSTRT2 - OUTPUT: INTERPOLATED ATMOSPHERE AND RESULTING C POPULATIONS TO BE USED AS START C APPROXIMATION FROM PARTIALLY CONVERGED C SOLUTION C LATHSE ATHSE - OUTPUT: NEW ATMOSPHERE RESULTING FROM HSE C INTEGRATIONS IN HSEINT C LDSCA2 DSCAL2 - OUTPUT: INTERPOLATED LG MASS SCALE SMOOTH IN C DELTA LG TAUNY C LWMAT WMAT - OUTPUT: COEFFICIENT MATRIX C LNIIT NIIT - OUTPUT: NI FOR EVERY ITERATION UNFORMATTED C LDUMS DUMS - : DUMMY FILE USED FOR STRIPPING INPUT C FILES FROM COMMENT LINES = LINES C BEGINNING WITH * C LDUMC DUMC - : DUMMY FILE USED IN ROUTINE COLRAT C LDUMI DUMI - : DUMMY FILE USED IN ROUTINE RINPUT C LOPC OPC - : FILE USED TO STORE OPACITIES C LJNY JNY - : FILE USED TO STORE JNY, SEE ROUTINE OPAC C LINIT INIT - : FILE USED TO STORE RHO, TOTN, XNORM, BH C LPHI PHI - : FILE USED TO STORE PHI C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLGMX' INCLUDE 'CLU' INCLUDE 'COPCL' C C INITIALISE C CALL CPUTIME('START ',0,0,3) CALL CPUTIME('INPUT ',0,0,2) C C OPEN GLOBAL FILES ALWAYS NEEDED C CALL OPEN(LOUT,'OUT',1,'NEW') CALL OPEN(LTIME,'TIME',1,'UNKNOWN') CALL OPEN(LJOBLO,'JOBLOG',1,'UNKNOWN') CALL OPEN(LDUMS,'DUMS',1,'UNKNOWN') CALL OPEN(LPHI,'PHI',0,'UNKNOWN') C C READ INPUT PARAMETERS C CALL RINPUT C C OPEN GLOBAL FILES REGULATED BY INPUT PARAMETERS C IF(IWWMAT.NE.0) CALL OPEN(LWMAT,'WMAT',0,'NEW') IF(IWNIIT.NE.0) CALL OPEN(LNIIT,'NIIT',0,'NEW') C CALL ATMOS CALL ATOM(1) CALL VALCHK CALL WHEAD CALL WATOM(1) CALL FREQ CALL OPAC(0) CALL DPCONV CALL WATMOS CALL LTEPOP CALL COLRAT CALL FIXRAD CALL WLTE CALL WNIIT(0) CALL DAMP XMU0=0. XMU1=1. CALL GAUSI(NMU,XMU0,XMU1,WMU,XMU) CALL CPUTIME('INPUT ',0,2,2) C CALL PROFIL C CALL LTEEQW CALL INITIA CALL WSTART CALL WNIIT(1) CALL CPUTIME('START ',0,2,3) C RETURN END C C********************************************************************** C SUBROUTINE STATEQ(ISUM,NETRAT) C C SOLVES THE EQUATIONS OF STATISTICAL EQUILIBRIUM FOR GIVEN RATES. C NETRAT=0 RADIATIVE DETAILED BALANCE C: C: STATEQ 87-12-10 MODIFICATIONS: (MATS CARLSSON) C: IF ISUM=0 ISUM IS SET TO THE LEVEL WITH THE LARGEST C: POPULATION AT EACH DEPTH-POINT. THIS MODIFICATION ALSO C: AFFECTS VALCHK AND ITER C: C: 88-04-13 MODIFICATIONS PGJ C: INCLUDES A CHECK ON AA(I,I) TO SEE IF THE NET RATES OUT OF A LEVEL C: ARE ZERO C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' C* INCLUDE 'CLU' C* C DIMENSION AA(MK,MK),Y(MK) LOGICAL LISUM C LISUM=ISUM.EQ.0 DO 200 K=1,NDEP DO 120 I=1,NK DO 110 J=1,NK AA(I,J)=0. 110 CONTINUE 120 CONTINUE C C RADIATION PART C IF(NETRAT.NE.0) THEN DO 130 KR=1,NRAD I=IRAD(KR) J=JRAD(KR) AA(I,J)=AA(I,J)+RJI(K,KR) AA(J,I)=AA(J,I)+RIJ(K,KR) AA(I,I)=AA(I,I)-AA(J,I) AA(J,J)=AA(J,J)-AA(I,J) 130 CONTINUE ENDIF C C FIXED TRANSITIONS PART C DO 150 I=1,NK DO 140 J=1,NK IF(J.EQ.I) GOTO 140 AA(I,J)=AA(I,J)+C(J,I,K) AA(I,I)=AA(I,I)-C(I,J,K) 140 CONTINUE 150 CONTINUE IF(LISUM) THEN POPMAX=N(1,K) ISUM=1 DO 152 I=2,NK IF(N(I,K).GT.POPMAX) THEN POPMAX=N(I,K) ISUM=I ENDIF 152 CONTINUE ENDIF DO 155 I=1,NK AA(ISUM,I)=1. Y(I)=0. C* IF (AA(I,I) .EQ. 0.0) THEN WRITE (LJOBLO,1000) I,K 1000 FORMAT (' STATEQ: LEVEL',I3, * ' HAS NET OUTWARD RATE ZERO',' AT DEPTH',I3) END IF C* 155 CONTINUE Y(ISUM)=TOTN(K) C CALL EQSYST(NK,MK,AA,Y,.TRUE.) C DO 160 I=1,NK N(I,K)=Y(I) 160 CONTINUE 200 CONTINUE IF(LISUM) ISUM=0 RETURN END C C*********************************************************************** C SUBROUTINE STOP(TEXT) C C CAUSES A STOP WITH A TEXT WRITTEN TO UNIT LJOBLO C CLOSES ALL OPEN FILES C INCLUDE 'PREC' INCLUDE 'CLU' INCLUDE 'COPCL' C CHARACTER*(*) TEXT C WRITE(LJOBLO,100) TEXT 100 FORMAT(1X,A) DO 200 I=1,MAXLU IF(LU(I)) CLOSE(I) 200 CONTINUE C STOP C END C C*********************************************************************** C SUBROUTINE TAUNYQ(NY,KR) C C CALCULATES TAUNY FOR MU=1 C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CCONST' INCLUDE 'CTAUQQ' C I=IRAD(KR) J=JRAD(KR) TAU(0)=-TAU(1) TAUQ(0)=0.0 X(0)=0.0 IF(KR.LE.NLINE) THEN DO 100 K=1,NDEP V=(Q(NY,KR)-VEL(K))/DNYD(K) CALL VOIGT(ADAMP(K,KR),V,H) PHI(K)=H/(DNYD(K)*1.772453851) Z(K)=N(I,K)-G(I)/G(J)*N(J,K) ALFA(K)=B(I,J)*PHI(K)*HNY4P X(K)=Z(K)*ALFA(K)/XNORM(K)+XCONT(K) TAUQ(K)=TAUQ(K-1)+0.5*(X(K-1)+X(K))*(TAU(K)-TAU(K-1)) TAUQQ(NY,K)=TAUQ(K) 100 CONTINUE ELSE KT=KTRANS(KR) DO 200 K=1,NDEP GIJ(K)=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FRQ(NY,KT)/BK/TEMP(K)) ALFA(K)=ALFAC(NY,KR-NLINE) Z(K)=N(I,K)-GIJ(K)*N(J,K) X(K)=Z(K)*ALFA(K)/XNORM(K)+XCONT(K) TAUQ(K)=TAUQ(K-1)+0.5*(X(K-1)+X(K))*(TAU(K)-TAU(K-1)) TAUQQ(NY,K)=TAUQ(K) 200 CONTINUE ENDIF RETURN END C C*********************************************************************** C SUBROUTINE TRCONT C C SOLVES THE EQUATION OF RADIATIVE TRANSFER FOR THE BACKGROUND C: C: TRCONT 88-01-21 MODIFICATIONS: (MATS CARLSSON) C: STORES CONTINUUM INTENSITY FOR USE IN CONTRIBUTION FUNCTIONS C: C: 93-01-21 MODIFICATIONS: (JO BRULS, MATS CARLSSON) C: READ OPACITY FILES AT LINE WING FOR IWIDE=1. CORRECTED. C: C: 95-03-11 MODIFICATIONS: (MATS CARLSSON) C: MODIFIED TO WORK WITH TRANSC (SCATTERING FEAUTRIER) C: ABSORPTION PART OF SOURCE FUNCTION IS GIVEN BY C: S(K)-RNY(K)*SBCK(K)+RNY(K)*SC(K) C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CLU' INCLUDE 'CCNTRB' C CALL CPUTIME('TRCONT',0,0,1) CALL REWIND(LOPC) IREC=0 DO 400 KR=1,NLINE FLUX(0,KR)=0. C C FIND NY FOR LINECENTER AND LINEWING C NY0=1 NY1=1 QMIN=ABS(Q(1,KR)) QMX=QMIN DO 50 NY=1,NQ(KR) IF(ABS(Q(NY,KR)).LT.QMIN) THEN NY0=NY QMIN=ABS(Q(NY,KR)) ENDIF IF(ABS(Q(NY,KR)).GT.QMX) THEN NY1=NY QMX=ABS(Q(NY,KR)) ENDIF 50 CONTINUE C C READ OPACITY FILES AT LINE-CENTER, JNY FILE AT LINE WING C DO 200 NY=1,NY0 CALL READX 200 CONTINUE CALL READJ(IREC+NY1) DO 210 K=1,NDEP X(K)=XCONT(K) RNY(K)=1.0 SBCK(K)=SC(K)+SCAT(K)*JNY(K) S(K)=SBCK(K) 210 CONTINUE DO 300 MU=1,NMU XMY=XMU(MU) C CALL TRANSP IF(MU.EQ.NMU) THEN DO 215 K=1,NDEP ICONT(KR,K)=IPLUS(K) 215 CONTINUE ENDIF C OUTINT(0,MU,KR)=IPLUS(0) FLUX(0,KR)=FLUX(0,KR)+WMU(MU)*XMY*IPLUS(0) 300 CONTINUE FLUX(0,KR)=2.*PI*FLUX(0,KR) C C POSITION FILES FOR NEXT TRANSITION C DO 350 NY=NY0+1,NQ(KR) CALL READX 350 CONTINUE IREC=IREC+NQ(KR) 400 CONTINUE CALL CPUTIME('TRCONT',0,1,1) RETURN END C C*********************************************************************** C SUBROUTINE TRPT C C SOLVES THE EQUATION OF RADIATIVE TRANSFER FOR GIVEN POPULATIONS C GIVES FLUXES AND INTENSITIES AND RADIATIVE RATES C JNY IS CALCULATED AND USED TO UPDATE BACKGROUND SOURCE FUNCTION, C SCATTERING TERM C C: TRPT 92-10-08 MODIFICATIONS: (MATS CARLSSON) C: INCIDENT RADIATION FIELD IS SET IF ITRAN.GE.10 C: C: 95-08-16 MODIFICATIONS: (MATS CARLSSON) C: BOTH OLD PRINTOUT ROUTINES (TO FILE OUT) AND IDL PRINTOUT C: ROUTINES ARE CALLED C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION SWITCHED ON BY INCRAD.NE.0 INSTEAD OF C: ITRAN=10-14 C: C: 95-08-21 MODIFICATIONS: (MATS CARLSSON) C: WARNING MESSAGES REGULATED BY IWARN C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CCONST' INCLUDE 'CLGMX' INCLUDE 'CLU' INCLUDE 'CIMIN' C LOGICAL CONT C CALL CPUTIME('TRPT ',0,0,1) CALL REWIND(LOPC) CALL REWIND(LPHI) IREC=0 EMAXJ=0.0 IF(ICONV.EQ.1) CALL MXDLG(0) DO 400 KR=1,NRAD I=IRAD(KR) J=JRAD(KR) C C CALCULATE SOME CONSTANTS C DO 270 K=1,NDEP RIJ(K,KR)=0. RJI(K,KR)=0. COOL(K,KR)=0. EJ(KR,K)=0. 270 CONTINUE IF(KR.GT.NLINE) THEN CONT=.TRUE. DO 275 K=1,NDEP WPHI(K,KR)=1.0 275 CONTINUE ELSE CONT=.FALSE. GIJK=G(I)/G(J) HN3C2=A(KR)/B(J,I) DO 280 K=1,NDEP Z(K)=N(I,K)-GIJK*N(J,K) IF(Z(K).LE.0.0) THEN IF(ICONV.EQ.1 .AND. IWARN.GE.2) THEN WRITE(LJOBLO,278) KR,K 278 FORMAT(' WARNING IN TRPT: NEGATIVE OPACITIES', * ' IN LINE',I3,' DEPTH',I4) ENDIF ENDIF SL(K,KR)=HN3C2*N(J,K)*GIJK/Z(K) GIJ(K)=GIJK 280 CONTINUE ENDIF IF(IWIDE(KR)) KT=KTRANS(KR) C DO 380 NY=1,NQ(KR) IREC=IREC+1 CALL READX CALL READJ(IREC) FLUX(NY,KR)=0. DO 285 K=1,NDEP SBCK(K)=SC(K)+SCAT(K)*JNY(K) JNYOLD(K)=JNY(K) JNY(K)=0.0 285 CONTINUE IF(CONT) THEN HN3C2=2.*HH*FRQ(NY,KT)/CC*FRQ(NY,KT)/CC*FRQ(NY,KT) DO 290 K=1,NDEP GIJ(K)=NSTAR(I,K)/NSTAR(J,K)* * EXP(-HH*FRQ(NY,KT)/BK/TEMP(K)) ALFA(K)=ALFAC(NY,KR-NLINE) Z(K)=N(I,K)-GIJ(K)*N(J,K) SL(K,KR)=HN3C2*N(J,K)*GIJ(K)/Z(K) 290 CONTINUE ENDIF C DO 370 MU=1,NMU XMY=XMU(MU) WQMU=WQ(NY,KR)*WMU(MU) IF(.NOT.CONT .AND. (IND(KR).EQ.2 .OR. MU.EQ.1)) THEN CALL READP DO 300 K=1,NDEP ALFA(K)=B(I,J)*PHI(K)*HNY4P 300 CONTINUE ENDIF IF(IND(KR).EQ.2 .OR. MU.EQ.1) THEN DO 310 K=1,NDEP X(K)=Z(K)*ALFA(K)/XNORM(K)+XCONT(K) RNY(K)=XCONT(K)/X(K) S(K)=(1.-RNY(K))*SL(K,KR)+RNY(K)*SBCK(K) 310 CONTINUE ENDIF C IF(INCRAD.NE.0) IMINUS(0)=XIMIN(NY,MU,KR) CALL TRANSP C IF(ICONV.EQ.1) CALL MXDLG(1) OUTINT(NY,MU,KR)=IPLUS(0) FLUX(NY,KR)=FLUX(NY,KR)+WMU(MU)*XMY*IPLUS(0) IF(IWIDE(KR)) THEN WQMUH=WQMU/HNY4P/FRQ(NY,KT)*FRQ(0,KT) ELSE WQMUH=WQMU/HNY4P ENDIF DO 360 K=1,NDEP RIJ(K,KR)=RIJ(K,KR)+WQMUH*WPHI(K,KR)*ALFA(K)*P(K) RJI(K,KR)=RJI(K,KR)+WQMUH*WPHI(K,KR)*ALFA(K)* * GIJ(K)*(P(K)+HN3C2) JNY(K)=JNY(K)+WMU(MU)*P(K) COOL(K,KR)=COOL(K,KR)-WQMU*X(K)*PMS(K) 360 CONTINUE 370 CONTINUE IF(ICONV.EQ.1) CALL WTEST(KR,NY) IF(ICONV.EQ.1) CALL WIDLNY(KR,NY) FLUX(NY,KR)=2.*PI*FLUX(NY,KR) DO 375 K=1,NDEP IF(JNY(K).GT.0.0) THEN DJ=(JNY(K)-JNYOLD(K))/JNY(K) ELSE DJ=0.0 ENDIF EMAXJ=MAX(ABS(DJ),EMAXJ) IF(ABS(DJ).GT.ABS(EJ(KR,K))) EJ(KR,K)=DJ 375 CONTINUE CALL WRITEJ(IREC) 380 CONTINUE DO 390 K=1,NDEP COOL(K,KR)=4.0*PI*COOL(K,KR)*XNORM(K)*QNORM*1.E13/ALAMB(KR) 390 CONTINUE 400 CONTINUE CALL CPUTIME('TRPT ',0,1,1) RETURN END C C*********************************************************************** C SUBROUTINE VALCHK C C CHECKS INPUT PARAMETER VALIDITY C: C: VALCHK 88-05-03 MODIFICATIONS: (MATS CARLSSON) C: ALLOWS ISUM=0, IOPAC=2,3,4 C: NO CHECK ON IWLINE, IWRAD, ICONV C: C: 95-08-23 MODIFICATIONS: (MATS CARLSSON) C: NEW CHECKS INTRODUCED C: C: 95-10-13 MODIFICATIONS: (MATS CARLSSON) C: CHECKS FOR NON-ALLOWED COMBINATION OF NON-SYMMETRIC C: FREQUENCY QUADRATURE AND VELOCITY FIELDS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CTERM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CLU' C LOGICAL ERROR C ERROR=.FALSE. IF(QNORM.LE.0.0 .OR. QNORM.GT.100.) THEN WRITE(LJOBLO,100) 'QNORM ',QNORM ERROR=.TRUE. ENDIF IF(DIFF.LE.0.0 .OR. DIFF.GT.1000.) THEN WRITE(LJOBLO,100) 'DIFF ',DIFF ERROR=.TRUE. ENDIF IF(ELIM1.LE.0.0 .OR. ELIM1.GT.1.0) THEN WRITE(LJOBLO,100) 'ELIM1 ',ELIM1 ERROR=.TRUE. ENDIF IF(ELIM2.LE.0.0 .OR. ELIM2.GT.1.0) THEN WRITE(LJOBLO,100) 'ELIM2 ',ELIM2 ERROR=.TRUE. ENDIF IF(THIN.LT.0.0 .OR. THIN.GT.10.0) THEN WRITE(LJOBLO,100) 'THIN ',THIN ERROR=.TRUE. ENDIF IF(NMU.LE.0 .OR. NMU.GT.MMU) THEN WRITE(LJOBLO,200) 'NMU ',NMU ERROR=.TRUE. ENDIF IF(ITMAX.LT.0) THEN WRITE(LJOBLO,200) 'ITMAX ',ITMAX ERROR=.TRUE. ENDIF IF(ISUM.LT.0 .OR. ISUM.GT.NK) THEN WRITE(LJOBLO,200) 'ISUM ',ISUM ERROR=.TRUE. ENDIF IF(ILAMBD.LT.0) THEN WRITE(LJOBLO,200) 'ILAMBD',ILAMBD ERROR=.TRUE. ENDIF IF(ICONV.GT.1) THEN WRITE(LJOBLO,200) 'ICONV ',ICONV WRITE(LJOBLO,400) 'COLLISIONAL RADIATIVE SWITCHING NOW'// * ' SWITCHED ON BY ICRSW' ERROR=.TRUE. ENDIF IF(IOPAC.GT.10) THEN WRITE(LJOBLO,200) 'IOPAC ',IOPAC WRITE(LJOBLO,400) 'LINE ABSORPTION NOW SWITCHED ON BY IOPACL' ERROR=.TRUE. ENDIF IF(IOPAC.LT.0 .OR. IOPAC.GT.4) THEN WRITE(LJOBLO,200) 'IOPAC ',IOPAC ERROR=.TRUE. ENDIF IF(ITRAN.GE.10) THEN WRITE(LJOBLO,200) 'ITRAN ',ITRAN WRITE(LJOBLO,400) 'INCIDENT RADIATION NOW SWITCHED ON BY INCRAD' ERROR=.TRUE. ENDIF IF(ITRAN.LT.0 .OR. ITRAN.GT.4) THEN WRITE(LJOBLO,200) 'ITRAN ',ITRAN ERROR=.TRUE. ENDIF IF(IWATMS.LT.0) THEN WRITE(LJOBLO,200) 'IWATMS',IWATMS ERROR=.TRUE. ENDIF IF(IWLTE.LT.0) THEN WRITE(LJOBLO,200) 'IWLTE ',IWLTE ERROR=.TRUE. ENDIF IF(IWDAMP.LT.0) THEN WRITE(LJOBLO,200) 'IWDAMP',IWDAMP ERROR=.TRUE. ENDIF IF(IWSTRT.LT.0) THEN WRITE(LJOBLO,200) 'IWSTRT',IWSTRT ERROR=.TRUE. ENDIF IF(IWN.LT.0) THEN WRITE(LJOBLO,200) 'IWN ',IWN ERROR=.TRUE. ENDIF IF(IWTAUQ.LT.0) THEN WRITE(LJOBLO,200) 'IWTAUQ',IWTAUQ ERROR=.TRUE. ENDIF IF(IWLGMX.LT.0) THEN WRITE(LJOBLO,200) 'IWLGMX',IWLGMX ERROR=.TRUE. ENDIF IF(IWRATE.LT.0) THEN WRITE(LJOBLO,200) 'IWRATE',IWRATE ERROR=.TRUE. ENDIF IF(IWEVEC.LT.0) THEN WRITE(LJOBLO,200) 'IWEVEC',IWEVEC ERROR=.TRUE. ENDIF IF(IWNIIT.LT.0) THEN WRITE(LJOBLO,200) 'IWNIIT',IWNIIT ERROR=.TRUE. ENDIF IF(IWHSE.LT.0) THEN WRITE(LJOBLO,200) 'IWHSE ',IWHSE ERROR=.TRUE. ENDIF IF(IWRAD.LT.0) THEN WRITE(LJOBLO,200) 'IWRAD ',IWRAD WRITE(LJOBLO,400) 'JFIX WRITING NOW SWITCHED ON BY IWJFIX.NE.0' ERROR=.TRUE. ENDIF IF(ATOMID.EQ.'H ' .AND. IHSE.NE.0 .AND. DPTYPE.EQ.'T') THEN WRITE(LJOBLO,300) ERROR=.TRUE. ENDIF DO 50 KR=1,NLINE IF(NTERM(KR).GT.1) THEN DO 40 K=1,NDEP IF(ABS(VEL(K)).GT.1.E-6) THEN CALL STOP('VALCHK: BLENDS + VELOCITY FIELDS NOT ALLOWED') ENDIF 40 CONTINUE ENDIF 50 CONTINUE C 100 FORMAT(' VALUE OUT OF RANGE ',A,' = ',1P,E10.2) 200 FORMAT(' VALUE OUT OF RANGE ',A,' = ',I10) 300 FORMAT(' DEPTH SCALE HAS TO BE A COLUMN MASS SCALE'/ * 'IF HSE INTEGRATIONS ARE TO BE DONE') 400 FORMAT(A) C IF(ERROR) CALL STOP('VALCHK: VALUE(S) OUT OF RANGE') C RETURN END C C*********************************************************************** C SUBROUTINE VOIGT(A,VS,H) C C COMPUTES A VOIGT FUNCTION H = H(A,V) C A=GAMMA/(4*PI*DNUD) AND V=(NU-NU0)/DNUD. THIS IS DONE AFTER C TRAVING (LANDOLT-B\RNSTEIN, P. 449). C C CODED BY: KJELL ERIKSSON (APR-1974). REVISED: (JULY-1975). C INCLUDE 'PREC' DIMENSION AK(19),A1(5) DATA AK /-1.12470432, -0.15516677, 3.28867591, -2.34357915, , 0.42139162, -4.48480194, 9.39456063, -6.61487486, 1.98919585, , -0.22041650, 0.554153432, 0.278711796,-0.188325687, 0.042991293, ,-0.003278278, 0.979895023,-0.962846325, 0.532770573,-0.122727278/ DATA SQP/1.772453851/,SQ2/1.414213562/ C V = ABS(VS) U = A + V V2 = V*V IF (A.EQ.0.0) GO TO 140 IF (A.GT.0.2) GO TO 120 IF (V.GE.5.0) GO TO 121 C EX=0. IF(V2.LT.100.)EX = EXP(-V2) K = 1 C 100 QUO = 1. IF (V.LT.2.4) GO TO 101 QUO = 1./(V2 - 1.5) M = 11 GO TO 102 C 101 M = 6 IF (V.LT.1.3) M = 1 102 DO 103 I=1,5 A1(I) = AK(M) M = M + 1 103 CONTINUE H1 = QUO*(A1(1) + V*(A1(2) + V*(A1(3) + V*(A1(4) + V*A1(5))))) IF (K.GT.1) GO TO 110 C C A LE 0.2 AND V LT 5. C H = H1*A + EX*(1. + A*A*(1. - 2.*V2)) RETURN C 110 PQS = 2./SQP H1P = H1 + PQS*EX H2P = PQS*H1P - 2.*V2*EX H3P = (PQS*(1. - EX*(1. - 2.*V2)) - 2.*V2*H1P)/3. + PQS*H2P H4P = (2.*V2*V2*EX - PQS*H1P)/3. + PQS*H3P PSI = AK(16) + A*(AK(17) + A*(AK(18) + A*AK(19))) C C 0.2 LT A LE 1.4 AND A + V LE 3.2 C H = PSI*(EX + A*(H1P + A*(H2P + A*(H3P + A*H4P)))) RETURN C 120 IF (A.GT.1.4.OR.U.GT.3.2) GO TO 130 EX=0. IF(V2.LT.100.)EX = EXP(-V2) K = 2 GO TO 100 C C A LE 0.2 AND V GE 5. C 121 H = A*(15. + 6.*V2 + 4.*V2*V2)/(4.*V2*V2*V2*SQP) RETURN C 130 A2 = A*A U = SQ2*(A2 + V2) U2 = 1./(U*U) C C A GT 1.4 OR A + V GT 3.2 C H = SQ2/SQP*A/U*(1. + U2*(3.*V2 - A2) + , U2*U2*(15.*V2*V2 - 30.*V2*A2 + 3.*A2*A2)) RETURN C C A EQ 0. C 140 H=0. IF(V2.LT.100.)H=EXP(-V2) RETURN END C C******************************************************************** C SUBROUTINE WATHSE C C WRITES THE ATMOSPHERE AND HYDROGEN POPULATIONS TO FILE ATHSE C TO BE USED AS NEW ATMOSPHERE AFTER HSE INTEGRATIONS C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IHSE.EQ.0) RETURN C CALL OPEN(LATHSE,'ATHSE',1,'UNKNOWN') WRITE(LATHSE,40) 40 FORMAT('* ATMOSPHERE FROM HSE INTEGRATIONS') WRITE(LATHSE,50) ATMOID 50 FORMAT(1X,A) IF(DPTYPE.EQ.'M') THEN WRITE(LATHSE,70) 70 FORMAT(' MASS SCALE') ELSE WRITE(LATHSE,80) 80 FORMAT(' TAU5000 SCALE') ENDIF WRITE(LATHSE,90) LOG10(GRAV) 90 FORMAT(F10.7) WRITE(LATHSE,100) NDEP 100 FORMAT(1X,I3) IF(DPTYPE.EQ.'M') THEN WRITE(LATHSE,200) (LOG10(CMASS(K)),TEMP(K),NE(K),VEL(K)*QNORM, * VTURB(K)*QNORM,K=1,NDEP) ELSE WRITE(LATHSE,200) (LOG10(TAU(K)),TEMP(K),NE(K),VEL(K)*QNORM, * VTURB(K)*QNORM,K=1,NDEP) ENDIF 200 FORMAT(1P,5E15.6) DO 400 K=1,NDEP WRITE(LATHSE,300) (NH(I,K),I=1,6) 300 FORMAT(1P,6E12.4) 400 CONTINUE CALL CLOSE(LATHSE) C RETURN END C C******************************************************************** C SUBROUTINE WATMOS C C WRITES ATMOSPHERIC PARAMETERS TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CINPUT' INCLUDE 'CLU' C* C* 1987 P. JUDGE. ALTERATION INCLUDE 'CCONST' C* LOCALS C* DIMENSION H1COL(MDEP),HKCOL(MDEP) C* END ALTERATION C* DATA XMH/1.67352E-24/,XMP/1.672661E-24/ C IF(IWATMS.EQ.0) RETURN C CALL CPUTIME('WATMOS',0,0,1) WRITE(LOUT,100) ATMOID 100 FORMAT('1 ATMOSPHERIC MODEL',3X,A) IF(DPTYPE.EQ.'M') THEN WRITE(LOUT,110) 110 FORMAT(' MASS SCALE WAS USED AS INPUT') ELSE WRITE(LOUT,120) 120 FORMAT(' TAU5000 SCALE WAS USED AS INPUT') ENDIF C* C* ALTERED BY PGJ (MAY 1987) TO OUTPUT THE HYDROGEN POPULATIONS C* ALTERED BY PGJ (JAN 1988) TO OUTPUT THE TEMPERATURE GRADIENT AND C* GAS PRESURES C* WRITE (LOUT,150) LOG10(GRAV) 150 FORMAT (' LOG GRAVITY=',F12.5,/,/,2X,'K',5X,'LG M',6X,'LG TAU', * 5X,'T',7X,'NHI',6X,'NHII',8X,'NE',7X,'VEL',4X,'VTURB',6X, * 'BH 1',7X,'HEIGHT',7X,'RHO',/,67X,'KM/S',29X,'KM',7X, * 'G/CM3') WRITE (LOUT,200) (K,LOG10(CMASS(K)),LOG10(TAU(K)),TEMP(K), * (NH(1,K)+NH(2,K)+NH(3,K)+NH(4,K)+NH(5,K)),NH(6,K),NE(K), * VEL(K)*QNORM,VTURB(K)*QNORM,BH(1,K),HEIGHT(K),RHO(K),K=1,NDEP, * IWATMS) 200 FORMAT (1X,I3,0P,2F10.5,F10.1,1P,3E10.3,0P,2F8.3,1P,E10.3,1P, * E13.5,1P,E10.3) C C GAS DATA C WRITE (LOUT,160) 160 FORMAT (/,/,' GAS DATA (ASSUMING HE=0.1*H, METALS = 0)',/, * ' DP/DH = P/HSCALE WHERE HSCALE = PRESSURE SCALE HEIGHT',/, * /,2X,'K LG M T DT/DH HSCALE ', * 'P=MG P(H+HE+E) ','PHYD PE PTURB PGAS/PTOT', * ' H(N=1) H(CONT)',/,2X, * ' G/CM2 K K/KM KM ', * 'DYNE/CM2... ',40X,'NUMBER COLUMN (/CM2)') C* C* TAKE TOP POINT AND MAKE INITIAL NUMBER COLUMN DENSITIES IN RATIO C* OF TOP IONIZATION FRACTION: C* H1COL(1) = CMASS(1)*NH(1,1)/ (NH(1,1)+NH(6,1))/1.36250/XMH HKCOL(1) = CMASS(1)*NH(6,1)/ (NH(1,1)+NH(6,1))/1.36250/XMP DO 1010 K = 2,NDEP H1COL(K) = H1COL(K-1) + (NH(1,K)+NH(1,K-1))*0.5E5* * (HEIGHT(K-1)-HEIGHT(K)) HKCOL(K) = HKCOL(K-1) + (NH(6,K)+NH(6,K-1))*0.5E5* * (HEIGHT(K-1)-HEIGHT(K)) 1010 CONTINUE C DO 101 K = 1,NDEP,IWATMS PHYD = (NH(1,K)+NH(2,K)+NH(3,K)+NH(4,K)+NH(5,K)+NH(6,K))*BK* * TEMP(K) PHE = 0.1*PHYD PELEC = NE(K)*BK*TEMP(K) PV = 0.5*RHO(K)*VTURB(K)*VTURB(K)*QNORM*QNORM*1.E10 PTOT = PHYD + PHE + PELEC + PV PTOT1 = GRAV*CMASS(K) DTDH = DIFF1(K,TEMP,HEIGHT,NDEP) SCHT = PTOT/GRAV/RHO(K)/1.E5 WRITE (LOUT,210) K,LOG10(CMASS(K)),TEMP(K),DTDH,SCHT, * PTOT1,PTOT,PHYD,PELEC,PV, (PHYD+PELEC+PHE)/PTOT,H1COL(K), * HKCOL(K) 210 FORMAT (1X,I3,0P,F10.5,F10.1,2G10.3,1P,5E9.2,0P,F8.3,1P, * 2 (1X,E9.2)) 101 CONTINUE C CALL CPUTIME('WATMOS',0,1,1) RETURN END C C************************************************************************ C FUNCTION DIFF1(K,Y,X,N) C C FIRST DIFFERENTIAL OF ARRAY Y WRT X AT INDEX K C INCLUDE 'PREC' INTEGER K,N DIMENSION X(N),Y(N) DIFF1 = 0.00 IF (N.EQ.1) RETURN IF (K.EQ.1) THEN DIFF1 = (Y(2)-Y(1))/ (X(2)-X(1)) ELSE IF (K.EQ.N) THEN DIFF1 = (Y(N)-Y(N-1))/ (X(N)-X(N-1)) ELSE DIFF1 = 0.50* ((Y(K)-Y(K-1))/ (X(K)-X(K-1))+ * (Y(K+1)-Y(K))/ (X(K+1)-X(K))) END IF RETURN END SUBROUTINE WATOM(ICALL) C C WRITES ATOMIC PARAMETERS TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CSLINE' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWATOM.EQ.0) RETURN C CALL CPUTIME('WATOM ',0,0,1) WRITE(LOUT,50) ATOMID 50 FORMAT('1 ATOMIC PARAMETERS ',A//) WRITE(LOUT,100) ABND, 10.**(ABND-12.) 100 FORMAT(' ABUNDANCE=',F7.4,1P,E11.4//7X,'E(CM-1)', * 6X,'E(EV)',4X,'G',24X,'ION'/) WRITE(LOUT,200) (I,EV(I)/CC/HH*EE,EV(I),G(I),LABEL(I), * ION(I),I=1,NK) 200 FORMAT(1X,I3,F11.3,F11.6,F5.0,2X,A,I3) C QN=QNORM*1.E8/CC WRITE(LOUT,300) 300 FORMAT(/' RADIATIVE TRANSITIONS') WRITE(LOUT,400) (JRAD(KR),IRAD(KR),CONVL(ALAMB(KR)),F(KR),A(KR), * NQ(KR),QMAX(KR),Q0(KR),QN*ALAMB(KR),IWIDE(KR),KR=1,NLINE) 400 FORMAT(' J I LAMBDA',8X,'F',10X,'A',6X,'NQ',6X,'QMAX', * 8X,'Q0',3X,'QNORM MA',' IWIDE'/ * (1X,I2,I3,0P,F10.3,1P,2E11.3,I4,0P,2F10.2,F11.3,4X,L1)) DO 550 KR=NLINE+1,NRAD WRITE(LOUT,500) JRAD(KR),IRAD(KR),CONVL(ALAMB(KR)),F(KR), * NQ(KR),QMAX(KR),Q0(KR),QN*ALAMB(KR) 500 FORMAT(/1X,I2,I3,0P,F10.3,1P,E11.3,11X,I4,0P,3F10.2) IF(Q0(KR).LT.0.0 .OR. QMAX(KR).LT.0.0) THEN WRITE(LOUT,540) (Q(NY,KR),ALFAC(NY,KR-NLINE),NY=1,NQ(KR)) 540 FORMAT(0P,F16.3,1P,E11.3) ENDIF 550 CONTINUE C IF(ICALL.EQ.1) THEN IF(NRFIX.GT.0) THEN WRITE(LOUT,600) 600 FORMAT(/' FIXED TRANSITIONS'/' J I LAMBDA',7X,'A0',8X, * 'TRAD',5X,'ITRAD') ELSE WRITE(LOUT,610) 610 FORMAT(/' NO FIXED TRANSITIONS') ENDIF C DO 800 KF=1,NRFIX CLAM=HCE/(EV(JFX(KF))-EV(IFX(KF))) IF(IPHO(KF).EQ.0) THEN WRITE(LOUT,700) JFX(KF),IFX(KF),CLAM,A0(KF), * TRAD(KF),ITRAD(KF) 700 FORMAT(2I3,F10.3,1P,E11.3,0P,F11.3,I5,' BOUND-BOUND') ELSE WRITE(LOUT,710) JFX(KF),IFX(KF),CLAM,A0(KF), * TRAD(KF),ITRAD(KF) 710 FORMAT(2I3,F10.3,1P,E11.3,0P,F11.3,I5,' PHOTOIONIZATION') ENDIF 800 CONTINUE ENDIF WRITE(LOUT,900) 900 FORMAT(////1X) CALL CPUTIME('WATOM ',0,1,1) RETURN END C C**************************************************************** C SUBROUTINE WCONT(KR) C C WRITES CONTINUUM INTENSITIES AND FLUXES TO OUTPUT LISTING C C: C: WCONT 88-11-01 MODIFICATION: (PHILIP JUDGE) C: ALSO PRINTS RADIATION TEMPERATURES AND C: OUTPUTS FULL CHARACTER STRING LABEL C: C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' C PGJ ADDITION IF(IWLINE.EQ.0) RETURN C KT=KTRANS(KR) I=IRAD(KR) J=JRAD(KR) WRITE(LOUT,100) J,I,LABEL(J),LABEL(I) 100 FORMAT('1 CONTINUUM TRANSITION',10X,I2,'-',I2, * ' (',A,'-',A,')') C C PRINT FLUXES AND INTENSITIES C C PGJ ADDITION: WRITE(LOUT,400) 400 FORMAT(//' PHYSICAL FLUXES IN ERGS/CM**2/S/HZ LAMBDA', * ' IN ANGSTROM'/' INTENSITIES IN ERGS/CM**2/S/HZ/STERADIAN'/ * ' EQUIVALENT RADIATION (BLACK BODY) TEMPERATURES IN K') DO 600 M=NQ(KR),1,-12 NYMIN=MAX(1,M-11) WRITE(LOUT,500) (CC/FRQ(NY,KT)*1.E8,NY=M,NYMIN,-1) WRITE(LOUT,510) (FLUX(NY,KR),NY=M,NYMIN,-1) C PGJ ADDITION WRITE(LOUT,540) (PGJRT(FLUX(NY,KR),CC/FRQ(NY,KT)*1.E8) * ,NY=M,NYMIN,-1) WRITE(LOUT,520) DO 450 MU=1,NMU WRITE(LOUT,530) XMU(MU),(OUTINT(NY,MU,KR),NY=M,NYMIN,-1) 450 CONTINUE 500 FORMAT(//' LAMBDA ',1P,12E9.2) 510 FORMAT(' FLUX ',1P,12E9.2) 520 FORMAT(' MY/ I ') 530 FORMAT(1X,0P,F8.5,1P,12E9.2) C PGJ ADDITION: 540 FORMAT(' RAD T ',0P,12F9.1) 600 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE WDAMP(K,KR,DDP,GAMMA,GR,GV,GS) C C WRITES DAMPING PARAMETERS TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWDAMP.EQ.0) RETURN C IF(K.EQ.0) THEN WRITE(LOUT,100) 100 FORMAT('1 DOPPLER WIDTH, ADAMP AND CONTRIBUTIONS TO ADAMP', * ' IN PERCENT'//3X,'K',2X,'DNYD',3X,'J I',3X,'DNYD',4X,'ADAMP', * 5X,'RAD',5X,'VW',5X,'STARK'/6X,'KM/S',10X,'MA') ELSE IF(MOD(K-1,IWDAMP).EQ.0) THEN J=JRAD(KR) I=IRAD(KR) IF(GAMMA.GT.0.0) THEN GAM=GAMMA ELSE GAM=1.0 ENDIF IF(KR.EQ.1) THEN WRITE(LOUT,200) K,DNYD(K)*QNORM,J,I,DDP,ADAMP(K,KR),GR/GAM, * GV/GAM,GS/GAM 200 FORMAT(1X,I3,F7.3,2I3,F8.3,1P,E10.3,2P,3F8.3) ELSE WRITE(LOUT,300) J,I,DDP,ADAMP(K,KR),GR/GAM,GV/GAM, * GS/GAM 300 FORMAT(11X,2I3,F8.3,1P,E10.3,2P,3F8.3) ENDIF ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE WEMAX(TEXT,VALUE) C C WRITES MAX CHANGE FROM PREVIOUS ITERATION TO OUTPUT LISTING C C: C: WEMAX 86-06-28 MODIFICATIONS: (MATS CARLSSON) C: WRITES BOTH TO OUTPUT FILE AND TO TERMINAL C: 89-03-18 MODIFICATIONS: (MATS CARLSSON) C: USES ONE MORE DIGIT IN E FORMAT TO MAKE ROOM FOR SIGN C: TAKES THE LOG OF THE ABSOLUTE VALUE C: INCLUDE 'PREC' INCLUDE 'CINPUT' INCLUDE 'CLU' C CHARACTER*(*) TEXT C IF(IWEMAX.EQ.0) RETURN C IF(TEXT.EQ.' ') THEN WRITE(LOUT,100) WRITE(*,100) 100 FORMAT(1X) ELSE IF(VALUE.NE.0.0) THEN VALLG=LOG10(ABS(VALUE)) ELSE VALLG=-99. ENDIF WRITE(LOUT,200) IT, TEXT,VALLG,TEXT,VALUE WRITE(*,200) IT, TEXT,VALLG,TEXT,VALUE 200 FORMAT(I4,1X,'LOG(',A,')=',F10.4,2X,A,'=',1P,E11.4) ENDIF RETURN END C C*********************************************************************** C SUBROUTINE WEQW C C WRITES THE EQUIVALENT WIDTHS C C: C: WEQW 88-01-19 MODIFICATIONS: (MATS CARLSSON) C: WEQ/WEQLTE IS SET TO 0.0 IF WEQLTE=0.0 C: 86-10-01 MODIFICATIONS: (PHILIP JUDGE) C: OUTPUT INTEGRATED FLUXES C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWEQW.EQ.0) RETURN C WRITE(LOUT,100) 100 FORMAT ('1 EQUIVALENT WIDTHS',/,/,3X,'KR',5X,'LAMBDA',5X,'J',4X, * 'I',7X,'W(MA)',5X,'W/W*',5X,' INTEGRATED FLUX (ERG/CM2/S)', * /) DO 300 KR=1,NLINE IF(WEQLTE(KR).EQ.0.0) THEN BWEQ=0.0 ELSE BWEQ=WEQ(KR)/WEQLTE(KR) ENDIF C PGJ ADDED INTEGRATED FLUX WRITE (LOUT,FMT=200) KR,CONVL(ALAMB(KR)),JRAD(KR),IRAD(KR), * WEQ(KR)*1.E3,BWEQ, * -WEQ(KR)*FLUX(0,KR)*2.9979246E18/ALAMB(KR)/ALAMB(KR) 200 FORMAT (I5,F12.3,2I5,F12.3,F10.3,4X,1P,E9.2) 300 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE WHSE(OLDNE) C C WRITES RELATIVE CHANGE IN THE ELECTRON DENSITIES C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION OLDNE(MDEP) C IF(IWHSE.EQ.0) RETURN C WRITE(LOUT,100) (J,J=1,10) 100 FORMAT(/' HYDROSTATIC EQUILIBRIUM INTEGRATION'/ * ' RELATIVE CHANGE IN ELECTRON DENSITIES IN PERCENT'/ * 5X,10I7) DO 300 M=1,NDEP,10 KMAX=MIN(NDEP,M+9) WRITE(LOUT,200) M-1,(NE(K)/OLDNE(K)-1.0,K=M,KMAX) 200 FORMAT(1X,I3,2X,10(1X,2P,F6.0)) 300 CONTINUE C RETURN END C C************************************************************************ C SUBROUTINE WINPUT C C WRITES INPUT PARAMETERS C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' INCLUDE 'CLU' C WRITE(LOUT,100) 'DIFF ',DIFF, 'ELIM1 ',ELIM1, * 'ELIM2 ',ELIM2, 'QNORM ',QNORM, 'THIN ',THIN WRITE(LOUT,200) 'IATOM2',IATOM2, 'ICONV ',ICONV, * 'IHSE ',IHSE, 'ILAMBD',ILAMBD, 'IOPAC ',IOPAC, * 'IOPACL',IOPACL, 'ISTART',ISTART, 'ISUM ',ISUM, * 'ITMAX ',ITMAX, 'ITRAN ',ITRAN, 'ISCAT ',ISCAT, * 'INCRAD',INCRAD, 'INGACC',INGACC, 'ICRSW ',ICRSW, * 'NMU ',NMU, * 'IWABND',IWABND, 'IWATMS',IWATMS, 'IWATOM',IWATOM WRITE(LOUT,200) * 'IWCHAN',IWCHAN, 'IWDAMP',IWDAMP, 'IWEMAX',IWEMAX, * 'IWEQW ',IWEQW, 'IWEVEC',IWEVEC, 'IWHEAD',IWHEAD, * 'IWHSE ',IWHSE, 'IWLGMX',IWLGMX, 'IWLINE',IWLINE, * 'IWLTE ',IWLTE, 'IWN ',IWN, 'IWNIIT',IWNIIT, * 'IWOPAC',IWOPAC, 'IWRAD ',IWRAD, 'IWRATE',IWRATE WRITE(LOUT,200) * 'IWSTRT',IWSTRT, 'IWTAUQ',IWTAUQ, 'IWTEST',IWTEST, * 'IWWMAT',IWWMAT, 'IWJFIX',IWJFIX, 'IWARN ',IWARN, * 'IDL1 ',IDL1, 'IDLNY ',IDLNY, 'IDLCNT',IDLCNT 100 FORMAT(1X,A,5X,1P,E15.7) 200 FORMAT(1X,A,5X,I7) C RETURN END C C*********************************************************************** C SUBROUTINE WLGMX C C WRITES MAX DELTA LG TAUQ, S, PMS TO OUTPUT LISTING C: C: WLGMX 88-01-26 MODIFICATIONS: (MATS CARLSSON) C: WRITES MAX AND MIN TAUQ INSTEAD OF MAX DELTA LG S, PMS C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLGMX' INCLUDE 'CLU' C IF(IWLGMX.EQ.0) RETURN C WRITE(LOUT,100) 100 FORMAT('1 D LG TAUNY IS MAXIMUM LG(TAUNY(K))-LG(TAUNY(K-1))'/ * ' K LG DEPTH D LG TAUNY MAX TAUNY MIN TAUNY') IF(DPTYPE.EQ.'M') THEN WRITE(LOUT,200) (K,LOG10(CMASS(K)),DLGTMX(K),(TAUMAX(K)), * (TAUMIN(K)),K=1,NDEP,IWLGMX) ELSE WRITE(LOUT,200) (K,LOG10(TAU(K)),DLGTMX(K),(TAUMAX(K)), * (TAUMIN(K)),K=1,NDEP,IWLGMX) ENDIF 200 FORMAT(1X,I3,0P,2F11.6,1P,2E11.3) RETURN END C C*********************************************************************** C SUBROUTINE WLINE(KR) C C WRITES LINE INTENSITIES AND FLUXES TO OUTPUT LISTING C C: WLINE 90-01-10 MODIFICATIONS: (MATS CARLSSON) C: WRITE -99 INSTEAD OF LOG(X) IF X.LE.0 C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWLINE.EQ.0) RETURN C TAUQ(0)=TAUQ(1) TAUC0=TAU(1)*XCONT(1) TAUC1=TAUC0 J=JRAD(KR) I=IRAD(KR) WRITE(LOUT,100) CONVL(ALAMB(KR)),J,I,LABEL(J),LABEL(I), * Q(NQ(KR),KR) 100 FORMAT('1 LINE CENTER, LAMBDA=',F10.3,10X,I2,'-',I2,' (', * A,'-',A,')',19X,'LINE WING,NY=', * F10.3//3X,'K',1X,'LG TAU0',3X,'LG TAUNY',1X,'D LG TAUNY', * 1X,'XLINE',7X,'LG TAUC',3X,'D LG TAUC XCONT',9X,'XCONT/XTOT', * 3X,'COOLING') C C PRINT TAU-SCALES AND OPACITIES C DO 300 K=1,NDEP IF(MOD(K-1,IABS(IWLINE)).EQ.0) THEN WRITE(LOUT,200) K,ALOGP(TAU(K)),ALOGP(TAUQ(K)), * ALOGP(TAUQ(K))-ALOGP(TAUQ(K-1)),X(K)-XCONT(K), * ALOGP(TAUC1),ALOGP(TAUC1)-ALOGP(TAUC0),XCONT(K), * RNY(K),COOL(K,KR) 200 FORMAT(1X,I3,0P,3F10.6,1P,E10.3,2X,0P,2F10.6,1P,E10.3, * 2(4X,1P,E10.3)) ENDIF TAUC0=TAUC1 IF(K.LT.NDEP) TAUC1=TAUC1+0.5*(XCONT(K+1)+XCONT(K))* * (TAU(K+1)-TAU(K)) 300 CONTINUE C C PRINT FLUXES AND INTENSITIES C WRITE(LOUT,310) 310 FORMAT(//' FREQUENCY QUADRATURE') DO 330 M=1,NQ(KR),10 NYMAX=MIN(NQ(KR),M+9) WRITE(LOUT,320) (Q(NY,KR),NY=M,NYMAX) WRITE(LOUT,325) (DLAMB(Q(NY,KR),KR),NY=M,NYMAX) 320 FORMAT(/8X,'Q',1P,10E11.3) 325 FORMAT(' ANGSTROM',1P,10E11.3) 330 CONTINUE WRITE(LOUT,400) 400 FORMAT(//' REL FLUX RELATIVE TO CONTINUUM FLUX, DELTA LAMBDA', * ' IN ANGSTROM'/' REL INT. RELATIVE TO CONTINUUM INTENSITY') DO 600 M=NQ(KR),1,-12 NYMIN=MAX(1,M-11) WRITE(LOUT,500) (DLAMB(Q(NY,KR),KR),NY=M,NYMIN,-1) WRITE(LOUT,505) (FLUX(NY,KR),NY=M,NYMIN,-1) WRITE(LOUT,510) (FLUX(NY,KR)/FLUX(0,KR),NY=M,NYMIN,-1) WRITE(LOUT,515) DO 440 MU=1,NMU WRITE(LOUT,530) XMU(MU),(OUTINT(NY,MU,KR),NY=M,NYMIN,-1) 440 CONTINUE WRITE(LOUT,520) DO 450 MU=1,NMU WRITE(LOUT,530) XMU(MU),(OUTINT(NY,MU,KR)/OUTINT(0,MU,KR), * NY=M,NYMIN,-1) 450 CONTINUE 500 FORMAT(//' D LAMBDA',1P,12E9.2) 505 FORMAT(' FLUX ',1P,12E9.2) 510 FORMAT(' REL FLUX',1P,12E9.2) 515 FORMAT(' MY/ I') 520 FORMAT(' MY/REL I') 530 FORMAT(1X,0P,F8.5,1P,12E9.2) 600 CONTINUE WRITE(LOUT,610) FLUX(0,KR) 610 FORMAT(//' CONTINUUM FLUX=',1P,E11.4) WRITE(LOUT,620) (XMU(MU),OUTINT(0,MU,KR),MU=1,NMU) 620 FORMAT(' MY',6X,' CONT. INT.'/(1X,0P,F8.5,1P,E11.4)) C RETURN END C C*********************************************************************** C FUNCTION ALOGP(X) C C RETURNS LOG10(X) IF X IS POSITIVE -99 IF X IS NEGATIVE OR ZERO C INCLUDE 'PREC' C IF(X.GT.0.) THEN ALOGP=LOG10(X) ELSE ALOGP=-99. ENDIF C RETURN END SUBROUTINE WLTE C C WRITES LTE POPULATIONS TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWLTE.EQ.0) RETURN C CALL CPUTIME('WLTE ',0,0,1) WRITE(LOUT,50) 50 FORMAT('1 LTE POPULATIONS') DO 400 M=1,NK,12 IMAX=MIN(NK,M+11) WRITE(LOUT,100) (I,I=M,IMAX) 100 FORMAT(//2X,'K',4X,12(I3,7X)/) DO 300 K=1,NDEP,IWLTE WRITE(LOUT,200) K,(NSTAR(I,K),I=M,IMAX) 200 FORMAT(1X,I3,1P,12E10.3) 300 CONTINUE 400 CONTINUE CALL CPUTIME('WLTE ',0,1,1) RETURN END C C*********************************************************************** C SUBROUTINE WNIIT(INIT) C C WRITES POPULATIONS TO FILE NIIT UNFORMATTED. C INIT = 0 BOLTZMANN POPULATIONS C INIT.GT.0 POPULATIONS C: C: WNIIT 89-03-19 MODIFICATIONS: (MATS CARLSSON) C: WRITES UNFORMATTED FILE IN SINGLE PRECISION C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWNIIT.EQ.0) RETURN C CALL CPUTIME('WNIIT ',0,0,1) IF(INIT.EQ.0) THEN WRITE(LNIIT) ((REAL(NSTAR(I,K)),I=1,NK),K=1,NDEP) ELSE WRITE(LNIIT) ((REAL(N(I,K)),I=1,NK),K=1,NDEP) ENDIF CALL CPUTIME('WNIIT ',0,1,1) RETURN END C C*********************************************************************** C SUBROUTINE WN C C WRITES POPULATION NUMBERS AND IONIZATION FRACTIONS TO OUTPUT LISTING C C: C: WN 88-02-04 MODIFICATIONS: (PHILIP JUDGE, MATS CARLSSON) C: POPULATION DENSITIES AND DEPARTURE COEFFICIENTS WRITTEN C: SEPARATELY. DEPARTURE COEFFICIENTS DEFINED AS N/NSTAR C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION ANION(5),ANIONS(5) CHARACTER*5 SYMB(10) DATA SYMB/'I','II','III','IV','V','VI','VII','VIII','IX','X'/ C IF(IWN.EQ.0) RETURN C C MODIFIED BY PGJ TO OUTPUT N'S AND B'S SEPARATELY C SMALL=1.E-10 DO 400 M=1,NK,10 IMAX=MIN(NK,M+9) WRITE(LOUT,100) (' ',I,I=M,IMAX) 100 FORMAT('1 SOLUTION- NUMBER DENSITIES'// * ' K',10(4X,A,'N(',I2,')')) DO 300 K=1,NDEP,IWN WRITE(LOUT,200) K,(N(I,K), I=M,IMAX) 200 FORMAT(1X,I3,10(1P,E10.3,1X)) 300 CONTINUE 400 CONTINUE C C DEPARTURE COEFFS. C DO 401 M=1,NK,10 IMAX=MIN(NK,M+9) WRITE(LOUT,101) (' ',I,I=M,IMAX) 101 FORMAT('1 SOLUTION- DEPARTURE COEFFICIENTS', * ' B(I)=N(I)/NSTAR(I)'// * ' K',10(4X,A,'B(',I2,')')) DO 301 K=1,NDEP,IWN WRITE(LOUT,201) K,(N(I,K)/NSTAR(I,K),I=M,IMAX) 201 FORMAT(1X,I3,10(1P,E10.3,1X)) 301 CONTINUE 401 CONTINUE C C FIND LOWEST AND HIGHEST IONIZATION STAGE C MINION=100 MAXION=0 DO 450 I=1,NK MINION=MIN(MINION,ION(I)) MAXION=MAX(MAXION,ION(I)) 450 CONTINUE IF(MAXION.GT.10) THEN DO 480 I=MINION,MAXION J=I-MINION+1 WRITE(SYMB(J),460) I,' ' 460 FORMAT(I2,A) 480 CONTINUE ENDIF C WRITE(LOUT,500) 500 FORMAT('1 IONIZATION FRACTIONS IN PERCENT AND DEVIATIONS FROM', * ' LTE'//) WRITE(LOUT,510) (SYMB(I),I=MINION,MAXION) 510 FORMAT(5(17X,A)) WRITE(LOUT,520) (' ',I=MINION,MAXION) 520 FORMAT(4X,5(A,7X,'PERC N/N*')/) DO 900 K=1,NDEP,IWN DO 600 I=MINION,MAXION ANION(I)=0.0 ANIONS(I)=0.0 600 CONTINUE DO 700 J=1,NK I=ION(J) ANION(I)=ANION(I)+N(J,K) ANIONS(I)=ANIONS(I)+NSTAR(J,K) 700 CONTINUE DO 750 I=MINION,MAXION ANIONS(I)=MAX(ANIONS(I),SMALL) ANIONS(I)=MAX(ANIONS(I),ANION(I)/99999.0) 750 CONTINUE WRITE(LOUT,800) K,(ANION(I)/TOTN(K),ANION(I)/ANIONS(I), * I=MINION,MAXION) 800 FORMAT(1X,I3,5(2P,F12.3,0P,F10.3)) 900 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE WOPAC(K,KR,WAVE) C C WRITES BACKGROUND OPACITIES TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMO2' INCLUDE 'CINPUT' INCLUDE 'CLU' C COMMON/ABSBID/ OPH,OPHE,OPM,OPSC,CHI,ABSLTE,HYDB(6),HELB(6), * AMTB(16),SCATB(4),ABH(5),ABHLTE(5) C IF(IWOPAC.EQ.0) RETURN C IF(KR.EQ.0 .AND. K.EQ.0) THEN WRITE(LOUT,100) 100 FORMAT('1 O P A C I T I E S'/24X,'PERCENTAGES OF ', * 'CONTINUUS OPACITY'/' K LAMBDA OPACITY',3X, * ' HBF HFF H-BF H-FF H2- H2+ HE SI MG AL', * ' FE C ELECT R(H) R(H2) R(HE)',3X, * 'ABS SCAT') ELSE IF(MOD(K-1,IWOPAC).EQ.0) THEN OP=CHI/RHO(K) IF(KR.EQ.0) THEN WRITE(LOUT,200) K,5000.,OP,(HYDB(I)/CHI,I=1,6),OPHE/CHI, * (AMTB(I)/CHI,I=1,5),(SCATB(I)/CHI,I=1,4), * 1.-OPSC/CHI,OPSC/CHI 200 FORMAT(/1X,I3,1X,F8.0,1P,E10.3,2X,2P,12F5.1,3X,4F5.1, * 3X,2F6.1) ELSE WRITE(LOUT,300) CONVL(WAVE),OP,(HYDB(I)/CHI,I=1,6),OPHE/CHI, * (AMTB(I)/CHI,I=1,5),(SCATB(I)/CHI,I=1,4), * 1.-OPSC/CHI,OPSC/CHI 300 FORMAT(5X,F8.0,1P,E10.3,2X,2P,12F5.1,3X,4F5.1, * 3X,2F6.1) ENDIF ENDIF ENDIF RETURN END C C**************************************************************** C SUBROUTINE WRATE C C WRITES RADIATIVE, COLLISIONAL AND NET RATES TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION R(MK,MK),P(MK,MK) CHARACTER CH1(MK,MK) C IF(IWRATE.EQ.0) RETURN SMALL=1.E-33 C WRITE(LOUT,100) 100 FORMAT('1 L O G A R I T H M I C R A T E S '// * ' R= RADIATIVE RATES PER ATOM AND SEC'/ * ' C= COLLISIONAL RATES PER ATOM AND SEC'/ * ' T= UPPER HALF: TOTAL NET RATE PER CM3 AND SEC'/ * ' I.E. ABS(UP-DOWN)'/ * ' LOWER HALF: RATIO OF UP TO DOWN RATES PER'/ * ' CM3'/ * ' DIAGONAL ELEMENTS CONTAIN POPULATIONS'/ * ' ZERO INDICATE ZERO RATE'//) C CH1(NK,NK)=' ' DO 800 K=1,NDEP,IWRATE READ(LDUMC,REC=K) ((P(I,J),I=1,NK),J=1,NK) DO 300 I=1,NK-1 CH1(I,I)=' ' DO 200 J=I+1,NK R(I,J)=SMALL R(J,I)=SMALL KR=KRAD(I,J) IF(KR.NE.0) THEN R(I,J)=RIJ(K,KR) R(J,I)=RJI(K,KR) END IF UP=N(I,K)*(R(I,J)+C(I,J,K)) DOWN=N(J,K)*(R(J,I)+C(J,I,K)) R(I,J)=R(I,J)+C(I,J,K)-P(I,J) R(J,I)=R(J,I)+C(J,I,K)-P(J,I) C(I,J,K)=P(I,J) C(J,I,K)=P(J,I) P(I,J)=ABS(UP-DOWN) IF(DOWN.NE.0.0) THEN P(J,I)=UP/DOWN ELSE P(J,I)=0.0 ENDIF IF(UP.GT.DOWN) THEN CH1(I,J)='U' ELSE CH1(I,J)='D' END IF CH1(J,I)=' ' IF(R(I,J).LE.SMALL) THEN R(I,J)=1. R(J,I)=1. END IF R(J,I)=MAX(R(J,I),SMALL) P(I,J)=MAX(P(I,J),SMALL) P(J,I)=MAX(P(J,I),SMALL) C(I,J,K)=MAX(C(I,J,K),SMALL) C(J,I,K)=MAX(C(J,I,K),SMALL) 200 CONTINUE R(I,I)=MAX(R(I,I),SMALL) P(I,I)=MAX(P(I,I),SMALL) C(I,I,K)=MAX(C(I,I,K),SMALL) 300 CONTINUE C WRITE(LOUT,400) K,TEMP(K),NE(K) 400 FORMAT(///' K=',I3,' T=',F7.0,' NE=',1P,E12.3) DO 700 JMIN=1,NK,15 JMAX=MIN(JMIN+14,NK) WRITE(LOUT,500) (J,J=JMIN,JMAX) 500 FORMAT(/' FROM/TO',I4,14I7) DO 600 I=1,NK WRITE(LOUT,510) 510 FORMAT(/) IF(I.LT.JMIN.OR.I.GT.JMAX) THEN WRITE(LOUT,520) I,(LOG10(R(I,J)),J=JMIN,JMAX) WRITE(LOUT,530) (LOG10(C(I,J,K)),J=JMIN,JMAX) WRITE(LOUT,540) (LOG10(P(I,J)),CH1(I,J),J=JMIN,JMAX) ELSE WRITE(LOUT,520) I,(LOG10(R(I,J)),J=JMIN,I-1), * LOG10(N(I,K)),(LOG10(R(I,J)),J=I+1,JMAX) WRITE(LOUT,530) (LOG10(C(I,J,K)),J=JMIN,I-1), * LOG10(N(I,K)),(LOG10(C(I,J,K)),J=I+1,JMAX) WRITE(LOUT,540) (LOG10(P(I,J)),CH1(I,J),J=JMIN,I-1), * LOG10(N(I,K)),CH1(I,I),(LOG10(P(I,J)),CH1(I,J), * J=I+1,JMAX) 520 FORMAT(I3,' R: ',15F7.2) 530 FORMAT(4X,'C. ',15F7.2) 540 FORMAT(4X,'T: ',1X,15(F6.2,A1)) END IF 600 CONTINUE 700 CONTINUE 800 CONTINUE C RETURN END C C********************************************************************* C SUBROUTINE WRSTRT C C WRITES THE SOLUTION TO FILE RSTRT2 C TO BE USED AS INITIAL SOLUTION FOR RESTART FROM PARTIALLY CONVERGED C SOLUTION C SAME FORMAT AS ATMOS C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CATMO2' INCLUDE 'CLU' C CALL OPEN(LRSTRT,'RSTRT2',1,'NEW') WRITE(LRSTRT,40) 40 FORMAT('* RESTART FILE') WRITE(LRSTRT,50) ATMOID 50 FORMAT(A) IF(DPTYPE.EQ.'M') THEN WRITE(LRSTRT,70) 70 FORMAT(' MASS SCALE') ELSE WRITE(LRSTRT,80) 80 FORMAT(' TAU5000 SCALE') ENDIF WRITE(LRSTRT,90) LOG10(GRAV) 90 FORMAT(F10.7) WRITE(LRSTRT,100) NDEP 100 FORMAT(1X,I3) IF(DPTYPE.EQ.'M') THEN WRITE(LRSTRT,200) (LOG10(CMASS(K)),TEMP(K),NE(K),VEL(K)*QNORM, * VTURB(K)*QNORM,K=1,NDEP) ELSE WRITE(LRSTRT,200) (LOG10(TAU(K)),TEMP(K),NE(K),VEL(K)*QNORM, * VTURB(K)*QNORM,K=1,NDEP) ENDIF 200 FORMAT(1P,5E15.6) DO 400 K=1,NDEP WRITE(LRSTRT,300) (N(I,K),I=1,NK) 300 FORMAT(1P,6E12.4) 400 CONTINUE CALL CLOSE(LRSTRT) C RETURN END C C******************************************************************** C SUBROUTINE WSTART C C WRITES STARTING SOLUTION TO OUTPUT LISTING C C: C: WSTART 88-02-04 MODIFICATIONS: (PHILIP JUDGE, MATS CARLSSON) C: POPULATION DENSITIES AND DEPARTURE COEFFICIENTS WRITTEN C: SEPARATELY. DEPARTURE COEFFICIENTS DEFINED AS N/NSTAR C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C IF(IWSTRT.EQ.0) RETURN C C ALTERED BY PGJ TO OUTPUT NUMBER DENSITIES & DEP COEFFS. SEPARATELY C DO 400 M=1,NK,10 IMAX=MIN(NK,M+9) WRITE(LOUT,100) (' ',I,I=M,IMAX) 100 FORMAT('1 STARTING APPROXIMATION - NUMBER DENSITIES'// * ' K',10(4X,A,'N(',I2,')')) DO 300 K=1,NDEP,IWN WRITE(LOUT,200) K,(N(I,K), I=M,IMAX) 200 FORMAT(1X,I3,10(1P,E10.3,1X)) 300 CONTINUE 400 CONTINUE C C DEPARTURE COEFFS. C DO 401 M=1,NK,10 IMAX=MIN(NK,M+9) WRITE(LOUT,101) (' ',I,I=M,IMAX) 101 FORMAT('1 STARTING APPROXIMATION- DEPARTURE COEFFICIENTS', * ' B(I)=N(I)/NSTAR(I)'// * ' K',10(4X,A,'B(',I2,')')) DO 301 K=1,NDEP,IWN WRITE(LOUT,201) K,(N(I,K)/NSTAR(I,K),I=M,IMAX) 201 FORMAT(1X,I3,10(1P,E10.3,1X)) 301 CONTINUE 401 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE WTAUQQ(KR) C C WRITES LOG TAUNY SCALES AND DELTA LOG TAUNY TO OUTPUT LISTING C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CTAUQQ' INCLUDE 'CLU' C IF(IWTAUQ.EQ.0) RETURN C I=IRAD(KR) J=JRAD(KR) WRITE(LOUT,100) J,I,CONVL(ALAMB(KR)) 100 FORMAT('1 TRANSITION ',I2,'-',I2,F13.3,' ANGSTROM') IF(KR.GT.NLINE) KT=KTRANS(KR) WRITE(LOUT,110) 110 FORMAT(/' LG TAUNY FOR LAMBDA IN ANGSTROM') C C PRINT LG TAUQQ AND STORE LG(TAUQQ(K)) - LG(TAUQQ(K-1)) IN C TAUQQ(K-1) C DO 300 M=NQ(KR),1,-12 NYMIN=MAX(1,M-11) IF(KR.LE.NLINE) THEN WRITE(LOUT,200) (DLAMB(Q(NY,KR),KR),NY=M,NYMIN,-1) 200 FORMAT(//' K/LAMBDA',1P,12E9.2/) ELSE WRITE(LOUT,201) (CC/FRQ(NY,KT)*1.E8,NY=M,NYMIN,-1) 201 FORMAT(//' K/LAMBDA',0P,12F9.0/) ENDIF DO 250 K=1,NDEP IF(MOD(K-1,IWTAUQ).EQ.0) THEN WRITE(LOUT,220) K,(LOG10(TAUQQ(NY,K)),NY=M,NYMIN,-1) 220 FORMAT(1X,I3,5X,12F9.5) ENDIF IF(K.GT.1) THEN DO 230 NY=M,NYMIN,-1 TAUQQ(NY,K-1)=LOG10(TAUQQ(NY,K))-LOG10(TAUQQ(NY,K-1)) 230 CONTINUE ENDIF 250 CONTINUE 300 CONTINUE C C PRINT D LGTAUNY C WRITE(LOUT,310) 310 FORMAT(///' D LG TAUNY FOR LAMBDA IN ANGSTROM') DO 350 M=NQ(KR),1,-12 NYMIN=MAX(1,M-11) IF(KR.LE.NLINE) THEN WRITE(LOUT,200) (DLAMB(Q(NY,KR),KR),NY=M,NYMIN,-1) ELSE WRITE(LOUT,201) (CC/FRQ(NY,KT)*1.E8,NY=M,NYMIN,-1) ENDIF DO 340 K=2,NDEP,IWTAUQ WRITE(LOUT,330) K,(TAUQQ(NY,K-1),NY=M,NYMIN,-1) 330 FORMAT(1X,I3,5X,12F9.5) 340 CONTINUE 350 CONTINUE C RETURN END C C********************************************************************* C SUBROUTINE COLRAT C C CHOOSES COLLISIONAL ROUTINE C IN FIXRAD FIXED RATES ARE ADDED INTO C. COLLISION PART IS THEREFORE C STORED TO ENABLE SEPARATION IN ROUTINE WRATE C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CCONST' INCLUDE 'CLU' C C INITIALIZE COLLISIONAL RATES TO ZERO C DO 300 K=1,NDEP DO 200 J=1,NK DO 100 I=1,NK C(I,J,K)=0.0 100 CONTINUE 200 CONTINUE 300 CONTINUE C C CHOOSE COLLISIONAL SUBROUTINE C IF(CROUT.EQ.'CA2COL') THEN CALL CA2COL ELSE IF(CROUT.EQ.'HCOL ') THEN CALL HCOL ELSE IF(CROUT.EQ.'COCOL ') THEN CALL COCOL ELSE IF(CROUT.EQ.'KCOL ') THEN CALL KCOL ELSE IF(CROUT.EQ.'GENCOL') THEN CALL GENCOL ELSE CALL STOP('COLRAT: COLLISIONAL ROUTINE UNKNOWN') ENDIF C C OUTPUT TO DISK C CALL CPUTIME('CDISK ',0,0,1) CALL OPEN(LDUMC,'DUMC',NK*NK,'UNKNOWN') DO 600 K=1,NDEP WRITE(LDUMC,REC=K) ((C(I,J,K),I=1,NK),J=1,NK) 600 CONTINUE CALL CPUTIME('CDISK ',0,1,1) RETURN END C C ********************************************************************** C SUBROUTINE KCOL C C COMPUTES COLLISIONAL RATES FOR ALKALI ELEMENTS SODIUM AND POTASSIUM C THE ROUTINE IS WRITTEN FOR BOUND LEVELS PLUS ONE CONTINUUM LEVEL C C OPTION TO CALCULATE C BOUND-FREE C OPTION -2.0 : ALLEN, 1973, ASTROPHYSICAL QUANTITIES, P.42 (NEUTRAL ATOM C APPROXIMATION, THIRD EQUATION OF P. 42; FOR EXPLANATION C OF SYMBOLS REFER TO PAGE 41). NO ADDITIONAL INPUT C DATA REQUIRED FOR THIS FORMULA. C OPTION -3.0 : INTERPOLATION OF TABLES OF KUNC AND ZGORZELSKI, ATOMIC C DATA AND NUCLEAR DATA TABLES, 19, 1-21, 1971. C FIRST INPUT ITEM AFTER OPTION NUMBER IS THE NUMBER OF C TEMPERATURE POINTS. AFTER THAT ALL THE TEMPERATURES ARE READ. C THEY HAVE TO BE IN ASCENDING ORDER, OTHERWISE A FATAL C ERROR IS ISSUED. C FINALLY THE COLLISIONAL RATES DATA FOR ALL LEVELS ARE READ C IN ORDER OF ASCENDING TEMPERATURE. C C KCOL: 90-01-01 NEW ROUTINE: (JO BRULS) C THE FORMULAE AND DATA USED ARE FROM PARK, JQSRT, 11, 7-36 C (1971). FOR A TRANSITION BETWEEN STATES THAT HAVE THE C SAME EFFECTIVE PRINCIPAL QUANTUM NUMBER THE STANDARD C FORMULA IS USED FOR A CONSTANT COLLISION STRENGTH OF C UNITY (SEE NOTES FROM NATASHA SHCHUKINA). C THIS IS INDICATED IN THE INPUT FILE BY CE = -1.0 C C 91-01-01 MODIFICATIONS: (JO BRULS) C THE OPTION -1 IS FOR CONSTANT COLLISION STRENGTH UNITY. C THE COLLISIONAL IONIZATION CROSS SECTIONS HAVE SEVERAL OPTIONS. C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CCONST' INCLUDE 'CLU' C PARAMETER(MT=20) DIMENSION CE(MK-1),RE(MK-1),QN(MK-1) DIMENSION CI(MK-1,MT),TEM(MT) C C FIRST CALCULATE ALL EFFECTIVE PRINCIPAL QUANTUM NUMBERS C RYD=13.595 DO 50 J=1,NK-1 QN(J)=(RYD/(EV(NK)-EV(J)))**0.5 50 CONTINUE C C CALCULATE C BOUND-BOUND C C USE EQ. 8 OF PARK (1971) TO CALCULATE THE RATES FOR ALL TRANSITIONS C THAT INVOLVE A CHANGE OF (EFFECTIVE) PRINCIPAL QUANTUM NUMBER, C I.E. FOR ALMOST ALL TRANSITIONS, EXCEPT BETWEEN THE LEVELS OF ONE TERM C (INDICATED BY CE = -1.0 AND ARBITRARY RE IN THE INPUT FILE) C DO 300 J=2,NK-1 IMAX=J-1 READ(LDUMS,*) (CE(I),I=1,IMAX) READ(LDUMS,*) (RE(I),I=1,IMAX) DO 200 I=1,J-1 C C CALCULATE THE DEPTH INDEPENDENT PART HERE C IF (CE(I).NE.-1.0) THEN CE(I)=CE(I)*G(J)/((QN(I)**5.0)*(QN(J)**5.0))/ * (((QN(I)**(-2.0))-(QN(J)**(-2.0)))**4) ELSE CE(I)=8.63E-8/G(I) RE(I)=-0.5 ENDIF C DO 100 K=1,NDEP C(I,J,K)=NE(K)*CE(I)*(TEMP(K)/10000.)**RE(I)* * EXP(-EK*(EV(J)-EV(I))/TEMP(K)) C(J,I,K)=C(I,J,K)*NSTAR(I,K)/NSTAR(J,K) 100 CONTINUE 200 CONTINUE 300 CONTINUE C READ(LDUMS,*) OPT IF (OPT.EQ.-2.0) THEN C C SEATON: UNITY COLLLISION STRENGTH. NO FURTHER INPUT REQUIRED. C DO 410 K=1,NDEP DO 400 I=1,NK-1 C(I,NK,K)=1.1E-8*(TEMP(K)**0.5)*((EV(NK)-EV(I))**(-2.0))* * 10.0**(-5040.0*(EV(NK)-EV(I))/TEMP(K))*NE(K) C(NK,I,K)=C(I,NK,K)*NSTAR(I,K)/NSTAR(NK,K) 400 CONTINUE 410 CONTINUE ELSE IF (OPT.EQ.-3.0) THEN C C KUNC AND ZGORZELSKI, ATDNDT TABLES INTERPOLATION. C READ(LDUMS,*) NT IF (NT.GT.MT) CALL STOP('KCOL: NT GREATER THAN MT') READ(LDUMS,*) (TEM(K),K=1,NT) DO 420 K=1,NT READ(LDUMS,*) (CI(L,K),L=1,NK-1) 420 CONTINUE C C DO A LINEAR INTERPOLATION ON CI VS. TEMP AFTER CHECKING THAT THE C INPUT TEMPERATURES ARE IN ASCENDING ORDER. C DO 430 K=2,NT IF (TEM(K).LE.TEM(K-1)) THEN CALL STOP('KCOL: INPUT TEMPERATURES NOT IN ASCENDING ORDER') ENDIF 430 CONTINUE DO 470 K=1,NDEP T=TEMP(K) L=1 440 L=L+1 IF (T.GT.TEM(L) .AND. L.LT.NT) GOTO 440 APOL=(T-TEM(L-1))/(TEM(L)-TEM(L-1)) DO 450 I=1,NK-1 C(I,NK,K)=(CI(I,L-1)+APOL*(CI(I,L)-CI(I,L-1)))*NE(K) C(NK,I,K)=C(I,NK,K)*NSTAR(I,K)/NSTAR(NK,K) 450 CONTINUE 470 CONTINUE ELSE CALL STOP('KCOL: COLLISIONAL IONIZATION OPTION UNKNOWN') ENDIF RETURN END C C ********************************************************************** C SUBROUTINE TRANSP C C SOLVES THE RADIATIVE TRANSFER EQUATION WITH GIVEN SOURCE FUNCTION. C SOLVES FOR P-S WHERE P IS FEAUTRIERS P; I.E. C C P = 0.5*(I(XMU)+I(-XMU)) C C IN THE PRESENCE OF VELOCITY FIELDS: C C P = 0.5*(I(NY,XMU) + I(-NY,-XMU)) C C ITRAN DETERMINES THE MODE OF FORMAL SOLUTION C ITRAN = 0 FEAUTRIER SOLUTION C 1 FEAUTRIER SOLUTION TO CUBIC SPLINE ACCURACY, REF: C KUNASZ, HUMMER, 1974, MNRAS 166,19 C MIHALAS, 1974, APJ SUPPL 28,343 C 2 FEAUTRIER SOLUTION HERMITE, REF: C AUER, 1976, JQSRT 16,931 C C 3 INTEGRAL CUBIC SPLINE METHOD ACCORDING TO SCHARMER C 4 INTEGRAL CUBIC SPLINE METHOD WITH LOCALLY DETERMINED C FIRST DERIVATIVE C C VARIABLES: C C XMU : ANGULAR COSINE (IN) C NDEP : NUMBER OF DEPTH POINTS (IN) C TAU : STANDARD OPTICAL DEPTH SCALE (IN) C X : RATIO OF MONOCHROMATIC TO STANDARD OPACITY (IN) C S : MONOCHROMATIC SOURCE FUNCTION (IN) C ITRAN : DETERMINES MODE OF FORMAL SOLUTION (IN) C C P : MEAN BIDIRECTIONAL INTENSITY (CF. ABOVE) (OUT) C PMS : P-S (OUT) C IPLUS : RADIATION INTENSITY, OUTGOING RAYS (OUT) C IMINUS: RADIATION INTENSITY, INGOING RAYS (OUT) C TAUQ : MONOCHROMATIC OPTICAL DEPTH (OUT) C DTAUQ : DTAUQ(K)=TAUQ(K)-TAUQ(K-1) (OUT) C A1 : 1/(DTAUQ(K+0.5)*DTAUQ(K)) (OUT) C C1 : 1/(DTAUQ(K+0.5)*DTAUQ(K+1)) (OUT) C C: C: TRANSP 86-09-04 MODIFICATIONS: (GORAN SCHARMER, MATS CARLSSON) C: INCLUDES CALL TO CUBIC SPLINE INTEGRAL FORMAL SOLVER C: C: 92-10-08 MODIFICATIONS: (MATS CARLSSON) C: ITRAN=10-14 ALLOWED (ITRAN.GT.10 INCLUDES INCOMING RADIATION) C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION SWITCHED ON BY INCRAD.NE.0 INSTEAD OF C: ITRAN=10-14 C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CINPUT' C C K=1: UPPER BOUNDARY C CMU=0.5/XMY DTAUQ(2)=(X(1)+X(2))*(TAU(2)-TAU(1))*CMU A1(1)=1./DTAUQ(2) T=TAU(1)*X(1)*2.0*CMU TAUQ(1)=T DTAUQ(1)=T C C CALCULATE DTAUQ C IF(DPTYPE.EQ.'H') THEN DO 100 K=2,NDEP DTAUQ(K)=(X(K)*XNORM(K)+X(K-1)*XNORM(K-1))* * (HEIGHT(K-1)-HEIGHT(K))*1.E5*CMU 100 CONTINUE DO 102 K=2,3 TAUQ(K)=TAUQ(K-1)+DTAUQ(K) 102 CONTINUE TAUQ(1)=EXP(2*LOG(TAUQ(2))-LOG(TAUQ(3))) ELSE DO 105 K=2,NDEP DTAUQ(K)=(X(K)+X(K-1))*(TAU(K)-TAU(K-1))*CMU 105 CONTINUE ENDIF C C CALCULATE TAUQ C DO 110 K=2,NDEP TAUQ(K)=TAUQ(K-1)+DTAUQ(K) 110 CONTINUE C C CALCULATE A1 AND C1 C DO 120 K=2,NDEP-1 A1(K)=2./(DTAUQ(K)+DTAUQ(K+1))/DTAUQ(K) C1(K)=2./(DTAUQ(K)+DTAUQ(K+1))/DTAUQ(K+1) 120 CONTINUE C C CHOOSE THE METHOD OF FORMAL SOLUTION C IF(ISCAT.EQ.1) THEN CALL TRANSC ELSE IF(ITRAN.LE.2) THEN CALL TRANF ELSE CALL TRANI ENDIF C C CALCULATE P(K)-S(K) C PMS(1)=P(1)-S(1) DO 130 K=2,NDEP-1 IF(ABS(A1(K)).GT.1.0) THEN PMS(K)=P(K)-S(K) ELSE PMS(K)=C1(K)*(P(K+1)-P(K))-A1(K)*(P(K)-P(K-1)) END IF 130 CONTINUE PMS(NDEP)=(P(NDEP-1)-P(NDEP)+S(NDEP)-S(NDEP-1)) * /(DTAUQ(NDEP)+0.5*DTAUQ(NDEP)**2) C RETURN END C C********************************************************************** C SUBROUTINE TRANF C C FORMAL SOLVER FOR THE TRANSFER EQUATION USING FEAUTRIER TECHNIQUE C C ITRAN DETERMINES THE MODE OF FORMAL SOLUTION C = 0 FEAUTRIER SOLUTION C 1 FEAUTRIER SOLUTION TO CUBIC SPLINE ACCURACY, REF: C AUER, 1976, JQSRT 16,931 C 2 FEAUTRIER SOLUTION HERMITE C C THE UPPER AND LOWER BOUNDARY CONDITIONS ARE DERIVED FROM SECOND C ORDER TAYLOR EXPANSIONS, WITH ESTIMATES OF THE INCIDENT RADIATION C FIELDS: C C P(2) = P(1) + DTAU*P'(1) + 0.5*DTAU**2*P''(1) C P'(1) = P(1) - I(-XMU) C I(-XMU) = S(1)*(1.-EXP(-TAUQ(1))) C C P(N-1) = P(N) - DTAU*P'(N) + 0.5*DTAU**2*P''(N) C P'(N) = I(XMU) - P(N) C I(XMU) = S(N) + (S(N)-S(N-1))/DTAU C C STEINS TRICK C (STORING THE SUM OF ELEMENTS INSTEAD OF THE DIAGONAL) IS USED TO C AVOID NUMERICAL ROUND-OFF PROBLEMS AT SMALL OPTICAL DEPTHS. C C CODED BY: A.NORDLUND (OCT-1981). REVISED (MAR-1982). C C: TRANF 92-10-08 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION FIELD IS SPECIFIED IN IMINUS(0) C: THIS IS INDICATED BY ITRAN.GT.10 C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION FIELD INDICATED BY INCRAD.NE.0 C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CINPUT' C DIMENSION SP1(MDEP),SP2(MDEP),SP3(MDEP) C T=TAUQ(1) SP1(1)=0. SP2(1)=1.+2.*A1(1) SP3(1)=-2.*A1(1)*A1(1) IF (T.LT.0.01) THEN EX1=T*(1.-T*(0.5-T*(0.1666667-T*0.041666667))) ELSE IF(T.LT.20.) THEN EX1=1.-EXP(-T) ELSE EX1=1. END IF END IF EX=1.-EX1 FACT=1.+2.*A1(1)*EX1 SP2(1)=SP2(1)/FACT SP3(1)=SP3(1)/FACT IF(INCRAD.EQ.0) THEN IMINUS(0)=0.0 P(1)=S(1) ELSE P(1)=S(1)+2.*A1(1)/FACT*EXP(-T)*IMINUS(0) ENDIF C C CALCULATE TRIDIAGONAL COEFFICIENTS C IF(ITRAN.EQ.0) THEN DO 200 K=2,NDEP-1 SP1(K)=-A1(K) SP2(K)=1. SP3(K)=-C1(K) P(K)=S(K) 200 CONTINUE ELSE IF(ITRAN.EQ.1) THEN DO 201 K=2,NDEP-1 AD=.166666666*DTAUQ(K)*2./(DTAUQ(K)+DTAUQ(K+1)) CD=.166666666*DTAUQ(K+1)*2./(DTAUQ(K)+DTAUQ(K+1)) SP1(K)=-A1(K)+AD SP2(K)=1. SP3(K)=-C1(K)+CD P(K)=S(K)+AD*(S(K-1)-S(K))+CD*(S(K+1)-S(K)) 201 CONTINUE ELSE IF(ITRAN.EQ.2) THEN DO 202 K=2,NDEP-1 AD=.166666666-0.083333333*DTAUQ(K+1)**2/DTAUQ(K)* * 2./(DTAUQ(K)+DTAUQ(K+1)) CD=.166666666-0.083333333*DTAUQ(K)**2/DTAUQ(K+1)* * 2./(DTAUQ(K)+DTAUQ(K+1)) SP1(K)=-A1(K)+AD SP2(K)=1. SP3(K)=-C1(K)+CD P(K)=S(K)+AD*(S(K-1)-S(K))+CD*(S(K+1)-S(K)) 202 CONTINUE ELSE CALL STOP(' TRANF: ITRAN OUTSIDE RANGE') ENDIF C C K=NDEP: LOWER BOUNDARY C SP1(NDEP)=-1. SP2(NDEP)=DTAUQ(NDEP)+0.5*DTAUQ(NDEP)**2 SP3(NDEP)=0. P(NDEP)=S(NDEP)*(DTAUQ(NDEP)+0.5*DTAUQ(NDEP)**2)+ * (S(NDEP)-S(NDEP-1)) C C ELIMINATE SUBDIAGONAL C DO 300 K=1,NDEP-1 F=-SP1(K+1)/(SP2(K)-SP3(K)) P(K+1)=P(K+1)+F*P(K) SP2(K+1)=SP2(K+1)+F*SP2(K) SP2(K)=SP2(K)-SP3(K) 300 CONTINUE SP2(NDEP)=SP2(NDEP)-SP3(NDEP) C C BACKSUBSTITUTE C P(NDEP)=P(NDEP)/SP2(NDEP) DO 320 K=NDEP-1,1,-1 P(K)=(P(K)-SP3(K)*P(K+1))/SP2(K) 320 CONTINUE C C SURFACE INTENSITY C IPLUS(0)=2.*(EX*P(1)+S(1)*0.5*EX1**2) - IMINUS(0)*EX*EX C C INTENSITIES, CALCULATED FROM A WEIGHTED AVERAGE OF C DP/DTAUNY(K+0.5) AND DP/DTAUNY(K-0.5) C IMINUS(1)=EX1*S(1) + EXP(-T)*IMINUS(0) IPLUS(1)=2.*P(1)-IMINUS(1) DPDT2=(P(2)-P(1))/DTAUQ(2) DO 340 K=2,NDEP-1 DPDT1=DPDT2 DPDT2=(P(K+1)-P(K))/DTAUQ(K+1) PPRIMK=(DTAUQ(K)*DPDT2+DTAUQ(K+1)*DPDT1)/ * (DTAUQ(K)+DTAUQ(K+1)) IPLUS(K)=P(K)+PPRIMK IMINUS(K)=P(K)-PPRIMK 340 CONTINUE IPLUS(NDEP)=S(NDEP)+(S(NDEP)-S(NDEP-1))/DTAUQ(NDEP) IMINUS(NDEP)=2.0*P(NDEP)-IPLUS(NDEP) C RETURN END C C********************************************************************** C SUBROUTINE TRANSC C C FORMAL SOLVER FOR THE TRANSFER EQUATION USING FEAUTRIER TECHNIQUE C VERSION SOLVING FOR ALL ANGLES AT ONCE TO BE USED IF C SCATTERING DOMINATES C FOR THIS ROUTINE TO WORK IT MUST BE POSSIBLE TO SEPARATE THE C ABSORPTION AND SCATTERING PARTS OF THE SOURCE FUNCTION C THIS IS DONE BY WRITING THE SOURCE FUNCTION AS C S = SABS + SSCAT*JNY C SABS AND SSCAT ARE EXTRACTED FROM THE VARIABLES SBCK, SC, SCAT C AND RNY FROM: C SSCAT(K)=RNY(K)*SCAT(K) C SABS(K)=S(K)-RNY(K)*SBCK(K)+RNY(K)*SC(K) C THESE VARIABLES THUS HAVE TO BE LOADED WITH THE CORRECT VALUES C C ROUTINE HAS NOT BEEN OPTIMIZED, MAJOR SAVINGS PROBABLY POSSIBLE C ITRAN=0 IS THE ONLY FEAUTRIER OPTION IMPLEMENTED C C: TRANSC 95-02-22 NEW ROUTINE: (MATS CARLSSON) C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CINPUT' INCLUDE 'CGAUSI' INCLUDE 'CLU' C DIMENSION SP1(MMU,MMU,MDEP),SP2(MMU,MMU,MDEP),SP3(MMU,MMU,MDEP) DIMENSION PSC(MMU,MDEP),SCJ(MDEP),SABS(MDEP),SSCAT(MDEP) DIMENSION WORK(MMU,MMU),FWORK(MMU,MMU),PWORK(MMU),PWORK2(MMU) LOGICAL LWARN DATA LWARN/.FALSE./ SAVE PSC,SCJ,SABS,SSCAT,LWARN C C SOLVE FOR ALL ANGLES WHEN MU=1 AND STORE RESULTS IN LOCAL VARIABLE C PSC C FOR ALL MU-VALUES, COPY FROM LOCAL VARIABLE INTO MU-DEPENDENT C P C QUANTITIES THAT INCLUDE XMY HAS TO BE ADJUSTED C C IF ITRAN HAS A NON-IMPLEMENTED VALUE, ISSUE A WARNING BUT ONLY ONCE C IF(ITRAN.NE.0 .AND. .NOT.LWARN) THEN WRITE(LJOBLO,*) 'TRANSC: WARNING - ONLY ITRAN=0 IMPLEMENTED' LWARN=.TRUE. ENDIF DO 100 MU=1,NMU IF(XMY.EQ.XMU(MU)) GOTO 110 100 CONTINUE CALL STOP('TRANSC: MU NOT FOUND') 110 CONTINUE MU0=MU IF(MU0.EQ.1) THEN DO 120 K=1,NDEP SSCAT(K)=RNY(K)*SCAT(K) SABS(K)=S(K)-RNY(K)*SBCK(K)+RNY(K)*SC(K) 120 CONTINUE DO 280 MU=1,NMU DO 150 MU2=1,NMU DO 140 K=1,NDEP SP1(MU,MU2,K)=0.0 SP2(MU,MU2,K)=0.0 SP3(MU,MU2,K)=0.0 140 CONTINUE 150 CONTINUE XMU1=XMU(MU)/XMU(1) XMU2=XMU1*XMU1 T=TAUQ(1)/XMU1 SP1(MU,MU,1)=0. SP2(MU,MU,1)=1.+2.*XMU1*A1(1) SP3(MU,MU,1)=-2.*XMU2*A1(1)*A1(1) IF (T.LT.0.01) THEN EX1=T*(1.-T*(0.5-T*(0.1666667-T*0.041666667))) ELSE IF(T.LT.20.) THEN EX1=1.-EXP(-T) ELSE EX1=1. END IF END IF EX=1.-EX1 FACT=1.+2.*XMU1*A1(1)*EX1 SP2(MU,MU,1)=SP2(MU,MU,1)/FACT SP3(MU,MU,1)=SP3(MU,MU,1)/FACT IF(INCRAD.EQ.0) THEN IMINUS(0)=0.0 PSC(MU,1)=SABS(1) ELSE PSC(MU,1)=SABS(1)+2.*XMU1*A1(1)/FACT*EXP(-T)*IMINUS(0) ENDIF C C CALCULATE TRIDIAGONAL COEFFICIENTS C DO 200 K=2,NDEP-1 SP1(MU,MU,K)=-XMU2*A1(K) SP2(MU,MU,K)=1. SP3(MU,MU,K)=-XMU2*C1(K) PSC(MU,K)=SABS(K) 200 CONTINUE C C K=NDEP: LOWER BOUNDARY C SP1(MU,MU,NDEP)=-1. SP2(MU,MU,NDEP)=DTAUQ(NDEP)/XMU1+0.5*DTAUQ(NDEP)**2/XMU2 SP3(MU,MU,NDEP)=0. SK=(DTAUQ(NDEP)/XMU1+0.5*DTAUQ(NDEP)**2/XMU2)+1.0 PSC(MU,NDEP)=SABS(NDEP)*SK - SABS(NDEP-1) C C NON-DIAGONAL SCATTERING ELEMENTS C DO 250 MU2=1,NMU SP2(MU,MU2,1)=SP2(MU,MU2,1)-SSCAT(1)*WMU(MU2) DO 220 K=2,NDEP-1 SP2(MU,MU2,K)=SP2(MU,MU2,K)-SSCAT(K)*WMU(MU2) 220 CONTINUE SP2(MU,MU2,NDEP)=SP2(MU,MU2,NDEP)-SSCAT(NDEP)*WMU(MU2)*SK+ * SSCAT(NDEP-1)*WMU(MU2) SP1(MU,MU2,NDEP)=SP1(MU,MU2,NDEP)+SSCAT(NDEP-2)*WMU(MU2) 250 CONTINUE 280 CONTINUE C C ELIMINATE SUBDIAGONAL C DO 390 K=1,NDEP-1 CTM CALL MATSUB(SP2(1,1,K),SP3(1,1,K),WORK,NMU,NMU,MMU,MMU) CALL MATSUB(SP2(:,:,K),SP3(:,:,K),WORK,NMU,NMU,MMU,MMU) CALL MATINV(WORK,NMU,MMU) CALL MATMUL(SP1(1,1,K+1),WORK,FWORK,NMU,NMU,NMU,MMU,MMU,MMU) CALL MATMUL(FWORK,PSC(1,K),PWORK,NMU,NMU,1,MMU,MMU,1) CALL MATMUL(FWORK,SP2(1,1,K),WORK,NMU,NMU,NMU,MMU,MMU,MMU) DO 380 MU=1,NMU PSC(MU,K+1)=PSC(MU,K+1)-PWORK(MU) DO 370 MU2=1,NMU SP2(MU,MU2,K+1)=SP2(MU,MU2,K+1)-WORK(MU,MU2) SP2(MU,MU2,K)=SP2(MU,MU2,K)-SP3(MU,MU2,K) 370 CONTINUE 380 CONTINUE 390 CONTINUE DO 420 MU=1,NMU DO 410 MU2=1,NMU SP2(MU,MU2,NDEP)=SP2(MU,MU2,NDEP)-SP3(MU,MU2,NDEP) 410 CONTINUE 420 CONTINUE C C BACKSUBSTITUTE C CALL MATINV(SP2(1,1,NDEP),NMU,MMU) CALL MATMUL(SP2(1,1,NDEP),PSC(1,NDEP),PWORK,NMU,NMU,1,MMU,MMU,1) SCJ(NDEP)=0.0 DO 500 MU=1,NMU PSC(MU,NDEP)=PWORK(MU) SCJ(NDEP)=SCJ(NDEP)+WMU(MU)*PSC(MU,NDEP) 500 CONTINUE DO 690 K=NDEP-1,1,-1 CALL MATMUL(SP3(1,1,K),PSC(1,K+1),PWORK,NMU,NMU,1,MMU,MMU,1) CALL MATSUB(PSC(1,K),PWORK,PWORK2,NMU,1,MMU,1) CALL MATINV(SP2(1,1,K),NMU,MMU) CALL MATMUL(SP2(1,1,K),PWORK2,PWORK,NMU,NMU,1,MMU,MMU,1) SCJ(K)=0.0 DO 650 MU=1,NMU PSC(MU,K)=PWORK(MU) SCJ(K)=SCJ(K)+WMU(MU)*PSC(MU,K) 650 CONTINUE 690 CONTINUE ENDIF C C MOVE LOCAL VARIABLES INTO MULTI VARIABLES C DO 700 K=1,NDEP P(K)=PSC(MU0,K) S(K)=SABS(K)+SSCAT(K)*SCJ(K) 700 CONTINUE C T=TAUQ(1) IF (T.LT.0.01) THEN EX1=T*(1.-T*(0.5-T*(0.1666667-T*0.041666667))) ELSE IF(T.LT.20.) THEN EX1=1.-EXP(-T) ELSE EX1=1. END IF END IF EX=1.-EX1 C C SURFACE INTENSITY C IPLUS(0)=2.*(EX*P(1)+S(1)*0.5*EX1**2) - IMINUS(0)*EX*EX C C INTENSITIES, CALCULATED FROM A WEIGHTED AVERAGE OF C DP/DTAUNY(K+0.5) AND DP/DTAUNY(K-0.5) C IMINUS(1)=EX1*S(1) + EXP(-T)*IMINUS(0) IPLUS(1)=2.*P(1)-IMINUS(1) DPDT2=(P(2)-P(1))/DTAUQ(2) DO 800 K=2,NDEP-1 DPDT1=DPDT2 DPDT2=(P(K+1)-P(K))/DTAUQ(K+1) PPRIMK=(DTAUQ(K)*DPDT2+DTAUQ(K+1)*DPDT1)/ * (DTAUQ(K)+DTAUQ(K+1)) IPLUS(K)=P(K)+PPRIMK IMINUS(K)=P(K)-PPRIMK 800 CONTINUE IPLUS(NDEP)=S(NDEP)+(S(NDEP)-S(NDEP-1))/DTAUQ(NDEP) IMINUS(NDEP)=2.0*P(NDEP)-IPLUS(NDEP) C RETURN END C C********************************************************************** C SUBROUTINE TRANI C C FORMAL SOLVER FOR THE TRANSFER EQUATION USING DIRECT INTEGRATION C FOR CUBIC SPLINE SOURCES C C ITRAN DETERMINES THE MODE OF FORMAL SOLUTION C = 3 INTEGRAL CUBIC SPLINE METHOD ACCORDING TO SCHARMER C 4 INTEGRAL CUBIC SPLINE METHOD WITH LOCALLY DETERMINED C FIRST DERIVATIVE C C REF G.B.SCHARMER C C: C: TRANI 86-09-04 NEW ROUTINE: (GORAN SCHARMER) C: INTEGRAL CUBIC SPLINE FORMAL SOLVER C: C: 92-10-08 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION FIELD IS SPECIFIED IN IMINUS(0) C: THIS IS INDICATED BY ITRAN.GT.10 C: C: 95-08-17 MODIFICATIONS: (MATS CARLSSON) C: INCOMING RADIATION FIELD INDICATED BY INCRAD.NE.0 C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CINPUT' C DIMENSION SPRIM(MDEP),SBISS(MDEP),STRISS(MDEP), * EXPD(MDEP),WW(MDEP) INTEGER KP(MDEP) DIMENSION AIP(0:MDEP),AIM(0:MDEP) DATA DT1/1.E-2/,DT2/1.E2/ C C THE PARAMETER DT1 DETERMINES WHEN I HAS TO BE CALCULATED IN A C THIN WAY DUE TO NUMERICAL REASONS. THE PARAMETER SHOULD BE SET C TO ABOUT THE THIRD ROOT OF THE MACHINE ACCURACY. C IF(INCRAD.EQ.0) IMINUS(0)=0.0 T=TAUQ(1) C IF(ITRAN.EQ.3) THEN CALL SPLIN0(NDEP,DTAUQ,S,SPRIM,SBISS,STRISS,WW) ELSE CALL SPLIN1(NDEP,DTAUQ,S,SPRIM,SBISS,STRISS) ENDIF C C CALCULATE ONLY NECESSARY EXPONENTIALS, AND DO IT IN A SEPARATE C LOOP TO ALLOW THE CRAY TO VECTORIZE THE CALLS C NM=0 DO 400 K=1,NDEP-1 IF(DTAUQ(K+1).GT.DT1 .AND. DTAUQ(K+1).LE.DT2) THEN NM=NM+1 WW(NM)=DTAUQ(K+1) KP(NM)=K ENDIF 400 CONTINUE DO 410 M=1,NM WW(M)=EXP(-WW(M)) 410 CONTINUE DO 420 M=1,NM EXPD(KP(M))=WW(M) 420 CONTINUE C C CALCULATE OUTGOING INTENSITY IPLUS(K) C AIP(NDEP)=S(NDEP)+SPRIM(NDEP)+SBISS(NDEP)+STRISS(NDEP) IPLUS(NDEP)=AIP(NDEP) DO 500 K=NDEP-1,1,-1 IF(DTAUQ(K+1).GT.DT2) THEN AIP(K)=S(K)+SPRIM(K)+SBISS(K)+STRISS(K) ELSE IF (DTAUQ(K+1).GT.DT1) THEN SBISS1=SBISS(K)+DTAUQ(K+1)*STRISS(K) AIP(K)=S(K)+SPRIM(K)+SBISS(K)+STRISS(K)+EXPD(K)* * (AIP(K+1)-S(K+1)-SPRIM(K+1)-SBISS1-STRISS(K)) ELSE SBISS1=SBISS(K)+DTAUQ(K+1)*STRISS(K) AIP(K)=AIP(K+1)-DTAUQ(K+1)*(AIP(K+1)-S(K+1)- * 0.5*DTAUQ(K+1)*(AIP(K+1)-S(K+1)-SPRIM(K+1)- * 0.333333333*DTAUQ(K+1)*(AIP(K+1)-S(K+1)-SPRIM(K+1)- * SBISS1-0.25*DTAUQ(K+1)*(IPLUS(K+1)-SPRIM(K+1)-SBISS1- * STRISS(K))))) ENDIF IPLUS(K)=AIP(K) 500 CONTINUE C C SURFACE INTENSITY C IPLUS(0)=IPLUS(1)+T*(1.-T*(.5-T*(1./6.-T/24.)))*(S(1)-IPLUS(1)) C C CALCULATE INCOMING INTENSITY C IF(T.LT.0.01) THEN AIM(1)=IMINUS(0)+T*(1.-T*(.5-T*(1./6.-T/24.)))* * (S(1)-IMINUS(0)) ELSE IF(T.LT.20.) THEN AIM(1)=IMINUS(0)+(1.-EXP(-T))*(S(1)-IMINUS(0)) ELSE AIM(1)=S(1) ENDIF IMINUS(1)=AIM(1) DO 600 K=1,NDEP-1 IF(DTAUQ(K+1).LE.DT1) THEN AIM(K+1)=AIM(K)-DTAUQ(K+1)*(AIM(K)-S(K)- * .5*DTAUQ(K+1)*(AIM(K)-S(K)+SPRIM(K)-.333333333* * DTAUQ(K+1)*(AIM(K)-S(K)+SPRIM(K)-SBISS(K)-.25* * DTAUQ(K+1)*(AIM(K)+SPRIM(K)-SBISS(K)+STRISS(K))))) ELSE IF(DTAUQ(K+1).LE.DT2) THEN SBISS1=SBISS(K)+DTAUQ(K+1)*STRISS(K) AIM(K+1)=S(K+1)-SPRIM(K+1)+SBISS1-STRISS(K)+EXPD(K)* * (AIM(K)-S(K)+SPRIM(K)-SBISS(K)+STRISS(K)) ELSE SBISS1=SBISS(K)+DTAUQ(K+1)*STRISS(K) AIM(K+1)=S(K+1)-SPRIM(K+1)+SBISS1-STRISS(K) ENDIF IMINUS(K+1)=AIM(K+1) 600 CONTINUE C C FEAUTRIERS P C DO 700 K=1,NDEP P(K)=0.5*(IPLUS(K)+IMINUS(K)) 700 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE SPLIN0(N,DX,F,D,D2,D3,WW) C C CALCULATES STANDARD CUBIC SPLINE WITH CONTINOUS SECOND DERIVATIVE C C CODED BY G.SCHARMER, 1982 C C DX(K)=X(K)-X(K-1) C WW IS A WORKING ARRAY C INCLUDE 'PREC' DIMENSION DX(N),F(N) DIMENSION D(N),D2(N),D3(N),WW(N) C FAC=-DX(2)/DX(3) D(2)=(2.-FAC)*(DX(2)+DX(3)) C2=DX(3)+FAC*DX(2) WW(2)=(F(2)-F(1))/DX(2) WW(3)=(F(3)-F(2))/DX(3) D3(2)=6.*(WW(3)-WW(2)) DO 100 K=3,N-1 D(K)=2.*(DX(K)+DX(K+1)) WW(K+1)=(F(K+1)-F(K))/DX(K+1) FAC=-DX(K)/D(K-1) D(K)=D(K)+FAC*DX(K) IF(K.EQ.3) D(K)=D(K)+FAC*(C2-DX(K)) D3(K)=6.*(WW(K+1)-WW(K))+FAC*D3(K-1) 100 CONTINUE FAC=-DX(N)/D(N-2) AN=-DX(N-1)-DX(N)+FAC*DX(N-1) D3(N)=FAC*D3(N-2) FAC=-AN/D(N-1) D2(N)=(D3(N)+FAC*D3(N-1))/(DX(N-1)+FAC*DX(N)) DO 150 K=N-1,3,-1 D2(K)=(D3(K)-DX(K+1)*D2(K+1))/D(K) 150 CONTINUE D2(2)=(D3(2)-C2*D2(3))/D(2) D2(1)=((DX(2)+DX(3))*D2(2)-DX(2)*D2(3))/DX(3) DO 180 K=1,N-1 D(K)=WW(K+1)-(D2(K+1)+D2(K)+D2(K))*DX(K+1)/6. D3(K)=(D2(K+1)-D2(K))/DX(K+1) 180 CONTINUE D3(N)=D3(N-1) D(N)=D(N-1)+DX(N)*(D2(N-1)+0.5*DX(N)*D3(N-1)) C RETURN END C C********************************************************************** C SUBROUTINE SPLIN1(N,DX,F,D,D2,D3) C C CALCULATES CUBIC SPLINES WITH LOCALLY DETERMINED FIRST DERIVATIVE C STABILITY BETTER THAN FOR STANDARD CUBIC SPLINE, AT COST OF C DISCONTINOUS SECOND DERIVATIVE. C C CODED BY AA.NORDLUND, MAR-83 C C DX(K)=X(K)-X(K-1) C INCLUDE 'PREC' DIMENSION DX(N),F(N) DIMENSION D(N),D2(N),D3(N) C C FIRST DERIVATIVE BY CENTERED DIFFERENCE C DO 100 K=2,N-1 D(K)=(F(K+1)-F(K-1))/(DX(K+1)+DX(K)) 100 CONTINUE D(1)=(F(2)-F(1))/DX(2) D(N)=(F(N)-F(N-1))/DX(N) C C SECOND AND THIRD DERIVATIVE FROM SPLINE CONDITIONS C DO 110 K=1,N-1 CX=1.0/DX(K+1) DFDX=(F(K+1)-F(K))*CX D2(K)=(6.*DFDX-4.*D(K)-2.*D(K+1))*CX D3(K)=6.*(D(K)+D(K+1)-2.*DFDX)*CX*CX 110 CONTINUE CXN=1.0/DX(N) DFDXN=(F(N)-F(N-1))*CXN D2(N)=(4.*D(N)+2.*D(N-1)-6.*DFDXN)*CXN D3(N)=D3(N-1) C RETURN END C C*********************************************************************** C SUBROUTINE MATMUL(A,B,C,N1,N2,N3,M1,M2,M3) C C GIVES A#B=C WHERE DIMENSIONS ARE C A: N1 ROWS BY N2 COLUMNS C B: N2 ROWS BY N3 COLUMNS C C: N1 ROWS BY N3 COLUMNS C C DIMENSIONS OF ARRAYS ARE M1,M2,M3, ACUTAL USED DIMENSIONS N1,N2,N3 C INCLUDE 'PREC' C DIMENSION A(M1,M2),B(M2,M3),C(M1,M3) C DO 300 I3=1,N3 DO 200 I1=1,N1 C(I1,I3)=0.0 DO 100 I2=1,N2 C(I1,I3)=C(I1,I3)+A(I1,I2)*B(I2,I3) 100 CONTINUE 200 CONTINUE 300 CONTINUE C END C C**************************************************************** C SUBROUTINE MATADD(A,B,C,N1,N2,M1,M2) C C GIVES A+B=C WHERE DIMENSIONS ARE C A: N1 ROWS BY N2 COLUMNS C B: N1 ROWS BY N2 COLUMNS C C: N1 ROWS BY N2 COLUMNS C C DIMENSIONS OF ARRAYS ARE M1,M2 ACUTAL USED DIMENSIONS N1,N2 C INCLUDE 'PREC' C DIMENSION A(M1,M2),B(M1,M2),C(M1,M2) C DO 200 I2=1,N2 DO 100 I1=1,N1 C(I1,I2)=A(I1,I2)+B(I1,I2) 100 CONTINUE 200 CONTINUE C END C C**************************************************************** C SUBROUTINE MATSUB(A,B,C,N1,N2,M1,M2) C C GIVES A-B=C WHERE DIMENSIONS ARE C A: N1 ROWS BY N2 COLUMNS C B: N1 ROWS BY N2 COLUMNS C C: N1 ROWS BY N2 COLUMNS C C DIMENSIONS OF ARRAYS ARE M1,M2 ACUTAL USED DIMENSIONS N1,N2 C INCLUDE 'PREC' C DIMENSION A(M1,M2),B(M1,M2),C(M1,M2) C DO 200 I2=1,N2 DO 100 I1=1,N1 C(I1,I2)=A(I1,I2)-B(I1,I2) 100 CONTINUE 200 CONTINUE C END C C**************************************************************** C SUBROUTINE matinv(a,n,np) c c adapted from gaussj of numerical recipes c include 'PREC' c INTEGER n,np,NMAX dimension a(np,np) C ADD TM C PARAMETER (NMAX=600) C PARAMETER (NMAX=50) C END ADD TM C INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX) INTEGER i,icol,irow,j,k,l,ll,indxc(np),indxr(np),ipiv(np) c do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0. do 13 j=1,n if(ipiv(j).ne.1)then do 12 k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big)then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then STOP 'singular matrix in matinv' endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.) STOP 'singular matrix in matinv' pivinv=1./a(icol,icol) a(icol,icol)=1. do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue do 21 ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0. do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue endif 21 continue 22 continue do 24 l=n,1,-1 if(indxr(l).ne.indxc(l))then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return END C (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw#