SUBROUTINE BMAT(W,NK1,NDEP1,J,I) C C FILLS THE DELTA LAMBDA ELEMENTS IN THE PRECONDITIONED MATRIX W C IN THE EQUATION C C W*DN/N = E C C DN IS THE APPROXIMATE CORRECTION, REQUIRED TO MAKE THE POPULATIONS C SATISFY THE EQUATIONS OF STATISTICAL EQUILIBRIUM. C C DELTA LAMBDA IS CALCULATED AS 1-LAMBDA WHERE 1 IS THE IDENTITY C OPERATOR. THIS IS DONE WHEN TAUQ IS GREATER THAN THIN AND C DELTA TAUQ IS LARGER THAN DIFF. IN THE OPTICALLY THIN PART, C LAMBDA IS SET TO ZERO, IN THE THICK PARTS THE DIFFUSION APPROXIMATION C IS USED. C C VARIABLES: C C K1 : FIRST DEPTH POINT FOR WHICH TAUQ .GT. THIN C C TAUQ : MONOCHROMATIC OPTICAL DEPTH C DTAUQ : DTAUQ(K)=TAUQ(K)-TAUQ(K-1) C TP : QUADRATURE POINT SUCH THAT I(TAUQ)=S(TP) (OUTGOING RAYS) C TM : DITTO SUCH THAT I(TAUQ)=WM*S(TM) (INGOING RAYS) C WM : WEIGHT FOR INGOING RAYS C KP : DEPTH INDEX SUCH THAT TAUQ(KP) .LT. TP .LT. TAUQ(KP+1) C KM : DEPTH INDEX SUCH THAT TAUQ(KM) .LT. TM .LT. TAUQ(KM+1) C AP0 : INTERPOLATION COEFFICIENT FOR OUTGOING RAYS C AM0 : INTERPOLATION COEFFICIENT FOR INGOING RAYS C RNY : CONTINUUM OPACITY RELATIVE TO TOTAL OPACITY C HN3C2 : 2*H*NY**3/C **2 C HNY4P : H*NY/4PI NY IN UNITS OF A TYPICAL DOPPLER WIDTH C ALFK : ALFA/(TOTAL OPACITY) C CIP : DSNY/DNI OUTGOING RADIATION C CIM : INGOING RADIATION C CJP : DSNY/DNJ OUTGOING RADIATION C CJM : INGOING RADIATION C A1 : DIFFUSION TERM FROM TRANSP C C1 : DIFFUSION TERM FROM TRANSP C C CODED BY: M.CARLSSON (FEB.1983). C C: C: BMAT 88-02-04 MODIFICATIONS: (MATS CARLSSON) C: BOUNDARY CONDITION FOR DTAUQ.GT.DIFF AND K=1 CORRECTED C: C: 90-07-31 MODIFICATIONS: (MATS CARLSSON) C: CHANGED TO CONTAIN DECLARATION OF W C: W 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 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CINPUT' C DIMENSION W(NK1,NDEP1,NK1,NDEP1) C DIMENSION CIP(0:MDEP),CIM(0:MDEP),CJP(0:MDEP),CJM(0:MDEP) DATA CIP(0),CIM(0),CJP(0),CJM(0)/4*0./ C C CALCULATE SOME CONSTANTS AND DEPTH INDICES C C ON A VECTOR-TYPE MACHINE, THE IF TEST SHOULD BE IN A C SEPARATE LOOP TO ALLOW VECTORIZATION OF LOOP 100 C WC=-0.5*WQMU/HNY4P K1=0 DO 100 K=1,NDEP ALFK=ALFA(K)/(X(K)*XNORM(K)) CIP(K)=-ALFK*IPLUS(K) CIM(K)=-ALFK*IMINUS(K) CJP(K)=(HN3C2+IPLUS(K))*GIJ(K)*ALFK CJM(K)=(HN3C2+IMINUS(K))*GIJ(K)*ALFK IF(TAUQ(K).LT.THIN) K1=K 100 CONTINUE KP=1 KM=1 C C TAUQ(K) .LT. THIN . SET LAMBDA=0 (I.E. LAMBDA ITERATE) C THIS LOOP VECTORIZES C DO 120 K=1,K1 DW=WC *ALFA(K)*Z(K) C C ADD IDENTITY MATRIX ELEMENTS C W(I,K,J,K)=W(I,K,J,K)-DW*(CJP(K)+CJM(K)) W(J,K,I,K)=W(J,K,I,K)+DW*(CIP(K)+CIM(K)) 120 CONTINUE C C TAUQ(K) .GT. THIN C COMMENTS ON POSSIBILITIES OF INCREASING VECTORIZATION: C THIS LOOP DOES NOT VECTORIZE. THERE ARE TWO DIFFICULTIES: C 1) THE DETERMINATION OF QUADRATURE POINTS FOR THE SCHARMER OPERATOR C INVOLVES AN INNER DO-LOOP WITH AN IF-TEST. C 2) THERE ARE TWO REGIMES: THE SCHARMER OPERATOR REGIME C AND THE DIFFUSION REGIME. IN GENERAL THEY HAVE SOME OVERLAP C IT IS POSSIBLE TO INCREASE VECTORIZATION BY: C A) SETTING UP THREE INTERVALS: C 1) SCHARMER OPERATOR INTERVAL C 2) DIFFUSION APPROXIMATION INTERVAL C 3) INTERVAL WHERE THEY OVERLAP C B) HAVE THE TEST OF THE LIMITS OF THE THREE INTERVALS IN SEPARATE C LOOPS FIRST C C) DETERMINE THE SCHARMER OPERATOR QUADRATURE POINTS KM AND KP IN C SEPARATE LOOPS AND STORE THEM IN ARRAYS C D) TREAT INTERVALS 1, 2 AND 3 IN SEPARATE LOOPS, ONLY INTERVAL 3 C DO NOT VECTORIZE. WITH INDIRECT ADDRESSING ALSO THIS LOOP CAN C GET VECTORIZED. C THE LOOPS IN B) AND C) WILL NORMALLY NOT VECTORIZE. ON THE AMDAHL C (FUJITSU) MACHINE A LOOP WITH A GOTO CAN VECTORIZE AND THUS ALSO C LOOPS UNDER B) AND C). ALMOST 100 PERCENT OF THE ROUTINE WILL THEN C BE IN A VECTORIZABLE FORM. C DO 150 K=K1+1,NDEP-1 T=TAUQ(K) DW=WC *ALFA(K)*Z(K) IF(DTAUQ(K).LT.DIFF .AND. DTAUQ(K+1).LT.DIFF) THEN C C INCOMING RADIATION C WM=1.-EXP(-T) TM=T/WM-1. DO 130 L=KM,NDEP-1 IF (TAUQ(L+1).GT.TM) GO TO 131 130 CONTINUE 131 KM=L AM0=(TAUQ(KM+1)-TM)/DTAUQ(KM+1) IF(K.EQ.1) THEN KM=1 AM0=1. ENDIF C W(I,K,J,KM)=W(I,K,J,KM)+DW*CJM(KM)*AM0*WM W(I,K,J,KM+1)=W(I,K,J,KM+1)+DW*CJM(KM+1)*(1.-AM0)*WM C W(J,K,I,KM)=W(J,K,I,KM)-DW*CIM(KM)*AM0*WM W(J,K,I,KM+1)=W(J,K,I,KM+1)-DW*CIM(KM+1)*(1.-AM0)*WM C C OUTGOING RADIATION C TP=T+1 DO 140 L=KP,NDEP-1 IF (TAUQ(L+1).GT.TP) GO TO 141 140 CONTINUE L=NDEP-2 141 KP=L AP0=(TAUQ(KP+1)-TP)/DTAUQ(KP+1) C W(I,K,J,KP)=W(I,K,J,KP)+DW*CJP(KP)*AP0 W(I,K,J,KP+1)=W(I,K,J,KP+1)+DW*CJP(KP+1)*(1.-AP0) C W(J,K,I,KP)=W(J,K,I,KP)-DW*CIP(KP)*AP0 W(J,K,I,KP+1)=W(J,K,I,KP+1)-DW*CIP(KP+1)*(1.-AP0) C C ADD IDENTITY MATRIX ELEMENTS C W(I,K,J,K)=W(I,K,J,K)-DW*(CJP(K)+CJM(K)) W(J,K,I,K)=W(J,K,I,K)+DW*(CIP(K)+CIM(K)) C ELSE IF(K.EQ.1) THEN C C BOUNDARY - DIFFUSION APPROXIMATION CAN NOT BE USED C FOR TAUQ.GE.20 SET I=S (DELTA LAMBDA=0, NO ELEMENTS FILLED IN) C FOR TAUQ.LT.20 SET I- = S*(1-EXP(-TAUQ)) C I+ = S(1) + (S(2)-S(1))/(TAUQ(2)-TAUQ(1)) IF(T.LT.20.) THEN W(I,K,J,K)=W(I,K,J,K)-DW*EXP(-T)*CJM(K) W(J,K,I,K)=W(J,K,I,K)+DW*EXP(-T)*CIM(K) W(I,K,J,K)=W(I,K,J,K)-DW*CJP(K)/DTAUQ(K+1) W(J,K,I,K)=W(J,K,I,K)+DW*CIP(K)/DTAUQ(K+1) W(I,K,J,K+1)=W(I,K,J,K+1)+DW*CJP(K+1)/DTAUQ(K+1) W(J,K,I,K+1)=W(J,K,I,K+1)-DW*CIP(K+1)/DTAUQ(K+1) ENDIF ELSE C C MAKE DIFFUSION APPROXIMATION WHEN DTAU.GT.DIFF C AM0=A1(K) AP0=C1(K) C W(I,K,J,K-1)=W(I,K,J,K-1)+DW*(CJP(K-1)+CJM(K-1))*AM0 W(I,K,J,K)=W(I,K,J,K)-DW*(CJP(K)+CJM(K))*(AM0+AP0) W(I,K,J,K+1)=W(I,K,J,K+1)+DW*(CJP(K+1)+CJM(K+1))*AP0 C W(J,K,I,K-1)=W(J,K,I,K-1)-DW*(CIP(K-1)+CIM(K-1))*AM0 W(J,K,I,K)=W(J,K,I,K)+DW*(CIP(K)+CIM(K))*(AM0+AP0) W(J,K,I,K+1)=W(J,K,I,K+1)-DW*(CIP(K+1)+CIM(K+1))*AP0 ENDIF 150 CONTINUE RETURN END C C*********************************************************************** C SUBROUTINE ITER(W,E,NK1,NDEP1) C C ADMINISTERS ITERATION. C FILLS THE MATRICES AND SOLVES THE SYSTEM OF EQUATIONS. C: C: ITER 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 STATEQ C: C: 88-04-21 MODIFICATIONS: (MATS CARLSSON) C: VERSION THAT IMPLEMENTS COLLISIONAL-RADIATIVE SWITCHING, C: HUMMER,D.G., VOELS,S.A.: 1988, ASTRON. ASTROPHYS. 192,279 C: SWITCHING IS REGULATED BY ICONV. C: ICONV=2 GIVES INTERACTIVE MODE C: 3 GIVES AUTOMATIC MODE C: C: 89-03-18 MODIFICATIONS: (MATS CARLSSON) C: EMAX PRINTOUT CHANGED TO INCLUDE SIGN AND C: LEVEL/TRANSITION AND DEPTH C: C: 89-06-06 MODIFICATIONS: (MATS CARLSSON) C: HSEINT ONLY CALLED IF SWITCH=1.0 C: C: 89-06-07 MODIFICATIONS: (MATS CARLSSON) C: WHEN SWITCHING WAS USED COLLISIONAL RATES WERE TAKEN C: FROM FIRST CALCULATED VALUES THROUGH VARIABLE COL READ C: FROM DUMC FILE. THIS GAVE ERRORS FOR HYDROGEN WHEN HSE C: INTEGRATIONS WERE PERFORMED SINCE THE UPDATED COLLISIONAL C: RATES WERE OVERWRITTEN WITH THE OLD VALUES. FIXED SO C: THAT COL IS ONLY USED WHEN SWITCH.GT.1.0 AND THE FIRST C: TIME SWITCH=1.0. SECOND CASE CONTROLLED BY VARIABLE OLDSW C: C: 90-01-11 MODIFICATIONS: (MATS CARLSSON) C: X(K) MAY BE 0 IF STIMULATED EMISSION EXCEEDS ABSORPTION C: X(K) IS SET TO SIGN(1,X(K))*MAX(ABS(X(K)),1.E-4*XCONT(K)) C: C: 90-07-31 MODIFICATIONS: (MATS CARLSSON) C: CHANGED TO CONTAIN DECLARATION OF W AND E C: THESE PASSED IN ARGUMENT LIST INSTEAD OF IN COMMON TO C: AVOID NK=MK1 AND NDEP=MDEP1 REQUIREMENT C: C: 90-10-16 MODIFICATIONS: (MATS CARLSSON) C: IMPLEMENTED NG ACCELERATION C: C: 90-10-18 MODIFICATIONS: (MATS CARLSSON) C: NG ACCELERATION ALSO INCLUDED IN HSE INTEGRATION. C: NG ACCELERATION IS SWITCHED ON WHEN HSE STEPS ARE C: PERFORMED FOR EACH ITERATION (EACH HSE STEP GIVES RISE C: TO EMAX.LT.ELIM1). THIS IS ACHIEVED THROUGH ITERATION C: COUNT VARIABLE ITHSE C: C: 92-06-05 MODIFICATIONS: (MATS CARLSSON) C: TEST OF KREJ AND KEJ INCLUDED TO AVOID ACCESSING C: EJ(0,0) WHEN NRAD=0 C: C: 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: C: 95-08-22 MODIFICATIONS: (MATS CARLSSON) C: COLLISIONAL RADIATIVE SWITCHING NOW SWITCHED ON WITH ICRSW C: ICRSW=-2 GIVES AUTOMATIC MODE C: ICRSW=-1 GIVES INTERACTIVE MODE C: ICRSW.GT.0 GIVES DECREASE OF C: THE SWITCHING PARAMETER WITH ICRSW STEPS PER DECADE C: IN THE INTERACTIVE MODE, A REDO OF THE PREVIOUS ITERATION IS ASKED C: FOR WITH A NEGATIVE VALUE OF SWITCH. ZERO CAUSES SWITCH TO C: AUTOMATIC MODE C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CTRAN' INCLUDE 'CSLINE' INCLUDE 'CGAUSI' INCLUDE 'CCONST' INCLUDE 'CINPUT' INCLUDE 'CLU' INCLUDE 'CIMIN' COMMON /CNGACC/ NGACC LOGICAL NGACC C DIMENSION W(NK1,NDEP1,NK1,NDEP1),E(NK1,NDEP1) CHARACTER*14 ETEXT DIMENSION COL(MK,MK,MDEP) DIMENSION FIX(MK,MK,MDEP) DIMENSION DEIK(MDEP) DIMENSION AN(MK,MDEP) LOGICAL NEWMAT,CONT,HSE LOGICAL LISUM,NGON,NGHSE,LWARN INTEGER ISUMA(MDEP) C CALL CPUTIME('SUM IT',0,0,3) ONE=1.0 HSE=ATOMID.EQ.'H ' .AND. IHSE.NE.0 EMAX=10. SWITCH=1. OLDSW=SWITCH EREDO=0.0 FS0=9. DFSMAX=1.3 NGON=.FALSE. NGHSE=.FALSE. ITHSE=0 C C SWITCHING INITIALIZATION C IF(ICRSW.NE.0) THEN C C CALL TRPT TO GET RIJ AND RJI C CALL TRPT C C FIND INITIAL SWITCHING PARAMETER VALUE C CALL CLOSE(LDUMC) CALL OPEN(LDUMC,'DUMC',NK*NK,'OLD') DO 90 K=1,NDEP READ(LDUMC,REC=K) ((COL(I,J,K),I=1,NK),J=1,NK) DO 80 J=1,NK DO 70 I=1,NK FIX(I,J,K)=C(I,J,K)-COL(I,J,K) 70 CONTINUE 80 CONTINUE 90 CONTINUE DO 150 KR=1,NLINE I=IRAD(KR) J=JRAD(KR) DO 100 K=1,NDEP SWITCH=MAX(SWITCH,RIJ(K,KR)/COL(I,J,K)) SWITCH=MAX(SWITCH,RJI(K,KR)/COL(J,I,K)) 100 CONTINUE 150 CONTINUE IF(ICRSW.EQ.-1) THEN WRITE(*,160) ' MAX(Rij/Cij) IS:',SWITCH 160 FORMAT(A,1P,E10.2) WRITE(*,*) ' GIVE INITIAL SWITCH VALUE TO BE USED' READ(*,*) SWITCH ELSE SWFAC=10.**(-1./ICRSW) SWITCH=SWITCH*5. ENDIF ENDIF C DO 9000 IT=1,ITMAX ITHSE=ITHSE+1 IF(SWITCH.GT.1.0 .OR. OLDSW.NE.1.0) THEN C C CALCULATE SWITCH SCALED C C DO 200 K=1,NDEP DO 190 J=1,NK DO 180 I=1,NK C(I,J,K)=COL(I,J,K)*SWITCH + FIX(I,J,K) 180 CONTINUE 190 CONTINUE 200 CONTINUE ENDIF C CALL CPUTIME('ITER ',0,0,2) EMAXJ=0.0 NEWMAT=EMAX.GT.ELIM1 .OR. SWITCH.NE.1. CALL CPUTIME('ZERO W',0,0,1) C C ZEROING OF MATRICES E AND W IN ROUTINE ZMAT C IF(NEWMAT) THEN CALL ZMAT(W,NK*NDEP*NK*NDEP) ENDIF CALL ZMAT(E,NK*NDEP) CALL CPUTIME('ZERO W',0,1,1) C C **** START OF DETAILED TRANSITIONS **** C CALCULATE THE ERROR TERM E BY SOLVING RADIATIVE TRANSFER IN DETAIL C CALCULATE MATRIX W C CALL CPUTIME('BMATTR',0,0,1) CALL REWIND(LOPC) CALL REWIND(LPHI) IREC=0 DO 400 KR=1,NRAD LWARN=.FALSE. I=IRAD(KR) J=JRAD(KR) C C CALCULATE SOURCE FUNCTIONS, LINE OPACITIES AND CROSS-SECTIONS C DO 270 K=1,NDEP DEIK(K)=0. EJ(KR,K)=0.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).LT.0.0) THEN IF(.NOT.LWARN) THEN LWARN=.TRUE. IF(IWARN.GE.2) WRITE(LJOBLO,278) KR,IT 278 FORMAT(' NEGATIVE OPACITIES, KR=',I3,' IT=',I3) ENDIF ENDIF SL(K,KR)=HN3C2*N(J,K)*GIJK/Z(K) GIJ(K)=GIJK 280 CONTINUE ENDIF IF(IWIDE(KR)) KT=KTRANS(KR) 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) JNYOLD(K)=JNY(K) JNY(K)=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 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) X(K)=SIGN(ONE,X(K))*MAX(ABS(X(K)),1.E-4*XCONT(K)) S(K)=(1.-RNY(K))*SL(K,KR)+RNY(K)*SBCK(K) 310 CONTINUE ENDIF C C CALCULATE INTENSITIES AND FILL NON-LOCAL ELEMENTS C IF(INCRAD.NE.0) IMINUS(0)=XIMIN(NY,MU,KR) CALL TRANSP IF(NEWMAT) CALL BMAT(W,NK1,NDEP1,J,I) C C RADIATION PART TO DIAGONAL IN W C IF(CONT) THEN WQMUH=WQMU/HNY4P/FRQ(NY,KT)*FRQ(0,KT) ELSE WQMUH=WQMU/HNY4P ENDIF IF(NEWMAT) THEN DO 340 K=1,NDEP W(I,K,J,K)=W(I,K,J,K)+WQMUH*ALFA(K)*GIJ(K)* * RNY(K)*(P(K)+HN3C2) W(J,K,I,K)=W(J,K,I,K)+WQMUH*ALFA(K)* * RNY(K)*P(K) 340 CONTINUE ENDIF C C RADIATION PART TO ERROR E C DO 360 K=1,NDEP DEIK(K)=DEIK(K)+WQMUH*ALFA(K)*(Z(K)*(PMS(K)+ * RNY(K)*SBCK(K))-N(J,K)*HN3C2*GIJ(K)*RNY(K)) JNY(K)=JNY(K)+WMU(MU)*P(K) 360 CONTINUE 370 CONTINUE IF(IWTEST.LT.0) CALL WTEST(KR,NY) IF(IDLNY.LT.0) CALL WIDLNY(KR,NY) C C CHECK CONVERGENCE OF JNY AND WRITE NEW JNY TO FILE C DO 375 K=1,NDEP CPGJ IF-LOOP ADDED BY PGJ: IF(JNY(K) .GT. 0.)THEN DJ=(JNY(K)-JNYOLD(K))/JNY(K) ELSE DJ=0. ENDIF IF(ABS(DJ).GT.ABS(EJ(KR,K))) EJ(KR,K)=DJ 375 CONTINUE CALL WRITEJ(IREC) C 380 CONTINUE DO 390 K=1,NDEP E(I,K)=E(I,K)+DEIK(K)*WPHI(K,KR) E(J,K)=E(J,K)-DEIK(K)*WPHI(K,KR) 390 CONTINUE IF(IWTEST.LT.0 .AND. IWEVEC.GT.0) CALL WEVEC(E,NK1,NDEP1) 400 CONTINUE CALL CPUTIME('BMATTR',0,1,1) C* C* 89-03-18 START MODIFICATION C STORE KR INDEX FOR MAX CHANGE IN KREJ, K INDEX IN KEJ C KREJ=0 KEJ=0 DO 404 KR=1,NRAD DO 402 K=1,NDEP IF(ABS(EJ(KR,K)).GT.EMAXJ) THEN KREJ=KR KEJ=K EMAXJ=ABS(EJ(KR,K)) ENDIF 402 CONTINUE 404 CONTINUE IF(KREJ.GT.0 .AND. KEJ.GT.0) EMAXJ=EJ(KREJ,KEJ) WRITE(ETEXT,405) KREJ,KEJ 405 FORMAT('EMAXJ(',I3,',',I3,')') C IF(IT.GT.1) THEN CALL WCHANG(E,NK1,NDEP1,2) CALL WEMAX(ETEXT,EMAXJ) ENDIF EMAXJ=ABS(EMAXJ) C* 89-03-18 END MODIFICATION C* C* CALCULATE ISUM C* LISUM=ISUM.EQ.0 IF(LISUM) THEN DO 408 K=1,NDEP POPMAX=N(1,K) ISUMA(K)=1 DO 406 I=2,NK IF(N(I,K).GT.POPMAX) THEN POPMAX=N(I,K) ISUMA(K)=I ENDIF 406 CONTINUE 408 CONTINUE ENDIF C C NORMALISE AND ADD INTO SMALL DIAGONALS C IF(NEWMAT) THEN DO 440 I=1,NK DO 430 J=1,NK IF(KRAD(I,J).EQ.0) GOTO 430 KR=KRAD(I,J) DO 420 L=1,NDEP DO 410 K=1,NDEP W(I,K,J,L)=W(I,K,J,L)*WPHI(K,KR)*N(J,L) W(J,K,J,L)=W(J,K,J,L)-W(I,K,J,L) 410 CONTINUE 420 CONTINUE 430 CONTINUE 440 CONTINUE CALL CPUTIME('NORM ',0,1,1) C C **** END OF DETAILED TRANSITIONS **** C C FIXED TRANSITIONS C W MATRIX C C(I,I,K) HAS BEEN INITIALIZED TO ZERO IN COLRAT C AN IF-STATEMENT TO SKIP THE CASES WHEN I=J IS THUS C UNNECESSARY C DO 575 K=1,NDEP DO 570 J=1,NK DO 560 I=1,NK W(I,K,I,K)=W(I,K,I,K)-C(I,J,K)*N(I,K) W(I,K,J,K)=W(I,K,J,K)+C(J,I,K)*N(J,K) 560 CONTINUE 570 CONTINUE 575 CONTINUE CALL CPUTIME('WCOLL ',0,1,1) C C PARTICLE CONSERVATION EQUATION C DO 600 K=1,NDEP IF(LISUM) ISUM=ISUMA(K) DO 590 J=1,NK DO 580 L=1,NDEP W(ISUM,K,J,L)=0. 580 CONTINUE W(ISUM,K,J,K)=N(J,K) 590 CONTINUE 600 CONTINUE ENDIF C C ERROR VECTOR. FOR I=J SEE COMMENT ABOVE C DO 700 K=1,NDEP DO 670 J=1,NK DO 660 I=1,NK E(I,K)=E(I,K)+N(I,K)*C(I,J,K)-N(J,K)*C(J,I,K) 660 CONTINUE 670 CONTINUE IF(LISUM) ISUM=ISUMA(K) E(ISUM,K)=TOTN(K) DO 690 J=1,NK E(ISUM,K)=E(ISUM,K)-N(J,K) 690 CONTINUE 700 CONTINUE IF(LISUM) ISUM=0 CALL CPUTIME('ERCOLL',0,1,1) IF(IWEVEC.GT.0) CALL WEVEC(E,NK1,NDEP1) C C OUTPUT TO FILE WMAT C IF(NEWMAT) CALL WWMAT(W,E,NK1,NDEP1,EMAX) C C CALCULATE THE NEW POPULATION NUMBERS C PLOT LOGARITHMIC RELATIVE CHANGE C NDIM=NDEP*NK MDIM=NDEP1*NK1 IF(NDIM.NE.MDIM) CALL STOP('ITER:NDIM.NE.MDIM') CALL CPUTIME('EQSYST',0,0,1) CALL EQSYST(NDIM,MDIM,W,E,NEWMAT) CALL CPUTIME('EQSYST',0,1,1) C CALL WCHANG(E,NK1,NDEP1,1) C C STORE I INDEX FOR MAX CHANGE IN IE, K INDEX IN KE C EMAX=0.0 IE=0 KE=0 DO 800 K=1,NDEP DO 780 I=1,NK IF(ABS(E(I,K)).GT.EMAX) THEN IE=I KE=K EMAX=ABS(E(I,K)) ENDIF 780 CONTINUE 800 CONTINUE EMAX=E(IE,KE) WRITE(ETEXT,810) IE,KE 810 FORMAT('EMAX (',I3,',',I3,')') CALL WNIIT(1) CALL WEMAX(' ',EMAX) IF(ICRSW.NE.0) CALL WEMAX('SWITCH ',SWITCH) CALL WEMAX(ETEXT,EMAX) EMAX=ABS(EMAX) C C IF SWITCHING IS ENABLED IN AUTOMATIC MODE: C CHECK TO SEE IF REDO IS TO BE DONE C IF(ICRSW.EQ.-2) THEN IF(LOG10(EMAX).GT.EREDO .AND. OLDSW.NE.1.0) THEN DIVFAC=MAX(ONE+0.1,(1.1+(LOG10(EMAX)-EREDO)*2.)) FS0=FS0/DIVFAC IF(FS0.LT.1.E-4) THEN CALL STOP('ITER: SWITCHING NOT CONVERGING') ENDIF SWITCH=MAX(ONE,OLDSW/(1.+FS0)) IF(IWEMAX.NE.0) THEN WRITE(*,815) FS0 WRITE(LOUT,815) FS0 815 FORMAT(' REDOING ITERATION, FS0= ',1P,E10.2,0P) ENDIF GOTO 9000 ENDIF ENDIF C IF(ICRSW.NE.0) OLDSW=SWITCH IF(ICRSW.EQ.-2) THEN FS0=MIN(DFSMAX,(1.+(-LOG10(EMAX)+EREDO)*0.2))*FS0 SWITCH=MAX(ONE,SWITCH/(1.+FS0)) ELSE IF(ICRSW.EQ.-1) THEN WRITE(*,*) ' GIVE NEW SWITCH VALUE (<0 TO REDO ITERATION)' READ(*,*) SWITCH IF(SWITCH.LT.0.) THEN SWITCH=-SWITCH GOTO 9000 ELSE IF(SWITCH.EQ.0.) THEN WRITE(*,*) ' SWITCHING TO AUTOMATIC MODE' WRITE(LOUT,*) ' SWITCHING TO AUTOMATIC MODE' SWITCH=OLDSW ICRSW=-2 ENDIF ELSE IF(ICRSW.GT.0) THEN IF(ABS(EMAX).LT.10.) THEN SWITCH=MAX(SWITCH*SWFAC,ONE) ELSE IF(ABS(EMAX).LT.100.) THEN SWITCH=MAX(SWITCH/100.*ABS(EMAX),ONE) ENDIF ENDIF C C UPDATE POPULATIONS C DO 830 K=1,NDEP DO 820 I=1,NK N(I,K)=N(I,K)*(1.0+E(I,K)) 820 CONTINUE 830 CONTINUE C IF(HSE .AND. EMAX.LT.ELIM1 .AND. SWITCH.EQ.1.) THEN CALL HSEINT(HSE) IF(ITHSE.EQ.1 .AND. NK+1.LE.MK) THEN C C NG ACCELERATION FOR HSE ITERATIONS. C SWITCHED ON IF HSE STEP EACH ITERATION (ITHSE.EQ.1) C ZERO COUNT IN NG IF NGHSE=.FALSE. C ACCELERATION ON BOTH N AND NE THROUGH LOCAL ARRAY AN C TO HAVE ROOM FOR NE, NK+1.LE.MK IS REQUIRED C IF(.NOT.NGHSE) CALL NG(N,NK,NDEP,NGHSE) NGHSE=.TRUE. DO 910 I=1,NK DO 900 K=1,NDEP AN(I,K)=N(I,K) 900 CONTINUE 910 CONTINUE DO 920 K=1,NDEP AN(NK+1,K)=NE(K) 920 CONTINUE CALL NG(AN,NK+1,NDEP,NGHSE) DO 940 I=1,NK DO 930 K=1,NDEP N(I,K)=AN(I,K) 930 CONTINUE 940 CONTINUE DO 950 K=1,NDEP NE(K)=AN(NK+1,K) 950 CONTINUE ELSE NGHSE=.FALSE. ENDIF EMAX=10. ITHSE=0 ENDIF CALL CPUTIME('ITER ',0,2,2) IF(EMAX.LT.ELIM2 .AND. SWITCH.EQ.1.) THEN ICONV=1 GOTO 9100 ENDIF C C NG ACCELERATION STEP C START NG SCHEME AS SOON AS CORRECTIONS FALL BELOW 100 PERCENT C AND SWITCH=1.0. NOTE THAT NG COUNT IS RESET TO 1 AT EACH HSE C ITERATION C IF(NGACC .AND. .NOT.NGHSE) THEN NGON=EMAX.LT.1.0 .AND. SWITCH.EQ.ONE CALL NG(N,NK,NDEP,NGON) ENDIF C 9000 CONTINUE C 9100 CALL CPUTIME('SUM IT',0,2,3) RETURN END C C********************************************************************** C SUBROUTINE ZMAT(W,N) C C SETS ARRAY W(N) TO ZERO C C: ZMAT 90-07-31 NEW ROUTINE: (MATS CARLSSON) C: INCLUDE 'PREC' DIMENSION W(N) C DO 100 I=1,N W(I)=0.0 100 CONTINUE C RETURN END C C*********************************************************************** C SUBROUTINE NG(AN,NK,NDEP,NGON) C C PERFORMS NG ACCELERATION C SEE L.H. AUER , P. 101 C IN KALKOFEN, ED., "NUMERICAL RADIATIVE TRANSFER", C CAMBRIDGE UNIVERSITY PRESS 1987 C C: NG 90-09-11 NEW ROUTINE: (MARTIN J. STIFT, MATS CARLSSON) C: C INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION AN(MK,MDEP) LOGICAL NGON C DIMENSION YS(5,MK,MDEP),AA(3,3),BB(3) DATA ICALL/0/ SAVE YS,ICALL C C IF NGON=.FALSE., RESET ITERATION COUNT AND RETURN C IF(.NOT.NGON) THEN ICALL=0 RETURN ENDIF ICALL=ICALL+1 C C STORE ITERATED LEVEL POPULATIONS FOR NG ACCELERATION C IS0 = MOD(ICALL-1,5) + 1 DO 910 I = 1,NK DO 900 K = 1,NDEP YS(IS0,I,K) = AN(I,K) 900 CONTINUE 910 CONTINUE C IF (IS0.EQ.5) THEN C IF(IWEMAX.NE.0) THEN WRITE(*,*) ' NG ACCELERATION' WRITE(LOUT,*) ' NG ACCELERATION' ENDIF DO 960 I = 1,NK C DO 930 K1 = 1,3 DO 920 K2 = 1,3 AA(K1,K2) = 0. 920 CONTINUE BB(K1) = 0. 930 CONTINUE C DO 940 K = 1,NDEP WT = 1. / YS(5,I,K)**2 D0 = YS(5,I,K) - YS(4,I,K) D1 = D0 - YS(4,I,K) + YS(3,I,K) D2 = D0 - YS(3,I,K) + YS(2,I,K) D3 = D0 - YS(2,I,K) + YS(1,I,K) AA(1,1) = AA(1,1) + WT * D1 * D1 AA(1,2) = AA(1,2) + WT * D1 * D2 AA(1,3) = AA(1,3) + WT * D1 * D3 AA(2,2) = AA(2,2) + WT * D2 * D2 AA(2,3) = AA(2,3) + WT * D2 * D3 AA(3,3) = AA(3,3) + WT * D3 * D3 AA(2,1) = AA(1,2) AA(3,1) = AA(1,3) AA(3,2) = AA(2,3) BB(1) = BB(1) + WT * D0 * D1 BB(2) = BB(2) + WT * D0 * D2 BB(3) = BB(3) + WT * D0 * D3 940 CONTINUE C CALL LINEQ (AA,BB,3) DO 950 K = 1,NDEP AN(I,K) = (1. - BB(1) - BB(2) - BB(3)) * YS(5,I,K) + * BB(1) * YS(4,I,K) + * BB(2) * YS(3,I,K) + * BB(3) * YS(2,I,K) 950 CONTINUE C 960 CONTINUE ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE LINEQ(A,B,N) C C FINDS SOLUTION OF SYSTEM OF LINEAR EQUATIONS C WITH GAUSSIAN ELIMINATION WITH PIVOTING C C: LINEQ 90-06-05 NEW ROUTINE: (MARTIN J STIFT) C: INCLUDE 'PREC' C DIMENSION A(3,3),B(3),C(3),ICOL(3) C C INITIALIZE COLUMN COUNT AND STARTING POINTS C DO 10 I = 1,N ICOL(I) = I 10 CONTINUE C IBEG = 1 JBEG = 1 C C DETERMINE PIVOT C 20 AMA = ABS(A(IBEG,JBEG)) IMA = IBEG JMA = JBEG DO 30 I = IBEG,N DO 30 J = JBEG,N IF(ABS(A(I,J)).LE.AMA) GOTO 30 AMA = ABS(A(I,J)) IMA = I JMA = J 30 CONTINUE C C ORDER MATRIX DEPENDING ON PIVOT C DO 40 I = 1,N TEMP = A(I,JMA) A(I,JMA) = A(I,JBEG) A(I,JBEG) = TEMP 40 CONTINUE C DO 50 J = JBEG,N TEMP = A(IMA,J) A(IMA,J) = A(IBEG,J) A(IBEG,J) = TEMP 50 CONTINUE C TEMP = B(IMA) B(IMA) = B(IBEG) B(IBEG) = TEMP IT = ICOL(JBEG) ICOL(JBEG) = ICOL(JMA) ICOL(JMA)= IT IBEG = IBEG + 1 JBEG = JBEG + 1 C C ELIMINATE C IMIN = IBEG - 1 JMIN = JBEG - 1 DO 70 I = IBEG,N QUOT = A(I,JMIN) / A(IMIN,JMIN) DO 60 J = JBEG,N A(I,J) = A(I,J) - QUOT * A(IMIN,J) 60 CONTINUE B(I) = B(I) - QUOT * B(IMIN) 70 CONTINUE IF (IBEG.LT.N) GOTO 20 C C DETERMINE COEFFICIENTS C B(N) = B(N) / A(N,N) N1 = N - 1 DO 90 I = N1,1,-1 I1 = I + 1 DO 80 J = N,I1,-1 B(I) = B(I) - A(I,J) * B(J) 80 CONTINUE B(I) = B(I) / A(I,I) 90 CONTINUE C C REORDER COEFFICIENTS C DO 100 I = 1,N C(ICOL(I))=B(I) 100 CONTINUE C DO 110 I = 1,N B(I) = C(I) 110 CONTINUE C RETURN END C C********************************************************************** C SUBROUTINE WWMAT(W0,E0,NK1,NDEP1,EMAX) C C WRITES MATRIX W TO UNFORMATTED FILE WMAT C THE MATRIX IS OUTPUT COLUMN BY COLUMN WITH ONE RECORD C PER COLUMN. C THE FIRST RECORD IS AN ID RECORD CONTAINING: C NDEP,NK,NLINE,NRAD,(IRAD(KR),KR=1,NRAD),(JRAD(KR),KR=1,NRAD), C ISUM,ATOMID,ATMOID,DPID C C IF IWWMAT.GT.0 THE FILE IS WRITTEN ONLY WHEN EMAX.LT.ELIM1*100 C TO AVOID TOO MANY WRITE OPERATIONS. IF EMAX GOES FROM A VALUE C LARGER THAN ELIM1*100 TO A VALUE SMALLER THAN ELIM1 IN ONE C ITERATION THIS MEANS THAT THERE WILL BE NO WMAT OUTPUT. C: C: WWMAT 88-02-26 MODIFICATIONS: (MATS CARLSSON) C: WRITES RIGHT HAND SIDE TO FILE IN LAST RECORD C: C: 90-07-31 MODIFICATIONS: (MATS CARLSSON) C: CHANGED TO CONTAIN DECLARATION OF W AND E C: THESE PASSED IN ARGUMENT LIST INSTEAD OF IN COMMON TO C: AVOID NK=MK1 AND NDEP=MDEP1 REQUIREMENT C: C: 90-10-19 MODIFICATIONS: (MATS CARLSSON) C: DOUBLE PRECISION VERSION. OUTPUTS MATRIX IN SINGLE C: PRECISION C: INCLUDE 'PREC' INCLUDE 'PARAM' INCLUDE 'CATOM' INCLUDE 'CATMOS' INCLUDE 'CINPUT' INCLUDE 'CLU' C DIMENSION W0(NK1*NDEP1,NK1*NDEP1),E0(NK1*NDEP1) C IF(IWWMAT.EQ.0) RETURN C CALL CPUTIME('WWMAT ',0,0,1) IF(IWWMAT.GT.0) THEN IF(EMAX.LT.ELIM1*100.) THEN CALL REWIND(LWMAT) WRITE(LWMAT) NDEP,NK,NLINE,NRAD,ISUM,(IRAD(KR),KR=1,NRAD), * (JRAD(KR),KR=1,NRAD),ATOMID,ATMOID,DPID DO 100 JL=1,NK*NDEP WRITE(LWMAT) (REAL(W0(IK,JL)),IK=1,NK*NDEP) 100 CONTINUE WRITE(LWMAT) (REAL(E0(IK)),IK=1,NK*NDEP) ENDIF ELSE IF(IT.EQ.1) THEN WRITE(LWMAT) NDEP,NK,NLINE,NRAD,ISUM,(IRAD(KR),KR=1,NRAD), * (JRAD(KR),KR=1,NRAD),ATOMID,ATMOID,DPID ENDIF DO 200 JL=1,NK*NDEP WRITE(LWMAT) (REAL(W0(IK,JL)),IK=1,NK*NDEP) 200 CONTINUE WRITE(LWMAT) (REAL(E0(IK)),IK=1,NK*NDEP) ENDIF CALL CPUTIME('WWMAT ',0,1,1) C RETURN END C C*********************************************************************** C